Skip to content

Commit 1410e29

Browse files
jbreue16jbreue16
committed
Move computation of film diffusion bulk term to particle model (#435)
1 parent 1315adf commit 1410e29

8 files changed

+103
-137
lines changed

‎src/libcadet/model/GeneralRateModelDG.cpp

Lines changed: 8 additions & 49 deletions
Original file line numberDiff line numberDiff line change
@@ -886,14 +886,20 @@ int GeneralRateModelDG::residualImpl(double t, unsigned int secIdx, StateType co
886886
const unsigned int colNode = pblk % _disc.nPoints;
887887

888888
linalg::BandedEigenSparseRowIterator jacIt(_globalJac, idxr.offsetCp(ParticleTypeIndex{ parType }, ParticleIndex{ colNode }));
889-
ColumnPosition colPos{ _convDispOp.relativeCoordinate(colNode), 0.0, 0.0 }; // Relative position of current node - needed in externally dependent adsorption kinetic
889+
model::columnPackingParameters packing
890+
{
891+
&_parTypeVolFrac[parType] + _disc.nParType * colNode,
892+
_colPorosity,
893+
ColumnPosition{ _convDispOp.relativeCoordinate(colNode), 0.0, 0.0 }
894+
};
890895

891896
_particle[parType].residual(t, secIdx,
892897
y + idxr.offsetCp(ParticleTypeIndex{ parType }, ParticleIndex{ colNode }),
893898
y + idxr.offsetC() + colNode * idxr.strideColNode(),
894899
yDot ? yDot + idxr.offsetCp(ParticleTypeIndex{ parType }, ParticleIndex{ colNode }) : nullptr,
895900
res ? res + idxr.offsetCp(ParticleTypeIndex{ parType }, ParticleIndex{ colNode }) : nullptr,
896-
colPos, jacIt, tlmAlloc,
901+
res ? res + idxr.offsetC() + colNode * idxr.strideColNode() : nullptr,
902+
packing, jacIt, tlmAlloc,
897903
typename cadet::ParamSens<ParamType>::enabled()
898904
);
899905
}
@@ -903,8 +909,6 @@ int GeneralRateModelDG::residualImpl(double t, unsigned int secIdx, StateType co
903909

904910
BENCH_STOP(_timerResidualPar);
905911

906-
residualFlux<StateType, ResidualType, ParamType>(t, secIdx, y, yDot, res);
907-
908912
// Handle inlet DOFs, which are simply copied to the residual
909913
for (unsigned int i = 0; i < _disc.nComp; ++i)
910914
{
@@ -964,51 +968,6 @@ int GeneralRateModelDG::residualBulk(double t, unsigned int secIdx, StateType co
964968
return 0;
965969
}
966970

967-
template <typename StateType, typename ResidualType, typename ParamType>
968-
int GeneralRateModelDG::residualFlux(double t, unsigned int secIdx, StateType const* yBase, double const* yDotBase, ResidualType* resBase)
969-
{
970-
Indexer idxr(_disc);
971-
972-
const ParamType invBetaC = 1.0 / static_cast<ParamType>(_colPorosity) - 1.0;
973-
974-
// Get offsets
975-
ResidualType* const resCol = resBase + idxr.offsetC();
976-
StateType const* const yCol = yBase + idxr.offsetC();
977-
978-
for (unsigned int type = 0; type < _disc.nParType; ++type)
979-
{
980-
ResidualType* const resParType = resBase + idxr.offsetCp(ParticleTypeIndex{type});
981-
StateType const* const yParType = yBase + idxr.offsetCp(ParticleTypeIndex{type});
982-
983-
const ParamType epsP = static_cast<ParamType>(_particle[type].getPorosity());
984-
985-
// Ordering of diffusion:
986-
// sec0type0comp0, sec0type0comp1, sec0type0comp2, sec0type1comp0, sec0type1comp1, sec0type1comp2,
987-
// sec1type0comp0, sec1type0comp1, sec1type0comp2, sec1type1comp0, sec1type1comp1, sec1type1comp2, ...
988-
active const* const filmDiff = getSectionDependentSlice(_filmDiffusion, _disc.nComp * _disc.nParType, secIdx) + type * _disc.nComp;
989-
990-
const ParamType surfaceToVolumeRatio = static_cast<ParamType>(_particle[type].surfaceToVolumeRatio());
991-
992-
const ParamType jacCF_val = invBetaC * surfaceToVolumeRatio;
993-
const ParamType jacPF_val = -1.0 / epsP;
994-
995-
// Add flux to column void / bulk volume
996-
for (unsigned int i = 0; i < _disc.nPoints * _disc.nComp; ++i)
997-
{
998-
const unsigned int colNode = i / _disc.nComp;
999-
const unsigned int comp = i - colNode * _disc.nComp;
1000-
// + 1/Beta_c * (surfaceToVolumeRatio_{p,j}) * d_j * (k_f * [c_l - c_p])
1001-
resCol[i] += static_cast<ParamType>(filmDiff[comp]) * jacCF_val * static_cast<ParamType>(_parTypeVolFrac[type + colNode * _disc.nParType])
1002-
* (yCol[i] - yParType[colNode * idxr.strideParBlock(type) + (_disc.nParPoints[type] - 1) * idxr.strideParNode(type) + comp]);
1003-
}
1004-
1005-
// Bead boundary condition is computed in particle residual.
1006-
1007-
}
1008-
1009-
return 0;
1010-
}
1011-
1012971
parts::cell::CellParameters GeneralRateModelDG::makeCellResidualParams(unsigned int parType, int const* qsReaction) const
1013972
{
1014973
return parts::cell::CellParameters

‎src/libcadet/model/LumpedRateModelWithPoresDG.cpp

Lines changed: 10 additions & 38 deletions
Original file line numberDiff line numberDiff line change
@@ -960,17 +960,25 @@ int LumpedRateModelWithPoresDG::residualImpl(double t, unsigned int secIdx, Stat
960960

961961
for (unsigned int parType = 0; parType < _disc.nParType; parType++)
962962
{
963+
963964
for (unsigned int colNode = 0; colNode < _disc.nPoints; colNode++)
964965
{
965-
const double z = _convDispOp.relativeCoordinate(colNode);
966+
model::columnPackingParameters packing
967+
{
968+
&_parTypeVolFrac[parType] + _disc.nParType * colNode,
969+
_colPorosity,
970+
ColumnPosition{ _convDispOp.relativeCoordinate(colNode), 0.0, 0.0 }
971+
};
972+
966973
linalg::BandedEigenSparseRowIterator jacIt(_globalJac, idxr.offsetCp(ParticleTypeIndex{ parType }, ParticleIndex{ colNode }) - idxr.offsetC());
967974

968975
_particle[parType].residual(t, secIdx,
969976
y + idxr.offsetCp(ParticleTypeIndex{ parType }, ParticleIndex{ colNode }),
970977
y + idxr.offsetC() + colNode * idxr.strideColNode(),
971978
yDot ? yDot + idxr.offsetCp(ParticleTypeIndex{ parType }, ParticleIndex{ colNode }) : nullptr,
972979
res ? res + idxr.offsetCp(ParticleTypeIndex{ parType }, ParticleIndex{ colNode }) : nullptr,
973-
ColumnPosition{ z, 0.0, 0.0 }, jacIt, tlmAlloc,
980+
res ? res + idxr.offsetC() + colNode * idxr.strideColNode() : nullptr,
981+
packing, jacIt, tlmAlloc,
974982
typename cadet::ParamSens<ParamType>::enabled()
975983
);
976984
}
@@ -981,8 +989,6 @@ int LumpedRateModelWithPoresDG::residualImpl(double t, unsigned int secIdx, Stat
981989
if (!wantRes)
982990
return 0;
983991

984-
residualFlux<StateType, ResidualType, ParamType>(t, secIdx, y, yDot, res);
985-
986992
// Handle inlet DOFs, which are simply copied to res
987993
for (unsigned int i = 0; i < _disc.nComp; ++i)
988994
{
@@ -1037,40 +1043,6 @@ int LumpedRateModelWithPoresDG::residualBulk(double t, unsigned int secIdx, Stat
10371043
return 0;
10381044
}
10391045

1040-
template <typename StateType, typename ResidualType, typename ParamType>
1041-
int LumpedRateModelWithPoresDG::residualFlux(double t, unsigned int secIdx, StateType const* yBase, double const* yDotBase, ResidualType* resBase)
1042-
{
1043-
Indexer idxr(_disc);
1044-
1045-
const ParamType invBetaC = 1.0 / static_cast<ParamType>(_colPorosity) - 1.0;
1046-
1047-
// Get offsets
1048-
ResidualType* const resCol = resBase + idxr.offsetC();
1049-
StateType const* const yCol = yBase + idxr.offsetC();
1050-
1051-
for (unsigned int type = 0; type < _disc.nParType; ++type)
1052-
{
1053-
ResidualType* const resParType = resBase + idxr.offsetCp(ParticleTypeIndex{ type });
1054-
StateType const* const yParType = yBase + idxr.offsetCp(ParticleTypeIndex{ type });
1055-
1056-
const ParamType surfToVolRatio = static_cast<ParamType>(_particle[type].surfaceToVolumeRatio());
1057-
active const* const filmDiff = getSectionDependentSlice(_filmDiffusion, _disc.nComp * _disc.nParType, secIdx) + type * _disc.nComp;
1058-
active const* const poreAccFactor = _poreAccessFactor.data() + type * _disc.nComp;
1059-
1060-
const ParamType jacCF_val = invBetaC * surfToVolRatio;
1061-
1062-
// Add flux to column void / bulk volume equations
1063-
for (unsigned int i = 0; i < _disc.nPoints * _disc.nComp; ++i)
1064-
{
1065-
const unsigned int colNode = i / _disc.nComp;
1066-
const unsigned int comp = i % _disc.nComp;
1067-
resCol[i] += jacCF_val * static_cast<ParamType>(filmDiff[comp]) * static_cast<ParamType>(_parTypeVolFrac[type + _disc.nParType * colNode]) * (yCol[i] - yParType[colNode * idxr.strideParBlock(type) + comp]);
1068-
}
1069-
}
1070-
1071-
return 0;
1072-
}
1073-
10741046
void LumpedRateModelWithPoresDG::assembleFluxJacobian(double t, unsigned int secIdx)
10751047
{
10761048
calcFluxJacobians(secIdx, false);

‎src/libcadet/model/particle/GeneralRateParticle.cpp

Lines changed: 37 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -246,44 +246,44 @@ namespace model
246246
};
247247
}
248248

249-
int GeneralRateParticle::residual(double t, unsigned int secIdx, double const* yPar, double const* yBulk, double const* yDotPar, double* resPar, ColumnPosition colPos, linalg::BandedEigenSparseRowIterator& jacIt, LinearBufferAllocator tlmAlloc, WithoutParamSensitivity)
249+
int GeneralRateParticle::residual(double t, unsigned int secIdx, double const* yPar, double const* yBulk, double const* yDotPar, double* resPar, double* resBulk, columnPackingParameters packing, linalg::BandedEigenSparseRowIterator& jacIt, LinearBufferAllocator tlmAlloc, WithoutParamSensitivity)
250250
{
251251
if (resPar)
252252
{
253253
if (jacIt.data())
254-
return residualImpl<double, double, double, true, true>(t, secIdx, yPar, yBulk, yDotPar, resPar, colPos, jacIt, tlmAlloc);
254+
return residualImpl<double, double, double, true, true>(t, secIdx, yPar, yBulk, yDotPar, resPar, resBulk, packing, jacIt, tlmAlloc);
255255
else
256-
return residualImpl<double, double, double, false, true>(t, secIdx, yPar, yBulk, yDotPar, resPar, colPos, jacIt, tlmAlloc);
256+
return residualImpl<double, double, double, false, true>(t, secIdx, yPar, yBulk, yDotPar, resPar, resBulk, packing, jacIt, tlmAlloc);
257257
}
258258
else if (jacIt.data())
259-
return residualImpl<double, double, double, true, false>(t, secIdx, yPar, yBulk, yDotPar, resPar, colPos, jacIt, tlmAlloc);
259+
return residualImpl<double, double, double, true, false>(t, secIdx, yPar, yBulk, yDotPar, resPar, resBulk, packing, jacIt, tlmAlloc);
260260
else
261261
return -1;
262262
}
263-
int GeneralRateParticle::residual(double t, unsigned int secIdx, double const* yPar, double const* yBulk, double const* yDotPar, active* resPar, ColumnPosition colPos, linalg::BandedEigenSparseRowIterator& jacIt, LinearBufferAllocator tlmAlloc, WithParamSensitivity)
263+
int GeneralRateParticle::residual(double t, unsigned int secIdx, double const* yPar, double const* yBulk, double const* yDotPar, active* resPar, active* resBulk, columnPackingParameters packing, linalg::BandedEigenSparseRowIterator& jacIt, LinearBufferAllocator tlmAlloc, WithParamSensitivity)
264264
{
265265
if (jacIt.data())
266-
return residualImpl<double, active, active, true, true>(t, secIdx, yPar, yBulk, yDotPar, resPar, colPos, jacIt, tlmAlloc);
266+
return residualImpl<double, active, active, true, true>(t, secIdx, yPar, yBulk, yDotPar, resPar, resBulk, packing, jacIt, tlmAlloc);
267267
else
268-
return residualImpl<double, active, active, false, true>(t, secIdx, yPar, yBulk, yDotPar, resPar, colPos, jacIt, tlmAlloc);
268+
return residualImpl<double, active, active, false, true>(t, secIdx, yPar, yBulk, yDotPar, resPar, resBulk, packing, jacIt, tlmAlloc);
269269
}
270-
int GeneralRateParticle::residual(double t, unsigned int secIdx, active const* yPar, active const* yBulk, double const* yDotPar, active* resPar, ColumnPosition colPos, linalg::BandedEigenSparseRowIterator& jacIt, LinearBufferAllocator tlmAlloc, WithoutParamSensitivity)
270+
int GeneralRateParticle::residual(double t, unsigned int secIdx, active const* yPar, active const* yBulk, double const* yDotPar, active* resPar, active* resBulk, columnPackingParameters packing, linalg::BandedEigenSparseRowIterator& jacIt, LinearBufferAllocator tlmAlloc, WithoutParamSensitivity)
271271
{
272272
if (jacIt.data())
273-
return residualImpl<active, active, double, true, true>(t, secIdx, yPar, yBulk, yDotPar, resPar, colPos, jacIt, tlmAlloc);
273+
return residualImpl<active, active, double, true, true>(t, secIdx, yPar, yBulk, yDotPar, resPar, resBulk, packing, jacIt, tlmAlloc);
274274
else
275-
return residualImpl<active, active, double, false, true>(t, secIdx, yPar, yBulk, yDotPar, resPar, colPos, jacIt, tlmAlloc);
275+
return residualImpl<active, active, double, false, true>(t, secIdx, yPar, yBulk, yDotPar, resPar, resBulk, packing, jacIt, tlmAlloc);
276276
}
277-
int GeneralRateParticle::residual(double t, unsigned int secIdx, active const* yPar, active const* yBulk, double const* yDotPar, active* resPar, ColumnPosition colPos, linalg::BandedEigenSparseRowIterator& jacIt, LinearBufferAllocator tlmAlloc, WithParamSensitivity)
277+
int GeneralRateParticle::residual(double t, unsigned int secIdx, active const* yPar, active const* yBulk, double const* yDotPar, active* resPar, active* resBulk, columnPackingParameters packing, linalg::BandedEigenSparseRowIterator& jacIt, LinearBufferAllocator tlmAlloc, WithParamSensitivity)
278278
{
279279
if (jacIt.data())
280-
return residualImpl<active, active, active, true, true>(t, secIdx, yPar, yBulk, yDotPar, resPar, colPos, jacIt, tlmAlloc);
280+
return residualImpl<active, active, active, true, true>(t, secIdx, yPar, yBulk, yDotPar, resPar, resBulk, packing, jacIt, tlmAlloc);
281281
else
282-
return residualImpl<active, active, active, false, true>(t, secIdx, yPar, yBulk, yDotPar, resPar, colPos, jacIt, tlmAlloc);
282+
return residualImpl<active, active, active, false, true>(t, secIdx, yPar, yBulk, yDotPar, resPar, resBulk, packing, jacIt, tlmAlloc);
283283
}
284284

285285
template <typename StateType, typename ResidualType, typename ParamType, bool wantNonLinJac, bool wantRes>
286-
int GeneralRateParticle::residualImpl(double t, unsigned int secIdx, StateType const* yPar, StateType const* yBulk, double const* yDotPar, ResidualType* resPar, ColumnPosition colPos, linalg::BandedEigenSparseRowIterator& jacIt, LinearBufferAllocator tlmAlloc)
286+
int GeneralRateParticle::residualImpl(double t, unsigned int secIdx, StateType const* yPar, StateType const* yBulk, double const* yDotPar, ResidualType* resPar, ResidualType* resBulk, columnPackingParameters packing, linalg::BandedEigenSparseRowIterator& jacIt, LinearBufferAllocator tlmAlloc)
287287
{
288288
int const* const qsBinding = _binding->reactionQuasiStationarity();
289289
const parts::cell::CellParameters cellResParams = makeCellResidualParams(qsBinding, _parDiffOp->nBound());
@@ -300,24 +300,43 @@ namespace model
300300
ResidualType* local_res = resPar ? resPar + par * stridePoint() : nullptr;
301301

302302
// r (particle) coordinate of current node (particle radius normed to 1) - needed in externally dependent adsorption kinetic
303-
colPos.particle = relativeCoordinate(par);
303+
packing.colPos.particle = relativeCoordinate(par);
304304

305305
if (wantRes)
306306
parts::cell::residualKernel<StateType, ResidualType, ParamType, parts::cell::CellParameters, linalg::BandedEigenSparseRowIterator, wantNonLinJac, true>(
307-
t, secIdx, colPos, local_y, local_yDot, local_res, jacIt, cellResParams, tlmAlloc
307+
t, secIdx, packing.colPos, local_y, local_yDot, local_res, jacIt, cellResParams, tlmAlloc
308308
);
309309
else
310310
parts::cell::residualKernel<StateType, ResidualType, ParamType, parts::cell::CellParameters, linalg::BandedEigenSparseRowIterator, wantNonLinJac, false, false>(
311-
t, secIdx, colPos, local_y, local_yDot, local_res, jacIt, cellResParams, tlmAlloc
311+
t, secIdx, packing.colPos, local_y, local_yDot, local_res, jacIt, cellResParams, tlmAlloc
312312
);
313313

314314
// Move rowiterator to next particle node
315315
jacIt += stridePoint();
316316
}
317317

318+
// particle diffusion, including film diffusion boundary condition
318319
ResidualType* wantResPtr = wantRes ? resPar : nullptr;
319320
linalg::BandedEigenSparseRowIterator jacSafe = wantNonLinJac ? jacBase : linalg::BandedEigenSparseRowIterator{};
320-
return _parDiffOp->residual(t, secIdx, yPar, yBulk, yDotPar, wantResPtr, jacSafe, typename ParamSens<ParamType>::enabled());
321+
_parDiffOp->residual(t, secIdx, yPar, yBulk, yDotPar, wantResPtr, jacSafe, typename ParamSens<ParamType>::enabled());
322+
323+
if (wantRes)
324+
{
325+
// film diffusion bulk eq. term
326+
active const* const filmDiff = _parDiffOp->getFilmDiffusion(secIdx);
327+
const ParamType invBetaC = 1.0 / static_cast<ParamType>(packing.colPorosity) - 1.0;
328+
const ParamType jacCF_val = invBetaC * static_cast<ParamType>(surfaceToVolumeRatio());
329+
const ParamType jacPF_val = -1.0 / static_cast<ParamType>(getPorosity());
330+
331+
// Add flux to column void / bulk volume
332+
for (unsigned int comp = 0; comp < _nComp; ++comp)
333+
{
334+
// + 1/Beta_c * (surfaceToVolumeRatio_{p,j}) * d_j * (k_f * [c_l - c_p])
335+
resBulk[0] += static_cast<ParamType>(filmDiff[comp]) * jacCF_val * static_cast<ParamType>(packing.parTypeVolFrac[0]) * (yBulk[0] - yPar[(nDiscPoints() - 1) * stridePoint() + comp]);
336+
}
337+
}
338+
339+
return 0;
321340
}
322341

323342
unsigned int GeneralRateParticle::jacobianNNZperParticle() const

‎src/libcadet/model/particle/GeneralRateParticle.hpp

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -94,10 +94,10 @@ namespace parts
9494

9595
bool notifyDiscontinuousSectionTransition(double t, unsigned int secIdx, active const* const filmDiff, active const* const poreAccessFactor) override;
9696

97-
int residual(double t, unsigned int secIdx, double const* yPar, double const* yBulk, double const* yDotPar, double* resPar, ColumnPosition colPos, linalg::BandedEigenSparseRowIterator& jacIt, LinearBufferAllocator tlmAlloc, WithoutParamSensitivity) override;
98-
int residual(double t, unsigned int secIdx, double const* yPar, double const* yBulk, double const* yDotPar, active* resPar, ColumnPosition colPos, linalg::BandedEigenSparseRowIterator& jacIt, LinearBufferAllocator tlmAlloc, WithParamSensitivity) override;
99-
int residual(double t, unsigned int secIdx, active const* yPar, active const* yBulk, double const* yDotPar, active* resPar, ColumnPosition colPos, linalg::BandedEigenSparseRowIterator& jacIt, LinearBufferAllocator tlmAlloc, WithoutParamSensitivity) override;
100-
int residual(double t, unsigned int secIdx, active const* yPar, active const* yBulk, double const* yDotPar, active* resPar, ColumnPosition colPos, linalg::BandedEigenSparseRowIterator& jacIt, LinearBufferAllocator tlmAlloc, WithParamSensitivity) override;
97+
int residual(double t, unsigned int secIdx, double const* yPar, double const* yBulk, double const* yDotPar, double* resPar, double* resBulk, columnPackingParameters packing, linalg::BandedEigenSparseRowIterator& jacIt, LinearBufferAllocator tlmAlloc, WithoutParamSensitivity) override;
98+
int residual(double t, unsigned int secIdx, double const* yPar, double const* yBulk, double const* yDotPar, active* resPar, active* resBulk, columnPackingParameters packing, linalg::BandedEigenSparseRowIterator& jacIt, LinearBufferAllocator tlmAlloc, WithParamSensitivity) override;
99+
int residual(double t, unsigned int secIdx, active const* yPar, active const* yBulk, double const* yDotPar, active* resPar, active* resBulk, columnPackingParameters packing, linalg::BandedEigenSparseRowIterator& jacIt, LinearBufferAllocator tlmAlloc, WithoutParamSensitivity) override;
100+
int residual(double t, unsigned int secIdx, active const* yPar, active const* yBulk, double const* yDotPar, active* resPar, active* resBulk, columnPackingParameters packing, linalg::BandedEigenSparseRowIterator& jacIt, LinearBufferAllocator tlmAlloc, WithParamSensitivity) override;
101101

102102
double relativeCoordinate(const unsigned int nodeIdx) const CADET_NOEXCEPT override { return _parDiffOp->relativeCoordinate(nodeIdx); }
103103

@@ -165,7 +165,7 @@ namespace parts
165165
inline int stridePoint() const CADET_NOEXCEPT { return static_cast<int>(_nComp) + _parDiffOp->strideBound(); }
166166

167167
template <typename StateType, typename ResidualType, typename ParamType, bool wantJac, bool wantRes>
168-
int residualImpl(double t, unsigned int secIdx, StateType const* yPar, StateType const* yBulk, double const* yDotPar, ResidualType* resPar, ColumnPosition colPos, linalg::BandedEigenSparseRowIterator& jacIt, LinearBufferAllocator tlmAlloc);
168+
int residualImpl(double t, unsigned int secIdx, StateType const* yPar, StateType const* yBulk, double const* yDotPar, ResidualType* resPar, ResidualType* resBulk, columnPackingParameters packing, linalg::BandedEigenSparseRowIterator& jacIt, LinearBufferAllocator tlmAlloc);
169169
};
170170

171171
} // namespace model

0 commit comments

Comments
 (0)