Refactor particle operator to be particle type specific (#423)

Author: jbreue16

src/libcadet/model/GeneralRateModelDG-InitialConditions.cpp

@@ -386,7 +386,7 @@ void GeneralRateModelDG::consistentInitialState(const SimulationTime& simTime, d
 				active* const localAdY = adJac.adY ? adJac.adY + localOffsetToParticle + localOffsetInParticle : nullptr;
 
 				// r (particle) coordinate of current node
-				const double r = _parDiffOp.relativeCoordinate(type, node);
+				const double r = _parDiffOp[type].relativeCoordinate(node);
 				const ColumnPosition colPos{ z, 0.0, r };
 
 				// Determine whether nonlinear solver is required
@@ -397,8 +397,8 @@ void GeneralRateModelDG::consistentInitialState(const SimulationTime& simTime, d
 				linalg::selectVectorSubset(qShell - _disc.nComp, mask, solution);
 
 				// Save values of conserved moieties
-				const double epsQ = 1.0 - static_cast<double>(_parDiffOp._parPorosity[type]);
-				linalg::conservedMoietiesFromPartitionedMask(mask, _disc.nBound + type * _disc.nComp, _disc.nComp, qShell - _disc.nComp, conservedQuants, static_cast<double>(_parDiffOp._parPorosity[type]), epsQ);
+				const double epsQ = 1.0 - static_cast<double>(_parDiffOp[type]._parPorosity);
+				linalg::conservedMoietiesFromPartitionedMask(mask, _disc.nBound + type * _disc.nComp, _disc.nComp, qShell - _disc.nComp, conservedQuants, static_cast<double>(_parDiffOp[type]._parPorosity), epsQ);
 
 				std::function<bool(double const* const, linalg::detail::DenseMatrixBase&)> jacFunc;
 
@@ -477,7 +477,7 @@ void GeneralRateModelDG::consistentInitialState(const SimulationTime& simTime, d
 								continue;
 							}
 
-							mat.native(rIdx, rIdx) = static_cast<double>(_parDiffOp._parPorosity[type]);
+							mat.native(rIdx, rIdx) = static_cast<double>(_parDiffOp[type]._parPorosity);
 
 							for (unsigned int bnd = 0; bnd < _disc.nBound[_disc.nComp * type + comp]; ++bnd, ++bndIdx)
 							{
@@ -525,7 +525,7 @@ void GeneralRateModelDG::consistentInitialState(const SimulationTime& simTime, d
 								continue;
 							}
 
-							mat.native(rIdx, rIdx) = static_cast<double>(_parDiffOp._parPorosity[type]);
+							mat.native(rIdx, rIdx) = static_cast<double>(_parDiffOp[type]._parPorosity);
 
 							for (unsigned int bnd = 0; bnd < _disc.nBound[_disc.nComp * type + comp]; ++bnd, ++bndIdx)
 							{
@@ -572,7 +572,7 @@ void GeneralRateModelDG::consistentInitialState(const SimulationTime& simTime, d
 								continue;
 							}
 
-							r[rIdx] = static_cast<double>(_parDiffOp._parPorosity[type]) * x[rIdx] - conservedQuants[rIdx];
+							r[rIdx] = static_cast<double>(_parDiffOp[type]._parPorosity) * x[rIdx] - conservedQuants[rIdx];
 
 							for (unsigned int bnd = 0; bnd < _disc.nBound[_disc.nComp * type + comp]; ++bnd, ++bndIdx)
 							{
@@ -713,7 +713,7 @@ void GeneralRateModelDG::consistentInitialTimeDerivative(const SimulationTime& s
 			if (_binding[type]->dependsOnTime())
 			{
 				// r (particle) coordinate of current node (particle radius normed to 1) - needed in externally dependent adsorption kinetic
-				const double r = _parDiffOp.relativeCoordinate(type, j);
+				const double r = _parDiffOp[type].relativeCoordinate(j);
 
 				_binding[type]->timeDerivativeQuasiStationaryFluxes(simTime.t, simTime.secIdx,
 					ColumnPosition{ z, 0.0, r },
@@ -797,10 +797,13 @@ void GeneralRateModelDG::consistentInitialTimeDerivative(const SimulationTime& s
  */
 void GeneralRateModelDG::leanConsistentInitialState(const SimulationTime& simTime, double* const vecStateY, const AdJacobianParams& adJac, double errorTol, util::ThreadLocalStorage& threadLocalMem)
 {
-	if (isSectionDependent(_parDiffOp._parDiffusionMode) || isSectionDependent(_parDiffOp._parSurfDiffusionMode))
-		LOG(Warning) << "Lean consistent initialization is not appropriate for section-dependent pore and surface diffusion";
+	for (int parType = 0; parType < _disc.nParType; parType++)
+	{
+		if (isSectionDependent(_parDiffOp[parType]._parDiffusionMode) || isSectionDependent(_parDiffOp[parType]._parSurfDiffusionMode))
+			LOG(Warning) << "Lean consistent initialization is not appropriate for section-dependent pore and surface diffusion";
 
-	BENCH_SCOPE(_timerConsistentInit);
+		BENCH_SCOPE(_timerConsistentInit);
+	}
 
 	Indexer idxr(_disc);
 
@@ -830,7 +833,7 @@ void GeneralRateModelDG::leanConsistentInitialState(const SimulationTime& simTim
 					const int localOffsetInParticle = static_cast<int>(shell) * idxr.strideParNode(type) + idxr.strideParLiquid();
 					double* const qShell = vecStateY + localOffsetToParticle + localOffsetInParticle;
 					// r (particle) coordinate of current node
-					const double r = _parDiffOp.relativeCoordinate(type, shell);
+					const double r = _parDiffOp[type].relativeCoordinate(shell);
 					const ColumnPosition colPos{ z, 0.0, r};
 
 					// Perform consistent initialization that does not require a full fledged nonlinear solver (that may fail or damage the current state vector)
@@ -885,8 +888,13 @@ void GeneralRateModelDG::leanConsistentInitialState(const SimulationTime& simTim
  */
 void GeneralRateModelDG::leanConsistentInitialTimeDerivative(double t, double const* const vecStateY, double* const vecStateYdot, double* const res, util::ThreadLocalStorage& threadLocalMem)
 {
-	if (isSectionDependent(_parDiffOp._parDiffusionMode) || isSectionDependent(_parDiffOp._parSurfDiffusionMode))
-		LOG(Warning) << "Lean consistent initialization is not appropriate for section-dependent pore and surface diffusion";
+	for (int parType = 0; parType < _disc.nParType; parType++)
+	{
+		if (isSectionDependent(_parDiffOp[parType]._parDiffusionMode) || isSectionDependent(_parDiffOp[parType]._parSurfDiffusionMode))
+			LOG(Warning) << "Lean consistent initialization is not appropriate for section-dependent pore and surface diffusion";
+
+		BENCH_SCOPE(_timerConsistentInit);
+	}
 
 	BENCH_SCOPE(_timerConsistentInit);
 
@@ -1268,8 +1276,13 @@ void GeneralRateModelDG::solveBulkTimeDerivativeSystem(const SimulationTime& sim
 void GeneralRateModelDG::leanConsistentInitialSensitivity(const SimulationTime& simTime, const ConstSimulationState& simState,
 	std::vector<double*>& vecSensY, std::vector<double*>& vecSensYdot, active const* const adRes, util::ThreadLocalStorage& threadLocalMem)
 {
-	if (isSectionDependent(_parDiffOp._parDiffusionMode) || isSectionDependent(_parDiffOp._parSurfDiffusionMode))
-		LOG(Warning) << "Lean consistent initialization is not appropriate for section-dependent pore and surface diffusion";
+	for (int parType = 0; parType < _disc.nParType; parType++)
+	{
+		if (isSectionDependent(_parDiffOp[parType]._parDiffusionMode) || isSectionDependent(_parDiffOp[parType]._parSurfDiffusionMode))
+			LOG(Warning) << "Lean consistent initialization is not appropriate for section-dependent pore and surface diffusion";
+
+		BENCH_SCOPE(_timerConsistentInit);
+	}
 
 	BENCH_SCOPE(_timerConsistentInit);
 

src/libcadet/model/GeneralRateModelDG-LinearSolver.cpp

@@ -220,7 +220,7 @@ void GeneralRateModelDG::assembleDiscretizedGlobalJacobian(double alpha, Indexer
  */
 void GeneralRateModelDG::addTimeDerivativeToJacobianParticleShell(linalg::BandedEigenSparseRowIterator& jac, const Indexer& idxr, double alpha, unsigned int parType)
 {
-	parts::cell::addTimeDerivativeToJacobianParticleShell<linalg::BandedEigenSparseRowIterator, true>(jac, alpha, static_cast<double>(_parDiffOp._parPorosity[parType]), _disc.nComp, _disc.nBound + _disc.nComp * parType,
+	parts::cell::addTimeDerivativeToJacobianParticleShell<linalg::BandedEigenSparseRowIterator, true>(jac, alpha, static_cast<double>(_parDiffOp[parType]._parPorosity), _disc.nComp, _disc.nBound + _disc.nComp * parType,
 		_poreAccessFactor.data() + _disc.nComp * parType, _disc.strideBound[parType], _disc.boundOffset + _disc.nComp * parType, _binding[parType]->reactionQuasiStationarity());
 }