Skip to content

Commit

Permalink
Fix bugs in DG units
Browse files Browse the repository at this point in the history 
* LRM inlet DOF treatment

* LRMP flux Jacobian analytical and AD fix

* GRM allocation of deltaR

* Clean up

* Fix bug in LRMP DG sensitivity init
  • Loading branch information
@jbreue16
jbreue16 committed Jun 18, 2024
1 parent a11c07a commit cbda1a2
Show file tree
Hide file tree
Showing 6 changed files with 99 additions and 112 deletions.
27 changes: 12 additions & 15 deletions src/libcadet/model/GeneralRateModelDG.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -210,7 +210,7 @@ bool GeneralRateModelDG::configureModelDiscretization(IParameterProvider& paramP

std::vector<int> parPolyDeg(_disc.nParType);
std::vector<int> ParNelem(_disc.nParType);
std::vector<bool> parExactInt(_disc.nParType);
std::vector<bool> parExactInt(_disc.nParType, true);
_disc.parPolyDeg = new unsigned int[_disc.nParType];
_disc.nParCell = new unsigned int[_disc.nParType];
_disc.parExactInt = new bool[_disc.nParType];
Expand All @@ -219,12 +219,17 @@ bool GeneralRateModelDG::configureModelDiscretization(IParameterProvider& paramP
if (paramProvider.exists("PAR_POLYDEG"))
{
parPolyDeg = paramProvider.getIntArray("PAR_POLYDEG");

if ((std::any_of(parPolyDeg.begin(), parPolyDeg.end(), [](int value) { return value < 1; })))
throw InvalidParameterException("Particle polynomial degrees must be at least 1!");
ParNelem = paramProvider.getIntArray("PAR_NELEM");

if ((std::any_of(ParNelem.begin(), ParNelem.end(), [](int value) { return value < 1; })))
throw InvalidParameterException("Particle number of elements must be at least 1!");
parExactInt = paramProvider.getBoolArray("PAR_EXACT_INTEGRATION");

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!";

Expand Down Expand Up @@ -766,6 +771,7 @@ bool GeneralRateModelDG::configure(IParameterProvider& paramProvider)
registerParam2DArray(_parameters, _parTypeVolFrac, [=](bool multi, unsigned cell, unsigned int type) { return makeParamId(hashString("PAR_TYPE_VOLFRAC"), _unitOpIdx, CompIndep, type, BoundStateIndep, ReactionIndep, cell); }, _disc.nParType);

// Calculate the particle radial discretization variables (_parCellSize, _parCenterRadius, etc.)
_disc.deltaR = new active[_disc.offsetMetric[_disc.nParType]];
updateRadialDisc();

// Register initial conditions parameters
Expand Down Expand Up @@ -954,18 +960,10 @@ void GeneralRateModelDG::notifyDiscontinuousSectionTransition(double t, unsigned
setJacobianPattern_GRM(_globalJac, 0, _dynReactionBulk);
_globalJacDisc = _globalJac;

// ConvectionDispersionOperator tells us whether flow direction has changed
if (!_convDispOp.notifyDiscontinuousSectionTransition(t, secIdx, _jacInlet)) {
// (re)compute DG particle Jacobian blocks (can only be done after notify)
updateSection(secIdx);
_disc.initializeDGjac(_parGeomSurfToVol);
return;
}
else {
// (re)compute DG particle Jacobian blocks
updateSection(secIdx);
_disc.initializeDGjac(_parGeomSurfToVol);
}
_convDispOp.notifyDiscontinuousSectionTransition(t, secIdx, _jacInlet);

updateSection(secIdx);
_disc.initializeDGjac(_parGeomSurfToVol);
}

void GeneralRateModelDG::setFlowRates(active const* in, active const* out) CADET_NOEXCEPT
Expand Down Expand Up @@ -2087,7 +2085,6 @@ void GeneralRateModelDG::setUserdefinedRadialDisc(unsigned int parType)
// This approach should only be used when necessary, i.e. solely when particle radius parameter sensitivity is required.
void GeneralRateModelDG::updateRadialDisc()
{
_disc.deltaR = new active[_disc.offsetMetric[_disc.nParType]];

for (unsigned int parType = 0; parType < _disc.nParType; ++parType)
{
Expand Down
116 changes: 58 additions & 58 deletions src/libcadet/model/GeneralRateModelDG.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -997,64 +997,6 @@ class GeneralRateModelDG : public UnitOperationBase
// ======================================== DG particle Jacobian ============================================= //
// ========================================================================================================================================================== //

/**
* @brief computes the jacobian via finite differences (testing purpose)
*/
MatrixXd calcFDJacobian(const double* y_, const double* yDot_, const SimulationTime simTime, util::ThreadLocalStorage& threadLocalMem, double alpha) {

// create solution vectors
Eigen::Map<const VectorXd> hmpf(y_, numDofs());
VectorXd y = hmpf;
VectorXd yDot;
if (yDot_) {
Eigen::Map<const VectorXd> hmpf2(yDot_, numDofs());
yDot = hmpf2;
}
else {
return MatrixXd::Zero(numDofs(), numDofs());
}
VectorXd res = VectorXd::Zero(numDofs());
const double* yPtr = &y[0];
const double* yDotPtr = &yDot[0];
double* resPtr = &res[0];
// create FD jacobian
MatrixXd Jacobian = MatrixXd::Zero(numDofs(), numDofs());
// set FD step
double epsilon = 0.01;

residualImpl<double, double, double, false>(simTime.t, simTime.secIdx, yPtr, yDotPtr, resPtr, threadLocalMem);

for (int col = 0; col < Jacobian.cols(); col++) {
Jacobian.col(col) = -(1.0 + alpha) * res;
}
/* Residual(y+h) */
// state DOFs
for (int dof = 0; dof < Jacobian.cols(); dof++) {
y[dof] += epsilon;
residualImpl<double, double, double, false>(simTime.t, simTime.secIdx, yPtr, yDotPtr, resPtr, threadLocalMem);
y[dof] -= epsilon;
Jacobian.col(dof) += res;
}

// state derivative Jacobian
for (int dof = 0; dof < Jacobian.cols(); dof++) {
yDot[dof] += epsilon;
residualImpl<double, double, double, false>(simTime.t, simTime.secIdx, yPtr, yDotPtr, resPtr, threadLocalMem);
yDot[dof] -= epsilon;
Jacobian.col(dof) += alpha * res;
}

///* exterminate numerical noise */
//for (int i = 0; i < Jacobian.rows(); i++) {
// for (int j = 0; j < Jacobian.cols(); j++) {
// if (std::abs(Jacobian(i, j)) < 1e-10) Jacobian(i, j) = 0.0;
// }
//}
Jacobian /= epsilon;

return Jacobian;
}

typedef Eigen::Triplet<double> T;
/**
* @brief calculates the particle dispersion jacobian Pattern of the exact/inexact integration DG scheme for the given particle type and bead
Expand Down Expand Up @@ -2594,6 +2536,64 @@ class GeneralRateModelDG : public UnitOperationBase
return _globalJac.isCompressed(); // check if the jacobian estimation fits the pattern
}

/**
* @brief computes the jacobian via finite differences (testing purpose)
*/
MatrixXd calcFDJacobian(const double* y_, const double* yDot_, const SimulationTime simTime, util::ThreadLocalStorage& threadLocalMem, double alpha) {

// create solution vectors
Eigen::Map<const VectorXd> hmpf(y_, numDofs());
VectorXd y = hmpf;
VectorXd yDot;
if (yDot_) {
Eigen::Map<const VectorXd> hmpf2(yDot_, numDofs());
yDot = hmpf2;
}
else {
return MatrixXd::Zero(numDofs(), numDofs());
}
VectorXd res = VectorXd::Zero(numDofs());
const double* yPtr = &y[0];
const double* yDotPtr = &yDot[0];
double* resPtr = &res[0];
// create FD jacobian
MatrixXd Jacobian = MatrixXd::Zero(numDofs(), numDofs());
// set FD step
double epsilon = 0.01;

residualImpl<double, double, double, false>(simTime.t, simTime.secIdx, yPtr, yDotPtr, resPtr, threadLocalMem);

for (int col = 0; col < Jacobian.cols(); col++) {
Jacobian.col(col) = -(1.0 + alpha) * res;
}
/* Residual(y+h) */
// state DOFs
for (int dof = 0; dof < Jacobian.cols(); dof++) {
y[dof] += epsilon;
residualImpl<double, double, double, false>(simTime.t, simTime.secIdx, yPtr, yDotPtr, resPtr, threadLocalMem);
y[dof] -= epsilon;
Jacobian.col(dof) += res;
}

// state derivative Jacobian
for (int dof = 0; dof < Jacobian.cols(); dof++) {
yDot[dof] += epsilon;
residualImpl<double, double, double, false>(simTime.t, simTime.secIdx, yPtr, yDotPtr, resPtr, threadLocalMem);
yDot[dof] -= epsilon;
Jacobian.col(dof) += alpha * res;
}

///* exterminate numerical noise */
//for (int i = 0; i < Jacobian.rows(); i++) {
// for (int j = 0; j < Jacobian.cols(); j++) {
// if (std::abs(Jacobian(i, j)) < 1e-10) Jacobian(i, j) = 0.0;
// }
//}
Jacobian /= epsilon;

return Jacobian;
}

};

} // namespace model
Expand Down
18 changes: 7 additions & 11 deletions src/libcadet/model/LumpedRateModelWithPoresDG-InitialConditions.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -343,7 +343,7 @@ namespace cadet
LinearBufferAllocator tlmAlloc = threadLocalMem.get();

// Reuse memory of band matrix for dense matrix
linalg::DenseMatrixView fullJacobianMatrix(_globalJacDisc.valuePtr() + _globalJacDisc.outerIndexPtr()[idxr.offsetCp(ParticleTypeIndex{ type }) - idxr.offsetC() + pblk], nullptr, mask.len, mask.len);
linalg::DenseMatrixView fullJacobianMatrix(_globalJacDisc.valuePtr() + _globalJacDisc.outerIndexPtr()[idxr.offsetCp(ParticleTypeIndex{ type }, ParticleIndex{ pblk }) - idxr.offsetC()], nullptr, mask.len, mask.len);

// z coordinate (column length normed to 1) of current node - needed in externally dependent adsorption kinetic
const double z = _convDispOp.relativeCoordinate(pblk);
Expand Down Expand Up @@ -402,7 +402,7 @@ namespace cadet
linalg::conservedMoietiesFromPartitionedMask(mask, _disc.nBound + type * _disc.nComp, _disc.nComp, qShell - _disc.nComp, conservedQuants, static_cast<double>(_parPorosity[type]), epsQ);

std::function<bool(double const* const, linalg::detail::DenseMatrixBase&)> jacFunc;
if (localAdY && localAdRes)
if (localAdY && localAdRes) // todo fix AD consistent initialization with req. binding
{
jacFunc = [&](double const* const x, linalg::detail::DenseMatrixBase& mat)
{
Expand All @@ -420,7 +420,7 @@ namespace cadet
simTime.t, simTime.secIdx, colPos, localAdY, nullptr, localAdRes, fullJacobianMatrix.row(0), cellResParams, tlmAlloc
);

// todo ? check analytical jacobian?
// // todo check analytical jacobian
//#ifdef CADET_CHECK_ANALYTIC_JACOBIAN
// std::copy_n(qShell - _disc.nComp, mask.len, fullX);
// linalg::applyVectorSubset(x, mask, fullX);
Expand All @@ -447,7 +447,7 @@ namespace cadet
for (unsigned int par = 0; par < _disc.nPoints; par++)
{
const int eqOffset_res = idxr.offsetCp(ParticleTypeIndex{ type }, ParticleIndex{ par });
const int eqOffset_mat = idxr.offsetCp(ParticleTypeIndex{ type }, ParticleIndex{ par });
const int eqOffset_mat = idxr.offsetCp(ParticleTypeIndex{ type }, ParticleIndex{ par }) - idxr.offsetC();
for (unsigned int phase = 0; phase < idxr.strideParBlock(type); phase++)
{
for (unsigned int phaseTo = 0; phaseTo < idxr.strideParBlock(type); phaseTo++)
Expand Down Expand Up @@ -912,7 +912,7 @@ namespace cadet
{
// Loop over components in cell
for (unsigned comp = 0; comp < _disc.nComp; ++comp)
stateYbulk[point * idxr.strideColCell() + comp * idxr.strideColComp()] = _initC[comp].getADValue(param);
stateYbulk[point * idxr.strideColNode() + comp * idxr.strideColComp()] = _initC[comp].getADValue(param);
}

// Loop over particles
Expand Down Expand Up @@ -1003,9 +1003,7 @@ namespace cadet
for (unsigned int i = _disc.nComp; i < numDofs(); ++i)
sensYdot[i] = -adRes[i].getADValue(param);

// Step 1: Solve algebraic equations

// Step 1a: Compute quasi-stationary binding model state
// Step 1: Compute quasi-stationary binding model state
for (unsigned int type = 0; type < _disc.nParType; ++type)
{
if (!_binding[type]->hasQuasiStationaryReactions())
Expand Down Expand Up @@ -1073,9 +1071,7 @@ namespace cadet
} CADET_PARFOR_END;
}

// Step 1: Compute the correct time derivative of the state vector

// Step 1a: Assemble, factorize, and solve diagonal blocks of linear system
// Step 2: Compute the correct time derivative of the state vector: Assemble, factorize, and solve diagonal blocks of linear system

// Compute right hand side by adding -dF / dy * s = -J * s to -dF / dp which is already stored in sensYdot
multiplyWithJacobian(simTime, simState, sensY, -1.0, 1.0, sensYdot);
Expand Down
19 changes: 8 additions & 11 deletions src/libcadet/model/LumpedRateModelWithPoresDG.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -181,7 +181,7 @@ bool LumpedRateModelWithPoresDG::configureModelDiscretization(IParameterProvider
_disc.nPoints = _disc.nNodes * _disc.nCol;

int polynomial_integration_mode = 0;
if(paramProvider.exists("EXACT_INTEGRATION"))
if (paramProvider.exists("EXACT_INTEGRATION"))
polynomial_integration_mode = paramProvider.getInt("EXACT_INTEGRATION");
_disc.exactInt = static_cast<bool>(polynomial_integration_mode); // only integration mode 0 applies the inexact collocated diagonal LGL mass matrix

Expand Down Expand Up @@ -667,11 +667,6 @@ void LumpedRateModelWithPoresDG::useAnalyticJacobian(const bool analyticJac)

void LumpedRateModelWithPoresDG::notifyDiscontinuousSectionTransition(double t, unsigned int secIdx, const ConstSimulationState& simState, const AdJacobianParams& adJac)
{
// Setup flux Jacobian blocks at the beginning of the simulation or in case of
// section dependent film or particle diffusion coefficients
if ((secIdx == 0) || isSectionDependent(_filmDiffusionMode))
assembleFluxJacobian(t, secIdx);

updateSection(secIdx);

_convDispOp.notifyDiscontinuousSectionTransition(t, secIdx, _jacInlet);
Expand Down Expand Up @@ -774,8 +769,6 @@ void LumpedRateModelWithPoresDG::extractJacobianFromAD(active const* const adRes
int eqOffset = 0;
ad::extractBandedBlockEigenJacobianFromAd(adVec, adDirOffset, diagDir, lowerBandwidth, upperBandwidth, eqOffset, bulkDoFs, _globalJac);

/* Film diffusion flux entries are handled analytically (todo: point to where that happens) */

/* Handle particle liquid and solid phase equations entries */
// Read particle Jacobian enries from dedicated AD directions
int offsetParticleTypeDirs = adDirOffset + _convDispOp.requiredADdirs();
Expand All @@ -796,6 +789,10 @@ void LumpedRateModelWithPoresDG::extractJacobianFromAD(active const* const adRes
}
offsetParticleTypeDirs += idxr.strideParBlock(type);
}

/* Film diffusion flux entries are handled analytically (only cross dependent entries) */
calcFluxJacobians(_disc.curSection, true);

}

#ifdef CADET_CHECK_ANALYTIC_JACOBIAN
Expand Down Expand Up @@ -955,11 +952,11 @@ int LumpedRateModelWithPoresDG::residualImpl(double t, unsigned int secIdx, Stat
// check for section switch
updateSection(secIdx);
bool success = 1;

if (wantJac)
{
if (_disc.newStaticJac) { // static (per section) transport Jacobian

if (_disc.newStaticJac) // static (per section) transport Jacobian
{
success = calcStaticAnaGlobalJacobian(secIdx);
_disc.newStaticJac = false;
if (cadet_unlikely(!success))
Expand Down
14 changes: 7 additions & 7 deletions src/libcadet/model/LumpedRateModelWithPoresDG.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -528,15 +528,15 @@ class LumpedRateModelWithPoresDG : public UnitOperationBase
* @brief analytically calculates the static (per section) bulk jacobian (inlet DOFs included!)
* @return 1 if jacobain estimation fits the predefined pattern of the jacobian, 0 if not.
*/
int calcStaticAnaGlobalJacobian(unsigned int secIdx) {
int calcStaticAnaGlobalJacobian(const unsigned int secIdx) {

bool success = _convDispOp.calcStaticAnaJacobian(_globalJac, _jacInlet);
success = success && calcFluxJacobians(secIdx);

return success;
}

int calcFluxJacobians(unsigned int secIdx) {
int calcFluxJacobians(const unsigned int secIdx, const bool crossDepsOnly = false) {

Indexer idxr(_disc);

Expand Down Expand Up @@ -565,7 +565,8 @@ class LumpedRateModelWithPoresDG : public UnitOperationBase
// add Cl on Cl entries (added since already set in bulk jacobian)
// row: already at bulk phase. already at current node and component.
// col: already at bulk phase. already at current node and component.
jacC[0] += jacCF_val * static_cast<double>(filmDiff[comp]) * static_cast<double>(_parTypeVolFrac[type + _disc.nParType * colNode]);
if(!crossDepsOnly)
jacC[0] += jacCF_val * static_cast<double>(filmDiff[comp]) * static_cast<double>(_parTypeVolFrac[type + _disc.nParType * colNode]);
// add Cl on Cp entries
// row: already at bulk phase. already at current node and component.
// col: already at bulk phase. already at current node and component.
Expand All @@ -574,13 +575,12 @@ class LumpedRateModelWithPoresDG : public UnitOperationBase
// add Cp on Cp entries
// row: already at particle. already at current node and liquid state.
// col: go to flux of current parType and adjust for offsetC. jump over previous colNodes and add component offset
jacP[0]
= -jacPF_val / static_cast<double>(poreAccFactor[comp]) * static_cast<double>(filmDiff[comp]);
if (!crossDepsOnly)
jacP[0] = -jacPF_val / static_cast<double>(poreAccFactor[comp]) * static_cast<double>(filmDiff[comp]);
// add Cp on Cl entries
// row: already at particle. already at current node and liquid state.
// col: go to flux of current parType and adjust for offsetC. jump over previous colNodes and add component offset
jacP[jacC.row() - jacP.row()]
= jacPF_val / static_cast<double>(poreAccFactor[comp]) * static_cast<double>(filmDiff[comp]);
jacP[jacC.row() - jacP.row()] = jacPF_val / static_cast<double>(poreAccFactor[comp]) * static_cast<double>(filmDiff[comp]);
}
}
}
Expand Down
Loading

0 comments on commit cbda1a2

Please sign in to comment.