add moity matrix to cstr

src/libcadet/model/StirredTankModel.cpp

@@ -75,6 +75,7 @@ CSTRModel::~CSTRModel() CADET_NOEXCEPT
 	delete[] _nBound;
 	delete[] _strideBound;
 	delete[] _offsetParType;
+	delete[] _tmp;
 
 	delete _dynReactionBulk;
 }
@@ -151,6 +152,14 @@ bool CSTRModel::configureModelDiscretization(IParameterProvider& paramProvider,
 		_offsetParType = nullptr;
 	}
 
+	if (paramProvider.exists("QuaseStationaryReactions"))
+	{
+		_tmp = new active[_nComp]; // temp
+		int nMoities = _nComp - _dynReactionBulk->numQSReactionLiquid();
+		_MconvMoityBulk = Eigen::MatrixXd::Zero(nMoities, _nComp);
+		_dynReactionBulk->fillConservedMoietiesBulk(_MconvMoityBulk, _qsComponentsBulk);
+	}
+
 	if ((_nParType > 1) && (_totalBound == 0))
 		throw InvalidParameterException("Multiple particle types but not a single bound state detected");
 
@@ -1309,39 +1318,36 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons
 
 		if (_dynReactionBulk->hasQuasiStationaryReactionsLiquid())
 		{
-			int nmoities = _nComp - _dynReactionBulk->numQSReactionLiquid();
-			Eigen::MatrixXd MconMoiBulk = Eigen::MatrixXd::Zero(nmoities, _nComp);
-			std::vector<int> qsComponent;
-			_dynReactionBulk->fillConservedMoietiesBulk(MconMoiBulk, qsComponent);
-
-			Eigen::VectorXd resCMoities = Eigen::VectorXd::Zero(_nComp);
-			resCMoities = Eigen::Map<Eigen::VectorXd>(resC, _nComp);
-
-
+			Eigen::Map<Vector<ResidualType, Dynamic>> resCMoities(_tmp, _nComp);
+			//Eigen::Map<Vector<double, Dynamic>> resCMoities(reinterpret_cast<double>(_tmp),_nComp);
+			resCMoities.setZero();
+
+			MoityIdx = 0
 			for (unsigned int comp = 0; comp < _nComp; ++comp)
 			{
-				if (qsComponent[comp])
+				if (_qsComponentsBulk[comp])
 				{
-					resCMoities[comp] = MconMoiBulk.row(comp) * resC[comp];
+					resCMoities[comp] = _MconvMoityBulk.row(MoityIdx) * resC[comp];
+					MoityIdx++;
 				}
-				if (comp < _nComp - nmoities)
+				else if (comp < _nComp - nmoities)
 				{
 					resCMoities[comp] += v * flux[comp];
 				}
-				if(comp > _nComp - nmoities)
+				else if(comp > _nComp - nmoities)
 				{
 					resCMoities[comp] = v * flux[comp];
 				}
 			}
+			Eigen::Map<Vector<ResidualType, Dynamic>> mapResC(resC, _nComp);
+			mapResC = resCMoities;
 		}
 		else
 		{
 			for (unsigned int comp = 0; comp < _nComp; ++comp)
 				resC[comp] += v * flux[comp];
 		}
 
-		//AB resC = resCMoities;
-
 		if (wantJac)
 		{
 			for (unsigned int comp = 0; comp < _nComp; ++comp)
@@ -1374,7 +1380,7 @@ int CSTRModel::residualImpl(double t, unsigned int secIdx, StateType const* cons
 						continue;
 
 					// Add time derivative to solid phase
-					resQ[idx] += qDot[idx]; //AB why here "+" ?
+					resQ[idx] += qDot[idx];
 				}
 			}
 		}

src/libcadet/model/StirredTankModel.hpp

@@ -169,7 +169,11 @@ class CSTRModel : public UnitOperationBase
 	std::vector<active> _initConditions; //!< Initial conditions, ordering: Liquid phase concentration, solid phase concentration, liquid volume
 	std::vector<double> _initConditionsDot; //!< Initial conditions for time derivative
 
+	active* _tmp;
+
 	IDynamicReactionModel* _dynReactionBulk; //!< Dynamic reactions in the bulk volume
+	Eigen::MatrixXd _MconvMoityBulk; //!<  Matrix with conservation of moieties in the bulk volume
+	std::vector<int> _qsComponentsBulk; //!< Indices of components that are not conserved in the bulk volume
 
 	class Exporter : public ISolutionExporter
 	{

src/libcadet/model/reaction/MassActionLawReaction.cpp

@@ -867,12 +867,12 @@ class MassActionLawReactionBase : public DynamicReactionModelBase
 		}
 
 		// Entferne Nullzeilen in-place
-		std::vector<int> nonZeroRowIndices;
+		std::vector<int> nonQSComponents;
 
 		for (Eigen::Index i = 0; i < QSS.rows(); ++i)
 		{
 			if (QSS.row(i).norm() > 1e-10) {
-				nonZeroRowIndices.push_back(i);
+				nonQSComponents.push_back(i);
 			}
 			else
 			{
@@ -881,21 +881,16 @@ class MassActionLawReactionBase : public DynamicReactionModelBase
 		}
 
 		// Redimensioniere QSS und kopiere nur die nicht-null Zeilen
-		if (!nonZeroRowIndices.empty())
+		if (!nonQSComponents.empty())
 		{
-			Eigen::MatrixXd QSSCompressed(nonZeroRowIndices.size(), QSS.cols());
-			for (std::size_t i = 0; i < nonZeroRowIndices.size(); ++i)
+			Eigen::MatrixXd QSSCompressed(nonQSComponents.size(), QSS.cols());
+			for (std::size_t i = 0; i < nonQSComponents.size(); ++i)
 			{
-				QSSCompressed.row(i) = QSS.row(nonZeroRowIndices[i]);
+				QSSCompressed.row(i) = QSS.row(nonQSComponents[i]);
 			}
-			QSS.swap(QSSCompressed); // 
-		}
-		else
-		{
-			QSS.resize(0, countQS); // Leere Matrix, falls alle Zeilen Null sind
+			QSS.swap(QSSCompressed); // Swap the compressed matrix back into QSS
 		}
 
-
 		//3. Test if the matrix is full rank
 		int rang = QSS.fullPivLu().rank();
 		if (rang != countQS)
@@ -906,8 +901,19 @@ class MassActionLawReactionBase : public DynamicReactionModelBase
 		//3. Calculate the null space of the matrix
 		Eigen::MatrixXd leftZeroSpace = QSS.transpose().fullPivLu().kernel().transpose();
 
+		if (!nonZeroRowIndices.empty())
+		{
+			M = leftZeroSpace;
+		}
+		else // add zero rows for each component that is not in the quasi stationary
+		{
+			for (int idx : nonQSComponents) {
+				leftZeroSpace.conservativeResize(NoChange, leftZeroSpace.cols() + 1);
+				leftZeroSpace.rightCols(leftZeroSpace.cols() - idx - 1) = leftZeroSpace.middleCols(idx, leftZeroSpace.cols() - idx - 1);
+				leftZeroSpace.col(idx).setZero();
+			}
+		}
 		M = leftZeroSpace;
-
 	}
 };