Add a competitive exchange model to MCT

src/libcadet/CMakeLists.txt

@@ -103,6 +103,7 @@ set(LIBCADET_BINDINGMODEL_SOURCES
 )
 set(LIBCADET_EXCHANGEMODEL_SOURCES
 	${CMAKE_SOURCE_DIR}/src/libcadet/model/exchange/LinearExchange.cpp
+	${CMAKE_SOURCE_DIR}/src/libcadet/model/exchange/LangumirExchange.cpp
 )
 
 

src/libcadet/ConfigurationHelper.hpp

@@ -59,6 +59,12 @@ class IConfigHelper
 	 */
 	virtual model::IBindingModel* createBindingModel(const std::string& name) const = 0;
 
+	/**
+	 * @brief Creates an IExchange object of the given @p name
+	 * @details The caller owns the returned IExchange object.
+	 * @param [in] name Name of the IExchange object
+	 * @return Object of the given IExchange @p name or @c nullptr if that name does not exist
+	 */
 	virtual model::IExchangeModel* createExchangeModel(const std::string& name) const = 0;
 
 

src/libcadet/ExchangeModelFactory.cpp

@@ -21,13 +21,15 @@ namespace cadet
 		namespace exchange
 		{
 			void registerLinearExModel(std::unordered_map<std::string, std::function<model::IExchangeModel* ()>>& exchange);
+			void registerLangumirExModel(std::unordered_map<std::string, std::function<model::IExchangeModel* ()>>& exchange);
 		}
 	}
 
 	ExchangeModelFactory::ExchangeModelFactory()
 	{
 		// Register all ExchangeModels here
 		model::exchange::registerLinearExModel(_exchangeModels);
+		model::exchange::registerLangumirExModel(_exchangeModels);
 
 	}
 

src/libcadet/model/MultiChannelTransportModel.cpp

@@ -11,6 +11,7 @@
 // =============================================================================
 
 #include "model/MultiChannelTransportModel.hpp"
+#include "ExchangeModelFactory.hpp"
 #include "ParamReaderHelper.hpp"
 #include "ParamReaderScopes.hpp"
 #include "cadet/Exceptions.hpp"
@@ -411,7 +412,6 @@ bool MultiChannelTransportModel::configureModelDiscretization(IParameterProvider
 
 
 	clearExchangeModels();
-
 	_exchange.push_back(nullptr);
 
 	if (paramProvider.exists("EXCHANGE_MODEL"))

src/libcadet/model/exchange/LangumirExchange.cpp

@@ -0,0 +1,303 @@
+// =============================================================================
+//  CADET
+//  
+//  Copyright © 2008-present: The CADET-Core Authors
+//            Please see the AUTHORS.md file.
+//  
+//  All rights reserved. This program and the accompanying materials
+//  are made available under the terms of the GNU Public License v3.0 (or, at
+//  your option, any later version) which accompanies this distribution, and
+//  is available at http://www.gnu.org/licenses/gpl.html
+// =============================================================================
+
+
+
+#include "model/ExternalFunctionSupport.hpp"
+#include "ParamIdUtil.hpp"
+#include "model/ModelUtils.hpp"
+#include "cadet/Exceptions.hpp"
+#include "model/Parameters.hpp"
+#include "LocalVector.hpp"
+#include "SimulationTypes.hpp"
+#include "Memory.hpp"
+#include "model/ExchangeModel.hpp"
+#include "ParamReaderHelper.hpp"
+#include "model/parts/MultiChannelConvectionDispersionOperator.hpp"
+
+#include <vector>
+#include <unordered_map>
+#include <functional>
+#include <algorithm>
+#include <iterator>
+#include <tuple>
+
+namespace cadet
+{
+
+namespace model
+{
+/**
+ * @brief Defines the linear exchange model
+ * @details Implements the linear exchange model for the MCT model
+ * The exchange is given by a matrix of exchange coefficients eij for each component i and j and the chross Ai section from channel i.
+ * The exchange from channel i to all other channel j is given by \f$ \frac{\mathrm{d}ci}{\mathrm{d}t} = \sum_j eij cj Aj/Ai - eji ci \f$.
+ */
+class LangumirExchangeBase : public IExchangeModel
+{
+public:
+
+	LangumirExchangeBase() : _nComp(0), _nChannel(0), _nCol(0) { }
+	virtual ~LangumirExchangeBase() CADET_NOEXCEPT { }
+
+	static const char* identifier() { return "LANGMUIR_EX"; }
+	virtual const char* name() const CADET_NOEXCEPT { return "LANGMUIR_EX"; }
+	virtual bool requiresConfiguration() const CADET_NOEXCEPT { return true; }
+	virtual bool usesParamProviderInDiscretizationConfig() const CADET_NOEXCEPT { return true; }
+
+	virtual bool configureModelDiscretization(IParameterProvider& paramProvider, unsigned int nComp, unsigned int nChannel, unsigned int nCol)
+	{
+		_nComp = nComp;
+		_nChannel = nChannel;
+		_nCol = nCol;
+
+		return true;
+	}
+
+	virtual bool configure(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx)
+	{
+		_parameters.clear();
+		readParameterMatrix(_exchangeMatrix, paramProvider, "EXCHANGE_MATRIX", _nChannel * _nChannel * _nComp, 1); // include parameterPeaderHelp in exchange modul
+		_crossSections = paramProvider.getDoubleArray("CHANNEL_CROSS_SECTION_AREAS");
+		readParameterMatrix(_capacityMatrix, paramProvider, "CAPACITY_MATRIX", _nComp * _nChannel, 1);
+
+		return true;
+	}
+
+	virtual void fillExchangeInitialParameters(ParameterId* params, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx) const CADET_NOEXCEPT
+	{
+		unsigned int ctr = 0;
+		for (int c = 0; c < _nComp; ++c)
+		{
+			for (unsigned int bp = 0; bp < _nChannel; ++bp, ++ctr)
+				params[ctr] = makeParamId(hashString("INIT_C"), unitOpIdx, c, parTypeIdx, bp, ReactionIndep, SectionIndep);
+		}
+	}
+	// ----------------------------------------------------------------------------------------------------------------------------// 
+	virtual std::unordered_map<ParameterId, double> getAllParameterValues() const
+	{
+		std::unordered_map<ParameterId, double> data;
+		std::transform(_parameters.begin(), _parameters.end(), std::inserter(data, data.end()),
+			[](const std::pair<const ParameterId, active*>& p) { return std::make_pair(p.first, static_cast<double>(*p.second)); });
+
+		return data;
+	}
+
+	virtual bool hasParameter(const ParameterId& pId) const
+	{
+		return _parameters.find(pId) != _parameters.end();
+	}
+
+	virtual bool setParameter(const ParameterId& pId, int value)
+	{
+		return false;
+	}
+
+	virtual bool setParameter(const ParameterId& pId, double value)
+	{
+		auto paramHandle = _parameters.find(pId);
+		if (paramHandle != _parameters.end())
+		{
+			paramHandle->second->setValue(value);
+			return true;
+		}
+
+		return false;
+	}
+
+	virtual bool setParameter(const ParameterId& pId, bool value)
+	{
+		return false;
+	}
+
+	virtual active* getParameter(const ParameterId& pId) {
+		auto paramHandle = _parameters.find(pId);
+		if (paramHandle != _parameters.end())
+		{
+			return paramHandle->second;
+		}
+
+		return nullptr;
+	}
+
+	virtual unsigned int workspaceSize(unsigned int nComp, unsigned int totalNumBoundStates, unsigned int const* nChannel) const CADET_NOEXCEPT
+	{
+		// return _paramHandler.cacheSize(nComp, totalNumBoundStates, nChannel);
+		return 0;
+	}
+
+	virtual void configure(IExternalFunction** extFuns, unsigned int size) {}
+
+	virtual int residual(active const* y, active* res, WithParamSensitivity, bool wantJac, linalg::BandedSparseRowIterator jacBegin) const
+	{
+		if (wantJac)
+			return residualImpl<double, active, active, true>(reinterpret_cast<const double*>(y), res, jacBegin);
+		else
+			return residualImpl<active, active, active, false>(y, res, jacBegin);
+	}
+
+	virtual int residual(active const* y, active* res, WithoutParamSensitivity, bool wantJac, linalg::BandedSparseRowIterator jacBegin) const
+	{	
+		if (wantJac)
+			return residualImpl<double, active, double, false>(reinterpret_cast<const double*>(y), res, jacBegin);
+		else
+			return residualImpl<active, active, double, false>(y, res, jacBegin);
+	}
+
+	virtual int residual(double const* y, active* res, WithParamSensitivity, bool wantJac, linalg::BandedSparseRowIterator jacBegin) const
+	{
+		if(wantJac)
+			return residualImpl<double, active, active, true>(y, res, jacBegin);
+		else
+			return residualImpl<double, active, active, false>(y, res, jacBegin);
+	}
+
+	virtual int residual(double const* y, double* res, WithoutParamSensitivity, bool wantJac, linalg::BandedSparseRowIterator jacBegin) const
+	{
+		if (wantJac)
+			return residualImpl<double, double, double, true>(y, res, jacBegin);
+		else
+			return residualImpl<double, double, double, false>(y, res, jacBegin);
+	}
+
+
+protected:
+	int _nComp; //!< Number of components
+	unsigned int _nChannel; //!< Total number of bound states
+	unsigned int _nCol; //!< Number of columns
+
+	std::vector<active> _exchangeMatrix; //!< Matrix of exchange coeffs for the linear inter-channel transport
+	std::vector<double> _crossSections; //!< Cross sections of the channels
+	std::vector<active> _capacityMatrix; //!< Capacity of the channels -> double  ncomp x nchannel
+	parts::MultiChannelConvectionDispersionOperator _conDis; //!< Convection dispersion operator
+
+	std::unordered_map<ParameterId, active*> _parameters; //!< Map used to translate ParameterIds to actual variables
+
+	virtual bool implementsAnalyticJacobian() const CADET_NOEXCEPT { return true; }
+
+
+	template <typename StateType, typename ResidualType, typename ParamType, bool wantJac>
+	int residualImpl(StateType const* y, ResidualType* res, linalg::BandedSparseRowIterator jacBegin) const
+	{	
+
+		const unsigned int offsetC = _nChannel * _nComp;
+		for (unsigned int col = 0; col < _nCol; ++col) // Collum
+		{
+			const unsigned int offsetColBlock = col * _nChannel * _nComp;
+			ResidualType* const resColBlock = res + offsetC + offsetColBlock;
+			StateType const* const yColBlock = y + offsetC + offsetColBlock;
+
+			for (unsigned int rad_orig = 0; rad_orig < _nChannel; ++rad_orig)  // Channel orig
+			{
+				const unsigned int offsetToRadOrigBlock = rad_orig * _nComp;
+				const unsigned int offsetColRadOrigBlock = offsetColBlock + offsetToRadOrigBlock;
+				ResidualType* const resColRadOrigBlock = resColBlock + offsetToRadOrigBlock;
+				StateType const* const yColRadOrigBlock = yColBlock + offsetToRadOrigBlock;
+
+				for (unsigned int rad_dest = 0; rad_dest < _nChannel; ++rad_dest) // Channel dest
+				{
+					if (rad_orig == rad_dest)
+						continue;
+
+					const unsigned int offsetToRadDestBlock = rad_dest * _nComp;
+					const unsigned int offsetColRadDestBlock = offsetColBlock + offsetToRadDestBlock;
+					ResidualType* const resColRadDestBlock = resColBlock + offsetToRadDestBlock;
+					StateType const* const yColRadDestBlock = yColBlock + offsetToRadDestBlock;
+
+					for (unsigned int comp = 0; comp < _nComp; ++comp) // Component
+					{
+						const unsigned int offsetCur_orig = offsetColRadOrigBlock + comp; // component in orig channel
+						const unsigned int offsetCur_dest = offsetColRadDestBlock + comp; // component in dest channel
+						StateType const* const yCur_orig = yColRadOrigBlock + comp; // state of component in orig channel
+						//StateType const* const yCur_dest = yColRadDestBlock + comp; // state of component in dest channel
+						ResidualType* const resCur_orig = resColRadOrigBlock + comp; // residual of component in orig channel
+						ResidualType* const resCur_dest = resColRadDestBlock + comp; // residual of component in dest channel
+
+						const ParamType exchange_orig_dest_comp = static_cast<ParamType>(_exchangeMatrix[rad_orig * _nChannel * _nComp + rad_dest * _nComp + comp]);
+						const ParamType capacity_orig_comp = static_cast<ParamType>(_capacityMatrix[rad_dest * _nComp + comp]);
+						//const ParamType capacity_dest_comp = static_cast<ParamType>(_capacityMatrix[rad_dest * _nComp + comp]);
+						if (cadet_likely(exchange_orig_dest_comp > 0.0))
+						{
+							double compSum_orig_dest = 0.0;
+							for (unsigned int component = 0; component < _nComp; component++)
+							{
+
+								StateType const* const compSum = yColRadDestBlock + component;
+								double compSum_dest = static_cast<double>(compSum[0]);
+
+								double cMax_destComp  = static_cast<double>(_capacityMatrix[rad_dest * _nComp + component]);
+
+								if (cMax_destComp < 1e-16)
+								{
+									continue;
+								}
+
+								compSum_orig_dest += compSum_dest / cMax_destComp;
+							 }
+
+							*resCur_orig += exchange_orig_dest_comp * yCur_orig[0] * capacity_orig_comp * (1 - compSum_orig_dest);
+							if (capacity_orig_comp < 1e-16)
+							{
+								*resCur_orig += exchange_orig_dest_comp * yCur_orig[0];
+							}
+
+							*resCur_dest -= exchange_orig_dest_comp * yCur_orig[0] * static_cast<ParamType>(_crossSections[rad_orig]) / static_cast<ParamType>(_crossSections[rad_dest]) * capacity_orig_comp * (1 - compSum_orig_dest);
+							if (capacity_orig_comp < 1e-16)
+							{
+								*resCur_dest -= exchange_orig_dest_comp * yCur_orig[0] * static_cast<ParamType>(_crossSections[rad_orig]) / static_cast<ParamType>(_crossSections[rad_dest]);
+							}
+
+
+							if (wantJac) {
+
+								linalg::BandedSparseRowIterator jacorig;
+								jacorig = jacBegin + offsetCur_orig;
+								jacorig[0] += static_cast<double>(exchange_orig_dest_comp) * static_cast<double>(capacity_orig_comp * (1 - compSum_orig_dest));
+								if (capacity_orig_comp < 1e-16)
+								{
+									jacorig[0] -= static_cast<double>(exchange_orig_dest_comp);
+								}
+
+								linalg::BandedSparseRowIterator jacdest;
+								jacdest = jacBegin + offsetCur_dest;
+								jacdest[static_cast<int>(offsetCur_orig) - static_cast<int>(offsetCur_dest)] -= static_cast<double>(exchange_orig_dest_comp) * static_cast<double>(capacity_orig_comp) * (1 - compSum_orig_dest);
+								if (capacity_orig_comp < 1e-16)
+								{
+									jacdest[static_cast<int>(offsetCur_orig) - static_cast<int>(offsetCur_dest)] += static_cast<double>(exchange_orig_dest_comp);
+								}
+							}
+
+						}
+
+					}
+
+				}
+
+			}
+		}
+		return 0;
+	}
+};
+
+typedef LangumirExchangeBase  LangumirExchange;
+
+namespace exchange
+{
+	void registerLangumirExModel(std::unordered_map<std::string, std::function<model::IExchangeModel*()>>& exchange)
+	{
+		exchange[LangumirExchange::identifier()] = []() { return new  LangumirExchange(); };
+	}
+}  // namespace exchange
+
+} // namespace model
+
+}// namespace cadet