Modularize particle DG discretization (#406)

Author: jbreue16

src/libcadet/CMakeLists.txt

@@ -214,6 +214,7 @@ if (ENABLE_DG)
 	list (APPEND LIBCADET_SOURCES
 		${CMAKE_SOURCE_DIR}/src/libcadet/model/parts/DGToolbox.cpp
 		${CMAKE_SOURCE_DIR}/src/libcadet/model/parts/ConvectionDispersionOperatorDG.cpp
+		${CMAKE_SOURCE_DIR}/src/libcadet/model/parts/ParticleDiffusionOperatorDG.cpp
 		${CMAKE_SOURCE_DIR}/src/libcadet/model/LumpedRateModelWithPoresDG.cpp
 		${CMAKE_SOURCE_DIR}/src/libcadet/model/LumpedRateModelWithoutPoresDG.cpp
 		${CMAKE_SOURCE_DIR}/src/libcadet/model/GeneralRateModelDG.cpp

src/libcadet/model/GeneralRateModelDG-InitialConditions.cpp

@@ -386,8 +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 = static_cast<double>(_disc.deltaR[type]) * std::floor(node / _disc.nParNode[type])
-					+ 0.5 * static_cast<double>(_disc.deltaR[type]) * (1 + _disc.parNodes[type][node % _disc.nParNode[type]]);
+				const double r = _parDiffOp.relativeCoordinate(type, node);
 				const ColumnPosition colPos{ z, 0.0, r };
 
 				// Determine whether nonlinear solver is required
@@ -398,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>(_parPorosity[type]);
-				linalg::conservedMoietiesFromPartitionedMask(mask, _disc.nBound + type * _disc.nComp, _disc.nComp, qShell - _disc.nComp, conservedQuants, static_cast<double>(_parPorosity[type]), epsQ);
+				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);
 
 				std::function<bool(double const* const, linalg::detail::DenseMatrixBase&)> jacFunc;
 
@@ -478,7 +477,7 @@ void GeneralRateModelDG::consistentInitialState(const SimulationTime& simTime, d
 								continue;
 							}
 
-							mat.native(rIdx, rIdx) = static_cast<double>(_parPorosity[type]);
+							mat.native(rIdx, rIdx) = static_cast<double>(_parDiffOp._parPorosity[type]);
 
 							for (unsigned int bnd = 0; bnd < _disc.nBound[_disc.nComp * type + comp]; ++bnd, ++bndIdx)
 							{
@@ -526,7 +525,7 @@ void GeneralRateModelDG::consistentInitialState(const SimulationTime& simTime, d
 								continue;
 							}
 
-							mat.native(rIdx, rIdx) = static_cast<double>(_parPorosity[type]);
+							mat.native(rIdx, rIdx) = static_cast<double>(_parDiffOp._parPorosity[type]);
 
 							for (unsigned int bnd = 0; bnd < _disc.nBound[_disc.nComp * type + comp]; ++bnd, ++bndIdx)
 							{
@@ -573,7 +572,7 @@ void GeneralRateModelDG::consistentInitialState(const SimulationTime& simTime, d
 								continue;
 							}
 
-							r[rIdx] = static_cast<double>(_parPorosity[type]) * x[rIdx] - conservedQuants[rIdx];
+							r[rIdx] = static_cast<double>(_parDiffOp._parPorosity[type]) * x[rIdx] - conservedQuants[rIdx];
 
 							for (unsigned int bnd = 0; bnd < _disc.nBound[_disc.nComp * type + comp]; ++bnd, ++bndIdx)
 							{
@@ -714,9 +713,8 @@ 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 = (static_cast<double>(_disc.deltaR[type]) * std::floor(j / _disc.nParNode[type])
-					+ 0.5 * static_cast<double>(_disc.deltaR[type]) * (1 + _disc.parNodes[type][j % _disc.nParNode[type]]))
-					/ (static_cast<double>(_parRadius[type]) - static_cast<double>(_parCoreRadius[type]));
+				const double r = _parDiffOp.relativeCoordinate(type, j);
+
 				_binding[type]->timeDerivativeQuasiStationaryFluxes(simTime.t, simTime.secIdx,
 					ColumnPosition{ z, 0.0, r },
 					qShellDot - _disc.nComp, qShellDot, dFluxDt, tlmAlloc);
@@ -799,7 +797,7 @@ 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(_parDiffusionMode) || isSectionDependent(_parSurfDiffusionMode))
+	if (isSectionDependent(_parDiffOp._parDiffusionMode) || isSectionDependent(_parDiffOp._parSurfDiffusionMode))
 		LOG(Warning) << "Lean consistent initialization is not appropriate for section-dependent pore and surface diffusion";
 
 	BENCH_SCOPE(_timerConsistentInit);
@@ -832,8 +830,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 = static_cast<double>(_disc.deltaR[type]) * std::floor(shell / _disc.nParNode[type])
-						+ 0.5 * static_cast<double>(_disc.deltaR[type]) * (1 + _disc.parNodes[type][shell % _disc.nParNode[type]]);
+					const double r = _parDiffOp.relativeCoordinate(type, 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)
@@ -843,7 +840,6 @@ void GeneralRateModelDG::leanConsistentInitialState(const SimulationTime& simTim
 			} CADET_PARFOR_END;
 		}
 	}
-
 }
 
 /**
@@ -889,7 +885,7 @@ 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(_parDiffusionMode) || isSectionDependent(_parSurfDiffusionMode))
+	if (isSectionDependent(_parDiffOp._parDiffusionMode) || isSectionDependent(_parDiffOp._parSurfDiffusionMode))
 		LOG(Warning) << "Lean consistent initialization is not appropriate for section-dependent pore and surface diffusion";
 
 	BENCH_SCOPE(_timerConsistentInit);
@@ -1272,7 +1268,7 @@ 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(_parDiffusionMode) || isSectionDependent(_parSurfDiffusionMode))
+	if (isSectionDependent(_parDiffOp._parDiffusionMode) || isSectionDependent(_parDiffOp._parSurfDiffusionMode))
 		LOG(Warning) << "Lean consistent initialization is not appropriate for section-dependent pore and surface diffusion";
 
 	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>(_parPorosity[parType]), _disc.nComp, _disc.nBound + _disc.nComp * parType,
+	parts::cell::addTimeDerivativeToJacobianParticleShell<linalg::BandedEigenSparseRowIterator, true>(jac, alpha, static_cast<double>(_parDiffOp._parPorosity[parType]), _disc.nComp, _disc.nBound + _disc.nComp * parType,
 		_poreAccessFactor.data() + _disc.nComp * parType, _disc.strideBound[parType], _disc.boundOffset + _disc.nComp * parType, _binding[parType]->reactionQuasiStationarity());
 }