Skip to content

Commit

Permalink
WIP add parameter dependencies to chromatography models
Browse files Browse the repository at this point in the history 
Support radial flow and param deps in createLWE
Add support for radial flow models and parameter dependence (i.e.,
dispersion and film diffusion coefficients) to createLWE

Add axial dispersion and film diffusion parameter dependence

Axial dispersion for all models, film diffusion parameter dependence for LRMP and GRM

Add tests for radial flow convection dispersion operators

Enable radial dispersion coeff dependency in operators

Co-authored-by: jbreue16 <[email protected]>
  • Loading branch information
@sleweke-bayer @jbreue16
sleweke-bayer and jbreue16 committed Apr 19, 2024
1 parent a5e8904 commit 5c5c55c
Show file tree
Hide file tree
Showing 24 changed files with 875 additions and 267 deletions.
2 changes: 1 addition & 1 deletion src/libcadet/model/GeneralRateModel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1263,7 +1263,7 @@ template <typename ConvDispOperator>
template <typename StateType, typename ResidualType, typename ParamType, bool wantJac>
int GeneralRateModel<ConvDispOperator>::residualBulk(double t, unsigned int secIdx, StateType const* yBase, double const* yDotBase, ResidualType* resBase, util::ThreadLocalStorage& threadLocalMem)
{
_convDispOp.residual(t, secIdx, yBase, yDotBase, resBase, wantJac, typename ParamSens<ParamType>::enabled());
_convDispOp.residual(*this, t, secIdx, yBase, yDotBase, resBase, wantJac, typename ParamSens<ParamType>::enabled());
if (!_dynReactionBulk || (_dynReactionBulk->numReactionsLiquid() == 0))
return 0;

Expand Down
4 changes: 2 additions & 2 deletions src/libcadet/model/GeneralRateModel2D.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -723,7 +723,7 @@ bool GeneralRateModel2D::configureModelDiscretization(IParameterProvider& paramP
}
}

const bool transportSuccess = _convDispOp.configureModelDiscretization(paramProvider, _disc.nComp, _disc.nCol, _disc.nRad, _dynReactionBulk);
const bool transportSuccess = _convDispOp.configureModelDiscretization(paramProvider, helper, _disc.nComp, _disc.nCol, _disc.nRad, _dynReactionBulk);

// Setup the memory for tempState based on state vector
_tempState = new double[numDofs()];
Expand Down Expand Up @@ -1336,7 +1336,7 @@ int GeneralRateModel2D::residualImpl(double t, unsigned int secIdx, StateType co
template <typename StateType, typename ResidualType, typename ParamType, bool wantJac>
int GeneralRateModel2D::residualBulk(double t, unsigned int secIdx, StateType const* yBase, double const* yDotBase, ResidualType* resBase, util::ThreadLocalStorage& threadLocalMem)
{
_convDispOp.residual(t, secIdx, yBase, yDotBase, resBase, wantJac, typename ParamSens<ParamType>::enabled());
_convDispOp.residual(*this, t, secIdx, yBase, yDotBase, resBase, wantJac, typename ParamSens<ParamType>::enabled());
if (!_dynReactionBulk || (_dynReactionBulk->numReactionsLiquid() == 0))
return 0;

Expand Down
2 changes: 1 addition & 1 deletion src/libcadet/model/GeneralRateModelDGFV.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -434,7 +434,7 @@ bool GeneralRateModelDGFV::configureModelDiscretization(IParameterProvider& para

_disc.dispersion = Eigen::VectorXd::Zero(_disc.nComp); // fill later on with convDispOpB (section and component dependent)

_disc.velocity = static_cast<double>(_convDispOp.currentVelocity()); // updated later on with convDispOpB (section dependent)
_disc.velocity = static_cast<double>(_convDispOp.currentVelocity(0)); // updated later on with convDispOpB (section dependent)
_disc.curSection = -1;

_disc.length_ = paramProvider.getDouble("COL_LENGTH");
Expand Down
81 changes: 75 additions & 6 deletions src/libcadet/model/LumpedRateModelWithPores.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@
#include "model/BindingModel.hpp"
#include "model/ReactionModel.hpp"
#include "model/parts/BindingCellKernel.hpp"
#include "model/ParameterDependence.hpp"
#include "SimulationTypes.hpp"
#include "linalg/DenseMatrix.hpp"
#include "linalg/BandMatrix.hpp"
Expand Down Expand Up @@ -63,7 +64,7 @@ int schurComplementMultiplierLRMPores(void* userData, double const* x, double* z

template <typename ConvDispOperator>
LumpedRateModelWithPores<ConvDispOperator>::LumpedRateModelWithPores(UnitOpIdx unitOpIdx) : UnitOperationBase(unitOpIdx),
_dynReactionBulk(nullptr), _jacP(0), _jacPdisc(0), _jacPF(0), _jacFP(0), _jacInlet(), _analyticJac(true),
_dynReactionBulk(nullptr), _filmDiffDep(nullptr), _jacP(0), _jacPdisc(0), _jacPF(0), _jacFP(0), _jacInlet(), _analyticJac(true),
_jacobianAdDirs(0), _factorizeJacobian(false), _tempState(nullptr), _initC(0), _initCp(0), _initQ(0),
_initState(0), _initStateDot(0)
{
Expand All @@ -75,6 +76,7 @@ LumpedRateModelWithPores<ConvDispOperator>::~LumpedRateModelWithPores() CADET_NO
delete[] _tempState;

delete _dynReactionBulk;
delete _filmDiffDep;

delete[] _disc.parTypeOffset;
delete[] _disc.nBound;
Expand Down Expand Up @@ -263,6 +265,18 @@ bool LumpedRateModelWithPores<ConvDispOperator>::configureModelDiscretization(IP

const bool transportSuccess = _convDispOp.configureModelDiscretization(paramProvider, helper, _disc.nComp, _disc.nCol);

if (paramProvider.exists("FILM_DIFFUSION_DEP"))
{
const std::string paramDepName = paramProvider.getString("FILM_DIFFUSION_DEP");
_filmDiffDep = helper.createParameterParameterDependence(paramDepName);
if (!_filmDiffDep)
throw InvalidParameterException("Unknown parameter dependence " + paramDepName + " in FILM_DIFFUSION_DEP");

_filmDiffDep->configureModelDiscretization(paramProvider);
}
else
_filmDiffDep = helper.createParameterParameterDependence("CONSTANT_ONE");

// Allocate memory
Indexer idxr(_disc);

Expand Down Expand Up @@ -402,6 +416,12 @@ bool LumpedRateModelWithPores<ConvDispOperator>::configure(IParameterProvider& p

const bool transportSuccess = _convDispOp.configure(_unitOpIdx, paramProvider, _parameters);

if (_filmDiffDep)
{
if (!_filmDiffDep->configure(paramProvider, _unitOpIdx, ParTypeIndep, BoundStateIndep, "FILM_DIFFUSION_DEP"))
throw InvalidParameterException("Failed to configure film diffusion parameter dependency (FILM_DIFFUSION_DEP)");
}

// Read geometry parameters
_colPorosity = paramProvider.getDouble("COL_POROSITY");
_singleParRadius = readAndRegisterMultiplexTypeParam(paramProvider, _parameters, _parRadius, "PAR_RADIUS", _disc.nParType, _unitOpIdx);
Expand Down Expand Up @@ -964,7 +984,7 @@ template <typename ConvDispOperator>
template <typename StateType, typename ResidualType, typename ParamType, bool wantJac>
int LumpedRateModelWithPores<ConvDispOperator>::residualBulk(double t, unsigned int secIdx, StateType const* yBase, double const* yDotBase, ResidualType* resBase, util::ThreadLocalStorage& threadLocalMem)
{
_convDispOp.residual(t, secIdx, yBase, yDotBase, resBase, wantJac, typename ParamSens<ParamType>::enabled());
_convDispOp.residual(*this, t, secIdx, yBase, yDotBase, resBase, wantJac, typename ParamSens<ParamType>::enabled());
if (!_dynReactionBulk || (_dynReactionBulk->numReactionsLiquid() == 0))
return 0;

Expand Down Expand Up @@ -1068,7 +1088,12 @@ int LumpedRateModelWithPores<ConvDispOperator>::residualFlux(double t, unsigned
{
const unsigned int colCell = i / _disc.nComp;
const unsigned int comp = i % _disc.nComp;
resCol[i] += jacCF_val * static_cast<ParamType>(filmDiff[comp]) * static_cast<ParamType>(_parTypeVolFrac[type + _disc.nParType * colCell]) * yFluxType[i];

const double relPos = _convDispOp.relativeCoordinate(colCell);
const active curVelocity = _convDispOp.currentVelocity(relPos);
const active modifier = _filmDiffDep->getValue(*this, ColumnPosition{relPos, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, curVelocity);

resCol[i] += jacCF_val * static_cast<ParamType>(filmDiff[comp]) * static_cast<ParamType>(modifier) * static_cast<ParamType>(_parTypeVolFrac[type + _disc.nParType * colCell]) * yFluxType[i];
}

// J_{f,0} block, adds bulk volume state c_i to flux equation
Expand All @@ -1084,10 +1109,15 @@ int LumpedRateModelWithPores<ConvDispOperator>::residualFlux(double t, unsigned
// J_{p,f} block, adds flux to particle / bead volume equations
for (unsigned int pblk = 0; pblk < _disc.nCol; ++pblk)
{
const double relPos = _convDispOp.relativeCoordinate(pblk);
const active curVelocity = _convDispOp.currentVelocity(relPos);

for (unsigned int comp = 0; comp < _disc.nComp; ++comp)
{
const active modifier = _filmDiffDep->getValue(*this, ColumnPosition{relPos, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, curVelocity);

const unsigned int eq = pblk * idxr.strideColCell() + comp * idxr.strideColComp();
resParType[pblk * idxr.strideParBlock(type) + comp] += jacPF_val / static_cast<ParamType>(poreAccFactor[comp]) * static_cast<ParamType>(filmDiff[comp]) * yFluxType[eq];
resParType[pblk * idxr.strideParBlock(type) + comp] += jacPF_val / static_cast<ParamType>(poreAccFactor[comp]) * static_cast<ParamType>(filmDiff[comp]) * static_cast<ParamType>(modifier) * yFluxType[eq];
}
}

Expand Down Expand Up @@ -1148,8 +1178,12 @@ void LumpedRateModelWithPores<ConvDispOperator>::assembleOffdiagJac(double t, un
const unsigned int colCell = eq / _disc.nComp;
const unsigned int comp = eq % _disc.nComp;

const double relPos = _convDispOp.relativeCoordinate(colCell);
const double curVelocity = static_cast<double>(_convDispOp.currentVelocity(relPos));
const double modifier = _filmDiffDep->getValue(*this, ColumnPosition{relPos, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, curVelocity);

// Main diagonal corresponds to j_{f,i} (flux) state variable
_jacCF.addElement(eq, eq + typeOffset, jacCF_val * static_cast<double>(filmDiff[comp]) * static_cast<double>(_parTypeVolFrac[type + _disc.nParType * colCell]));
_jacCF.addElement(eq, eq + typeOffset, jacCF_val * static_cast<double>(filmDiff[comp]) * modifier * static_cast<double>(_parTypeVolFrac[type + _disc.nParType * colCell]));
}

// J_{f,0} block, adds bulk volume state c_i to flux equation
Expand All @@ -1161,11 +1195,16 @@ void LumpedRateModelWithPores<ConvDispOperator>::assembleOffdiagJac(double t, un
// J_{p,f} block, implements bead boundary condition in outer bead shell equation
for (unsigned int pblk = 0; pblk < _disc.nCol; ++pblk)
{
const double relPos = _convDispOp.relativeCoordinate(pblk);
const double curVelocity = static_cast<double>(_convDispOp.currentVelocity(relPos));

for (unsigned int comp = 0; comp < _disc.nComp; ++comp)
{
const unsigned int eq = typeOffset + pblk * idxr.strideColCell() + comp * idxr.strideColComp();
const unsigned int col = pblk * idxr.strideParBlock(type) + comp;
_jacPF[type].addElement(col, eq, jacPF_val / static_cast<double>(poreAccFactor[comp]) * static_cast<double>(filmDiff[comp]));

const double modifier = _filmDiffDep->getValue(*this, ColumnPosition{relPos, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, curVelocity);
_jacPF[type].addElement(col, eq, jacPF_val / static_cast<double>(poreAccFactor[comp]) * static_cast<double>(filmDiff[comp]) * modifier);
}
}

Expand Down Expand Up @@ -1468,6 +1507,15 @@ bool LumpedRateModelWithPores<ConvDispOperator>::setParameter(const ParameterId&

if (_convDispOp.setParameter(pId, value))
return true;

if (_filmDiffDep)
{
if (_filmDiffDep->hasParameter(pId))
{
_filmDiffDep->setParameter(pId, value);
return true;
}
}
}

return UnitOperationBase::setParameter(pId, value);
Expand Down Expand Up @@ -1509,6 +1557,16 @@ void LumpedRateModelWithPores<ConvDispOperator>::setSensitiveParameterValue(cons

if (_convDispOp.setSensitiveParameterValue(_sensParams, pId, value))
return;

if (_filmDiffDep)
{
active* const param = _filmDiffDep->getParameter(pId);
if (param)
{
param->setValue(value);
return;
}
}
}

UnitOperationBase::setSensitiveParameterValue(pId, value);
Expand Down Expand Up @@ -1575,6 +1633,17 @@ bool LumpedRateModelWithPores<ConvDispOperator>::setSensitiveParameter(const Par
LOG(Debug) << "Found parameter " << pId << ": Dir " << adDirection << " is set to " << adValue;
return true;
}

if (_filmDiffDep)
{
active* const param = _filmDiffDep->getParameter(pId);
if (param)
{
LOG(Debug) << "Found parameter " << pId << ": Dir " << adDirection << " is set to " << adValue;
param->setADValue(adDirection, adValue);
return true;
}
}
}

return UnitOperationBase::setSensitiveParameter(pId, adDirection, adValue);
Expand Down
1 change: 1 addition & 0 deletions src/libcadet/model/LumpedRateModelWithPores.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -264,6 +264,7 @@ class LumpedRateModelWithPores : public UnitOperationBase

ConvDispOperator _convDispOp; //!< Convection dispersion operator for interstitial volume transport
IDynamicReactionModel* _dynReactionBulk; //!< Dynamic reactions in the bulk volume
IParameterParameterDependence* _filmDiffDep; //!< Film diffusion dependency on local velocity

std::vector<linalg::BandMatrix> _jacP; //!< Particle jacobian diagonal blocks (all of them for each particle type)
std::vector<linalg::FactorizableBandMatrix> _jacPdisc; //!< Particle jacobian diagonal blocks (all of them for each particle type) with time derivatives from BDF method
Expand Down
18 changes: 9 additions & 9 deletions src/libcadet/model/LumpedRateModelWithoutPores.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@ namespace
class ConvOpResidual
{
public:
static inline void call(ConvDispOperator& op, double t, unsigned int secIdx, StateType const* const y, double const* const yDot, ResidualType* const res, cadet::linalg::BandMatrix& jac)
static inline void call(cadet::IModel* const model, ConvDispOperator& op, double t, unsigned int secIdx, StateType const* const y, double const* const yDot, ResidualType* const res, cadet::linalg::BandMatrix& jac)
{
// This should not be reached
cadet_assert(false);
Expand All @@ -57,29 +57,29 @@ namespace
class ConvOpResidual<ConvDispOperator, double, ResidualType, ParamType, true>
{
public:
static inline void call(ConvDispOperator& op, double t, unsigned int secIdx, double const* const y, double const* const yDot, ResidualType* const res, cadet::linalg::BandMatrix& jac)
static inline void call(cadet::IModel* const model, ConvDispOperator& op, double t, unsigned int secIdx, double const* const y, double const* const yDot, ResidualType* const res, cadet::linalg::BandMatrix& jac)
{
op.residual(t, secIdx, y, yDot, res, jac);
op.residual(*model, t, secIdx, y, yDot, res, jac);
}
};

template <typename ConvDispOperator, typename ResidualType, typename ParamType>
class ConvOpResidual<ConvDispOperator, double, ResidualType, ParamType, false>
{
public:
static inline void call(ConvDispOperator& op, double t, unsigned int secIdx, double const* const y, double const* const yDot, ResidualType* const res, cadet::linalg::BandMatrix& jac)
static inline void call(cadet::IModel* const model, ConvDispOperator& op, double t, unsigned int secIdx, double const* const y, double const* const yDot, ResidualType* const res, cadet::linalg::BandMatrix& jac)
{
op.residual(t, secIdx, y, yDot, res, typename cadet::ParamSens<ParamType>::enabled());
op.residual(*model, t, secIdx, y, yDot, res, typename cadet::ParamSens<ParamType>::enabled());
}
};

template <typename ConvDispOperator, typename ResidualType, typename ParamType>
class ConvOpResidual<ConvDispOperator, cadet::active, ResidualType, ParamType, false>
{
public:
static inline void call(ConvDispOperator& op, double t, unsigned int secIdx, cadet::active const* const y, double const* const yDot, ResidualType* const res, cadet::linalg::BandMatrix& jac)
static inline void call(cadet::IModel* const model, ConvDispOperator& op, double t, unsigned int secIdx, cadet::active const* const y, double const* const yDot, ResidualType* const res, cadet::linalg::BandMatrix& jac)
{
op.residual(t, secIdx, y, yDot, res, typename cadet::ParamSens<ParamType>::enabled());
op.residual(*model, t, secIdx, y, yDot, res, typename cadet::ParamSens<ParamType>::enabled());
}
};
}
Expand Down Expand Up @@ -609,7 +609,7 @@ template <typename ConvDispOperator>
template <typename StateType, typename ResidualType, typename ParamType, bool wantJac>
int LumpedRateModelWithoutPores<ConvDispOperator>::residualImpl(double t, unsigned int secIdx, StateType const* const y, double const* const yDot, ResidualType* const res, util::ThreadLocalStorage& threadLocalMem)
{
ConvOpResidual<ConvDispOperator, StateType, ResidualType, ParamType, wantJac>::call(_convDispOp, t, secIdx, y, yDot, res, _jac);
ConvOpResidual<ConvDispOperator, StateType, ResidualType, ParamType, wantJac>::call(this, _convDispOp, t, secIdx, y, yDot, res, _jac);

Indexer idxr(_disc);

Expand Down Expand Up @@ -785,7 +785,7 @@ template <typename ConvDispOperator>
unsigned int LumpedRateModelWithoutPores<ConvDispOperator>::localOutletComponentIndex(unsigned int port) const CADET_NOEXCEPT
{
// Inlets are duplicated so need to be accounted for
if (static_cast<double>(_convDispOp.currentVelocity()) >= 0.0)
if (_convDispOp.forwardFlow())
// Forward Flow: outlet is last cell
return _disc.nComp + (_disc.nCol - 1) * (_disc.nComp + _disc.strideBound);
else
Expand Down
Loading

0 comments on commit 5c5c55c

Please sign in to comment.