Add FV particle diffusion operator (#465)

Also makes field `SPATIAL_METHOD` in particle discretization mandatory

Author: jbreue16

.github/workflows/ci.yml

@@ -224,6 +224,7 @@ jobs:
           ${BUILD_DIR}/test/testRunner [CI_sens22]
           ${BUILD_DIR}/test/testRunner [CI_sens23]
           ${BUILD_DIR}/test/testRunner [CI_sens24]
+          ${BUILD_DIR}/test/testRunner [CI_sens25]
       - name: Run CI test set III - crystallization
         run: |
           set -e

.github/workflows/coverage.yml

@@ -87,6 +87,7 @@ jobs:
           ${BUILD_DIR}/test/testRunner [CI_sens22]
           ${BUILD_DIR}/test/testRunner [CI_sens23]
           ${BUILD_DIR}/test/testRunner [CI_sens24]
+          ${BUILD_DIR}/test/testRunner [CI_sens25]
       - name: Run CI test set III - crystallization
         run: |
           ${BUILD_DIR}/test/testRunner [CI_cry1]

src/libcadet/CMakeLists.txt

@@ -221,6 +221,7 @@ if (ENABLE_DG)
 		${CMAKE_SOURCE_DIR}/src/libcadet/model/particle/HomogeneousParticle.cpp
 		${CMAKE_SOURCE_DIR}/src/libcadet/model/parts/ParticleDiffusionOperatorBase.cpp
 		${CMAKE_SOURCE_DIR}/src/libcadet/model/parts/ParticleDiffusionOperatorDG.cpp
+		${CMAKE_SOURCE_DIR}/src/libcadet/model/parts/ParticleDiffusionOperatorFV.cpp
 		${CMAKE_SOURCE_DIR}/src/libcadet/model/LumpedRateModelWithoutPoresDG.cpp
 		${CMAKE_SOURCE_DIR}/src/libcadet/model/ColumnModel1D.cpp
 	)

src/libcadet/model/particle/GeneralRateParticle.cpp

@@ -12,6 +12,7 @@
 
 #include "model/particle/GeneralRateParticle.hpp"
 #include "model/parts/ParticleDiffusionOperatorDG.hpp"
+#include "model/parts/ParticleDiffusionOperatorFV.hpp"
 
 #include "cadet/Exceptions.hpp"
 #include "BindingModelFactory.hpp"
@@ -115,16 +116,13 @@ namespace model
 
 		paramProvider.pushScope("discretization");
 
-		if (paramProvider.exists("SPATIAL_METHOD"))
-		{
-			const std::string parSpatialMethod = paramProvider.getString("SPATIAL_METHOD");
-			if (parSpatialMethod != "DG")
-				throw InvalidParameterException("Unsupported SPATIAL_METHOD '" + parSpatialMethod + "' for GeneralRateParticle. Only 'DG' is supported for now.");
-
+		const std::string parSpatialMethod = paramProvider.getString("SPATIAL_METHOD");
+		if (parSpatialMethod == "DG")
 			_parDiffOp = new parts::ParticleDiffusionOperatorDG();
-		}
+		else if (parSpatialMethod == "FV")
+			_parDiffOp = new parts::ParticleDiffusionOperatorFV();
 		else
-			_parDiffOp = new parts::ParticleDiffusionOperatorDG();
+			throw InvalidParameterException("Unsupported SPATIAL_METHOD '" + parSpatialMethod + "' for GeneralRateParticle. Only 'DG' and 'FV' are supported.");
 
 		paramProvider.popScope();
 
@@ -330,7 +328,7 @@ namespace model
 			// Add flux to column void / bulk volume
 			for (unsigned int comp = 0; comp < _nComp; ++comp)
 			{
-				// + 1/Beta_c * (surfaceToVolumeRatio_{p,j}) * d_j * (k_f * [c_l - c_p])
+				// + 1/Beta^c * (surfaceToVolumeRatio^p_j) * d_j * (k_f * [c^b - c^p])
 				resBulk[comp] += static_cast<ParamType>(filmDiff[comp]) * jacCF_val * static_cast<ParamType>(packing.parTypeVolFrac) * (yBulk[comp] - yPar[(nDiscPoints() - 1) * stridePoint() + comp]);
 			}
 		}

src/libcadet/model/parts/ParticleDiffusionOperatorDG.cpp

@@ -1610,7 +1610,7 @@ namespace parts
 					jacCl[0] += static_cast<double>(filmDiff[comp]) * (1.0 - colPorosity) / colPorosity
 					* _parGeomSurfToVol / static_cast<double>(_parRadius)
 					* static_cast<double>(parTypeVolFrac[_parTypeIdx + blk * nParType]);
-				// add Cl on Cp entries (added since these entries are also touched by bulk jacobian)
+				// add Cl on Cp entries
 				// row: already at bulk phase. already at current node and component.
 				// col: go to current particle phase entry.
 				jacCl[jacCp.row() - jacCl.row()] = -static_cast<double>(filmDiff[comp]) * (1.0 - colPorosity) / colPorosity
@@ -1622,7 +1622,7 @@ namespace parts
 				for (int node = _parPolyDeg; node >= 0; node--, jacCp -= strideParNode()) {
 					// row: already at particle. Already at current node and liquid state.
 					// col: original entry at outer node.
-					if (!outliersOnly) // Cp on Cb
+					if (!outliersOnly) // Cp on Cp
 						jacCp[entry - jacCp.row()]
 						+= static_cast<double>(filmDiff[comp]) * 2.0 / static_cast<double>(_deltaR[0]) * _parInvMM[_nParElem - 1](node, _nParNode - 1) * exIntLiftContribution / static_cast<double>(_parPorosity) / static_cast<double>(_poreAccessFactor[comp]);
 					// row: already at particle. Already at current node and liquid state.

src/libcadet/model/parts/ParticleDiffusionOperatorDG.hpp

@@ -200,7 +200,7 @@ namespace parts
 
 		/* DG specific operators */
 
-		active* _deltaR; //!< equidistant particle element spacing
+		active* _deltaR; //!< particle element spacing
 		Eigen::VectorXd _parNodes; //!< Array with positions of nodes in radial reference element for each particle
 		Eigen::MatrixXd _parPolyDerM; //!< Array with polynomial derivative Matrix for each particle
 		Eigen::MatrixXd* _minus_InvMM_ST; //!< equals minus inverse mass matrix times transposed stiffness matrix.