Skip to content

Commit 1d033f5

Browse files
AntoniaBergerAntoniaBerger
committed
WIP undependent of old calculation
1 parent 4a1f39d commit 1d033f5

8 files changed

+236
-327
lines changed

‎src/libcadet/model/ReactionModel.hpp

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -265,8 +265,8 @@ class IDynamicReactionModel
265265
* @param [in,out] QsCompBulk Vector to be filled with the components involved in qs reactions
266266
*/
267267
virtual void fillConservedMoietiesBulk(Eigen::MatrixXd& M, std::vector<int>& QsCompBulk) {}
268-
virtual void fillConservedMoietiesBulk2(Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>& M, int& QsCompBulk) {}
269-
virtual void fillConservedMoietiesBulk2(Eigen::Matrix<active, Eigen::Dynamic, Eigen::Dynamic>& M, 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) {}
270270

271271

272272
/**

‎src/libcadet/model/StirredTankModel.cpp

Lines changed: 20 additions & 80 deletions
Original file line numberDiff line numberDiff line change
@@ -372,33 +372,10 @@ bool CSTRModel::configure(IParameterProvider& paramProvider)
372372

373373
if (_qsReactionBulk != nullptr)
374374
{
375-
_dynReactionBulk->fillConservedMoietiesBulk(_MconvMoityBulk, _QsCompBulk); // fill conserved moities matrix
376-
_dynReactionBulk->fillConservedMoietiesBulk2(_MconvMoityBulk2, _nConservedQuants);
377-
int nMoitiesBulk = _MconvMoityBulk.rows();
378-
if (nMoitiesBulk != 0)
379-
{
380-
int mIdx = 0;
381-
for (int state = 0; state < _nComp; state++)
382-
{
383-
std::bitset<3> c;
384-
if (_QsCompBulk[state] == 0) // state is dynamic
385-
{
386-
c.set(0);
387-
_stateMap.push_back(c);
388-
}
389-
else if (mIdx < nMoitiesBulk) // state is dynamic and calculated with conserved moities
390-
{
391-
c.set(1);
392-
_stateMap.push_back(c);
393-
mIdx++;
394-
}
395-
else if (mIdx >= nMoitiesBulk) // state algebraic
396-
{
397-
c.set(2);
398-
_stateMap.push_back(c);
399-
}
400-
}
401-
}
375+
//_dynReactionBulk->fillConservedMoietiesBulk(_MconvMoityBulk, _QsCompBulk); // fill conserved moities matrix
376+
_QsCompBulk.resize(_nComp);
377+
_dynReactionBulk->fillConservedMoietiesBulk2(_MconvMoityBulk2, _nConservedQuants, _QsCompBulk);
378+
402379

403380
bool hasQSBinding = false;
404381
for (int i = 0; i < _nParType; i++)
@@ -732,7 +709,7 @@ void CSTRModel::consistentInitialState(const SimulationTime& simTime, double* co
732709
continue;
733710
double dotProduct = 0.0;
734711
for (unsigned int j = 0; j < _MconvMoityBulk2.cols(); ++j)
735-
dotProduct += _MconvMoityBulk2(i, j) * c[j];
712+
dotProduct += static_cast<double>(_MconvMoityBulk2(i, j)) * c[j];
736713

737714
conservedQuants[mIdx] = dotProduct;
738715
mIdx++;
@@ -795,14 +772,19 @@ void CSTRModel::consistentInitialState(const SimulationTime& simTime, double* co
795772
int mIdx = 0;
796773
for (unsigned int i = 0; i < _nComp; ++i)
797774
{
798-
if (_QsCompBulk[i]== 0)
775+
if (_QsCompBulk[i] == 0)
799776
continue;
800777
if (mIdx >= _nConservedQuants)
801778
continue;
802-
for (unsigned int j = 0; j < _MconvMoityBulk2.cols(); ++i)
779+
780+
int jIdx = 0;
781+
for (unsigned int j = 0; j < _MconvMoityBulk2.cols(); ++j)
803782
{
804-
mat.native(mIdx, j) = _MconvMoityBulk2(i, j);
805-
j++;
783+
if (_QsCompBulk[j] == 0)
784+
continue;
785+
786+
mat.native(mIdx, jIdx) = static_cast<double>(_MconvMoityBulk2(i, j));
787+
jIdx++;
806788
}
807789
mIdx++;
808790
}
@@ -840,7 +822,7 @@ void CSTRModel::consistentInitialState(const SimulationTime& simTime, double* co
840822
if(_QsCompBulk[j] == 0)
841823
continue;
842824

843-
mat.native(mIdx, jIdx) = _MconvMoityBulk2(i, j);
825+
mat.native(mIdx, jIdx) = static_cast<double>(_MconvMoityBulk2(i, j));
844826
jIdx++;
845827
}
846828
mIdx++;
@@ -877,7 +859,7 @@ void CSTRModel::consistentInitialState(const SimulationTime& simTime, double* co
877859
{
878860
if (_QsCompBulk[j] == 0)
879861
continue;
880-
dotProduct += _MconvMoityBulk2(i, j) * x[jIdx];
862+
dotProduct += static_cast<double>(_MconvMoityBulk2(i, j)) * x[jIdx];
881863
jIdx++;
882864
}
883865
r[mIdx] = dotProduct - conservedQuants[mIdx];
@@ -1734,17 +1716,18 @@ void CSTRModel::applyConservedMoitiesBulk2(double t, unsigned int secIdx, const
17341716

17351717
// calculate conserved moities
17361718
Eigen::Matrix<ResidualType, Eigen::Dynamic, Eigen::Dynamic> M(_nComp, _nComp);
1737-
_dynReactionBulk->fillConservedMoietiesBulk2(M, _nConservedQuants); // fill conserved moities matrix (alternative method)
1719+
_dynReactionBulk->fillConservedMoietiesBulk2(M, _nConservedQuants, _QsCompBulk); // fill conserved moities matrix (alternative method)
17381720

17391721

17401722
// buffer memory for transformed residual
17411723
BufferedArray<ResidualType> temp = subAlloc.array<ResidualType>(_nComp);
17421724
Eigen::Map<Eigen::Vector<ResidualType, Eigen::Dynamic>> resCWithMoities(static_cast<ResidualType*>(temp), _nComp);
17431725
resCWithMoities.setZero();
17441726

1727+
Eigen::Matrix<ResidualType, Eigen::Dynamic, Eigen::Dynamic> MconvMoityBulk2Cast = _MconvMoityBulk2.template cast<ResidualType>();
17451728

17461729
// multiply conserved moities matrix with residual
1747-
resCWithMoities = M * mapResC;
1730+
resCWithMoities = MconvMoityBulk2Cast * mapResC;
17481731

17491732
// add quasi stationary reaction to residium
17501733
const int nQsReac = _dynReactionBulk->numReactionQuasiStationary();
@@ -1768,7 +1751,7 @@ void CSTRModel::applyConservedMoitiesBulk2(double t, unsigned int secIdx, const
17681751

17691752
if (wantJac)
17701753
{
1771-
EigenMatrixTimesDemseMatrix(M, _jac);
1754+
EigenMatrixTimesDemseMatrix(MconvMoityBulk2Cast, _jac);
17721755
int mIdx = 0;
17731756
int rIdx = 0;
17741757
for (int i = 0; i < _nComp; i++)
@@ -1787,52 +1770,9 @@ void CSTRModel::applyConservedMoitiesBulk2(double t, unsigned int secIdx, const
17871770

17881771
}
17891772

1790-
//std::cout << "Jacobian with conserved moities" << std::endl;
1791-
//std::cout << "Jacobian before" << std::endl;
1792-
//for (int i = 0; i < _nComp + 1; i++)
1793-
//{
1794-
// for (int j = 0; j < _nComp+ 1; j++)
1795-
// {
1796-
// std::cout << _jac.native(i, j) << " ";
1797-
// }
1798-
// std::cout << std::endl;
1799-
//}
1800-
// multiply conserved moities matrix with jacobian
1801-
//if (wantJac)
1802-
//{
1803-
// EigenMatrixTimesDemseMatrix(M, _jac);
1804-
1805-
// int state = _nComp - _dynReactionBulk->numReactionQuasiStationary();
1806-
// for (int qsReac = 0; qsReac < nQsReac; ++qsReac) // todo this in a function
1807-
// {
1808-
// if (state < _nComp)
1809-
// {
1810-
// _dynReactionBulk->analyticJacobianQuasiStationaryReaction(t, secIdx, colPos, reinterpret_cast<double const*>(c), state, qsReac, _jac.row(state), subAlloc);
1811-
// _jac.native(state, _nComp + _totalBound) = 0.0; // dF_{ci}/dvliquid = 0
1812-
// state++;
1813-
// }
1814-
// else
1815-
// throw InvalidParameterException(
1816-
// "Jacobian implementation with conserved moities: Too many quasi stationary reactions detected. "
1817-
// "Please check the implementation of the model."
1818-
// );
1819-
// }
1820-
1821-
//std::cout << "Jacobian with conserved moities" << std::endl;
1822-
//std::cout << "Jacobian after" << std::endl;
1823-
//for (int i = 0; i < _nComp+ 1; i++)
1824-
//{
1825-
// for (int j = 0; j < _nComp+1; j++)
1826-
//{
1827-
// std::cout << _jac.native(i, j) << " ";
1828-
// }
1829-
//std::cout << std::endl;
1830-
//}
1831-
18321773
}
18331774

18341775

1835-
18361776
template <typename StateType, typename ResidualType, typename ParamType, bool wantJac>
18371777
int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* const y, double const* const yDot, ResidualType* const res, LinearBufferAllocator tlmAlloc)
18381778
{

‎src/libcadet/model/StirredTankModel.hpp

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -31,7 +31,7 @@
3131

3232
namespace cadet
3333
{
34-
struct ColumnPosition;
34+
struct ColumnPosition;
3535

3636
namespace model
3737
{
@@ -183,7 +183,7 @@ class CSTRModel : public UnitOperationBase
183183
IDynamicReactionModel* _dynReactionBulk; //!< Dynamic reactions in the bulk volume
184184

185185
Eigen::MatrixXd _MconvMoityBulk; //!< Matrix with conservation of moieties in the bulk volume
186-
Eigen::MatrixXd _MconvMoityBulk2; //!< 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
187187
int const* _qsReactionBulk; //!< Indices of reactions that are not conserved in the bulk volume
188188
bool _hasQuasiStationaryReactionBulk; //!< Flag that determines whether there are quasi-stationary reactions in the bulk volume
189189
std::vector<int> _QsCompBulk; //!< Indices of components that are conserved in the bulk volume

‎src/libcadet/model/reaction/CrystallizationReaction.cpp

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -198,14 +198,14 @@ namespace cadet
198198
virtual bool requiresConfiguration() const CADET_NOEXCEPT { return true; }
199199
virtual bool usesParamProviderInDiscretizationConfig() const CADET_NOEXCEPT { return false; }
200200
template <typename RowIterator>
201-
void jacobianQuasiSteadyLiquidImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, int reaction, const RowIterator& jac, LinearBufferAllocator workSpace) const { }
201+
void jacobianQuasiStationaryBulkImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, int reaction, const RowIterator& jac, LinearBufferAllocator workSpace) const { }
202202

203203
template <typename RowIterator>
204204
void jacobianSingleFluxImpl(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, int state, int reaction, const RowIterator& jac, LinearBufferAllocator workSpace) const { }
205205

206206
template<typename StateType, typename ResidualType>
207207
int quasiStationaryFlux(double t, unsigned int secIdx, const ColumnPosition& colPos, StateType const* y,
208-
Eigen::Map<Eigen::Vector<ResidualType, Eigen::Dynamic>> fluxes, int const* mapQSReac, LinearBufferAllocator workSpace){return 0;}
208+
Eigen::Map<Eigen::Vector<ResidualType, Eigen::Dynamic>> fluxes, LinearBufferAllocator workSpace){return 0;}
209209

210210
virtual int const* reactionQuasiStationarity() const CADET_NOEXCEPT { return nullptr; }
211211
virtual void timeDerivativeQuasiStationaryReaction(double t, unsigned int secIdx, const ColumnPosition& colPos, double const* y, double* dReacDt, LinearBufferAllocator workSpace){ }
@@ -347,8 +347,8 @@ namespace cadet
347347
virtual unsigned int numReactionsLiquid() const CADET_NOEXCEPT { return 1; }
348348
virtual unsigned int numReactionsCombined() const CADET_NOEXCEPT { return 1; }
349349
virtual unsigned int numReactionQuasiStationary() const CADET_NOEXCEPT { return 0; }
350-
template <typename ResidualType>
351-
void fillConservedMoietiesBulk21(Eigen::Matrix<ResidualType, Eigen::Dynamic, Eigen::Dynamic>& M, int& QsCompBulk) {}
350+
template<typename ResidualType>
351+
void fillConservedMoietiesBulk(Eigen::Matrix<ResidualType, Eigen::Dynamic, Eigen::Dynamic>& M, int& conservedState, std::vector<int>& QsCompBulk) {}
352352

353353
CADET_DYNAMICREACTIONMODEL_BOILERPLATE
354354

‎src/libcadet/model/reaction/DummyReaction.cpp

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -129,8 +129,8 @@ class DummyDynamicReaction : public IDynamicReactionModel
129129
virtual unsigned int numReactionsLiquid() const CADET_NOEXCEPT { return 0; }
130130
virtual unsigned int numReactionsCombined() const CADET_NOEXCEPT { return 0; }
131131
virtual unsigned int numReactionQuasiStationary() const CADET_NOEXCEPT { return 0; }
132-
template <typename ResidualType>
133-
void fillConservedMoietiesBulk21(Eigen::Matrix<ResidualType, Eigen::Dynamic, Eigen::Dynamic>& M, int& QsCompBulk) {}
132+
template<typename ResidualType>
133+
void fillConservedMoietiesBulk(Eigen::Matrix<ResidualType, Eigen::Dynamic, Eigen::Dynamic>& M, int& conservedState, std::vector<int>& QsCompBulk) {}
134134
protected:
135135
};
136136

0 commit comments

Comments
 (0)