Skip to content

Commit b1dd002

Browse files
jbreue16jbreue16
committed
Outsource film diffusion Jacobian
1 parent 137d7fe commit b1dd002

File tree

4 files changed

+65
-71
lines changed

4 files changed

+65
-71
lines changed

‎src/libcadet/model/GeneralRateModelDG.cpp‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1050,7 +1050,10 @@ void GeneralRateModelDG::extractJacobianFromAD(active const* const adRes, unsign
10501050

10511051
/* Add analytically derived flux entries (only those that are part of the outlier bands) */
10521052
// todo extract these entries instead of analytical calculation?
1053-
calcFluxJacobians(_disc.curSection, true);
1053+
for (unsigned int parType = 0; parType < _disc.nParType; parType++)
1054+
{
1055+
_parDiffOp.calcFilmDiffJacobian(_disc.curSection, parType, idxr.offsetCp(ParticleTypeIndex{ static_cast<unsigned int>(parType) }), idxr.offsetC(), _disc.nPoints, static_cast<double>(_colPorosity), &_parTypeVolFrac[0], _globalJac, true);
1056+
}
10541057
}
10551058

10561059
#ifdef CADET_CHECK_ANALYTIC_JACOBIAN

‎src/libcadet/model/GeneralRateModelDG.hpp‎

Lines changed: 4 additions & 70 deletions
Original file line numberDiff line numberDiff line change
@@ -625,91 +625,25 @@ class GeneralRateModelDG : public UnitOperationBase
625625
globalJ.setFromTriplets(tripletList.begin(), tripletList.end());
626626
}
627627

628-
int calcFluxJacobians(unsigned int secIdx, bool outliersOnly = false)
629-
{
630-
Indexer idxr(_disc);
631-
632-
for (unsigned int type = 0; type < _disc.nParType; type++) {
633-
634-
// lifting matrix entry for exact integration scheme depends on metrics for sphere and cylinder
635-
double exIntLiftContribution = static_cast<double>(_parDiffOp._Ir[_parDiffOp._offsetMetric[type] + _disc.nParCell[type] - 1][_disc.nParNode[type] - 1]);
636-
if (_parGeomSurfToVol[type] == _disc.SurfVolRatioSlab)
637-
exIntLiftContribution = 1.0;
638-
639-
// Ordering of diffusion:
640-
// sec0type0comp0, sec0type0comp1, sec0type0comp2, sec0type1comp0, sec0type1comp1, sec0type1comp2,
641-
// sec1type0comp0, sec1type0comp1, sec1type0comp2, sec1type1comp0, sec1type1comp1, sec1type1comp2, ...
642-
active const* const filmDiff = getSectionDependentSlice(_filmDiffusion, _disc.nComp * _disc.nParType, secIdx) + type * _disc.nComp;
643-
644-
linalg::BandedEigenSparseRowIterator jacCl(_globalJac, idxr.offsetC());
645-
linalg::BandedEigenSparseRowIterator jacCp(_globalJac, idxr.offsetCp(ParticleTypeIndex{ type }) + (_disc.nParPoints[type] - 1) * idxr.strideParNode(type));
646-
647-
for (unsigned int colNode = 0; colNode < _disc.nPoints; colNode++)
648-
{
649-
for (unsigned int comp = 0; comp < _disc.nComp; comp++, ++jacCp, ++jacCl) {
650-
// add Cl on Cl entries (added since these entries are also touched by bulk jacobian)
651-
// row: already at bulk phase. already at current node and component.
652-
// col: already at bulk phase. already at current node and component.
653-
if (!outliersOnly)
654-
jacCl[0] += static_cast<double>(filmDiff[comp]) * (1.0 - static_cast<double>(_colPorosity)) / static_cast<double>(_colPorosity)
655-
* _parGeomSurfToVol[type] / static_cast<double>(_parRadius[type])
656-
* _parTypeVolFrac[type + colNode * _disc.nParType].getValue();
657-
// add Cl on Cp entries (added since these entries are also touched by bulk jacobian)
658-
// row: already at bulk phase. already at current node and component.
659-
// col: go to current particle phase entry.
660-
jacCl[jacCp.row() - jacCl.row()] = -static_cast<double>(filmDiff[comp]) * (1.0 - static_cast<double>(_colPorosity)) / static_cast<double>(_colPorosity)
661-
* _parGeomSurfToVol[type] / static_cast<double>(_parRadius[type])
662-
* _parTypeVolFrac[type + colNode * _disc.nParType].getValue();
663-
664-
665-
unsigned int entry = jacCp.row();
666-
for (int node = _disc.parPolyDeg[type]; node >= 0; node--, jacCp -= idxr.strideParNode(type)) {
667-
// row: already at particle. Already at current node and liquid state.
668-
// col: original entry at outer node.
669-
if (!outliersOnly) // Cp on Cb
670-
jacCp[entry - jacCp.row()]
671-
+= static_cast<double>(filmDiff[comp]) * 2.0 / static_cast<double>(_parDiffOp._deltaR[_parDiffOp._offsetMetric[type]]) * _parDiffOp._parInvMM[_parDiffOp._offsetMetric[type] + _disc.nParCell[type] - 1](node, _disc.nParNode[type] - 1) * exIntLiftContribution / static_cast<double>(_parPorosity[type]) / static_cast<double>(_poreAccessFactor[type * _disc.nComp + comp]);
672-
// row: already at particle. Already at current node and liquid state.
673-
// col: go to current bulk phase.
674-
jacCp[jacCl.row() - jacCp.row()]
675-
= -static_cast<double>(filmDiff[comp]) * 2.0 / static_cast<double>(_parDiffOp._deltaR[_parDiffOp._offsetMetric[type]]) * _parDiffOp._parInvMM[_parDiffOp._offsetMetric[type] + _disc.nParCell[type] - 1](node, _disc.nParNode[type] - 1) * exIntLiftContribution / static_cast<double>(_parPorosity[type]) / static_cast<double>(_poreAccessFactor[type * _disc.nComp + comp]);
676-
}
677-
// set back iterator to first node as required by component loop
678-
jacCp += _disc.nParNode[type] * idxr.strideParNode(type);
679-
}
680-
if (colNode < _disc.nPoints - 1) // execute iteration statement only when condition is true in next loop.
681-
jacCp += _disc.strideBound[type] + (_disc.nParPoints[type] - 1) * idxr.strideParNode(type);
682-
}
683-
}
684-
685-
return 1;
686-
}
687-
688628
int calcStaticAnaJacobian_GRM(unsigned int secIdx)
689629
{
690630
Indexer idxr(_disc);
691631
// inlet and bulk jacobian
692632
_convDispOp.calcStaticAnaJacobian(_globalJac, _jacInlet, idxr.offsetC());
693633

694-
double* invBetaP = new double[_disc.nComp];
695-
696634
// particle jacobian (without isotherm, which is handled in residualKernel)
697635
for (int colNode = 0; colNode < _disc.nPoints; colNode++)
698636
{
699637
for (int type = 0; type < _disc.nParType; type++)
700638
{
701-
for (int comp = 0; comp < _disc.nComp; comp++)
702-
{
703-
invBetaP[comp] = (1.0 - static_cast<double>(_parPorosity[type])) / (static_cast<double>(_poreAccessFactor[_disc.nComp * type + comp]) * static_cast<double>(_parPorosity[type]));
704-
}
705-
706639
_parDiffOp.calcStaticAnaParticleDiffJacobian(secIdx, type, colNode, _binding[type]->reactionQuasiStationarity(), idxr.offsetCp(ParticleTypeIndex{ static_cast<unsigned int>(type) }, ParticleIndex{ static_cast<unsigned int>(colNode) }), _globalJac);
707640
}
708641
}
709642

710-
delete[] invBetaP;
711-
712-
calcFluxJacobians(secIdx);
643+
for (unsigned int parType = 0; parType < _disc.nParType; parType++)
644+
{
645+
_parDiffOp.calcFilmDiffJacobian(secIdx, parType, idxr.offsetCp(ParticleTypeIndex{ static_cast<unsigned int>(parType) }), idxr.offsetC(), _disc.nPoints, static_cast<double>(_colPorosity), &_parTypeVolFrac[0], _globalJac);
646+
}
713647

714648
return _globalJac.isCompressed(); // check if the jacobian estimation fits the pattern
715649
}

‎src/libcadet/model/parts/ParticleDiffusionOperatorDG.cpp‎

Lines changed: 55 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2297,6 +2297,61 @@ namespace parts
22972297
return true;
22982298
}
22992299

2300+
int ParticleDiffusionOperatorDG::calcFilmDiffJacobian(unsigned int secIdx, const int parType, const int offsetCp, const int offsetC, const int nBulkPoints, const double colPorosity, const active* const parTypeVolFrac, Eigen::SparseMatrix<double, RowMajor>& globalJac, bool outliersOnly)
2301+
{
2302+
// lifting matrix entry for exact integration scheme depends on metrics for sphere and cylinder
2303+
double exIntLiftContribution = static_cast<double>(_Ir[_offsetMetric[parType] + _nParElem[parType] - 1][_nParNode[parType] - 1]);
2304+
if (_parGeomSurfToVol[parType] == _SurfVolRatioSlab)
2305+
exIntLiftContribution = 1.0;
2306+
2307+
// Ordering of diffusion:
2308+
// sec0type0comp0, sec0type0comp1, sec0type0comp2, sec0type1comp0, sec0type1comp1, sec0type1comp2,
2309+
// sec1type0comp0, sec1type0comp1, sec1type0comp2, sec1type1comp0, sec1type1comp1, sec1type1comp2, ...
2310+
active const* const filmDiff = getSectionDependentSlice(_filmDiffusion, _nComp * _nParType, secIdx) + parType * _nComp;
2311+
2312+
linalg::BandedEigenSparseRowIterator jacCl(globalJac, offsetC);
2313+
linalg::BandedEigenSparseRowIterator jacCp(globalJac, offsetCp + (_nParPoints[parType] - 1) * strideParNode(parType)); // iterator at the outer particle boundary
2314+
2315+
for (unsigned int blk = 0; blk < nBulkPoints; blk++)
2316+
{
2317+
for (unsigned int comp = 0; comp < _nComp; comp++, ++jacCp, ++jacCl) {
2318+
// add Cl on Cl entries (added since these entries are also touched by bulk jacobian)
2319+
// row: already at bulk phase. already at current node and component.
2320+
// col: already at bulk phase. already at current node and component.
2321+
if (!outliersOnly)
2322+
jacCl[0] += static_cast<double>(filmDiff[comp]) * (1.0 - colPorosity) / colPorosity
2323+
* _parGeomSurfToVol[parType] / static_cast<double>(_parRadius[parType])
2324+
* static_cast<double>(parTypeVolFrac[parType + blk * _nParType]);
2325+
// add Cl on Cp entries (added since these entries are also touched by bulk jacobian)
2326+
// row: already at bulk phase. already at current node and component.
2327+
// col: go to current particle phase entry.
2328+
jacCl[jacCp.row() - jacCl.row()] = -static_cast<double>(filmDiff[comp]) * (1.0 - colPorosity) / colPorosity
2329+
* _parGeomSurfToVol[parType] / static_cast<double>(_parRadius[parType])
2330+
* static_cast<double>(parTypeVolFrac[parType + blk * _nParType]);
2331+
2332+
2333+
unsigned int entry = jacCp.row();
2334+
for (int node = _parPolyDeg[parType]; node >= 0; node--, jacCp -= strideParNode(parType)) {
2335+
// row: already at particle. Already at current node and liquid state.
2336+
// col: original entry at outer node.
2337+
if (!outliersOnly) // Cp on Cb
2338+
jacCp[entry - jacCp.row()]
2339+
+= static_cast<double>(filmDiff[comp]) * 2.0 / static_cast<double>(_deltaR[_offsetMetric[parType]]) * _parInvMM[_offsetMetric[parType] + _nParElem[parType] - 1](node, _nParNode[parType] - 1) * exIntLiftContribution / static_cast<double>(_parPorosity[parType]) / static_cast<double>(_poreAccessFactor[parType * _nComp + comp]);
2340+
// row: already at particle. Already at current node and liquid state.
2341+
// col: go to current bulk phase.
2342+
jacCp[jacCl.row() - jacCp.row()]
2343+
= -static_cast<double>(filmDiff[comp]) * 2.0 / static_cast<double>(_deltaR[_offsetMetric[parType]]) * _parInvMM[_offsetMetric[parType] + _nParElem[parType] - 1](node, _nParNode[parType] - 1) * exIntLiftContribution / static_cast<double>(_parPorosity[parType]) / static_cast<double>(_poreAccessFactor[parType * _nComp + comp]);
2344+
}
2345+
// set back iterator to first node as required by component loop
2346+
jacCp += _nParNode[parType] * strideParNode(parType);
2347+
}
2348+
if (blk < nBulkPoints - 1) // execute iteration statement only when condition is true in next loop.
2349+
jacCp += _strideBound[parType] + (_nParPoints[parType] - 1) * strideParNode(parType);
2350+
}
2351+
2352+
return 1;
2353+
}
2354+
23002355
} // namespace parts
23012356

23022357
} // namespace model

‎src/libcadet/model/parts/ParticleDiffusionOperatorDG.hpp‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -223,6 +223,8 @@ namespace parts
223223

224224
void initializeDGjac(std::vector<double> parGeomSurfToVol);
225225

226+
int calcFilmDiffJacobian(unsigned int secIdx, const int parType, const int offsetCp, const int offsetC, const int nBulkPoints, const double colPorosity, const active* const parTypeVolFrac, Eigen::SparseMatrix<double, RowMajor>& globalJac, bool outliersOnly = false);
227+
226228
int getParticleCoordinates(unsigned int parType, double* coords) const;
227229

228230
typedef Eigen::Triplet<double> T;

0 commit comments

Comments
 (0)