diff --git a/src/libcadet/BindingModelFactory.cpp b/src/libcadet/BindingModelFactory.cpp index ecfe15bf4..6091e4b9a 100644 --- a/src/libcadet/BindingModelFactory.cpp +++ b/src/libcadet/BindingModelFactory.cpp @@ -38,6 +38,7 @@ namespace cadet void registerMultiComponentSpreadingModel(std::unordered_map>& bindings); void registerGeneralizedIonExchangeModel(std::unordered_map>& bindings); void registerColloidalModel(std::unordered_map>& bindings); + void registerConstantWaterActivityModel(std::unordered_map>& bindings); } } @@ -61,6 +62,7 @@ namespace cadet model::binding::registerMultiComponentSpreadingModel(_bindingModels); model::binding::registerGeneralizedIonExchangeModel(_bindingModels); model::binding::registerColloidalModel(_bindingModels); + model::binding::registerConstantWaterActivityModel(_bindingModels); registerModel(); } diff --git a/src/libcadet/CMakeLists.txt b/src/libcadet/CMakeLists.txt index 9615bafaf..eda2e100c 100644 --- a/src/libcadet/CMakeLists.txt +++ b/src/libcadet/CMakeLists.txt @@ -94,6 +94,7 @@ set(LIBCADET_BINDINGMODEL_SOURCES ${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/MultiComponentSpreadingBinding.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/GeneralizedIonExchangeBinding.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/ColloidalBinding.cpp + ${CMAKE_SOURCE_DIR}/src/libcadet/model/binding/ConstantWaterActivityBinding.cpp ) # LIBCADET_REACTIONMODEL_SOURCES holds all source files of reaction models diff --git a/src/libcadet/model/binding/ConstantWaterActivityBinding.cpp b/src/libcadet/model/binding/ConstantWaterActivityBinding.cpp new file mode 100644 index 000000000..14e1a884e --- /dev/null +++ b/src/libcadet/model/binding/ConstantWaterActivityBinding.cpp @@ -0,0 +1,287 @@ +// ============================================================================= +// CADET +// +// Copyright © 2008-2021: 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/binding/BindingModelMacros.hpp" +#include "model/binding/RefConcentrationSupport.hpp" +#include "model/ModelUtils.hpp" +#include "cadet/Exceptions.hpp" +#include "model/Parameters.hpp" +#include "LocalVector.hpp" +#include "SimulationTypes.hpp" + +#include +#include +#include +#include + +/* +{ + "name": "CWAParamHandler", + "externalName": "ExtCWAParamHandler", + "parameters": + [ + { "type": "ScalarParameter", "varName": "beta0", "confName": "BETA0"}, + { "type": "ScalarParameter", "varName": "beta1", "confName": "BETA1"}, + { "type": "ScalarComponentDependentParameter", "varName": "kA", "confName": "KA"}, + { "type": "ScalarComponentDependentParameter", "varName": "kD", "confName": "KD"}, + { "type": "ScalarComponentDependentParameter", "varName": "nu", "confName": "NU"}, + { "type": "ScalarComponentDependentParameter", "varName": "qMax", "confName": "QMAX"} + ] +} +*/ + +/* Parameter description + ------------------------ + beta0 = bulk-like ordered water molecules at infinitely diluted salt concentration + beta1 = influence of salt concentration on bulk-like ordered water molecules + kA = Adsorption constant + kD = Desorption constant + nu = Number of binding sites + qMax = Maximum binding capacity +*/ + +namespace cadet +{ + + namespace model + { + + inline const char* CWAParamHandler::identifier() CADET_NOEXCEPT { return "CONSTANT_WATER_ACTIVITY"; } + + inline bool CWAParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) + { + if ((_kA.size() != _kD.size()) || (_kA.size() != _nu.size()) || (_kA.size() != _qMax.size()) || (_kA.size() < nComp)) + throw InvalidParameterException("KA, KD, NU, and QMAX have to have the same size"); + + return true; + } + + inline const char* ExtCWAParamHandler::identifier() CADET_NOEXCEPT { return "EXT_CONSTANT_WATER_ACTIVITY"; } + + inline bool ExtCWAParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) + { + if ((_kA.size() != _kD.size()) || (_kA.size() != _nu.size()) || (_kA.size() != _qMax.size()) || (_kA.size() < nComp)) + throw InvalidParameterException("KA, KD, NU, and QMAX have to have the same size"); + + return true; + } + + + /** + * @brief Defines the HIC Isotherm assuming a constant water activity as described by Jäpel and Buyel 2022 + * @details Implements the the HIC Isotherm assuming a constant water activity: \f[ \begin{align} + * \beta&=\beta_0 e^{c_{p,0}\beta_1}\\ + * \frac{\mathrm{d}q_i}{\mathrm{d}t} &= k_{a,i} c_{p,i} \left( 1 - \sum_j \frac{q_j}{q_{max,j}} \right)^{\nu_i} - k_{d,i} q_i 0.1^{\nu_i \beta} + * \end{align} \f] + * Component @c 0 is assumed to be salt without a bound state. Multiple bound states are not supported. + * Components without bound state (i.e., salt and non-binding components) are supported. + * + * See @cite Jäpel and Buyel 2022. + * @tparam ParamHandler_t Type that can add support for external function dependence + */ + template + class ConstantWaterActivityBindingBase : public ParamHandlerBindingModelBase + { + public: + + ConstantWaterActivityBindingBase() { } + virtual ~ConstantWaterActivityBindingBase() CADET_NOEXCEPT { } + + static const char* identifier() { return ParamHandler_t::identifier(); } + + virtual bool configureModelDiscretization(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBound, unsigned int const* boundOffset) + { + const bool res = BindingModelBase::configureModelDiscretization(paramProvider, nComp, nBound, boundOffset); + + // Guarantee that salt has no bound state + if (nBound[0] != 0) + throw InvalidParameterException("CWA binding model requires exactly zero bound states for salt component"); + + // First flux is salt, which is always quasi-stationary + _reactionQuasistationarity[0] = false; + + return res; + } + + virtual bool hasSalt() const CADET_NOEXCEPT { return true; } + virtual bool supportsMultistate() const CADET_NOEXCEPT { return false; } + virtual bool supportsNonBinding() const CADET_NOEXCEPT { return true; } + virtual bool hasQuasiStationaryReactions() const CADET_NOEXCEPT { return false; } + virtual bool implementsAnalyticJacobian() const CADET_NOEXCEPT { return true; } + + virtual bool preConsistentInitialState(double t, unsigned int secIdx, const ColumnPosition& colPos, double* y, double const* yCp, LinearBufferAllocator workSpace) const + { + typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); + return true; + } + + virtual void postConsistentInitialState(double t, unsigned int secIdx, const ColumnPosition& colPos, double* y, double const* yCp, LinearBufferAllocator workSpace) const + { + preConsistentInitialState(t, secIdx, colPos, y, yCp, workSpace); + } + + + 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 CpStateParamTypeParamType = typename DoubleActivePromoter::type; + using StateParamType = typename DoubleActivePromoter::type; + + typename ParamHandler_t::ParamsHandle const _p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); + + + auto beta0 = static_cast(_p->beta0); + auto beta1 = static_cast(_p->beta1); + auto water_activity = static_cast(0.1); + + auto beta = static_cast(beta0 * exp(beta1 * yCp[0])); + + ResidualType qSum = 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; + + qSum -= y[bndIdx] / 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; + + auto kD = static_cast(_p->kD[i]); + auto kA = static_cast(_p->kA[i]); + auto nu = static_cast(_p->nu[i]); + + ResidualType bulk_water = pow(water_activity, static_cast(_p->nu[i]) * beta); + + res[bndIdx] = kD * y[bndIdx] * static_cast(bulk_water) - kA * static_cast(pow(qSum, nu)) * yCp[i]; + + // 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); + + auto beta0 = static_cast(_p->beta0); + auto beta1 = static_cast(_p->beta1); + auto water_activity = static_cast(0.1); + + auto beta = static_cast(beta0 * exp(beta1 * yCp[0])); + + double qSum = 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; + + qSum -= y[bndIdx] / 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; + + const double kD = static_cast(_p->kD[i]); + const double kA = static_cast(_p->kA[i]); + const double nu = static_cast(_p->nu[i]); + + const double bulk_water = pow(water_activity, nu * beta); + + // dres_i / dc_{p,0} + jac[-bndIdx - offsetCp] = kD * y[bndIdx] * static_cast(bulk_water) * static_cast(std::log(water_activity)) * nu * beta * beta1; + + + // dres_i / dc_{p,i} + jac[i - bndIdx - offsetCp] = -kA * static_cast(pow(qSum, nu)); + // 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 / dq_j + int bndIdx2 = 0; + for (int j = 1; j < _nComp; ++j) + { + // Skip components without bound states (bound state index bndIdx is not advanced) + if (_nBoundStates[j] == 0) + continue; + + // dres_i / dq_j + jac[bndIdx2 - bndIdx] = -kA * yCp[i] * nu * static_cast(pow(qSum, nu - 1)) / (-static_cast(_p->qMax[j])); + // Getting to q_j: -bndIdx takes us to q_0, another +bndIdx2 to q_j. This means jac[bndIdx2 - bndIdx] corresponds to q_j. + + ++bndIdx2; + } + + // Add to dres_i / dq_i + jac[0] += kD * bulk_water; + + // Advance to next flux and Jacobian row + ++bndIdx; + ++jac; + } + } + }; + + + typedef ConstantWaterActivityBindingBase ConstantWaterActivityBinding; + typedef ConstantWaterActivityBindingBase ExternalConstantWaterActivityBinding; + + namespace binding + { + void registerConstantWaterActivityModel(std::unordered_map>& bindings) + { + bindings[ConstantWaterActivityBinding::identifier()] = []() { return new ConstantWaterActivityBinding(); }; + bindings[ExternalConstantWaterActivityBinding::identifier()] = []() { return new ExternalConstantWaterActivityBinding(); }; + } + } // namespace binding + + } // namespace model + +} // namespace cadet