From eb9fde17e34a76eca6553c983127c581cc60d2ea Mon Sep 17 00:00:00 2001 From: Samuel Leweke Date: Wed, 21 Jul 2021 01:42:44 +0200 Subject: [PATCH 1/8] Add vector valued parameters Adds a vector valued and a vector valued component dependent parameter class. The ordering of the latter is component-major (i.e., component index changes the slowest). The reaction index is used for the vector element index in the ParameterId. --- src/libcadet/model/Parameters.hpp | 329 ++++++++++++++++++++++++++++++ 1 file changed, 329 insertions(+) diff --git a/src/libcadet/model/Parameters.hpp b/src/libcadet/model/Parameters.hpp index 808075c5c..ebc69d466 100644 --- a/src/libcadet/model/Parameters.hpp +++ b/src/libcadet/model/Parameters.hpp @@ -186,6 +186,42 @@ class ScalarParameter }; +/** + * @brief Vector parameter + * @details A single value per component. + */ +class VectorParameter +{ +public: + + typedef std::vector storage_t; + + VectorParameter(std::vector& p) : _p(&p) { } + VectorParameter(std::vector* p) : _p(p) { } + + inline void configure(const std::string& varName, IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBoundStates) + { + readScalarParameterOrArray(*_p, paramProvider, varName, 1); + } + + inline void registerParam(const std::string& varName, std::unordered_map& parameters, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, unsigned int nComp, unsigned int const* nBoundStates) + { + const StringHash hashVar = hashStringRuntime(varName); + registerParam1DArray(parameters, *_p, [=](bool multi, unsigned int rid) { return makeParamId(hashVar, unitOpIdx, CompIndep, parTypeIdx, BoundStateIndep, rid, SectionIndep); }); + } + + inline void reserve(unsigned int numElem, unsigned int numSlices, unsigned int nComp, unsigned int const* nBoundStates) { } + + inline std::size_t size() const CADET_NOEXCEPT { return _p->size(); } + + inline const std::vector& get() const CADET_NOEXCEPT { return *_p; } + inline std::vector& get() CADET_NOEXCEPT { return *_p; } + +protected: + std::vector* _p; +}; + + /** * @brief Component dependent scalar external parameter * @details A single value per component. @@ -304,6 +340,42 @@ class ScalarReactionDependentParameter }; +/** + * @brief Component dependent vectorial external parameter + * @details A vector per component, the component index changes the slowest. + */ +class VectorComponentDependentParameter +{ +public: + + typedef std::vector storage_t; + + VectorComponentDependentParameter(std::vector& p) : _p(&p) { } + VectorComponentDependentParameter(std::vector* p) : _p(p) { } + + inline void configure(const std::string& varName, IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBoundStates) + { + readScalarParameterOrArray(*_p, paramProvider, varName, 1); + } + + inline void registerParam(const std::string& varName, std::unordered_map& parameters, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, unsigned int nComp, unsigned int const* nBoundStates) + { + const StringHash nameHash = hashStringRuntime(varName); + registerParam2DArray(parameters, *_p, [=](bool multi, unsigned int comp, unsigned int rid) { return makeParamId(nameHash, unitOpIdx, comp, parTypeIdx, BoundStateIndep, multi ? rid : ReactionIndep, SectionIndep); }, _p->size() / nComp); + } + + inline void reserve(unsigned int numElem, unsigned int numSlices, unsigned int nComp, unsigned int const* nBoundStates) { } + + inline std::size_t size() const CADET_NOEXCEPT { return _p->size(); } + + inline const std::vector& get() const CADET_NOEXCEPT { return *_p; } + inline std::vector& get() CADET_NOEXCEPT { return *_p; } + +protected: + std::vector* _p; +}; + + /** * @brief Component and bound state dependent external parameter * @details A single value per component and bound state / binding site type. @@ -749,6 +821,174 @@ class ExternalScalarParameter }; +/** + * @brief Vector external parameter + * @details A vector. + */ +class ExternalVectorParameter +{ +public: + + /** + * @brief Underlying type + */ + typedef std::vector storage_t; + + /** + * @brief Reads parameters and verifies them + * @details See IBindingModel::configure() for details. + * @param [in] varName Name of the parameter + * @param [in] paramProvider IParameterProvider used for reading parameters + * @param [in] nComp Number of components + * @param [in] nBoundStates Array with number of bound states for each component + * @return @c true if the parameters were read and validated successfully, otherwise @c false + */ + inline void configure(const std::string& varName, IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBoundStates) + { + readScalarParameterOrArray(_base, paramProvider, "EXT_" + varName, 1); + readScalarParameterOrArray(_linear, paramProvider, "EXT_" + varName + "_T", 1); + readScalarParameterOrArray(_quad, paramProvider, "EXT_" + varName + "_TT", 1); + readScalarParameterOrArray(_cube, paramProvider, "EXT_" + varName + "_TTT", 1); + } + + /** + * @brief Registers the parameters in a map for further use + * @param [in] varName Name of the parameter + * @param [in,out] parameters Map in which the parameters are stored + * @param [in] unitOpIdx Index of the unit operation used for registering the parameters + * @param [in] nComp Number of components + * @param [in] nBoundStates Array with number of bound states for each component + */ + inline void registerParam(const std::string& varName, std::unordered_map& parameters, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, unsigned int nComp, unsigned int const* nBoundStates) + { + const StringHash hashConst = hashStringRuntime("EXT_" + varName); + registerParam1DArray(parameters, _base, [=](bool multi, unsigned int rid) { return makeParamId(hashConst, unitOpIdx, CompIndep, parTypeIdx, BoundStateIndep, rid, SectionIndep); }); + + const StringHash hashLinear = hashStringRuntime("EXT_" + varName + "_T"); + registerParam1DArray(parameters, _linear, [=](bool multi, unsigned int rid) { return makeParamId(hashLinear, unitOpIdx, CompIndep, parTypeIdx, BoundStateIndep, rid, SectionIndep); }); + + const StringHash hashQuad = hashStringRuntime("EXT_" + varName + "_TT"); + registerParam1DArray(parameters, _quad, [=](bool multi, unsigned int rid) { return makeParamId(hashQuad, unitOpIdx, CompIndep, parTypeIdx, BoundStateIndep, rid, SectionIndep); }); + + const StringHash hashCube = hashStringRuntime("EXT_" + varName + "_TTT"); + registerParam1DArray(parameters, _cube, [=](bool multi, unsigned int rid) { return makeParamId(hashCube, unitOpIdx, CompIndep, parTypeIdx, BoundStateIndep, rid, SectionIndep); }); + } + + /** + * @brief Reserves space in the storage of the parameters + * @param [in] numElem Total number of components in all slices / binding site types + * @param [in] numSlices Number of slices / binding site types + * @param [in] nComp Number of components + * @param [in] nBoundStates Array with number of bound states for each component + */ + inline void reserve(unsigned int numElem, unsigned int numSlices, unsigned int nComp, unsigned int const* nBoundStates) { } + + /** + * @brief Calculates a parameter in order to take the external profile into account + * @param [out] result Stores the result of the paramter + * @param [in] extVal Value of the external function + * @param [in] nComp Number of components + * @param [in] nBoundStates Array with number of bound states for each component + */ + inline void update(std::vector& result, double extVal, unsigned int nComp, unsigned int const* nBoundStates) const + { + update(result.data(), extVal, nComp, nBoundStates); + } + + inline void update(active* result, double extVal, unsigned int nComp, unsigned int const* nBoundStates) const + { + for (std::size_t i = 0; i < _base.size(); ++i) + result[i] = _base[i] + extVal * (_linear[i] + extVal * (_quad[i] + extVal * _cube[i])); + } + + /** + * @brief Calculates time derivative of parameter in case of external dependence + * @param [out] result Stores the result of the paramter + * @param [in] extVal Value of the external function + * @param [in] extTimeDiff Time derivative of the external function + * @param [in] nComp Number of components + * @param [in] nBoundStates Array with number of bound states for each component + */ + inline void updateTimeDerivative(active& result, double extVal, double extTimeDiff, unsigned int nComp, unsigned int const* nBoundStates) const + { + updateTimeDerivative(&result, extVal, extTimeDiff, nComp, nBoundStates); + } + + /** + * @brief Calculates time derivative of parameter in case of external dependence + * @param [out] result Stores the result of the paramter + * @param [in] extVal Value of the external function + * @param [in] extTimeDiff Time derivative of the external function + * @param [in] nComp Number of components + * @param [in] nBoundStates Array with number of bound states for each component + */ + inline void updateTimeDerivative(std::vector& result, double extVal, double extTimeDiff, unsigned int nComp, unsigned int const* nBoundStates) const + { + updateTimeDerivative(result.data(), extVal, extTimeDiff, nComp, nBoundStates); + } + + inline void updateTimeDerivative(active* result, double extVal, double extTimeDiff, unsigned int nComp, unsigned int const* nBoundStates) const + { + for (std::size_t i = 0; i < _base.size(); ++i) + result[i] = extTimeDiff * (static_cast(_linear[i]) + extVal * (2.0 * static_cast(_quad[i]) + 3.0 * extVal * static_cast(_cube[i]))); + } + + /** + * @brief Returns the base value that does not depend on an external value + * @return Constant base value + */ + inline storage_t& base() CADET_NOEXCEPT { return _base; } + inline const storage_t& base() const CADET_NOEXCEPT { return _base; } + + /** + * @brief Returns the amount of additional memory (usually dynamically allocated by containers) for storing the final parameters + * @details In a model, externally dependent parameters are stored in a struct, usually called + * VariableParams. This is sufficient for "static" parameter types that do + * not require additional memory (which is usually allocated dynamically). + * For containers using additional dynamic memory, only the container itself + * is stored in the struct. Memory for the content of the container (i.e., the + * elements) is still required. This function computes this amount of additional + * memory. + * + * @param [in] nComp Number of components + * @param [in] totalNumBoundStates Total number of bound states + * @param [in] nBoundStates Array with number of bound states for each component + * @return Amount of additional memory in bytes + */ + inline std::size_t additionalDynamicMemory(unsigned int nComp, unsigned int totalNumBoundStates, unsigned int const* nBoundStates) const CADET_NOEXCEPT { return 0; } + + /** + * @brief Prepares the cache for the updated values + * @details The cache is a local version of storage_t (e.g., LocalVector). + * @param [in,out] cache Cache object to be prepared + * @param [in] ptr Pointer to cache buffer + */ + template + inline void prepareCache(T& cache, LinearBufferAllocator& buffer) const + { + cache.fromTemplate(buffer, _base); + } + + /** + * @brief Returns the number of elements in the parameter + * @return Number of elements in the parameter + */ + inline std::size_t size() const CADET_NOEXCEPT { return _base.size(); } + + /** + * @brief Returns whether base, linear, quadratic, and cubic arrays have the same size + * @return @c true if the arrays are of the same size, otherwise @c false + */ + inline bool allSameSize() const CADET_NOEXCEPT { return (_base.size() == _linear.size()) && (_base.size() == _quad.size()) && (_base.size() == _cube.size()); } + +protected: + std::vector _base; + std::vector _linear; + std::vector _quad; + std::vector _cube; +}; + + /** * @brief Component dependent scalar external parameter * @details A single value per component. @@ -998,6 +1238,95 @@ class ExternalScalarReactionDependentParameter }; +/** + * @brief Component dependent vectorial external parameter + * @details A vector per component, the component index changes the slowest. + */ +class ExternalVectorComponentDependentParameter +{ +public: + + typedef std::vector storage_t; + + inline void configure(const std::string& varName, IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBoundStates) + { + readScalarParameterOrArray(_base, paramProvider, "EXT_" + varName, 1); + readScalarParameterOrArray(_linear, paramProvider, "EXT_" + varName + "_T", 1); + readScalarParameterOrArray(_quad, paramProvider, "EXT_" + varName + "_TT", 1); + readScalarParameterOrArray(_cube, paramProvider, "EXT_" + varName + "_TTT", 1); + } + + inline void registerParam(const std::string& varName, std::unordered_map& parameters, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, unsigned int nComp, unsigned int const* nBoundStates) + { + const unsigned int stride = _base.size() / nComp; + const StringHash hashConst = hashStringRuntime("EXT_" + varName); + registerParam2DArray(parameters, _base, [=](bool multi, unsigned int comp, unsigned int rid) { return makeParamId(hashConst, unitOpIdx, comp, parTypeIdx, BoundStateIndep, multi ? rid : ReactionIndep, SectionIndep); }, stride); + + const StringHash hashLinear = hashStringRuntime("EXT_" + varName + "_T"); + registerParam2DArray(parameters, _linear, [=](bool multi, unsigned int comp, unsigned int rid) { return makeParamId(hashLinear, unitOpIdx, comp, parTypeIdx, BoundStateIndep, multi ? rid : ReactionIndep, SectionIndep); }, stride); + + const StringHash hashQuad = hashStringRuntime("EXT_" + varName + "_TT"); + registerParam2DArray(parameters, _quad, [=](bool multi, unsigned int comp, unsigned int rid) { return makeParamId(hashQuad, unitOpIdx, comp, parTypeIdx, BoundStateIndep, multi ? rid : ReactionIndep, SectionIndep); }, stride); + + const StringHash hashCube = hashStringRuntime("EXT_" + varName + "_TTT"); + registerParam2DArray(parameters, _cube, [=](bool multi, unsigned int comp, unsigned int rid) { return makeParamId(hashCube, unitOpIdx, comp, parTypeIdx, BoundStateIndep, multi ? rid : ReactionIndep, SectionIndep); }, stride); + } + + inline void reserve(unsigned int numElem, unsigned int numSlices, unsigned int nComp, unsigned int const* nBoundStates) { } + + inline void update(std::vector& result, double extVal, unsigned int nComp, unsigned int const* nBoundStates) const + { + update(result.data(), extVal, nComp, nBoundStates); + } + + inline void update(active* result, double extVal, unsigned int nComp, unsigned int const* nBoundStates) const + { + for (std::size_t i = 0; i < _base.size(); ++i) + result[i] = _base[i] + extVal * (_linear[i] + extVal * (_quad[i] + extVal * _cube[i])); + } + + inline void updateTimeDerivative(std::vector& result, double extVal, double extTimeDiff, unsigned int nComp, unsigned int const* nBoundStates) const + { + updateTimeDerivative(result.data(), extVal, extTimeDiff, nComp, nBoundStates); + } + + inline void updateTimeDerivative(active* result, double extVal, double extTimeDiff, unsigned int nComp, unsigned int const* nBoundStates) const + { + for (std::size_t i = 0; i < _base.size(); ++i) + result[i] = extTimeDiff * (static_cast(_linear[i]) + extVal * (2.0 * static_cast(_quad[i]) + 3.0 * extVal * static_cast(_cube[i]))); + } + + inline std::size_t size() const CADET_NOEXCEPT { return _base.size(); } + inline bool allSameSize() const CADET_NOEXCEPT { return (_base.size() == _linear.size()) && (_base.size() == _quad.size()) && (_base.size() == _cube.size()); } + + inline std::size_t additionalDynamicMemory(unsigned int nComp, unsigned int totalNumBoundStates, unsigned int const* nBoundStates) const CADET_NOEXCEPT { return nComp * sizeof(active) + alignof(active); } + + inline storage_t& base() CADET_NOEXCEPT { return _base; } + inline const storage_t& base() const CADET_NOEXCEPT { return _base; } + + inline storage_t& linear() CADET_NOEXCEPT { return _linear; } + inline const storage_t& linear() const CADET_NOEXCEPT { return _linear; } + + inline storage_t& quadratic() CADET_NOEXCEPT { return _quad; } + inline const storage_t& quadratic() const CADET_NOEXCEPT { return _quad; } + + inline storage_t& cubic() CADET_NOEXCEPT { return _cube; } + inline const storage_t& cubic() const CADET_NOEXCEPT { return _cube; } + + template + inline void prepareCache(T& cache, LinearBufferAllocator& buffer) const + { + cache.fromTemplate(buffer, _base); + } + +protected: + std::vector _base; + std::vector _linear; + std::vector _quad; + std::vector _cube; +}; + + /** * @brief Component and bound state dependent external parameter * @details A single value per component and bound state / binding site type. From 5fc9b89d427992b0b1a42d7cc594f5e249f83405 Mon Sep 17 00:00:00 2001 From: Samuel Leweke Date: Wed, 21 Jul 2021 02:07:38 +0200 Subject: [PATCH 2/8] Convert nu to piecewise cubic poly in GIEX model The characteristic charge nu is extended to a piecewise cubic polynomial in the GIEX adsorption model. The nu of each component can have different polynomial coefficients and breaks, but their amount must be the same over all components. If the breaks are not provided, a single global polynomial is applied. This maintains backwards compatibility. --- .../binding/generalized_ion_exchange.rst | 58 +++-- src/libcadet/Spline.hpp | 175 +++++++++++++++ .../binding/GeneralizedIonExchangeBinding.cpp | 207 ++++++++++++++++-- 3 files changed, 405 insertions(+), 35 deletions(-) create mode 100644 src/libcadet/Spline.hpp diff --git a/doc/interface/binding/generalized_ion_exchange.rst b/doc/interface/binding/generalized_ion_exchange.rst index c12b2e8ca..98fa74733 100644 --- a/doc/interface/binding/generalized_ion_exchange.rst +++ b/doc/interface/binding/generalized_ion_exchange.rst @@ -106,38 +106,70 @@ Generalized Ion Exchange **Unit:** :math:`m_{MP}^{3} mol^{-1}` -=================== ========================= +=================== ========================= **Type:** double **Length:** NCOMP =================== ========================= ``GIEX_NU`` Base value for characteristic charges of the protein; The number of sites :math:`\nu` that the protein interacts with on the resin - surface + surface. If more than NCOMP values are provided, :math:`\nu` is given + by a piecewise cubic polynomial. Data is expected in component-major + ordering. -=================== ========================= -**Type:** double **Length:** NCOMP -=================== ========================= +=================== =============================== +**Type:** double **Length:** multiples of NCOMP +=================== =============================== ``GIEX_NU_LIN`` Coefficient of linear dependence of characteristic charge on modifier - component + component. If more than NCOMP values are provided, :math:`\nu` is given + by a piecewise cubic polynomial. This parameter is optional, it defaults + to all 0. Data is expected in component-major ordering. **Unit:** :math:`\text{[Mod]}^{-1}` -=================== ========================= -**Type:** double **Length:** NCOMP -=================== ========================= +=================== =============================== +**Type:** double **Length:** same as GIEX_NU +=================== =============================== ``GIEX_NU_QUAD`` Coefficient of quadratic dependence of characteristic charge on - modifier component + modifier component. If more than NCOMP values are provided, + :math:`\nu` is given by a piecewise cubic polynomial. This parameter + is optional, it defaults to all 0. Data is expected in component-major + ordering. **Unit:** :math:`\text{[Mod]}^{-2}` -=================== ========================= -**Type:** double **Length:** NCOMP -=================== ========================= +=================== =============================== +**Type:** double **Length:** same as GIEX_NU +=================== =============================== + +``GIEX_NU_CUBE`` + Coefficient of cubic dependence of characteristic charge on + modifier component. If more than NCOMP values are provided, + :math:`\nu` is given by a piecewise cubic polynomial. This + parameter is optional, it defaults to all 0. Data is expected + in component-major ordering. + +**Unit:** :math:`\text{[Mod]}^{-3}` + +=================== =============================== +**Type:** double **Length:** same as GIEX_NU +=================== =============================== + +``GIEX_NU_BREAKS`` + Optional, only required if a piecewise cubic polynomial is + used for :math:`\nu`. Contains the breaks of the pieces + in component-major ordering. Note that :math:`n` pieces + have :math:`n+1` breaks. + +**Unit:** :math:`\text{[Mod]}` + +=================== ====================================== +**Type:** double **Length:** length of GIEX_NU + NCOMP +=================== ====================================== ``GIEX_SIGMA`` Steric factors of the protein; The number of sites :math:`\sigma` on diff --git a/src/libcadet/Spline.hpp b/src/libcadet/Spline.hpp new file mode 100644 index 000000000..469e3cb3a --- /dev/null +++ b/src/libcadet/Spline.hpp @@ -0,0 +1,175 @@ +// ============================================================================= +// 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 +// ============================================================================= + +#ifndef LIBCADET_SPLINE_HPP_ +#define LIBCADET_SPLINE_HPP_ + +#include +#include + +#include "cadet/cadetCompilerInfo.hpp" +#include "AutoDiff.hpp" + +namespace cadet +{ + + /** + * @brief Evaluates a piecewise cubic polynomial + * @details If the evaluation point @p x is outside the domain of the piecewise polynomial, + * constant extrapolation from the first and last value of the polynomial on its + * domain is applied. + * @param[in] x Evaluation position + * @param[in] breaks Array with @c nPieces+1 elements (strictly increasing) that defines the domain pieces + * @param[in] constCoeff Array with constant coefficients of size @p nPieces + * @param[in] linCoeff Array with linear coefficients of size @p nPieces + * @param[in] quadCoeff Array with quadratic coefficients of size @p nPieces + * @param[in] cubeCoeff Array with cubic coefficients of size @p nPieces + * @param[in] nPieces Number of pieces (at least 1) + * @tparam value_t Type of the evaluation point + * @tparam result_t Type of the result + * @return Value of the piecewise cubic polynomial at the given point @p x + */ + template + result_t evaluateCubicPiecewisePolynomial(const value_t& x, active const* breaks, active const* constCoeff, + active const* linCoeff, active const* quadCoeff, active const* cubeCoeff, int nPieces) CADET_NOEXCEPT + { + // Test if outside of domain, apply constant extrapolation + if (x < breaks[0]) + { + return constCoeff[0]; + } + + if (x >= breaks[nPieces]) + { + // Return the value at the right of the domain + const result_t y = static_cast(breaks[nPieces]) - static_cast(breaks[nPieces - 1]); + const int idx = nPieces - 1; + return static_cast(constCoeff[idx]) + y * ( + static_cast(linCoeff[idx]) + y * ( + static_cast(quadCoeff[idx]) + y * static_cast(cubeCoeff[idx]) + ) + ); + } + + // Find correct piece + int idx = 0; + for (; idx < nPieces; ++idx) + { + if (breaks[idx] >= x) + break; + } + + // Evaluate polynomial + const typename DoubleActivePromoter::type y = x - static_cast(breaks[idx]); + return static_cast(constCoeff[idx]) + y * ( + static_cast(linCoeff[idx]) + y * ( + static_cast(quadCoeff[idx]) + y * static_cast(cubeCoeff[idx]) + ) + ); + } + + /** + * @brief Evaluates a piecewise cubic polynomial and returns base and dependent values + * @details If the evaluation point @p x is outside the domain of the piecewise polynomial, + * constant extrapolation from the first and last value of the polynomial on its + * domain is applied. + * + * The function returns the constant base value and the variable part depending on @p x. + * @param[in] x Evaluation position + * @param[in] breaks Array with @c nPieces+1 elements (strictly increasing) that defines the domain pieces + * @param[in] constCoeff Array with constant coefficients of size @p nPieces + * @param[in] linCoeff Array with linear coefficients of size @p nPieces + * @param[in] quadCoeff Array with quadratic coefficients of size @p nPieces + * @param[in] cubeCoeff Array with cubic coefficients of size @p nPieces + * @param[in] nPieces Number of pieces (at least 1) + * @tparam value_t Type of the evaluation point + * @tparam result_t Type of the result + * @return Constant and dynamic value of the piecewise cubic polynomial at the given point @p x + */ + template + std::tuple evaluateCubicPiecewisePolynomialSplit(const value_t& x, active const* breaks, active const* constCoeff, + active const* linCoeff, active const* quadCoeff, active const* cubeCoeff, int nPieces) CADET_NOEXCEPT + { + // Test if outside of domain, apply constant extrapolation + if (x < breaks[0]) + { + return {static_cast(constCoeff[0]), result_t(0.0)}; + } + + if (x >= breaks[nPieces]) + { + // Return the value at the right of the domain + const result_t y = static_cast(breaks[nPieces]) - static_cast(breaks[nPieces - 1]); + const int idx = nPieces - 1; + return {static_cast(constCoeff[idx]) + y * ( + static_cast(linCoeff[idx]) + y * ( + static_cast(quadCoeff[idx]) + y * static_cast(cubeCoeff[idx]) + ) + ), result_t(0.0)}; + } + + // Find correct piece + int idx = 0; + for (; idx < nPieces; ++idx) + { + if (breaks[idx] >= x) + break; + } + + // Evaluate polynomial + const typename DoubleActivePromoter::type y = x - static_cast(breaks[idx]); + return {static_cast(constCoeff[idx]), y * ( + static_cast(linCoeff[idx]) + y * ( + static_cast(quadCoeff[idx]) + y * static_cast(cubeCoeff[idx]) + ) + )}; + } + + /** + * @brief Evaluates the derivative of a piecewise cubic polynomial + * @details If the evaluation point @p x is outside the domain of the piecewise polynomial, + * constant extrapolation from the first and last value of the polynomial on its + * domain is applied. + * @param[in] x Evaluation position + * @param[in] breaks Array with @c nPieces+1 elements (strictly increasing) that defines the domain pieces + * @param[in] linCoeff Array with linear coefficients of size @p nPieces + * @param[in] quadCoeff Array with quadratic coefficients of size @p nPieces + * @param[in] cubeCoeff Array with cubic coefficients of size @p nPieces + * @param[in] nPieces Number of pieces (at least 1) + * @return Derivative of the piecewise cubic polynomial at the given point @p x + */ + inline double evaluateCubicPiecewisePolynomialDerivative(double x, active const* breaks, active const* linCoeff, + active const* quadCoeff, active const* cubeCoeff, int nPieces) CADET_NOEXCEPT + { + // Test if outside of domain, apply constant extrapolation + if (x < breaks[0]) + return 0.0; + + if (x >= breaks[nPieces]) + return 0.0; + + // Find correct piece + int idx = 0; + for (; idx < nPieces; ++idx) + { + if (breaks[idx] >= x) + break; + } + + // Evaluate polynomial + const double y = x - static_cast(breaks[idx]); + return static_cast(linCoeff[idx]) + y * (2.0 * static_cast(quadCoeff[idx]) + 3.0 * y * static_cast(cubeCoeff[idx])); + } + +} // namespace cadet + +#endif // LIBCADET_SPLINE_HPP_ diff --git a/src/libcadet/model/binding/GeneralizedIonExchangeBinding.cpp b/src/libcadet/model/binding/GeneralizedIonExchangeBinding.cpp index e98598bbf..97d6413be 100644 --- a/src/libcadet/model/binding/GeneralizedIonExchangeBinding.cpp +++ b/src/libcadet/model/binding/GeneralizedIonExchangeBinding.cpp @@ -19,6 +19,7 @@ #include "model/Parameters.hpp" #include "LocalVector.hpp" #include "SimulationTypes.hpp" +#include "Spline.hpp" #include #include @@ -42,9 +43,11 @@ { "type": "ScalarComponentDependentParameter", "varName": "kDQuad", "confName": "GIEX_KD_QUAD"}, { "type": "ScalarComponentDependentParameter", "varName": "kDSalt", "confName": "GIEX_KD_SALT"}, { "type": "ScalarComponentDependentParameter", "varName": "kDProt", "confName": "GIEX_KD_PROT"}, - { "type": "ScalarComponentDependentParameter", "varName": "nu", "confName": "GIEX_NU"}, - { "type": "ScalarComponentDependentParameter", "varName": "nuLin", "confName": "GIEX_NU_LIN"}, - { "type": "ScalarComponentDependentParameter", "varName": "nuQuad", "confName": "GIEX_NU_QUAD"}, + { "type": "VectorComponentDependentParameter", "varName": "nu", "confName": "GIEX_NU"}, + { "type": "VectorComponentDependentParameter", "varName": "nuLin", "confName": "GIEX_NU_LIN", "skipConfig": true}, + { "type": "VectorComponentDependentParameter", "varName": "nuQuad", "confName": "GIEX_NU_QUAD", "skipConfig": true}, + { "type": "VectorComponentDependentParameter", "varName": "nuCube", "confName": "GIEX_NU_CUBE", "skipConfig": true}, + { "type": "VectorComponentDependentParameter", "varName": "nuBreaks", "confName": "GIEX_NU_BREAKS", "skipConfig": true}, { "type": "ScalarComponentDependentParameter", "varName": "sigma", "confName": "GIEX_SIGMA"} ], "constantParameters": @@ -66,6 +69,82 @@ refPhC0,refPhQ = Reference concentrations for pH dependent powers */ +namespace +{ + void assignZeros(cadet::model::VectorComponentDependentParameter& p, unsigned int size) + { + p.get() = std::vector(size, 0.0); + } + + void assignZeros(cadet::model::ExternalVectorComponentDependentParameter& p, unsigned int size) + { + p.base() = std::vector(size, 0.0); + p.linear() = std::vector(size, 0.0); + p.quadratic() = std::vector(size, 0.0); + p.cubic() = std::vector(size, 0.0); + } + + void assignSinglePiece(cadet::model::VectorComponentDependentParameter& p) + { + p.get().clear(); + } + + void assignSinglePiece(cadet::model::ExternalVectorComponentDependentParameter& p) + { + p.base().clear(); + p.linear().clear(); + p.quadratic().clear(); + p.cubic().clear(); + } + + template + std::tuple cpQNuPowers(int comp, int nPieces, const params_t& p, ph_t pH, cp_state_t cpBase, cp_state_t cpVar, q_state_t qBase, q_state_t qVar) + { + if (p.nuBreaks.size() == 0) + { + const cp_state_t nu_i_0_over_nu0 = static_cast(p.nu[comp]) / static_cast(p.nu[0]); + const cp_state_t nu_i_pH_over_nu0 = pH * (static_cast(p.nuLin[comp]) + pH * (static_cast(p.nuQuad[comp]) + pH * static_cast(p.nuCube[comp]))) / static_cast(p.nu[0]); + return {pow(cpBase, nu_i_0_over_nu0) * pow(cpVar, nu_i_pH_over_nu0), pow(qBase, nu_i_0_over_nu0) * pow(qVar, nu_i_pH_over_nu0)}; + } + else + { + const int offset = comp * nPieces; + const auto [nuConst, nuVar] = cadet::evaluateCubicPiecewisePolynomialSplit::type>(pH, p.nuBreaks.data() + comp * (nPieces + 1), p.nu.data() + offset, p.nuLin.data() + offset, p.nuQuad.data() + offset, p.nuCube.data() + offset, nPieces); + const cp_state_t nu_i_0_over_nu0 = static_cast(nuConst) / static_cast(p.nu[0]); + const cp_state_t nu_i_pH_over_nu0 = nuVar / static_cast(p.nu[0]); + return {pow(cpBase, nu_i_0_over_nu0) * pow(cpVar, nu_i_pH_over_nu0), pow(qBase, nu_i_0_over_nu0) * pow(qVar, nu_i_pH_over_nu0)}; + } + } + + template + std::tuple evaluateNu(int comp, int nPieces, const params_t& p, ph_t pH) + { + if (p.nuBreaks.size() == 0) + { + return {static_cast(p.nu[comp]), pH * (static_cast(p.nuLin[comp]) + pH * (static_cast(p.nuQuad[comp]) + pH * static_cast(p.nuCube[comp])))}; + } + else + { + const int offset = comp * nPieces; + return cadet::evaluateCubicPiecewisePolynomialSplit(pH, p.nuBreaks.data() + comp * (nPieces + 1), p.nu.data() + offset, p.nuLin.data() + offset, p.nuQuad.data() + offset, p.nuCube.data() + offset, nPieces); + } + } + + template + double nuDerivative(int comp, int nPieces, const params_t& p, double pH) + { + if (p.nuBreaks.size() == 0) + { + return static_cast(p.nuLin[comp]) + pH * (2.0 * static_cast(p.nuQuad[comp]) + 3.0 * pH * static_cast(p.nuCube[comp])); + } + else + { + const int offset = comp * nPieces; + return cadet::evaluateCubicPiecewisePolynomialDerivative(pH, p.nuBreaks.data() + comp * (nPieces + 1), p.nuLin.data() + offset, p.nuQuad.data() + offset, p.nuCube.data() + offset, nPieces); + } + } +} + namespace cadet { @@ -84,13 +163,35 @@ inline bool GIEXParamHandler::validateConfig(unsigned int nComp, unsigned int co if (_kD.size() < nComp) throw InvalidParameterException("GIEX_KD requires NCOMP entries"); if (_nu.size() < nComp) - throw InvalidParameterException("GIEX_NU requires NCOMP entries"); + throw InvalidParameterException("GIEX_NU requires at least NCOMP entries"); if (_sigma.size() < nComp) throw InvalidParameterException("GIEX_SIGMA requires NCOMP entries"); + if (_nu.size() != _nuLin.size()) + throw InvalidParameterException("GIEX_NU and GIEX_NU_LIN do not have the same length"); + + if (_nu.size() != _nuQuad.size()) + throw InvalidParameterException("GIEX_NU and GIEX_NU_QUAD do not have the same length"); + + if (_nu.size() != _nuCube.size()) + throw InvalidParameterException("GIEX_NU and GIEX_NU_CUBE do not have the same length"); + + if ((_nu.size() % nComp) != 0) + throw InvalidParameterException("Length of GIEX_NU must be a multiple of NCOMP"); + + if ((_nu.size() + nComp != _nuBreaks.size()) && (_nuBreaks.size() > 0)) + throw InvalidParameterException("GIEX_NU_BREAKS must have one entry more than polynomial pieces in GIEX_NU"); + + if ((_nu.size() != nComp) && (_nuBreaks.size() == 0)) + throw InvalidParameterException("GIEX_NU is expected to have NCOMP entries"); + // Assume monovalent salt ions by default - if (_nu.get()[0] <= 0.0) - _nu.get()[0] = 1.0; + const int nPieces = _nu.size() / nComp; + for (int i = 0; i < nPieces; ++i) + { + if (_nu.get()[i] <= 0.0) + _nu.get()[i] = 1.0; + } return true; } @@ -107,13 +208,46 @@ inline bool ExtGIEXParamHandler::validateConfig(unsigned int nComp, unsigned int if (_kD.size() < nComp) throw InvalidParameterException("EXT_GIEX_KD requires NCOMP entries"); if (_nu.size() < nComp) - throw InvalidParameterException("EXT_GIEX_NU requires NCOMP entries"); + throw InvalidParameterException("EXT_GIEX_NU requires at least NCOMP entries"); if (_sigma.size() < nComp) throw InvalidParameterException("EXT_GIEX_SIGMA requires NCOMP entries"); + if (_nu.size() != _nuLin.size()) + throw InvalidParameterException("EXT_GIEX_NU and EXT_GIEX_NU_LIN do not have the same length"); + + if (_nu.size() != _nuQuad.size()) + throw InvalidParameterException("EXT_GIEX_NU and EXT_GIEX_NU_QUAD do not have the same length"); + + if (_nu.size() != _nuCube.size()) + throw InvalidParameterException("EXT_GIEX_NU and EXT_GIEX_NU_CUBE do not have the same length"); + + if ((_nu.size() % nComp) != 0) + throw InvalidParameterException("Length of EXT_GIEX_NU must be a multiple of NCOMP"); + + if ((_nu.size() + nComp != _nuBreaks.size()) && (_nuBreaks.size() > 0)) + throw InvalidParameterException("EXT_GIEX_NU_BREAKS must have one entry more than polynomial pieces in EXT_GIEX_NU"); + + if ((_nu.size() != nComp) && (_nuBreaks.size() == 0)) + throw InvalidParameterException("EXT_GIEX_NU is expected to have NCOMP entries"); + + if (!_nu.allSameSize()) + throw InvalidParameterException("EXT_GIEX_NU, EXT_GIEX_NU_T, EXT_GIEX_NU_TT, and EXT_GIEX_NU_TTT must have the same size"); + if (!_nuLin.allSameSize()) + throw InvalidParameterException("EXT_GIEX_NU_LIN, EXT_GIEX_NU_LIN_T, EXT_GIEX_NU_LIN_TT, and EXT_GIEX_NU_LIN_TTT must have the same size"); + if (!_nuQuad.allSameSize()) + throw InvalidParameterException("EXT_GIEX_NU_QUAD, EXT_GIEX_NU_QUAD_T, EXT_GIEX_NU_QUAD_TT, and EXT_GIEX_NU_QUAD_TTT must have the same size"); + if (!_nuCube.allSameSize()) + throw InvalidParameterException("EXT_GIEX_NU_CUBE, EXT_GIEX_NU_CUBE_T, EXT_GIEX_NU_CUBE_TT, and EXT_GIEX_NU_CUBE_TTT must have the same size"); + if (!_nuBreaks.allSameSize()) + throw InvalidParameterException("EXT_GIEX_NU_BREAKS, EXT_GIEX_NU_BREAKS_T, EXT_GIEX_NU_BREAKS_TT, and EXT_GIEX_NU_BREAKS_TTT must have the same size"); + // Assume monovalent salt ions by default - if (_nu.base()[0] <= 0.0) - _nu.base()[0] = 1.0; + const int nPieces = _nu.size() / nComp; + for (int i = 0; i < nPieces; ++i) + { + if ((_nu.base()[i] <= 0.0) && (_nu.linear()[i] <= 0.0) && (_nu.quadratic()[i] <= 0.0) && (_nu.cubic()[i] <= 0.0)) + _nu.base()[i] = 1.0; + } return true; } @@ -220,6 +354,26 @@ class GeneralizedIonExchangeBindingBase : public ParamHandlerBindingModelBase::configureImpl(paramProvider, unitOpIdx, parTypeIdx); + if (paramProvider.exists(std::string(_paramHandler.prefixInConfiguration()) + "GIEX_NU_LIN")) + _paramHandler.nuLin().configure("GIEX_PH", paramProvider, _nComp, _nBoundStates); + else + assignZeros(_paramHandler.nuLin(), _paramHandler.nu().size()); + + if (paramProvider.exists(std::string(_paramHandler.prefixInConfiguration()) + "GIEX_NU_QUAD")) + _paramHandler.nuQuad().configure("GIEX_PH", paramProvider, _nComp, _nBoundStates); + else + assignZeros(_paramHandler.nuQuad(), _paramHandler.nu().size()); + + if (paramProvider.exists(std::string(_paramHandler.prefixInConfiguration()) + "GIEX_NU_CUBE")) + _paramHandler.nuCube().configure("GIEX_PH", paramProvider, _nComp, _nBoundStates); + else + assignZeros(_paramHandler.nuCube(), _paramHandler.nu().size()); + + if (paramProvider.exists(std::string(_paramHandler.prefixInConfiguration()) + "GIEX_NU_BREAKS")) + _paramHandler.nuBreaks().configure("GIEX_PH", paramProvider, _nComp, _nBoundStates); + else + assignSinglePiece(_paramHandler.nuBreaks()); + if (paramProvider.exists(std::string(_paramHandler.prefixInConfiguration()) + "GIEX_PHREFC0") && paramProvider.exists(std::string(_paramHandler.prefixInConfiguration()) + "GIEX_PHREFQ")) { // Parameters are present, use them @@ -241,11 +395,14 @@ class GeneralizedIonExchangeBindingBase : public ParamHandlerBindingModelBase::type; using StateParamType = typename DoubleActivePromoter::type; + using PhType = typename ActiveRefOrDouble::type; typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); + const int nPieces = (p->nuBreaks.size() > 0) ? (p->nuBreaks.size() - 1) : 0; + // Pseudo component 1 is pH - const CpStateType pH = yCp[1]; + const PhType pH = yCp[1]; // Salt flux: nu_0 * q_0 - Lambda + Sum[nu_j * q_j, j] == 0 // <=> nu_0 * q_0 == Lambda - Sum[nu_j * q_j, j] @@ -260,9 +417,9 @@ class GeneralizedIonExchangeBindingBase : public ParamHandlerBindingModelBase(p->nu[j]) + pH * (static_cast(p->nuLin[j]) + pH * static_cast(p->nuQuad[j])); + const auto [nuConst, nuVar] = evaluateNu(j, nPieces, *p, pH); - res[0] += nu_j * y[bndIdx]; + res[0] += (nuConst + nuVar) * y[bndIdx]; q0_bar -= static_cast(p->sigma[j]) * y[bndIdx]; // Next bound component @@ -289,11 +446,14 @@ class GeneralizedIonExchangeBindingBase : public ParamHandlerBindingModelBase(p->nu[i]) + pH * (static_cast(p->nuLin[i]) + pH * static_cast(p->nuQuad[i]))) / static_cast(p->nu[0]); // const CpStateParamType c0_pow_nu = pow(yCp0_divRef, nu_i_over_nu0); // const StateParamType q0_bar_pow_nu = pow(q0_bar_divRef, nu_i_over_nu0); - const CpStateParamType nu_i_0_over_nu0 = static_cast(p->nu[i]) / static_cast(p->nu[0]); - const CpStateParamType nu_i_pH_over_nu0 = pH * (static_cast(p->nuLin[i]) + pH * static_cast(p->nuQuad[i])) / static_cast(p->nu[0]); - const CpStateParamType c0_pow_nu = pow(yCp0_divRef, nu_i_0_over_nu0) * pow(yCp0_ph_divRef, nu_i_pH_over_nu0); - const StateParamType q0_bar_pow_nu = pow(q0_bar_divRef, nu_i_0_over_nu0) * pow(q0_bar_ph_divRef, nu_i_pH_over_nu0); + +// const CpStateParamType nu_i_0_over_nu0 = static_cast(p->nu[i]) / static_cast(p->nu[0]); +// const CpStateParamType nu_i_pH_over_nu0 = pH * (static_cast(p->nuLin[i]) + pH * static_cast(p->nuQuad[i])) / static_cast(p->nu[0]); +// const CpStateParamType c0_pow_nu = pow(yCp0_divRef, nu_i_0_over_nu0) * pow(yCp0_ph_divRef, nu_i_pH_over_nu0); +// const StateParamType q0_bar_pow_nu = pow(q0_bar_divRef, nu_i_0_over_nu0) * pow(q0_bar_ph_divRef, nu_i_pH_over_nu0); + const auto [c0_pow_nu, q0_bar_pow_nu] = cpQNuPowers(i, nPieces, *p, pH, yCp0_divRef, yCp0_ph_divRef, q0_bar_divRef, q0_bar_ph_divRef); + // k_{a,i}(c_p, q, \mathrm{pH}) = k_{a,i,0} \exp(k_{a,i,1} \mathrm{pH} + k_{a,i,2} \mathrm{pH}^2 + k_{a,i,\mathrm{salt}} c_{p,0} + k_{a,i,\mathrm{prot}} c_{p,i}) const CpStateParamType ka_i = static_cast(p->kA[i]) * exp(pH * (static_cast(p->kALin[i]) + pH * static_cast(p->kAQuad[i])) @@ -318,6 +478,7 @@ class GeneralizedIonExchangeBindingBase : public ParamHandlerBindingModelBasenuBreaks.size() > 0) ? (p->nuBreaks.size() - 1) : 0; // Pseudo component 1 is pH const double pH = yCp[1]; @@ -337,10 +498,10 @@ class GeneralizedIonExchangeBindingBase : public ParamHandlerBindingModelBase(p->nu[j]) + pH * (static_cast(p->nuLin[j]) + pH * static_cast(p->nuQuad[j])); + const auto [nuConst, nuVar] = evaluateNu(j, nPieces, *p, pH); - jac[bndIdx] = nu_j; - jac[1 - offsetCp] += (static_cast(p->nuLin[j]) + 2.0 * pH * static_cast(p->nuQuad[j])) * y[bndIdx]; + jac[bndIdx] = nuConst + nuVar; + jac[1 - offsetCp] += nuDerivative(j, nPieces, *p, pH) * y[bndIdx]; // Calculate \bar{q}_0 = nu_0 * q_0 - Sum[sigma_j * q_j, j] q0_bar -= static_cast(p->sigma[j]) * y[bndIdx]; @@ -370,12 +531,14 @@ class GeneralizedIonExchangeBindingBase : public ParamHandlerBindingModelBase(i, nPieces, *p, pH); + const double ka = static_cast(p->kA[i]); const double kd = static_cast(p->kD[i]); - const double nu_0 = static_cast(p->nu[i]) / static_cast(p->nu[0]); - const double nu_pH = pH * (static_cast(p->nuLin[i]) + pH * static_cast(p->nuQuad[i])) / static_cast(p->nu[0]); + const double nu_0 = nuConst / static_cast(p->nu[0]); + const double nu_pH = nuVar / static_cast(p->nu[0]); const double nu = nu_0 + nu_pH; - const double dNuDpH = (static_cast(p->nuLin[i]) + 2.0 * pH * static_cast(p->nuQuad[i])) / static_cast(p->nu[0]); + const double dNuDpH = nuDerivative(i, nPieces, *p, pH) / static_cast(p->nu[0]); const double c0_pow_nu = pow(yCp0_divRef, nu_0) * pow(yCp0_ph_divRef, nu_pH); const double q0_bar_pow_nu = pow(q0_bar_divRef, nu_0) * pow(q0_bar_ph_divRef, nu_pH); From 3558279f1277fa573d17948d305e482a0952eebd Mon Sep 17 00:00:00 2001 From: Samuel Leweke Date: Thu, 22 Jul 2021 00:56:10 +0200 Subject: [PATCH 3/8] Decouple validation and config in binding base Renames the validateConfig() function in the parameter handler to validate() and makes it public. Splits configure() into configure() and configureAndValidate(). This allows more fine-grained control, especially in situations where adsorption models require custom configuration code. --- .../model/binding/AntiLangmuirBinding.cpp | 4 ++-- .../model/binding/BiLangmuirBinding.cpp | 4 ++-- .../binding/BiStericMassActionBinding.cpp | 4 ++-- .../model/binding/BindingModelBase.hpp | 4 ++-- .../model/binding/ColloidalBinding.cpp | 4 ++-- ...dedMobilePhaseModulatorLangmuirBinding.cpp | 4 ++-- .../binding/ExternalFunctionTemplate.cpp | 22 ++++++++++++----- .../binding/GeneralizedIonExchangeBinding.cpp | 13 ++++++---- .../model/binding/KumarLangmuirBinding.cpp | 4 ++-- .../model/binding/LangmuirBinding.cpp | 4 ++-- src/libcadet/model/binding/LinearBinding.cpp | 24 +++++++++---------- .../MobilePhaseModulatorLangmuirBinding.cpp | 4 ++-- .../MultiComponentSpreadingBinding.cpp | 4 ++-- .../MultiStateStericMassActionBinding.cpp | 4 ++-- src/libcadet/model/binding/SaskaBinding.cpp | 4 ++-- .../model/binding/SelfAssociationBinding.cpp | 4 ++-- .../model/binding/StericMassActionBinding.cpp | 4 ++-- 17 files changed, 65 insertions(+), 50 deletions(-) diff --git a/src/libcadet/model/binding/AntiLangmuirBinding.cpp b/src/libcadet/model/binding/AntiLangmuirBinding.cpp index b9e3ce6e5..9bb5e2036 100644 --- a/src/libcadet/model/binding/AntiLangmuirBinding.cpp +++ b/src/libcadet/model/binding/AntiLangmuirBinding.cpp @@ -52,7 +52,7 @@ namespace model inline const char* AntiLangmuirParamHandler::identifier() CADET_NOEXCEPT { return "MULTI_COMPONENT_ANTILANGMUIR"; } -inline bool AntiLangmuirParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +inline bool AntiLangmuirParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates) { if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() != _antiLangmuir.size()) || (_kA.size() < nComp)) throw InvalidParameterException("MCAL_KA, MCAL_KD, MCAL_QMAX, and MCAL_ANTILANGMUIR have to have the same size"); @@ -62,7 +62,7 @@ inline bool AntiLangmuirParamHandler::validateConfig(unsigned int nComp, unsigne inline const char* ExtAntiLangmuirParamHandler::identifier() CADET_NOEXCEPT { return "EXT_MULTI_COMPONENT_ANTILANGMUIR"; } -inline bool ExtAntiLangmuirParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +inline bool ExtAntiLangmuirParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates) { if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() != _antiLangmuir.size()) || (_kA.size() < nComp)) throw InvalidParameterException("MCAL_KA, MCAL_KD, MCAL_QMAX, and MCAL_ANTILANGMUIR have to have the same size"); diff --git a/src/libcadet/model/binding/BiLangmuirBinding.cpp b/src/libcadet/model/binding/BiLangmuirBinding.cpp index 61cfe4bdb..e7a1ce260 100644 --- a/src/libcadet/model/binding/BiLangmuirBinding.cpp +++ b/src/libcadet/model/binding/BiLangmuirBinding.cpp @@ -51,7 +51,7 @@ namespace model inline const char* BiLangmuirParamHandler::identifier() CADET_NOEXCEPT { return "MULTI_COMPONENT_BILANGMUIR"; } -inline bool BiLangmuirParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +inline bool BiLangmuirParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates) { if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() < nComp)) throw InvalidParameterException("MCBL_KA, MCBL_KD, and MCBL_QMAX have to have the same size"); @@ -61,7 +61,7 @@ inline bool BiLangmuirParamHandler::validateConfig(unsigned int nComp, unsigned inline const char* ExtBiLangmuirParamHandler::identifier() CADET_NOEXCEPT { return "EXT_MULTI_COMPONENT_BILANGMUIR"; } -inline bool ExtBiLangmuirParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +inline bool ExtBiLangmuirParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates) { if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() < nComp)) throw InvalidParameterException("EXT_MCBL_KA, EXT_MCBL_KD, and EXT_MCBL_QMAX have to have the same size"); diff --git a/src/libcadet/model/binding/BiStericMassActionBinding.cpp b/src/libcadet/model/binding/BiStericMassActionBinding.cpp index 209b7b663..7ac2c975b 100644 --- a/src/libcadet/model/binding/BiStericMassActionBinding.cpp +++ b/src/libcadet/model/binding/BiStericMassActionBinding.cpp @@ -62,7 +62,7 @@ namespace model inline const char* BiSMAParamHandler::identifier() CADET_NOEXCEPT { return "BI_STERIC_MASS_ACTION"; } -inline bool BiSMAParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +inline bool BiSMAParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates) { const unsigned int numStates = firstNonEmptyBoundStates(nBoundStates, nComp); @@ -83,7 +83,7 @@ inline bool BiSMAParamHandler::validateConfig(unsigned int nComp, unsigned int c inline const char* ExtBiSMAParamHandler::identifier() CADET_NOEXCEPT { return "EXT_BI_STERIC_MASS_ACTION"; } -inline bool ExtBiSMAParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +inline bool ExtBiSMAParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates) { const unsigned int numStates = firstNonEmptyBoundStates(nBoundStates, nComp); diff --git a/src/libcadet/model/binding/BindingModelBase.hpp b/src/libcadet/model/binding/BindingModelBase.hpp index 85ef09f72..6c2faa996 100644 --- a/src/libcadet/model/binding/BindingModelBase.hpp +++ b/src/libcadet/model/binding/BindingModelBase.hpp @@ -133,12 +133,12 @@ class ParamHandlerBindingModelBase : public BindingModelBase virtual bool configureImpl(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx) { // Read parameters - _paramHandler.configure(paramProvider, _nComp, _nBoundStates); + const bool valid = _paramHandler.configureAndValidate(paramProvider, _nComp, _nBoundStates); // Register parameters _paramHandler.registerParameters(_parameters, unitOpIdx, parTypeIdx, _nComp, _nBoundStates); - return true; + return valid; } }; diff --git a/src/libcadet/model/binding/ColloidalBinding.cpp b/src/libcadet/model/binding/ColloidalBinding.cpp index 8089050bb..04c582ba4 100644 --- a/src/libcadet/model/binding/ColloidalBinding.cpp +++ b/src/libcadet/model/binding/ColloidalBinding.cpp @@ -71,7 +71,7 @@ namespace model inline const char* ColloidalParamHandler::identifier() CADET_NOEXCEPT { return "MULTI_COMPONENT_COLLOIDAL"; } -inline bool ColloidalParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +inline bool ColloidalParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates) { if ((_kEqPhExp.size() != _kEqSaltPowerExp.size()) || (_kEqPhExp.size() != _kEqSaltPowerFact.size()) @@ -97,7 +97,7 @@ inline bool ColloidalParamHandler::validateConfig(unsigned int nComp, unsigned i inline const char* ExtColloidalParamHandler::identifier() CADET_NOEXCEPT { return "EXT_MULTI_COMPONENT_COLLOIDAL"; } -inline bool ExtColloidalParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +inline bool ExtColloidalParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates) { if ((_kEqPhExp.size() != _kEqSaltPowerExp.size()) || (_kEqPhExp.size() != _kEqSaltPowerFact.size()) diff --git a/src/libcadet/model/binding/ExtendedMobilePhaseModulatorLangmuirBinding.cpp b/src/libcadet/model/binding/ExtendedMobilePhaseModulatorLangmuirBinding.cpp index 4d573a01e..14290e643 100644 --- a/src/libcadet/model/binding/ExtendedMobilePhaseModulatorLangmuirBinding.cpp +++ b/src/libcadet/model/binding/ExtendedMobilePhaseModulatorLangmuirBinding.cpp @@ -55,7 +55,7 @@ namespace model inline const char* EMPMLangmuirParamHandler::identifier() CADET_NOEXCEPT { return "EXTENDED_MOBILE_PHASE_MODULATOR"; } -inline bool EMPMLangmuirParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +inline bool EMPMLangmuirParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates) { if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() != _gamma.size()) || (_kA.size() != _beta.size()) || (_kA.size() < nComp)) @@ -66,7 +66,7 @@ inline bool EMPMLangmuirParamHandler::validateConfig(unsigned int nComp, unsigne inline const char* ExtEMPMLangmuirParamHandler::identifier() CADET_NOEXCEPT { return "EXT_EXTENDED_MOBILE_PHASE_MODULATOR"; } -inline bool ExtEMPMLangmuirParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +inline bool ExtEMPMLangmuirParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates) { if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() != _gamma.size()) || (_kA.size() != _beta.size()) || (_kA.size() < nComp)) diff --git a/src/libcadet/model/binding/ExternalFunctionTemplate.cpp b/src/libcadet/model/binding/ExternalFunctionTemplate.cpp index 7fc91b04f..688b57872 100644 --- a/src/libcadet/model/binding/ExternalFunctionTemplate.cpp +++ b/src/libcadet/model/binding/ExternalFunctionTemplate.cpp @@ -83,7 +83,7 @@ class {{ name }} : public cadet::model::ConstParamHandlerBase {% endif %} { } - inline bool configure(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBoundStates) + inline void configure(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBoundStates) { {% for p in parameters %} {% if not existsIn(p, "skipConfig") %} @@ -109,7 +109,12 @@ class {{ name }} : public cadet::model::ConstParamHandlerBase {% endif %} {% endfor %} {% endif %} - return validateConfig(nComp, nBoundStates); + } + + inline bool configureAndValidate(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBoundStates) + { + configure(paramProvider, nComp, nBoundStates); + return validate(nComp, nBoundStates); } inline void registerParameters(std::unordered_map& parameters, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, unsigned int nComp, unsigned int const* nBoundStates) @@ -171,9 +176,9 @@ class {{ name }} : public cadet::model::ConstParamHandlerBase {% endif %} inline char const* prefixInConfiguration() const CADET_NOEXCEPT { return ""; } + inline bool validate(unsigned int nComp, unsigned int const* nBoundStates); protected: - inline bool validateConfig(unsigned int nComp, unsigned int const* nBoundStates); ConstParams _localParams; @@ -262,7 +267,7 @@ class {{ externalName }} : public cadet::model::ExternalParamHandlerBase {% endif %} { } - inline bool configure(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBoundStates) + inline void configure(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBoundStates) { {% for p in parameters %} {% if not existsIn(p, "skipConfig") %} @@ -289,7 +294,12 @@ class {{ externalName }} : public cadet::model::ExternalParamHandlerBase {% endfor %} {% endif %} ExternalParamHandlerBase::configure(paramProvider, {{ length(parameters) }}); - return validateConfig(nComp, nBoundStates); + } + + inline bool configureAndValidate(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBoundStates) + { + configure(paramProvider, nComp, nBoundStates); + return validate(nComp, nBoundStates); } inline void registerParameters(std::unordered_map& parameters, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, unsigned int nComp, unsigned int const* nBoundStates) @@ -409,9 +419,9 @@ class {{ externalName }} : public cadet::model::ExternalParamHandlerBase {% endif %} inline char const* prefixInConfiguration() const CADET_NOEXCEPT { return "EXT_"; } + inline bool validate(unsigned int nComp, unsigned int const* nBoundStates); protected: - inline bool validateConfig(unsigned int nComp, unsigned int const* nBoundStates); ConstParams _constParams; diff --git a/src/libcadet/model/binding/GeneralizedIonExchangeBinding.cpp b/src/libcadet/model/binding/GeneralizedIonExchangeBinding.cpp index 97d6413be..c4a81d03e 100644 --- a/src/libcadet/model/binding/GeneralizedIonExchangeBinding.cpp +++ b/src/libcadet/model/binding/GeneralizedIonExchangeBinding.cpp @@ -153,7 +153,7 @@ namespace model inline const char* GIEXParamHandler::identifier() CADET_NOEXCEPT { return "GENERALIZED_ION_EXCHANGE"; } -inline bool GIEXParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +inline bool GIEXParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates) { if (nComp <= 2) throw InvalidParameterException("GENERALIZED_ION_EXCHANGE requires at least 3 components"); @@ -198,7 +198,7 @@ inline bool GIEXParamHandler::validateConfig(unsigned int nComp, unsigned int co inline const char* ExtGIEXParamHandler::identifier() CADET_NOEXCEPT { return "EXT_GENERALIZED_ION_EXCHANGE"; } -inline bool ExtGIEXParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +inline bool ExtGIEXParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates) { if (nComp <= 2) throw InvalidParameterException("EXT_GENERALIZED_ION_EXCHANGE requires at least 3 components"); @@ -349,10 +349,12 @@ class GeneralizedIonExchangeBindingBase : public ParamHandlerBindingModelBase::_reactionQuasistationarity; using ParamHandlerBindingModelBase::_nComp; using ParamHandlerBindingModelBase::_nBoundStates; + using ParamHandlerBindingModelBase::_parameters; virtual bool configureImpl(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx) { - const bool valid = ParamHandlerBindingModelBase::configureImpl(paramProvider, unitOpIdx, parTypeIdx); + // Read parameters + _paramHandler.configure(paramProvider, _nComp, _nBoundStates); if (paramProvider.exists(std::string(_paramHandler.prefixInConfiguration()) + "GIEX_NU_LIN")) _paramHandler.nuLin().configure("GIEX_PH", paramProvider, _nComp, _nBoundStates); @@ -386,7 +388,10 @@ class GeneralizedIonExchangeBindingBase : public ParamHandlerBindingModelBase diff --git a/src/libcadet/model/binding/KumarLangmuirBinding.cpp b/src/libcadet/model/binding/KumarLangmuirBinding.cpp index 93fefaa9c..6f6e4576d 100644 --- a/src/libcadet/model/binding/KumarLangmuirBinding.cpp +++ b/src/libcadet/model/binding/KumarLangmuirBinding.cpp @@ -58,7 +58,7 @@ namespace model inline const char* KumarLangmuirParamHandler::identifier() CADET_NOEXCEPT { return "KUMAR_MULTI_COMPONENT_LANGMUIR"; } -inline bool KumarLangmuirParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +inline bool KumarLangmuirParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates) { if ((_kA.size() != _kAct.size()) || (_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() != _nu.size()) || (_kA.size() < nComp)) throw InvalidParameterException("KMCL_KA, KMCL_KD, KMCL_KACT, KMCL_NU, and KMCL_QMAX have to have the same size"); @@ -68,7 +68,7 @@ inline bool KumarLangmuirParamHandler::validateConfig(unsigned int nComp, unsign inline const char* ExtKumarLangmuirParamHandler::identifier() CADET_NOEXCEPT { return "EXT_KUMAR_MULTI_COMPONENT_LANGMUIR"; } -inline bool ExtKumarLangmuirParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +inline bool ExtKumarLangmuirParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates) { if ((_kA.size() != _kAct.size()) || (_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() != _nu.size()) || (_kA.size() < nComp)) throw InvalidParameterException("KMCL_KA, KMCL_KD, KMCL_KACT, KMCL_NU, and KMCL_QMAX have to have the same size"); diff --git a/src/libcadet/model/binding/LangmuirBinding.cpp b/src/libcadet/model/binding/LangmuirBinding.cpp index dcc475a13..650fefc90 100644 --- a/src/libcadet/model/binding/LangmuirBinding.cpp +++ b/src/libcadet/model/binding/LangmuirBinding.cpp @@ -51,7 +51,7 @@ namespace model inline const char* LangmuirParamHandler::identifier() CADET_NOEXCEPT { return "MULTI_COMPONENT_LANGMUIR"; } -inline bool LangmuirParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +inline bool LangmuirParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates) { if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() < nComp)) throw InvalidParameterException("MCL_KA, MCL_KD, and MCL_QMAX have to have the same size"); @@ -61,7 +61,7 @@ inline bool LangmuirParamHandler::validateConfig(unsigned int nComp, unsigned in inline const char* ExtLangmuirParamHandler::identifier() CADET_NOEXCEPT { return "EXT_MULTI_COMPONENT_LANGMUIR"; } -inline bool ExtLangmuirParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +inline bool ExtLangmuirParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates) { if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() < nComp)) throw InvalidParameterException("EXT_MCL_KA, EXT_MCL_KD, and EXT_MCL_QMAX have to have the same size"); diff --git a/src/libcadet/model/binding/LinearBinding.cpp b/src/libcadet/model/binding/LinearBinding.cpp index 57dd58d46..6e03833fa 100644 --- a/src/libcadet/model/binding/LinearBinding.cpp +++ b/src/libcadet/model/binding/LinearBinding.cpp @@ -99,11 +99,11 @@ class LinearParamHandler : public ConstParamHandlerBase * @param [in] nBoundStates Array with number of bound states for each component * @return @c true if the parameters were read and validated successfully, otherwise @c false */ - inline bool configure(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBoundStates) + inline bool configureAndValidate(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBoundStates) { _kA.configure("LIN_KA", paramProvider, nComp, nBoundStates); _kD.configure("LIN_KD", paramProvider, nComp, nBoundStates); - return validateConfig(nComp, nBoundStates); + return validate(nComp, nBoundStates); } /** @@ -164,15 +164,13 @@ class LinearParamHandler : public ConstParamHandlerBase return std::make_tuple(&_localParams, nullptr); } -protected: - /** * @brief Validates recently read parameters * @param [in] nComp Number of components * @param [in] nBoundStates Array with number of bound states for each component * @return @c true if the parameters were validated successfully, otherwise @c false */ - inline bool validateConfig(unsigned int nComp, unsigned int const* nBoundStates) + inline bool validate(unsigned int nComp, unsigned int const* nBoundStates) { if ((_kA.size() != _kD.size()) || (_kA.size() < nComp)) throw InvalidParameterException("LIN_KA and LIN_KD have to have the same size"); @@ -180,6 +178,8 @@ class LinearParamHandler : public ConstParamHandlerBase return true; } +protected: + ConstParams _localParams; //!< Actual parameter data // Handlers provide configure(), reserve(), and registerParam() for parameters @@ -228,14 +228,14 @@ class ExtLinearParamHandler : public ExternalParamHandlerBase * @param [in] nBoundStates Array with number of bound states for each component * @return @c true if the parameters were read and validated successfully, otherwise @c false */ - inline bool configure(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBoundStates) + inline bool configureAndValidate(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBoundStates) { _kA.configure("LIN_KA", paramProvider, nComp, nBoundStates); _kD.configure("LIN_KD", paramProvider, nComp, nBoundStates); // Number of externally dependent parameters (2) needs to be given to ExternalParamHandlerBase::configure() ExternalParamHandlerBase::configure(paramProvider, 2); - return validateConfig(nComp, nBoundStates); + return validate(nComp, nBoundStates); } /** @@ -358,15 +358,13 @@ class ExtLinearParamHandler : public ExternalParamHandlerBase + 2 * (_kA.additionalDynamicMemory(nComp, totalNumBoundStates, nBoundStates) + _kD.additionalDynamicMemory(nComp, totalNumBoundStates, nBoundStates)); } -protected: - /** * @brief Validates recently read parameters * @param [in] nComp Number of components * @param [in] nBoundStates Array with number of bound states for each component * @return @c true if the parameters were validated successfully, otherwise @c false */ - inline bool validateConfig(unsigned int nComp, unsigned int const* nBoundStates) + inline bool validate(unsigned int nComp, unsigned int const* nBoundStates) { if ((_kA.size() != _kD.size()) || (_kA.size() < nComp)) throw InvalidParameterException("LIN_KA and LIN_KD have to have the same size"); @@ -374,6 +372,8 @@ class ExtLinearParamHandler : public ExternalParamHandlerBase return true; } +protected: + // Handlers provide configure(), reserve(), and registerParam() for parameters ExternalScalarComponentDependentParameter _kA; //!< Handler for adsorption rate ExternalScalarComponentDependentParameter _kD; //!< Handler for desorption rate @@ -446,12 +446,12 @@ class LinearBindingBase : public IBindingModel _parameters.clear(); // Read parameters (k_a and k_d) - _paramHandler.configure(paramProvider, _nComp, _nBoundStates); + const bool valid = _paramHandler.configureAndValidate(paramProvider, _nComp, _nBoundStates); // Register parameters _paramHandler.registerParameters(_parameters, unitOpIdx, parTypeIdx, _nComp, _nBoundStates); - return true; + return valid; } virtual void fillBoundPhaseInitialParameters(ParameterId* params, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx) const CADET_NOEXCEPT diff --git a/src/libcadet/model/binding/MobilePhaseModulatorLangmuirBinding.cpp b/src/libcadet/model/binding/MobilePhaseModulatorLangmuirBinding.cpp index 47ecb00a0..9c840c501 100644 --- a/src/libcadet/model/binding/MobilePhaseModulatorLangmuirBinding.cpp +++ b/src/libcadet/model/binding/MobilePhaseModulatorLangmuirBinding.cpp @@ -55,7 +55,7 @@ namespace model inline const char* MPMLangmuirParamHandler::identifier() CADET_NOEXCEPT { return "MOBILE_PHASE_MODULATOR"; } -inline bool MPMLangmuirParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +inline bool MPMLangmuirParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates) { if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() != _gamma.size()) || (_kA.size() != _beta.size()) || (_kA.size() < nComp)) @@ -66,7 +66,7 @@ inline bool MPMLangmuirParamHandler::validateConfig(unsigned int nComp, unsigned inline const char* ExtMPMLangmuirParamHandler::identifier() CADET_NOEXCEPT { return "EXT_MOBILE_PHASE_MODULATOR"; } -inline bool ExtMPMLangmuirParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +inline bool ExtMPMLangmuirParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates) { if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() != _gamma.size()) || (_kA.size() != _beta.size()) || (_kA.size() < nComp)) diff --git a/src/libcadet/model/binding/MultiComponentSpreadingBinding.cpp b/src/libcadet/model/binding/MultiComponentSpreadingBinding.cpp index dc324a47c..8fc12a84d 100644 --- a/src/libcadet/model/binding/MultiComponentSpreadingBinding.cpp +++ b/src/libcadet/model/binding/MultiComponentSpreadingBinding.cpp @@ -55,7 +55,7 @@ namespace model inline const char* SpreadingParamHandler::identifier() CADET_NOEXCEPT { return "MULTI_COMPONENT_SPREADING"; } -inline bool SpreadingParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +inline bool SpreadingParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates) { if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() < nComp * 2)) throw InvalidParameterException("MCSPR_KA, MCSPR_KD, and MCSPR_QMAX have to have the same size"); @@ -67,7 +67,7 @@ inline bool SpreadingParamHandler::validateConfig(unsigned int nComp, unsigned i inline const char* ExtSpreadingParamHandler::identifier() CADET_NOEXCEPT { return "EXT_MULTI_COMPONENT_SPREADING"; } -inline bool ExtSpreadingParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +inline bool ExtSpreadingParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates) { if ((_kA.size() != _kD.size()) || (_kA.size() != _qMax.size()) || (_kA.size() < nComp * 2)) throw InvalidParameterException("EXT_MCSPR_KA, EXT_MCSPR_KD, and EXT_MCSPR_QMAX have to have the same size"); diff --git a/src/libcadet/model/binding/MultiStateStericMassActionBinding.cpp b/src/libcadet/model/binding/MultiStateStericMassActionBinding.cpp index 118273cfe..b921fd739 100644 --- a/src/libcadet/model/binding/MultiStateStericMassActionBinding.cpp +++ b/src/libcadet/model/binding/MultiStateStericMassActionBinding.cpp @@ -64,7 +64,7 @@ namespace model inline const char* MSSMAParamHandler::identifier() CADET_NOEXCEPT { return "MULTISTATE_STERIC_MASS_ACTION"; } -inline bool MSSMAParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +inline bool MSSMAParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates) { const unsigned int totalBoundStates = numBoundStates(nBoundStates, nComp); @@ -91,7 +91,7 @@ inline bool MSSMAParamHandler::validateConfig(unsigned int nComp, unsigned int c inline const char* ExtMSSMAParamHandler::identifier() CADET_NOEXCEPT { return "EXT_MULTISTATE_STERIC_MASS_ACTION"; } -inline bool ExtMSSMAParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +inline bool ExtMSSMAParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates) { const unsigned int totalBoundStates = numBoundStates(nBoundStates, nComp); diff --git a/src/libcadet/model/binding/SaskaBinding.cpp b/src/libcadet/model/binding/SaskaBinding.cpp index a49c67212..89559e8ec 100644 --- a/src/libcadet/model/binding/SaskaBinding.cpp +++ b/src/libcadet/model/binding/SaskaBinding.cpp @@ -50,7 +50,7 @@ namespace model inline const char* SaskaParamHandler::identifier() CADET_NOEXCEPT { return "SASKA"; } -inline bool SaskaParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +inline bool SaskaParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates) { if ((_h.size() != _k.slices()) || (_k.size() != nComp * nComp) || (_h.size() < nComp)) throw InvalidParameterException("SASKA_K has to have NCOMP^2 and SASKA_H NCOMP many elements"); @@ -60,7 +60,7 @@ inline bool SaskaParamHandler::validateConfig(unsigned int nComp, unsigned int c inline const char* ExtSaskaParamHandler::identifier() CADET_NOEXCEPT { return "EXT_SASKA"; } -inline bool ExtSaskaParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +inline bool ExtSaskaParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates) { if ((_h.size() != _k.slices()) || (_k.size() != nComp * nComp) || (_h.size() < nComp)) throw InvalidParameterException("EXT_SASKA_K has to have NCOMP^2 and EXT_SASKA_H NCOMP many elements"); diff --git a/src/libcadet/model/binding/SelfAssociationBinding.cpp b/src/libcadet/model/binding/SelfAssociationBinding.cpp index 90aa8e365..885e21eec 100644 --- a/src/libcadet/model/binding/SelfAssociationBinding.cpp +++ b/src/libcadet/model/binding/SelfAssociationBinding.cpp @@ -64,7 +64,7 @@ namespace model inline const char* SelfAssociationParamHandler::identifier() CADET_NOEXCEPT { return "SELF_ASSOCIATION"; } -inline bool SelfAssociationParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +inline bool SelfAssociationParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates) { if ((_kA.size() != _kA2.size()) || (_kA.size() != _kD.size()) || (_kA.size() != _nu.size()) || (_kA.size() != _sigma.size()) || (_kA.size() < nComp)) @@ -79,7 +79,7 @@ inline bool SelfAssociationParamHandler::validateConfig(unsigned int nComp, unsi inline const char* ExtSelfAssociationParamHandler::identifier() CADET_NOEXCEPT { return "EXT_SELF_ASSOCIATION"; } -inline bool ExtSelfAssociationParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +inline bool ExtSelfAssociationParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates) { if ((_kA.size() != _kA2.size()) || (_kA.size() != _kD.size()) || (_kA.size() != _nu.size()) || (_kA.size() != _sigma.size()) || (_kA.size() < nComp)) diff --git a/src/libcadet/model/binding/StericMassActionBinding.cpp b/src/libcadet/model/binding/StericMassActionBinding.cpp index b783cd6b9..0ea5aeaae 100644 --- a/src/libcadet/model/binding/StericMassActionBinding.cpp +++ b/src/libcadet/model/binding/StericMassActionBinding.cpp @@ -62,7 +62,7 @@ namespace model inline const char* SMAParamHandler::identifier() CADET_NOEXCEPT { return "STERIC_MASS_ACTION"; } -inline bool SMAParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +inline bool SMAParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates) { if ((_kA.size() != _kD.size()) || (_kA.size() != _nu.size()) || (_kA.size() != _sigma.size()) || (_kA.size() < nComp)) throw InvalidParameterException("SMA_KA, SMA_KD, SMA_NU, and SMA_SIGMA have to have the same size"); @@ -76,7 +76,7 @@ inline bool SMAParamHandler::validateConfig(unsigned int nComp, unsigned int con inline const char* ExtSMAParamHandler::identifier() CADET_NOEXCEPT { return "EXT_STERIC_MASS_ACTION"; } -inline bool ExtSMAParamHandler::validateConfig(unsigned int nComp, unsigned int const* nBoundStates) +inline bool ExtSMAParamHandler::validate(unsigned int nComp, unsigned int const* nBoundStates) { if ((_kA.size() != _kD.size()) || (_kA.size() != _nu.size()) || (_kA.size() != _sigma.size()) || (_kA.size() < nComp)) throw InvalidParameterException("EXT_SMA_KA, EXT_SMA_KD, EXT_SMA_NU, and EXT_SMA_SIGMA have to have the same size"); From 4b69dc016a95a768d560d96b7d82b29248586cf6 Mon Sep 17 00:00:00 2001 From: Samuel Leweke Date: Thu, 22 Jul 2021 00:59:18 +0200 Subject: [PATCH 4/8] Decouple validation and config in reaction base Renames the validateConfig() function in the parameter handler to validate() and makes it public. Splits configure() into configure() and configureAndValidate(). This allows more fine-grained control, especially in situations where reaction models require custom configuration code. --- .../reaction/ExternalFunctionTemplate.cpp | 22 ++++++++++++++----- .../model/reaction/MassActionLawReaction.cpp | 6 ++--- 2 files changed, 19 insertions(+), 9 deletions(-) diff --git a/src/libcadet/model/reaction/ExternalFunctionTemplate.cpp b/src/libcadet/model/reaction/ExternalFunctionTemplate.cpp index d0ac4e16b..37f19b530 100644 --- a/src/libcadet/model/reaction/ExternalFunctionTemplate.cpp +++ b/src/libcadet/model/reaction/ExternalFunctionTemplate.cpp @@ -83,7 +83,7 @@ class {{ name }} : public cadet::model::ConstParamHandlerBase {% endif %} { } - inline bool configure(IParameterProvider& paramProvider, unsigned int nReactions, unsigned int nComp, unsigned int const* nBoundStates) + inline void configure(IParameterProvider& paramProvider, unsigned int nReactions, unsigned int nComp, unsigned int const* nBoundStates) { {% for p in parameters %} {% if not existsIn(p, "skipConfig") %} @@ -109,7 +109,12 @@ class {{ name }} : public cadet::model::ConstParamHandlerBase {% endif %} {% endfor %} {% endif %} - return validateConfig(nReactions, nComp, nBoundStates); + } + + inline bool configureAndValidate(IParameterProvider& paramProvider, unsigned int nReactions, unsigned int nComp, unsigned int const* nBoundStates) + { + configure(paramProvider, nReactions, nComp, nBoundStates); + return validate(nReactions, nComp, nBoundStates); } inline void registerParameters(std::unordered_map& parameters, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, unsigned int nComp, unsigned int const* nBoundStates) @@ -171,9 +176,9 @@ class {{ name }} : public cadet::model::ConstParamHandlerBase {% endif %} inline char const* prefixInConfiguration() const CADET_NOEXCEPT { return ""; } + inline bool validate(unsigned int nReactions, unsigned int nComp, unsigned int const* nBoundStates); protected: - inline bool validateConfig(unsigned int nReactions, unsigned int nComp, unsigned int const* nBoundStates); ConstParams _localParams; @@ -262,7 +267,7 @@ class {{ externalName }} : public cadet::model::ExternalParamHandlerBase {% endif %} { } - inline bool configure(IParameterProvider& paramProvider, unsigned int nReactions, unsigned int nComp, unsigned int const* nBoundStates) + inline void configure(IParameterProvider& paramProvider, unsigned int nReactions, unsigned int nComp, unsigned int const* nBoundStates) { {% for p in parameters %} {% if not existsIn(p, "skipConfig") %} @@ -289,7 +294,12 @@ class {{ externalName }} : public cadet::model::ExternalParamHandlerBase {% endfor %} {% endif %} ExternalParamHandlerBase::configure(paramProvider, {{ length(parameters) }}); - return validateConfig(nReactions, nComp, nBoundStates); + } + + inline bool configureAndValidate(IParameterProvider& paramProvider, unsigned int nReactions, unsigned int nComp, unsigned int const* nBoundStates) + { + configure(paramProvider, nReactions, nComp, nBoundStates); + return validate(nReactions, nComp, nBoundStates); } inline void registerParameters(std::unordered_map& parameters, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, unsigned int nComp, unsigned int const* nBoundStates) @@ -409,9 +419,9 @@ class {{ externalName }} : public cadet::model::ExternalParamHandlerBase {% endif %} inline char const* prefixInConfiguration() const CADET_NOEXCEPT { return "EXT_"; } + inline bool validate(unsigned int nReactions, unsigned int nComp, unsigned int const* nBoundStates); protected: - inline bool validateConfig(unsigned int nReactions, unsigned int nComp, unsigned int const* nBoundStates); ConstParams _constParams; diff --git a/src/libcadet/model/reaction/MassActionLawReaction.cpp b/src/libcadet/model/reaction/MassActionLawReaction.cpp index 3a550ae21..e1658802f 100644 --- a/src/libcadet/model/reaction/MassActionLawReaction.cpp +++ b/src/libcadet/model/reaction/MassActionLawReaction.cpp @@ -60,14 +60,14 @@ namespace model inline const char* MassActionLawParamHandler::identifier() CADET_NOEXCEPT { return "MASS_ACTION_LAW"; } -inline bool MassActionLawParamHandler::validateConfig(unsigned int nReactions, unsigned int nComp, unsigned int const* nBoundStates) +inline bool MassActionLawParamHandler::validate(unsigned int nReactions, unsigned int nComp, unsigned int const* nBoundStates) { return true; } inline const char* ExtMassActionLawParamHandler::identifier() CADET_NOEXCEPT { return "EXT_MASS_ACTION_LAW"; } -inline bool ExtMassActionLawParamHandler::validateConfig(unsigned int nReactions, unsigned int nComp, unsigned int const* nBoundStates) +inline bool ExtMassActionLawParamHandler::validate(unsigned int nReactions, unsigned int nComp, unsigned int const* nBoundStates) { return true; } @@ -546,7 +546,7 @@ class MassActionLawReactionBase : public DynamicReactionModelBase readAndRegisterExponents(paramProvider, _parameters, unitOpIdx, parTypeIdx, "MAL_EXPONENTS_SOLID_FWD_MODLIQUID", _expSolidFwdLiquid, _nComp, nullptr); readAndRegisterExponents(paramProvider, _parameters, unitOpIdx, parTypeIdx, "MAL_EXPONENTS_SOLID_BWD_MODLIQUID", _expSolidBwdLiquid, _nComp, nullptr); - return true; + return _paramHandler.validate(maxNumReactions(), _nComp, _nBoundStates); } template From 25c2ea796850cde927801c2f2b8765bdf41c1f60 Mon Sep 17 00:00:00 2001 From: Samuel Leweke Date: Thu, 22 Jul 2021 00:56:10 +0200 Subject: [PATCH 5/8] Fix bugs in GIEX parameter configuration Fixes copy and paste bugs in reading optional parameters of the GIEX adsorption model. Also adds a check for strict monotonicity of the polynomial breaks. --- .../binding/GeneralizedIonExchangeBinding.cpp | 27 ++++++++++++++----- 1 file changed, 20 insertions(+), 7 deletions(-) diff --git a/src/libcadet/model/binding/GeneralizedIonExchangeBinding.cpp b/src/libcadet/model/binding/GeneralizedIonExchangeBinding.cpp index c4a81d03e..351143068 100644 --- a/src/libcadet/model/binding/GeneralizedIonExchangeBinding.cpp +++ b/src/libcadet/model/binding/GeneralizedIonExchangeBinding.cpp @@ -193,6 +193,19 @@ inline bool GIEXParamHandler::validate(unsigned int nComp, unsigned int const* n _nu.get()[i] = 1.0; } + // Check breaks + for (int i = 0; i < nComp; ++i) + { + cadet::active const* const b = _nuBreaks.get().data() + nPieces * i; + for (int j = 0; j < nPieces; ++j) + { + if (b[j] >= b[j+1]) + { + throw InvalidParameterException("GIEX_NU_BREAKS must be strictly increasing for each component"); + } + } + } + return true; } @@ -311,6 +324,7 @@ class GeneralizedIonExchangeBindingBase : public ParamHandlerBindingModelBasenuBreaks.size() > 0) ? (p->nuBreaks.size() - 1) : 0; // Salt equation: nu_0 * q_0 - Lambda + Sum[nu_j(pH) * q_j, j] == 0 // <=> q_0 == (Lambda - Sum[nu_j(pH) * q_j, j]) / nu_0 @@ -323,9 +337,8 @@ class GeneralizedIonExchangeBindingBase : public ParamHandlerBindingModelBase(p->nu[j]) + pH * (static_cast(p->nuLin[j]) + pH * static_cast(p->nuQuad[j])); - - y[0] -= nu_j * y[bndIdx]; + const auto [nuConst, nuVar] = evaluateNu(j, nPieces, *p, pH); + y[0] -= (nuConst + nuVar) * y[bndIdx]; // Next bound component ++bndIdx; @@ -357,22 +370,22 @@ class GeneralizedIonExchangeBindingBase : public ParamHandlerBindingModelBase Date: Thu, 22 Jul 2021 13:33:23 +0200 Subject: [PATCH 6/8] Fix more bugs in GIEX adsorption model Fixes some more bugs in the GIEX adsorption model. --- .../model/binding/GeneralizedIonExchangeBinding.cpp | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/src/libcadet/model/binding/GeneralizedIonExchangeBinding.cpp b/src/libcadet/model/binding/GeneralizedIonExchangeBinding.cpp index 351143068..09843f3bc 100644 --- a/src/libcadet/model/binding/GeneralizedIonExchangeBinding.cpp +++ b/src/libcadet/model/binding/GeneralizedIonExchangeBinding.cpp @@ -194,14 +194,17 @@ inline bool GIEXParamHandler::validate(unsigned int nComp, unsigned int const* n } // Check breaks - for (int i = 0; i < nComp; ++i) + if (_nuBreaks.size() > 1) { - cadet::active const* const b = _nuBreaks.get().data() + nPieces * i; - for (int j = 0; j < nPieces; ++j) + for (int i = 0; i < nComp; ++i) { - if (b[j] >= b[j+1]) + cadet::active const* const b = _nuBreaks.get().data() + (nPieces + 1) * i; + for (int j = 0; j < nPieces; ++j) { - throw InvalidParameterException("GIEX_NU_BREAKS must be strictly increasing for each component"); + if (b[j] >= b[j+1]) + { + throw InvalidParameterException("GIEX_NU_BREAKS must be strictly increasing for each component"); + } } } } From b834e8c1c382decaec600edf382382ba1b6a4b8a Mon Sep 17 00:00:00 2001 From: Samuel Leweke Date: Thu, 22 Jul 2021 21:11:16 +0200 Subject: [PATCH 7/8] Fix bug in computing number of GIEX spline pieces Number of spline pieces was fundamentally wrong. --- .../model/binding/GeneralizedIonExchangeBinding.cpp | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/libcadet/model/binding/GeneralizedIonExchangeBinding.cpp b/src/libcadet/model/binding/GeneralizedIonExchangeBinding.cpp index 09843f3bc..34ccd7995 100644 --- a/src/libcadet/model/binding/GeneralizedIonExchangeBinding.cpp +++ b/src/libcadet/model/binding/GeneralizedIonExchangeBinding.cpp @@ -110,6 +110,7 @@ namespace { const int offset = comp * nPieces; const auto [nuConst, nuVar] = cadet::evaluateCubicPiecewisePolynomialSplit::type>(pH, p.nuBreaks.data() + comp * (nPieces + 1), p.nu.data() + offset, p.nuLin.data() + offset, p.nuQuad.data() + offset, p.nuCube.data() + offset, nPieces); + LOG(Debug) << "cpQNuPowers " << comp << ": at " << static_cast(pH) << ":" << static_cast(nuConst) << " + " << static_cast(nuVar) << " = " << static_cast(nuConst) + static_cast(nuVar); const cp_state_t nu_i_0_over_nu0 = static_cast(nuConst) / static_cast(p.nu[0]); const cp_state_t nu_i_pH_over_nu0 = nuVar / static_cast(p.nu[0]); return {pow(cpBase, nu_i_0_over_nu0) * pow(cpVar, nu_i_pH_over_nu0), pow(qBase, nu_i_0_over_nu0) * pow(qVar, nu_i_pH_over_nu0)}; @@ -126,7 +127,10 @@ namespace else { const int offset = comp * nPieces; - return cadet::evaluateCubicPiecewisePolynomialSplit(pH, p.nuBreaks.data() + comp * (nPieces + 1), p.nu.data() + offset, p.nuLin.data() + offset, p.nuQuad.data() + offset, p.nuCube.data() + offset, nPieces); +// return cadet::evaluateCubicPiecewisePolynomialSplit(pH, p.nuBreaks.data() + comp * (nPieces + 1), p.nu.data() + offset, p.nuLin.data() + offset, p.nuQuad.data() + offset, p.nuCube.data() + offset, nPieces); + const auto [nuConst, nuVar] = cadet::evaluateCubicPiecewisePolynomialSplit(pH, p.nuBreaks.data() + comp * (nPieces + 1), p.nu.data() + offset, p.nuLin.data() + offset, p.nuQuad.data() + offset, p.nuCube.data() + offset, nPieces); + LOG(Debug) << "evalNu " << comp << ": at " << static_cast(pH) << ": " << static_cast(nuConst) << " + " << static_cast(nuVar) << " = " << static_cast(nuConst) + static_cast(nuVar); + return std::make_tuple(nuConst, nuVar); } } @@ -327,7 +331,7 @@ class GeneralizedIonExchangeBindingBase : public ParamHandlerBindingModelBasenuBreaks.size() > 0) ? (p->nuBreaks.size() - 1) : 0; + const int nPieces = (p->nuBreaks.size() > 0) ? (p->nu.size() / _nComp) : 0; // Salt equation: nu_0 * q_0 - Lambda + Sum[nu_j(pH) * q_j, j] == 0 // <=> q_0 == (Lambda - Sum[nu_j(pH) * q_j, j]) / nu_0 @@ -420,7 +424,7 @@ class GeneralizedIonExchangeBindingBase : public ParamHandlerBindingModelBasenuBreaks.size() > 0) ? (p->nuBreaks.size() - 1) : 0; + const int nPieces = (p->nuBreaks.size() > 0) ? (p->nu.size() / _nComp) : 0; // Pseudo component 1 is pH const PhType pH = yCp[1]; @@ -499,7 +503,7 @@ class GeneralizedIonExchangeBindingBase : public ParamHandlerBindingModelBasenuBreaks.size() > 0) ? (p->nuBreaks.size() - 1) : 0; + const int nPieces = (p->nuBreaks.size() > 0) ? (p->nu.size() / _nComp) : 0; // Pseudo component 1 is pH const double pH = yCp[1]; From d804b30ec25bcc053dd626f2618ed3ff587cd92d Mon Sep 17 00:00:00 2001 From: Samuel Leweke Date: Thu, 22 Jul 2021 21:21:04 +0200 Subject: [PATCH 8/8] Remove debug logging --- .../model/binding/GeneralizedIonExchangeBinding.cpp | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/src/libcadet/model/binding/GeneralizedIonExchangeBinding.cpp b/src/libcadet/model/binding/GeneralizedIonExchangeBinding.cpp index 34ccd7995..3daec43da 100644 --- a/src/libcadet/model/binding/GeneralizedIonExchangeBinding.cpp +++ b/src/libcadet/model/binding/GeneralizedIonExchangeBinding.cpp @@ -110,7 +110,6 @@ namespace { const int offset = comp * nPieces; const auto [nuConst, nuVar] = cadet::evaluateCubicPiecewisePolynomialSplit::type>(pH, p.nuBreaks.data() + comp * (nPieces + 1), p.nu.data() + offset, p.nuLin.data() + offset, p.nuQuad.data() + offset, p.nuCube.data() + offset, nPieces); - LOG(Debug) << "cpQNuPowers " << comp << ": at " << static_cast(pH) << ":" << static_cast(nuConst) << " + " << static_cast(nuVar) << " = " << static_cast(nuConst) + static_cast(nuVar); const cp_state_t nu_i_0_over_nu0 = static_cast(nuConst) / static_cast(p.nu[0]); const cp_state_t nu_i_pH_over_nu0 = nuVar / static_cast(p.nu[0]); return {pow(cpBase, nu_i_0_over_nu0) * pow(cpVar, nu_i_pH_over_nu0), pow(qBase, nu_i_0_over_nu0) * pow(qVar, nu_i_pH_over_nu0)}; @@ -127,10 +126,7 @@ namespace else { const int offset = comp * nPieces; -// return cadet::evaluateCubicPiecewisePolynomialSplit(pH, p.nuBreaks.data() + comp * (nPieces + 1), p.nu.data() + offset, p.nuLin.data() + offset, p.nuQuad.data() + offset, p.nuCube.data() + offset, nPieces); - const auto [nuConst, nuVar] = cadet::evaluateCubicPiecewisePolynomialSplit(pH, p.nuBreaks.data() + comp * (nPieces + 1), p.nu.data() + offset, p.nuLin.data() + offset, p.nuQuad.data() + offset, p.nuCube.data() + offset, nPieces); - LOG(Debug) << "evalNu " << comp << ": at " << static_cast(pH) << ": " << static_cast(nuConst) << " + " << static_cast(nuVar) << " = " << static_cast(nuConst) + static_cast(nuVar); - return std::make_tuple(nuConst, nuVar); + return cadet::evaluateCubicPiecewisePolynomialSplit(pH, p.nuBreaks.data() + comp * (nPieces + 1), p.nu.data() + offset, p.nuLin.data() + offset, p.nuQuad.data() + offset, p.nuCube.data() + offset, nPieces); } }