@@ -709,7 +709,7 @@ void CSTRModel::consistentInitialState(const SimulationTime& simTime, double* co
709709 continue ;
710710 double dotProduct = 0.0 ;
711711 for (unsigned int j = 0 ; j < _MconvMoityBulk2.cols (); ++j)
712- dotProduct += _MconvMoityBulk2 (i, j) * c[j];
712+ dotProduct += static_cast < double >( _MconvMoityBulk2 (i, j) ) * c[j];
713713
714714 conservedQuants[mIdx ] = dotProduct;
715715 mIdx ++;
@@ -772,14 +772,19 @@ void CSTRModel::consistentInitialState(const SimulationTime& simTime, double* co
772772 int mIdx = 0 ;
773773 for (unsigned int i = 0 ; i < _nComp; ++i)
774774 {
775- if (_QsCompBulk[i]== 0 )
775+ if (_QsCompBulk[i] == 0 )
776776 continue ;
777777 if (mIdx >= _nConservedQuants)
778778 continue ;
779- for (unsigned int j = 0 ; j < _MconvMoityBulk2.cols (); ++i)
779+
780+ int jIdx = 0 ;
781+ for (unsigned int j = 0 ; j < _MconvMoityBulk2.cols (); ++j)
780782 {
781- mat.native (mIdx , j) = _MconvMoityBulk2 (i, j);
782- j++;
783+ if (_QsCompBulk[j] == 0 )
784+ continue ;
785+
786+ mat.native (mIdx , jIdx) = static_cast <double >(_MconvMoityBulk2 (i, j));
787+ jIdx++;
783788 }
784789 mIdx ++;
785790 }
@@ -817,7 +822,7 @@ void CSTRModel::consistentInitialState(const SimulationTime& simTime, double* co
817822 if (_QsCompBulk[j] == 0 )
818823 continue ;
819824
820- mat.native (mIdx , jIdx) = _MconvMoityBulk2 (i, j);
825+ mat.native (mIdx , jIdx) = static_cast < double >( _MconvMoityBulk2 (i, j) );
821826 jIdx++;
822827 }
823828 mIdx ++;
@@ -854,7 +859,7 @@ void CSTRModel::consistentInitialState(const SimulationTime& simTime, double* co
854859 {
855860 if (_QsCompBulk[j] == 0 )
856861 continue ;
857- dotProduct += _MconvMoityBulk2 (i, j) * x[jIdx];
862+ dotProduct += static_cast < double >( _MconvMoityBulk2 (i, j) ) * x[jIdx];
858863 jIdx++;
859864 }
860865 r[mIdx ] = dotProduct - conservedQuants[mIdx ];
@@ -1719,9 +1724,10 @@ void CSTRModel::applyConservedMoitiesBulk2(double t, unsigned int secIdx, const
17191724 Eigen::Map<Eigen::Vector<ResidualType, Eigen::Dynamic>> resCWithMoities (static_cast <ResidualType*>(temp), _nComp);
17201725 resCWithMoities.setZero ();
17211726
1727+ Eigen::Matrix<ResidualType, Eigen::Dynamic, Eigen::Dynamic> MconvMoityBulk2Cast = _MconvMoityBulk2.template cast <ResidualType>();
17221728
17231729 // multiply conserved moities matrix with residual
1724- resCWithMoities = M * mapResC;
1730+ resCWithMoities = MconvMoityBulk2Cast * mapResC;
17251731
17261732 // add quasi stationary reaction to residium
17271733 const int nQsReac = _dynReactionBulk->numReactionQuasiStationary ();
@@ -1745,7 +1751,7 @@ void CSTRModel::applyConservedMoitiesBulk2(double t, unsigned int secIdx, const
17451751
17461752 if (wantJac)
17471753 {
1748- EigenMatrixTimesDemseMatrix (M , _jac);
1754+ EigenMatrixTimesDemseMatrix (MconvMoityBulk2Cast , _jac);
17491755 int mIdx = 0 ;
17501756 int rIdx = 0 ;
17511757 for (int i = 0 ; i < _nComp; i++)
@@ -1764,52 +1770,9 @@ void CSTRModel::applyConservedMoitiesBulk2(double t, unsigned int secIdx, const
17641770
17651771 }
17661772
1767- // std::cout << "Jacobian with conserved moities" << std::endl;
1768- // std::cout << "Jacobian before" << std::endl;
1769- // for (int i = 0; i < _nComp + 1; i++)
1770- // {
1771- // for (int j = 0; j < _nComp+ 1; j++)
1772- // {
1773- // std::cout << _jac.native(i, j) << " ";
1774- // }
1775- // std::cout << std::endl;
1776- // }
1777- // multiply conserved moities matrix with jacobian
1778- // if (wantJac)
1779- // {
1780- // EigenMatrixTimesDemseMatrix(M, _jac);
1781-
1782- // int state = _nComp - _dynReactionBulk->numReactionQuasiStationary();
1783- // for (int qsReac = 0; qsReac < nQsReac; ++qsReac) // todo this in a function
1784- // {
1785- // if (state < _nComp)
1786- // {
1787- // _dynReactionBulk->analyticJacobianQuasiStationaryReaction(t, secIdx, colPos, reinterpret_cast<double const*>(c), state, qsReac, _jac.row(state), subAlloc);
1788- // _jac.native(state, _nComp + _totalBound) = 0.0; // dF_{ci}/dvliquid = 0
1789- // state++;
1790- // }
1791- // else
1792- // throw InvalidParameterException(
1793- // "Jacobian implementation with conserved moities: Too many quasi stationary reactions detected. "
1794- // "Please check the implementation of the model."
1795- // );
1796- // }
1797-
1798- // std::cout << "Jacobian with conserved moities" << std::endl;
1799- // std::cout << "Jacobian after" << std::endl;
1800- // for (int i = 0; i < _nComp+ 1; i++)
1801- // {
1802- // for (int j = 0; j < _nComp+1; j++)
1803- // {
1804- // std::cout << _jac.native(i, j) << " ";
1805- // }
1806- // std::cout << std::endl;
1807- // }
1808-
18091773}
18101774
18111775
1812-
18131776template <typename StateType, typename ResidualType, typename ParamType, bool wantJac>
18141777int CSTRModel::residualImpl (double t, unsigned int secIdx, StateType const * const y, double const * const yDot, ResidualType* const res, LinearBufferAllocator tlmAlloc)
18151778{
0 commit comments