From fc8a3cefdbe2ee9e0859a2e4e790340ffbcaed01 Mon Sep 17 00:00:00 2001 From: unknown Date: Fri, 23 Feb 2024 17:42:51 +0000 Subject: [PATCH 1/3] Changed file LangmuirLDFBinding --- .../model/binding/LangmuirLDFBinding.cpp | 52 +++++++++++++------ 1 file changed, 37 insertions(+), 15 deletions(-) diff --git a/src/libcadet/model/binding/LangmuirLDFBinding.cpp b/src/libcadet/model/binding/LangmuirLDFBinding.cpp index 3a48a3f23..ce639a6a7 100644 --- a/src/libcadet/model/binding/LangmuirLDFBinding.cpp +++ b/src/libcadet/model/binding/LangmuirLDFBinding.cpp @@ -29,20 +29,28 @@ "externalName": "ExtLangmuirLDFParamHandler", "parameters": [ - { "type": "ScalarComponentDependentParameter", "varName": "keq", "confName": "MCLLDF_KEQ"}, - { "type": "ScalarComponentDependentParameter", "varName": "kkin", "confName": "MCLLDF_KKIN"}, - { "type": "ScalarComponentDependentParameter", "varName": "qMax", "confName": "MCLLDF_QMAX"} + { "type": "ScalarComponentDependentParameter", "varName": "keq", "confName": "KEQ"}, + { "type": "ScalarComponentDependentParameter", "varName": "kkin", "confName": "KKIN"}, + { "type": "ScalarComponentDependentParameter", "varName": "qMax", "confName": "QMAX"}, + { "type": "ScalarComponentDependentParameter", "varName": "SolSA", "confName": "SOLSA"}, + { "type": "ScalarComponentDependentParameter", "varName": "SolSB", "confName": "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 { @@ -165,7 +173,7 @@ class LangmuirLDFBindingBase : public ParamHandlerBindingModelBase(p->keq[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; @@ -179,7 +187,7 @@ class LangmuirLDFBindingBase : public ParamHandlerBindingModelBase(p->kkin[i]) * (y[bndIdx] - static_cast(p->keq[i]) * yCp[i] * static_cast(p->qMax[i]) / cpSum); + 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; @@ -194,7 +202,7 @@ class LangmuirLDFBindingBase : public ParamHandlerBindingModelBase(p->keq[i]); + 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] * keq + qMax * exp(solsb * yCp[0]); // Next bound component ++bndIdx; @@ -216,39 +230,47 @@ class LangmuirLDFBindingBase : public ParamHandlerBindingModelBase(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]); + //(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; + jac[i - bndIdx - offsetCp] = -kkin * keq * qMax * exp((-solsa + solsb) * 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 = keq * kkin * static_cast(p->qMax[i]) * yCp[i] / (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; - + const double keq = static_cast(p->keq[j]); + const double qMax = static_cast(p->qMax[j]); + const double kkin = static_cast(p->kkin[j]); + const double solsa = static_cast(p->SolSA[j]); + const double solsb = static_cast(p->SolSB[j]); // dres_i / dc_j - jac[j - bndIdx - offsetCp] += static_cast(p->keq[j]) * commonFactor; + jac[j - bndIdx - offsetCp] += pow(keq, 2.0) * kkin * qMax * yCp[j] * exp((-solsa + solsb) * yCp[0]) / pow(cpSum, 2.0); // 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] += kkin; + jac[0] += static_cast(p->kkin[i]); // Advance to next flux and Jacobian row ++bndIdx; ++jac; } - } + } }; typedef LangmuirLDFBindingBase LangmuirLDFBinding; @@ -256,7 +278,7 @@ typedef LangmuirLDFBindingBase ExternalLangmuirLDFBi namespace binding { - void registerLangmuirLDFModel(std::unordered_map>& bindings) + void registerLangmuirLDFModel(std::unordered_map>& bindings) { bindings[LangmuirLDFBinding::identifier()] = []() { return new LangmuirLDFBinding(); }; bindings[ExternalLangmuirLDFBinding::identifier()] = []() { return new ExternalLangmuirLDFBinding(); }; From b0c85a579406769e56a13974317f54d179287047 Mon Sep 17 00:00:00 2001 From: unknown Date: Mon, 8 Apr 2024 18:18:44 +0100 Subject: [PATCH 2/3] Fixed analytic jacobian and parameter prefixes --- .../model/binding/LangmuirLDFBinding.cpp | 423 +++++++++--------- 1 file changed, 207 insertions(+), 216 deletions(-) diff --git a/src/libcadet/model/binding/LangmuirLDFBinding.cpp b/src/libcadet/model/binding/LangmuirLDFBinding.cpp index ce639a6a7..9ba276c49 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 @@ -29,11 +30,11 @@ "externalName": "ExtLangmuirLDFParamHandler", "parameters": [ - { "type": "ScalarComponentDependentParameter", "varName": "keq", "confName": "KEQ"}, - { "type": "ScalarComponentDependentParameter", "varName": "kkin", "confName": "KKIN"}, - { "type": "ScalarComponentDependentParameter", "varName": "qMax", "confName": "QMAX"}, - { "type": "ScalarComponentDependentParameter", "varName": "SolSA", "confName": "SOLSA"}, - { "type": "ScalarComponentDependentParameter", "varName": "SolSB", "confName": "SOLSB"} + { "type": "ScalarComponentDependentParameter", "varName": "keq", "confName": "MCLLDF_KEQ"}, + { "type": "ScalarComponentDependentParameter", "varName": "kkin", "confName": "MCLLDF_KKIN"}, + { "type": "ScalarComponentDependentParameter", "varName": "qMax", "confName": "MCLLDF_QMAX"}, + { "type": "ScalarComponentDependentParameter", "varName": "SolSA", "confName": "MCLLDF_SOLSA"}, + { "type": "ScalarComponentDependentParameter", "varName": "SolSB", "confName": "MCLLDF_SOLSB"} ] } */ @@ -50,241 +51,231 @@ 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]); - - // 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); + inline const char* LangmuirLDFParamHandler::identifier() CADET_NOEXCEPT { return "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 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; + 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"); - 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; + 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; + inline const char* ExtLangmuirLDFParamHandler::identifier() CADET_NOEXCEPT { return "EXT_MULTI_COMPONENT_LANGMUIR_LDF"; } - // 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); + 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"); - // Next bound component - ++bndIdx; + return true; } - 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); - - // 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 = 0.0; - int bndIdx = 0; - for (int i = 0; i < _nComp; ++i) + /** + * @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 { - // Skip components without bound states (bound state index bndIdx is not advanced) - if (_nBoundStates[i] == 0) - continue; + public: - 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]); + LangmuirLDFBindingBase() { } + virtual ~LangmuirLDFBindingBase() CADET_NOEXCEPT { } - cpSum += yCp[i] * keq + qMax * exp(solsb * yCp[0]); + static const char* identifier() { return ParamHandler_t::identifier(); } - // Next bound component - ++bndIdx; - } + virtual bool implementsAnalyticJacobian() const CADET_NOEXCEPT { 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; - - 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]); - - //(keq *keq *static_cast(p->qMax[i]) * yCp[i] / (cpSum*cpSum)) - // dres_i / dc_{p,i} - jac[i - bndIdx - offsetCp] = -kkin * keq * qMax * exp((-solsa + solsb) * 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 - - int bndIdx2 = 0; - for (int j = 0; j < _nComp; ++j) + /*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 { - // Skip components without bound states (bound state index bndIdx is not advanced) - if (_nBoundStates[j] == 0) - continue; - const double keq = static_cast(p->keq[j]); - const double qMax = static_cast(p->qMax[j]); - const double kkin = static_cast(p->kkin[j]); - const double solsa = static_cast(p->SolSA[j]); - const double solsb = static_cast(p->SolSB[j]); - // dres_i / dc_j - jac[j - bndIdx - offsetCp] += pow(keq, 2.0) * kkin * qMax * yCp[j] * exp((-solsa + solsb) * yCp[0]) / pow(cpSum, 2.0); - // 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; + 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; } - // Add to dres_i / dq_i - jac[0] += static_cast(p->kkin[i]); - - // Advance to next flux and Jacobian row - ++bndIdx; - ++jac; - } - } -}; + 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; + } + } + }; -typedef LangmuirLDFBindingBase LangmuirLDFBinding; -typedef LangmuirLDFBindingBase ExternalLangmuirLDFBinding; + 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 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 From f86ea581e8802df1003127358e575379926f7b15 Mon Sep 17 00:00:00 2001 From: Kats1247 <151026874+Kats1247@users.noreply.github.com> Date: Fri, 17 May 2024 11:29:10 +0000 Subject: [PATCH 3/3] Fixed parenthesis --- src/libcadet/model/binding/LangmuirLDFBinding.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/libcadet/model/binding/LangmuirLDFBinding.cpp b/src/libcadet/model/binding/LangmuirLDFBinding.cpp index 9ba276c49..7c243e622 100644 --- a/src/libcadet/model/binding/LangmuirLDFBinding.cpp +++ b/src/libcadet/model/binding/LangmuirLDFBinding.cpp @@ -187,7 +187,7 @@ namespace cadet 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); + 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;