From 1ec5abc5bfcd7f527903b0dd6eb961c378964bc5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Schm=C3=B6lder?= Date: Mon, 22 Nov 2021 13:02:22 +0100 Subject: [PATCH] Fix consistent initial state for linear isotherm. Previously, it was assumed that the linear algebraic equations are easily solved, not requiring a full nonlinear solver. However, it was shown that this can lead to situations in which mass is not conserved. To fix this, even this simple model needs to be fully initialized. --- src/libcadet/model/binding/LinearBinding.cpp | 33 ++------------------ 1 file changed, 3 insertions(+), 30 deletions(-) diff --git a/src/libcadet/model/binding/LinearBinding.cpp b/src/libcadet/model/binding/LinearBinding.cpp index 4e311383b..0016827f6 100644 --- a/src/libcadet/model/binding/LinearBinding.cpp +++ b/src/libcadet/model/binding/LinearBinding.cpp @@ -584,36 +584,9 @@ class LinearBindingBase : public IBindingModel } } - virtual bool preConsistentInitialState(double t, unsigned int secIdx, const ColumnPosition& colPos, double* y, double const* yCp, LinearBufferAllocator workSpace) const - { - // The linear algebraic equations are easily solved, that is, we don't require a full nonlinear solver - - typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); - - 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; - - // Skip dynamic binding reactions - if (!_reactionQuasistationarity[i]) - { - // Next bound component - ++bndIdx; - continue; - } - - // q = k_a / k_d * c_p - y[bndIdx] = static_cast(p->kA[i]) / static_cast(p->kD[i]) * yCp[i]; - - // Next bound component - ++bndIdx; - } - - // Consistent initialization is complete, don't call a nonlinear solver - return false; + virtual bool preConsistentInitialState(double t, unsigned int secIdx, const ColumnPosition& colPos, double* y, double const* yCp, LinearBufferAllocator workSpace) const { return true; } + { + // For a proper consistent initial state which conserves mass, nonlinear solver needs to be called. } virtual void postConsistentInitialState(double t, unsigned int secIdx, const ColumnPosition& colPos, double* y, double const* yCp, LinearBufferAllocator workSpace) const