Skip to content

Commit 837b3b8

Browse files
jbreue16jbreue16
committed
fixup!
1 parent 6ae3000 commit 837b3b8

8 files changed

+65
-128
lines changed

‎src/libcadet/model/GeneralRateModelDG.cpp

Lines changed: 1 addition & 47 deletions
Original file line numberDiff line numberDiff line change
@@ -898,6 +898,7 @@ int GeneralRateModelDG::residualImpl(double t, unsigned int secIdx, StateType co
898898
y + idxr.offsetC() + colNode * idxr.strideColNode(),
899899
yDot ? yDot + idxr.offsetCp(ParticleTypeIndex{ parType }, ParticleIndex{ colNode }) : nullptr,
900900
res ? res + idxr.offsetCp(ParticleTypeIndex{ parType }, ParticleIndex{ colNode }) : nullptr,
901+
res ? res + idxr.offsetC() + colNode * idxr.strideColNode() : nullptr,
901902
packing, jacIt, tlmAlloc,
902903
typename cadet::ParamSens<ParamType>::enabled()
903904
);
@@ -908,8 +909,6 @@ int GeneralRateModelDG::residualImpl(double t, unsigned int secIdx, StateType co
908909

909910
BENCH_STOP(_timerResidualPar);
910911

911-
residualFlux<StateType, ResidualType, ParamType>(t, secIdx, y, yDot, res);
912-
913912
// Handle inlet DOFs, which are simply copied to the residual
914913
for (unsigned int i = 0; i < _disc.nComp; ++i)
915914
{
@@ -969,51 +968,6 @@ int GeneralRateModelDG::residualBulk(double t, unsigned int secIdx, StateType co
969968
return 0;
970969
}
971970

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

‎src/libcadet/model/LumpedRateModelWithPoresDG.cpp

Lines changed: 1 addition & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -977,6 +977,7 @@ int LumpedRateModelWithPoresDG::residualImpl(double t, unsigned int secIdx, Stat
977977
y + idxr.offsetC() + colNode * idxr.strideColNode(),
978978
yDot ? yDot + idxr.offsetCp(ParticleTypeIndex{ parType }, ParticleIndex{ colNode }) : nullptr,
979979
res ? res + idxr.offsetCp(ParticleTypeIndex{ parType }, ParticleIndex{ colNode }) : nullptr,
980+
res ? res + idxr.offsetC() + colNode * idxr.strideColNode() : nullptr,
980981
packing, jacIt, tlmAlloc,
981982
typename cadet::ParamSens<ParamType>::enabled()
982983
);
@@ -988,8 +989,6 @@ int LumpedRateModelWithPoresDG::residualImpl(double t, unsigned int secIdx, Stat
988989
if (!wantRes)
989990
return 0;
990991

991-
residualFlux<StateType, ResidualType, ParamType>(t, secIdx, y, yDot, res);
992-
993992
// Handle inlet DOFs, which are simply copied to res
994993
for (unsigned int i = 0; i < _disc.nComp; ++i)
995994
{
@@ -1044,40 +1043,6 @@ int LumpedRateModelWithPoresDG::residualBulk(double t, unsigned int secIdx, Stat
10441043
return 0;
10451044
}
10461045

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

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

Lines changed: 31 additions & 15 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, columnPackingParameters packing, 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, packing, 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, packing, 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, packing, 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, columnPackingParameters packing, 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, packing, 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, packing, 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, columnPackingParameters packing, 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, packing, 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, packing, 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, columnPackingParameters packing, 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, packing, 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, packing, 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, columnPackingParameters packing, 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());
@@ -315,9 +315,25 @@ namespace model
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+
// film diffusion bulk eq. term
324+
active const* const filmDiff = _parDiffOp->getFilmDiffusion(secIdx);
325+
const ParamType invBetaC = 1.0 / static_cast<ParamType>(packing.colPorosity) - 1.0;
326+
const ParamType jacCF_val = invBetaC * static_cast<ParamType>(surfaceToVolumeRatio());
327+
const ParamType jacPF_val = -1.0 / static_cast<ParamType>(getPorosity());
328+
329+
// Add flux to column void / bulk volume
330+
for (unsigned int comp = 0; comp < _nComp; ++comp)
331+
{
332+
// + 1/Beta_c * (surfaceToVolumeRatio_{p,j}) * d_j * (k_f * [c_l - c_p])
333+
resBulk[0] += static_cast<ParamType>(filmDiff[comp]) * jacCF_val * static_cast<ParamType>(packing.parTypeVolFrac[0]) * (yBulk[0] - yPar[(nDiscPoints() - 1) * stridePoint() + comp]);
334+
}
335+
336+
return 0;
321337
}
322338

323339
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, 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, 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, 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, columnPackingParameters packing, 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, columnPackingParameters packing, 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)