Outsource film diffusion Jacobian

Author: jbreue16

src/libcadet/model/GeneralRateModelDG.hpp

@@ -709,7 +709,7 @@ class GeneralRateModelDG : public UnitOperationBase
 
 		delete[] invBetaP;
 
-		calcFluxJacobians(secIdx);
+		_parDiffOp.calcFilmDiffJacobian(secIdx, idxr.offsetCp(), idxr.offsetC(), _disc.nPoints, static_cast<double>(_colPorosity), &_parTypeVolFrac[0], _globalJac);
 
 		return _globalJac.isCompressed(); // check if the jacobian estimation fits the pattern
 	}

src/libcadet/model/parts/ParticleDiffusionOperatorDG.cpp

@@ -2297,6 +2297,64 @@ namespace parts
 		return true;
 	}
 
+	int ParticleDiffusionOperatorDG::calcFilmDiffJacobian(unsigned int secIdx, const int offsetCp, const int offsetC, const int nBulkPoints, const double colPorosity, const active* const parTypeVolFrac, Eigen::SparseMatrix<double, RowMajor>& globalJac, bool outliersOnly)
+	{
+		for (unsigned int type = 0; type < _nParType; type++)
+		{
+			// lifting matrix entry for exact integration scheme depends on metrics for sphere and cylinder
+			double exIntLiftContribution = static_cast<double>(_Ir[_offsetMetric[type] + _nParElem[type] - 1][_nParNode[type] - 1]);
+			if (_parGeomSurfToVol[type] == _SurfVolRatioSlab)
+				exIntLiftContribution = 1.0;
+
+			// Ordering of diffusion:
+			// sec0type0comp0, sec0type0comp1, sec0type0comp2, sec0type1comp0, sec0type1comp1, sec0type1comp2,
+			// sec1type0comp0, sec1type0comp1, sec1type0comp2, sec1type1comp0, sec1type1comp1, sec1type1comp2, ...
+			active const* const filmDiff = getSectionDependentSlice(_filmDiffusion, _nComp * _nParType, secIdx) + type * _nComp;
+
+			linalg::BandedEigenSparseRowIterator jacCl(globalJac, offsetC);
+			linalg::BandedEigenSparseRowIterator jacCp(globalJac, offsetCp + (_nParPoints[type] + _nParPoints[type] - 1) * strideParNode(type));
+
+			for (unsigned int blk = 0; blk < nBulkPoints; blk++)
+			{
+				for (unsigned int comp = 0; comp < _nComp; comp++, ++jacCp, ++jacCl) {
+					// add Cl on Cl entries (added since these entries are also touched by bulk jacobian)
+					// row: already at bulk phase. already at current node and component.
+					// col: already at bulk phase. already at current node and component.
+					if (!outliersOnly)
+						jacCl[0] += static_cast<double>(filmDiff[comp]) * (1.0 - colPorosity) / colPorosity
+						* _parGeomSurfToVol[type] / static_cast<double>(_parRadius[type])
+						* static_cast<double>(parTypeVolFrac[type + blk * _nParType]);
+					// add Cl on Cp entries (added since these entries are also touched by bulk jacobian)
+					// 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
+						* _parGeomSurfToVol[type] / static_cast<double>(_parRadius[type])
+						* static_cast<double>(parTypeVolFrac[type + blk * _nParType]);
+
+
+					unsigned int entry = jacCp.row();
+					for (int node = _parPolyDeg[type]; node >= 0; node--, jacCp -= strideParNode(type)) {
+						// row: already at particle. Already at current node and liquid state.
+						// col: original entry at outer node.
+						if (!outliersOnly) // Cp on Cb
+							jacCp[entry - jacCp.row()]
+							+= static_cast<double>(filmDiff[comp]) * 2.0 / static_cast<double>(_deltaR[_offsetMetric[type]]) * _parInvMM[_offsetMetric[type] + _nParElem[type] - 1](node, _nParNode[type] - 1) * exIntLiftContribution / static_cast<double>(_parPorosity[type]) / static_cast<double>(_poreAccessFactor[type * _nComp + comp]);
+						// row: already at particle. Already at current node and liquid state.
+						// col: go to current bulk phase.
+						jacCp[jacCl.row() - jacCp.row()]
+							= -static_cast<double>(filmDiff[comp]) * 2.0 / static_cast<double>(_deltaR[_offsetMetric[type]]) * _parInvMM[_offsetMetric[type] + _nParElem[type] - 1](node, _nParNode[type] - 1) * exIntLiftContribution / static_cast<double>(_parPorosity[type]) / static_cast<double>(_poreAccessFactor[type * _nComp + comp]);
+					}
+					// set back iterator to first node as required by component loop
+					jacCp += _nParNode[type] * strideParNode(type);
+				}
+				if (blk < nBulkPoints - 1) // execute iteration statement only when condition is true in next loop.
+					jacCp += _strideBound[type] + (_nParPoints[type] - 1) * strideParNode(type);
+			}
+		}
+
+		return 1;
+	}
+
 }  // namespace parts
 
 }  // namespace model

src/libcadet/model/parts/ParticleDiffusionOperatorDG.hpp

@@ -223,6 +223,8 @@ namespace parts
 
 		void initializeDGjac(std::vector<double> parGeomSurfToVol);
 
+		int calcFilmDiffJacobian(unsigned int secIdx, const int offsetCp, const int offsetC, const int nBulkPoints, const double colPorosity, const active* const parTypeVolFrac, Eigen::SparseMatrix<double, RowMajor>& globalJac, bool outliersOnly = false);
+
 		int getParticleCoordinates(unsigned int parType, double* coords) const;
 
 		typedef Eigen::Triplet<double> T;