Skip to content

Modularize particle DG discretization #406

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
merged 2 commits into from
May 19, 2025
Merged

Modularize particle DG discretization #406

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
d93381a
Modularize particle DG discretization (#406)
jbreue16 May 15, 2025
ab452f2
Documentation: Fix indentation in linux install guide
jbreue16 May 19, 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
1 change: 1 addition & 0 deletions BUILD-LINUX.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,7 @@ sudo apt -y install liblapack3 liblapack-dev libblas3 libblas-dev
- Open a terminal and change to `CADET-Core/build`
- If using MKL, execute `export MKLROOT=/opt/intel/mkl`
- To compile:

- Using standard LAPACK: Execute `cmake -DCMAKE_INSTALL_PREFIX="../install" ../`

- Using MKL (sequential): Execute `cmake -DCMAKE_INSTALL_PREFIX="../install" -DBLA_VENDOR=Intel10_64lp_seq ../`
Expand Down
1 change: 1 addition & 0 deletions src/libcadet/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -214,6 +214,7 @@ if (ENABLE_DG)
list (APPEND LIBCADET_SOURCES
${CMAKE_SOURCE_DIR}/src/libcadet/model/parts/DGToolbox.cpp
${CMAKE_SOURCE_DIR}/src/libcadet/model/parts/ConvectionDispersionOperatorDG.cpp
${CMAKE_SOURCE_DIR}/src/libcadet/model/parts/ParticleDiffusionOperatorDG.cpp
${CMAKE_SOURCE_DIR}/src/libcadet/model/LumpedRateModelWithPoresDG.cpp
${CMAKE_SOURCE_DIR}/src/libcadet/model/LumpedRateModelWithoutPoresDG.cpp
${CMAKE_SOURCE_DIR}/src/libcadet/model/GeneralRateModelDG.cpp
Expand Down
28 changes: 12 additions & 16 deletions src/libcadet/model/GeneralRateModelDG-InitialConditions.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -386,8 +386,7 @@ void GeneralRateModelDG::consistentInitialState(const SimulationTime& simTime, d
active* const localAdY = adJac.adY ? adJac.adY + localOffsetToParticle + localOffsetInParticle : nullptr;

// r (particle) coordinate of current node
const double r = static_cast<double>(_disc.deltaR[type]) * std::floor(node / _disc.nParNode[type])
+ 0.5 * static_cast<double>(_disc.deltaR[type]) * (1 + _disc.parNodes[type][node % _disc.nParNode[type]]);
const double r = _parDiffOp.relativeCoordinate(type, node);
const ColumnPosition colPos{ z, 0.0, r };

// Determine whether nonlinear solver is required
Expand All @@ -398,8 +397,8 @@ void GeneralRateModelDG::consistentInitialState(const SimulationTime& simTime, d
linalg::selectVectorSubset(qShell - _disc.nComp, mask, solution);

// Save values of conserved moieties
const double epsQ = 1.0 - static_cast<double>(_parPorosity[type]);
linalg::conservedMoietiesFromPartitionedMask(mask, _disc.nBound + type * _disc.nComp, _disc.nComp, qShell - _disc.nComp, conservedQuants, static_cast<double>(_parPorosity[type]), epsQ);
const double epsQ = 1.0 - static_cast<double>(_parDiffOp._parPorosity[type]);
linalg::conservedMoietiesFromPartitionedMask(mask, _disc.nBound + type * _disc.nComp, _disc.nComp, qShell - _disc.nComp, conservedQuants, static_cast<double>(_parDiffOp._parPorosity[type]), epsQ);

std::function<bool(double const* const, linalg::detail::DenseMatrixBase&)> jacFunc;

Expand Down Expand Up @@ -478,7 +477,7 @@ void GeneralRateModelDG::consistentInitialState(const SimulationTime& simTime, d
continue;
}

mat.native(rIdx, rIdx) = static_cast<double>(_parPorosity[type]);
mat.native(rIdx, rIdx) = static_cast<double>(_parDiffOp._parPorosity[type]);

for (unsigned int bnd = 0; bnd < _disc.nBound[_disc.nComp * type + comp]; ++bnd, ++bndIdx)
{
Expand Down Expand Up @@ -526,7 +525,7 @@ void GeneralRateModelDG::consistentInitialState(const SimulationTime& simTime, d
continue;
}

mat.native(rIdx, rIdx) = static_cast<double>(_parPorosity[type]);
mat.native(rIdx, rIdx) = static_cast<double>(_parDiffOp._parPorosity[type]);

for (unsigned int bnd = 0; bnd < _disc.nBound[_disc.nComp * type + comp]; ++bnd, ++bndIdx)
{
Expand Down Expand Up @@ -573,7 +572,7 @@ void GeneralRateModelDG::consistentInitialState(const SimulationTime& simTime, d
continue;
}

r[rIdx] = static_cast<double>(_parPorosity[type]) * x[rIdx] - conservedQuants[rIdx];
r[rIdx] = static_cast<double>(_parDiffOp._parPorosity[type]) * x[rIdx] - conservedQuants[rIdx];

for (unsigned int bnd = 0; bnd < _disc.nBound[_disc.nComp * type + comp]; ++bnd, ++bndIdx)
{
Expand Down Expand Up @@ -714,9 +713,8 @@ void GeneralRateModelDG::consistentInitialTimeDerivative(const SimulationTime& s
if (_binding[type]->dependsOnTime())
{
// r (particle) coordinate of current node (particle radius normed to 1) - needed in externally dependent adsorption kinetic
const double r = (static_cast<double>(_disc.deltaR[type]) * std::floor(j / _disc.nParNode[type])
+ 0.5 * static_cast<double>(_disc.deltaR[type]) * (1 + _disc.parNodes[type][j % _disc.nParNode[type]]))
/ (static_cast<double>(_parRadius[type]) - static_cast<double>(_parCoreRadius[type]));
const double r = _parDiffOp.relativeCoordinate(type, j);

_binding[type]->timeDerivativeQuasiStationaryFluxes(simTime.t, simTime.secIdx,
ColumnPosition{ z, 0.0, r },
qShellDot - _disc.nComp, qShellDot, dFluxDt, tlmAlloc);
Expand Down Expand Up @@ -799,7 +797,7 @@ void GeneralRateModelDG::consistentInitialTimeDerivative(const SimulationTime& s
*/
void GeneralRateModelDG::leanConsistentInitialState(const SimulationTime& simTime, double* const vecStateY, const AdJacobianParams& adJac, double errorTol, util::ThreadLocalStorage& threadLocalMem)
{
if (isSectionDependent(_parDiffusionMode) || isSectionDependent(_parSurfDiffusionMode))
if (isSectionDependent(_parDiffOp._parDiffusionMode) || isSectionDependent(_parDiffOp._parSurfDiffusionMode))
LOG(Warning) << "Lean consistent initialization is not appropriate for section-dependent pore and surface diffusion";

BENCH_SCOPE(_timerConsistentInit);
Expand Down Expand Up @@ -832,8 +830,7 @@ void GeneralRateModelDG::leanConsistentInitialState(const SimulationTime& simTim
const int localOffsetInParticle = static_cast<int>(shell) * idxr.strideParNode(type) + idxr.strideParLiquid();
double* const qShell = vecStateY + localOffsetToParticle + localOffsetInParticle;
// r (particle) coordinate of current node
const double r = static_cast<double>(_disc.deltaR[type]) * std::floor(shell / _disc.nParNode[type])
+ 0.5 * static_cast<double>(_disc.deltaR[type]) * (1 + _disc.parNodes[type][shell % _disc.nParNode[type]]);
const double r = _parDiffOp.relativeCoordinate(type, shell);
const ColumnPosition colPos{ z, 0.0, r};

// Perform consistent initialization that does not require a full fledged nonlinear solver (that may fail or damage the current state vector)
Expand All @@ -843,7 +840,6 @@ void GeneralRateModelDG::leanConsistentInitialState(const SimulationTime& simTim
} CADET_PARFOR_END;
}
}

}

/**
Expand Down Expand Up @@ -889,7 +885,7 @@ void GeneralRateModelDG::leanConsistentInitialState(const SimulationTime& simTim
*/
void GeneralRateModelDG::leanConsistentInitialTimeDerivative(double t, double const* const vecStateY, double* const vecStateYdot, double* const res, util::ThreadLocalStorage& threadLocalMem)
{
if (isSectionDependent(_parDiffusionMode) || isSectionDependent(_parSurfDiffusionMode))
if (isSectionDependent(_parDiffOp._parDiffusionMode) || isSectionDependent(_parDiffOp._parSurfDiffusionMode))
LOG(Warning) << "Lean consistent initialization is not appropriate for section-dependent pore and surface diffusion";

BENCH_SCOPE(_timerConsistentInit);
Expand Down Expand Up @@ -1272,7 +1268,7 @@ void GeneralRateModelDG::solveBulkTimeDerivativeSystem(const SimulationTime& sim
void GeneralRateModelDG::leanConsistentInitialSensitivity(const SimulationTime& simTime, const ConstSimulationState& simState,
std::vector<double*>& vecSensY, std::vector<double*>& vecSensYdot, active const* const adRes, util::ThreadLocalStorage& threadLocalMem)
{
if (isSectionDependent(_parDiffusionMode) || isSectionDependent(_parSurfDiffusionMode))
if (isSectionDependent(_parDiffOp._parDiffusionMode) || isSectionDependent(_parDiffOp._parSurfDiffusionMode))
LOG(Warning) << "Lean consistent initialization is not appropriate for section-dependent pore and surface diffusion";

BENCH_SCOPE(_timerConsistentInit);
Expand Down
2 changes: 1 addition & 1 deletion src/libcadet/model/GeneralRateModelDG-LinearSolver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -220,7 +220,7 @@ void GeneralRateModelDG::assembleDiscretizedGlobalJacobian(double alpha, Indexer
*/
void GeneralRateModelDG::addTimeDerivativeToJacobianParticleShell(linalg::BandedEigenSparseRowIterator& jac, const Indexer& idxr, double alpha, unsigned int parType)
{
parts::cell::addTimeDerivativeToJacobianParticleShell<linalg::BandedEigenSparseRowIterator, true>(jac, alpha, static_cast<double>(_parPorosity[parType]), _disc.nComp, _disc.nBound + _disc.nComp * parType,
parts::cell::addTimeDerivativeToJacobianParticleShell<linalg::BandedEigenSparseRowIterator, true>(jac, alpha, static_cast<double>(_parDiffOp._parPorosity[parType]), _disc.nComp, _disc.nBound + _disc.nComp * parType,
_poreAccessFactor.data() + _disc.nComp * parType, _disc.strideBound[parType], _disc.boundOffset + _disc.nComp * parType, _binding[parType]->reactionQuasiStationarity());
}

Expand Down
Loading
Loading