Skip to content

Commit

Permalink
Restart implementation
Browse files Browse the repository at this point in the history
  • Loading branch information
@AntoniaBerger
AntoniaBerger committed Jan 15, 2025
1 parent 5333e26 commit 6a4367e
Show file tree
Hide file tree
Showing 9 changed files with 303 additions and 18 deletions.
19 changes: 17 additions & 2 deletions src/libcadet/model/ReactionModel.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<int>& 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
Expand All @@ -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<int>& QSReaction, std::vector<int>& _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
Expand Down Expand Up @@ -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:
};
Expand Down
110 changes: 98 additions & 12 deletions src/libcadet/model/StirredTankModel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,10 @@
#include <algorithm>
#include <functional>

#include <Eigen/Dense>

#include<iostream>

namespace cadet
{

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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<IDynamicReactionModel*>(_nParType, nullptr);
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
{
Expand Down Expand Up @@ -1247,7 +1276,7 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons

const ParamType flowIn = static_cast<ParamType>(_flowRateIn);
const ParamType flowOut = static_cast<ParamType>(_flowRateOut);

// Inlet DOF
for (unsigned int i = 0; i < _nComp; ++i)
{
Expand All @@ -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<ParamType>(_constSolidVolume);
ResidualType* const resC = res + _nComp;
ResidualType* resC = res + _nComp;
for (unsigned int i = 0; i < _nComp; ++i)
{
resC[i] = 0.0;
Expand Down Expand Up @@ -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<ResidualType> flux = subAlloc.array<ResidualType>(_nComp);
LinearBufferAllocator subAlloc = tlmAlloc.manageRemainingMemory(); // todo nachgucken ob das wirklich so geht

BufferedArray<ResidualType> flux = subAlloc.array<ResidualType>(_nComp);
std::fill_n(static_cast<ResidualType*>(flux), _nComp, 0.0);

BufferedArray<ResidualType> qsflux = subAlloc.array<ResidualType>(_nQsReacBulk);
std::fill_n(static_cast<ResidualType*>(qsflux), _nQsReacBulk, 0.0);

BufferedArray<ResidualType> temp = subAlloc.array<ResidualType>(_nComp);
std::fill_n(static_cast<ResidualType*>(temp), _nComp, 0.0);

_dynReactionBulk->residualLiquidAdd(t, secIdx, colPos, c, static_cast<ResidualType*>(flux), -1.0, subAlloc);
_dynReactionBulk->quasiStationaryFlux(t, secIdx, colPos, c, static_cast<ResidualType*>(qsflux), _qsReacBulk, subAlloc);



Eigen::Map<Eigen::Vector<ResidualType, Eigen::Dynamic>> resCMoities(static_cast<ResidualType*>(temp), _nComp);
resCMoities.setZero();

Eigen::Map<Eigen::Vector<ResidualType, Eigen::Dynamic>> mapResC(resC, _nComp);
std::vector<int> 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<ResidualType>(_MconvMoityBulk(MoityIdx, i)) * (mapResC[i]);

if (wantJac)
_jac.native(i, state) = vDotTimeFactor + static_cast<double>(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<double>(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];

Expand All @@ -1331,7 +1419,7 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons
_dynReactionBulk->analyticJacobianLiquidAdd(t, secIdx, colPos, reinterpret_cast<double const*>(c), -static_cast<double>(v), _jac.row(0), subAlloc);
}
}

// Bound states
for (unsigned int type = 0; type < _nParType; ++type)
{
Expand Down Expand Up @@ -1401,7 +1489,6 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons
}
}
}

if (wantJac)
{
// Assemble Jacobian: Reaction
Expand Down Expand Up @@ -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<ParamType>(_curFlowRateFilter);

return 0;
}

Expand Down
7 changes: 7 additions & 0 deletions src/libcadet/model/StirredTankModel.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<int> _qsReacBulk; //!< Indices of components that are not conserved in the bulk volume
std::vector<int> _QsCompBulk;
unsigned int _nQsReacBulk;
unsigned int _nMoitiesBulk;

class Exporter : public ISolutionExporter
{
public:
Expand Down
5 changes: 4 additions & 1 deletion src/libcadet/model/reaction/CrystallizationReaction.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<int>& QSReaction, std::vector<int>& QSComponent) { }
bool hasQuasiStationaryReactionsLiquid() { return false; }
virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y,
active* fluxes, std::vector<int>& mapQSReac, LinearBufferAllocator workSpace) const { return 0; }
virtual bool configure(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx)
{
readScalarParameterOrArray(_bins, paramProvider, "CRY_BINS", 1);
Expand Down
6 changes: 6 additions & 0 deletions src/libcadet/model/reaction/DummyReaction.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<int>& QSReaction, std::vector<int>& QSComponent) {}
bool hasQuasiStationaryReactionsLiquid() { return false;}

virtual int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y,
active* fluxes, std::vector<int>& mapQSReac, LinearBufferAllocator workSpace) const {return 0;}

protected:
};
Expand Down
Loading

0 comments on commit 6a4367e

Please sign in to comment.