Skip to content

Commit becf645

Browse files
jbreue16jbreue16
authored
Move film diffusion to particle model (#435)
1 parent 1315adf commit becf645

17 files changed

+229
-300
lines changed

‎src/libcadet/model/GeneralRateModelDG-LinearSolver.cpp

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -221,7 +221,7 @@ void GeneralRateModelDG::assembleDiscretizedGlobalJacobian(double alpha, Indexer
221221
void GeneralRateModelDG::addTimeDerivativeToJacobianParticleShell(linalg::BandedEigenSparseRowIterator& jac, const Indexer& idxr, double alpha, unsigned int parType)
222222
{
223223
parts::cell::addTimeDerivativeToJacobianParticleShell<linalg::BandedEigenSparseRowIterator, true>(jac, alpha, static_cast<double>(_particle[parType].getPorosity()), _disc.nComp, _disc.nBound + _disc.nComp * parType,
224-
_poreAccessFactor.data() + _disc.nComp * parType, _disc.strideBound[parType], _disc.boundOffset + _disc.nComp * parType, _binding[parType]->reactionQuasiStationarity());
224+
_particle[parType].getPoreAccessfactor(), _disc.strideBound[parType], _disc.boundOffset + _disc.nComp * parType, _binding[parType]->reactionQuasiStationarity());
225225
}
226226

227227
} // namespace model

‎src/libcadet/model/GeneralRateModelDG.cpp

Lines changed: 35 additions & 91 deletions
Original file line numberDiff line numberDiff line change
@@ -375,22 +375,6 @@ bool GeneralRateModelDG::configure(IParameterProvider& paramProvider)
375375
}
376376

377377
// Read vectorial parameters (which may also be section dependent; transport)
378-
_filmDiffusionMode = readAndRegisterMultiplexCompTypeSecParam(paramProvider, _parameters, _filmDiffusion, "FILM_DIFFUSION", _disc.nParType, _disc.nComp, _unitOpIdx);
379-
380-
if ((_filmDiffusion.size() < _disc.nComp * _disc.nParType) || (_filmDiffusion.size() % (_disc.nComp * _disc.nParType) != 0))
381-
throw InvalidParameterException("Number of elements in field FILM_DIFFUSION is not a positive multiple of NCOMP * NPARTYPE (" + std::to_string(_disc.nComp * _disc.nParType) + ")");
382-
383-
if (paramProvider.exists("PORE_ACCESSIBILITY"))
384-
_poreAccessFactorMode = readAndRegisterMultiplexCompTypeSecParam(paramProvider, _parameters, _poreAccessFactor, "PORE_ACCESSIBILITY", _disc.nParType, _disc.nComp, _unitOpIdx);
385-
else
386-
{
387-
_poreAccessFactorMode = MultiplexMode::ComponentType;
388-
_poreAccessFactor = std::vector<cadet::active>(_disc.nComp * _disc.nParType, 1.0);
389-
}
390-
391-
if (_disc.nComp * _disc.nParType != _poreAccessFactor.size())
392-
throw InvalidParameterException("Number of elements in field PORE_ACCESSIBILITY differs from NCOMP * NPARTYPE (" + std::to_string(_disc.nComp * _disc.nParType) + ")");
393-
394378
// Add parameters to map
395379
_parameters[makeParamId(hashString("COL_POROSITY"), _unitOpIdx, CompIndep, ParTypeIndep, BoundStateIndep, ReactionIndep, SectionIndep)] = &_colPorosity;
396380

@@ -552,7 +536,7 @@ void GeneralRateModelDG::notifyDiscontinuousSectionTransition(double t, unsigned
552536

553537
for (int parType = 0; parType < _disc.nParType; parType++)
554538
{
555-
_particle[parType].notifyDiscontinuousSectionTransition(t, secIdx, getSectionDependentSlice(_filmDiffusion, _disc.nComp * _disc.nParType, secIdx), &_poreAccessFactor[0]);
539+
_particle[parType].notifyDiscontinuousSectionTransition(t, secIdx);
556540
}
557541

558542
_disc.curSection = secIdx;
@@ -886,14 +870,20 @@ int GeneralRateModelDG::residualImpl(double t, unsigned int secIdx, StateType co
886870
const unsigned int colNode = pblk % _disc.nPoints;
887871

888872
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
873+
model::columnPackingParameters packing
874+
{
875+
_parTypeVolFrac[parType + _disc.nParType * colNode],
876+
_colPorosity,
877+
ColumnPosition{ _convDispOp.relativeCoordinate(colNode), 0.0, 0.0 }
878+
};
890879

891880
_particle[parType].residual(t, secIdx,
892881
y + idxr.offsetCp(ParticleTypeIndex{ parType }, ParticleIndex{ colNode }),
893882
y + idxr.offsetC() + colNode * idxr.strideColNode(),
894883
yDot ? yDot + idxr.offsetCp(ParticleTypeIndex{ parType }, ParticleIndex{ colNode }) : nullptr,
895884
res ? res + idxr.offsetCp(ParticleTypeIndex{ parType }, ParticleIndex{ colNode }) : nullptr,
896-
colPos, jacIt, tlmAlloc,
885+
res ? res + idxr.offsetC() + colNode * idxr.strideColNode() : nullptr,
886+
packing, jacIt, tlmAlloc,
897887
typename cadet::ParamSens<ParamType>::enabled()
898888
);
899889
}
@@ -903,8 +893,6 @@ int GeneralRateModelDG::residualImpl(double t, unsigned int secIdx, StateType co
903893

904894
BENCH_STOP(_timerResidualPar);
905895

906-
residualFlux<StateType, ResidualType, ParamType>(t, secIdx, y, yDot, res);
907-
908896
// Handle inlet DOFs, which are simply copied to the residual
909897
for (unsigned int i = 0; i < _disc.nComp; ++i)
910898
{
@@ -964,51 +952,6 @@ int GeneralRateModelDG::residualBulk(double t, unsigned int secIdx, StateType co
964952
return 0;
965953
}
966954

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-
1012955
parts::cell::CellParameters GeneralRateModelDG::makeCellResidualParams(unsigned int parType, int const* qsReaction) const
1013956
{
1014957
return parts::cell::CellParameters
@@ -1019,7 +962,7 @@ parts::cell::CellParameters GeneralRateModelDG::makeCellResidualParams(unsigned
1019962
_disc.strideBound[parType],
1020963
qsReaction,
1021964
_particle[parType].getPorosity(),
1022-
_poreAccessFactor.data() + _disc.nComp * parType,
965+
_particle[parType].getPoreAccessfactor(),
1023966
_binding[parType],
1024967
(_dynReaction[parType] && (_dynReaction[parType]->numReactionsCombined() > 0)) ? _dynReaction[parType] : nullptr
1025968
};
@@ -1214,10 +1157,6 @@ bool GeneralRateModelDG::setParameter(const ParameterId& pId, double value)
12141157
{
12151158
if (pId.unitOperation == _unitOpIdx)
12161159
{
1217-
if (multiplexCompTypeSecParameterValue(pId, hashString("PORE_ACCESSIBILITY"), _poreAccessFactorMode, _poreAccessFactor, _disc.nParType, _disc.nComp, value, nullptr))
1218-
return true;
1219-
if (multiplexCompTypeSecParameterValue(pId, hashString("FILM_DIFFUSION"), _filmDiffusionMode, _filmDiffusion, _disc.nParType, _disc.nComp, value, nullptr))
1220-
return true;
12211160
const int mpIc = multiplexInitialConditions(pId, value, false);
12221161
if (mpIc > 0)
12231162
return true;
@@ -1241,7 +1180,12 @@ bool GeneralRateModelDG::setParameter(const ParameterId& pId, double value)
12411180
for (int parType = 0; parType < _disc.nParType; parType++)
12421181
{
12431182
if (_particle[parType].setParameter(pId, value))
1244-
return true;
1183+
{ // continue loop for particle type independent parameters to set the respective parameter sensitivity in all particle types
1184+
if ((pId.particleType != ParTypeIndep && parType == pId.particleType) || (pId.particleType == ParTypeIndep && parType == _disc.nParType - 1))
1185+
{
1186+
return true;
1187+
}
1188+
}
12451189
}
12461190

12471191
if (_convDispOp.setParameter(pId, value))
@@ -1273,7 +1217,12 @@ bool GeneralRateModelDG::setParameter(const ParameterId& pId, int value)
12731217
for (int parType = 0; parType < _disc.nParType; parType++)
12741218
{
12751219
if (_particle[parType].setParameter(pId, value))
1276-
return true;
1220+
{ // continue loop for particle type independent parameters to set the respective parameter sensitivity in all particle types
1221+
if ((pId.particleType != ParTypeIndep && parType == pId.particleType) || (pId.particleType == ParTypeIndep && parType == _disc.nParType - 1))
1222+
{
1223+
return true;
1224+
}
1225+
}
12771226
}
12781227

12791228
if (pId.unitOperation == _unitOpIdx)
@@ -1293,7 +1242,12 @@ bool GeneralRateModelDG::setParameter(const ParameterId& pId, bool value)
12931242
for (int parType = 0; parType < _disc.nParType; parType++)
12941243
{
12951244
if (_particle[parType].setParameter(pId, value))
1296-
return true;
1245+
{ // continue loop for particle type independent parameters to set the respective parameter sensitivity in all particle types
1246+
if ((pId.particleType != ParTypeIndep && parType == pId.particleType) || (pId.particleType == ParTypeIndep && parType == _disc.nParType - 1))
1247+
{
1248+
return true;
1249+
}
1250+
}
12971251
}
12981252

12991253
if (pId.unitOperation == _unitOpIdx)
@@ -1309,10 +1263,6 @@ void GeneralRateModelDG::setSensitiveParameterValue(const ParameterId& pId, doub
13091263
{
13101264
if (pId.unitOperation == _unitOpIdx)
13111265
{
1312-
if (multiplexCompTypeSecParameterValue(pId, hashString("PORE_ACCESSIBILITY"), _poreAccessFactorMode, _poreAccessFactor, _disc.nParType, _disc.nComp, value, &_sensParams))
1313-
return;
1314-
if (multiplexCompTypeSecParameterValue(pId, hashString("FILM_DIFFUSION"), _filmDiffusionMode, _filmDiffusion, _disc.nParType, _disc.nComp, value, &_sensParams))
1315-
return;
13161266
if (multiplexInitialConditions(pId, value, true) != 0)
13171267
return;
13181268

@@ -1336,7 +1286,12 @@ void GeneralRateModelDG::setSensitiveParameterValue(const ParameterId& pId, doub
13361286
for (int parType = 0; parType < _disc.nParType; parType++)
13371287
{
13381288
if (_particle[parType].setSensitiveParameterValue(_sensParams, pId, value))
1339-
return;
1289+
{ // continue loop for particle type independent parameters to set the respective parameter sensitivity in all particle types
1290+
if ((pId.particleType != ParTypeIndep && parType == pId.particleType) || (pId.particleType == ParTypeIndep && parType == _disc.nParType - 1))
1291+
{
1292+
return;
1293+
}
1294+
}
13401295
}
13411296

13421297
if (_convDispOp.setSensitiveParameterValue(_sensParams, pId, value))
@@ -1362,18 +1317,6 @@ bool GeneralRateModelDG::setSensitiveParameter(const ParameterId& pId, unsigned
13621317
{
13631318
if (pId.unitOperation == _unitOpIdx)
13641319
{
1365-
if (multiplexCompTypeSecParameterAD(pId, hashString("PORE_ACCESSIBILITY"), _poreAccessFactorMode, _poreAccessFactor, _disc.nParType, _disc.nComp, adDirection, adValue, _sensParams))
1366-
{
1367-
LOG(Debug) << "Found parameter " << pId << ": Dir " << adDirection << " is set to " << adValue;
1368-
return true;
1369-
}
1370-
1371-
if (multiplexCompTypeSecParameterAD(pId, hashString("FILM_DIFFUSION"), _filmDiffusionMode, _filmDiffusion, _disc.nParType, _disc.nComp, adDirection, adValue, _sensParams))
1372-
{
1373-
LOG(Debug) << "Found parameter " << pId << ": Dir " << adDirection << " is set to " << adValue;
1374-
return true;
1375-
}
1376-
13771320
const int mpIc = multiplexInitialConditions(pId, adDirection, adValue);
13781321
if (mpIc > 0)
13791322
{
@@ -1405,6 +1348,7 @@ bool GeneralRateModelDG::setSensitiveParameter(const ParameterId& pId, unsigned
14051348
{
14061349
if (_particle[parType].setSensitiveParameter(_sensParams, pId, adDirection, adValue))
14071350
{
1351+
// continue loop for particle type independent parameters to set the respective parameter sensitivity in all particle types
14081352
if ((pId.particleType != ParTypeIndep && parType == pId.particleType) || (pId.particleType == ParTypeIndep && parType == _disc.nParType - 1))
14091353
{
14101354
LOG(Debug) << "Found parameter " << pId << ": Dir " << adDirection << " is set to " << adValue;

‎src/libcadet/model/GeneralRateModelDG.hpp

Lines changed: 0 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -311,12 +311,6 @@ class GeneralRateModelDG : public UnitOperationBase
311311
active _colPorosity; //!< Column porosity (external porosity) \f$ \varepsilon_c \f$
312312
std::vector<active> _parTypeVolFrac; //!< Volume fraction of each particle type
313313

314-
// Vectorial parameters
315-
std::vector<active> _filmDiffusion; //!< Film diffusion coefficient \f$ k_f \f$
316-
MultiplexMode _filmDiffusionMode;
317-
std::vector<active> _poreAccessFactor; //!< Pore accessibility factor \f$ F_{\text{acc}} \f$
318-
MultiplexMode _poreAccessFactorMode;
319-
320314
bool _axiallyConstantParTypeVolFrac; //!< Determines whether particle type volume fraction is homogeneous across axial coordinate
321315
bool _analyticJac; //!< Determines whether AD or analytic Jacobians are used
322316
unsigned int _jacobianAdDirs; //!< Number of AD seed vectors required for Jacobian computation

‎src/libcadet/model/LumpedRateModelWithPoresDG-InitialConditions.cpp

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -375,7 +375,7 @@ namespace cadet
375375
_disc.strideBound[type],
376376
_binding[type]->reactionQuasiStationarity(),
377377
_particle[type].getPorosity(),
378-
_poreAccessFactor.data() + _disc.nComp * type,
378+
_particle[type].getPoreAccessfactor(),
379379
_binding[type],
380380
(_dynReaction[type] && (_dynReaction[type]->numReactionsCombined() > 0)) ? _dynReaction[type] : nullptr
381381
};

‎src/libcadet/model/LumpedRateModelWithPoresDG-LinearSolver.cpp

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -202,7 +202,7 @@ void LumpedRateModelWithPoresDG::addTimeDerivativeToJacobianParticleBlock(linalg
202202
// Add derivative with respect to dc_p / dt to Jacobian
203203
jac[0] += alpha;
204204

205-
const double invBetaP = (1.0 - static_cast<double>(_particle[parType].getPorosity())) / (static_cast<double>(_poreAccessFactor[parType * _disc.nComp + comp]) * static_cast<double>(_particle[parType].getPorosity()));
205+
const double invBetaP = (1.0 - static_cast<double>(_particle[parType].getPorosity())) / (static_cast<double>(_particle[parType].getPoreAccessfactor()[comp]) * static_cast<double>(_particle[parType].getPorosity()));
206206

207207
// Add derivative with respect to dq / dt to Jacobian
208208
const int nBound = static_cast<int>(_disc.nBound[parType * _disc.nComp + comp]);

0 commit comments

Comments
 (0)