add multiplyMatrixWithDerivativeJacobian as an alternative for multip…

…lyWithDerivativeJacobian

src/libcadet/model/StirredTankModel.cpp

@@ -1791,6 +1791,77 @@ void CSTRModel::multiplyWithDerivativeJacobian(const SimulationTime& simTime, co
 	r[_nComp + _totalBound] = s[_nComp + _totalBound];
 }
 
+
+void CSTRModel::multiplyMatrixWithDerivativeJacobian(const ConstSimulationState& simState, double const* sDot, double* ret)
+{	
+
+	// Handle inlet DOFs (all algebraic -> dRes/dcDot = 0)
+	std::fill_n(ret, _nComp, 0.0);
+	double* const r = ret + _nComp;
+
+
+	// create zero matrix
+	// create here a zero matrix with the same size as the Jacobian
+	// ActiveDenseMatrix matrix = ActiveDenseMatrix(_nComp + _totalBound + 1, _nComp + _totalBound + 1);
+
+	addTimeDerivativeJacobian(0.0, 1.0, simState, matrix);
+
+	// Handle actual ODE DOFs
+	double* const mat = matrix + _nComp;
+	double const* const s = sDot + _nComp;
+
+	// Concentrations: V^l * \dot{c}_i \dot{V^l} * c_i + V^s * ([sum_j sum_m d_j \dot{q}_{j,i,m}]) - c_{in,i} * F_in + c_i * F_out == 0
+	for (unsigned int i = 0; i < _nComp; ++i)
+	{
+		r[i] = mat.native(i, i) *  s[i]; // dRes / dcDot
+
+		double qSum = 0.0;
+		for (unsigned int type = 0; type < _nParType; ++type)
+		{
+			double const* const qi = q + _offsetParType[type] + _boundOffset[type * _nComp + i];
+			const unsigned int localOffset = _nComp + _offsetParType[type] + _boundOffset[type * _nComp + i];
+			const double vSolidParVolFrac = vSolid * static_cast<double>(_parTypeVolFrac[type]);
+			double qSumType = 0.0;
+			for (unsigned int j = 0; j < _nBound[type * _nComp + i]; ++j)
+			{
+				r[i] += native(i, localOffset + j) * s[localOffset + j]; // dRes / d_qDot
+				// + _nComp: Moves over liquid phase components
+				// + _offsetParType[type]: Moves to particle type
+				// + _boundOffset[i]: Moves over bound states of previous components
+				// + j: Moves to current bound state j of component i
+			}
+		}
+		r[i] += mat.native(i, _nComp + _totalBound) *  s[_nComp + _totalBound]; // dRes / dvLiquidDot
+	}
+
+	// Bound states
+	for (unsigned int type = 0; type < _nParType; ++type)
+	{
+		// Jump to solid phase of current particle type
+		double* const rQ = r + _nComp + _offsetParType[type];
+		double const* const sQ = s + _nComp + _offsetParType[type];
+		std::fill_n(rQ, _strideBound[type], 0.0);
+
+		// Skip binding models without dynamic binding fluxes
+		IBindingModel* const binding = _binding[type];
+		if (!binding->hasDynamicReactions())
+			continue;
+
+		int const* const qsReaction = binding->reactionQuasiStationarity();
+		for (unsigned int idx = 0; idx < _strideBound[type]; ++idx)
+		{
+			// Skip quasi-stationary fluxes
+			if (qsReaction[idx])
+				continue;
+
+			rQ[idx] = mat.native(idx, idx) *sQ[idx];
+		}
+	}
+
+	// Volume: \dot{V} - F_{in} + F_{out} + F_{filter} == 0
+	r[_nComp + _totalBound] = mat.native(_nComp + _totalBound, _nComp + _totalBound) * s[_nComp + _totalBound];
+}
+
 int CSTRModel::linearSolve(double t, double alpha, double tol, double* const rhs, double const* const weight,
 	const ConstSimulationState& simState)
 {

src/libcadet/model/StirredTankModel.hpp

@@ -108,6 +108,7 @@ class CSTRModel : public UnitOperationBase
 
 	virtual void multiplyWithJacobian(const SimulationTime& simTime, const ConstSimulationState& simState, double const* yS, double alpha, double beta, double* ret);
 	virtual void multiplyWithDerivativeJacobian(const SimulationTime& simTime, const ConstSimulationState& simState, double const* sDot, double* ret);
+	void multiplyMatrixWithDerivativeJacobian(const ConstSimulationState& simState, double const* sDot, double* ret);
 
 	virtual bool hasInlet() const CADET_NOEXCEPT { return true; }
 	virtual bool hasOutlet() const CADET_NOEXCEPT { return true; }

test/CSTR-Residual.cpp

@@ -227,7 +227,7 @@ inline void checkTimeDerivativeJacobianFD(double flowRateIn, double flowRateOut,
 	cstr->notifyDiscontinuousSectionTransition(0.0, 0u, {y.data(), yDot.data()}, cadet::AdJacobianParams{nullptr, nullptr, 0u});
 
 	// Compare Jacobians
-	cadet::test::compareTimeDerivativeJacobianFD(cstr, cstr, y.data(), yDot.data(), jacDir.data(), jacCol1.data(), jacCol2.data(), tls);
+	cadet::test::compareTimeDerivativeJacobianFDCSTR(cstr, cstr, y.data(), yDot.data(), jacDir.data(), jacCol1.data(), jacCol2.data(), tls);
 
 	mb->destroyUnitOperation(cstr);
 	destroyModelBuilder(mb);

test/JacobianHelper.hpp

@@ -521,6 +521,15 @@ inline void compareTimeDerivativeJacobianFD(cadet::IUnitOperation* modelA, cadet
 		yDot, dir, colA, colB, modelA->numDofs(), modelA->numDofs(), h, absTol, relTol);
 }
 
+inline void compareTimeDerivativeJacobianFDCSTR(cadet::IUnitOperation* modelA, cadet::IUnitOperation* modelB, double const* y, double const* yDot, double* dir, double* colA, double* colB,
+	cadet::util::ThreadLocalStorage& tls, double h = 1e-6, double absTol = 0.0, double relTol = std::numeric_limits<float>::epsilon() * 100.0)
+{
+	compareJacobianFD(
+		[=, &tls](double const* lDir, double* res) -> void { modelA->residual(SimulationTime{ 0.0, 0u }, ConstSimulationState{ y, lDir }, res, tls); },
+		[=](double const* lDir, double* res) -> void { modelB->multiplyMatrixWithDerivativeJacobian(ConstSimulationState{ y, yDot },lDir, ,res); },
+		yDot, dir, colA, colB, modelA->numDofs(), modelA->numDofs(), h, absTol, relTol);
+}
+
 } // namespace test
 } // namespace cadet