diff --git a/src/libcadet/model/binding/LangmuirLDFBinding.cpp b/src/libcadet/model/binding/LangmuirLDFBinding.cpp index 3a48a3f23..7c243e622 100644 --- a/src/libcadet/model/binding/LangmuirLDFBinding.cpp +++ b/src/libcadet/model/binding/LangmuirLDFBinding.cpp @@ -1,7 +1,7 @@ // ============================================================================= // CADET // -// Copyright © 2008-2022: The CADET Authors +// Copyright © 2008-2024: The CADET Authors // Please see the AUTHORS and CONTRIBUTORS file. // // All rights reserved. This program and the accompanying materials @@ -18,6 +18,7 @@ #include "LocalVector.hpp" #include "SimulationTypes.hpp" +#include #include #include #include @@ -31,238 +32,250 @@ [ { "type": "ScalarComponentDependentParameter", "varName": "keq", "confName": "MCLLDF_KEQ"}, { "type": "ScalarComponentDependentParameter", "varName": "kkin", "confName": "MCLLDF_KKIN"}, - { "type": "ScalarComponentDependentParameter", "varName": "qMax", "confName": "MCLLDF_QMAX"} + { "type": "ScalarComponentDependentParameter", "varName": "qMax", "confName": "MCLLDF_QMAX"}, + { "type": "ScalarComponentDependentParameter", "varName": "SolSA", "confName": "MCLLDF_SOLSA"}, + { "type": "ScalarComponentDependentParameter", "varName": "SolSB", "confName": "MCLLDF_SOLSB"} ] } */ + /* Parameter description ------------------------ keq = Equillibrium constant kkin = Linear driving force qMax = Capacity + solsa = solvent strength numerator + solsb = solvent strength denominator + + q = keq * exp(-solsa * phi) * C_p/(1 + keq/qmax * exp(-solsb * phi) Cp) */ namespace cadet { -namespace model -{ - -inline const char* LangmuirLDFParamHandler::identifier() CADET_NOEXCEPT { return "MULTI_COMPONENT_LANGMUIR_LDF"; } - -inline bool LangmuirLDFParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) -{ - if ((_keq.size() != _kkin.size()) || (_keq.size() != _qMax.size()) || (_keq.size() < nComp)) - throw InvalidParameterException("MCLLDF_KEQ, MCLLDF_KKIN, and MCLLDF_QMAX have to have the same size"); - - return true; -} - -inline const char* ExtLangmuirLDFParamHandler::identifier() CADET_NOEXCEPT { return "EXT_MULTI_COMPONENT_LANGMUIR_LDF"; } - -inline bool ExtLangmuirLDFParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) -{ - if ((_keq.size() != _kkin.size()) || (_keq.size() != _qMax.size()) || (_keq.size() < nComp)) - throw InvalidParameterException("EXT_MCLLDF_KEQ, EXT_MCLLDF_KKIN, and EXT_MCLLDF_QMAX have to have the same size"); - - return true; -} - - -/** - * @brief Defines the multi component Langmuir binding model - * @details Implements the Langmuir adsorption model: \f[ \begin{align} - * \frac{\mathrm{d}q_i}{\mathrm{d}t} &= k_{a,i} c_{p,i} q_{\text{max},i} \left( 1 - \sum_j \frac{q_j}{q_{\text{max},j}} \right) - k_{d,i} q_i - * \end{align} \f] - * Multiple bound states are not supported. - * Components without bound state (i.e., non-binding components) are supported. - * - * See @cite Langmuir1916. - * @tparam ParamHandler_t Type that can add support for external function dependence - */ -template -class LangmuirLDFBindingBase : public ParamHandlerBindingModelBase -{ -public: - - LangmuirLDFBindingBase() { } - virtual ~LangmuirLDFBindingBase() CADET_NOEXCEPT { } - - static const char* identifier() { return ParamHandler_t::identifier(); } - - virtual bool implementsAnalyticJacobian() const CADET_NOEXCEPT { return true; } - - /*virtual void timeDerivativeQuasiStationaryFluxes(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* yCp, double const* y, double* dResDt, LinearBufferAllocator workSpace) const + namespace model { - if (!this->hasQuasiStationaryReactions()) - return; - - if (!ParamHandler_t::dependsOnTime()) - return; - - // Update external function and compute time derivative of parameters - typename ParamHandler_t::ParamsHandle p; - typename ParamHandler_t::ParamsHandle dpDt; - std::tie(p, dpDt) = _paramHandler.updateTimeDerivative(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); - - // Protein fluxes: k_{kin,i}q_i - k_{kin,i} \frac{q_{m,i}k_{eq,i}c_i}{1+\sum_{j=1}^{n_{comp}} k_{eq,j}c_j} - double cpSum = 1.0; - double cpSumT = 0.0; - 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; - const double summand = y[bndIdx] / static_cast(p->qMax[i]); - cpSum -= summand; - cpSumT += summand / static_cast(p->qMax[i]) * static_cast(dpDt->qMax[i]); + inline const char* LangmuirLDFParamHandler::identifier() CADET_NOEXCEPT { return "MULTI_COMPONENT_LANGMUIR_LDF"; } - // Next bound component - ++bndIdx; - } - - bndIdx = 0; - for (int i = 0; i < _nComp; ++i) + inline bool LangmuirLDFParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) { - // Skip components without bound states (bound state index bndIdx is not advanced) - if (_nBoundStates[i] == 0) - continue; - - // Residual - dResDt[bndIdx] = static_cast(dpDt->kD[i]) * y[bndIdx] - - yCp[i] * (static_cast(dpDt->kA[i]) * static_cast(p->qMax[i]) * cpSum - + static_cast(p->kA[i]) * static_cast(dpDt->qMax[i]) * cpSum - + static_cast(p->kA[i]) * static_cast(p->qMax[i]) * cpSumT); - - // Next bound component - ++bndIdx; - } - }*/ + if ((_keq.size() != _kkin.size()) || (_keq.size() != _qMax.size()) || (_keq.size() < nComp)) + throw InvalidParameterException("MCLLDF_KEQ, MCLLDF_KKIN, and MCLLDF_QMAX have to have the same size"); - CADET_BINDINGMODELBASE_BOILERPLATE - -protected: - using ParamHandlerBindingModelBase::_paramHandler; - using ParamHandlerBindingModelBase::_reactionQuasistationarity; - using ParamHandlerBindingModelBase::_nComp; - using ParamHandlerBindingModelBase::_nBoundStates; + return true; + } - template - int fluxImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, StateType const* y, - CpStateType const* yCp, ResidualType* res, LinearBufferAllocator workSpace) const - { - typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); + inline const char* ExtLangmuirLDFParamHandler::identifier() CADET_NOEXCEPT { return "EXT_MULTI_COMPONENT_LANGMUIR_LDF"; } - // Protein fluxes: k_{kin,i}q_i - k_{kin,i} \frac{q_{m,i}k_{eq,i}c_i}{1+\sum_{j=1}^{n_{comp}} k_{eq,j}c_j} - ResidualType cpSum = 1.0; - unsigned int bndIdx = 0; - for (int i = 0; i < _nComp; ++i) + inline bool ExtLangmuirLDFParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) { - // Skip components without bound states (bound state index bndIdx is not advanced) - if (_nBoundStates[i] == 0) - continue; + if ((_keq.size() != _kkin.size()) || (_keq.size() != _qMax.size()) || (_keq.size() < nComp)) + throw InvalidParameterException("EXT_MCLLDF_KEQ, EXT_MCLLDF_KKIN, and EXT_MCLLDF_QMAX have to have the same size"); - cpSum += yCp[i] * static_cast(p->keq[i]); - - // Next bound component - ++bndIdx; + return true; } - 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; - // Residual - res[bndIdx] = static_cast(p->kkin[i]) * (y[bndIdx] - static_cast(p->keq[i]) * yCp[i] * static_cast(p->qMax[i]) / cpSum); + /** + * @brief Defines the multi component Langmuir binding model + * @details Implements the Langmuir adsorption model: \f[ \begin{align} + * \frac{\mathrm{d}q_i}{\mathrm{d}t} &= k_{a,i} c_{p,i} q_{\text{max},i} \left( 1 - \sum_j \frac{q_j}{q_{\text{max},j}} \right) - k_{d,i} q_i + * \end{align} \f] + * Multiple bound states are not supported. + * Components without bound state (i.e., non-binding components) are supported. + * + * See @cite Langmuir1916. + * @tparam ParamHandler_t Type that can add support for external function dependence + */ + template + class LangmuirLDFBindingBase : public ParamHandlerBindingModelBase + { + public: - // Next bound component - ++bndIdx; - } + LangmuirLDFBindingBase() { } + virtual ~LangmuirLDFBindingBase() CADET_NOEXCEPT { } - return 0; - } + static const char* identifier() { return ParamHandler_t::identifier(); } - 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); + virtual bool implementsAnalyticJacobian() const CADET_NOEXCEPT { return true; } - // Protein fluxes: k_{kin,i}q_i - k_{kin,i} \frac{q_{m,i}k_{eq,i}c_i}{1+\sum_{j=1}^{n_{comp}} k_{eq,j}c_j} - double cpSum = 1.0; - 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; + /*virtual void timeDerivativeQuasiStationaryFluxes(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* yCp, double const* y, double* dResDt, LinearBufferAllocator workSpace) const + { + if (!this->hasQuasiStationaryReactions()) + return; + + if (!ParamHandler_t::dependsOnTime()) + return; + + // Update external function and compute time derivative of parameters + typename ParamHandler_t::ParamsHandle p; + typename ParamHandler_t::ParamsHandle dpDt; + std::tie(p, dpDt) = _paramHandler.updateTimeDerivative(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); + + // Protein fluxes: k_{kin,i}q_i - k_{kin,i} \frac{q_{m,i}k_{eq,i}c_i}{1+\sum_{j=1}^{n_{comp}} k_{eq,j}c_j} + double cpSum = 1.0; + double cpSumT = 0.0; + 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; + + const double summand = y[bndIdx] / static_cast(p->qMax[i]); + cpSum -= summand; + cpSumT += summand / static_cast(p->qMax[i]) * static_cast(dpDt->qMax[i]); + + // Next bound component + ++bndIdx; + } + + 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; + + // Residual + dResDt[bndIdx] = static_cast(dpDt->kD[i]) * y[bndIdx] + - yCp[i] * (static_cast(dpDt->kA[i]) * static_cast(p->qMax[i]) * cpSum + + static_cast(p->kA[i]) * static_cast(dpDt->qMax[i]) * cpSum + + static_cast(p->kA[i]) * static_cast(p->qMax[i]) * cpSumT); + + // Next bound component + ++bndIdx; + } + }*/ + + 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 + { + typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); + + // Protein fluxes: k_{kin,i}q_i - k_{kin,i} \frac{q_{m,i}k_{eq,i}c_i}{1+\sum_{j=1}^{n_{comp}} k_{eq,j}c_j} + ResidualType cpSum = 1.0; + 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; + + cpSum += yCp[i] * static_cast(p->keq[i]) * exp(-static_cast(p->SolSB[i]) * yCp[0]) / static_cast(p->qMax[i]); + + // Next bound component + ++bndIdx; + } + + 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; + + // Residual + res[bndIdx] = static_cast(p->kkin[i]) * (y[bndIdx] - static_cast(p->keq[i]) * exp(-static_cast(p->SolSA[i]) * yCp[0]) * yCp[i] / cpSum); + + // Next bound component + ++bndIdx; + } + + return 0; + } - cpSum += yCp[i] * static_cast(p->keq[i]); + 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); + + // Protein fluxes: k_{kin,i}q_i - k_{kin,i} \frac{q_{m,i}k_{eq,i}c_i}{1+\sum_{j=1}^{n_{comp}} k_{eq,j}c_j} + double cpSum = 1.0; + 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; + + const double keq = static_cast(p->keq[i]); + const double qMax = static_cast(p->qMax[i]); + const double kkin = static_cast(p->kkin[i]); + const double solsa = static_cast(p->SolSA[i]); + const double solsb = static_cast(p->SolSB[i]); + + cpSum += yCp[i] * static_cast(p->keq[i]) * exp(-static_cast(p->SolSB[i]) * yCp[0]) / static_cast(p->qMax[i]); + + // Next bound component + ++bndIdx; + } + + 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; + + //(keq *keq *static_cast(p->qMax[i]) * yCp[i] / (cpSum*cpSum)) + // dres_i / dc_{p,i} + jac[i - bndIdx - offsetCp] = -static_cast(p->kkin[i]) * static_cast(p->keq[i]) * exp(-static_cast(p->SolSA[i]) * yCp[0]) / cpSum; + // 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}. + + // Fill dres_i / dc_j + const double commonFactor = static_cast(p->keq[i]) * static_cast(p->kkin[i]) * yCp[i] * exp(-static_cast(p->SolSA[i]) * yCp[0]) / (cpSum * cpSum); + + int bndIdx2 = 0; + for (int j = 0; j < _nComp; ++j) + { + // Skip components without bound states (bound state index bndIdx is not advanced) + if (_nBoundStates[j] == 0) + continue; + + // dres_i / dc_j + jac[j - bndIdx - offsetCp] += static_cast(p->keq[j]) * exp(-static_cast(p->SolSB[j]) * yCp[0]) / static_cast(p->qMax[j]) * commonFactor; + // 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}. + + + ++bndIdx2; + } + + // Add to dres_i / dq_i + jac[0] += static_cast(p->kkin[i]); + + // Advance to next flux and Jacobian row + ++bndIdx; + ++jac; + } + } + }; - // Next bound component - ++bndIdx; - } + typedef LangmuirLDFBindingBase LangmuirLDFBinding; + typedef LangmuirLDFBindingBase ExternalLangmuirLDFBinding; - bndIdx = 0; - for (int i = 0; i < _nComp; ++i) + namespace binding { - // Skip components without bound states (bound state index bndIdx is not advanced) - if (_nBoundStates[i] == 0) - continue; - - const double keq = static_cast(p->keq[i]); - const double kkin = static_cast(p->kkin[i]); - //(keq *keq *static_cast(p->qMax[i]) * yCp[i] / (cpSum*cpSum)) - // dres_i / dc_{p,i} - jac[i - bndIdx - offsetCp] = -kkin * keq * static_cast(p->qMax[i]) / cpSum; - // 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}. - - // Fill dres_i / dc_j - const double commonFactor = keq * kkin * static_cast(p->qMax[i]) * yCp[i] / (cpSum * cpSum); - int bndIdx2 = 0; - for (int j = 0; j < _nComp; ++j) + void registerLangmuirLDFModel(std::unordered_map>& bindings) { - // Skip components without bound states (bound state index bndIdx is not advanced) - if (_nBoundStates[j] == 0) - continue; - - // dres_i / dc_j - jac[j - bndIdx - offsetCp] += static_cast(p->keq[j]) * commonFactor; - // 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}. - - - ++bndIdx2; + bindings[LangmuirLDFBinding::identifier()] = []() { return new LangmuirLDFBinding(); }; + bindings[ExternalLangmuirLDFBinding::identifier()] = []() { return new ExternalLangmuirLDFBinding(); }; } + } // namespace binding - // Add to dres_i / dq_i - jac[0] += kkin; - - // Advance to next flux and Jacobian row - ++bndIdx; - ++jac; - } - } -}; - -typedef LangmuirLDFBindingBase LangmuirLDFBinding; -typedef LangmuirLDFBindingBase ExternalLangmuirLDFBinding; - -namespace binding -{ - void registerLangmuirLDFModel(std::unordered_map>& bindings) - { - bindings[LangmuirLDFBinding::identifier()] = []() { return new LangmuirLDFBinding(); }; - bindings[ExternalLangmuirLDFBinding::identifier()] = []() { return new ExternalLangmuirLDFBinding(); }; - } -} // namespace binding - -} // namespace model + } // namespace model } // namespace cadet