Add Multicomponent Freundlich LDF isotherm. (#259)

Co-authored-by: Ronald Jäpel <[email protected]>

doc/interface/binding/index.rst

@@ -68,6 +68,7 @@ This group also takes precedence over a possibly existing ``/input/model/unit_XX
     saska
     self_association
     freundlich_ldf
+    multi_component_ldf_freundlich
     hic_water_on_hydrophobic_surfaces
     hic_constant_water_activity
     multi_component_colloidal

doc/interface/binding/multi_component_ldf_freundlich.rst

@@ -0,0 +1,63 @@
+.. _multi_component_ldf_freundlich_config:
+
+Multi Component Linear Driving Force Freundlich
+===============================================
+
+**Group /input/model/unit_XXX/adsorption – ADSORPTION_MODEL = MULTI_COMPONENT_LDF_FREUNDLICH**
+
+For information on model equations, refer to :ref:`multi_component_ldf_freundlich_model`.
+
+``IS_KINETIC``
+   Selects kinetic or quasi-stationary adsorption mode: 1 = kinetic, 0 =
+   quasi-stationary. If a single value is given, the mode is set for all
+   bound states. Otherwise, the adsorption mode is set for each bound
+   state separately.
+
+===================  =========================  =========================================
+**Type:** int        **Range:** {0,1}  		    **Length:** 1/NTOTALBND
+===================  =========================  =========================================
+
+``MCLDFFRL_KLDF``
+   Rate constants in linear driving force approach
+
+**Unit:** :math:`s^{-1}`
+
+===================  =========================  =========================================
+**Type:** double     **Range:** :math:`\ge 0`   **Length:** NCOMP
+===================  =========================  =========================================
+
+``MCLDFFRL_KF``
+   Proportionality constants
+
+**Unit:** :math:`m_{MP}^{3/n}~m_{SP}^{-3}~mol^{1-1/n}`
+
+===================  =========================  ==================================
+**Type:** double     **Range:** :math:`\ge 0`   **Length:** NCOMP
+===================  =========================  ================================== 
+
+``MCLDFFRL_EXP``
+   Freundlich exponent
+
+**Unit:** [-]
+
+===================  =========================  ==================================
+**Type:** double     **Range:** :math:`\gt 0`   **Length:** NCOMP
+===================  =========================  ================================== 
+
+``MCLDFFRL_A``
+   Component influences in row-major ordering
+
+**Unit:** [-]
+
+===================  =========================  ==================================
+**Type:** double     **Range:** :math:`\ge 0`   **Length:** :math:`\text{NCOMP}^2`
+===================  =========================  ================================== 
+
+``MCLDFFRL_TAU``
+   Small constant that ensures numerical stability
+
+**Unit:** :math:`mol~m_{MP}^{-3}`
+
+===================  =========================  =========================================
+**Type:** double     **Range:** :math:`\gt 0`   **Length:** 1
+===================  =========================  =========================================

doc/modelling/binding/freundlich_ldf.rst

@@ -14,7 +14,7 @@ This variant of the model is based on the linear driving force approximation (se
 No interaction between the components is considered when the model has multiple components. 
 One of the limitation of this isotherm is the first order Jacobian :math:`\left(\frac{dq^*}{dc_p}\right)` tends to infinity as :math:`c_{p} \rightarrow 0` for :math:`n>1`.
 To address this issue an approximation of isotherm is considered near the origin.
-This approximation matches the isotherm in such a way that  :math:`q=0` at :math:`c_p=0` and also matches the first derivative of the istotherm at :math:`c_p = \varepsilon`, where :math:`\varepsilon` is a very small number, for example :math:`1e-14`.
+This approximation matches the isotherm in such a way that  :math:`q=0` at :math:`c_p=0` and also matches the first derivative of the isotherm at :math:`c_p = \varepsilon`, where :math:`\varepsilon` is a very small number, for example :math:`1e-14`.
 The form of approximation and its derivative is given below for :math:`c_p < \varepsilon` and :math:`n>1`:
 
 .. math::
@@ -38,7 +38,7 @@ where :math:`\alpha_0=0` and :math:`\alpha_1` and :math:`\alpha_2` are determine
 		\alpha_2 = \frac{1-n}{n}k_f \varepsilon^{\frac{1-2 n}{n}}
 	\end{aligned}
 
-This approximation can be used for any pore phase cocentration :math:`c_p < \varepsilon` given :math:`n>1`.
-For the case, when :math:`n \le 1` no special treament near the origin is required.
+This approximation can be used for any pore phase concentration :math:`c_p < \varepsilon` given :math:`n>1`.
+For the case, when :math:`n \le 1` no special treatment near the origin is required.
 For more information on model parameters required to define in CADET file format, see :ref:`freundlich_ldf_config`.
 

doc/modelling/binding/index.rst

@@ -277,6 +277,11 @@ The models also differ in whether a mobile phase modifier (e.g., salt) is suppor
      - x
      - ✓
      - x
+   * - :ref:`multi_component_ldf_freundlich_model`
+     - ✓
+     - ×
+     - ✓
+     - ×
 
 
 .. toctree::

doc/modelling/binding/multi_component_ldf_freundlich.rst

@@ -0,0 +1,20 @@
+.. _multi_component_ldf_freundlich_model:
+
+Multi Component Linear Driving Force Freundlich
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+A multi-component extension to the classical Freundlich adsorption model.
+A linear driving force approach is applied to obtain a kinetic form.
+
+.. math::
+
+    \begin{aligned}
+        \frac{\mathrm{d} q_i}{\mathrm{d} t} &= k_{\text{ldf},i}\: \left(q_i^* - q_i \right) & i = 0, \dots, N_{\text{comp}} - 1 \\
+        q_i^* &= k_{F,i} c_{p,i} \left( \sum_j a_{ij} c_{p,j} + \tau \right)^{1 / n_i - 1}.
+    \end{aligned}
+
+Here, :math:`\tau > 0` is a small constant that ensures numerical stability.
+
+In a rapid-equilibrium setting with a diagonal matrix (i.e., :math:`a_{ii} = 1` and :math:`a_{ij} = 0` for :math:`j \neq i`), the traditional Freundlich isotherm is recovered.
+
+For more information on model parameters required to define in CADET file format, see :ref:`multi_component_ldf_freundlich_config`.

src/libcadet/BindingModelFactory.cpp

@@ -45,6 +45,7 @@ namespace cadet
 			void registerBiLangmuirLDFModel(std::unordered_map<std::string, std::function<model::IBindingModel* ()>>& bindings);
 			void registerHICWaterOnHydrophobicSurfacesModel(std::unordered_map<std::string, std::function<model::IBindingModel*()>>& bindings);
 			void registerHICConstantWaterActivityModel(std::unordered_map<std::string, std::function<model::IBindingModel*()>>& bindings);
+			void registerMultiComponentLDFFreundlichModel(std::unordered_map<std::string, std::function<model::IBindingModel*()>>& bindings);
 		}
 	}
 
@@ -75,6 +76,7 @@ namespace cadet
 		model::binding::registerBiLangmuirLDFModel(_bindingModels);
 		model::binding::registerHICWaterOnHydrophobicSurfacesModel(_bindingModels);
 		model::binding::registerHICConstantWaterActivityModel(_bindingModels);
+		model::binding::registerMultiComponentLDFFreundlichModel(_bindingModels);
 		registerModel<model::SimplifiedMultiStateStericMassActionBinding>();
 	}
 

src/libcadet/CMakeLists.txt

@@ -100,6 +100,7 @@ set(LIBCADET_BINDINGMODEL_SOURCES
 	${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/BiLangmuirLDFBinding.cpp
 	${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/HICWaterOnHydrophobicSurfacesBinding.cpp
 	${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/HICConstantWaterActivityBinding.cpp
+	${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/MultiComponentLDFFreundlichBinding.cpp
 )
 set(LIBCADET_EXCHANGEMODEL_SOURCES
 	${CMAKE_SOURCE_DIR}/src/libcadet/model/exchange/LinearExchange.cpp

src/libcadet/model/binding/MultiComponentLDFFreundlichBinding.cpp

@@ -0,0 +1,207 @@
+// =============================================================================
+//  CADET
+//  
+//  Copyright © 2008-2022: The CADET Authors
+//            Please see the AUTHORS and CONTRIBUTORS 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/binding/BindingModelBase.hpp"
+#include "model/ExternalFunctionSupport.hpp"
+#include "model/ModelUtils.hpp"
+#include "cadet/Exceptions.hpp"
+#include "model/Parameters.hpp"
+#include "LocalVector.hpp"
+#include "SimulationTypes.hpp"
+
+#include <functional>
+#include <unordered_map>
+#include <string>
+#include <vector>
+
+/*<codegen>
+{
+	"name": "MultiComponentLDFFreundlichParamHandler",
+	"externalName": "ExtMultiComponentLDFFreundlichParamHandler",
+	"parameters":
+		[
+			{ "type": "ScalarComponentDependentParameter", "varName": "kLDF", "confName": "MCLDFFRL_KLDF"},
+			{ "type": "ScalarComponentDependentParameter", "varName": "kF", "confName": "MCLDFFRL_KF"},
+			{ "type": "ScalarComponentDependentParameter", "varName": "n", "confName": "MCLDFFRL_EXP"},
+			{ "type": "ComponentDependentComponentVectorParameter", "varName": "a", "confName": "MCLDFFRL_A"}
+		],
+	"constantParameters":
+		[
+			{ "type": "ScalarParameter", "varName": "tau", "confName": "MCLDFFRL_TAU"}
+		]
+}
+</codegen>*/
+
+namespace cadet
+{
+
+namespace model
+{
+
+inline const char* MultiComponentLDFFreundlichParamHandler::identifier() CADET_NOEXCEPT { return "MULTI_COMPONENT_LDF_FREUNDLICH"; }
+
+inline bool MultiComponentLDFFreundlichParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
+{
+	if ((_kLDF.size() != _kF.size()) || (_kLDF.size() != _n.size()) || (_kLDF.size() < nComp))
+		throw InvalidParameterException("MCLDFFRL_KLDF, MCLDFFRL_KF, and MCLDFFRL_N have to have the same size");
+
+	if ((_a.slices() != nComp) || (_a.size() != nComp * nComp))
+		throw InvalidParameterException("MCLDFFRL_A has to have NCOMP^2 elements");
+
+	return true;
+}
+
+inline const char* ExtMultiComponentLDFFreundlichParamHandler::identifier() CADET_NOEXCEPT { return "EXT_MULTI_COMPONENT_LDF_FREUNDLICH"; }
+
+inline bool ExtMultiComponentLDFFreundlichParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates)
+{
+	if ((_kLDF.size() != _kF.size()) || (_kLDF.size() != _n.size()) || (_kLDF.size() < nComp))
+		throw InvalidParameterException("EXT_MCLDFFRL_KLDF, EXT_MCLDFFRL_KF, and EXT_MCLDFFRL_N have to have the same size");
+
+	if ((_a.slices() != nComp) || (_a.size() != nComp * nComp))
+		throw InvalidParameterException("EXT_MCLDFFRL_A has to have NCOMP^2 elements");
+
+	return true;
+}
+
+
+/**
+ * @brief Defines the multi component Freundlich binding model with linear driving force
+ * @details Implements the multi component linear driving force Freundlich adsorption model: \f[ \begin{align} 
+ *              \frac{\mathrm{d}q_i}{\mathrm{d}t} &= k_{ldf,i} \left( q_i^* - q_i \right) \\
+ *              q_i^* &= k_{f,i} c_{p,i} \left( \sum_j a_{ij} c_{p,j} \right)^{n_i-1}
+ *          \end{align} \f]
+ *          where @f$ a_{ii} = 1 @f$ for all i.
+ *          Multiple bound states are not supported. 
+ *          Components without bound state (i.e., non-binding components) are supported.
+ * @tparam ParamHandler_t Type that can add support for external function dependence
+ */
+template <class ParamHandler_t>
+class MultiComponentLDFFreundlichBindingBase : public ParamHandlerBindingModelBase<ParamHandler_t>
+{
+public:
+
+	MultiComponentLDFFreundlichBindingBase() { }
+	virtual ~MultiComponentLDFFreundlichBindingBase() CADET_NOEXCEPT { }
+
+	static const char* identifier() { return ParamHandler_t::identifier(); }
+
+	virtual bool implementsAnalyticJacobian() const CADET_NOEXCEPT { return true; }
+
+	CADET_BINDINGMODELBASE_BOILERPLATE
+
+protected:
+	using ParamHandlerBindingModelBase<ParamHandler_t>::_paramHandler;
+	using ParamHandlerBindingModelBase<ParamHandler_t>::_reactionQuasistationarity;
+	using ParamHandlerBindingModelBase<ParamHandler_t>::_nComp;
+	using ParamHandlerBindingModelBase<ParamHandler_t>::_nBoundStates;
+
+	template <typename StateType, typename CpStateType, typename ResidualType, typename ParamType>
+	int fluxImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, StateType const* y,
+		CpStateType const* yCp, ResidualType* res, LinearBufferAllocator workSpace) const
+	{
+		using CpStateParamType = typename DoubleActivePromoter<CpStateType, ParamType>::type;
+		using StateParamType = typename DoubleActivePromoter<StateType, ParamType>::type;
+		typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace);
+
+		// Fluxes: k_{ldf,i} ( q_i - k_{f,i} c_{p,i} \left( tau + \sum_j a_{ij} c_{p,j} \right)^{n_i-1} )
+		unsigned int bndIdx = 0;
+		for (int i = 0; i < _nComp; ++i)
+		{
+			// Skip components without bound states (bound state index bndIdx is not advanced)
+			if (_nBoundStates[i] == 0)
+				continue;
+
+			active const* const aSlice = p->a[i];
+			ResidualType sum = static_cast<double>(p->tau);
+			for (int j = 0; j < _nComp; ++j)
+			{
+				sum += static_cast<ParamType>(aSlice[j]) * yCp[j];
+			}
+
+			// Residual
+			res[bndIdx] = static_cast<ParamType>(p->kLDF[i]) * ( y[bndIdx] - static_cast<ParamType>(p->kF[i]) * yCp[i] * pow(sum, 1 / static_cast<ParamType>(p->n[i]) - 1.0) );
+
+			// Next bound component
+			++bndIdx;
+		}
+
+		return 0;
+	}
+
+	template <typename RowIterator>
+	void jacobianImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double const* yCp, int offsetCp, RowIterator jac, LinearBufferAllocator workSpace) const
+	{
+		typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace);
+
+		// Fluxes: k_{ldf,i} ( q_i - k_{f,i} c_{p,i} \left( tau + \sum_j a_{ij} c_{p,j} \right)^{n_i-1} )
+
+		int bndIdx = 0;
+		for (int i = 0; i < _nComp; ++i)
+		{
+			// Skip components without bound states (bound state index bndIdx is not advanced)
+			if (_nBoundStates[i] == 0)
+				continue;
+
+			active const* const aSlice = p->a[i];
+
+			double cpSum = static_cast<double>(p->tau);
+			for (int j = 0; j < _nComp; ++j)
+			{
+				cpSum += static_cast<double>(aSlice[j]) * yCp[j];
+			}
+
+			const double kldf = static_cast<double>(p->kLDF[i]);
+			const double kf_ldf = static_cast<double>(p->kF[i]) * kldf;
+			const double n = static_cast<double>(p->n[i]);
+
+			// dres_i / dc_{p,i}
+			jac[i - bndIdx - offsetCp] = -kf_ldf * pow(cpSum, 1 / n - 1.0);
+			// Getting to c_{p,i}: -bndIdx takes us to q_0, another -offsetCp to c_{p,0} and a +i to c_{p,i}.
+			//                     This means jac[i - bndIdx - offsetCp] corresponds to c_{p,i}.
+
+			const double factor = kf_ldf * yCp[i] * ( 1 / n - 1.0) * pow(cpSum, 1 / n - 2.0);
+
+			// Fill dres_i / dc_{p,j}
+			for (int j = 0; j < _nComp; ++j)
+			{
+				// dres_i / dc_{p,j}
+				jac[j - bndIdx - offsetCp] -= factor * static_cast<double>(aSlice[j]);
+				// Getting to c_{p,j}: -bndIdx takes us to q_0, another -offsetCp to c_{p,0} and a +j to c_{p,j}.
+				//                     This means jac[j - bndIdx - offsetCp] corresponds to c_{p,j}.
+			}
+
+			// Add to dres_i / dq_i
+			jac[0] += kldf;
+
+			// Advance to next flux and Jacobian row
+			++bndIdx;
+			++jac;
+		}
+	}	
+};
+
+typedef MultiComponentLDFFreundlichBindingBase<MultiComponentLDFFreundlichParamHandler> MultiComponentLDFFreundlichBinding;
+typedef MultiComponentLDFFreundlichBindingBase<ExtMultiComponentLDFFreundlichParamHandler> ExternalMultiComponentLDFFreundlichBinding;
+
+namespace binding
+{
+	void registerMultiComponentLDFFreundlichModel(std::unordered_map<std::string, std::function<model::IBindingModel*()>>& bindings)
+	{
+		bindings[MultiComponentLDFFreundlichBinding::identifier()] = []() { return new MultiComponentLDFFreundlichBinding(); };
+		bindings[ExternalMultiComponentLDFFreundlichBinding::identifier()] = []() { return new ExternalMultiComponentLDFFreundlichBinding(); };
+	}
+}  // namespace binding
+
+}  // namespace model
+
+}  // namespace cadet