@@ -709,7 +709,7 @@ void CSTRModel::consistentInitialState(const SimulationTime& simTime, double* co
709
709
continue ;
710
710
double dotProduct = 0.0 ;
711
711
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];
713
713
714
714
conservedQuants[mIdx ] = dotProduct;
715
715
mIdx ++;
@@ -772,14 +772,19 @@ void CSTRModel::consistentInitialState(const SimulationTime& simTime, double* co
772
772
int mIdx = 0 ;
773
773
for (unsigned int i = 0 ; i < _nComp; ++i)
774
774
{
775
- if (_QsCompBulk[i]== 0 )
775
+ if (_QsCompBulk[i] == 0 )
776
776
continue ;
777
777
if (mIdx >= _nConservedQuants)
778
778
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)
780
782
{
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++;
783
788
}
784
789
mIdx ++;
785
790
}
@@ -817,7 +822,7 @@ void CSTRModel::consistentInitialState(const SimulationTime& simTime, double* co
817
822
if (_QsCompBulk[j] == 0 )
818
823
continue ;
819
824
820
- mat.native (mIdx , jIdx) = _MconvMoityBulk2 (i, j);
825
+ mat.native (mIdx , jIdx) = static_cast < double >( _MconvMoityBulk2 (i, j) );
821
826
jIdx++;
822
827
}
823
828
mIdx ++;
@@ -854,7 +859,7 @@ void CSTRModel::consistentInitialState(const SimulationTime& simTime, double* co
854
859
{
855
860
if (_QsCompBulk[j] == 0 )
856
861
continue ;
857
- dotProduct += _MconvMoityBulk2 (i, j) * x[jIdx];
862
+ dotProduct += static_cast < double >( _MconvMoityBulk2 (i, j) ) * x[jIdx];
858
863
jIdx++;
859
864
}
860
865
r[mIdx ] = dotProduct - conservedQuants[mIdx ];
@@ -1719,9 +1724,10 @@ void CSTRModel::applyConservedMoitiesBulk2(double t, unsigned int secIdx, const
1719
1724
Eigen::Map<Eigen::Vector<ResidualType, Eigen::Dynamic>> resCWithMoities (static_cast <ResidualType*>(temp), _nComp);
1720
1725
resCWithMoities.setZero ();
1721
1726
1727
+ Eigen::Matrix<ResidualType, Eigen::Dynamic, Eigen::Dynamic> MconvMoityBulk2Cast = _MconvMoityBulk2.template cast <ResidualType>();
1722
1728
1723
1729
// multiply conserved moities matrix with residual
1724
- resCWithMoities = M * mapResC;
1730
+ resCWithMoities = MconvMoityBulk2Cast * mapResC;
1725
1731
1726
1732
// add quasi stationary reaction to residium
1727
1733
const int nQsReac = _dynReactionBulk->numReactionQuasiStationary ();
@@ -1745,7 +1751,7 @@ void CSTRModel::applyConservedMoitiesBulk2(double t, unsigned int secIdx, const
1745
1751
1746
1752
if (wantJac)
1747
1753
{
1748
- EigenMatrixTimesDemseMatrix (M , _jac);
1754
+ EigenMatrixTimesDemseMatrix (MconvMoityBulk2Cast , _jac);
1749
1755
int mIdx = 0 ;
1750
1756
int rIdx = 0 ;
1751
1757
for (int i = 0 ; i < _nComp; i++)
@@ -1764,52 +1770,9 @@ void CSTRModel::applyConservedMoitiesBulk2(double t, unsigned int secIdx, const
1764
1770
1765
1771
}
1766
1772
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
-
1809
1773
}
1810
1774
1811
1775
1812
-
1813
1776
template <typename StateType, typename ResidualType, typename ParamType, bool wantJac>
1814
1777
int CSTRModel::residualImpl (double t, unsigned int secIdx, StateType const * const y, double const * const yDot, ResidualType* const res, LinearBufferAllocator tlmAlloc)
1815
1778
{
0 commit comments