diff --git a/doc/interface/binding/index.rst b/doc/interface/binding/index.rst index a4e6e2c74..3593970ba 100644 --- a/doc/interface/binding/index.rst +++ b/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 diff --git a/doc/interface/binding/multi_component_ldf_freundlich.rst b/doc/interface/binding/multi_component_ldf_freundlich.rst new file mode 100644 index 000000000..b2adec56a --- /dev/null +++ b/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 +=================== ========================= ========================================= diff --git a/doc/modelling/binding/freundlich_ldf.rst b/doc/modelling/binding/freundlich_ldf.rst index b5a24a07d..e225f32bd 100644 --- a/doc/modelling/binding/freundlich_ldf.rst +++ b/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`. diff --git a/doc/modelling/binding/index.rst b/doc/modelling/binding/index.rst index 98211587f..617a6cc6e 100644 --- a/doc/modelling/binding/index.rst +++ b/doc/modelling/binding/index.rst @@ -272,6 +272,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:: diff --git a/doc/modelling/binding/multi_component_ldf_freundlich.rst b/doc/modelling/binding/multi_component_ldf_freundlich.rst new file mode 100644 index 000000000..53c819812 --- /dev/null +++ b/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`. diff --git a/src/libcadet/BindingModelFactory.cpp b/src/libcadet/BindingModelFactory.cpp index 2bdb14bc9..5e36f3863 100644 --- a/src/libcadet/BindingModelFactory.cpp +++ b/src/libcadet/BindingModelFactory.cpp @@ -45,6 +45,7 @@ namespace cadet void registerBiLangmuirLDFModel(std::unordered_map>& bindings); void registerHICWaterOnHydrophobicSurfacesModel(std::unordered_map>& bindings); void registerHICConstantWaterActivityModel(std::unordered_map>& bindings); + void registerMultiComponentLDFFreundlichModel(std::unordered_map>& bindings); } } @@ -75,6 +76,7 @@ namespace cadet model::binding::registerBiLangmuirLDFModel(_bindingModels); model::binding::registerHICWaterOnHydrophobicSurfacesModel(_bindingModels); model::binding::registerHICConstantWaterActivityModel(_bindingModels); + model::binding::registerMultiComponentLDFFreundlichModel(_bindingModels); registerModel(); } diff --git a/src/libcadet/CMakeLists.txt b/src/libcadet/CMakeLists.txt index 6cb580be8..874b150a4 100644 --- a/src/libcadet/CMakeLists.txt +++ b/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 ) # LIBCADET_REACTIONMODEL_SOURCES holds all source files of reaction models diff --git a/src/libcadet/model/binding/MultiComponentLDFFreundlichBinding.cpp b/src/libcadet/model/binding/MultiComponentLDFFreundlichBinding.cpp new file mode 100644 index 000000000..47f052ed4 --- /dev/null +++ b/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 +#include +#include +#include + +/* +{ + "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"} + ] +} +*/ + +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 MultiComponentLDFFreundlichBindingBase : public ParamHandlerBindingModelBase +{ +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; + using ParamHandlerBindingModelBase::_reactionQuasistationarity; + using ParamHandlerBindingModelBase::_nComp; + using ParamHandlerBindingModelBase::_nBoundStates; + + template + 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::type; + using StateParamType = typename DoubleActivePromoter::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(p->tau); + for (int j = 0; j < _nComp; ++j) + { + sum += static_cast(aSlice[j]) * yCp[j]; + } + + // Residual + res[bndIdx] = static_cast(p->kLDF[i]) * ( y[bndIdx] - static_cast(p->kF[i]) * yCp[i] * pow(sum, 1 / static_cast(p->n[i]) - 1.0) ); + + // Next bound component + ++bndIdx; + } + + return 0; + } + + template + 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(p->tau); + for (int j = 0; j < _nComp; ++j) + { + cpSum += static_cast(aSlice[j]) * yCp[j]; + } + + const double kldf = static_cast(p->kLDF[i]); + const double kf_ldf = static_cast(p->kF[i]) * kldf; + const double n = static_cast(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(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 MultiComponentLDFFreundlichBinding; +typedef MultiComponentLDFFreundlichBindingBase ExternalMultiComponentLDFFreundlichBinding; + +namespace binding +{ + void registerMultiComponentLDFFreundlichModel(std::unordered_map>& bindings) + { + bindings[MultiComponentLDFFreundlichBinding::identifier()] = []() { return new MultiComponentLDFFreundlichBinding(); }; + bindings[ExternalMultiComponentLDFFreundlichBinding::identifier()] = []() { return new ExternalMultiComponentLDFFreundlichBinding(); }; + } +} // namespace binding + +} // namespace model + +} // namespace cadet diff --git a/test/BindingModels.cpp b/test/BindingModels.cpp index c430f7b16..f2cf846b4 100644 --- a/test/BindingModels.cpp +++ b/test/BindingModels.cpp @@ -973,6 +973,57 @@ CADET_BINDINGTEST("SASKA", "EXT_SASKA", (1,1), (1,0,1), (1.0, 2.0, 0.0, 0.0), (1 1e-10, 1e-10, CADET_NONBINDING_LIQUIDPHASE_COMP_USED, CADET_COMPARE_BINDING_VS_NONBINDING) +CADET_BINDINGTEST("MULTI_COMPONENT_LDF_FREUNDLICH", "EXT_MULTI_COMPONENT_LDF_FREUNDLICH", (1,1), (1,0,1), (1.0, 2.0, 0.0, 0.0), (1.0, 3.0, 2.0, 0.0, 0.0), \ + R"json( "MCLDFFRL_KLDF": [1.14, 2.0], + "MCLDFFRL_KF": [1.3, 0.9], + "MCLDFFRL_EXP": [0.8, 3.5], + "MCLDFFRL_A": [3.0, 2.2, 1.5, 0.5], + "MCLDFFRL_TAU": 0.1 + )json", \ + R"json( "MCLDFFRL_KLDF": [1.14, 1.0, 2.0], + "MCLDFFRL_KF": [1.3, 2.0, 0.9], + "MCLDFFRL_EXP": [0.8, 3.0, 3.5], + "MCLDFFRL_A": [2.2, 1.1, 0.4, 0.1, 0.94, 2.8, 0.5, 1.2, 2.4], + "MCLDFFRL_TAU": 0.1 + )json", \ + R"json( "EXT_MCLDFFRL_KLDF": [0.0, 0.0], + "EXT_MCLDFFRL_KLDF_T": [1.14, 2.0], + "EXT_MCLDFFRL_KLDF_TT": [0.0, 0.0], + "EXT_MCLDFFRL_KLDF_TTT": [0.0, 0.0], + "EXT_MCLDFFRL_KF": [0.0, 0.0], + "EXT_MCLDFFRL_KF_T": [1.3, 0.9], + "EXT_MCLDFFRL_KF_TT": [0.0, 0.0], + "EXT_MCLDFFRL_KF_TTT": [0.0, 0.0], + "EXT_MCLDFFRL_EXP": [0.0, 0.0], + "EXT_MCLDFFRL_EXP_T": [0.8, 3.5], + "EXT_MCLDFFRL_EXP_TT": [0.0, 0.0], + "EXT_MCLDFFRL_EXP_TTT": [0.0, 0.0], + "EXT_MCLDFFRL_A": [0.0, 0.0, 0.0, 0.0], + "EXT_MCLDFFRL_A_T": [3.0, 2.2, 1.5, 0.5], + "EXT_MCLDFFRL_A_TT": [0.0, 0.0, 0.0, 0.0], + "EXT_MCLDFFRL_A_TTT": [0.0, 0.0, 0.0, 0.0], + "EXT_MCLDFFRL_TAU": 0.1 + )json", \ + R"json( "EXT_MCLDFFRL_KLDF": [0.0, 0.0, 0.0], + "EXT_MCLDFFRL_KLDF_T": [1.14, 1.0, 2.0], + "EXT_MCLDFFRL_KLDF_TT": [0.0, 0.0, 0.0], + "EXT_MCLDFFRL_KLDF_TTT": [0.0, 0.0, 0.0], + "EXT_MCLDFFRL_KF": [0.0, 0.0, 0.0], + "EXT_MCLDFFRL_KF_T": [1.3, 2.0, 0.9], + "EXT_MCLDFFRL_KF_TT": [0.0, 0.0, 0.0], + "EXT_MCLDFFRL_KF_TTT": [0.0, 0.0, 0.0], + "EXT_MCLDFFRL_EXP": [0.0, 0.0, 0.0], + "EXT_MCLDFFRL_EXP_T": [0.8, 3.0, 3.5], + "EXT_MCLDFFRL_EXP_TT": [0.0, 0.0, 0.0], + "EXT_MCLDFFRL_EXP_TTT": [0.0, 0.0, 0.0], + "EXT_MCLDFFRL_A": [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + "EXT_MCLDFFRL_A_T": [2.2, 1.1, 0.4, 0.1, 0.94, 2.8, 0.5, 1.2, 2.4], + "EXT_MCLDFFRL_A_TT": [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + "EXT_MCLDFFRL_A_TTT": [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], + "EXT_MCLDFFRL_TAU": 0.1 + )json", \ + 1e-10, 1e-10, CADET_NONBINDING_LIQUIDPHASE_COMP_USED, CADET_DONT_COMPARE_BINDING_VS_NONBINDING) + CADET_BINDINGTEST("MULTI_COMPONENT_SPREADING", "EXT_MULTI_COMPONENT_SPREADING", (2,2), (2,0,2), (1.2, 1.5, 0.1, 0.2, 0.3, 0.4), (1.2, 0.5, 1.5, 0.1, 0.2, 0.3, 0.4), \ R"json( "MCSPR_KA": [1.14, 2.0, 1.5, 1.9], "MCSPR_KD": [0.004, 0.008, 0.006, 0.002],