Skip to content

Deprecate collocation variant of DGSEM particle discretization #411

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
jbreue16 merged 5 commits into from
May 15, 2025
Merged

Deprecate collocation variant of DGSEM particle discretization #411

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
d4c6796
Add numerical reference test for 2DLRMP to CI
jbreue16 May 6, 2025
5542647
Rename nCol to nElem in DG units
jbreue16 May 6, 2025
dce44c5
Update dev guide: test flags
jbreue16 May 8, 2025
9ab8665
Improve documentation of sensitivity interface
jbreue16 May 13, 2025
dcca54e
Deprecate collocation variant of DGSEM particle discretization (#411)
jbreue16 May 12, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions doc/developer_guide/testing.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,8 @@ To debug specific tests (with flag [testHere]) from the Visual Studio IDE, you c
}

Select the testRunner.exe in the startup item dropdown menu and you can start debugging tests with the specified flag.
To see the options, run `testrunner --help`. To get the available test flags for instance, execute the testrunner with the corresponding option, i.e. `testrunner -t`.
You can provide multiple arguments like `"[testHere][testHere2]"` to run all test that have both flags, or with comma separation `"[testHere],[testHere2]"` to run all tests with the first and/or the second flag.

Adding tests for your model
---------------------------
Expand Down
6 changes: 3 additions & 3 deletions doc/interface/sensitivities.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,23 +39,23 @@ Group /input/sensitivity/param_XXX

``SENS_NAME``

Name of the parameter (Note that ``PAR_RADIUS`` and ``PAR_CORE_RADIUS`` sensitivities are only available for Finite Volume discretization)
Name of the parameter (Note that ``PAR_RADIUS`` and ``PAR_CORERADIUS`` sensitivities are only available for Finite Volume discretization)

================ ===========================
**Type:** string **Length:** :math:`\geq 1`
================ ===========================

``SENS_COMP``

Component index (:math:`-1` if parameter is independent of components)
Component index (:math:`-1` if parameter is independent of components, see the multiplexing of the corresponding parameter, e.g. `FILM_DIFFUSION_MULTIPLEX`)

============= ========================== ============================
**Type:** int **Range:** :math:`\geq -1` **Length:** :math:`\geq 1`
============= ========================== ============================

``SENS_PARTYPE``

Particle type index (:math:`-1` if parameter is independent of particle types)
Particle type index (:math:`-1` if parameter is independent of particle types, see the multiplexing of the corresponding parameter, e.g. `PAR_DIFFUSION_MULTIPLEX`)

============= ========================== ===========================
**Type:** int **Range:** :math:`\geq -1` **Length:** :math:`\geq 1`
Expand Down
9 changes: 0 additions & 9 deletions src/libcadet/model/GeneralRateModelDG-InitialConditions.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -303,15 +303,6 @@ void GeneralRateModelDG::consistentInitialState(const SimulationTime& simTime, d

Indexer idxr(_disc);

// initialization for inexact integration DG discretization of particle mass balance not supported. This would require the consideration of the additional algebraic constraints,
// but the general performance (stiffness due to the additional algebraic constraints) of this scheme does not justify the effort here.
for (unsigned int type = 0; type < _disc.nParType; type++) {
if (_disc.parExactInt[type] == false) {
LOG(Error) << "No consistent initialization for inexact integration DG discretization in particles (cf. par_exact_integration). If consistent initialization is required, change to exact integration.";
return;
}
}

// Step 1: Solve algebraic equations

// Step 1a: Compute quasi-stationary binding model state
Expand Down
30 changes: 2 additions & 28 deletions src/libcadet/model/GeneralRateModelDG-LinearSolver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -141,7 +141,7 @@ int GeneralRateModelDG::linearSolve(double t, double alpha, double outerTol, dou

// Handle inlet DOFs:
// Inlet at z = 0 for forward flow, at z = L for backward flow.
unsigned int offInlet = _convDispOp.forwardFlow() ? 0 : (_disc.nCol - 1u) * idxr.strideColCell();
unsigned int offInlet = _convDispOp.forwardFlow() ? 0 : (_disc.nElem - 1u) * idxr.strideColCell();

for (int comp = 0; comp < _disc.nComp; comp++) {
for (int node = 0; node < (_disc.exactInt ? _disc.nNodes : 1); node++) {
Expand Down Expand Up @@ -189,33 +189,7 @@ void GeneralRateModelDG::assembleDiscretizedGlobalJacobian(double alpha, Indexer

linalg::BandedEigenSparseRowIterator jac(_globalJacDisc, idxr.offsetCp(ParticleTypeIndex{ parType }, ParticleIndex{ colNode }));

// If special case for inexact integration DG scheme:
// Do not add time derivative to particle mass balance equation at inner particle boundary for mass balance(s)
if (!_disc.parExactInt[parType] && _parGeomSurfToVol[parType] != _disc.SurfVolRatioSlab && _parCoreRadius[parType] == 0.0) {
// we still need to add the derivative for the binding
// move iterator to solid phase
jac += idxr.strideParLiquid();
// Solid phase
for (unsigned int comp = 0; comp < _disc.nComp; ++comp) {
for (unsigned int bnd = 0; bnd < _disc.nBound[parType * _disc.nComp + comp]; ++bnd, ++jac)
{
// Add derivative with respect to dynamic states to Jacobian
if (_binding[parType]->reactionQuasiStationarity()[idxr.offsetBoundComp(ParticleTypeIndex{ parType }, ComponentIndex{ comp }) + bnd])
continue;

// surface diffusion + kinetic binding -> additional DG-discretized mass balance equation for solid, for which the (inexact integration) discretization special case also holds
else if (_hasSurfaceDiffusion[parType]
&& static_cast<double>(getSectionDependentSlice(_parSurfDiffusion, _disc.strideBound[_disc.nParType], _disc.curSection)[_disc.nBoundBeforeType[parType] + getOffsetSurfDiff(parType, comp, bnd)]) != 0.0)
continue;

// Add derivative with respect to dq / dt to Jacobian
jac[0] += alpha;
}
}
}
else { // else, treat boundary node "normally"
addTimeDerivativeToJacobianParticleShell(jac, idxr, alpha, parType);
}
addTimeDerivativeToJacobianParticleShell(jac, idxr, alpha, parType);

// compute time derivative of remaining points
// Iterator jac has already been advanced to next shell
Expand Down
95 changes: 13 additions & 82 deletions src/libcadet/model/GeneralRateModelDG.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -157,16 +157,16 @@ bool GeneralRateModelDG::configureModelDiscretization(IParameterProvider& paramP
_disc.nNodes = _disc.polyDeg + 1;

if (paramProvider.exists("NELEM"))
_disc.nCol = paramProvider.getInt("NELEM");
_disc.nElem = paramProvider.getInt("NELEM");
else if (paramProvider.exists("NCOL"))
_disc.nCol = std::max(1u, paramProvider.getInt("NCOL") / _disc.nNodes); // number of elements is rounded down
_disc.nElem = std::max(1u, paramProvider.getInt("NCOL") / _disc.nNodes); // number of elements is rounded down
else
throw InvalidParameterException("Specify field NELEM (or NCOL)");

if (_disc.nCol < 1)
if (_disc.nElem < 1)
throw InvalidParameterException("Number of column elements must be at least 1!");

_disc.nPoints = _disc.nNodes * _disc.nCol;
_disc.nPoints = _disc.nNodes * _disc.nElem;

int polynomial_integration_mode = 0;
if (paramProvider.exists("EXACT_INTEGRATION"))
Expand All @@ -186,12 +186,10 @@ bool GeneralRateModelDG::configureModelDiscretization(IParameterProvider& paramP

std::vector<int> parPolyDeg(_disc.nParType);
std::vector<int> ParNelem(_disc.nParType);
std::vector<bool> parExactInt(_disc.nParType, true);
if (firstConfigCall)
{
_disc.parPolyDeg = new unsigned int[_disc.nParType];
_disc.nParCell = new unsigned int[_disc.nParType];
_disc.parExactInt = new bool[_disc.nParType];
_disc.parGSM = new bool[_disc.nParType];
}

Expand All @@ -206,12 +204,6 @@ bool GeneralRateModelDG::configureModelDiscretization(IParameterProvider& paramP
if ((std::any_of(ParNelem.begin(), ParNelem.end(), [](int value) { return value < 1; })))
throw InvalidParameterException("Particle number of elements must be at least 1!");

if(paramProvider.exists("PAR_EXACT_INTEGRATION"))
parExactInt = paramProvider.getBoolArray("PAR_EXACT_INTEGRATION");

if ((std::any_of(parExactInt.begin(), parExactInt.end(), [](bool value) { return !value; })))
LOG(Warning) << "Inexact integration method (cf. PAR_EXACT_INTEGRATION) in particles might add severe! stiffness to the system and disables consistent initialization!";

if (parPolyDeg.size() == 1)
{
// Multiplex number of particle shells to all particle types
Expand All @@ -232,16 +224,6 @@ bool GeneralRateModelDG::configureModelDiscretization(IParameterProvider& paramP
throw InvalidParameterException("Field PAR_NELEM must have 1 or NPARTYPE (" + std::to_string(_disc.nParType) + ") entries");
else
std::copy_n(ParNelem.begin(), _disc.nParType, _disc.nParCell);
if (parExactInt.size() == 1)
{
// Multiplex number of particle shells to all particle types
for (unsigned int i = 0; i < _disc.nParType; ++i)
std::fill(_disc.parExactInt, _disc.parExactInt + _disc.nParType, parExactInt[0]);
}
else if (parExactInt.size() < _disc.nParType)
throw InvalidParameterException("Field PAR_EXACT_INTEGRATION must have 1 or NPARTYPE (" + std::to_string(_disc.nParType) + ") entries");
else
std::copy_n(parExactInt.begin(), _disc.nParType, _disc.parExactInt);
}
else if (paramProvider.exists("NPAR"))
{
Expand Down Expand Up @@ -507,7 +489,7 @@ bool GeneralRateModelDG::configureModelDiscretization(IParameterProvider& paramP
}

unsigned int strideColNode = _disc.nComp;
const bool transportSuccess = _convDispOp.configureModelDiscretization(paramProvider, helper, _disc.nComp, polynomial_integration_mode, _disc.nCol, _disc.polyDeg, strideColNode);
const bool transportSuccess = _convDispOp.configureModelDiscretization(paramProvider, helper, _disc.nComp, polynomial_integration_mode, _disc.nElem, _disc.polyDeg, strideColNode);

_disc.curSection = -1;

Expand Down Expand Up @@ -1366,9 +1348,6 @@ int GeneralRateModelDG::residualParticle(double t, unsigned int parType, unsigne

LinearBufferAllocator tlmAlloc = threadLocalMem.get();

// Special case: individual treatment of time derivatives in particle mass balance at inner particle boundary node
bool specialCase = !_disc.parExactInt[parType] && (_parGeomSurfToVol[parType] != _disc.SurfVolRatioSlab && _parCoreRadius[parType] == 0.0);

// Prepare parameters
active const* const parDiff = getSectionDependentSlice(_parDiffusion, _disc.nComp * _disc.nParType, secIdx) + parType * _disc.nComp;

Expand Down Expand Up @@ -1407,62 +1386,14 @@ int GeneralRateModelDG::residualParticle(double t, unsigned int parType, unsigne
);
const ColumnPosition colPos{ z, 0.0, r };

// Handle time derivatives, binding, dynamic reactions.
// Special case: Dont add time derivatives to inner boundary node for DG discretized mass balance equations.
// This can be achieved by setting yDot pointer to null before passing to residual kernel, and adding only the time derivative for dynamic binding
// TODO Check Treatment of reactions (do we need yDot then?)
if (cadet_unlikely(par == 0 && specialCase))
{
if (wantRes)
parts::cell::residualKernel<StateType, ResidualType, ParamType, parts::cell::CellParameters, linalg::BandedEigenSparseRowIterator, wantJac, true>(
t, secIdx, colPos, local_y, nullptr, local_res, jac, cellResParams, tlmAlloc // TODO Check Treatment of reactions (do we need yDot then?)
);
else
parts::cell::residualKernel<StateType, ResidualType, ParamType, parts::cell::CellParameters, linalg::BandedEigenSparseRowIterator, wantJac, false, false>(
t, secIdx, colPos, local_y, nullptr, local_res, jac, cellResParams, tlmAlloc // TODO Check Treatment of reactions (do we need yDot then?)
);

if (!wantRes)
continue;

if (cellResParams.binding->hasDynamicReactions() && local_yDot)
{
unsigned int idx = 0;
for (unsigned int comp = 0; comp < cellResParams.nComp; ++comp)
{
for (unsigned int state = 0; state < cellResParams.nBound[comp]; ++state, ++idx)
{
// Skip quasi-stationary fluxes
if (cellResParams.qsReaction[idx])
continue;

// For kinetic bindings and surface diffusion, we have an additional DG-discretized mass balance eq.
// -> add time derivate at inner bonudary node only without surface diffusion
else if (_hasSurfaceDiffusion[parType])
continue;
// Some bound states might still not be effected by surface diffusion
else if (parSurfDiff[idx] != 0.0)
continue;

// Add time derivative to solid phase
local_res[idxr.strideParLiquid() + idx] += local_yDot[idxr.strideParLiquid() + idx];
}
}
}
}
if (wantRes)
parts::cell::residualKernel<StateType, ResidualType, ParamType, parts::cell::CellParameters, linalg::BandedEigenSparseRowIterator, wantJac, true>(
t, secIdx, colPos, local_y, local_yDot, local_res, jac, cellResParams, tlmAlloc
);
else
{
if (wantRes)
parts::cell::residualKernel<StateType, ResidualType, ParamType, parts::cell::CellParameters, linalg::BandedEigenSparseRowIterator, wantJac, true>(
t, secIdx, colPos, local_y, local_yDot, local_res, jac, cellResParams, tlmAlloc
);
else
{
parts::cell::residualKernel<StateType, ResidualType, ParamType, parts::cell::CellParameters, linalg::BandedEigenSparseRowIterator, wantJac, false, false>(
t, secIdx, colPos, local_y, local_yDot, local_res, jac, cellResParams, tlmAlloc
);
}
}
parts::cell::residualKernel<StateType, ResidualType, ParamType, parts::cell::CellParameters, linalg::BandedEigenSparseRowIterator, wantJac, false, false>(
t, secIdx, colPos, local_y, local_yDot, local_res, jac, cellResParams, tlmAlloc
);

// Move rowiterator to next particle node
jac += idxr.strideParNode(parType);
Expand Down Expand Up @@ -1797,7 +1728,7 @@ void GeneralRateModelDG::multiplyWithJacobian(const SimulationTime& simTime, con

// Map inlet DOFs to the column inlet (first bulk cells)
// Inlet at z = 0 for forward flow, at z = L for backward flow.
unsigned int offInlet = _convDispOp.forwardFlow() ? 0 : (_disc.nCol - 1u) * idxr.strideColCell();
unsigned int offInlet = _convDispOp.forwardFlow() ? 0 : (_disc.nElem - 1u) * idxr.strideColCell();

for (unsigned int comp = 0; comp < _disc.nComp; comp++) {
for (unsigned int node = 0; node < (_disc.exactInt ? _disc.nNodes : 1); node++) {
Expand Down
Loading
Loading