wip fix reaction tests

Author: jbreue16

src/libcadet/model/ColumnModel1D.hpp

@@ -501,7 +501,7 @@ class ColumnModel1D : public UnitOperationBase
 			}
 		}
 
-		// particle jacobian (including film diffusion, isotherm and time derivative)
+		// particle jacobian (including film diffusion, isotherm, reaction and time derivative)
 		for (int colNode = 0; colNode < _disc.nPoints; colNode++)
 		{
 			for (int type = 0; type < _disc.nParType; type++)
@@ -519,7 +519,7 @@ class ColumnModel1D : public UnitOperationBase
 		// inlet and bulk jacobian
 		_convDispOp.calcTransportJacobian(_globalJac, _jacInlet, idxr.offsetC());
 
-		// particle jacobian (without isotherm, which is handled in residualKernel)
+		// particle transport jacobian
 		for (int colNode = 0; colNode < _disc.nPoints; colNode++)
 		{
 			for (int parType = 0; parType < _disc.nParType; parType++)

src/libcadet/model/particle/GeneralRateParticle.cpp

@@ -162,6 +162,7 @@ namespace model
 				paramProvider.popScope();
 			}
 
+			MultiplexedScopeSelector scopeGuard(paramProvider, "reaction", _dynReaction->usesParamProviderInDiscretizationConfig());
 			reactionConfSuccess = _dynReaction->configureModelDiscretization(paramProvider, _nComp, _nBound.get(), _parDiffOp->offsetBoundComp()) && reactionConfSuccess;
 		}
 

src/libcadet/model/particle/HomogeneousParticle.cpp

@@ -176,6 +176,7 @@ namespace model
 				paramProvider.popScope();
 			}
 
+			MultiplexedScopeSelector scopeGuard(paramProvider, "reaction", _dynReaction->usesParamProviderInDiscretizationConfig());
 			reactionConfSuccess = _dynReaction->configureModelDiscretization(paramProvider, _nComp, _nBound.get(), _boundOffset) && reactionConfSuccess;
 		}
 
@@ -404,24 +405,22 @@ namespace model
 
 	void HomogeneousParticle::setParJacPattern(std::vector<Eigen::Triplet<double>>& tripletList, const unsigned int offsetPar, const unsigned int offsetBulk, unsigned int colNode, unsigned int secIdx) const
 	{
-		// binding Jacobian pattern
-		// add dense nComp * nBound blocks, since all solid and liquid entries can be coupled through binding.
+		// binding and reaction Jacobian pattern
+		// add dense nComp * nBound blocks, since all solid and liquid entries can be coupled through reaction and or binding.
 		for (unsigned int parState = 0; parState < _nComp + _strideBound; parState++) {
 			for (unsigned int toParState = 0; toParState < _nComp + _strideBound; toParState++) {
 				tripletList.push_back(Eigen::Triplet<double>(offsetPar + parState, offsetPar + toParState, 0.0));
 			}
 		}
 
 		// flux Jacobian pattern
-
 		for (unsigned int comp = 0; comp < _nComp; comp++)
 		{
 			tripletList.push_back(Eigen::Triplet<double>(offsetPar + comp, offsetBulk + comp, 0.0));
 			tripletList.push_back(Eigen::Triplet<double>(offsetBulk + comp, offsetPar + comp, 0.0));
 		}
 	}
 
-
 	unsigned int HomogeneousParticle::jacobianNNZperParticle() const
 	{
 		return (_nComp + _strideBound) * (_nComp + _strideBound) + _nComp * 4; // reaction, binding patter + film diffusion pattern for one particle

src/libcadet/model/parts/ParticleDiffusionOperatorDG.cpp

@@ -912,12 +912,13 @@ namespace parts
 	// ========================================						DG particle Jacobian							=============================================  //
 	// ==========================================================================================================================================================  //
 
-	void ParticleDiffusionOperatorDG::parBindingPattern(std::vector<Eigen::Triplet<double>>& tripletList, const int offset, const unsigned int colNode)
+	void ParticleDiffusionOperatorDG::parBindingAndReactionPattern(std::vector<Eigen::Triplet<double>>& tripletList, const int offset, const unsigned int colNode)
 	{
-		// every bound state might depend on every bound and liquid state
 		for (int parNode = 0; parNode < _nParPoints; parNode++)
 		{
-			for (int state = 0; state < strideParPoint(); state++)
+			// binding: every bound state might depend on every bound and liquid state
+			// reaction: every liquid and bound state might depend on every liquid and bound state
+			for (int state = 0; state < _strideBound; state++)
 			{
 				for (int stateDep = 0; stateDep < strideParPoint(); stateDep++) {
 					// row: jump over previous nodes and liquid states and add current bound state offset
@@ -928,14 +929,14 @@ namespace parts
 			}
 		}
 	}
+
 	/**
 	 * @brief calculates the particle dispersion jacobian Pattern, including entries for the dependence of particle entries on bulk entries through film diffusion boundary condition
 	 * @detail Does NOT add film diffusion entries for the dependence of bulk conc. on particle conc.
 	*/
 	void ParticleDiffusionOperatorDG::setParticleJacobianPattern(std::vector<ParticleDiffusionOperatorDG::T>& tripletList, unsigned int offsetPar, unsigned int offsetBulk, unsigned int colNode, unsigned int secIdx)
 	{
-
-		parBindingPattern(tripletList, offsetPar, colNode);
+		parBindingAndReactionPattern(tripletList, offsetPar, colNode);
 
 		// Ordering of particle surface diffusion:
 		// bnd0comp0, bnd1comp0, bnd0comp1, bnd1comp1, bnd0comp2, bnd1comp2

src/libcadet/model/parts/ParticleDiffusionOperatorDG.hpp

@@ -143,7 +143,7 @@ namespace parts
 
 	protected:
 
-		void parBindingPattern(std::vector<Eigen::Triplet<double>>& tripletList, const int offset, const unsigned int colNode);
+		void parBindingAndReactionPattern(std::vector<Eigen::Triplet<double>>& tripletList, const int offset, const unsigned int colNode);
 
 		template <typename StateType, typename ResidualType, typename ParamType, bool wantJac, bool wantRes>
 		int residualImpl(double t, unsigned int secIdx, StateType const* yPar, StateType const* yBulk, double const* yDotPar, ResidualType* resPar, linalg::BandedEigenSparseRowIterator& jacBase);

test/ReactionModelTests.cpp

@@ -221,13 +221,18 @@ namespace reaction
 
 				std::vector<int> nBound = jpp.getIntArray("NBOUND");
 				const int nTotalBound = std::accumulate(nBound.begin(), nBound.end(), 0);
-				jpp.set("REACTION_MODEL", "MASS_ACTION_LAW");
 
 				std::string scope;
 				if (!isCSTR)
+				{
+					jpp.set("REACTION_MODEL", "MASS_ACTION_LAW");
 					scope = "reaction";
+				}
 				else
+				{
+					jpp.set("REACTION_MODEL_PARTICLE", "MASS_ACTION_LAW");
 					scope = "reaction_particle";
+				}
 
 				jpp.addScope(scope);
 				auto gs2 = util::makeGroupScope(jpp, scope);