Skip to content

Commit

Permalink
Change 2D DG from port per point to port per element
Browse files Browse the repository at this point in the history
  • Loading branch information
@jbreue16
jbreue16 committed Dec 9, 2024
1 parent 482273d commit b76f5d8
Show file tree
Hide file tree
Showing 4 changed files with 155 additions and 219 deletions.
57 changes: 30 additions & 27 deletions src/libcadet/model/LumpedRateModelWithPoresDG2D.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -449,6 +449,8 @@ bool LumpedRateModelWithPoresDG2D::configureModelDiscretization(IParameterProvid

_disc.axNPoints = _convDispOp.axNPoints();
_disc.radNPoints = _convDispOp.radNPoints();
_disc.radNNodes = _convDispOp.radNNodes();
_disc.radNElem = _convDispOp.radNElem();
_disc.nBulkPoints = _disc.axNPoints * _disc.radNPoints;

// Precompute offsets of particle type DOFs
Expand Down Expand Up @@ -626,30 +628,30 @@ bool LumpedRateModelWithPoresDG2D::configure(IParameterProvider& paramProvider)

// Let PAR_TYPE_VOLFRAC default to 1.0 for backwards compatibility
if (paramProvider.exists("PAR_TYPE_VOLFRAC"))
_parTypeVolFracMode = readAndRegisterMultiplexParam(paramProvider, _parameters, _parTypeVolFrac, "PAR_TYPE_VOLFRAC", _disc.axNPoints, _disc.radNPoints, _disc.nParType, _unitOpIdx);
_parTypeVolFracMode = readAndRegisterMultiplexParam(paramProvider, _parameters, _parTypeVolFrac, "PAR_TYPE_VOLFRAC", _disc.axNPoints, _disc.radNElem, _disc.nParType, _unitOpIdx);
else
{
// Only one particle type present
_parTypeVolFrac.resize(_disc.axNPoints * _disc.radNPoints, 1.0);
_parTypeVolFrac.resize(_disc.axNPoints * _disc.radNElem, 1.0);
_parTypeVolFracMode = MultiplexMode::Independent;
}

// Check whether all sizes are matched
// Check whether all sizes are matched
if (_disc.nParType != _parRadius.size())
throw InvalidParameterException("Number of elements in field PAR_RADIUS does not match number of particle types");
if (_disc.nParType * _disc.nBulkPoints != _parTypeVolFrac.size())
throw InvalidParameterException("Number of elements in field PAR_TYPE_VOLFRAC does not match number of particle types");
if (_disc.nParType * _disc.axNPoints * _disc.radNElem != _parTypeVolFrac.size())
throw InvalidParameterException("Number of elements in field PAR_TYPE_VOLFRAC does not match number of particle types times RAD_NELEM * (AX_POLYDEG+1)*AX_NELEM");
if (_disc.nParType != _parPorosity.size())
throw InvalidParameterException("Number of elements in field PAR_POROSITY does not match number of particle types");

// Check that particle volume fractions sum to 1.0
for (unsigned int i = 0; i < _disc.axNPoints * _disc.radNPoints; ++i)
for (unsigned int i = 0; i < _disc.axNPoints * _disc.radNElem; ++i)
{
const double volFracSum = std::accumulate(_parTypeVolFrac.begin() + i * _disc.nParType, _parTypeVolFrac.begin() + (i+1) * _disc.nParType, 0.0,
[](double a, const active& b) -> double { return a + static_cast<double>(b); });
if (std::abs(1.0 - volFracSum) > 1e-10)
throw InvalidParameterException("Sum of field PAR_TYPE_VOLFRAC differs from 1.0 (is " + std::to_string(volFracSum) + ") in axial cell " + std::to_string(i / _disc.radNPoints) + " radial cell " + std::to_string(i % _disc.radNPoints));
throw InvalidParameterException("Sum of field PAR_TYPE_VOLFRAC differs from 1.0 (is " + std::to_string(volFracSum) + ") in axial element " + std::to_string(i / _disc.radNElem) + " radial element " + std::to_string(i % _disc.radNElem));
}

// Read vectorial parameters (which may also be section dependent; transport)
Expand All @@ -671,17 +673,17 @@ bool LumpedRateModelWithPoresDG2D::configure(IParameterProvider& paramProvider)

// Register initial conditions parameters
registerParam1DArray(_parameters, _initC, [=](bool multi, unsigned int comp) { return makeParamId(hashString("INIT_C"), _unitOpIdx, comp, ParTypeIndep, BoundStateIndep, ReactionIndep, SectionIndep); });
if (_disc.radNPoints > 1)
if (_disc.radNElem > 1)
registerParam2DArray(_parameters, _initC, [=](bool multi, unsigned int rad, unsigned int comp) { return makeParamId(hashString("INIT_C"), _unitOpIdx, comp, ParTypeIndep, BoundStateIndep, rad, SectionIndep); }, _disc.nComp);

if (_singleBinding)
{
for (unsigned int c = 0; c < _disc.nComp; ++c)
_parameters[makeParamId(hashString("INIT_CP"), _unitOpIdx, c, ParTypeIndep, BoundStateIndep, ReactionIndep, SectionIndep)] = &_initCp[c];

if (_disc.radNPoints > 1)
if (_disc.radNElem > 1)
{
for (unsigned int r = 0; r < _disc.radNPoints; ++r)
for (unsigned int r = 0; r < _disc.radNElem; ++r)
{
for (unsigned int c = 0; c < _disc.nComp; ++c)
_parameters[makeParamId(hashString("INIT_CP"), _unitOpIdx, c, ParTypeIndep, BoundStateIndep, r, SectionIndep)] = &_initCp[r * _disc.nComp * _disc.nParType + c];
Expand All @@ -691,7 +693,7 @@ bool LumpedRateModelWithPoresDG2D::configure(IParameterProvider& paramProvider)
else
{
registerParam2DArray(_parameters, _initCp, [=](bool multi, unsigned int type, unsigned int comp) { return makeParamId(hashString("INIT_CP"), _unitOpIdx, comp, type, BoundStateIndep, ReactionIndep, SectionIndep); }, _disc.nComp);
if (_disc.radNPoints > 1)
if (_disc.radNElem > 1)
registerParam3DArray(_parameters, _initCp, [=](bool multi, unsigned int rad, unsigned int type, unsigned int comp) { return makeParamId(hashString("INIT_CP"), _unitOpIdx, comp, type, BoundStateIndep, rad, SectionIndep); }, _disc.nComp, _disc.nParType);
}

Expand Down Expand Up @@ -722,13 +724,13 @@ bool LumpedRateModelWithPoresDG2D::configure(IParameterProvider& paramProvider)
}

// Register radially dependent
if (_disc.radNPoints > 1)
if (_disc.radNElem > 1)
{
if (_singleBinding)
{
_binding[0]->fillBoundPhaseInitialParameters(initParams.data(), _unitOpIdx, ParTypeIndep);

for (unsigned int r = 0; r < _disc.radNPoints; ++r)
for (unsigned int r = 0; r < _disc.radNElem; ++r)
{
for (ParameterId& pId : initParams)
pId.reaction = r;
Expand All @@ -740,7 +742,7 @@ bool LumpedRateModelWithPoresDG2D::configure(IParameterProvider& paramProvider)
}
else
{
for (unsigned int r = 0; r < _disc.radNPoints; ++r)
for (unsigned int r = 0; r < _disc.radNElem; ++r)
{
for (unsigned int type = 0; type < _disc.nParType; ++type)
{
Expand Down Expand Up @@ -1303,7 +1305,8 @@ int LumpedRateModelWithPoresDG2D::residualFlux(double t, unsigned int secIdx, St

for (unsigned int j = 0; j < _disc.radNPoints; ++j)
{
invBetaC = 1.0 / static_cast<ParamType>(_convDispOp.columnPorosity(j)) - 1.0;
const int radZone = std::floor(j / _disc.radNNodes);
invBetaC = 1.0 / static_cast<ParamType>(_convDispOp.columnPorosity(radZone)) - 1.0;
jacCF_val = invBetaC * _parGeomSurfToVol[type] / radius;
jacPF_val = -_parGeomSurfToVol[type] / (epsP * radius);

Expand All @@ -1319,7 +1322,7 @@ int LumpedRateModelWithPoresDG2D::residualFlux(double t, unsigned int secIdx, St
if (wantRes)
{
// Add flux to column void / bulk volume equations
resCol[ClIdx] += jacCF_val * static_cast<ParamType>(filmDiff[comp]) * static_cast<ParamType>(_parTypeVolFrac[type + i * _disc.nParType * _disc.radNPoints + j * _disc.nParType]) * (yCol[ClIdx] - yCol[CpIdx]);
resCol[ClIdx] += jacCF_val * static_cast<ParamType>(filmDiff[comp]) * static_cast<ParamType>(_parTypeVolFrac[type + i * _disc.nParType * _disc.radNElem + radZone * _disc.nParType]) * (yCol[ClIdx] - yCol[CpIdx]);

// Add flux to particle / bead volume equations
resCol[CpIdx] += jacPF_val / static_cast<ParamType>(poreAccFactor[comp]) * static_cast<ParamType>(filmDiff[comp]) * (yCol[ClIdx] - yCol[CpIdx]);
Expand All @@ -1330,11 +1333,11 @@ int LumpedRateModelWithPoresDG2D::residualFlux(double t, unsigned int secIdx, St
// add Cl on Cl entries (added since entries already set in transport jacobian)
// row: already at current bulk node and component
// col: already at current bulk node and component
_globalJac.coeffRef(ClIdx, ClIdx) += static_cast<double>(jacCF_val) * static_cast<double>(filmDiff[comp]) * static_cast<double>(_parTypeVolFrac[type + i * _disc.nParType * _disc.radNPoints + j * _disc.nParType]);
_globalJac.coeffRef(ClIdx, ClIdx) += static_cast<double>(jacCF_val) * static_cast<double>(filmDiff[comp]) * static_cast<double>(_parTypeVolFrac[type + i * _disc.nParType * _disc.radNElem + radZone * _disc.nParType]);
// add Cl on Cp entries
// row: already at current bulk node and component
// col: jump to particle phase
_globalJac.coeffRef(ClIdx, CpIdx) = -static_cast<double>(jacCF_val) * static_cast<double>(filmDiff[comp]) * static_cast<double>(_parTypeVolFrac[type + i * _disc.nParType * _disc.radNPoints + j * _disc.nParType]);
_globalJac.coeffRef(ClIdx, CpIdx) = -static_cast<double>(jacCF_val) * static_cast<double>(filmDiff[comp]) * static_cast<double>(_parTypeVolFrac[type + i * _disc.nParType * _disc.radNElem + radZone * _disc.nParType]);

// add Cp on Cp entries (added since entries already set in particle jacobian)
// row: already at particle. already at current node and liquid state
Expand Down Expand Up @@ -1526,10 +1529,10 @@ unsigned int LumpedRateModelWithPoresDG2D::localOutletComponentIndex(unsigned in
{
// Inlets are duplicated so need to be accounted for
if (static_cast<double>(_convDispOp.currentVelocity(port)) >= 0.0)
// Forward Flow: outlet is last cell
// Forward Flow: outlet is last node
return _disc.nComp * _disc.radNPoints + (_disc.axNPoints - 1) * _disc.nComp * _disc.radNPoints + port * _disc.nComp;
else
// Backward flow: Outlet is first cell
// Backward flow: Outlet is first node
return _disc.nComp * _disc.radNPoints + _disc.nComp * port;
}

Expand Down Expand Up @@ -1558,7 +1561,7 @@ bool LumpedRateModelWithPoresDG2D::setParameter(const ParameterId& pId, double v
{
if (pId.unitOperation == _unitOpIdx)
{
if (multiplexParameterValue(pId, hashString("PAR_TYPE_VOLFRAC"), _parTypeVolFracMode, _parTypeVolFrac, _disc.axNPoints, _disc.radNPoints, _disc.nParType, value, nullptr))
if (multiplexParameterValue(pId, hashString("PAR_TYPE_VOLFRAC"), _parTypeVolFracMode, _parTypeVolFrac, _disc.axNPoints, _disc.radNElem, _disc.nParType, value, nullptr))
return true;
if (multiplexCompTypeSecParameterValue(pId, hashString("PORE_ACCESSIBILITY"), _poreAccessFactorMode, _poreAccessFactor, _disc.nParType, _disc.nComp, value, nullptr))
return true;
Expand Down Expand Up @@ -1588,7 +1591,7 @@ void LumpedRateModelWithPoresDG2D::setSensitiveParameterValue(const ParameterId&
{
if (pId.unitOperation == _unitOpIdx)
{
if (multiplexParameterValue(pId, hashString("PAR_TYPE_VOLFRAC"), _parTypeVolFracMode, _parTypeVolFrac, _disc.axNPoints, _disc.radNPoints, _disc.nParType, value, &_sensParams))
if (multiplexParameterValue(pId, hashString("PAR_TYPE_VOLFRAC"), _parTypeVolFracMode, _parTypeVolFrac, _disc.axNPoints, _disc.radNElem, _disc.nParType, value, &_sensParams))
return;
if (multiplexCompTypeSecParameterValue(pId, hashString("PORE_ACCESSIBILITY"), _poreAccessFactorMode, _poreAccessFactor, _disc.nParType, _disc.nComp, value, &_sensParams))
return;
Expand All @@ -1613,7 +1616,7 @@ bool LumpedRateModelWithPoresDG2D::setSensitiveParameter(const ParameterId& pId,
{
if (pId.unitOperation == _unitOpIdx)
{
if (multiplexParameterAD(pId, hashString("PAR_TYPE_VOLFRAC"), _parTypeVolFracMode, _parTypeVolFrac, _disc.axNPoints, _disc.radNPoints, _disc.nParType, adDirection, adValue, _sensParams))
if (multiplexParameterAD(pId, hashString("PAR_TYPE_VOLFRAC"), _parTypeVolFracMode, _parTypeVolFrac, _disc.axNPoints, _disc.radNElem, _disc.nParType, adDirection, adValue, _sensParams))
{
LOG(Debug) << "Found parameter " << pId << ": Dir " << adDirection << " is set to " << adValue;
return true;
Expand Down Expand Up @@ -1728,8 +1731,8 @@ int LumpedRateModelWithPoresDG2D::Exporter::writeParticleMobilePhase(unsigned in

int LumpedRateModelWithPoresDG2D::Exporter::writeInlet(unsigned int port, double* buffer) const
{
cadet_assert(port < _disc.radNPoints);
std::copy_n(_data + port * _disc.nComp, _disc.nComp, buffer);
cadet_assert(port < _disc.radNElem);
std::copy_n(_data + port * _disc.nComp * _disc.radNNodes, _disc.nComp * _disc.radNNodes, buffer);
return _disc.nComp;
}

Expand All @@ -1741,12 +1744,12 @@ int LumpedRateModelWithPoresDG2D::Exporter::writeInlet(double* buffer) const

int LumpedRateModelWithPoresDG2D::Exporter::writeOutlet(unsigned int port, double* buffer) const
{
cadet_assert(port < _disc.radNPoints);
cadet_assert(port < _disc.radNElem);

if (_model._convDispOp.currentVelocity(port) >= 0)
std::copy_n(&_idx.c(_data, _disc.axNPoints - 1, port, 0), _disc.nComp, buffer);
std::copy_n(&_idx.c(_data, _disc.axNPoints - 1, port * _disc.radNNodes, 0), _disc.nComp * _disc.radNNodes, buffer);
else
std::copy_n(&_idx.c(_data, 0, port, 0), _disc.nComp, buffer);
std::copy_n(&_idx.c(_data, 0, port * _disc.radNNodes, 0), _disc.nComp * _disc.radNNodes, buffer);

return _disc.nComp;
}
Expand Down
10 changes: 6 additions & 4 deletions src/libcadet/model/LumpedRateModelWithPoresDG2D.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -78,8 +78,8 @@ class LumpedRateModelWithPoresDG2D : public UnitOperationBase
virtual UnitOpIdx unitOperationId() const CADET_NOEXCEPT { return _unitOpIdx; }
virtual unsigned int numComponents() const CADET_NOEXCEPT { return _disc.nComp; }
virtual void setFlowRates(active const* in, active const* out) CADET_NOEXCEPT;
virtual unsigned int numInletPorts() const CADET_NOEXCEPT { return _disc.radNPoints; }
virtual unsigned int numOutletPorts() const CADET_NOEXCEPT { return _disc.radNPoints; }
virtual unsigned int numInletPorts() const CADET_NOEXCEPT { return _disc.radNElem; }
virtual unsigned int numOutletPorts() const CADET_NOEXCEPT { return _disc.radNElem; }
virtual bool canAccumulate() const CADET_NOEXCEPT { return false; }

static const char* identifier() { return "LUMPED_RATE_MODEL_WITH_PORES_2D"; }
Expand Down Expand Up @@ -230,6 +230,8 @@ class LumpedRateModelWithPoresDG2D : public UnitOperationBase
unsigned int nComp; //!< Number of components
unsigned int axNPoints; //!< Number of axial discrete points
unsigned int radNPoints; //!< Number of radial discrete points
unsigned int radNElem; //!< Number of radial DG elements, which can be used as radial zones for radial parameter inhomogeneities
unsigned int radNNodes; //!< Number of radial DG nodes per element
unsigned int nBulkPoints; //!< Number of total bulk discrete points
unsigned int nParType; //!< Number of particle types
unsigned int* parTypeOffset; //!< Array with offsets (in particle block) to particle type, additional last element contains total number of particle DOFs
Expand Down Expand Up @@ -375,8 +377,8 @@ class LumpedRateModelWithPoresDG2D : public UnitOperationBase
virtual unsigned int numComponents() const CADET_NOEXCEPT { return _disc.nComp; }
virtual unsigned int numPrimaryCoordinates() const CADET_NOEXCEPT { return _disc.axNPoints; }
virtual unsigned int numSecondaryCoordinates() const CADET_NOEXCEPT { return _disc.radNPoints; }
virtual unsigned int numInletPorts() const CADET_NOEXCEPT { return _disc.radNPoints; }
virtual unsigned int numOutletPorts() const CADET_NOEXCEPT { return _disc.radNPoints; }
virtual unsigned int numInletPorts() const CADET_NOEXCEPT { return _disc.radNElem; }
virtual unsigned int numOutletPorts() const CADET_NOEXCEPT { return _disc.radNElem; }
virtual unsigned int numParticleTypes() const CADET_NOEXCEPT { return _disc.nParType; }
virtual unsigned int numParticleShells(unsigned int parType) const CADET_NOEXCEPT { return 1; }
virtual unsigned int numBoundStates(unsigned int parType) const CADET_NOEXCEPT { return _disc.strideBound[parType]; }
Expand Down
Loading

0 comments on commit b76f5d8

Please sign in to comment.