From 6a4367e30a94050e25e87aecb2e414a984b2377a Mon Sep 17 00:00:00 2001 From: "a.berger" Date: Wed, 15 Jan 2025 15:24:36 +0100 Subject: [PATCH 01/11] Restart implementation --- src/libcadet/model/ReactionModel.hpp | 19 ++- src/libcadet/model/StirredTankModel.cpp | 110 ++++++++++-- src/libcadet/model/StirredTankModel.hpp | 7 + .../reaction/CrystallizationReaction.cpp | 5 +- src/libcadet/model/reaction/DummyReaction.cpp | 6 + .../model/reaction/MassActionLawReaction.cpp | 158 ++++++++++++++++++ .../reaction/MichaelisMentenReaction.cpp | 3 + .../model/reaction/ReactionModelBase.cpp | 1 + .../model/reaction/ReactionModelBase.hpp | 12 +- 9 files changed, 303 insertions(+), 18 deletions(-) diff --git a/src/libcadet/model/ReactionModel.hpp b/src/libcadet/model/ReactionModel.hpp index 0b2ab8397..861ea1589 100644 --- a/src/libcadet/model/ReactionModel.hpp +++ b/src/libcadet/model/ReactionModel.hpp @@ -213,6 +213,10 @@ class IDynamicReactionModel virtual int residualLiquidAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double* res, double factor, LinearBufferAllocator workSpace) const = 0; + virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, + active* fluxes, std::vector& mapQSReac, LinearBufferAllocator workSpace) const = 0; + + /** * @brief Adds the analytical Jacobian of the reaction terms for one liquid phase cell * @details Adds the Jacobian of the dynamic reaction terms for the given liquid phase @@ -237,9 +241,15 @@ class IDynamicReactionModel virtual void analyticJacobianLiquidAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, linalg::BandMatrix::RowIterator jac, LinearBufferAllocator workSpace) const = 0; virtual void analyticJacobianLiquidAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, linalg::DenseBandedRowIterator jac, LinearBufferAllocator workSpace) const = 0; virtual void analyticJacobianLiquidAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, linalg::BandedSparseRowIterator jac, LinearBufferAllocator workSpace) const = 0; - #ifdef ENABLE_DG + + + virtual void fillConservedMoietiesBulk(Eigen::MatrixXd& M, std::vector& QSReaction, std::vector& _QsCompBulk) = 0; + + +#ifdef ENABLE_DG virtual void analyticJacobianLiquidAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, linalg::BandedEigenSparseRowIterator jac, LinearBufferAllocator workSpace) const = 0; - #endif + +#endif /** * @brief Evaluates the residual for one combined phase cell @@ -314,6 +324,11 @@ class IDynamicReactionModel * @return Number of reactions */ virtual unsigned int numReactionsCombined() const CADET_NOEXCEPT = 0; + + /** + * @brief Returns wheather liquid reactions are in quasi-stationary state + * @return boolean + */ protected: }; diff --git a/src/libcadet/model/StirredTankModel.cpp b/src/libcadet/model/StirredTankModel.cpp index 671a59de1..15229e0b3 100644 --- a/src/libcadet/model/StirredTankModel.cpp +++ b/src/libcadet/model/StirredTankModel.cpp @@ -31,6 +31,10 @@ #include #include +#include + +#include + namespace cadet { @@ -76,7 +80,8 @@ CSTRModel::~CSTRModel() CADET_NOEXCEPT delete[] _strideBound; delete[] _offsetParType; - delete _dynReactionBulk; + delete[] _dynReactionBulk; + //delete[] _temp; -> Jan warum funktioniert das nicht? } unsigned int CSTRModel::numDofs() const CADET_NOEXCEPT @@ -235,6 +240,26 @@ bool CSTRModel::configureModelDiscretization(IParameterProvider& paramProvider, if (_dynReactionBulk->usesParamProviderInDiscretizationConfig()) paramProvider.popScope(); + + if (true) //paramProvider.exists("QUASI_STATIONARY_REACTION_BULK") + { + //_qsReacBulk = paramProvider.getIntArray("QUASI_STATIONARY_REACTION_BULK"); + _qsReacBulk = {1}; + _nQsReacBulk = _qsReacBulk.size(); + //_temp = new active[_nComp]; + + _nMoitiesBulk = (_nComp + _totalBound) - _nQsReacBulk; + _MconvMoityBulk = Eigen::MatrixXd::Zero(_nMoitiesBulk, _nComp); // matrix for conserved moities + } + else + { + _QsCompBulk.clear(); + _qsReacBulk.clear(); + _nMoitiesBulk = 0; + _nQsReacBulk = 0; + _MconvMoityBulk = Eigen::MatrixXd::Zero(0, 0); + + } } _dynReaction = std::vector(_nParType, nullptr); @@ -356,6 +381,12 @@ bool CSTRModel::configure(IParameterProvider& paramProvider) paramProvider.pushScope("reaction_bulk"); dynReactionConfSuccess = _dynReactionBulk->configure(paramProvider, _unitOpIdx, ParTypeIndep); paramProvider.popScope(); + + if (true) //paramProvider.exists("QUASI_STATIONARY_REACTION_BULK") + { + _dynReactionBulk->fillConservedMoietiesBulk(_MconvMoityBulk, _qsReacBulk, _QsCompBulk); // fill conserved moities matrix + int a = 0; + } } for (unsigned int type = 0; type < _nParType; ++type) @@ -422,9 +453,7 @@ unsigned int CSTRModel::threadLocalMemorySize() const CADET_NOEXCEPT return lms.bufferSize(); } -void CSTRModel::setSectionTimes(double const* secTimes, bool const* secContinuity, unsigned int nSections) -{ -} +void CSTRModel::setSectionTimes(double const* secTimes, bool const* secContinuity, unsigned int nSections){} void CSTRModel::useAnalyticJacobian(const bool analyticJac) { @@ -1247,7 +1276,7 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons const ParamType flowIn = static_cast(_flowRateIn); const ParamType flowOut = static_cast(_flowRateOut); - + // Inlet DOF for (unsigned int i = 0; i < _nComp; ++i) { @@ -1256,7 +1285,7 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons // Concentrations: \dot{V^l} * c + V^l * \dot{c} + V^s * sum_j sum_m \dot{q}_{j,m}] = c_in * F_in - c * F_out const ParamType vsolid = static_cast(_constSolidVolume); - ResidualType* const resC = res + _nComp; + ResidualType* resC = res + _nComp; for (unsigned int i = 0; i < _nComp; ++i) { resC[i] = 0.0; @@ -1312,14 +1341,73 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons // Reactions in liquid phase const ColumnPosition colPos{0.0, 0.0, 0.0}; - if (_dynReactionBulk && (_dynReactionBulk->numReactionsLiquid() > 0)) + if (_dynReactionBulk && (_nQsReacBulk > 0)) { - LinearBufferAllocator subAlloc = tlmAlloc.manageRemainingMemory(); - BufferedArray flux = subAlloc.array(_nComp); + LinearBufferAllocator subAlloc = tlmAlloc.manageRemainingMemory(); // todo nachgucken ob das wirklich so geht + BufferedArray flux = subAlloc.array(_nComp); std::fill_n(static_cast(flux), _nComp, 0.0); + + BufferedArray qsflux = subAlloc.array(_nQsReacBulk); + std::fill_n(static_cast(qsflux), _nQsReacBulk, 0.0); + + BufferedArray temp = subAlloc.array(_nComp); + std::fill_n(static_cast(temp), _nComp, 0.0); + _dynReactionBulk->residualLiquidAdd(t, secIdx, colPos, c, static_cast(flux), -1.0, subAlloc); + _dynReactionBulk->quasiStationaryFlux(t, secIdx, colPos, c, static_cast(qsflux), _qsReacBulk, subAlloc); + + + + Eigen::Map> resCMoities(static_cast(temp), _nComp); + resCMoities.setZero(); + + Eigen::Map> mapResC(resC, _nComp); + std::vector visitedQSComp(_nComp, 0); + int MoityIdx = 0; + + for (unsigned int state = 0; state < (_nComp - _nQsReacBulk); ++state) + { + if (_QsCompBulk[state] == 1) + { + ResidualType dotProduct = 0.0; + for (unsigned int i; i < _MconvMoityBulk.cols(); i++) // hier Optimierung durch Vermeidung von 0 Zeilen in MconvMoityBulk + { + dotProduct += static_cast(_MconvMoityBulk(MoityIdx, i)) * (mapResC[i]); + + if (wantJac) + _jac.native(i, state) = vDotTimeFactor + static_cast(flowOut); // dF_{ci}/dcj = v_liquidDot + F_out + } + resCMoities[state] = dotProduct; + + + MoityIdx++; + } + else if (_QsCompBulk[state] == 0) + { + resCMoities[state] += v * flux[state]; // hier sicher stellen, was in flux steht entweder res + flux oder nur flux + + if (wantJac) + _jac.native(state, _nComp + _totalBound) += static_cast(flux[state]); dF/dvliquid = flux + int a = 0; // add function that adds the jacobian for one state or change analyticJacobianLiquidAdd + } + } + + int state = (_nComp - _nQsReacBulk); + for (unsigned int qsreac = 0; qsreac < _nQsReacBulk; ++qsreac) + { + resCMoities[state++] = v * qsflux[qsreac]; + + if(wantJac) + int a = 0; // add function that adds the jacobian for single reactions -> maybe der klammer + } + + mapResC = resCMoities; + + } + else + { for (unsigned int comp = 0; comp < _nComp; ++comp) resC[comp] += v * flux[comp]; @@ -1331,7 +1419,7 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons _dynReactionBulk->analyticJacobianLiquidAdd(t, secIdx, colPos, reinterpret_cast(c), -static_cast(v), _jac.row(0), subAlloc); } } - + // Bound states for (unsigned int type = 0; type < _nParType; ++type) { @@ -1401,7 +1489,6 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons } } } - if (wantJac) { // Assemble Jacobian: Reaction @@ -1440,7 +1527,6 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons // Volume: \dot{V} = F_{in} - F_{out} - F_{filter} res[2 * _nComp + _totalBound] = vDot - flowIn + flowOut + static_cast(_curFlowRateFilter); - return 0; } diff --git a/src/libcadet/model/StirredTankModel.hpp b/src/libcadet/model/StirredTankModel.hpp index dc720b997..0fe5578ae 100644 --- a/src/libcadet/model/StirredTankModel.hpp +++ b/src/libcadet/model/StirredTankModel.hpp @@ -171,6 +171,13 @@ class CSTRModel : public UnitOperationBase IDynamicReactionModel* _dynReactionBulk; //!< Dynamic reactions in the bulk volume + // active* _temp; + Eigen::MatrixXd _MconvMoityBulk; //!< Matrix with conservation of moieties in the bulk volume + std::vector _qsReacBulk; //!< Indices of components that are not conserved in the bulk volume + std::vector _QsCompBulk; + unsigned int _nQsReacBulk; + unsigned int _nMoitiesBulk; + class Exporter : public ISolutionExporter { public: diff --git a/src/libcadet/model/reaction/CrystallizationReaction.cpp b/src/libcadet/model/reaction/CrystallizationReaction.cpp index 77eea953b..d9eadf942 100644 --- a/src/libcadet/model/reaction/CrystallizationReaction.cpp +++ b/src/libcadet/model/reaction/CrystallizationReaction.cpp @@ -197,7 +197,10 @@ namespace cadet virtual bool requiresConfiguration() const CADET_NOEXCEPT { return true; } virtual bool usesParamProviderInDiscretizationConfig() const CADET_NOEXCEPT { return false; } - + void fillConservedMoietiesBulk(Eigen::MatrixXd& M, std::vector& QSReaction, std::vector& QSComponent) { } + bool hasQuasiStationaryReactionsLiquid() { return false; } + virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, + active* fluxes, std::vector& mapQSReac, LinearBufferAllocator workSpace) const { return 0; } virtual bool configure(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx) { readScalarParameterOrArray(_bins, paramProvider, "CRY_BINS", 1); diff --git a/src/libcadet/model/reaction/DummyReaction.cpp b/src/libcadet/model/reaction/DummyReaction.cpp index 4b2581ede..1ac3786d1 100644 --- a/src/libcadet/model/reaction/DummyReaction.cpp +++ b/src/libcadet/model/reaction/DummyReaction.cpp @@ -99,6 +99,12 @@ class DummyDynamicReaction : public IDynamicReactionModel virtual unsigned int numReactionsLiquid() const CADET_NOEXCEPT { return 0; } virtual unsigned int numReactionsCombined() const CADET_NOEXCEPT { return 0; } + + void fillConservedMoietiesBulk(Eigen::MatrixXd& M, std::vector& QSReaction, std::vector& QSComponent) {} + bool hasQuasiStationaryReactionsLiquid() { return false;} + + virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, + active* fluxes, std::vector& mapQSReac, LinearBufferAllocator workSpace) const {return 0;} protected: }; diff --git a/src/libcadet/model/reaction/MassActionLawReaction.cpp b/src/libcadet/model/reaction/MassActionLawReaction.cpp index 2c11bff5d..498d6f3bd 100644 --- a/src/libcadet/model/reaction/MassActionLawReaction.cpp +++ b/src/libcadet/model/reaction/MassActionLawReaction.cpp @@ -25,7 +25,9 @@ #include #include #include +#include +using namespace Eigen; /* { "name": "MassActionLawParamHandler", @@ -312,11 +314,17 @@ class MassActionLawReactionBase : public DynamicReactionModelBase virtual void setExternalFunctions(IExternalFunction** extFuns, unsigned int size) { _paramHandler.setExternalFunctions(extFuns, size); } virtual bool dependsOnTime() const CADET_NOEXCEPT { return ParamHandler_t::dependsOnTime(); } virtual bool requiresWorkspace() const CADET_NOEXCEPT { return true; } + bool hasQuasiStationaryReactionsLiquid() { return true; } + virtual bool hasDynamicReactions() const CADET_NOEXCEPT { return true; } + virtual int const* reactionQuasiStationarity() const CADET_NOEXCEPT { return _reactionQuasistationarity.data(); } + virtual unsigned int workspaceSize(unsigned int nComp, unsigned int totalNumBoundStates, unsigned int const* nBoundStates) const CADET_NOEXCEPT { return _paramHandler.cacheSize(maxNumReactions(), nComp, totalNumBoundStates) + std::max(maxNumReactions() * sizeof(active), 2 * (_nComp + totalNumBoundStates) * sizeof(double)); } + + virtual bool configureModelDiscretization(IParameterProvider& paramProvider, unsigned int nComp, unsigned int const* nBound, unsigned int const* boundOffset) { DynamicReactionModelBase::configureModelDiscretization(paramProvider, nComp, nBound, boundOffset); @@ -371,6 +379,93 @@ class MassActionLawReactionBase : public DynamicReactionModelBase return true; } + + //void fillConservedMoietiesBulk(Eigen::MatrixXd M, std::vectorQSReaction, std::vectorQSComponent) { + // ConservedMoietiesLiquid(M, QSReaction,QSComponent); + //} + /* + * @brief Calculates the conserved moieties based on the stoichiometric matrix and reaction quasistationarity. + * @param S The stoichiometric matrix. + * @param _reactionQuasistationarity The reaction quasistationarity vector. + * @return The conserved moieties matrix. + */ + virtual void fillConservedMoietiesBulk(Eigen::MatrixXd& M, std::vector& mapQSReac, std::vector& _QsCompBulk) + { + + //1. get stoichmetic matrix with only reaction in quasi stationary + // S dim -> ncomp x nreac + + int nQS = std::count(mapQSReac.begin(), mapQSReac.end(), 1); + // Count the number of entries with value 1 in _reactionQuasistationarity + + + Eigen::MatrixXd QSS(_stoichiometryBulk.rows(), nQS); + int rowIndex = 0; + + // Fülle die Spalten basierend auf _reactionQuasistationarity + for (int i = 0; i < _stoichiometryBulk.columns(); ++i) + { + const double* rowData = reinterpret_cast(_stoichiometryBulk.data() + i * _stoichiometryBulk.columns()); + Eigen::Map Srow(rowData, _stoichiometryBulk.columns()); + QSS.row(rowIndex++) = Srow; + } + + //Remove zero rows from QSS + int nQScomp = 0; // Number of quasi stationary active components + for (int i = 0; i < QSS.rows(); ++i) + { + if (QSS.row(i).norm() < 1e-10) { + _QsCompBulk.push_back(0); + } + else + { + _QsCompBulk.push_back(1); + nQScomp++; + } + } + // Redimensioniere QSS und kopiere nur die nicht-null Zeilen + if (_nComp - nQScomp < 0.0 ) // if kinetic components exists we can resize QSS + { + Eigen::MatrixXd QSSCompressed(nQScomp, QSS.cols()); + for (std::size_t i = 0; i < _QsCompBulk.size(); ++i) + { + if (_QsCompBulk[i] == 1) + { + QSSCompressed.row(i) = QSS.row(_QsCompBulk[i]); + } + } + QSS.swap(QSSCompressed); // Swap the compressed matrix back into QSS + } + + //3. Test if the matrix is full rank + int rang = QSS.fullPivLu().rank(); + if (rang != nQS) + { + throw std::runtime_error("The matrix is not full rank"); + } + + //3. Calculate the null space of the matrix + Eigen::MatrixXd leftZeroSpace = QSS.transpose().fullPivLu().kernel().transpose(); + + if (_nComp - nQScomp < 1e-10) // if there are no kinetic components we can return the zero matrix + { + M = leftZeroSpace; + } + else // add zero rows for each component that is not in the quasi stationary + { + for (std::size_t i = 0; i < _QsCompBulk.size(); ++i) + { + if (_QsCompBulk[i] == 0) + { + leftZeroSpace.conservativeResize(NoChange, leftZeroSpace.cols() + 1); + leftZeroSpace.rightCols(leftZeroSpace.cols() - i - 1) = leftZeroSpace.middleCols(i, leftZeroSpace.cols() - i - 1); + leftZeroSpace.col(i).setZero(); + } + } + M = leftZeroSpace; + } + + } virtual unsigned int numReactionsLiquid() const CADET_NOEXCEPT { return _stoichiometryBulk.columns(); } virtual unsigned int numReactionsCombined() const CADET_NOEXCEPT { return _stoichiometryLiquid.columns() + _stoichiometrySolid.columns(); } @@ -396,6 +491,8 @@ class MassActionLawReactionBase : public DynamicReactionModelBase linalg::ActiveDenseMatrix _expSolidFwdLiquid; linalg::ActiveDenseMatrix _expSolidBwdLiquid; + std::vector _reactionQuasistationarity; + inline int maxNumReactions() const CADET_NOEXCEPT { return std::max(std::max(_stoichiometryBulk.columns(), _stoichiometryLiquid.columns()), _stoichiometrySolid.columns()); } virtual bool configureImpl(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx) @@ -549,6 +646,67 @@ class MassActionLawReactionBase : public DynamicReactionModelBase return true; } + int singleFlux(int r, active const* y, typename DoubleActivePromoter::type flow, typename ParamHandler_t::ParamsHandle const p) + { + typedef typename DoubleActivePromoter::type flux_t; + + flux_t fwd = rateConstantOrZero(static_cast::type>(p->kFwdBulk[r]), r, _expBulkFwd, _nComp); + for (int c = 0; c < _nComp; ++c) + { + if (_expBulkFwd.native(c, r) != 0.0) + { + if (static_cast(y[c]) > 0.0) + fwd *= pow(static_cast::type>(y[c]), + static_cast::type>(_expBulkFwd.native(c, r))); + else + { + fwd *= 0.0; + break; + } + } + } + + flux_t bwd = rateConstantOrZero(static_cast::type>(p->kBwdBulk[r]), r, _expBulkBwd, _nComp); + for (int c = 0; c < _nComp; ++c) + { + if (_expBulkBwd.native(c, r) != 0.0) + { + if (static_cast(y[c]) > 0.0) + bwd *= pow(static_cast::type>(y[c]), + static_cast::type>(_expBulkBwd.native(c, r))); + else + { + bwd *= 0.0; + break; + } + } + } + + flow = fwd - bwd; + return 0; + } + + + virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, + active* fluxes, std::vector& mapQSReac, LinearBufferAllocator workSpace) const + { + typedef typename DoubleActivePromoter::type flux_t; + typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); + + flux_t flow; + + for (int r = 0; r < _stoichiometryBulk.columns(); ++r) + { + if (mapQSReac[r] == 1) + { + //flow.clear() + singleFlux(r, y, flow, p); + fluxes[r] = flow; + } + return 0; + } + } + template int residualLiquidImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, StateType const* y, ResidualType* res, const FactorType& factor, LinearBufferAllocator workSpace) const diff --git a/src/libcadet/model/reaction/MichaelisMentenReaction.cpp b/src/libcadet/model/reaction/MichaelisMentenReaction.cpp index ee1f209b8..aad1ef500 100644 --- a/src/libcadet/model/reaction/MichaelisMentenReaction.cpp +++ b/src/libcadet/model/reaction/MichaelisMentenReaction.cpp @@ -152,6 +152,9 @@ class MichaelisMentenReactionBase : public DynamicReactionModelBase virtual unsigned int numReactionsLiquid() const CADET_NOEXCEPT { return _stoichiometryBulk.columns(); } virtual unsigned int numReactionsCombined() const CADET_NOEXCEPT { return 0; } + void fillConservedMoietiesBulk(Eigen::MatrixXd& M, std::vector& QSReaction, std::vector& QSComponent) {} + virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, + active* fluxes, std::vector& mapQSReac, LinearBufferAllocator workSpace) const { return 0; } CADET_DYNAMICREACTIONMODEL_BOILERPLATE diff --git a/src/libcadet/model/reaction/ReactionModelBase.cpp b/src/libcadet/model/reaction/ReactionModelBase.cpp index 71e59ca5f..e3d8f9979 100644 --- a/src/libcadet/model/reaction/ReactionModelBase.cpp +++ b/src/libcadet/model/reaction/ReactionModelBase.cpp @@ -87,6 +87,7 @@ bool DynamicReactionModelBase::setParameter(const ParameterId& pId, bool value) return false; } + active* DynamicReactionModelBase::getParameter(const ParameterId& pId) { auto paramHandle = _parameters.find(pId); diff --git a/src/libcadet/model/reaction/ReactionModelBase.hpp b/src/libcadet/model/reaction/ReactionModelBase.hpp index fe4f50724..83b54266c 100644 --- a/src/libcadet/model/reaction/ReactionModelBase.hpp +++ b/src/libcadet/model/reaction/ReactionModelBase.hpp @@ -19,8 +19,8 @@ #define LIBCADET_REACTIONMODELBASE_HPP_ #include "model/ReactionModel.hpp" +#include "linalg/ActiveDenseMatrix.hpp" #include "ParamIdUtil.hpp" - #include namespace cadet @@ -56,13 +56,18 @@ class DynamicReactionModelBase : public IDynamicReactionModel virtual active* getParameter(const ParameterId& pId); - virtual void setExternalFunctions(IExternalFunction** extFuns, unsigned int size) { } + virtual void setExternalFunctions(IExternalFunction** extFuns, unsigned int size) { }; + virtual void fillConservedMoietiesBulk(Eigen::MatrixXd& M, std::vector& QSReaction, std::vector& _QsCompBulk) = 0; + virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, + active* fluxes, std::vector& mapQSReac, LinearBufferAllocator workSpace) const = 0; + protected: int _nComp; //!< Number of components unsigned int const* _nBoundStates; //!< Array with number of bound states for each component unsigned int const* _boundOffset; //!< Array with offsets to the first bound state of each component int _nTotalBoundStates; + std::unordered_map _parameters; //!< Map used to translate ParameterIds to actual variables @@ -91,7 +96,7 @@ class DynamicReactionModelBase : public IDynamicReactionModel * The implementation is inserted inline in the class declaration. */ #ifdef ENABLE_DG -#define CADET_DYNAMICREACTIONMODEL_BOILERPLATE \ +#define CADET_DYNAMICREACTIONMODEL_BOILERPLATE \ virtual int residualLiquidAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, active const* y, \ active* res, const active& factor, LinearBufferAllocator workSpace) const \ { \ @@ -181,6 +186,7 @@ class DynamicReactionModelBase : public IDynamicReactionModel { \ jacobianCombinedImpl(t, secIdx, colPos, yLiquid, ySolid, factor, jacLiquid, jacSolid, workSpace); \ } + #else #define CADET_DYNAMICREACTIONMODEL_BOILERPLATE \ virtual int residualLiquidAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, active const* y, \ From e17e097986d862f899594f94769418ecc495d71d Mon Sep 17 00:00:00 2001 From: "a.berger" Date: Tue, 21 Jan 2025 12:13:48 +0100 Subject: [PATCH 02/11] changed residual, jacobi and jacobi derivativ to handle consvered moities --- src/libcadet/model/ReactionModel.hpp | 9 +- src/libcadet/model/StirredTankModel.cpp | 89 ++++++++++--------- src/libcadet/model/StirredTankModel.hpp | 4 +- .../reaction/CrystallizationReaction.cpp | 11 ++- src/libcadet/model/reaction/DummyReaction.cpp | 10 ++- .../model/reaction/MassActionLawReaction.cpp | 59 ++++++++---- .../reaction/MichaelisMentenReaction.cpp | 23 +++-- .../model/reaction/ReactionModelBase.hpp | 33 ++++++- 8 files changed, 164 insertions(+), 74 deletions(-) diff --git a/src/libcadet/model/ReactionModel.hpp b/src/libcadet/model/ReactionModel.hpp index 861ea1589..51fa4c6f7 100644 --- a/src/libcadet/model/ReactionModel.hpp +++ b/src/libcadet/model/ReactionModel.hpp @@ -212,9 +212,13 @@ class IDynamicReactionModel virtual int residualLiquidAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double* res, double factor, LinearBufferAllocator workSpace) const = 0; + + + virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, active const* y, + Eigen::Map fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) = 0; virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, - active* fluxes, std::vector& mapQSReac, LinearBufferAllocator workSpace) const = 0; + Eigen::Map fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) = 0; /** @@ -242,6 +246,9 @@ class IDynamicReactionModel virtual void analyticJacobianLiquidAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, linalg::DenseBandedRowIterator jac, LinearBufferAllocator workSpace) const = 0; virtual void analyticJacobianLiquidAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, linalg::BandedSparseRowIterator jac, LinearBufferAllocator workSpace) const = 0; + virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, linalg::BandMatrix::RowIterator jac, LinearBufferAllocator workSpace) const = 0; + virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, linalg::DenseBandedRowIterator jac, LinearBufferAllocator workSpace) const = 0; + virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, linalg::BandedSparseRowIterator jac, LinearBufferAllocator workSpace) const = 0; virtual void fillConservedMoietiesBulk(Eigen::MatrixXd& M, std::vector& QSReaction, std::vector& _QsCompBulk) = 0; diff --git a/src/libcadet/model/StirredTankModel.cpp b/src/libcadet/model/StirredTankModel.cpp index 15229e0b3..e5a975ebd 100644 --- a/src/libcadet/model/StirredTankModel.cpp +++ b/src/libcadet/model/StirredTankModel.cpp @@ -81,7 +81,8 @@ CSTRModel::~CSTRModel() CADET_NOEXCEPT delete[] _offsetParType; delete[] _dynReactionBulk; - //delete[] _temp; -> Jan warum funktioniert das nicht? + delete[] _temp; + delete[] _temp2; } unsigned int CSTRModel::numDofs() const CADET_NOEXCEPT @@ -241,19 +242,20 @@ bool CSTRModel::configureModelDiscretization(IParameterProvider& paramProvider, if (_dynReactionBulk->usesParamProviderInDiscretizationConfig()) paramProvider.popScope(); - if (true) //paramProvider.exists("QUASI_STATIONARY_REACTION_BULK") + if (paramProvider.exists("QUASI_STATIONARY_REACTION_BULK")) { - //_qsReacBulk = paramProvider.getIntArray("QUASI_STATIONARY_REACTION_BULK"); + _qsReacBulk = paramProvider.getIntArray("QUASI_STATIONARY_REACTION_BULK"); _qsReacBulk = {1}; _nQsReacBulk = _qsReacBulk.size(); - //_temp = new active[_nComp]; + _temp = new active[_nComp]; + _temp2 = new double[_nQsReacBulk]; _nMoitiesBulk = (_nComp + _totalBound) - _nQsReacBulk; _MconvMoityBulk = Eigen::MatrixXd::Zero(_nMoitiesBulk, _nComp); // matrix for conserved moities } else { - _QsCompBulk.clear(); + _QsCompBulk.clear(); _qsReacBulk.clear(); _nMoitiesBulk = 0; _nQsReacBulk = 0; @@ -382,10 +384,9 @@ bool CSTRModel::configure(IParameterProvider& paramProvider) dynReactionConfSuccess = _dynReactionBulk->configure(paramProvider, _unitOpIdx, ParTypeIndep); paramProvider.popScope(); - if (true) //paramProvider.exists("QUASI_STATIONARY_REACTION_BULK") + if (paramProvider.exists("QUASI_STATIONARY_REACTION_BULK")) { _dynReactionBulk->fillConservedMoietiesBulk(_MconvMoityBulk, _qsReacBulk, _QsCompBulk); // fill conserved moities matrix - int a = 0; } } @@ -718,7 +719,7 @@ void CSTRModel::consistentInitialState(const SimulationTime& simTime, double* co if (!qsMask[_nComp + bndIdx]) continue; - conservedQuants[idx] += factor * c[_nComp + bndIdx]; + conservedQuants[idx] += factor * c[_nComp + bndIdx]; } } @@ -1341,70 +1342,64 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons // Reactions in liquid phase const ColumnPosition colPos{0.0, 0.0, 0.0}; - if (_dynReactionBulk && (_nQsReacBulk > 0)) - { - LinearBufferAllocator subAlloc = tlmAlloc.manageRemainingMemory(); // todo nachgucken ob das wirklich so geht + LinearBufferAllocator subAlloc = tlmAlloc.manageRemainingMemory(); - BufferedArray flux = subAlloc.array(_nComp); - std::fill_n(static_cast(flux), _nComp, 0.0); + BufferedArray flux = subAlloc.array(_nComp); + std::fill_n(static_cast(flux), _nComp, 0.0); - BufferedArray qsflux = subAlloc.array(_nQsReacBulk); - std::fill_n(static_cast(qsflux), _nQsReacBulk, 0.0); + if (_dynReactionBulk && (_nQsReacBulk > 0)) + { + bool wantJac = true; + Eigen::Map> resCMoities(reinterpret_cast(_temp), _nComp); + resCMoities.setZero(); - BufferedArray temp = subAlloc.array(_nComp); - std::fill_n(static_cast(temp), _nComp, 0.0); + Eigen::Map> qsflux(_temp2, _nQsReacBulk); + qsflux.setZero(); _dynReactionBulk->residualLiquidAdd(t, secIdx, colPos, c, static_cast(flux), -1.0, subAlloc); - _dynReactionBulk->quasiStationaryFlux(t, secIdx, colPos, c, static_cast(qsflux), _qsReacBulk, subAlloc); - - - - Eigen::Map> resCMoities(static_cast(temp), _nComp); - resCMoities.setZero(); + _dynReactionBulk->quasiStationaryFlux(t, secIdx, colPos, c, qsflux, _qsReacBulk, subAlloc); Eigen::Map> mapResC(resC, _nComp); - std::vector visitedQSComp(_nComp, 0); int MoityIdx = 0; - for (unsigned int state = 0; state < (_nComp - _nQsReacBulk); ++state) { if (_QsCompBulk[state] == 1) { ResidualType dotProduct = 0.0; - for (unsigned int i; i < _MconvMoityBulk.cols(); i++) // hier Optimierung durch Vermeidung von 0 Zeilen in MconvMoityBulk + for (unsigned int i = 0 ; i < _MconvMoityBulk.cols(); ++i) // hier Optimierung durch Vermeidung von 0 Zeilen in MconvMoityBulk { dotProduct += static_cast(_MconvMoityBulk(MoityIdx, i)) * (mapResC[i]); - - if (wantJac) - _jac.native(i, state) = vDotTimeFactor + static_cast(flowOut); // dF_{ci}/dcj = v_liquidDot + F_out + if (wantJac && state != i) + { + _jac.native(state, i) = _MconvMoityBulk(MoityIdx, i)*(static_cast(vDot) + static_cast(flowOut)); // dF_{ci}/dcj = v_liquidDot + F_out + if (cadet_likely(yDot)) + _jac.native(state, _nComp + _totalBound) += cDot[i]; // ACHTUNG falsch fur M != 1 + } } resCMoities[state] = dotProduct; - - MoityIdx++; } else if (_QsCompBulk[state] == 0) { - resCMoities[state] += v * flux[state]; // hier sicher stellen, was in flux steht entweder res + flux oder nur flux - - if (wantJac) - _jac.native(state, _nComp + _totalBound) += static_cast(flux[state]); dF/dvliquid = flux - int a = 0; // add function that adds the jacobian for one state or change analyticJacobianLiquidAdd + resCMoities[state] = v * flux[state]; } } int state = (_nComp - _nQsReacBulk); for (unsigned int qsreac = 0; qsreac < _nQsReacBulk; ++qsreac) { - resCMoities[state++] = v * qsflux[qsreac]; + resCMoities[state] = v * qsflux[qsreac]; if(wantJac) - int a = 0; // add function that adds the jacobian for single reactions -> maybe der klammer + { + _dynReactionBulk->analyticQuasiSteadyJacobianLiquid(t, secIdx, colPos, reinterpret_cast(c), -static_cast(v), _jac.row(state), subAlloc); + _jac.native(state, _nComp + _totalBound) = 0.0; // dF_{ci}/dvliquid = 0 + } + state++; } mapResC = resCMoities; - } else { @@ -1919,10 +1914,22 @@ void CSTRModel::addTimeDerivativeJacobian(double t, double alpha, const ConstSim // Assemble Jacobian: dRes / dyDot // Concentrations: \dot{V^l} * c_i + V^l * \dot{c}_i + V^s * [sum_j sum_m d_j \dot{q}_{j,i,m}]) - c_{in,i} * F_in + c_i * F_out == 0 - for (unsigned int i = 0; i < _nComp; ++i) + for (unsigned int i = 0; i < _nComp - _nQsReacBulk; ++i) { mat.native(i, i) += timeV; // dRes / dcDot + if (_nQsReacBulk > 0) + { + for(unsigned int comp = 0; comp < _nComp; ++comp) + { + if(_QsCompBulk[comp] == 1 && i!=comp) + mat.native(i, comp) = timeV * static_cast(_MconvMoityBulk(i, comp)); // dRes / dcDot + mat.native(i, _nComp + _totalBound) += alpha * c[comp]; // dRes / dVlDot + } + + if(i >= _nComp - _nQsReacBulk) + break; + } double qSum = 0.0; for (unsigned int type = 0; type < _nParType; ++type) { @@ -1939,8 +1946,6 @@ void CSTRModel::addTimeDerivativeJacobian(double t, double alpha, const ConstSim } } - - mat.native(i, _nComp + _totalBound) += alpha * c[i]; // dRes / dVlDot } // Bound states diff --git a/src/libcadet/model/StirredTankModel.hpp b/src/libcadet/model/StirredTankModel.hpp index 0fe5578ae..fac1c0684 100644 --- a/src/libcadet/model/StirredTankModel.hpp +++ b/src/libcadet/model/StirredTankModel.hpp @@ -171,7 +171,9 @@ class CSTRModel : public UnitOperationBase IDynamicReactionModel* _dynReactionBulk; //!< Dynamic reactions in the bulk volume - // active* _temp; + active* _temp; + double* _temp2; + Eigen::MatrixXd _MconvMoityBulk; //!< Matrix with conservation of moieties in the bulk volume std::vector _qsReacBulk; //!< Indices of components that are not conserved in the bulk volume std::vector _QsCompBulk; diff --git a/src/libcadet/model/reaction/CrystallizationReaction.cpp b/src/libcadet/model/reaction/CrystallizationReaction.cpp index d9eadf942..a682246d7 100644 --- a/src/libcadet/model/reaction/CrystallizationReaction.cpp +++ b/src/libcadet/model/reaction/CrystallizationReaction.cpp @@ -199,8 +199,17 @@ namespace cadet virtual bool usesParamProviderInDiscretizationConfig() const CADET_NOEXCEPT { return false; } void fillConservedMoietiesBulk(Eigen::MatrixXd& M, std::vector& QSReaction, std::vector& QSComponent) { } bool hasQuasiStationaryReactionsLiquid() { return false; } + + template + void jacobianQuasiSteadyLiquidImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, const RowIterator& jac, LinearBufferAllocator workSpace) const { } + + virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, active const* y, + double* fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) { return 0; } + virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, - active* fluxes, std::vector& mapQSReac, LinearBufferAllocator workSpace) const { return 0; } + double* fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) { return 0; } + + virtual bool configure(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx) { readScalarParameterOrArray(_bins, paramProvider, "CRY_BINS", 1); diff --git a/src/libcadet/model/reaction/DummyReaction.cpp b/src/libcadet/model/reaction/DummyReaction.cpp index 1ac3786d1..3bd5ba325 100644 --- a/src/libcadet/model/reaction/DummyReaction.cpp +++ b/src/libcadet/model/reaction/DummyReaction.cpp @@ -77,6 +77,11 @@ class DummyDynamicReaction : public IDynamicReactionModel virtual void analyticJacobianLiquidAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, linalg::BandMatrix::RowIterator jac, LinearBufferAllocator workSpace) const { } virtual void analyticJacobianLiquidAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, linalg::DenseBandedRowIterator jac, LinearBufferAllocator workSpace) const { } virtual void analyticJacobianLiquidAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, linalg::BandedSparseRowIterator jac, LinearBufferAllocator workSpace) const { } + + virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, linalg::BandMatrix::RowIterator jac, LinearBufferAllocator workSpace) const { } + virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, linalg::DenseBandedRowIterator jac, LinearBufferAllocator workSpace) const { } + virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, linalg::BandedSparseRowIterator jac, LinearBufferAllocator workSpace) const { } + #ifdef ENABLE_DG virtual void analyticJacobianLiquidAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, linalg::BandedEigenSparseRowIterator jac, LinearBufferAllocator workSpace) const { } #endif @@ -103,8 +108,11 @@ class DummyDynamicReaction : public IDynamicReactionModel void fillConservedMoietiesBulk(Eigen::MatrixXd& M, std::vector& QSReaction, std::vector& QSComponent) {} bool hasQuasiStationaryReactionsLiquid() { return false;} + virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, active const* y, + Eigen::Map fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) { return 0; } + virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, - active* fluxes, std::vector& mapQSReac, LinearBufferAllocator workSpace) const {return 0;} + Eigen::Map fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) { return 0; } protected: }; diff --git a/src/libcadet/model/reaction/MassActionLawReaction.cpp b/src/libcadet/model/reaction/MassActionLawReaction.cpp index 498d6f3bd..6962f033e 100644 --- a/src/libcadet/model/reaction/MassActionLawReaction.cpp +++ b/src/libcadet/model/reaction/MassActionLawReaction.cpp @@ -403,7 +403,7 @@ class MassActionLawReactionBase : public DynamicReactionModelBase int rowIndex = 0; // Fülle die Spalten basierend auf _reactionQuasistationarity - for (int i = 0; i < _stoichiometryBulk.columns(); ++i) + for (int i = 0; i < _stoichiometryBulk.rows(); ++i) { const double* rowData = reinterpret_cast(_stoichiometryBulk.data() + i * _stoichiometryBulk.columns()); Eigen::Map Srow(rowData, _stoichiometryBulk.columns()); @@ -646,18 +646,15 @@ class MassActionLawReactionBase : public DynamicReactionModelBase return true; } - int singleFlux(int r, active const* y, typename DoubleActivePromoter::type flow, typename ParamHandler_t::ParamsHandle const p) + double singleFlux(int r, double const* y, double kFwdBulk_r, double kBwdBulk_r) { - typedef typename DoubleActivePromoter::type flux_t; - - flux_t fwd = rateConstantOrZero(static_cast::type>(p->kFwdBulk[r]), r, _expBulkFwd, _nComp); + double fwd = rateConstantOrZero(kFwdBulk_r, r, _expBulkFwd, _nComp); for (int c = 0; c < _nComp; ++c) { if (_expBulkFwd.native(c, r) != 0.0) { if (static_cast(y[c]) > 0.0) - fwd *= pow(static_cast::type>(y[c]), - static_cast::type>(_expBulkFwd.native(c, r))); + fwd *= pow(static_cast(y[c]), static_cast(_expBulkFwd.native(c, r))); else { fwd *= 0.0; @@ -666,14 +663,13 @@ class MassActionLawReactionBase : public DynamicReactionModelBase } } - flux_t bwd = rateConstantOrZero(static_cast::type>(p->kBwdBulk[r]), r, _expBulkBwd, _nComp); + double bwd = rateConstantOrZero(kBwdBulk_r, r, _expBulkBwd, _nComp); for (int c = 0; c < _nComp; ++c) { if (_expBulkBwd.native(c, r) != 0.0) { if (static_cast(y[c]) > 0.0) - bwd *= pow(static_cast::type>(y[c]), - static_cast::type>(_expBulkBwd.native(c, r))); + bwd *= pow(static_cast(y[c]), static_cast(_expBulkBwd.native(c, r))); else { bwd *= 0.0; @@ -681,26 +677,26 @@ class MassActionLawReactionBase : public DynamicReactionModelBase } } } - - flow = fwd - bwd; - return 0; + return fwd - bwd; } + virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, active const* y, + Eigen::Map fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) { return 0; } virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, - active* fluxes, std::vector& mapQSReac, LinearBufferAllocator workSpace) const + Eigen::Map fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) { - typedef typename DoubleActivePromoter::type flux_t; typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); - - flux_t flow; + for (int r = 0; r < _stoichiometryBulk.columns(); ++r) { + double kFwdBulk_r = static_cast(p->kFwdBulk[r]); + double kBwdBulk_r = static_cast(p->kBwdBulk[r]); + if (mapQSReac[r] == 1) { - //flow.clear() - singleFlux(r, y, flow, p); + double flow = singleFlux(r, y, kFwdBulk_r, kBwdBulk_r); fluxes[r] = flow; } return 0; @@ -936,6 +932,31 @@ class MassActionLawReactionBase : public DynamicReactionModelBase } } + template + void jacobianQuasiSteadyLiquidImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, const RowIterator& jac, LinearBufferAllocator workSpace) const + { + typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); + + BufferedArray fluxes = workSpace.array(2 * _nComp); + double* const fluxGradFwd = static_cast(fluxes); + double* const fluxGradBwd = fluxGradFwd + _nComp; + for (int r = 0; r < _stoichiometryBulk.columns(); ++r) + { + // Calculate gradients of forward and backward fluxes + fluxGradLiquid(fluxGradFwd, r, _nComp, static_cast(p->kFwdBulk[r]), _expBulkFwd, y); + fluxGradLiquid(fluxGradBwd, r, _nComp, static_cast(p->kBwdBulk[r]), _expBulkBwd, y); + + // Add gradients to Jacobian + RowIterator curJac = jac; // right row iterator + for (int row = 0; row < _nComp; ++row, ++curJac) + { + const double colFactor = static_cast(_stoichiometryBulk.native(row, r)) * factor; + for (int col = 0; col < _nComp; ++col) + curJac[0] = colFactor * (fluxGradFwd[col] - fluxGradBwd[col]); + } + } + } + template void jacobianCombinedImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* yLiquid, double const* ySolid, double factor, const RowIteratorLiquid& jacLiquid, const RowIteratorSolid& jacSolid, LinearBufferAllocator workSpace) const { diff --git a/src/libcadet/model/reaction/MichaelisMentenReaction.cpp b/src/libcadet/model/reaction/MichaelisMentenReaction.cpp index aad1ef500..9e733e004 100644 --- a/src/libcadet/model/reaction/MichaelisMentenReaction.cpp +++ b/src/libcadet/model/reaction/MichaelisMentenReaction.cpp @@ -153,8 +153,20 @@ class MichaelisMentenReactionBase : public DynamicReactionModelBase virtual unsigned int numReactionsLiquid() const CADET_NOEXCEPT { return _stoichiometryBulk.columns(); } virtual unsigned int numReactionsCombined() const CADET_NOEXCEPT { return 0; } void fillConservedMoietiesBulk(Eigen::MatrixXd& M, std::vector& QSReaction, std::vector& QSComponent) {} + + template + void jacobianQuasiSteadyLiquidImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, const RowIterator& jac, LinearBufferAllocator workSpace) const { } + + + virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, active const* y, + Eigen::Map fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) { + return 0; + } + virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, - active* fluxes, std::vector& mapQSReac, LinearBufferAllocator workSpace) const { return 0; } + Eigen::Map fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) { + return 0; + } CADET_DYNAMICREACTIONMODEL_BOILERPLATE @@ -171,7 +183,7 @@ class MichaelisMentenReactionBase : public DynamicReactionModelBase if ((_stoichiometryBulk.columns() > 0) && ((_paramHandler.vMax().size() < _stoichiometryBulk.columns()) || (_paramHandler.kMM().size() < _stoichiometryBulk.columns()))) throw InvalidParameterException("MM_VMAX and MM_KMM have to have the same size (number of reactions)"); - + if ((_stoichiometryBulk.columns() > 0) && (_paramHandler.kInhibit().size() < _stoichiometryBulk.columns() * _nComp)) throw InvalidParameterException("MM_KI have to have the size (number of reactions) x (number of components)"); @@ -222,7 +234,7 @@ class MichaelisMentenReactionBase : public DynamicReactionModelBase continue; } - fluxes[r] = static_cast::type>(p->vMax[r]) * y[idxSubs] / (static_cast::type>(p->kMM[r]) + y[idxSubs]); + fluxes[r] = static_cast::type>(p->vMax[r])* y[idxSubs] / (static_cast::type>(p->kMM[r]) + y[idxSubs]); for (int comp = 0; comp < _nComp; ++comp) { @@ -306,9 +318,8 @@ class MichaelisMentenReactionBase : public DynamicReactionModelBase } template - void jacobianCombinedImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* yLiquid, double const* ySolid, double factor, const RowIteratorLiquid& jacLiquid, const RowIteratorSolid& jacSolid, LinearBufferAllocator workSpace) const - { - } + void jacobianCombinedImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* yLiquid, double const* ySolid, double factor, const RowIteratorLiquid& jacLiquid, const RowIteratorSolid& jacSolid, LinearBufferAllocator workSpace) const {} + }; typedef MichaelisMentenReactionBase MichaelisMentenReaction; diff --git a/src/libcadet/model/reaction/ReactionModelBase.hpp b/src/libcadet/model/reaction/ReactionModelBase.hpp index 83b54266c..453ee71e4 100644 --- a/src/libcadet/model/reaction/ReactionModelBase.hpp +++ b/src/libcadet/model/reaction/ReactionModelBase.hpp @@ -59,8 +59,12 @@ class DynamicReactionModelBase : public IDynamicReactionModel virtual void setExternalFunctions(IExternalFunction** extFuns, unsigned int size) { }; virtual void fillConservedMoietiesBulk(Eigen::MatrixXd& M, std::vector& QSReaction, std::vector& _QsCompBulk) = 0; + + virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, active const* y, + Eigen::Map fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) = 0; + virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, - active* fluxes, std::vector& mapQSReac, LinearBufferAllocator workSpace) const = 0; + Eigen::Map fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) = 0; protected: int _nComp; //!< Number of components @@ -162,9 +166,32 @@ class DynamicReactionModelBase : public IDynamicReactionModel { \ jacobianLiquidImpl(t, secIdx, colPos, y, factor, jac, workSpace); \ } \ - \ + virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, \ + double factor, linalg::BandMatrix::RowIterator jac, LinearBufferAllocator workSpace) const \ + { \ + jacobianQuasiSteadyLiquidImpl(t, secIdx, colPos, y, factor, jac, workSpace); \ + } \ + \ + virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, \ + double factor, linalg::DenseBandedRowIterator jac, LinearBufferAllocator workSpace) const \ + { \ + jacobianQuasiSteadyLiquidImpl(t, secIdx, colPos, y, factor, jac, workSpace); \ + } \ + \ + virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, \ + double factor, linalg::BandedSparseRowIterator jac, LinearBufferAllocator workSpace) const \ + { \ + jacobianQuasiSteadyLiquidImpl(t, secIdx, colPos, y, factor, jac, workSpace); \ + } \ + \ + virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, \ + double factor, linalg::BandedEigenSparseRowIterator jac, LinearBufferAllocator workSpace) const \ + { \ + jacobianQuasiSteadyLiquidImpl(t, secIdx, colPos, y, factor, jac, workSpace); \ + } \ + \ virtual void analyticJacobianCombinedAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* yLiquid, double const* ySolid, \ - double factor, linalg::BandedEigenSparseRowIterator jacLiquid, linalg::DenseBandedRowIterator jacSolid, LinearBufferAllocator workSpace) const \ + double factor, linalg::BandedEigenSparseRowIterator jacLiquid, linalg::DenseBandedRowIterator jacSolid, LinearBufferAllocator workSpace) const \ { \ jacobianCombinedImpl(t, secIdx, colPos, yLiquid, ySolid, factor, jacLiquid, jacSolid, workSpace); \ } \ From bbb83bf94aa915502d6053fa929b6768a84adbbe Mon Sep 17 00:00:00 2001 From: "a.berger" Date: Fri, 24 Jan 2025 14:03:13 +0100 Subject: [PATCH 03/11] Progress of the week: add functionality for consisitent intialisation consisitent intitalisation does not work for singualr jacobian --- src/libcadet/Memory.hpp | 8 + src/libcadet/model/ReactionModel.hpp | 6 +- src/libcadet/model/StirredTankModel.cpp | 147 ++++++++++++++++-- src/libcadet/model/StirredTankModel.hpp | 2 +- .../reaction/CrystallizationReaction.cpp | 17 +- src/libcadet/model/reaction/DummyReaction.cpp | 6 +- .../model/reaction/MassActionLawReaction.cpp | 12 +- .../reaction/MichaelisMentenReaction.cpp | 4 +- .../model/reaction/ReactionModelBase.hpp | 16 +- 9 files changed, 175 insertions(+), 43 deletions(-) diff --git a/src/libcadet/Memory.hpp b/src/libcadet/Memory.hpp index 9125a7056..378e9a1c6 100644 --- a/src/libcadet/Memory.hpp +++ b/src/libcadet/Memory.hpp @@ -353,6 +353,14 @@ namespace cadet explicit operator T*() const { return _ptr; } + // method to copy elements from std::vector + void copyFromVector(const std::vector& vec) { + if (vec.size() > _numElements) { + throw std::out_of_range("Vector size exceeds BufferedArray size"); + } + std::copy(vec.begin(), vec.end(), _ptr); + } + private: BufferedArray(T* ptr, unsigned int numElements) : _ptr(ptr), _numElements(numElements) { } diff --git a/src/libcadet/model/ReactionModel.hpp b/src/libcadet/model/ReactionModel.hpp index 51fa4c6f7..6c57d22cd 100644 --- a/src/libcadet/model/ReactionModel.hpp +++ b/src/libcadet/model/ReactionModel.hpp @@ -246,9 +246,9 @@ class IDynamicReactionModel virtual void analyticJacobianLiquidAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, linalg::DenseBandedRowIterator jac, LinearBufferAllocator workSpace) const = 0; virtual void analyticJacobianLiquidAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, linalg::BandedSparseRowIterator jac, LinearBufferAllocator workSpace) const = 0; - virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, linalg::BandMatrix::RowIterator jac, LinearBufferAllocator workSpace) const = 0; - virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, linalg::DenseBandedRowIterator jac, LinearBufferAllocator workSpace) const = 0; - virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, linalg::BandedSparseRowIterator jac, LinearBufferAllocator workSpace) const = 0; + virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int nMoities, linalg::BandMatrix::RowIterator jac, LinearBufferAllocator workSpace) const = 0; + virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int nMoities, linalg::DenseBandedRowIterator jac, LinearBufferAllocator workSpace) const = 0; + virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int nMoities, linalg::BandedSparseRowIterator jac, LinearBufferAllocator workSpace) const = 0; virtual void fillConservedMoietiesBulk(Eigen::MatrixXd& M, std::vector& QSReaction, std::vector& _QsCompBulk) = 0; diff --git a/src/libcadet/model/StirredTankModel.cpp b/src/libcadet/model/StirredTankModel.cpp index e5a975ebd..414acf25d 100644 --- a/src/libcadet/model/StirredTankModel.cpp +++ b/src/libcadet/model/StirredTankModel.cpp @@ -242,9 +242,9 @@ bool CSTRModel::configureModelDiscretization(IParameterProvider& paramProvider, if (_dynReactionBulk->usesParamProviderInDiscretizationConfig()) paramProvider.popScope(); - if (paramProvider.exists("QUASI_STATIONARY_REACTION_BULK")) + if (true)//(paramProvider.exists("QUASI_STATIONARY_REACTION_BULK")) { - _qsReacBulk = paramProvider.getIntArray("QUASI_STATIONARY_REACTION_BULK"); + //_qsReacBulk = paramProvider.getIntArray("QUASI_STATIONARY_REACTION_BULK"); _qsReacBulk = {1}; _nQsReacBulk = _qsReacBulk.size(); _temp = new active[_nComp]; @@ -384,7 +384,7 @@ bool CSTRModel::configure(IParameterProvider& paramProvider) dynReactionConfSuccess = _dynReactionBulk->configure(paramProvider, _unitOpIdx, ParTypeIndep); paramProvider.popScope(); - if (paramProvider.exists("QUASI_STATIONARY_REACTION_BULK")) + if (true)//(paramProvider.exists("QUASI_STATIONARY_REACTION_BULK")) { _dynReactionBulk->fillConservedMoietiesBulk(_MconvMoityBulk, _qsReacBulk, _QsCompBulk); // fill conserved moities matrix } @@ -1075,11 +1075,11 @@ void CSTRModel::leanConsistentInitialState(const SimulationTime& simTime, double // It is assumed that the bound states are (approximately) correctly initialized. // Thus, only the liquid phase has to be initialized - double* const c = vecStateY + _nComp; - const double vLiquid = c[_nComp]; - const double vSolid = static_cast(_constSolidVolume); + double* const c = vecStateY + _nComp; + const double vLiquid = c[_nComp]; + const double vSolid = static_cast(_constSolidVolume); // Check if volume is 0 - if (vLiquid == 0.0) + if (vLiquid == 0.0) { const double flowIn = static_cast(_flowRateIn); const double flowOut = static_cast(_flowRateOut); @@ -1135,6 +1135,129 @@ void CSTRModel::leanConsistentInitialState(const SimulationTime& simTime, double } } } + if (_nQsReacBulk > 0) + { + LinearBufferAllocator tlmAlloc = threadLocalMem.get(); + BufferedArray qsMask = tlmAlloc.array(_nComp); + + qsMask.copyFromVector(_QsCompBulk); + + const linalg::ConstMaskArray mask{ static_cast(qsMask), static_cast(_nComp) }; + const int probSize = linalg::numMaskActive(mask); + + linalg::DenseMatrixView jacobianMatrix(_jacFact.data(), _jacFact.pivotData(), probSize, probSize); + BufferedArray baFullX = tlmAlloc.array(numDofs()); + double* const fullX = static_cast(baFullX); + + BufferedArray baFullResidual = tlmAlloc.array(numDofs()); + double* const fullResidual = static_cast(baFullResidual); + + BufferedArray baNonlinMem = tlmAlloc.array(_nonlinearSolver->workspaceSize(probSize)); + double* const nonlinMem = static_cast(baNonlinMem); + + // Extract initial values from current state + BufferedArray solution = tlmAlloc.array(probSize); + linalg::selectVectorSubset(c, mask, static_cast(solution)); + + std::function jacFunc; + if (adJac.adY && adJac.adRes) + { + ad::copyToAd(vecStateY, adJac.adY, _nComp); + + jacFunc = [&](double const* const x, linalg::detail::DenseMatrixBase& mat) + { + active* const localAdY = adJac.adY + _nComp; + active* const localAdRes = adJac.adRes + _nComp; + + // Copy over state vector to AD state vector (without changing directional values to keep seed vectors) + // and initialize residuals with zero (also resetting directional values) + ad::copyToAd(c, localAdY, mask.len); + // @todo Check if this is necessary + ad::resetAd(localAdRes, mask.len); + + // Prepare input vector by overwriting masked items + linalg::applyVectorSubset(x, mask, localAdY); + + // Call residual function + residualImpl(simTime.t, simTime.secIdx, adJac.adY, nullptr, adJac.adRes, tlmAlloc.manageRemainingMemory()); + +#ifdef CADET_CHECK_ANALYTIC_JACOBIAN + std::copy_n(c, mask.len, fullX); + linalg::applyVectorSubset(x, mask, fullX); + + // Compute analytic Jacobian + residualImpl(simTime.t, simTime.secIdx, fullX, nullptr, fullResidual, tlmAlloc.manageRemainingMemory()); + + // Compare + const double diff = ad::compareDenseJacobianWithAd(localAdRes, adJac.adDirOffset, _jac); + LOG(Debug) << "MaxDiff: " << diff; +#endif + + // Extract Jacobian from AD + ad::extractDenseJacobianFromAd(localAdRes, adJac.adDirOffset, _jac); + + // Extract Jacobian from full Jacobian + mat.setAll(0.0); + linalg::copyMatrixSubset(_jac, mask, mask, mat); + + return true; + }; + } + else + { + jacFunc = [&](double const* const x, linalg::detail::DenseMatrixBase& mat) + { + // Prepare input vector by overwriting masked items + std::copy_n(c, mask.len, fullX + _nComp); + fullX[(2 * _nComp + _totalBound)] = c[(_nComp + _totalBound)]; + + linalg::applyVectorSubset(x, mask, fullX + _nComp); + + // Call residual function to initialize jacobian + residualImpl(simTime.t, simTime.secIdx, fullX, nullptr, fullResidual, tlmAlloc.manageRemainingMemory()); + + // Extract Jacobian from full Jacobian + mat.setAll(0.0); + linalg::copyMatrixSubset(_jac, mask, mask, mat); + + return true; + }; + } + + _nonlinearSolver->solve( + [&](double const* const x, double* const r) + { + // Prepare input vector by overwriting masked items + std::copy_n(c, mask.len, fullX + _nComp); + fullX[(2*_nComp + _totalBound)] = c[(_nComp + _totalBound)]; + + linalg::applyVectorSubset(x, mask, fullX + _nComp); + + // Call residual function + residualImpl(simTime.t, simTime.secIdx, fullX, nullptr, fullResidual, tlmAlloc.manageRemainingMemory()); + + // Extract values from residual + linalg::selectVectorSubset(fullResidual + _nComp, mask, r); + + std::cout << "Residual: " << std::endl; + for (unsigned int i = 0; i < probSize; ++i){ + std::cout << r[i] << std::endl; + } + + return true; + }, + jacFunc, errorTol, static_cast(solution), nonlinMem, jacobianMatrix, probSize); + + // Apply solution + linalg::applyVectorSubset(static_cast(solution), mask, c); + std::cout << "Solution: " << std::endl; + for (unsigned int i = 0; i < _nComp; ++i) + std::cout << solution[i] << std::endl; + + int a = 1; + // Refine / correct solution + } + } void CSTRModel::leanConsistentInitialTimeDerivative(double t, double const* const vecStateY, double* const vecStateYdot, double* const res, util::ThreadLocalStorage& threadLocalMem) @@ -1277,7 +1400,7 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons const ParamType flowIn = static_cast(_flowRateIn); const ParamType flowOut = static_cast(_flowRateOut); - + bool wantJac = true; // just for debugg reasion // Inlet DOF for (unsigned int i = 0; i < _nComp; ++i) { @@ -1349,7 +1472,7 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons if (_dynReactionBulk && (_nQsReacBulk > 0)) { - bool wantJac = true; + Eigen::Map> resCMoities(reinterpret_cast(_temp), _nComp); resCMoities.setZero(); @@ -1374,7 +1497,7 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons { _jac.native(state, i) = _MconvMoityBulk(MoityIdx, i)*(static_cast(vDot) + static_cast(flowOut)); // dF_{ci}/dcj = v_liquidDot + F_out if (cadet_likely(yDot)) - _jac.native(state, _nComp + _totalBound) += cDot[i]; // ACHTUNG falsch fur M != 1 + _jac.native(state, _nComp + _totalBound) += _MconvMoityBulk(MoityIdx, i) * cDot[i]; // dF_{ci}/dVl = sum_i M[i,-] * dcidt } } resCMoities[state] = dotProduct; @@ -1393,8 +1516,9 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons if(wantJac) { - _dynReactionBulk->analyticQuasiSteadyJacobianLiquid(t, secIdx, colPos, reinterpret_cast(c), -static_cast(v), _jac.row(state), subAlloc); + _dynReactionBulk->analyticQuasiSteadyJacobianLiquid(t, secIdx, colPos, reinterpret_cast(c), state , _jac.row(state), subAlloc); _jac.native(state, _nComp + _totalBound) = 0.0; // dF_{ci}/dvliquid = 0 + } state++; } @@ -1920,6 +2044,7 @@ void CSTRModel::addTimeDerivativeJacobian(double t, double alpha, const ConstSim if (_nQsReacBulk > 0) { + mat.native(i, i) *= static_cast(_MconvMoityBulk(i, i)); // dRes / dcDot for(unsigned int comp = 0; comp < _nComp; ++comp) { if(_QsCompBulk[comp] == 1 && i!=comp) diff --git a/src/libcadet/model/StirredTankModel.hpp b/src/libcadet/model/StirredTankModel.hpp index fac1c0684..52e94378a 100644 --- a/src/libcadet/model/StirredTankModel.hpp +++ b/src/libcadet/model/StirredTankModel.hpp @@ -178,7 +178,7 @@ class CSTRModel : public UnitOperationBase std::vector _qsReacBulk; //!< Indices of components that are not conserved in the bulk volume std::vector _QsCompBulk; unsigned int _nQsReacBulk; - unsigned int _nMoitiesBulk; + int _nMoitiesBulk; class Exporter : public ISolutionExporter { diff --git a/src/libcadet/model/reaction/CrystallizationReaction.cpp b/src/libcadet/model/reaction/CrystallizationReaction.cpp index a682246d7..6f376891b 100644 --- a/src/libcadet/model/reaction/CrystallizationReaction.cpp +++ b/src/libcadet/model/reaction/CrystallizationReaction.cpp @@ -197,19 +197,20 @@ namespace cadet virtual bool requiresConfiguration() const CADET_NOEXCEPT { return true; } virtual bool usesParamProviderInDiscretizationConfig() const CADET_NOEXCEPT { return false; } - void fillConservedMoietiesBulk(Eigen::MatrixXd& M, std::vector& QSReaction, std::vector& QSComponent) { } - bool hasQuasiStationaryReactionsLiquid() { return false; } + void fillConservedMoietiesBulk(Eigen::MatrixXd& M, std::vector& QSReaction, std::vector& QSComponent) {} template - void jacobianQuasiSteadyLiquidImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, const RowIterator& jac, LinearBufferAllocator workSpace) const { } - - virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, active const* y, - double* fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) { return 0; } + void jacobianQuasiSteadyLiquidImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int nMoities, const RowIterator& jac, LinearBufferAllocator workSpace) const { } + virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, active const* y, + Eigen::Map fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) { + return 0; + } virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, - double* fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) { return 0; } + Eigen::Map fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) { + return 0; + } - virtual bool configure(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx) { readScalarParameterOrArray(_bins, paramProvider, "CRY_BINS", 1); diff --git a/src/libcadet/model/reaction/DummyReaction.cpp b/src/libcadet/model/reaction/DummyReaction.cpp index 3bd5ba325..17e8b8a83 100644 --- a/src/libcadet/model/reaction/DummyReaction.cpp +++ b/src/libcadet/model/reaction/DummyReaction.cpp @@ -78,9 +78,9 @@ class DummyDynamicReaction : public IDynamicReactionModel virtual void analyticJacobianLiquidAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, linalg::DenseBandedRowIterator jac, LinearBufferAllocator workSpace) const { } virtual void analyticJacobianLiquidAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, linalg::BandedSparseRowIterator jac, LinearBufferAllocator workSpace) const { } - virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, linalg::BandMatrix::RowIterator jac, LinearBufferAllocator workSpace) const { } - virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, linalg::DenseBandedRowIterator jac, LinearBufferAllocator workSpace) const { } - virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, linalg::BandedSparseRowIterator jac, LinearBufferAllocator workSpace) const { } + virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int nMoities , linalg::BandMatrix::RowIterator jac, LinearBufferAllocator workSpace) const { } + virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int nMoities, linalg::DenseBandedRowIterator jac, LinearBufferAllocator workSpace) const { } + virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int nMoities, linalg::BandedSparseRowIterator jac, LinearBufferAllocator workSpace) const { } #ifdef ENABLE_DG virtual void analyticJacobianLiquidAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, linalg::BandedEigenSparseRowIterator jac, LinearBufferAllocator workSpace) const { } diff --git a/src/libcadet/model/reaction/MassActionLawReaction.cpp b/src/libcadet/model/reaction/MassActionLawReaction.cpp index 6962f033e..a03384787 100644 --- a/src/libcadet/model/reaction/MassActionLawReaction.cpp +++ b/src/libcadet/model/reaction/MassActionLawReaction.cpp @@ -933,14 +933,14 @@ class MassActionLawReactionBase : public DynamicReactionModelBase } template - void jacobianQuasiSteadyLiquidImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, const RowIterator& jac, LinearBufferAllocator workSpace) const + void jacobianQuasiSteadyLiquidImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, const RowIterator& jac, LinearBufferAllocator workSpace) const { typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); BufferedArray fluxes = workSpace.array(2 * _nComp); double* const fluxGradFwd = static_cast(fluxes); double* const fluxGradBwd = fluxGradFwd + _nComp; - for (int r = 0; r < _stoichiometryBulk.columns(); ++r) + for (int r = 0; r < _stoichiometryBulk.columns(); ++r) // reactions { // Calculate gradients of forward and backward fluxes fluxGradLiquid(fluxGradFwd, r, _nComp, static_cast(p->kFwdBulk[r]), _expBulkFwd, y); @@ -948,11 +948,11 @@ class MassActionLawReactionBase : public DynamicReactionModelBase // Add gradients to Jacobian RowIterator curJac = jac; // right row iterator - for (int row = 0; row < _nComp; ++row, ++curJac) + + const double colFactor = static_cast(_stoichiometryBulk.native(state, r)); + for (int col = 0; col < _nComp; ++col) { - const double colFactor = static_cast(_stoichiometryBulk.native(row, r)) * factor; - for (int col = 0; col < _nComp; ++col) - curJac[0] = colFactor * (fluxGradFwd[col] - fluxGradBwd[col]); + curJac[col - state] = colFactor * (fluxGradFwd[col] - fluxGradBwd[col]); } } } diff --git a/src/libcadet/model/reaction/MichaelisMentenReaction.cpp b/src/libcadet/model/reaction/MichaelisMentenReaction.cpp index 9e733e004..b0525c9d7 100644 --- a/src/libcadet/model/reaction/MichaelisMentenReaction.cpp +++ b/src/libcadet/model/reaction/MichaelisMentenReaction.cpp @@ -155,14 +155,12 @@ class MichaelisMentenReactionBase : public DynamicReactionModelBase void fillConservedMoietiesBulk(Eigen::MatrixXd& M, std::vector& QSReaction, std::vector& QSComponent) {} template - void jacobianQuasiSteadyLiquidImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, const RowIterator& jac, LinearBufferAllocator workSpace) const { } - + void jacobianQuasiSteadyLiquidImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int nMoities, const RowIterator& jac, LinearBufferAllocator workSpace) const { } virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, active const* y, Eigen::Map fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) { return 0; } - virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, Eigen::Map fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) { return 0; diff --git a/src/libcadet/model/reaction/ReactionModelBase.hpp b/src/libcadet/model/reaction/ReactionModelBase.hpp index 453ee71e4..6e2edfd58 100644 --- a/src/libcadet/model/reaction/ReactionModelBase.hpp +++ b/src/libcadet/model/reaction/ReactionModelBase.hpp @@ -167,27 +167,27 @@ class DynamicReactionModelBase : public IDynamicReactionModel jacobianLiquidImpl(t, secIdx, colPos, y, factor, jac, workSpace); \ } \ virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, \ - double factor, linalg::BandMatrix::RowIterator jac, LinearBufferAllocator workSpace) const \ + int nMoities, linalg::BandMatrix::RowIterator jac, LinearBufferAllocator workSpace) const \ { \ - jacobianQuasiSteadyLiquidImpl(t, secIdx, colPos, y, factor, jac, workSpace); \ + jacobianQuasiSteadyLiquidImpl(t, secIdx, colPos, y, nMoities, jac, workSpace); \ } \ \ virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, \ - double factor, linalg::DenseBandedRowIterator jac, LinearBufferAllocator workSpace) const \ + int nMoities, linalg::DenseBandedRowIterator jac, LinearBufferAllocator workSpace) const \ { \ - jacobianQuasiSteadyLiquidImpl(t, secIdx, colPos, y, factor, jac, workSpace); \ + jacobianQuasiSteadyLiquidImpl(t, secIdx, colPos, y, nMoities, jac, workSpace); \ } \ \ virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, \ - double factor, linalg::BandedSparseRowIterator jac, LinearBufferAllocator workSpace) const \ + int nMoities, linalg::BandedSparseRowIterator jac, LinearBufferAllocator workSpace) const \ { \ - jacobianQuasiSteadyLiquidImpl(t, secIdx, colPos, y, factor, jac, workSpace); \ + jacobianQuasiSteadyLiquidImpl(t, secIdx, colPos, y, nMoities, jac, workSpace); \ } \ \ virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, \ - double factor, linalg::BandedEigenSparseRowIterator jac, LinearBufferAllocator workSpace) const \ + int nMoities, linalg::BandedEigenSparseRowIterator jac, LinearBufferAllocator workSpace) const \ { \ - jacobianQuasiSteadyLiquidImpl(t, secIdx, colPos, y, factor, jac, workSpace); \ + jacobianQuasiSteadyLiquidImpl(t, secIdx, colPos, y, nMoities, jac, workSpace); \ } \ \ virtual void analyticJacobianCombinedAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* yLiquid, double const* ySolid, \ From 0b3aaa47a4a93cac7a80cc611d32350cdf6b2f95 Mon Sep 17 00:00:00 2001 From: "a.berger" Date: Fri, 31 Jan 2025 16:35:48 +0100 Subject: [PATCH 04/11] edded leanConsistenInitalisationTimeDerivative, fixed some bugs --- src/libcadet/model/ReactionModel.hpp | 7 +- src/libcadet/model/StirredTankModel.cpp | 250 ++++++++++++++---- src/libcadet/model/StirredTankModel.hpp | 6 +- src/libcadet/model/UnitOperation.hpp | 1 + .../reaction/CrystallizationReaction.cpp | 4 +- src/libcadet/model/reaction/DummyReaction.cpp | 8 +- .../model/reaction/MassActionLawReaction.cpp | 53 +++- .../reaction/MichaelisMentenReaction.cpp | 4 +- .../model/reaction/ReactionModelBase.hpp | 18 +- 9 files changed, 267 insertions(+), 84 deletions(-) diff --git a/src/libcadet/model/ReactionModel.hpp b/src/libcadet/model/ReactionModel.hpp index 6c57d22cd..e89d1b564 100644 --- a/src/libcadet/model/ReactionModel.hpp +++ b/src/libcadet/model/ReactionModel.hpp @@ -246,12 +246,13 @@ class IDynamicReactionModel virtual void analyticJacobianLiquidAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, linalg::DenseBandedRowIterator jac, LinearBufferAllocator workSpace) const = 0; virtual void analyticJacobianLiquidAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, linalg::BandedSparseRowIterator jac, LinearBufferAllocator workSpace) const = 0; - virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int nMoities, linalg::BandMatrix::RowIterator jac, LinearBufferAllocator workSpace) const = 0; - virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int nMoities, linalg::DenseBandedRowIterator jac, LinearBufferAllocator workSpace) const = 0; - virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int nMoities, linalg::BandedSparseRowIterator jac, LinearBufferAllocator workSpace) const = 0; + virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, int reaction, linalg::BandMatrix::RowIterator jac, LinearBufferAllocator workSpace) const = 0; + virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, int reaction, linalg::DenseBandedRowIterator jac, LinearBufferAllocator workSpace) const = 0; + virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, int reaction, linalg::BandedSparseRowIterator jac, LinearBufferAllocator workSpace) const = 0; virtual void fillConservedMoietiesBulk(Eigen::MatrixXd& M, std::vector& QSReaction, std::vector& _QsCompBulk) = 0; + virtual void timeDerivativeQuasiStationaryReaction(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double* dReacDt, LinearBufferAllocator workSpace) = 0; #ifdef ENABLE_DG virtual void analyticJacobianLiquidAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, linalg::BandedEigenSparseRowIterator jac, LinearBufferAllocator workSpace) const = 0; diff --git a/src/libcadet/model/StirredTankModel.cpp b/src/libcadet/model/StirredTankModel.cpp index 414acf25d..c03e36adc 100644 --- a/src/libcadet/model/StirredTankModel.cpp +++ b/src/libcadet/model/StirredTankModel.cpp @@ -250,7 +250,7 @@ bool CSTRModel::configureModelDiscretization(IParameterProvider& paramProvider, _temp = new active[_nComp]; _temp2 = new double[_nQsReacBulk]; - _nMoitiesBulk = (_nComp + _totalBound) - _nQsReacBulk; + _nMoitiesBulk = 1; // Achtung nurs fuers debugging _MconvMoityBulk = Eigen::MatrixXd::Zero(_nMoitiesBulk, _nComp); // matrix for conserved moities } else @@ -451,6 +451,25 @@ unsigned int CSTRModel::threadLocalMemorySize() const CADET_NOEXCEPT lms.add(_totalBound); lms.commit(); + if (_nQsReacBulk > 0) + { + // Memory for leanConsistentInitialTimeDerivative() + lms.add(_nComp + _totalBound); + lms.add(_nComp + _totalBound); + lms.add(_nComp); + lms.add(_nComp + _totalBound); + lms.add(numDofs()); + lms.add(_nonlinearSolver->workspaceSize(_nComp + _totalBound) * sizeof(double)); + lms.addBlock(resImplSize); + lms.commit(); + + // Memory for leanConsistentInitialState() + lms.add(_totalBound); + lms.add(_nComp + 1); + lms.add(_nComp ); + lms.add(_totalBound); + lms.commit(); + } return lms.bufferSize(); } @@ -476,6 +495,8 @@ void CSTRModel::notifyDiscontinuousSectionTransition(double t, unsigned int secI { _curFlowRateFilter = _flowRateFilter[0]; } + + _curSecIdx = secIdx; } void CSTRModel::reportSolution(ISolutionRecorder& recorder, double const* const solution) const @@ -1078,6 +1099,9 @@ void CSTRModel::leanConsistentInitialState(const SimulationTime& simTime, double double* const c = vecStateY + _nComp; const double vLiquid = c[_nComp]; const double vSolid = static_cast(_constSolidVolume); + + if(_nQsReacBulk == 0) + return; // Check if volume is 0 if (vLiquid == 0.0) { @@ -1137,14 +1161,33 @@ void CSTRModel::leanConsistentInitialState(const SimulationTime& simTime, double } if (_nQsReacBulk > 0) { + //Todo check meomory -> something is wring with fullX LinearBufferAllocator tlmAlloc = threadLocalMem.get(); BufferedArray qsMask = tlmAlloc.array(_nComp); + // Mark components with quasi-stationary reactions partition qsMask.copyFromVector(_QsCompBulk); const linalg::ConstMaskArray mask{ static_cast(qsMask), static_cast(_nComp) }; const int probSize = linalg::numMaskActive(mask); + // Extract initial values from current state + BufferedArray solution = tlmAlloc.array(probSize); + linalg::selectVectorSubset(c, mask, static_cast(solution)); + + // Save values of conserved moieties; + const unsigned int numActiveComp = numMaskActive(mask, _nComp); + BufferedArray conservedQuants = tlmAlloc.array(_nMoitiesBulk); + + // Calculate conserved quantities for the inital state + for (unsigned int MoityIdx = 0; MoityIdx < _nMoitiesBulk; ++MoityIdx) + { + double dotProduct = 0.0; + for (unsigned int i = 0; i < _MconvMoityBulk.cols(); ++i) + dotProduct += _MconvMoityBulk(MoityIdx, i) * (c[i]); + conservedQuants[MoityIdx] = dotProduct; + } + linalg::DenseMatrixView jacobianMatrix(_jacFact.data(), _jacFact.pivotData(), probSize, probSize); BufferedArray baFullX = tlmAlloc.array(numDofs()); double* const fullX = static_cast(baFullX); @@ -1155,9 +1198,6 @@ void CSTRModel::leanConsistentInitialState(const SimulationTime& simTime, double BufferedArray baNonlinMem = tlmAlloc.array(_nonlinearSolver->workspaceSize(probSize)); double* const nonlinMem = static_cast(baNonlinMem); - // Extract initial values from current state - BufferedArray solution = tlmAlloc.array(probSize); - linalg::selectVectorSubset(c, mask, static_cast(solution)); std::function jacFunc; if (adJac.adY && adJac.adRes) @@ -1209,8 +1249,7 @@ void CSTRModel::leanConsistentInitialState(const SimulationTime& simTime, double { // Prepare input vector by overwriting masked items std::copy_n(c, mask.len, fullX + _nComp); - fullX[(2 * _nComp + _totalBound)] = c[(_nComp + _totalBound)]; - + linalg::applyVectorSubset(x, mask, fullX + _nComp); // Call residual function to initialize jacobian @@ -1219,18 +1258,27 @@ void CSTRModel::leanConsistentInitialState(const SimulationTime& simTime, double // Extract Jacobian from full Jacobian mat.setAll(0.0); linalg::copyMatrixSubset(_jac, mask, mask, mat); - + // Replace upper part with conservation relations + mat.submatrixSetAll(0.0, 0, 0, _nMoitiesBulk, probSize); + unsigned int rIdx = 0; + for (unsigned int MoityIdx = 0; MoityIdx < _nMoitiesBulk; ++MoityIdx) + { + double dotProduct = 0.0; + for (unsigned int i = 0; i < _MconvMoityBulk.cols(); ++i) // hier Optimierung durch Vermeidung von 0 Zeilen in MconvMoityBulk + mat.native(rIdx,i) = _MconvMoityBulk(MoityIdx, i); + rIdx++; + } return true; }; } - + + std::copy_n(vecStateY, _nComp, fullX); + fullX[(2 * _nComp + _totalBound)] = c[(_nComp + _totalBound)]; _nonlinearSolver->solve( [&](double const* const x, double* const r) { // Prepare input vector by overwriting masked items std::copy_n(c, mask.len, fullX + _nComp); - fullX[(2*_nComp + _totalBound)] = c[(_nComp + _totalBound)]; - linalg::applyVectorSubset(x, mask, fullX + _nComp); // Call residual function @@ -1238,12 +1286,26 @@ void CSTRModel::leanConsistentInitialState(const SimulationTime& simTime, double // Extract values from residual linalg::selectVectorSubset(fullResidual + _nComp, mask, r); + std::fill_n(r, _nMoitiesBulk, 0.0); + + int rIdx = 0; + for (unsigned int MoityIdx = 0; MoityIdx < _nMoitiesBulk; ++MoityIdx) + { + double dotProduct = 0.0; + for (unsigned int j = 0; j < _MconvMoityBulk.cols(); ++j) // hier Optimierung durch Vermeidung von 0 Zeilen in MconvMoityBulk + dotProduct += _MconvMoityBulk(MoityIdx, j) * x[j]; + r[rIdx] = dotProduct - conservedQuants[rIdx]; + rIdx++; + } std::cout << "Residual: " << std::endl; - for (unsigned int i = 0; i < probSize; ++i){ + for (unsigned int i = 0; i < probSize; ++i) { std::cout << r[i] << std::endl; } - + std::cout << "Solution: " << std::endl; + for (unsigned int i = 0; i < _nComp; ++i) + std::cout << x[i] << std::endl; + return true; }, jacFunc, errorTol, static_cast(solution), nonlinMem, jacobianMatrix, probSize); @@ -1253,8 +1315,7 @@ void CSTRModel::leanConsistentInitialState(const SimulationTime& simTime, double std::cout << "Solution: " << std::endl; for (unsigned int i = 0; i < _nComp; ++i) std::cout << solution[i] << std::endl; - - int a = 1; + // Refine / correct solution } @@ -1273,12 +1334,15 @@ void CSTRModel::leanConsistentInitialTimeDerivative(double t, double const* cons const double flowIn = static_cast(_flowRateIn); const double flowOut = static_cast(_flowRateOut); + LinearBufferAllocator tlmAlloc = threadLocalMem.get(); + // Volume: \dot{V} = F_{in} - F_{out} - F_{filter} const double vDot = flowIn - flowOut - static_cast(_curFlowRateFilter); cDot[_nComp] = vDot; - + + // Check if volume is 0 - if (vLiquid == 0.0) + if (vLiquid == 0.0 && _nQsReacBulk == 0) { // We have the equation // V^l * \dot{c}_i + \dot{V}^l * c_i + V^s * [sum_j sum_m d_j q_{i,m}]} = c_{in,i} * F_in - c_i * F_out @@ -1342,34 +1406,98 @@ void CSTRModel::leanConsistentInitialTimeDerivative(double t, double const* cons } } } - else + if (_nQsReacBulk > 0) { - // Concentrations: V^l * \dot{c} + \dot{V}^l * c + V^s * [sum_j sum_m d_j \dot{q}_{j,m}] = c_in * F_in + c * F_out - // <=> V^l * \dot{c} = c_in * F_in + c * F_out - V^s * [sum_j sum_m d_j \dot{q}_{j,m}] - \dot{V}^l * c - // = -res - V^s * [sum_j sum_m d_j \dot{q}_{j,m}] - \dot{V}^l * c - // => \dot{c} = (-res - V^s * [sum_j sum_m d_j \dot{q}_{j,m}] - \dot{V}^l * c) / V^l - for (unsigned int i = 0; i < _nComp; ++i) + if (vLiquid == 0) + return; // todo not implmenteted yet + + for (unsigned int i = 0; i < numDofs(); ++i) + vecStateYdot[i] = -vecStateYdot[i] - res[i]; + + + // Assemble time derivative Jacobian + _jacFact.setAll(0.0); + addTimeDerivativeJacobian(t, 1.0, ConstSimulationState{ vecStateY, nullptr }, _jacFact); + + LinearBufferAllocator tlmAlloc = threadLocalMem.get(); // -> todo hier aufpassen ob ich den mem auch benutzen darf + + // Overwrite rows corresponding to algebraic equations with the Jacobian and set right hand side to 0 + BufferedArray dReacDt = tlmAlloc.array(_nComp + 1); + + BufferedArray qsMask = tlmAlloc.array(_nComp); + qsMask.copyFromVector(_QsCompBulk); + + const linalg::ConstMaskArray mask{ static_cast(qsMask), static_cast(_nComp) }; + + // Obtain derivative of fluxes wrt. time + std::fill_n(static_cast(dReacDt), _nComp + 1, 0.0); + + _dynReactionBulk->timeDerivativeQuasiStationaryReaction(t, _curSecIdx, ColumnPosition{ 0.0, 0.0, 0.0 }, c, static_cast(dReacDt), tlmAlloc); + + // Copy row from original Jacobian and set right hand side + double* const cShellDot = cDot + _nComp; + int idx = 0; + for (unsigned int i = 0; i < _nComp; ++i, ++idx) { - double qSum = 0.0; - double qDotSum = 0.0; - for (unsigned int type = 0; type < _nParType; ++type) + if(i >= _nComp - _nMoitiesBulk) + _jacFact.copyRowFrom(_jac, idx, idx); + cShellDot[i] = -dReacDt[i]; + } + + // Factorize + const bool result = _jacFact.robustFactorize(static_cast(dReacDt)); + if (!result) + { + LOG(Error) << "Factorize() failed"; + } + + // Solve + const bool result2 = _jacFact.robustSolve(vecStateYdot + _nComp, static_cast(dReacDt)); + if (!result2) + { + LOG(Error) << "Solve() failed"; + } + + for (int i = 0; i < _nComp + 1; i++) + std::cout << " resC" << i << ": " << resC[i] << std::endl; + + for (int i = 0; i < _nComp + 1; i++) + std::cout <<" Ydot" << i << ": " << cDot[i] << std::endl; + + + for (int i = 0; i < _nComp + 1; i++) + std::cout << " Y" << i <<": " << c[i] << std::endl; + + } + else + { + // Concentrations: V^l * \dot{c} + \dot{V}^l * c + V^s * [sum_j sum_m d_j \dot{q}_{j,m}] = c_in * F_in + c * F_out + // <=> V^l * \dot{c} = c_in * F_in + c * F_out - V^s * [sum_j sum_m d_j \dot{q}_{j,m}] - \dot{V}^l * c + // = -res - V^s * [sum_j sum_m d_j \dot{q}_{j,m}] - \dot{V}^l * c + // => \dot{c} = (-res - V^s * [sum_j sum_m d_j \dot{q}_{j,m}] - \dot{V}^l * c) / V^l + for (unsigned int i = 0; i < _nComp; ++i) { - double const* const localQ = c + _nComp + _offsetParType[type] + _boundOffset[type * _nComp + i]; - double const* const localQdot = cDot + _nComp + _offsetParType[type] + _boundOffset[type * _nComp + i]; - double qSumType = 0.0; - double qDotSumType = 0.0; - for (unsigned int j = 0; j < _nBound[type * _nComp + i]; ++j) + double qSum = 0.0; + double qDotSum = 0.0; + for (unsigned int type = 0; type < _nParType; ++type) { - qSumType += localQ[j]; - qDotSumType += localQdot[j]; + double const* const localQ = c + _nComp + _offsetParType[type] + _boundOffset[type * _nComp + i]; + double const* const localQdot = cDot + _nComp + _offsetParType[type] + _boundOffset[type * _nComp + i]; + double qSumType = 0.0; + double qDotSumType = 0.0; + for (unsigned int j = 0; j < _nBound[type * _nComp + i]; ++j) + { + qSumType += localQ[j]; + qDotSumType += localQdot[j]; + } + qSum += static_cast(_parTypeVolFrac[type]) * qSumType; + qDotSum += static_cast(_parTypeVolFrac[type]) * qDotSumType; } - qSum += static_cast(_parTypeVolFrac[type]) * qSumType; - qDotSum += static_cast(_parTypeVolFrac[type]) * qDotSumType; - } - cDot[i] = (-resC[i] - vSolid * qDotSum - vLiquid * c[i]) / vLiquid; - } + cDot[i] = (-resC[i] - vSolid * qDotSum - vLiquid * c[i]) / vLiquid; + } } + } void CSTRModel::leanConsistentInitialSensitivity(const SimulationTime& simTime, const ConstSimulationState& simState, @@ -1439,10 +1567,9 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons // Divide by beta and add c_i and dc_i / dt resC[i] = cDot[i] * v + vDot * c[i] + vsolid * qDotSum; } - resC[i] += -flowIn * cIn[i] + flowOut * c[i]; } - + if (wantJac) { _jac.setAll(0.0); @@ -1485,44 +1612,54 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons Eigen::Map> mapResC(resC, _nComp); int MoityIdx = 0; - for (unsigned int state = 0; state < (_nComp - _nQsReacBulk); ++state) + int state = 0; + for (unsigned int comp = 0; comp < _nComp; comp++) { - if (_QsCompBulk[state] == 1) + if (_QsCompBulk[comp] == 1 && MoityIdx < _nMoitiesBulk) { ResidualType dotProduct = 0.0; for (unsigned int i = 0 ; i < _MconvMoityBulk.cols(); ++i) // hier Optimierung durch Vermeidung von 0 Zeilen in MconvMoityBulk { dotProduct += static_cast(_MconvMoityBulk(MoityIdx, i)) * (mapResC[i]); + + if (wantJac && i == state) + { + _jac.native(state, state) *= _MconvMoityBulk(MoityIdx, i); // dF_{ci}/dcj = v_liquidDot + F_out + if (cadet_likely(yDot)) + { + _jac.native(i, _nComp + _totalBound) *= _MconvMoityBulk(MoityIdx, i); // dF/dvliquid = cDot + } + } if (wantJac && state != i) { _jac.native(state, i) = _MconvMoityBulk(MoityIdx, i)*(static_cast(vDot) + static_cast(flowOut)); // dF_{ci}/dcj = v_liquidDot + F_out if (cadet_likely(yDot)) _jac.native(state, _nComp + _totalBound) += _MconvMoityBulk(MoityIdx, i) * cDot[i]; // dF_{ci}/dVl = sum_i M[i,-] * dcidt - } + } } resCMoities[state] = dotProduct; MoityIdx++; + state++; } - else if (_QsCompBulk[state] == 0) + else if (_QsCompBulk[comp] == 0) { resCMoities[state] = v * flux[state]; + state++; } } - int state = (_nComp - _nQsReacBulk); + state = (_nComp - _nQsReacBulk); for (unsigned int qsreac = 0; qsreac < _nQsReacBulk; ++qsreac) { - resCMoities[state] = v * qsflux[qsreac]; + resCMoities[state] = qsflux[qsreac]; if(wantJac) { - _dynReactionBulk->analyticQuasiSteadyJacobianLiquid(t, secIdx, colPos, reinterpret_cast(c), state , _jac.row(state), subAlloc); + _dynReactionBulk->analyticQuasiSteadyJacobianLiquid(t, secIdx, colPos, reinterpret_cast(c), state , qsreac, _jac.row(state), subAlloc); _jac.native(state, _nComp + _totalBound) = 0.0; // dF_{ci}/dvliquid = 0 - } state++; } - mapResC = resCMoities; } else @@ -1646,6 +1783,15 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons // Volume: \dot{V} = F_{in} - F_{out} - F_{filter} res[2 * _nComp + _totalBound] = vDot - flowIn + flowOut + static_cast(_curFlowRateFilter); + + + /*std::cout << "Jacobian: " << std::endl; + for (unsigned int i = 0; i < 4; ++i) + { + for (unsigned int j = 0; j < 4; ++j) + std::cout << _jac.native(i, j) << " "; + std::cout << std::endl; + }*/ return 0; } @@ -2040,16 +2186,16 @@ void CSTRModel::addTimeDerivativeJacobian(double t, double alpha, const ConstSim // Concentrations: \dot{V^l} * c_i + V^l * \dot{c}_i + V^s * [sum_j sum_m d_j \dot{q}_{j,i,m}]) - c_{in,i} * F_in + c_i * F_out == 0 for (unsigned int i = 0; i < _nComp - _nQsReacBulk; ++i) { - mat.native(i, i) += timeV; // dRes / dcDot if (_nQsReacBulk > 0) { - mat.native(i, i) *= static_cast(_MconvMoityBulk(i, i)); // dRes / dcDot + mat.native(i, i) += timeV*static_cast(_MconvMoityBulk(i, i)); // dRes / dcDot for(unsigned int comp = 0; comp < _nComp; ++comp) { if(_QsCompBulk[comp] == 1 && i!=comp) - mat.native(i, comp) = timeV * static_cast(_MconvMoityBulk(i, comp)); // dRes / dcDot - mat.native(i, _nComp + _totalBound) += alpha * c[comp]; // dRes / dVlDot + mat.native(i, comp) += timeV * static_cast(_MconvMoityBulk(i, comp)); // dRes / dcDot + // todo this is wrong + mat.native(i, _nComp + _totalBound) += alpha * static_cast(_MconvMoityBulk(i, comp))* c[comp]; // dRes / dVlDot } if(i >= _nComp - _nQsReacBulk) diff --git a/src/libcadet/model/StirredTankModel.hpp b/src/libcadet/model/StirredTankModel.hpp index 52e94378a..d6bd1c739 100644 --- a/src/libcadet/model/StirredTankModel.hpp +++ b/src/libcadet/model/StirredTankModel.hpp @@ -99,8 +99,9 @@ class CSTRModel : public UnitOperationBase std::vector& vecSensY, std::vector& vecSensYdot, active const* const adRes, util::ThreadLocalStorage& threadLocalMem); virtual void leanConsistentInitialState(const SimulationTime& simTime, double* const vecStateY, const AdJacobianParams& adJac, double errorTol, util::ThreadLocalStorage& threadLocalMem); - virtual void leanConsistentInitialTimeDerivative(double t, double const* const vecStateY, double* const vecStateYdot, double* const res, util::ThreadLocalStorage& threadLocalMem); - + virtual void leanConsistentInitialTimeDerivative(double time, double const* const vecStateY, double* const vecStateYdot, double* const res, util::ThreadLocalStorage& threadLocalMem); + virtual void leanConsistentInitialTimeDerivative(const SimulationTime& simTime, double const* const vecStateY, double* const vecStateYdot, double* const res, util::ThreadLocalStorage& threadLocalMem) {}; + virtual void leanConsistentInitialSensitivity(const SimulationTime& simTime, const ConstSimulationState& simState, std::vector& vecSensY, std::vector& vecSensYdot, active const* const adRes, util::ThreadLocalStorage& threadLocalMem); @@ -158,6 +159,7 @@ class CSTRModel : public UnitOperationBase active _flowRateIn; //!< Volumetric flow rate of incoming stream active _flowRateOut; //!< Volumetric flow rate of drawn outgoing stream active _curFlowRateFilter; //!< Current volumetric flow rate of liquid outtake stream for this section + unsigned int _curSecIdx; //!< Current section index std::vector _flowRateFilter; //!< Volumetric flow rate of liquid outtake stream std::vector _parTypeVolFrac; //!< Volume fraction of each particle type diff --git a/src/libcadet/model/UnitOperation.hpp b/src/libcadet/model/UnitOperation.hpp index 96e803476..6919af80c 100644 --- a/src/libcadet/model/UnitOperation.hpp +++ b/src/libcadet/model/UnitOperation.hpp @@ -498,6 +498,7 @@ class IUnitOperation : public IModel * @param [in] res On entry, residual without taking time derivatives into account. The data is overwritten during execution of the function. */ virtual void leanConsistentInitialTimeDerivative(double t, double const* const vecStateY, double* const vecStateYdot, double* const res, util::ThreadLocalStorage& threadLocalMem) = 0; + virtual void leanConsistentInitialTimeDerivative(const SimulationTime& simTime, double const* const vecStateY, double* const vecStateYdot, double* const res, util::ThreadLocalStorage& threadLocalMem) = 0; /** * @brief Computes consistent initial conditions for all sensitivity subsystems diff --git a/src/libcadet/model/reaction/CrystallizationReaction.cpp b/src/libcadet/model/reaction/CrystallizationReaction.cpp index 6f376891b..7d4522a1d 100644 --- a/src/libcadet/model/reaction/CrystallizationReaction.cpp +++ b/src/libcadet/model/reaction/CrystallizationReaction.cpp @@ -200,7 +200,7 @@ namespace cadet void fillConservedMoietiesBulk(Eigen::MatrixXd& M, std::vector& QSReaction, std::vector& QSComponent) {} template - void jacobianQuasiSteadyLiquidImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int nMoities, const RowIterator& jac, LinearBufferAllocator workSpace) const { } + void jacobianQuasiSteadyLiquidImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, int reaction, const RowIterator& jac, LinearBufferAllocator workSpace) const { } virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, active const* y, Eigen::Map fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) { @@ -210,7 +210,7 @@ namespace cadet Eigen::Map fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) { return 0; } - + virtual void timeDerivativeQuasiStationaryReaction(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double* dReacDt, LinearBufferAllocator workSpace){ } virtual bool configure(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx) { readScalarParameterOrArray(_bins, paramProvider, "CRY_BINS", 1); diff --git a/src/libcadet/model/reaction/DummyReaction.cpp b/src/libcadet/model/reaction/DummyReaction.cpp index 17e8b8a83..48ebc03c9 100644 --- a/src/libcadet/model/reaction/DummyReaction.cpp +++ b/src/libcadet/model/reaction/DummyReaction.cpp @@ -78,10 +78,12 @@ class DummyDynamicReaction : public IDynamicReactionModel virtual void analyticJacobianLiquidAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, linalg::DenseBandedRowIterator jac, LinearBufferAllocator workSpace) const { } virtual void analyticJacobianLiquidAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, linalg::BandedSparseRowIterator jac, LinearBufferAllocator workSpace) const { } - virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int nMoities , linalg::BandMatrix::RowIterator jac, LinearBufferAllocator workSpace) const { } - virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int nMoities, linalg::DenseBandedRowIterator jac, LinearBufferAllocator workSpace) const { } - virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int nMoities, linalg::BandedSparseRowIterator jac, LinearBufferAllocator workSpace) const { } + virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, int reaction, linalg::BandMatrix::RowIterator jac, LinearBufferAllocator workSpace) const { } + virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, int reaction, linalg::DenseBandedRowIterator jac, LinearBufferAllocator workSpace) const { } + virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, int reaction, linalg::BandedSparseRowIterator jac, LinearBufferAllocator workSpace) const { } + virtual void timeDerivativeQuasiStationaryReaction(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double* dReacDt, LinearBufferAllocator workSpace) { } + #ifdef ENABLE_DG virtual void analyticJacobianLiquidAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, linalg::BandedEigenSparseRowIterator jac, LinearBufferAllocator workSpace) const { } #endif diff --git a/src/libcadet/model/reaction/MassActionLawReaction.cpp b/src/libcadet/model/reaction/MassActionLawReaction.cpp index a03384787..8855d516f 100644 --- a/src/libcadet/model/reaction/MassActionLawReaction.cpp +++ b/src/libcadet/model/reaction/MassActionLawReaction.cpp @@ -400,16 +400,31 @@ class MassActionLawReactionBase : public DynamicReactionModelBase Eigen::MatrixXd QSS(_stoichiometryBulk.rows(), nQS); - int rowIndex = 0; + // Fülle die Spalten basierend auf _reactionQuasistationarity - for (int i = 0; i < _stoichiometryBulk.rows(); ++i) + /*for (int i = 0; i < _stoichiometryBulk.rows(); ++i) { const double* rowData = reinterpret_cast(_stoichiometryBulk.data() + i * _stoichiometryBulk.columns()); Eigen::Map Srow(rowData, _stoichiometryBulk.columns()); QSS.row(rowIndex++) = Srow; } + */ + for (int i = 0; i < _stoichiometryBulk.rows(); ++i) + { + for(int j = 0; j < _stoichiometryBulk.columns(); j++) + { + int rowIndex = 0; + for (int j = 0; j < _stoichiometryBulk.columns(); j++) + { + if (mapQSReac[j] == 0) + continue; + QSS(i, rowIndex) = static_cast(_stoichiometryBulk.native(i, j)); + rowIndex++; + } + } + } //Remove zero rows from QSS int nQScomp = 0; // Number of quasi stationary active components for (int i = 0; i < QSS.rows(); ++i) @@ -424,7 +439,7 @@ class MassActionLawReactionBase : public DynamicReactionModelBase } } // Redimensioniere QSS und kopiere nur die nicht-null Zeilen - if (_nComp - nQScomp < 0.0 ) // if kinetic components exists we can resize QSS + if (_nComp - nQScomp > 0 ) // if kinetic components exists we can resize QSS { Eigen::MatrixXd QSSCompressed(nQScomp, QSS.cols()); for (std::size_t i = 0; i < _QsCompBulk.size(); ++i) @@ -441,13 +456,13 @@ class MassActionLawReactionBase : public DynamicReactionModelBase int rang = QSS.fullPivLu().rank(); if (rang != nQS) { - throw std::runtime_error("The matrix is not full rank"); + throw std::runtime_error("Calculation Conserved Moities Matrix: The Stoichemetric Matrix is sigular"); } //3. Calculate the null space of the matrix Eigen::MatrixXd leftZeroSpace = QSS.transpose().fullPivLu().kernel().transpose(); - if (_nComp - nQScomp < 1e-10) // if there are no kinetic components we can return the zero matrix + if (_nComp - nQScomp < 0) // if there are no kinetic components we can return the zero matrix { M = leftZeroSpace; } @@ -933,28 +948,27 @@ class MassActionLawReactionBase : public DynamicReactionModelBase } template - void jacobianQuasiSteadyLiquidImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, const RowIterator& jac, LinearBufferAllocator workSpace) const + void jacobianQuasiSteadyLiquidImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, int reaction ,const RowIterator& jac, LinearBufferAllocator workSpace) const { typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); BufferedArray fluxes = workSpace.array(2 * _nComp); double* const fluxGradFwd = static_cast(fluxes); double* const fluxGradBwd = fluxGradFwd + _nComp; - for (int r = 0; r < _stoichiometryBulk.columns(); ++r) // reactions - { + // Calculate gradients of forward and backward fluxes - fluxGradLiquid(fluxGradFwd, r, _nComp, static_cast(p->kFwdBulk[r]), _expBulkFwd, y); - fluxGradLiquid(fluxGradBwd, r, _nComp, static_cast(p->kBwdBulk[r]), _expBulkBwd, y); + fluxGradLiquid(fluxGradFwd, reaction, _nComp, static_cast(p->kFwdBulk[reaction]), _expBulkFwd, y); + fluxGradLiquid(fluxGradBwd, reaction, _nComp, static_cast(p->kBwdBulk[reaction]), _expBulkBwd, y); // Add gradients to Jacobian RowIterator curJac = jac; // right row iterator - const double colFactor = static_cast(_stoichiometryBulk.native(state, r)); + const double colFactor = static_cast(_stoichiometryBulk.native(state, reaction)); for (int col = 0; col < _nComp; ++col) { curJac[col - state] = colFactor * (fluxGradFwd[col] - fluxGradBwd[col]); } - } + } template @@ -1001,6 +1015,21 @@ class MassActionLawReactionBase : public DynamicReactionModelBase } } } + + virtual void timeDerivativeQuasiStationaryReaction(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double* dY, LinearBufferAllocator workSpace) + { + typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); + + for (int r = 0; r < _stoichiometryBulk.columns(); ++r) + { + double kFwdBulk_r = static_cast(p->kFwdBulk[r]); + double kBwdBulk_r = static_cast(p->kBwdBulk[r]); + + double flow = singleFlux(r, y, kFwdBulk_r, kBwdBulk_r); + for (int c = 0; c < _nComp; ++c) + dY[c] += static_cast(_stoichiometryBulk.native(c, r)) * flow; + } + } }; typedef MassActionLawReactionBase MassActionLawReaction; diff --git a/src/libcadet/model/reaction/MichaelisMentenReaction.cpp b/src/libcadet/model/reaction/MichaelisMentenReaction.cpp index b0525c9d7..f8d146aad 100644 --- a/src/libcadet/model/reaction/MichaelisMentenReaction.cpp +++ b/src/libcadet/model/reaction/MichaelisMentenReaction.cpp @@ -155,8 +155,8 @@ class MichaelisMentenReactionBase : public DynamicReactionModelBase void fillConservedMoietiesBulk(Eigen::MatrixXd& M, std::vector& QSReaction, std::vector& QSComponent) {} template - void jacobianQuasiSteadyLiquidImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int nMoities, const RowIterator& jac, LinearBufferAllocator workSpace) const { } - + void jacobianQuasiSteadyLiquidImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, int reaction ,const RowIterator& jac, LinearBufferAllocator workSpace) const { } + virtual void timeDerivativeQuasiStationaryReaction(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double* dReacDt, LinearBufferAllocator workSpace) { } virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, active const* y, Eigen::Map fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) { return 0; diff --git a/src/libcadet/model/reaction/ReactionModelBase.hpp b/src/libcadet/model/reaction/ReactionModelBase.hpp index 6e2edfd58..665e96a84 100644 --- a/src/libcadet/model/reaction/ReactionModelBase.hpp +++ b/src/libcadet/model/reaction/ReactionModelBase.hpp @@ -65,6 +65,8 @@ class DynamicReactionModelBase : public IDynamicReactionModel virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, Eigen::Map fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) = 0; + + virtual void timeDerivativeQuasiStationaryReaction(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double* dReacDt, LinearBufferAllocator workSpace) = 0; protected: int _nComp; //!< Number of components @@ -167,27 +169,27 @@ class DynamicReactionModelBase : public IDynamicReactionModel jacobianLiquidImpl(t, secIdx, colPos, y, factor, jac, workSpace); \ } \ virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, \ - int nMoities, linalg::BandMatrix::RowIterator jac, LinearBufferAllocator workSpace) const \ + int state,int reaction, linalg::BandMatrix::RowIterator jac, LinearBufferAllocator workSpace) const \ { \ - jacobianQuasiSteadyLiquidImpl(t, secIdx, colPos, y, nMoities, jac, workSpace); \ + jacobianQuasiSteadyLiquidImpl(t, secIdx, colPos, y, state, reaction,jac, workSpace); \ } \ \ virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, \ - int nMoities, linalg::DenseBandedRowIterator jac, LinearBufferAllocator workSpace) const \ + int state, int reaction,linalg::DenseBandedRowIterator jac, LinearBufferAllocator workSpace) const \ { \ - jacobianQuasiSteadyLiquidImpl(t, secIdx, colPos, y, nMoities, jac, workSpace); \ + jacobianQuasiSteadyLiquidImpl(t, secIdx, colPos, y, state, reaction,jac, workSpace); \ } \ \ virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, \ - int nMoities, linalg::BandedSparseRowIterator jac, LinearBufferAllocator workSpace) const \ + int state, int reaction,linalg::BandedSparseRowIterator jac, LinearBufferAllocator workSpace) const \ { \ - jacobianQuasiSteadyLiquidImpl(t, secIdx, colPos, y, nMoities, jac, workSpace); \ + jacobianQuasiSteadyLiquidImpl(t, secIdx, colPos, y, state, reaction, jac, workSpace); \ } \ \ virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, \ - int nMoities, linalg::BandedEigenSparseRowIterator jac, LinearBufferAllocator workSpace) const \ + int state,int reaction, linalg::BandedEigenSparseRowIterator jac, LinearBufferAllocator workSpace) const \ { \ - jacobianQuasiSteadyLiquidImpl(t, secIdx, colPos, y, nMoities, jac, workSpace); \ + jacobianQuasiSteadyLiquidImpl(t, secIdx, colPos, y, state, reaction,jac, workSpace); \ } \ \ virtual void analyticJacobianCombinedAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* yLiquid, double const* ySolid, \ From 69f453c096fe17c7be2f8af7af96741cd11a7486 Mon Sep 17 00:00:00 2001 From: "a.berger" Date: Tue, 4 Feb 2025 10:23:42 +0100 Subject: [PATCH 05/11] Jacobian with CM is now independet from jacobian in generell case --- src/libcadet/model/ReactionModel.hpp | 4 + src/libcadet/model/StirredTankModel.cpp | 153 ++++++++++++------ src/libcadet/model/StirredTankModel.hpp | 1 - src/libcadet/model/UnitOperation.hpp | 2 - .../reaction/CrystallizationReaction.cpp | 3 + src/libcadet/model/reaction/DummyReaction.cpp | 4 + .../model/reaction/MassActionLawReaction.cpp | 29 +++- .../reaction/MichaelisMentenReaction.cpp | 4 + .../model/reaction/ReactionModelBase.hpp | 43 +++-- 9 files changed, 175 insertions(+), 68 deletions(-) diff --git a/src/libcadet/model/ReactionModel.hpp b/src/libcadet/model/ReactionModel.hpp index e89d1b564..dcadf7399 100644 --- a/src/libcadet/model/ReactionModel.hpp +++ b/src/libcadet/model/ReactionModel.hpp @@ -249,6 +249,10 @@ class IDynamicReactionModel virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, int reaction, linalg::BandMatrix::RowIterator jac, LinearBufferAllocator workSpace) const = 0; virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, int reaction, linalg::DenseBandedRowIterator jac, LinearBufferAllocator workSpace) const = 0; virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, int reaction, linalg::BandedSparseRowIterator jac, LinearBufferAllocator workSpace) const = 0; + + virtual void analyticJacobianLiquidSingleFluxAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, int reaction, linalg::BandMatrix::RowIterator jac, LinearBufferAllocator workSpace) const = 0; + virtual void analyticJacobianLiquidSingleFluxAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, int reaction, linalg::DenseBandedRowIterator jac, LinearBufferAllocator workSpace) const = 0; + virtual void analyticJacobianLiquidSingleFluxAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, int reaction, linalg::BandedSparseRowIterator jac, LinearBufferAllocator workSpace) const = 0; virtual void fillConservedMoietiesBulk(Eigen::MatrixXd& M, std::vector& QSReaction, std::vector& _QsCompBulk) = 0; diff --git a/src/libcadet/model/StirredTankModel.cpp b/src/libcadet/model/StirredTankModel.cpp index c03e36adc..9bd9688f4 100644 --- a/src/libcadet/model/StirredTankModel.cpp +++ b/src/libcadet/model/StirredTankModel.cpp @@ -1264,9 +1264,13 @@ void CSTRModel::leanConsistentInitialState(const SimulationTime& simTime, double for (unsigned int MoityIdx = 0; MoityIdx < _nMoitiesBulk; ++MoityIdx) { double dotProduct = 0.0; - for (unsigned int i = 0; i < _MconvMoityBulk.cols(); ++i) // hier Optimierung durch Vermeidung von 0 Zeilen in MconvMoityBulk - mat.native(rIdx,i) = _MconvMoityBulk(MoityIdx, i); - rIdx++; + for (unsigned int i = 0; i < _MconvMoityBulk.cols(); ++i) + {// hier Optimierung durch Vermeidung von 0 Zeilen in MconvMoityBulk + if (!mask.mask[i]) + continue; + mat.native(rIdx, i) = _MconvMoityBulk(MoityIdx, i); + rIdx++; + } } return true; }; @@ -1293,9 +1297,13 @@ void CSTRModel::leanConsistentInitialState(const SimulationTime& simTime, double { double dotProduct = 0.0; for (unsigned int j = 0; j < _MconvMoityBulk.cols(); ++j) // hier Optimierung durch Vermeidung von 0 Zeilen in MconvMoityBulk + { + if(!mask.mask[j]) + continue; dotProduct += _MconvMoityBulk(MoityIdx, j) * x[j]; - r[rIdx] = dotProduct - conservedQuants[rIdx]; - rIdx++; + r[rIdx] = dotProduct - conservedQuants[rIdx]; + rIdx++; + } } std::cout << "Residual: " << std::endl; @@ -1570,7 +1578,7 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons resC[i] += -flowIn * cIn[i] + flowOut * c[i]; } - if (wantJac) + if (wantJac && _nQsReacBulk == 0) { _jac.setAll(0.0); @@ -1599,6 +1607,7 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons if (_dynReactionBulk && (_nQsReacBulk > 0)) { + _jac.setAll(0.0); Eigen::Map> resCMoities(reinterpret_cast(_temp), _nComp); resCMoities.setZero(); @@ -1624,15 +1633,15 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons if (wantJac && i == state) { - _jac.native(state, state) *= _MconvMoityBulk(MoityIdx, i); // dF_{ci}/dcj = v_liquidDot + F_out + _jac.native(state, state) += _MconvMoityBulk(MoityIdx, i) * (static_cast(vDot) + static_cast(flowOut)); // dF_{ci}/dcj = v_liquidDot + F_out if (cadet_likely(yDot)) { - _jac.native(i, _nComp + _totalBound) *= _MconvMoityBulk(MoityIdx, i); // dF/dvliquid = cDot + _jac.native(i, _nComp + _totalBound) += _MconvMoityBulk(MoityIdx, i) * cDot[i]; // dF/dvliquid = cDot } } if (wantJac && state != i) { - _jac.native(state, i) = _MconvMoityBulk(MoityIdx, i)*(static_cast(vDot) + static_cast(flowOut)); // dF_{ci}/dcj = v_liquidDot + F_out + _jac.native(state, i) += _MconvMoityBulk(MoityIdx, i)*(static_cast(vDot) + static_cast(flowOut)); // dF_{ci}/dcj = v_liquidDot + F_out if (cadet_likely(yDot)) _jac.native(state, _nComp + _totalBound) += _MconvMoityBulk(MoityIdx, i) * cDot[i]; // dF_{ci}/dVl = sum_i M[i,-] * dcidt } @@ -1644,6 +1653,14 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons else if (_QsCompBulk[comp] == 0) { resCMoities[state] = v * flux[state]; + + + if (wantJac) + { + + _jac.native(state,comp) += (static_cast(vDot) + static_cast(flowOut)); // dF_{ci}/dcj = v_liquidDot + F_out + _dynReactionBulk->analyticJacobianLiquidSingleFluxAdd(t, secIdx, colPos, reinterpret_cast(c), state, comp, _jac.row(state), subAlloc); + } state++; } } @@ -1785,13 +1802,13 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons res[2 * _nComp + _totalBound] = vDot - flowIn + flowOut + static_cast(_curFlowRateFilter); - /*std::cout << "Jacobian: " << std::endl; + std::cout << "Jacobian: " << std::endl; for (unsigned int i = 0; i < 4; ++i) { for (unsigned int j = 0; j < 4; ++j) std::cout << _jac.native(i, j) << " "; std::cout << std::endl; - }*/ + } return 0; } @@ -2174,7 +2191,7 @@ int CSTRModel::linearSolve(double t, double alpha, double tol, double* const rhs template void CSTRModel::addTimeDerivativeJacobian(double t, double alpha, const ConstSimulationState& simState, MatrixType& mat) -{ +{ double const* const c = simState.vecStateY + _nComp; double const* const q = simState.vecStateY + 2 * _nComp; const double v = simState.vecStateY[2 * _nComp + _totalBound]; @@ -2182,70 +2199,102 @@ void CSTRModel::addTimeDerivativeJacobian(double t, double alpha, const ConstSim const double timeVSolid = static_cast(_constSolidVolume) * alpha; // Assemble Jacobian: dRes / dyDot - + // Concentrations: \dot{V^l} * c_i + V^l * \dot{c}_i + V^s * [sum_j sum_m d_j \dot{q}_{j,i,m}]) - c_{in,i} * F_in + c_i * F_out == 0 - for (unsigned int i = 0; i < _nComp - _nQsReacBulk; ++i) + + if (_nQsReacBulk > 0) { - if (_nQsReacBulk > 0) + int MoityIdx = 0; + int state = 0; + for (unsigned int i = 0; i < _nComp; i++) { - mat.native(i, i) += timeV*static_cast(_MconvMoityBulk(i, i)); // dRes / dcDot - for(unsigned int comp = 0; comp < _nComp; ++comp) + if (_QsCompBulk[i] == 1 && MoityIdx < _nMoitiesBulk) { - if(_QsCompBulk[comp] == 1 && i!=comp) - mat.native(i, comp) += timeV * static_cast(_MconvMoityBulk(i, comp)); // dRes / dcDot + + mat.native(i, i) += timeV * (static_cast(_MconvMoityBulk(MoityIdx, MoityIdx))); // dRes / dcDot + for (int comp = 0; comp < _nComp; comp++) + { + if (_QsCompBulk[comp] == 1 && i != comp) + mat.native(i, comp) += timeV * static_cast(_MconvMoityBulk(MoityIdx, comp)); // dRes / dcDot // todo this is wrong - mat.native(i, _nComp + _totalBound) += alpha * static_cast(_MconvMoityBulk(i, comp))* c[comp]; // dRes / dVlDot - } + mat.native(i, _nComp + _totalBound) += alpha * static_cast(_MconvMoityBulk(MoityIdx, comp)) * c[comp]; // dRes / dVlDot + } - if(i >= _nComp - _nQsReacBulk) - break; + MoityIdx++; + state++; + } + else if (_QsCompBulk[i] == 0) + { + mat.native(state, i) += timeV; // dRes / dcDot + state++; + } } - double qSum = 0.0; - for (unsigned int type = 0; type < _nParType; ++type) + } + else + { + + for (unsigned int i = 0; i < _nComp; ++i) { - double const* const qi = q + _offsetParType[type] + _boundOffset[type * _nComp + i]; - const unsigned int localOffset = _nComp + _offsetParType[type] + _boundOffset[type * _nComp + i]; - const double vSolidParVolFrac = timeVSolid * static_cast(_parTypeVolFrac[type]); - for (unsigned int j = 0; j < _nBound[type * _nComp + i]; ++j) + mat.native(i, i) += timeV; // dRes / dcDot + + double qSum = 0.0; + for (unsigned int type = 0; type < _nParType; ++type) { - mat.native(i, localOffset + j) += vSolidParVolFrac; // dRes / dqDot - // + _nComp: Moves over liquid phase components - // + _offsetParType[type]: Moves to particle type - // + _boundOffset[i]: Moves over bound states of previous components - // + j: Moves to current bound state j of component i + double const* const qi = q + _offsetParType[type] + _boundOffset[type * _nComp + i]; + const unsigned int localOffset = _nComp + _offsetParType[type] + _boundOffset[type * _nComp + i]; + const double vSolidParVolFrac = timeVSolid * static_cast(_parTypeVolFrac[type]); + for (unsigned int j = 0; j < _nBound[type * _nComp + i]; ++j) + { + mat.native(i, localOffset + j) += vSolidParVolFrac; // dRes / dqDot + // + _nComp: Moves over liquid phase components + // + _offsetParType[type]: Moves to particle type + // + _boundOffset[i]: Moves over bound states of previous components + // + j: Moves to current bound state j of component i + } } - } - } - // Bound states - unsigned int globalIdx = _nComp; - for (unsigned int type = 0; type < _nParType; ++type) - { - IBindingModel* const binding = _binding[type]; - if (!binding->hasDynamicReactions()) - { - // Skip binding models without dynamic binding fluxes - globalIdx += _strideBound[type]; - continue; + mat.native(i, _nComp + _totalBound) += alpha * c[i]; // dRes / dVlDot } - int const* const qsReaction = binding->reactionQuasiStationarity(); - for (unsigned int idx = 0; idx < _strideBound[type]; ++idx, ++globalIdx) + + // Bound states + unsigned int globalIdx = _nComp; + for (unsigned int type = 0; type < _nParType; ++type) { - // Skip quasi-stationary fluxes - if (qsReaction[idx]) + IBindingModel* const binding = _binding[type]; + if (!binding->hasDynamicReactions()) + { + // Skip binding models without dynamic binding fluxes + globalIdx += _strideBound[type]; continue; + } + + int const* const qsReaction = binding->reactionQuasiStationarity(); + for (unsigned int idx = 0; idx < _strideBound[type]; ++idx, ++globalIdx) + { + // Skip quasi-stationary fluxes + if (qsReaction[idx]) + continue; - mat.native(globalIdx, globalIdx) += alpha; + mat.native(globalIdx, globalIdx) += alpha; + } } } - // Volume: \dot{V} - F_{in} + F_{out} + F_{filter} == 0 mat.native(_nComp + _totalBound, _nComp + _totalBound) += alpha; + + std::cout << "Time derivative Jacobian: " << std::endl; + for (unsigned int i = 0; i < 4; ++i) + { + for (unsigned int j = 0; j < 4; ++j) + std::cout << mat.native(i, j) << " "; + std::cout << std::endl; + } } + /** * @brief Extracts the system Jacobian from AD seed vectors * @param [in] adRes Residual vector of AD datatypes with band compressed seed vectors diff --git a/src/libcadet/model/StirredTankModel.hpp b/src/libcadet/model/StirredTankModel.hpp index d6bd1c739..0b5d05834 100644 --- a/src/libcadet/model/StirredTankModel.hpp +++ b/src/libcadet/model/StirredTankModel.hpp @@ -100,7 +100,6 @@ class CSTRModel : public UnitOperationBase virtual void leanConsistentInitialState(const SimulationTime& simTime, double* const vecStateY, const AdJacobianParams& adJac, double errorTol, util::ThreadLocalStorage& threadLocalMem); virtual void leanConsistentInitialTimeDerivative(double time, double const* const vecStateY, double* const vecStateYdot, double* const res, util::ThreadLocalStorage& threadLocalMem); - virtual void leanConsistentInitialTimeDerivative(const SimulationTime& simTime, double const* const vecStateY, double* const vecStateYdot, double* const res, util::ThreadLocalStorage& threadLocalMem) {}; virtual void leanConsistentInitialSensitivity(const SimulationTime& simTime, const ConstSimulationState& simState, std::vector& vecSensY, std::vector& vecSensYdot, active const* const adRes, util::ThreadLocalStorage& threadLocalMem); diff --git a/src/libcadet/model/UnitOperation.hpp b/src/libcadet/model/UnitOperation.hpp index 6919af80c..fe07406d0 100644 --- a/src/libcadet/model/UnitOperation.hpp +++ b/src/libcadet/model/UnitOperation.hpp @@ -498,8 +498,6 @@ class IUnitOperation : public IModel * @param [in] res On entry, residual without taking time derivatives into account. The data is overwritten during execution of the function. */ virtual void leanConsistentInitialTimeDerivative(double t, double const* const vecStateY, double* const vecStateYdot, double* const res, util::ThreadLocalStorage& threadLocalMem) = 0; - virtual void leanConsistentInitialTimeDerivative(const SimulationTime& simTime, double const* const vecStateY, double* const vecStateYdot, double* const res, util::ThreadLocalStorage& threadLocalMem) = 0; - /** * @brief Computes consistent initial conditions for all sensitivity subsystems * @details Given the DAE \f[ F(t, y, \dot{y}, p) = 0, \f] the corresponding (linear) forward sensitivity diff --git a/src/libcadet/model/reaction/CrystallizationReaction.cpp b/src/libcadet/model/reaction/CrystallizationReaction.cpp index 7d4522a1d..d39e26954 100644 --- a/src/libcadet/model/reaction/CrystallizationReaction.cpp +++ b/src/libcadet/model/reaction/CrystallizationReaction.cpp @@ -202,6 +202,9 @@ namespace cadet template void jacobianQuasiSteadyLiquidImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, int reaction, const RowIterator& jac, LinearBufferAllocator workSpace) const { } + template + void jacobianSingleFluxImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, int reaction, const RowIterator& jac, LinearBufferAllocator workSpace) const { } + virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, active const* y, Eigen::Map fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) { return 0; diff --git a/src/libcadet/model/reaction/DummyReaction.cpp b/src/libcadet/model/reaction/DummyReaction.cpp index 48ebc03c9..0e5ab5133 100644 --- a/src/libcadet/model/reaction/DummyReaction.cpp +++ b/src/libcadet/model/reaction/DummyReaction.cpp @@ -82,6 +82,10 @@ class DummyDynamicReaction : public IDynamicReactionModel virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, int reaction, linalg::DenseBandedRowIterator jac, LinearBufferAllocator workSpace) const { } virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, int reaction, linalg::BandedSparseRowIterator jac, LinearBufferAllocator workSpace) const { } + virtual void analyticJacobianLiquidSingleFluxAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, int reaction, linalg::BandMatrix::RowIterator jac, LinearBufferAllocator workSpace) const { } + virtual void analyticJacobianLiquidSingleFluxAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, int reaction, linalg::DenseBandedRowIterator jac, LinearBufferAllocator workSpace) const { } + virtual void analyticJacobianLiquidSingleFluxAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, int reaction, linalg::BandedSparseRowIterator jac, LinearBufferAllocator workSpace) const { } + virtual void timeDerivativeQuasiStationaryReaction(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double* dReacDt, LinearBufferAllocator workSpace) { } #ifdef ENABLE_DG diff --git a/src/libcadet/model/reaction/MassActionLawReaction.cpp b/src/libcadet/model/reaction/MassActionLawReaction.cpp index 8855d516f..ed00fe13b 100644 --- a/src/libcadet/model/reaction/MassActionLawReaction.cpp +++ b/src/libcadet/model/reaction/MassActionLawReaction.cpp @@ -446,7 +446,7 @@ class MassActionLawReactionBase : public DynamicReactionModelBase { if (_QsCompBulk[i] == 1) { - QSSCompressed.row(i) = QSS.row(_QsCompBulk[i]); + QSSCompressed.row(i) = QSS.row(i); } } QSS.swap(QSSCompressed); // Swap the compressed matrix back into QSS @@ -947,6 +947,29 @@ class MassActionLawReactionBase : public DynamicReactionModelBase } } + template + void jacobianSingleFluxImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, int reaction, const RowIterator& jac, LinearBufferAllocator workSpace) const + { + typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); + + BufferedArray fluxes = workSpace.array(2 * _nComp); + double* const fluxGradFwd = static_cast(fluxes); + double* const fluxGradBwd = fluxGradFwd + _nComp; + for (int r = 0; r < _stoichiometryBulk.columns(); ++r) + { + // Calculate gradients of forward and backward fluxes + fluxGradLiquid(fluxGradFwd, r, _nComp, static_cast(p->kFwdBulk[r]), _expBulkFwd, y); + fluxGradLiquid(fluxGradBwd, r, _nComp, static_cast(p->kBwdBulk[r]), _expBulkBwd, y); + + // Add gradients to Jacobian + RowIterator curJac = jac; + const double colFactor = static_cast(_stoichiometryBulk.native(state, r)) * 1.0; + for (int col = 0; col < _nComp; ++col) + curJac[col - static_cast(state)] += colFactor * (fluxGradFwd[col] - fluxGradBwd[col]); + + } + } + template void jacobianQuasiSteadyLiquidImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, int reaction ,const RowIterator& jac, LinearBufferAllocator workSpace) const { @@ -963,10 +986,10 @@ class MassActionLawReactionBase : public DynamicReactionModelBase // Add gradients to Jacobian RowIterator curJac = jac; // right row iterator - const double colFactor = static_cast(_stoichiometryBulk.native(state, reaction)); + //const double colFactor = static_cast(_stoichiometryBulk.native(state, reaction)); for (int col = 0; col < _nComp; ++col) { - curJac[col - state] = colFactor * (fluxGradFwd[col] - fluxGradBwd[col]); + curJac[col - state] = (fluxGradFwd[col] - fluxGradBwd[col]); } } diff --git a/src/libcadet/model/reaction/MichaelisMentenReaction.cpp b/src/libcadet/model/reaction/MichaelisMentenReaction.cpp index f8d146aad..f1517863b 100644 --- a/src/libcadet/model/reaction/MichaelisMentenReaction.cpp +++ b/src/libcadet/model/reaction/MichaelisMentenReaction.cpp @@ -156,6 +156,10 @@ class MichaelisMentenReactionBase : public DynamicReactionModelBase template void jacobianQuasiSteadyLiquidImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, int reaction ,const RowIterator& jac, LinearBufferAllocator workSpace) const { } + + template + void jacobianSingleFluxImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, int reaction, const RowIterator& jac, LinearBufferAllocator workSpace) const { } + virtual void timeDerivativeQuasiStationaryReaction(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double* dReacDt, LinearBufferAllocator workSpace) { } virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, active const* y, Eigen::Map fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) { diff --git a/src/libcadet/model/reaction/ReactionModelBase.hpp b/src/libcadet/model/reaction/ReactionModelBase.hpp index 665e96a84..420a1d254 100644 --- a/src/libcadet/model/reaction/ReactionModelBase.hpp +++ b/src/libcadet/model/reaction/ReactionModelBase.hpp @@ -169,27 +169,27 @@ class DynamicReactionModelBase : public IDynamicReactionModel jacobianLiquidImpl(t, secIdx, colPos, y, factor, jac, workSpace); \ } \ virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, \ - int state,int reaction, linalg::BandMatrix::RowIterator jac, LinearBufferAllocator workSpace) const \ + int state,int reaction, linalg::BandMatrix::RowIterator jac, LinearBufferAllocator workSpace) const \ { \ - jacobianQuasiSteadyLiquidImpl(t, secIdx, colPos, y, state, reaction,jac, workSpace); \ + jacobianQuasiSteadyLiquidImpl(t, secIdx, colPos, y, state, reaction,jac, workSpace); \ } \ \ virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, \ - int state, int reaction,linalg::DenseBandedRowIterator jac, LinearBufferAllocator workSpace) const \ + int state, int reaction,linalg::DenseBandedRowIterator jac, LinearBufferAllocator workSpace) const \ { \ - jacobianQuasiSteadyLiquidImpl(t, secIdx, colPos, y, state, reaction,jac, workSpace); \ + jacobianQuasiSteadyLiquidImpl(t, secIdx, colPos, y, state, reaction,jac, workSpace); \ } \ \ virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, \ - int state, int reaction,linalg::BandedSparseRowIterator jac, LinearBufferAllocator workSpace) const \ + int state, int reaction,linalg::BandedSparseRowIterator jac, LinearBufferAllocator workSpace) const \ { \ - jacobianQuasiSteadyLiquidImpl(t, secIdx, colPos, y, state, reaction, jac, workSpace); \ + jacobianQuasiSteadyLiquidImpl(t, secIdx, colPos, y, state, reaction, jac, workSpace); \ } \ \ virtual void analyticQuasiSteadyJacobianLiquid(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, \ - int state,int reaction, linalg::BandedEigenSparseRowIterator jac, LinearBufferAllocator workSpace) const \ + int state,int reaction, linalg::BandedEigenSparseRowIterator jac, LinearBufferAllocator workSpace) const \ { \ - jacobianQuasiSteadyLiquidImpl(t, secIdx, colPos, y, state, reaction,jac, workSpace); \ + jacobianQuasiSteadyLiquidImpl(t, secIdx, colPos, y, state, reaction,jac, workSpace); \ } \ \ virtual void analyticJacobianCombinedAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* yLiquid, double const* ySolid, \ @@ -214,8 +214,31 @@ class DynamicReactionModelBase : public IDynamicReactionModel double factor, linalg::BandMatrix::RowIterator jacLiquid, linalg::DenseBandedRowIterator jacSolid, LinearBufferAllocator workSpace) const \ { \ jacobianCombinedImpl(t, secIdx, colPos, yLiquid, ySolid, factor, jacLiquid, jacSolid, workSpace); \ - } - + } \ + virtual void analyticJacobianLiquidSingleFluxAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, \ + int state, int reaction, linalg::BandMatrix::RowIterator jac, LinearBufferAllocator workSpace) const \ + { \ + jacobianSingleFluxImpl(t, secIdx, colPos, y, state, reaction, jac, workSpace); \ + } \ + \ + virtual void analyticJacobianLiquidSingleFluxAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, \ + int state, int reaction, linalg::DenseBandedRowIterator jac, LinearBufferAllocator workSpace) const \ + { \ + jacobianSingleFluxImpl(t, secIdx, colPos, y, state, reaction, jac, workSpace); \ + } \ + \ + virtual void analyticJacobianLiquidSingleFluxAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, \ + int state, int reaction, linalg::BandedSparseRowIterator jac, LinearBufferAllocator workSpace) const \ + { \ + jacobianSingleFluxImpl(t, secIdx, colPos, y, state, reaction, jac, workSpace); \ + } \ + /*\ + virtual void analyticJacobianLiquidSingleFluxAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, \ + int state, int reaction, linalg::BandedEigenSparseRowIterator jac, LinearBufferAllocator workSpace) const \ + { \ + jacobianSingleFluxImpl(t, secIdx, colPos, y, state, reaction, jac, workSpace); \ + } \ + \*/ #else #define CADET_DYNAMICREACTIONMODEL_BOILERPLATE \ virtual int residualLiquidAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, active const* y, \ From 98dd395f58d70c79208c7f3e5fbfc8bf38cbe27c Mon Sep 17 00:00:00 2001 From: "a.berger" Date: Tue, 4 Feb 2025 16:53:22 +0100 Subject: [PATCH 06/11] first version of interface --- src/libcadet/model/ReactionModel.hpp | 7 +- src/libcadet/model/StirredTankModel.cpp | 66 +++++++++---------- src/libcadet/model/StirredTankModel.hpp | 6 +- .../reaction/CrystallizationReaction.cpp | 7 +- src/libcadet/model/reaction/DummyReaction.cpp | 8 +-- .../model/reaction/MassActionLawReaction.cpp | 50 ++++++++++---- .../reaction/MichaelisMentenReaction.cpp | 8 +-- .../model/reaction/ReactionModelBase.hpp | 8 +-- 8 files changed, 93 insertions(+), 67 deletions(-) diff --git a/src/libcadet/model/ReactionModel.hpp b/src/libcadet/model/ReactionModel.hpp index dcadf7399..4fab0278b 100644 --- a/src/libcadet/model/ReactionModel.hpp +++ b/src/libcadet/model/ReactionModel.hpp @@ -215,10 +215,10 @@ class IDynamicReactionModel virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, active const* y, - Eigen::Map fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) = 0; + Eigen::Map fluxes, int const* mapQSReac, LinearBufferAllocator workSpace) = 0; virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, - Eigen::Map fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) = 0; + Eigen::Map fluxes, int const* mapQSReac, LinearBufferAllocator workSpace) = 0; /** @@ -254,10 +254,11 @@ class IDynamicReactionModel virtual void analyticJacobianLiquidSingleFluxAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, int reaction, linalg::DenseBandedRowIterator jac, LinearBufferAllocator workSpace) const = 0; virtual void analyticJacobianLiquidSingleFluxAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, int reaction, linalg::BandedSparseRowIterator jac, LinearBufferAllocator workSpace) const = 0; - virtual void fillConservedMoietiesBulk(Eigen::MatrixXd& M, std::vector& QSReaction, std::vector& _QsCompBulk) = 0; + virtual void fillConservedMoietiesBulk(Eigen::MatrixXd& M, unsigned int & QSReaction, std::vector& _QsCompBulk) = 0; virtual void timeDerivativeQuasiStationaryReaction(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double* dReacDt, LinearBufferAllocator workSpace) = 0; + virtual int const* reactionQuasiStationarity() const = 0; #ifdef ENABLE_DG virtual void analyticJacobianLiquidAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, linalg::BandedEigenSparseRowIterator jac, LinearBufferAllocator workSpace) const = 0; diff --git a/src/libcadet/model/StirredTankModel.cpp b/src/libcadet/model/StirredTankModel.cpp index 9bd9688f4..a446d42bb 100644 --- a/src/libcadet/model/StirredTankModel.cpp +++ b/src/libcadet/model/StirredTankModel.cpp @@ -242,26 +242,16 @@ bool CSTRModel::configureModelDiscretization(IParameterProvider& paramProvider, if (_dynReactionBulk->usesParamProviderInDiscretizationConfig()) paramProvider.popScope(); - if (true)//(paramProvider.exists("QUASI_STATIONARY_REACTION_BULK")) - { - //_qsReacBulk = paramProvider.getIntArray("QUASI_STATIONARY_REACTION_BULK"); - _qsReacBulk = {1}; - _nQsReacBulk = _qsReacBulk.size(); + + _temp = new active[_nComp]; - _temp2 = new double[_nQsReacBulk]; + _temp2 = new double[_nComp]; - _nMoitiesBulk = 1; // Achtung nurs fuers debugging - _MconvMoityBulk = Eigen::MatrixXd::Zero(_nMoitiesBulk, _nComp); // matrix for conserved moities - } - else - { + _MconvMoityBulk = Eigen::MatrixXd::Zero(0,0); // matrix for conserved moities _QsCompBulk.clear(); - _qsReacBulk.clear(); _nMoitiesBulk = 0; - _nQsReacBulk = 0; - _MconvMoityBulk = Eigen::MatrixXd::Zero(0, 0); - - } + _nqsReactionBulk = 0; + _qsReactionBulk = nullptr; } _dynReaction = std::vector(_nParType, nullptr); @@ -384,9 +374,15 @@ bool CSTRModel::configure(IParameterProvider& paramProvider) dynReactionConfSuccess = _dynReactionBulk->configure(paramProvider, _unitOpIdx, ParTypeIndep); paramProvider.popScope(); - if (true)//(paramProvider.exists("QUASI_STATIONARY_REACTION_BULK")) - { - _dynReactionBulk->fillConservedMoietiesBulk(_MconvMoityBulk, _qsReacBulk, _QsCompBulk); // fill conserved moities matrix + // Get quasi stationary reactions information vector + _qsReactionBulk = _dynReactionBulk->reactionQuasiStationarity(); + + if (_qsReactionBulk != nullptr) + { + // resize Conserved moities matrix + + _dynReactionBulk->fillConservedMoietiesBulk(_MconvMoityBulk, _nqsReactionBulk, _QsCompBulk); // fill conserved moities matrix + _nMoitiesBulk = _MconvMoityBulk.rows(); } } @@ -451,7 +447,7 @@ unsigned int CSTRModel::threadLocalMemorySize() const CADET_NOEXCEPT lms.add(_totalBound); lms.commit(); - if (_nQsReacBulk > 0) + if (_nqsReactionBulk > 0) { // Memory for leanConsistentInitialTimeDerivative() lms.add(_nComp + _totalBound); @@ -1100,7 +1096,7 @@ void CSTRModel::leanConsistentInitialState(const SimulationTime& simTime, double const double vLiquid = c[_nComp]; const double vSolid = static_cast(_constSolidVolume); - if(_nQsReacBulk == 0) + if(_nqsReactionBulk == 0) return; // Check if volume is 0 if (vLiquid == 0.0) @@ -1159,7 +1155,7 @@ void CSTRModel::leanConsistentInitialState(const SimulationTime& simTime, double } } } - if (_nQsReacBulk > 0) + if (_nqsReactionBulk > 0) { //Todo check meomory -> something is wring with fullX LinearBufferAllocator tlmAlloc = threadLocalMem.get(); @@ -1350,7 +1346,7 @@ void CSTRModel::leanConsistentInitialTimeDerivative(double t, double const* cons // Check if volume is 0 - if (vLiquid == 0.0 && _nQsReacBulk == 0) + if (vLiquid == 0.0 && _nqsReactionBulk == 0) { // We have the equation // V^l * \dot{c}_i + \dot{V}^l * c_i + V^s * [sum_j sum_m d_j q_{i,m}]} = c_{in,i} * F_in - c_i * F_out @@ -1414,7 +1410,7 @@ void CSTRModel::leanConsistentInitialTimeDerivative(double t, double const* cons } } } - if (_nQsReacBulk > 0) + if (_nqsReactionBulk > 0) { if (vLiquid == 0) return; // todo not implmenteted yet @@ -1578,7 +1574,7 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons resC[i] += -flowIn * cIn[i] + flowOut * c[i]; } - if (wantJac && _nQsReacBulk == 0) + if (wantJac && _nqsReactionBulk == 0) { _jac.setAll(0.0); @@ -1605,18 +1601,18 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons BufferedArray flux = subAlloc.array(_nComp); std::fill_n(static_cast(flux), _nComp, 0.0); - if (_dynReactionBulk && (_nQsReacBulk > 0)) + if (_dynReactionBulk && (_nqsReactionBulk > 0)) { _jac.setAll(0.0); Eigen::Map> resCMoities(reinterpret_cast(_temp), _nComp); resCMoities.setZero(); - Eigen::Map> qsflux(_temp2, _nQsReacBulk); + Eigen::Map> qsflux(_temp2, _nqsReactionBulk); qsflux.setZero(); _dynReactionBulk->residualLiquidAdd(t, secIdx, colPos, c, static_cast(flux), -1.0, subAlloc); - _dynReactionBulk->quasiStationaryFlux(t, secIdx, colPos, c, qsflux, _qsReacBulk, subAlloc); + _dynReactionBulk->quasiStationaryFlux(t, secIdx, colPos, c, qsflux, _qsReactionBulk, subAlloc); Eigen::Map> mapResC(resC, _nComp); @@ -1665,8 +1661,8 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons } } - state = (_nComp - _nQsReacBulk); - for (unsigned int qsreac = 0; qsreac < _nQsReacBulk; ++qsreac) + state = (_nComp - _nqsReactionBulk); + for (unsigned int qsreac = 0; qsreac < _nqsReactionBulk; ++qsreac) { resCMoities[state] = qsflux[qsreac]; @@ -1801,14 +1797,14 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons // Volume: \dot{V} = F_{in} - F_{out} - F_{filter} res[2 * _nComp + _totalBound] = vDot - flowIn + flowOut + static_cast(_curFlowRateFilter); - + /* std::cout << "Jacobian: " << std::endl; for (unsigned int i = 0; i < 4; ++i) { for (unsigned int j = 0; j < 4; ++j) std::cout << _jac.native(i, j) << " "; std::cout << std::endl; - } + }*/ return 0; } @@ -2202,7 +2198,7 @@ void CSTRModel::addTimeDerivativeJacobian(double t, double alpha, const ConstSim // Concentrations: \dot{V^l} * c_i + V^l * \dot{c}_i + V^s * [sum_j sum_m d_j \dot{q}_{j,i,m}]) - c_{in,i} * F_in + c_i * F_out == 0 - if (_nQsReacBulk > 0) + if (_nqsReactionBulk > 0) { int MoityIdx = 0; @@ -2284,14 +2280,14 @@ void CSTRModel::addTimeDerivativeJacobian(double t, double alpha, const ConstSim } // Volume: \dot{V} - F_{in} + F_{out} + F_{filter} == 0 mat.native(_nComp + _totalBound, _nComp + _totalBound) += alpha; - + /* std::cout << "Time derivative Jacobian: " << std::endl; for (unsigned int i = 0; i < 4; ++i) { for (unsigned int j = 0; j < 4; ++j) std::cout << mat.native(i, j) << " "; std::cout << std::endl; - } + }*/ } diff --git a/src/libcadet/model/StirredTankModel.hpp b/src/libcadet/model/StirredTankModel.hpp index 0b5d05834..51c605302 100644 --- a/src/libcadet/model/StirredTankModel.hpp +++ b/src/libcadet/model/StirredTankModel.hpp @@ -176,9 +176,9 @@ class CSTRModel : public UnitOperationBase double* _temp2; Eigen::MatrixXd _MconvMoityBulk; //!< Matrix with conservation of moieties in the bulk volume - std::vector _qsReacBulk; //!< Indices of components that are not conserved in the bulk volume - std::vector _QsCompBulk; - unsigned int _nQsReacBulk; + int const* _qsReactionBulk; //!< Indices of components that are not conserved in the bulk volume + std::vector _QsCompBulk; //!< Indices of components that are conserved in the bulk volume + unsigned int _nqsReactionBulk; int _nMoitiesBulk; class Exporter : public ISolutionExporter diff --git a/src/libcadet/model/reaction/CrystallizationReaction.cpp b/src/libcadet/model/reaction/CrystallizationReaction.cpp index d39e26954..bb7bb6c5d 100644 --- a/src/libcadet/model/reaction/CrystallizationReaction.cpp +++ b/src/libcadet/model/reaction/CrystallizationReaction.cpp @@ -197,7 +197,7 @@ namespace cadet virtual bool requiresConfiguration() const CADET_NOEXCEPT { return true; } virtual bool usesParamProviderInDiscretizationConfig() const CADET_NOEXCEPT { return false; } - void fillConservedMoietiesBulk(Eigen::MatrixXd& M, std::vector& QSReaction, std::vector& QSComponent) {} + void fillConservedMoietiesBulk(Eigen::MatrixXd& M, unsigned int& QSReaction, std::vector& QSComponent) {} template void jacobianQuasiSteadyLiquidImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, int reaction, const RowIterator& jac, LinearBufferAllocator workSpace) const { } @@ -206,13 +206,14 @@ namespace cadet void jacobianSingleFluxImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, int reaction, const RowIterator& jac, LinearBufferAllocator workSpace) const { } virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, active const* y, - Eigen::Map fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) { + Eigen::Map fluxes, int const* mapQSReac, LinearBufferAllocator workSpace) { return 0; } virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, - Eigen::Map fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) { + Eigen::Map fluxes, int const* mapQSReac, LinearBufferAllocator workSpace) { return 0; } + virtual int const* reactionQuasiStationarity() const CADET_NOEXCEPT { return nullptr; } virtual void timeDerivativeQuasiStationaryReaction(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double* dReacDt, LinearBufferAllocator workSpace){ } virtual bool configure(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx) { diff --git a/src/libcadet/model/reaction/DummyReaction.cpp b/src/libcadet/model/reaction/DummyReaction.cpp index 0e5ab5133..543edc5f1 100644 --- a/src/libcadet/model/reaction/DummyReaction.cpp +++ b/src/libcadet/model/reaction/DummyReaction.cpp @@ -87,7 +87,7 @@ class DummyDynamicReaction : public IDynamicReactionModel virtual void analyticJacobianLiquidSingleFluxAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, int reaction, linalg::BandedSparseRowIterator jac, LinearBufferAllocator workSpace) const { } virtual void timeDerivativeQuasiStationaryReaction(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double* dReacDt, LinearBufferAllocator workSpace) { } - + virtual int const* reactionQuasiStationarity() const CADET_NOEXCEPT { return nullptr; } #ifdef ENABLE_DG virtual void analyticJacobianLiquidAdd(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double factor, linalg::BandedEigenSparseRowIterator jac, LinearBufferAllocator workSpace) const { } #endif @@ -111,14 +111,14 @@ class DummyDynamicReaction : public IDynamicReactionModel virtual unsigned int numReactionsLiquid() const CADET_NOEXCEPT { return 0; } virtual unsigned int numReactionsCombined() const CADET_NOEXCEPT { return 0; } - void fillConservedMoietiesBulk(Eigen::MatrixXd& M, std::vector& QSReaction, std::vector& QSComponent) {} + void fillConservedMoietiesBulk(Eigen::MatrixXd& M, unsigned int & QSReaction, std::vector& QSComponent) {} bool hasQuasiStationaryReactionsLiquid() { return false;} virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, active const* y, - Eigen::Map fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) { return 0; } + Eigen::Map fluxes, int const* mapQSReac, LinearBufferAllocator workSpace) { return 0; } virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, - Eigen::Map fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) { return 0; } + Eigen::Map fluxes, int const* mapQSReac, LinearBufferAllocator workSpace) { return 0; } protected: }; diff --git a/src/libcadet/model/reaction/MassActionLawReaction.cpp b/src/libcadet/model/reaction/MassActionLawReaction.cpp index ed00fe13b..17384350f 100644 --- a/src/libcadet/model/reaction/MassActionLawReaction.cpp +++ b/src/libcadet/model/reaction/MassActionLawReaction.cpp @@ -341,6 +341,32 @@ class MassActionLawReactionBase : public DynamicReactionModelBase _expBulkFwd.resize(nComp, nReactions); _expBulkBwd.resize(nComp, nReactions); } + if (paramProvider.exists("IS_KINETIC")) + { + _reactionQuasistationarity.resize(_stoichiometryBulk.columns(), false); + + const std::vector vecKin = paramProvider.getIntArray("IS_KINETIC"); + if (vecKin.size() == 1) + { + // Treat an array with a single element as scalar + std::fill(_reactionQuasistationarity.begin(), _reactionQuasistationarity.end(), !static_cast(vecKin[0])); + } + else if (vecKin.size() < _reactionQuasistationarity.size()) + { + // Error on too few elements + throw InvalidParameterException("IS_KINETIC has to have at least " + std::to_string(_reactionQuasistationarity.size()) + " elements"); + } + else + { + // Copy what we need (ignore excess values) + std::transform(vecKin.begin(), vecKin.begin() + _reactionQuasistationarity.size(), _reactionQuasistationarity.begin(), [](int val) { return !static_cast(val); }); + } + } + else + { + const bool kineticBinding = paramProvider.getInt("IS_KINETIC"); + std::fill(_reactionQuasistationarity.begin(), _reactionQuasistationarity.end(), !kineticBinding); + } if (!nBound || !boundOffset) return true; @@ -377,6 +403,7 @@ class MassActionLawReactionBase : public DynamicReactionModelBase _expSolidBwdLiquid.resize(nComp, nReactions); } + return true; } @@ -389,17 +416,17 @@ class MassActionLawReactionBase : public DynamicReactionModelBase * @param _reactionQuasistationarity The reaction quasistationarity vector. * @return The conserved moieties matrix. */ - virtual void fillConservedMoietiesBulk(Eigen::MatrixXd& M, std::vector& mapQSReac, std::vector& _QsCompBulk) + virtual void fillConservedMoietiesBulk(Eigen::MatrixXd& M, unsigned int& numQSReac, std::vector& _QsCompBulk) { //1. get stoichmetic matrix with only reaction in quasi stationary // S dim -> ncomp x nreac - - int nQS = std::count(mapQSReac.begin(), mapQSReac.end(), 1); + _reactionQuasistationarity = {1}; + numQSReac = std::count(_reactionQuasistationarity.begin(), _reactionQuasistationarity.end(), 1); // Count the number of entries with value 1 in _reactionQuasistationarity - - - Eigen::MatrixXd QSS(_stoichiometryBulk.rows(), nQS); + + M.resize(numQSReac,_nComp); + Eigen::MatrixXd QSS(_stoichiometryBulk.rows(), numQSReac); // Fülle die Spalten basierend auf _reactionQuasistationarity @@ -418,7 +445,7 @@ class MassActionLawReactionBase : public DynamicReactionModelBase int rowIndex = 0; for (int j = 0; j < _stoichiometryBulk.columns(); j++) { - if (mapQSReac[j] == 0) + if (_reactionQuasistationarity[j] == 0) continue; QSS(i, rowIndex) = static_cast(_stoichiometryBulk.native(i, j)); rowIndex++; @@ -454,7 +481,7 @@ class MassActionLawReactionBase : public DynamicReactionModelBase //3. Test if the matrix is full rank int rang = QSS.fullPivLu().rank(); - if (rang != nQS) + if (rang != numQSReac) { throw std::runtime_error("Calculation Conserved Moities Matrix: The Stoichemetric Matrix is sigular"); } @@ -507,6 +534,7 @@ class MassActionLawReactionBase : public DynamicReactionModelBase linalg::ActiveDenseMatrix _expSolidBwdLiquid; std::vector _reactionQuasistationarity; + int _nQuasiStationaryReactionsLiquid; inline int maxNumReactions() const CADET_NOEXCEPT { return std::max(std::max(_stoichiometryBulk.columns(), _stoichiometryLiquid.columns()), _stoichiometrySolid.columns()); } @@ -696,10 +724,10 @@ class MassActionLawReactionBase : public DynamicReactionModelBase } virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, active const* y, - Eigen::Map fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) { return 0; } + Eigen::Map fluxes, int const* mapQSReac, LinearBufferAllocator workSpace) { return 0; } virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, - Eigen::Map fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) + Eigen::Map fluxes, int const* mapQSReac, LinearBufferAllocator workSpace) { typename ParamHandler_t::ParamsHandle const p = _paramHandler.update(t, secIdx, colPos, _nComp, _nBoundStates, workSpace); @@ -709,7 +737,7 @@ class MassActionLawReactionBase : public DynamicReactionModelBase double kFwdBulk_r = static_cast(p->kFwdBulk[r]); double kBwdBulk_r = static_cast(p->kBwdBulk[r]); - if (mapQSReac[r] == 1) + if (_reactionQuasistationarity[r] == 1) { double flow = singleFlux(r, y, kFwdBulk_r, kBwdBulk_r); fluxes[r] = flow; diff --git a/src/libcadet/model/reaction/MichaelisMentenReaction.cpp b/src/libcadet/model/reaction/MichaelisMentenReaction.cpp index f1517863b..d504d5ebc 100644 --- a/src/libcadet/model/reaction/MichaelisMentenReaction.cpp +++ b/src/libcadet/model/reaction/MichaelisMentenReaction.cpp @@ -149,10 +149,10 @@ class MichaelisMentenReactionBase : public DynamicReactionModelBase return true; } - + virtual int const* reactionQuasiStationarity() const CADET_NOEXCEPT { return nullptr; } virtual unsigned int numReactionsLiquid() const CADET_NOEXCEPT { return _stoichiometryBulk.columns(); } virtual unsigned int numReactionsCombined() const CADET_NOEXCEPT { return 0; } - void fillConservedMoietiesBulk(Eigen::MatrixXd& M, std::vector& QSReaction, std::vector& QSComponent) {} + void fillConservedMoietiesBulk(Eigen::MatrixXd& M, unsigned int& QSReaction, std::vector& QSComponent) {} template void jacobianQuasiSteadyLiquidImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, int reaction ,const RowIterator& jac, LinearBufferAllocator workSpace) const { } @@ -162,11 +162,11 @@ class MichaelisMentenReactionBase : public DynamicReactionModelBase virtual void timeDerivativeQuasiStationaryReaction(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double* dReacDt, LinearBufferAllocator workSpace) { } virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, active const* y, - Eigen::Map fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) { + Eigen::Map fluxes, int const* mapQSReac, LinearBufferAllocator workSpace) { return 0; } virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, - Eigen::Map fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) { + Eigen::Map fluxes, int const* mapQSReac, LinearBufferAllocator workSpace) { return 0; } diff --git a/src/libcadet/model/reaction/ReactionModelBase.hpp b/src/libcadet/model/reaction/ReactionModelBase.hpp index 420a1d254..27a9215f4 100644 --- a/src/libcadet/model/reaction/ReactionModelBase.hpp +++ b/src/libcadet/model/reaction/ReactionModelBase.hpp @@ -57,14 +57,14 @@ class DynamicReactionModelBase : public IDynamicReactionModel virtual active* getParameter(const ParameterId& pId); virtual void setExternalFunctions(IExternalFunction** extFuns, unsigned int size) { }; - virtual void fillConservedMoietiesBulk(Eigen::MatrixXd& M, std::vector& QSReaction, std::vector& _QsCompBulk) = 0; - + virtual void fillConservedMoietiesBulk(Eigen::MatrixXd& M, unsigned int& QSReaction, std::vector& _QsCompBulk) = 0; + virtual int const* reactionQuasiStationarity() const = 0; virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, active const* y, - Eigen::Map fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) = 0; + Eigen::Map fluxes, int const* mapQSReac, LinearBufferAllocator workSpace) = 0; virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, - Eigen::Map fluxes, std::vector mapQSReac, LinearBufferAllocator workSpace) = 0; + Eigen::Map fluxes, int const* mapQSReac, LinearBufferAllocator workSpace) = 0; virtual void timeDerivativeQuasiStationaryReaction(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double* dReacDt, LinearBufferAllocator workSpace) = 0; From c7a8d80ba7ed9dfbe42c4d0b58ae06ef06dfb789 Mon Sep 17 00:00:00 2001 From: "a.berger" Date: Wed, 5 Feb 2025 17:44:09 +0100 Subject: [PATCH 07/11] little changes --- src/libcadet/model/StirredTankModel.cpp | 37 +++++++++---------- .../model/reaction/MassActionLawReaction.cpp | 21 +++++------ 2 files changed, 27 insertions(+), 31 deletions(-) diff --git a/src/libcadet/model/StirredTankModel.cpp b/src/libcadet/model/StirredTankModel.cpp index a446d42bb..ace7f9c03 100644 --- a/src/libcadet/model/StirredTankModel.cpp +++ b/src/libcadet/model/StirredTankModel.cpp @@ -379,8 +379,6 @@ bool CSTRModel::configure(IParameterProvider& paramProvider) if (_qsReactionBulk != nullptr) { - // resize Conserved moities matrix - _dynReactionBulk->fillConservedMoietiesBulk(_MconvMoityBulk, _nqsReactionBulk, _QsCompBulk); // fill conserved moities matrix _nMoitiesBulk = _MconvMoityBulk.rows(); } @@ -449,21 +447,23 @@ unsigned int CSTRModel::threadLocalMemorySize() const CADET_NOEXCEPT if (_nqsReactionBulk > 0) { - // Memory for leanConsistentInitialTimeDerivative() - lms.add(_nComp + _totalBound); - lms.add(_nComp + _totalBound); - lms.add(_nComp); - lms.add(_nComp + _totalBound); - lms.add(numDofs()); - lms.add(_nonlinearSolver->workspaceSize(_nComp + _totalBound) * sizeof(double)); + // Memory for leanConsistentInitialTimeDerivative() + lms.add(_nComp); // map of kinetic components + lms.add(_nComp); // solution + lms.add(_nComp); // conserved quantities + //lms.add(_nComp + _totalBound); + lms.add(numDofs()); // fullX + lms.add(numDofs()); // fullRes + lms.add(_nonlinearSolver->workspaceSize(_nComp + _totalBound) * sizeof(double)); //non linear solver workspace lms.addBlock(resImplSize); lms.commit(); // Memory for leanConsistentInitialState() - lms.add(_totalBound); - lms.add(_nComp + 1); - lms.add(_nComp ); - lms.add(_totalBound); + //lms.add(_totalBound); + lms.add(_nComp); // dReacDt + lms.add(_nComp ); //map + lms.add(_nComp); //additional puffer + //lms.add(_totalBound); lms.commit(); } return lms.bufferSize(); @@ -1423,7 +1423,7 @@ void CSTRModel::leanConsistentInitialTimeDerivative(double t, double const* cons _jacFact.setAll(0.0); addTimeDerivativeJacobian(t, 1.0, ConstSimulationState{ vecStateY, nullptr }, _jacFact); - LinearBufferAllocator tlmAlloc = threadLocalMem.get(); // -> todo hier aufpassen ob ich den mem auch benutzen darf + //LinearBufferAllocator tlmAlloc = threadLocalMem.get(); // -> todo hier aufpassen ob ich den mem auch benutzen darf // Overwrite rows corresponding to algebraic equations with the Jacobian and set right hand side to 0 BufferedArray dReacDt = tlmAlloc.array(_nComp + 1); @@ -1642,7 +1642,7 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons _jac.native(state, _nComp + _totalBound) += _MconvMoityBulk(MoityIdx, i) * cDot[i]; // dF_{ci}/dVl = sum_i M[i,-] * dcidt } } - resCMoities[state] = dotProduct; + resCMoities[state] = dotProduct + v * flux[state]; MoityIdx++; state++; } @@ -1650,7 +1650,6 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons { resCMoities[state] = v * flux[state]; - if (wantJac) { @@ -1664,7 +1663,7 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons state = (_nComp - _nqsReactionBulk); for (unsigned int qsreac = 0; qsreac < _nqsReactionBulk; ++qsreac) { - resCMoities[state] = qsflux[qsreac]; + resCMoities[state] += qsflux[qsreac]; if(wantJac) { @@ -1797,14 +1796,14 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons // Volume: \dot{V} = F_{in} - F_{out} - F_{filter} res[2 * _nComp + _totalBound] = vDot - flowIn + flowOut + static_cast(_curFlowRateFilter); - /* + std::cout << "Jacobian: " << std::endl; for (unsigned int i = 0; i < 4; ++i) { for (unsigned int j = 0; j < 4; ++j) std::cout << _jac.native(i, j) << " "; std::cout << std::endl; - }*/ + } return 0; } diff --git a/src/libcadet/model/reaction/MassActionLawReaction.cpp b/src/libcadet/model/reaction/MassActionLawReaction.cpp index 17384350f..7a4a36e49 100644 --- a/src/libcadet/model/reaction/MassActionLawReaction.cpp +++ b/src/libcadet/model/reaction/MassActionLawReaction.cpp @@ -26,6 +26,7 @@ #include #include #include +#include using namespace Eigen; /* @@ -362,11 +363,6 @@ class MassActionLawReactionBase : public DynamicReactionModelBase std::transform(vecKin.begin(), vecKin.begin() + _reactionQuasistationarity.size(), _reactionQuasistationarity.begin(), [](int val) { return !static_cast(val); }); } } - else - { - const bool kineticBinding = paramProvider.getInt("IS_KINETIC"); - std::fill(_reactionQuasistationarity.begin(), _reactionQuasistationarity.end(), !kineticBinding); - } if (!nBound || !boundOffset) return true; @@ -421,13 +417,13 @@ class MassActionLawReactionBase : public DynamicReactionModelBase //1. get stoichmetic matrix with only reaction in quasi stationary // S dim -> ncomp x nreac - _reactionQuasistationarity = {1}; - numQSReac = std::count(_reactionQuasistationarity.begin(), _reactionQuasistationarity.end(), 1); + //_reactionQuasistationarity = {1}; + numQSReac = std::count(_reactionQuasistationarity.begin(), _reactionQuasistationarity.end(), true); // Count the number of entries with value 1 in _reactionQuasistationarity M.resize(numQSReac,_nComp); Eigen::MatrixXd QSS(_stoichiometryBulk.rows(), numQSReac); - + // Fülle die Spalten basierend auf _reactionQuasistationarity /*for (int i = 0; i < _stoichiometryBulk.rows(); ++i) @@ -440,16 +436,13 @@ class MassActionLawReactionBase : public DynamicReactionModelBase for (int i = 0; i < _stoichiometryBulk.rows(); ++i) { + int rowIndex = 0; for(int j = 0; j < _stoichiometryBulk.columns(); j++) { - int rowIndex = 0; - for (int j = 0; j < _stoichiometryBulk.columns(); j++) - { if (_reactionQuasistationarity[j] == 0) continue; QSS(i, rowIndex) = static_cast(_stoichiometryBulk.native(i, j)); rowIndex++; - } } } //Remove zero rows from QSS @@ -742,6 +735,10 @@ class MassActionLawReactionBase : public DynamicReactionModelBase double flow = singleFlux(r, y, kFwdBulk_r, kBwdBulk_r); fluxes[r] = flow; } + else + { + fluxes[r] = 0.0; + } return 0; } } From cc4c251d9b201710684a1c7d2a01d31a1e828af1 Mon Sep 17 00:00:00 2001 From: "a.berger" Date: Fri, 7 Feb 2025 12:40:28 +0100 Subject: [PATCH 08/11] restructre residum --- src/libcadet/model/StirredTankModel.cpp | 68 ++++++++----------- .../model/reaction/MassActionLawReaction.cpp | 4 ++ 2 files changed, 32 insertions(+), 40 deletions(-) diff --git a/src/libcadet/model/StirredTankModel.cpp b/src/libcadet/model/StirredTankModel.cpp index ace7f9c03..69406c016 100644 --- a/src/libcadet/model/StirredTankModel.cpp +++ b/src/libcadet/model/StirredTankModel.cpp @@ -1574,7 +1574,7 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons resC[i] += -flowIn * cIn[i] + flowOut * c[i]; } - if (wantJac && _nqsReactionBulk == 0) + if (wantJac) { _jac.setAll(0.0); @@ -1600,8 +1600,20 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons BufferedArray flux = subAlloc.array(_nComp); std::fill_n(static_cast(flux), _nComp, 0.0); + _dynReactionBulk->residualLiquidAdd(t, secIdx, colPos, c, static_cast(flux), -1.0, subAlloc); + + for (unsigned int comp = 0; comp < _nComp; ++comp) + resC[comp] += v * flux[comp]; + + if (wantJac) + { + for (unsigned int comp = 0; comp < _nComp; ++comp) + _jac.native(comp, _nComp + _totalBound) += static_cast(flux[comp]); // dF/dvliquid = flux + + _dynReactionBulk->analyticJacobianLiquidAdd(t, secIdx, colPos, reinterpret_cast(c), -static_cast(v), _jac.row(0), subAlloc); + } - if (_dynReactionBulk && (_nqsReactionBulk > 0)) + if (_nqsReactionBulk > 0) { _jac.setAll(0.0); @@ -1611,15 +1623,23 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons Eigen::Map> qsflux(_temp2, _nqsReactionBulk); qsflux.setZero(); - _dynReactionBulk->residualLiquidAdd(t, secIdx, colPos, c, static_cast(flux), -1.0, subAlloc); _dynReactionBulk->quasiStationaryFlux(t, secIdx, colPos, c, qsflux, _qsReactionBulk, subAlloc); - Eigen::Map> mapResC(resC, _nComp); int MoityIdx = 0; int state = 0; for (unsigned int comp = 0; comp < _nComp; comp++) { + if (_QsCompBulk[comp] == 0) + { + resCMoities[state] = resC[comp]; + if (wantJac) + { + _jac.native(state, comp) += (static_cast(vDot) + static_cast(flowOut)); // dF_{ci}/dcj = v_liquidDot + F_out + _dynReactionBulk->analyticJacobianLiquidSingleFluxAdd(t, secIdx, colPos, reinterpret_cast(c), state, comp, _jac.row(state), subAlloc); + } + state++; + } if (_QsCompBulk[comp] == 1 && MoityIdx < _nMoitiesBulk) { ResidualType dotProduct = 0.0; @@ -1627,37 +1647,18 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons { dotProduct += static_cast(_MconvMoityBulk(MoityIdx, i)) * (mapResC[i]); - if (wantJac && i == state) + if (wantJac) { - _jac.native(state, state) += _MconvMoityBulk(MoityIdx, i) * (static_cast(vDot) + static_cast(flowOut)); // dF_{ci}/dcj = v_liquidDot + F_out + _jac.native(state, i) += _MconvMoityBulk(MoityIdx, i) * (static_cast(vDot) + static_cast(flowOut)); // dF_{ci}/dcj = v_liquidDot + F_out if (cadet_likely(yDot)) - { _jac.native(i, _nComp + _totalBound) += _MconvMoityBulk(MoityIdx, i) * cDot[i]; // dF/dvliquid = cDot - } - } - if (wantJac && state != i) - { - _jac.native(state, i) += _MconvMoityBulk(MoityIdx, i)*(static_cast(vDot) + static_cast(flowOut)); // dF_{ci}/dcj = v_liquidDot + F_out - if (cadet_likely(yDot)) - _jac.native(state, _nComp + _totalBound) += _MconvMoityBulk(MoityIdx, i) * cDot[i]; // dF_{ci}/dVl = sum_i M[i,-] * dcidt + } } - resCMoities[state] = dotProduct + v * flux[state]; + resCMoities[state] = dotProduct; MoityIdx++; state++; } - else if (_QsCompBulk[comp] == 0) - { - resCMoities[state] = v * flux[state]; - - if (wantJac) - { - - _jac.native(state,comp) += (static_cast(vDot) + static_cast(flowOut)); // dF_{ci}/dcj = v_liquidDot + F_out - _dynReactionBulk->analyticJacobianLiquidSingleFluxAdd(t, secIdx, colPos, reinterpret_cast(c), state, comp, _jac.row(state), subAlloc); - } - state++; - } } state = (_nComp - _nqsReactionBulk); @@ -1674,19 +1675,6 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons } mapResC = resCMoities; } - else - { - for (unsigned int comp = 0; comp < _nComp; ++comp) - resC[comp] += v * flux[comp]; - - if (wantJac) - { - for (unsigned int comp = 0; comp < _nComp; ++comp) - _jac.native(comp, _nComp + _totalBound) += static_cast(flux[comp]); // dF/dvliquid = flux - - _dynReactionBulk->analyticJacobianLiquidAdd(t, secIdx, colPos, reinterpret_cast(c), -static_cast(v), _jac.row(0), subAlloc); - } - } // Bound states for (unsigned int type = 0; type < _nParType; ++type) diff --git a/src/libcadet/model/reaction/MassActionLawReaction.cpp b/src/libcadet/model/reaction/MassActionLawReaction.cpp index 7a4a36e49..a7bf3dee8 100644 --- a/src/libcadet/model/reaction/MassActionLawReaction.cpp +++ b/src/libcadet/model/reaction/MassActionLawReaction.cpp @@ -458,6 +458,8 @@ class MassActionLawReactionBase : public DynamicReactionModelBase nQScomp++; } } + if (nQScomp == 0) + return; // Redimensioniere QSS und kopiere nur die nicht-null Zeilen if (_nComp - nQScomp > 0 ) // if kinetic components exists we can resize QSS { @@ -982,6 +984,8 @@ class MassActionLawReactionBase : public DynamicReactionModelBase double* const fluxGradBwd = fluxGradFwd + _nComp; for (int r = 0; r < _stoichiometryBulk.columns(); ++r) { + if (_reactionQuasistationarity[r]) + continue; // Calculate gradients of forward and backward fluxes fluxGradLiquid(fluxGradFwd, r, _nComp, static_cast(p->kFwdBulk[r]), _expBulkFwd, y); fluxGradLiquid(fluxGradBwd, r, _nComp, static_cast(p->kBwdBulk[r]), _expBulkBwd, y); From debd5cb8d2f22dc4f76dd7453978e9e50b0ccc48 Mon Sep 17 00:00:00 2001 From: "a.berger" Date: Fri, 7 Feb 2025 17:48:53 +0100 Subject: [PATCH 09/11] fix jacobian bug --- src/libcadet/model/StirredTankModel.cpp | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/src/libcadet/model/StirredTankModel.cpp b/src/libcadet/model/StirredTankModel.cpp index 69406c016..7eb24f43f 100644 --- a/src/libcadet/model/StirredTankModel.cpp +++ b/src/libcadet/model/StirredTankModel.cpp @@ -1615,7 +1615,7 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons if (_nqsReactionBulk > 0) { - _jac.setAll(0.0); + //_jac.setAll(0.0); Eigen::Map> resCMoities(reinterpret_cast(_temp), _nComp); resCMoities.setZero(); @@ -1646,13 +1646,12 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons for (unsigned int i = 0 ; i < _MconvMoityBulk.cols(); ++i) // hier Optimierung durch Vermeidung von 0 Zeilen in MconvMoityBulk { dotProduct += static_cast(_MconvMoityBulk(MoityIdx, i)) * (mapResC[i]); - if (wantJac) { - _jac.native(state, i) += _MconvMoityBulk(MoityIdx, i) * (static_cast(vDot) + static_cast(flowOut)); // dF_{ci}/dcj = v_liquidDot + F_out + _jac.native(state, i) += (static_cast(vDot) + static_cast(flowOut)); // dF_{ci}/dcj = v_liquidDot + F_out + _jac.native(state, i) *= _MconvMoityBulk(MoityIdx, i) ; // dF_{ci}/dcj = v_liquidDot + F_out if (cadet_likely(yDot)) _jac.native(i, _nComp + _totalBound) += _MconvMoityBulk(MoityIdx, i) * cDot[i]; // dF/dvliquid = cDot - } } resCMoities[state] = dotProduct; @@ -1784,14 +1783,14 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons // Volume: \dot{V} = F_{in} - F_{out} - F_{filter} res[2 * _nComp + _totalBound] = vDot - flowIn + flowOut + static_cast(_curFlowRateFilter); - + /* std::cout << "Jacobian: " << std::endl; for (unsigned int i = 0; i < 4; ++i) { for (unsigned int j = 0; j < 4; ++j) std::cout << _jac.native(i, j) << " "; std::cout << std::endl; - } + }*/ return 0; } From 8b64fc1bb21e32ccc580eb4e7c205a9c7a88e9d6 Mon Sep 17 00:00:00 2001 From: "a.berger" Date: Tue, 11 Feb 2025 14:19:20 +0100 Subject: [PATCH 10/11] fix Consistent initalisation --- src/libcadet/model/StirredTankModel.cpp | 35 ++++++++++++++----------- 1 file changed, 19 insertions(+), 16 deletions(-) diff --git a/src/libcadet/model/StirredTankModel.cpp b/src/libcadet/model/StirredTankModel.cpp index 7eb24f43f..267028b40 100644 --- a/src/libcadet/model/StirredTankModel.cpp +++ b/src/libcadet/model/StirredTankModel.cpp @@ -1157,7 +1157,7 @@ void CSTRModel::leanConsistentInitialState(const SimulationTime& simTime, double } if (_nqsReactionBulk > 0) { - //Todo check meomory -> something is wring with fullX + LinearBufferAllocator tlmAlloc = threadLocalMem.get(); BufferedArray qsMask = tlmAlloc.array(_nComp); @@ -1180,7 +1180,7 @@ void CSTRModel::leanConsistentInitialState(const SimulationTime& simTime, double { double dotProduct = 0.0; for (unsigned int i = 0; i < _MconvMoityBulk.cols(); ++i) - dotProduct += _MconvMoityBulk(MoityIdx, i) * (c[i]); + dotProduct += _MconvMoityBulk(MoityIdx, i) * c[i]; conservedQuants[MoityIdx] = dotProduct; } @@ -1258,15 +1258,17 @@ void CSTRModel::leanConsistentInitialState(const SimulationTime& simTime, double mat.submatrixSetAll(0.0, 0, 0, _nMoitiesBulk, probSize); unsigned int rIdx = 0; for (unsigned int MoityIdx = 0; MoityIdx < _nMoitiesBulk; ++MoityIdx) - { + { + double dotProduct = 0.0; - for (unsigned int i = 0; i < _MconvMoityBulk.cols(); ++i) - {// hier Optimierung durch Vermeidung von 0 Zeilen in MconvMoityBulk + int j = 0; + for (unsigned int i = 0; i < _MconvMoityBulk.cols(); ++i) { if (!mask.mask[i]) continue; - mat.native(rIdx, i) = _MconvMoityBulk(MoityIdx, i); - rIdx++; + mat.native(rIdx, j) += _MconvMoityBulk(MoityIdx, i); + j++; } + rIdx++; } return true; }; @@ -1279,7 +1281,7 @@ void CSTRModel::leanConsistentInitialState(const SimulationTime& simTime, double { // Prepare input vector by overwriting masked items std::copy_n(c, mask.len, fullX + _nComp); - linalg::applyVectorSubset(x, mask, fullX + _nComp); + linalg::applyVectorSubset(x, mask, fullX + _nComp); //Todo hier ist ein Fehler ! mask genauer anschauen! // Call residual function residualImpl(simTime.t, simTime.secIdx, fullX, nullptr, fullResidual, tlmAlloc.manageRemainingMemory()); @@ -1292,16 +1294,16 @@ void CSTRModel::leanConsistentInitialState(const SimulationTime& simTime, double for (unsigned int MoityIdx = 0; MoityIdx < _nMoitiesBulk; ++MoityIdx) { double dotProduct = 0.0; - for (unsigned int j = 0; j < _MconvMoityBulk.cols(); ++j) // hier Optimierung durch Vermeidung von 0 Zeilen in MconvMoityBulk + for (unsigned int j = 0; j < _MconvMoityBulk.cols(); ++j) { if(!mask.mask[j]) continue; dotProduct += _MconvMoityBulk(MoityIdx, j) * x[j]; - r[rIdx] = dotProduct - conservedQuants[rIdx]; - rIdx++; } + r[rIdx] = dotProduct - conservedQuants[rIdx]; + rIdx++; } - + std::cout << "Residual: " << std::endl; for (unsigned int i = 0; i < probSize; ++i) { std::cout << r[i] << std::endl; @@ -1316,10 +1318,11 @@ void CSTRModel::leanConsistentInitialState(const SimulationTime& simTime, double // Apply solution linalg::applyVectorSubset(static_cast(solution), mask, c); + std::cout << "Solution: " << std::endl; for (unsigned int i = 0; i < _nComp; ++i) std::cout << solution[i] << std::endl; - + // Refine / correct solution } @@ -1416,7 +1419,7 @@ void CSTRModel::leanConsistentInitialTimeDerivative(double t, double const* cons return; // todo not implmenteted yet for (unsigned int i = 0; i < numDofs(); ++i) - vecStateYdot[i] = -vecStateYdot[i] - res[i]; + vecStateYdot[i] = -vecStateYdot[i]- res[i]; // Assemble time derivative Jacobian @@ -1461,7 +1464,7 @@ void CSTRModel::leanConsistentInitialTimeDerivative(double t, double const* cons { LOG(Error) << "Solve() failed"; } - + /* for (int i = 0; i < _nComp + 1; i++) std::cout << " resC" << i << ": " << resC[i] << std::endl; @@ -1471,7 +1474,7 @@ void CSTRModel::leanConsistentInitialTimeDerivative(double t, double const* cons for (int i = 0; i < _nComp + 1; i++) std::cout << " Y" << i <<": " << c[i] << std::endl; - + */ } else { From be69201191b2f6f54b19ab62c159965517014f60 Mon Sep 17 00:00:00 2001 From: "a.berger" Date: Thu, 13 Feb 2025 17:09:11 +0100 Subject: [PATCH 11/11] add characterisation map for dynamic, conserv and algebraic state entries --- src/libcadet/model/StirredTankModel.cpp | 205 +++++++++++------- src/libcadet/model/StirredTankModel.hpp | 3 +- .../model/reaction/MassActionLawReaction.cpp | 59 +++-- 3 files changed, 153 insertions(+), 114 deletions(-) diff --git a/src/libcadet/model/StirredTankModel.cpp b/src/libcadet/model/StirredTankModel.cpp index 267028b40..9fcfdd1e8 100644 --- a/src/libcadet/model/StirredTankModel.cpp +++ b/src/libcadet/model/StirredTankModel.cpp @@ -30,7 +30,7 @@ #include #include - +#include #include #include @@ -243,7 +243,6 @@ bool CSTRModel::configureModelDiscretization(IParameterProvider& paramProvider, paramProvider.popScope(); - _temp = new active[_nComp]; _temp2 = new double[_nComp]; @@ -377,10 +376,32 @@ bool CSTRModel::configure(IParameterProvider& paramProvider) // Get quasi stationary reactions information vector _qsReactionBulk = _dynReactionBulk->reactionQuasiStationarity(); - if (_qsReactionBulk != nullptr) + if (false) { _dynReactionBulk->fillConservedMoietiesBulk(_MconvMoityBulk, _nqsReactionBulk, _QsCompBulk); // fill conserved moities matrix _nMoitiesBulk = _MconvMoityBulk.rows(); + + int mIdx = 0; + for (int state = 0; state < _nComp; state++) + { + std::bitset<3> c; + if (_QsCompBulk[state] == 0) // state is dynamic + { + c.set(0); + stateMap.push_back(c); + } + else if(mIdx < _nMoitiesBulk) // state is dynamic and calculated with conserved moities + { + c.set(1); + stateMap.push_back(c); + mIdx++; + } + else if(mIdx >= _nMoitiesBulk) // state algebraic + { + c.set(2); + stateMap.push_back(c); + } + } } } @@ -1161,8 +1182,16 @@ void CSTRModel::leanConsistentInitialState(const SimulationTime& simTime, double LinearBufferAllocator tlmAlloc = threadLocalMem.get(); BufferedArray qsMask = tlmAlloc.array(_nComp); + for(int state = 0 ; state < _nComp; state++) + { + if(stateMap[state].test(1) || stateMap[state].test(2)) + qsMask[state] = 1; + else + qsMask[state] = 0; + } + // Mark components with quasi-stationary reactions partition - qsMask.copyFromVector(_QsCompBulk); + //qsMask.copyFromVector(_QsCompBulk); const linalg::ConstMaskArray mask{ static_cast(qsMask), static_cast(_nComp) }; const int probSize = linalg::numMaskActive(mask); @@ -1179,8 +1208,9 @@ void CSTRModel::leanConsistentInitialState(const SimulationTime& simTime, double for (unsigned int MoityIdx = 0; MoityIdx < _nMoitiesBulk; ++MoityIdx) { double dotProduct = 0.0; - for (unsigned int i = 0; i < _MconvMoityBulk.cols(); ++i) + for (unsigned int i = 0; i < _MconvMoityBulk.cols(); ++i) dotProduct += _MconvMoityBulk(MoityIdx, i) * c[i]; + conservedQuants[MoityIdx] = dotProduct; } @@ -1255,25 +1285,27 @@ void CSTRModel::leanConsistentInitialState(const SimulationTime& simTime, double mat.setAll(0.0); linalg::copyMatrixSubset(_jac, mask, mask, mat); // Replace upper part with conservation relations - mat.submatrixSetAll(0.0, 0, 0, _nMoitiesBulk, probSize); + //mat.submatrixSetAll(0.0, 0, 0, _nMoitiesBulk, probSize); unsigned int rIdx = 0; - for (unsigned int MoityIdx = 0; MoityIdx < _nMoitiesBulk; ++MoityIdx) + unsigned int MoityIdx = 0; + for (unsigned int state = 0; state < _nComp; ++state) { - - double dotProduct = 0.0; - int j = 0; - for (unsigned int i = 0; i < _MconvMoityBulk.cols(); ++i) { - if (!mask.mask[i]) - continue; - mat.native(rIdx, j) += _MconvMoityBulk(MoityIdx, i); - j++; + if (stateMap[state].test(1)) + { + double dotProduct = 0.0; + int j = 0; + for (unsigned int i = 0; i < _MconvMoityBulk.cols(); ++i) { + if (_QsCompBulk[i]==0) + continue; + mat.native(rIdx, j) = _MconvMoityBulk(MoityIdx, i); + j++; + } + rIdx++; } - rIdx++; } return true; }; } - std::copy_n(vecStateY, _nComp, fullX); fullX[(2 * _nComp + _totalBound)] = c[(_nComp + _totalBound)]; _nonlinearSolver->solve( @@ -1281,27 +1313,28 @@ void CSTRModel::leanConsistentInitialState(const SimulationTime& simTime, double { // Prepare input vector by overwriting masked items std::copy_n(c, mask.len, fullX + _nComp); - linalg::applyVectorSubset(x, mask, fullX + _nComp); //Todo hier ist ein Fehler ! mask genauer anschauen! + linalg::applyVectorSubset(x, mask, fullX + _nComp); // Call residual function residualImpl(simTime.t, simTime.secIdx, fullX, nullptr, fullResidual, tlmAlloc.manageRemainingMemory()); // Extract values from residual linalg::selectVectorSubset(fullResidual + _nComp, mask, r); - std::fill_n(r, _nMoitiesBulk, 0.0); + //std::fill_n(r, _nMoitiesBulk, 0.0); - int rIdx = 0; - for (unsigned int MoityIdx = 0; MoityIdx < _nMoitiesBulk; ++MoityIdx) + int MoityIdx = 0; + int rIdx = 0; + for (unsigned int state = 0; state < _nComp; ++state) { - double dotProduct = 0.0; - for (unsigned int j = 0; j < _MconvMoityBulk.cols(); ++j) - { - if(!mask.mask[j]) - continue; - dotProduct += _MconvMoityBulk(MoityIdx, j) * x[j]; + if (stateMap[state].test(1)) + { + double dotProduct = 0.0; + for (unsigned int j = 0; j < _MconvMoityBulk.cols(); ++j) + dotProduct += _MconvMoityBulk(MoityIdx, j) * x[j]; + r[rIdx] = dotProduct - conservedQuants[MoityIdx]; + MoityIdx++; + rIdx++; } - r[rIdx] = dotProduct - conservedQuants[rIdx]; - rIdx++; } std::cout << "Residual: " << std::endl; @@ -1309,9 +1342,8 @@ void CSTRModel::leanConsistentInitialState(const SimulationTime& simTime, double std::cout << r[i] << std::endl; } std::cout << "Solution: " << std::endl; - for (unsigned int i = 0; i < _nComp; ++i) + for (unsigned int i = 0; i < probSize; ++i) std::cout << x[i] << std::endl; - return true; }, jacFunc, errorTol, static_cast(solution), nonlinMem, jacobianMatrix, probSize); @@ -1419,7 +1451,7 @@ void CSTRModel::leanConsistentInitialTimeDerivative(double t, double const* cons return; // todo not implmenteted yet for (unsigned int i = 0; i < numDofs(); ++i) - vecStateYdot[i] = -vecStateYdot[i]- res[i]; + vecStateYdot[i] = -vecStateYdot[i]; // Assemble time derivative Jacobian @@ -1464,17 +1496,7 @@ void CSTRModel::leanConsistentInitialTimeDerivative(double t, double const* cons { LOG(Error) << "Solve() failed"; } - /* - for (int i = 0; i < _nComp + 1; i++) - std::cout << " resC" << i << ": " << resC[i] << std::endl; - for (int i = 0; i < _nComp + 1; i++) - std::cout <<" Ydot" << i << ": " << cDot[i] << std::endl; - - - for (int i = 0; i < _nComp + 1; i++) - std::cout << " Y" << i <<": " << c[i] << std::endl; - */ } else { @@ -1535,7 +1557,21 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons const ParamType flowIn = static_cast(_flowRateIn); const ParamType flowOut = static_cast(_flowRateOut); - bool wantJac = true; // just for debugg reasion + //bool wantJac = true; // just for debugg reasion + + /* testcase + std::bitset<3> c1; + std::bitset<3> c2; + std::bitset<3> c3; + + c1.set(0); // state 0 is dynamic + c2.set(1); // state 1 is conserved + c3.set(2); // state 2 is algebraic + + stateMap.push_back(c1); + stateMap.push_back(c2); + stateMap.push_back(c3); */ + // Inlet DOF for (unsigned int i = 0; i < _nComp; ++i) { @@ -1630,23 +1666,30 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons Eigen::Map> mapResC(resC, _nComp); int MoityIdx = 0; - int state = 0; - for (unsigned int comp = 0; comp < _nComp; comp++) + int comp = 0; + int qsreac = 0; + for (unsigned int state = 0; state < _nComp; state++) { - if (_QsCompBulk[comp] == 0) - { - resCMoities[state] = resC[comp]; - if (wantJac) + if (stateMap[state].test(0)) // dynamic + { + for (comp = 0; comp < _nComp; comp++) { - _jac.native(state, comp) += (static_cast(vDot) + static_cast(flowOut)); // dF_{ci}/dcj = v_liquidDot + F_out - _dynReactionBulk->analyticJacobianLiquidSingleFluxAdd(t, secIdx, colPos, reinterpret_cast(c), state, comp, _jac.row(state), subAlloc); + if(_QsCompBulk[comp] == 1) + continue; + resCMoities[state] = resC[comp]; + if (wantJac) + { + _jac.native(state, state) = 0.0; + _jac.native(state, comp) = (static_cast(vDot) + static_cast(flowOut)); // dF_{ci}/dcj = v_liquidDot + F_out + _dynReactionBulk->analyticJacobianLiquidSingleFluxAdd(t, secIdx, colPos, reinterpret_cast(c), state, comp, _jac.row(state), subAlloc); + } } - state++; } - if (_QsCompBulk[comp] == 1 && MoityIdx < _nMoitiesBulk) - { + else if (stateMap[state].test(1)) // conserved + { + // todo test of matrix times vector faster ResidualType dotProduct = 0.0; - for (unsigned int i = 0 ; i < _MconvMoityBulk.cols(); ++i) // hier Optimierung durch Vermeidung von 0 Zeilen in MconvMoityBulk + for (unsigned int i = 0 ; i < _MconvMoityBulk.cols(); ++i) { dotProduct += static_cast(_MconvMoityBulk(MoityIdx, i)) * (mapResC[i]); if (wantJac) @@ -1657,24 +1700,22 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons _jac.native(i, _nComp + _totalBound) += _MconvMoityBulk(MoityIdx, i) * cDot[i]; // dF/dvliquid = cDot } } + resCMoities[state] = dotProduct; MoityIdx++; - state++; } - } + else if (stateMap[state].test(2)) // algebraic + { + resCMoities[state] = qsflux[qsreac]; - state = (_nComp - _nqsReactionBulk); - for (unsigned int qsreac = 0; qsreac < _nqsReactionBulk; ++qsreac) - { - resCMoities[state] += qsflux[qsreac]; - - if(wantJac) - { - _dynReactionBulk->analyticQuasiSteadyJacobianLiquid(t, secIdx, colPos, reinterpret_cast(c), state , qsreac, _jac.row(state), subAlloc); - _jac.native(state, _nComp + _totalBound) = 0.0; // dF_{ci}/dvliquid = 0 + if (wantJac) + { + _dynReactionBulk->analyticQuasiSteadyJacobianLiquid(t, secIdx, colPos, reinterpret_cast(c), state, qsreac, _jac.row(state), subAlloc); + _jac.native(state, _nComp + _totalBound) = 0.0; // dF_{ci}/dvliquid = 0 + } + qsreac++; } - state++; - } + } mapResC = resCMoities; } @@ -2192,27 +2233,29 @@ void CSTRModel::addTimeDerivativeJacobian(double t, double alpha, const ConstSim int MoityIdx = 0; int state = 0; - for (unsigned int i = 0; i < _nComp; i++) + for (unsigned int state = 0; state < _nComp; state++) { - if (_QsCompBulk[i] == 1 && MoityIdx < _nMoitiesBulk) + if (stateMap[state].test(0)) { + for (int comp = 0; comp < _nComp; comp++) + { + if (_QsCompBulk[comp] == 1) + continue; + mat.native(state, comp) += timeV; // dRes / dcDot - mat.native(i, i) += timeV * (static_cast(_MconvMoityBulk(MoityIdx, MoityIdx))); // dRes / dcDot + } + } + if (stateMap[state].test(1)) + { + mat.native(state, state) += timeV * (static_cast(_MconvMoityBulk(MoityIdx, MoityIdx))); // dRes / dcDot for (int comp = 0; comp < _nComp; comp++) { - if (_QsCompBulk[comp] == 1 && i != comp) - mat.native(i, comp) += timeV * static_cast(_MconvMoityBulk(MoityIdx, comp)); // dRes / dcDot - // todo this is wrong - mat.native(i, _nComp + _totalBound) += alpha * static_cast(_MconvMoityBulk(MoityIdx, comp)) * c[comp]; // dRes / dVlDot + if (_QsCompBulk[comp] == 1 && state != comp) + mat.native(state, comp) += timeV * static_cast(_MconvMoityBulk(MoityIdx, comp)); // dRes / dcDot + mat.native(state, _nComp + _totalBound) += alpha * static_cast(_MconvMoityBulk(MoityIdx, comp)) * c[comp]; // dRes / dVlDot } MoityIdx++; - state++; - } - else if (_QsCompBulk[i] == 0) - { - mat.native(state, i) += timeV; // dRes / dcDot - state++; } } } diff --git a/src/libcadet/model/StirredTankModel.hpp b/src/libcadet/model/StirredTankModel.hpp index 51c605302..2b1d8368e 100644 --- a/src/libcadet/model/StirredTankModel.hpp +++ b/src/libcadet/model/StirredTankModel.hpp @@ -24,7 +24,7 @@ #include "linalg/DenseMatrix.hpp" #include "model/ModelUtils.hpp" #include "Memory.hpp" - +#include #include #include @@ -180,6 +180,7 @@ class CSTRModel : public UnitOperationBase std::vector _QsCompBulk; //!< Indices of components that are conserved in the bulk volume unsigned int _nqsReactionBulk; int _nMoitiesBulk; + std::vector> stateMap; // !< Bitset that tracks the properties of the state vector (dynamic, moities, algebraic) class Exporter : public ISolutionExporter { diff --git a/src/libcadet/model/reaction/MassActionLawReaction.cpp b/src/libcadet/model/reaction/MassActionLawReaction.cpp index a7bf3dee8..264a5ec0b 100644 --- a/src/libcadet/model/reaction/MassActionLawReaction.cpp +++ b/src/libcadet/model/reaction/MassActionLawReaction.cpp @@ -347,8 +347,13 @@ class MassActionLawReactionBase : public DynamicReactionModelBase _reactionQuasistationarity.resize(_stoichiometryBulk.columns(), false); const std::vector vecKin = paramProvider.getIntArray("IS_KINETIC"); - if (vecKin.size() == 1) + int numqsReaction = std::count(vecKin.begin(), vecKin.end(), 1); + if (numqsReaction == 0) { + _reactionQuasistationarity = {}; + } + else if (vecKin.size() == 1) + { // Treat an array with a single element as scalar std::fill(_reactionQuasistationarity.begin(), _reactionQuasistationarity.end(), !static_cast(vecKin[0])); } @@ -357,10 +362,11 @@ class MassActionLawReactionBase : public DynamicReactionModelBase // Error on too few elements throw InvalidParameterException("IS_KINETIC has to have at least " + std::to_string(_reactionQuasistationarity.size()) + " elements"); } - else + else if (numqsReaction != 0) { // Copy what we need (ignore excess values) std::transform(vecKin.begin(), vecKin.begin() + _reactionQuasistationarity.size(), _reactionQuasistationarity.begin(), [](int val) { return !static_cast(val); }); + } } @@ -417,23 +423,19 @@ class MassActionLawReactionBase : public DynamicReactionModelBase //1. get stoichmetic matrix with only reaction in quasi stationary // S dim -> ncomp x nreac - //_reactionQuasistationarity = {1}; - numQSReac = std::count(_reactionQuasistationarity.begin(), _reactionQuasistationarity.end(), true); + // Count the number of entries with value 1 in _reactionQuasistationarity + numQSReac = std::count(_reactionQuasistationarity.begin(), _reactionQuasistationarity.end(), true); - M.resize(numQSReac,_nComp); - Eigen::MatrixXd QSS(_stoichiometryBulk.rows(), numQSReac); - + if (numQSReac == 0) + return; - // Fülle die Spalten basierend auf _reactionQuasistationarity - /*for (int i = 0; i < _stoichiometryBulk.rows(); ++i) - { - const double* rowData = reinterpret_cast(_stoichiometryBulk.data() + i * _stoichiometryBulk.columns()); - Eigen::Map Srow(rowData, _stoichiometryBulk.columns()); - QSS.row(rowIndex++) = Srow; - } - */ + M.resize(numQSReac,_nComp); + // QSS quasistationary stoichiometry matrix contains only coloums with qs reactions + Eigen::MatrixXd QSS(_stoichiometryBulk.rows(), numQSReac); + + // Fill QSS with the stoichiometry of the quasi stationary reactions for (int i = 0; i < _stoichiometryBulk.rows(); ++i) { int rowIndex = 0; @@ -445,11 +447,11 @@ class MassActionLawReactionBase : public DynamicReactionModelBase rowIndex++; } } - //Remove zero rows from QSS + int nQScomp = 0; // Number of quasi stationary active components for (int i = 0; i < QSS.rows(); ++i) { - if (QSS.row(i).norm() < 1e-10) { + if (QSS.row(i).norm() < 1e-15) { _QsCompBulk.push_back(0); } else @@ -458,30 +460,27 @@ class MassActionLawReactionBase : public DynamicReactionModelBase nQScomp++; } } - if (nQScomp == 0) - return; - // Redimensioniere QSS und kopiere nur die nicht-null Zeilen + + // optional scip this step if (_nComp - nQScomp > 0 ) // if kinetic components exists we can resize QSS { Eigen::MatrixXd QSSCompressed(nQScomp, QSS.cols()); for (std::size_t i = 0; i < _QsCompBulk.size(); ++i) { if (_QsCompBulk[i] == 1) - { QSSCompressed.row(i) = QSS.row(i); - } } QSS.swap(QSSCompressed); // Swap the compressed matrix back into QSS } - //3. Test if the matrix is full rank - int rang = QSS.fullPivLu().rank(); - if (rang != numQSReac) + //Test if the matrix is full rank + int rank = QSS.fullPivLu().rank(); + if (rank != numQSReac) { throw std::runtime_error("Calculation Conserved Moities Matrix: The Stoichemetric Matrix is sigular"); } - //3. Calculate the null space of the matrix + //Calculate the null space of the matrix Eigen::MatrixXd leftZeroSpace = QSS.transpose().fullPivLu().kernel().transpose(); if (_nComp - nQScomp < 0) // if there are no kinetic components we can return the zero matrix @@ -501,7 +500,6 @@ class MassActionLawReactionBase : public DynamicReactionModelBase } M = leftZeroSpace; } - } virtual unsigned int numReactionsLiquid() const CADET_NOEXCEPT { return _stoichiometryBulk.columns(); } @@ -737,12 +735,9 @@ class MassActionLawReactionBase : public DynamicReactionModelBase double flow = singleFlux(r, y, kFwdBulk_r, kBwdBulk_r); fluxes[r] = flow; } - else - { - fluxes[r] = 0.0; - } - return 0; + } + return 0; } template