Skip to content

Commit 6ddcd08

Browse files
AntoniaBergerAntoniaBerger
committed
fix2
1 parent 41e7ff3 commit 6ddcd08

File tree

3 files changed

+23
-124
lines changed

3 files changed

+23
-124
lines changed

‎src/libcadet/model/ReactionModel.hpp

Lines changed: 2 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -264,9 +264,8 @@ class IDynamicReactionModel
264264
* @param [in,out] M Matrix to be filled with the conserved moieties
265265
* @param [in,out] QsCompBulk Vector to be filled with the components involved in qs reactions
266266
*/
267-
virtual void fillConservedMoietiesBulk(Eigen::MatrixXd& M, std::vector<int>& QsCompBulk) {}
268-
virtual void fillConservedMoietiesBulk2(Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>& M, int& nConservedQuants, std::vector<int>& QsCompBulk) {}
269-
virtual void fillConservedMoietiesBulk2(Eigen::Matrix<active, Eigen::Dynamic, Eigen::Dynamic>& M, int& nConservedQuants, std::vector<int>& QsCompBulk) {}
267+
virtual void fillConservedMoietiesBulk(Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>& M, int& nConservedQuants, std::vector<int>& QsCompBulk) {}
268+
virtual void fillConservedMoietiesBulk(Eigen::Matrix<active, Eigen::Dynamic, Eigen::Dynamic>& M, int& nConservedQuants, std::vector<int>& QsCompBulk) {}
270269

271270

272271
/**

‎src/libcadet/model/StirredTankModel.cpp

Lines changed: 20 additions & 115 deletions
Original file line numberDiff line numberDiff line change
@@ -34,7 +34,7 @@
3434
#include <Eigen/Dense>
3535
#include <Eigen/Sparse>
3636

37-
#include<iostream>
37+
#include<iostream>//todo delete
3838

3939
namespace cadet
4040
{
@@ -82,8 +82,7 @@ CSTRModel::~CSTRModel() CADET_NOEXCEPT
8282
delete[] _offsetParType;
8383

8484
delete[] _dynReactionBulk;
85-
//delete[] _temp;
86-
//delete[] _temp2;
85+
8786
}
8887

8988
unsigned int CSTRModel::numDofs() const CADET_NOEXCEPT
@@ -372,9 +371,8 @@ bool CSTRModel::configure(IParameterProvider& paramProvider)
372371

373372
if (_qsReactionBulk != nullptr)
374373
{
375-
//_dynReactionBulk->fillConservedMoietiesBulk(_MconvMoityBulk, _QsCompBulk); // fill conserved moities matrix
376374
_QsCompBulk.resize(_nComp);
377-
_dynReactionBulk->fillConservedMoietiesBulk2(_MconvMoityBulk2, _nConservedQuants, _QsCompBulk);
375+
_dynReactionBulk->fillConservedMoietiesBulk(_MconvMoityBulk, _nConservedQuants, _QsCompBulk);
378376

379377

380378
bool hasQSBinding = false;
@@ -395,8 +393,10 @@ bool CSTRModel::configure(IParameterProvider& paramProvider)
395393
}
396394
else
397395
{
398-
_MconvMoityBulk = Eigen::MatrixXd::Zero(0, 0); // matrix for conserved moities
396+
_MconvMoityBulk = Eigen::Matrix<active, Eigen::Dynamic, Eigen::Dynamic>::Zero(0, 0); // matrix for conserved moities
399397
_QsCompBulk.clear();
398+
_qsReactionBulk = nullptr;
399+
_nConservedQuants = 0;
400400
}
401401

402402
for (unsigned int type = 0; type < _nParType; ++type)
@@ -708,8 +708,8 @@ void CSTRModel::consistentInitialState(const SimulationTime& simTime, double* co
708708
if (mIdx >= _nConservedQuants)
709709
continue;
710710
double dotProduct = 0.0;
711-
for (unsigned int j = 0; j < _MconvMoityBulk2.cols(); ++j)
712-
dotProduct += static_cast<double>(_MconvMoityBulk2(i, j)) * c[j];
711+
for (unsigned int j = 0; j < _MconvMoityBulk.cols(); ++j)
712+
dotProduct += static_cast<double>(_MconvMoityBulk(i, j)) * c[j];
713713

714714
conservedQuants[mIdx] = dotProduct;
715715
mIdx++;
@@ -778,12 +778,12 @@ void CSTRModel::consistentInitialState(const SimulationTime& simTime, double* co
778778
continue;
779779

780780
int jIdx = 0;
781-
for (unsigned int j = 0; j < _MconvMoityBulk2.cols(); ++j)
781+
for (unsigned int j = 0; j < _MconvMoityBulk.cols(); ++j)
782782
{
783783
if (_QsCompBulk[j] == 0)
784784
continue;
785785

786-
mat.native(mIdx, jIdx) = static_cast<double>(_MconvMoityBulk2(i, j));
786+
mat.native(mIdx, jIdx) = static_cast<double>(_MconvMoityBulk(i, j));
787787
jIdx++;
788788
}
789789
mIdx++;
@@ -817,12 +817,12 @@ void CSTRModel::consistentInitialState(const SimulationTime& simTime, double* co
817817
continue;
818818

819819
int jIdx = 0;
820-
for (unsigned int j = 0; j < _MconvMoityBulk2.cols(); ++j)
820+
for (unsigned int j = 0; j < _MconvMoityBulk.cols(); ++j)
821821
{
822822
if(_QsCompBulk[j] == 0)
823823
continue;
824824

825-
mat.native(mIdx, jIdx) = static_cast<double>(_MconvMoityBulk2(i, j));
825+
mat.native(mIdx, jIdx) = static_cast<double>(_MconvMoityBulk(i, j));
826826
jIdx++;
827827
}
828828
mIdx++;
@@ -855,25 +855,17 @@ void CSTRModel::consistentInitialState(const SimulationTime& simTime, double* co
855855
continue;
856856
int jIdx = 0;
857857
double dotProduct = 0.0;
858-
for (unsigned int j = 0; j < _MconvMoityBulk2.cols(); ++j)
858+
for (unsigned int j = 0; j < _MconvMoityBulk.cols(); ++j)
859859
{
860860
if (_QsCompBulk[j] == 0)
861861
continue;
862-
dotProduct += static_cast<double>(_MconvMoityBulk2(i, j)) * x[jIdx];
862+
dotProduct += static_cast<double>(_MconvMoityBulk(i, j)) * x[jIdx];
863863
jIdx++;
864864
}
865865
r[mIdx] = dotProduct - conservedQuants[mIdx];
866866
mIdx++;
867867
}
868868

869-
std::cout << "Residual: " << std::endl;
870-
for (unsigned int i = 0; i < probSize; ++i) {
871-
std::cout << r[i] << std::endl;
872-
}
873-
std::cout << "Solution: " << std::endl;
874-
for (unsigned int i = 0; i < probSize; ++i)
875-
std::cout << x[i] << std::endl;
876-
877869
return true;
878870
},
879871
jacFunc, errorTol, static_cast<double*>(solution), nonlinMem, jacobianMatrix, probSize);
@@ -1566,88 +1558,7 @@ int CSTRModel::residual(const SimulationTime& simTime, const ConstSimulationStat
15661558
{
15671559
return residualImpl<double, double, double, false>(simTime.t, simTime.secIdx, simState.vecStateY, simState.vecStateYdot, res, threadLocalMem.get());
15681560
}
1569-
template <typename StateType, typename ResidualType, typename ParamType, bool wantJac>
1570-
void CSTRModel::applyConservedMoitiesBulk(double t, unsigned int secIdx, const ColumnPosition& colPos, StateType const* const c, double const* const yDot, ResidualType* const resC, LinearBufferAllocator tlmAlloc)
1571-
{
1572-
const ParamType flowIn = static_cast<ParamType>(_flowRateIn);
1573-
const ParamType flowOut = static_cast<ParamType>(_flowRateOut);
1574-
double const* const cDot = yDot ? yDot + _nComp : nullptr;
1575-
const double vDot = yDot ? yDot[2 * _nComp + _totalBound] : 0.0;
1576-
1577-
LinearBufferAllocator subAlloc = tlmAlloc.manageRemainingMemory();
1578-
1579-
BufferedArray<ResidualType> temp = subAlloc.array<ResidualType>(_nComp);
1580-
Eigen::Map<Eigen::Vector<ResidualType, Eigen::Dynamic>> resCWithMoities(static_cast<ResidualType*>(temp), _nComp);
1581-
resCWithMoities.setZero();
1582-
1583-
BufferedArray<ResidualType> temp2 = subAlloc.array<ResidualType>(_nComp);
1584-
Eigen::Map<Eigen::Vector<ResidualType, Eigen::Dynamic>> qsFlux(static_cast<ResidualType*>(temp2), _dynReactionBulk->numReactionsLiquid());
1585-
qsFlux.setZero();
1586-
1587-
_dynReactionBulk->computeQuasiStationaryReactionFlux(t, secIdx, colPos, c, qsFlux, _qsReactionBulk, subAlloc);
1588-
Eigen::Map<Eigen::Vector<ResidualType, Eigen::Dynamic>> mapResC(resC, _nComp);
1589-
1590-
int MoityIdx = 0;
1591-
int qsreac = 0;
1592-
int comp = 0;
1593-
1594-
for (unsigned int state = 0; state < _nComp; state++)
1595-
{
1596-
if (_stateMap[state].test(0)) // dynamic
1597-
{
1598-
resCWithMoities[state] = resC[state];
1599-
if (wantJac)
1600-
{
1601-
_jac.native(state, state) += (static_cast<double>(vDot) + static_cast<double>(flowOut)); // dF_{ci}/dcj = v_liquidDot + F_out
1602-
_dynReactionBulk->analyticJacobianLiquidSingleFluxAdd(t, secIdx, colPos, reinterpret_cast<double const*>(c), state, 0, _jac.row(state), subAlloc);
1603-
}
1604-
comp++;
1605-
}
1606-
else if (_stateMap[state].test(1)) // conserved
1607-
{
1608-
// todo test of matrix times vector isfaster
1609-
ResidualType dotProduct = 0.0;
1610-
for (unsigned int i = 0; i < _MconvMoityBulk.cols(); ++i)
1611-
{
1612-
dotProduct += _MconvMoityBulk(MoityIdx, i) * (mapResC[i]);
1613-
if (wantJac)
1614-
{
1615-
_jac.native(state, i) += (static_cast<double>(vDot) + static_cast<double>(flowOut)) * _MconvMoityBulk(MoityIdx, i); // dF_{ci}/dcj = v_liquidDot + F_out
1616-
if (cadet_likely(yDot))
1617-
_jac.native(i, _nComp + _totalBound) += _MconvMoityBulk(MoityIdx, i) * cDot[i] ; // dF/dvliquid = cDot
1618-
}
1619-
}
1620-
if (wantJac)
1621-
_dynReactionBulk->analyticJacobianLiquidSingleFluxAdd(t, secIdx, colPos, reinterpret_cast<double const*>(c), state, 0, _jac.row(state), subAlloc);
1622-
resCWithMoities[state] = dotProduct;
1623-
MoityIdx++;
1624-
comp++;
1625-
}
1626-
else if (_stateMap[state].test(2)) // algebraic
1627-
{
1628-
resCWithMoities[state] = qsFlux[qsreac];
16291561

1630-
if (wantJac)
1631-
{
1632-
_dynReactionBulk->analyticJacobianQuasiStationaryReaction(t, secIdx, colPos, reinterpret_cast<double const*>(c), state, qsreac, _jac.row(state), subAlloc);
1633-
_jac.native(state, _nComp + _totalBound) = 0.0; // dF_{ci}/dvliquid = 0
1634-
}
1635-
qsreac++;
1636-
}
1637-
}
1638-
mapResC = resCWithMoities;
1639-
1640-
std::cout << "Jacobian with conserved moities" << std::endl;
1641-
std::cout << "Jacobian after" << std::endl;
1642-
for (int i = 0; i < _nComp + 1; i++)
1643-
{
1644-
for (int j = 0; j < _nComp + 1; j++)
1645-
{
1646-
std::cout << _jac.native(i, j) << " ";
1647-
}
1648-
std::cout << std::endl;
1649-
}
1650-
}
16511562
template <typename ResidualType>
16521563
void CSTRModel::EigenMatrixTimesDemseMatrix(Eigen::Matrix<ResidualType, Eigen::Dynamic, Eigen::Dynamic> A, linalg::DenseMatrix& B)
16531564
{
@@ -1699,7 +1610,7 @@ void CSTRModel::EigenMatrixTimesDemseMatrix(Eigen::Matrix<ResidualType, Eigen::D
16991610
}
17001611

17011612
template <typename StateType, typename ResidualType, typename ParamType, bool wantJac>
1702-
void CSTRModel::applyConservedMoitiesBulk2(double t, unsigned int secIdx, const ColumnPosition& colPos, StateType const* const c, double const* const yDot, ResidualType* const resC, LinearBufferAllocator tlmAlloc)
1613+
void CSTRModel::applyConservedMoitiesBulk(double t, unsigned int secIdx, const ColumnPosition& colPos, StateType const* const c, double const* const yDot, ResidualType* const resC, LinearBufferAllocator tlmAlloc)
17031614
{
17041615
// prepare memory
17051616
LinearBufferAllocator subAlloc = tlmAlloc.manageRemainingMemory();
@@ -1714,17 +1625,12 @@ void CSTRModel::applyConservedMoitiesBulk2(double t, unsigned int secIdx, const
17141625

17151626
_dynReactionBulk->computeQuasiStationaryReactionFlux(t, secIdx, colPos, c, qsFlux, _qsReactionBulk, subAlloc);
17161627

1717-
// calculate conserved moities
1718-
Eigen::Matrix<ResidualType, Eigen::Dynamic, Eigen::Dynamic> M(_nComp, _nComp);
1719-
_dynReactionBulk->fillConservedMoietiesBulk2(M, _nConservedQuants, _QsCompBulk); // fill conserved moities matrix (alternative method)
1720-
1721-
17221628
// buffer memory for transformed residual
17231629
BufferedArray<ResidualType> temp = subAlloc.array<ResidualType>(_nComp);
17241630
Eigen::Map<Eigen::Vector<ResidualType, Eigen::Dynamic>> resCWithMoities(static_cast<ResidualType*>(temp), _nComp);
17251631
resCWithMoities.setZero();
17261632

1727-
Eigen::Matrix<ResidualType, Eigen::Dynamic, Eigen::Dynamic> MconvMoityBulk2Cast = _MconvMoityBulk2.template cast<ResidualType>();
1633+
Eigen::Matrix<ResidualType, Eigen::Dynamic, Eigen::Dynamic> MconvMoityBulk2Cast = _MconvMoityBulk.template cast<ResidualType>();
17281634

17291635
// multiply conserved moities matrix with residual
17301636
resCWithMoities = MconvMoityBulk2Cast * mapResC;
@@ -1867,7 +1773,7 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons
18671773

18681774
if (_hasQuasiStationaryReactionBulk)
18691775
{
1870-
applyConservedMoitiesBulk2<StateType, ResidualType, ParamType, wantJac>(t, secIdx, colPos, c, yDot, resC, subAlloc);
1776+
applyConservedMoitiesBulk<StateType, ResidualType, ParamType, wantJac>(t, secIdx, colPos, c, yDot, resC, subAlloc);
18711777

18721778
}
18731779
}
@@ -2401,11 +2307,10 @@ void CSTRModel::addTimeDerivativeJacobian(double t, double alpha, const ConstSim
24012307
continue;
24022308
else
24032309
{
2404-
//mat.native(state, state) += timeV * (static_cast<double>(_MconvMoityBulk(MoityIdx, MoityIdx))); // dRes / dcDot
24052310
for (int j = 0; j < _nComp; j++)
24062311
{
2407-
mat.native(i, j) += timeV * static_cast<double>(_MconvMoityBulk2(i, j)); // dRes / dcDot
2408-
mat.native(i, _nComp + _totalBound) += alpha * static_cast<double>(_MconvMoityBulk2(i, j)) * c[j]; // dRes / dVlDot
2312+
mat.native(i, j) += timeV * static_cast<double>(_MconvMoityBulk(i, j)); // dRes / dcDot
2313+
mat.native(i, _nComp + _totalBound) += alpha * static_cast<double>(_MconvMoityBulk(i, j)) * c[j]; // dRes / dVlDot
24092314
}
24102315
MoityIdx++;
24112316
}
@@ -2475,7 +2380,7 @@ void CSTRModel::addTimeDerivativeJacobian(double t, double alpha, const ConstSim
24752380

24762381
/*if (_hasQuasiStationaryReactionBulk)
24772382
{
2478-
EigenMatrixTimesDemseMatrix(_MconvMoityBulk2, mat);
2383+
EigenMatrixTimesDemseMatrix(_MconvMoityBulk, mat);
24792384
}*/
24802385

24812386
std::cout << "Jacobian Derivative: " << std::endl;

‎src/libcadet/model/StirredTankModel.hpp

Lines changed: 1 addition & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -85,8 +85,6 @@ class CSTRModel : public UnitOperationBase
8585

8686
template <typename StateType, typename ResidualType, typename ParamType, bool wantJac>
8787
void applyConservedMoitiesBulk(double t, unsigned int secIdx, const ColumnPosition& colPos, StateType const* const y, double const* const yDot, ResidualType* const resC, LinearBufferAllocator tlmAlloc);
88-
template <typename StateType, typename ResidualType, typename ParamType, bool wantJac>
89-
void applyConservedMoitiesBulk2(double t, unsigned int secIdx, const ColumnPosition& colPos, StateType const* const y, double const* const yDot, ResidualType* const resC, LinearBufferAllocator tlmAlloc);
9088

9189
virtual int jacobian(const SimulationTime& simTime, const ConstSimulationState& simState, double* const res, const AdJacobianParams& adJac, util::ThreadLocalStorage& threadLocalMem);
9290
virtual int residualWithJacobian(const SimulationTime& simTime, const ConstSimulationState& simState, double* const res, const AdJacobianParams& adJac, util::ThreadLocalStorage& threadLocalMem);
@@ -103,7 +101,6 @@ class CSTRModel : public UnitOperationBase
103101

104102
virtual void consistentInitialState(const SimulationTime& simTime, double* const vecStateY, const AdJacobianParams& adJac, double errorTol, util::ThreadLocalStorage& threadLocalMem);
105103
virtual void consistentInitialTimeDerivative(const SimulationTime& simTime, double const* vecStateY, double* const vecStateYdot, util::ThreadLocalStorage& threadLocalMem);
106-
107104
virtual void initializeSensitivityStates(const std::vector<double*>& vecSensY) const;
108105
virtual void consistentInitialSensitivity(const SimulationTime& simTime, const ConstSimulationState& simState,
109106
std::vector<double*>& vecSensY, std::vector<double*>& vecSensYdot, active const* const adRes, util::ThreadLocalStorage& threadLocalMem);
@@ -182,13 +179,11 @@ class CSTRModel : public UnitOperationBase
182179

183180
IDynamicReactionModel* _dynReactionBulk; //!< Dynamic reactions in the bulk volume
184181

185-
Eigen::MatrixXd _MconvMoityBulk; //!< Matrix with conservation of moieties in the bulk volume
186-
Eigen::Matrix<active, Eigen::Dynamic, Eigen::Dynamic> _MconvMoityBulk2; //!< Matrix with conservation of moieties in the bulk volume
182+
Eigen::Matrix<active, Eigen::Dynamic, Eigen::Dynamic> _MconvMoityBulk; //!< Matrix with conservation of moieties in the bulk volume
187183
int const* _qsReactionBulk; //!< Indices of reactions that are not conserved in the bulk volume
188184
bool _hasQuasiStationaryReactionBulk; //!< Flag that determines whether there are quasi-stationary reactions in the bulk volume
189185
std::vector<int> _QsCompBulk; //!< Indices of components that are conserved in the bulk volume
190186
int _nConservedQuants;
191-
std::vector<std::bitset<3>> _stateMap; // !< Bitset that tracks the properties of the state vector (dynamic, moities, algebraic)
192187

193188
class Exporter : public ISolutionExporter
194189
{

0 commit comments

Comments
 (0)