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 22, 2024
1 parent 40aba28 commit f8f20c0
Show file tree
Hide file tree
Showing 23 changed files with 874 additions and 266 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
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
21 changes: 9 additions & 12 deletions src/libcadet/model/ParameterDependence.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,7 @@ namespace cadet

class IParameterProvider;
struct ColumnPosition;
class IModel;

namespace model
{
Expand Down Expand Up @@ -454,59 +455,55 @@ class IParameterParameterDependence
* @brief Evaluates the parameter
* @details This function is called simultaneously from multiple threads.
*
* @param [in] parTypeIdx Index of the particle type this parameter dependence belongs to
* @param [in] params Parameters of the unit operation
* @param [in] model Model that owns this parameter dependence
* @param [in] colPos Position in normalized coordinates (column inlet = 0, column outlet = 1; outer shell = 1, inner center = 0)
* @param [in] comp Index of the component the parameter belongs to (or @c -1 if independent of components)
* @param [in] parType Index of the particle type the parameter belongs to (or @c -1 if independent of particle types)
* @param [in] bnd Index of the bound state the parameter belongs to (or @c -1 if independent of bound states)
* @return Actual parameter value
*/
virtual double getValue(UnitOpIdx unitOpIdx, const std::unordered_map<ParameterId, active*>& params, const ColumnPosition& colPos, int comp, int parType, int bnd) const = 0;
virtual double getValue(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd) const = 0;

/**
* @brief Evaluates the parameter
* @details This function is called simultaneously from multiple threads.
*
* @param [in] parTypeIdx Index of the particle type this parameter dependence belongs to
* @param [in] params Parameters of the unit operation
* @param [in] model Model that owns this parameter dependence
* @param [in] colPos Position in normalized coordinates (column inlet = 0, column outlet = 1; outer shell = 1, inner center = 0)
* @param [in] comp Index of the component the parameter belongs to (or @c -1 if independent of components)
* @param [in] parType Index of the particle type the parameter belongs to (or @c -1 if independent of particle types)
* @param [in] bnd Index of the bound state the parameter belongs to (or @c -1 if independent of bound states)
* @return Actual parameter value
*/
virtual active getValueActive(UnitOpIdx unitOpIdx, const std::unordered_map<ParameterId, active*>& params, const ColumnPosition& colPos, int comp, int parType, int bnd) const = 0;
virtual active getValueActive(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd) const = 0;

/**
* @brief Evaluates the parameter
* @details This function is called simultaneously from multiple threads.
*
* @param [in] parTypeIdx Index of the particle type this parameter dependence belongs to
* @param [in] params Parameters of the unit operation
* @param [in] model Model that owns this parameter dependence
* @param [in] colPos Position in normalized coordinates (column inlet = 0, column outlet = 1; outer shell = 1, inner center = 0)
* @param [in] comp Index of the component the parameter belongs to (or @c -1 if independent of components)
* @param [in] parType Index of the particle type the parameter belongs to (or @c -1 if independent of particle types)
* @param [in] bnd Index of the bound state the parameter belongs to (or @c -1 if independent of bound states)
* @param [in] val Additional parameter-dependent value
* @return Actual parameter value
*/
virtual double getValue(UnitOpIdx unitOpIdx, const std::unordered_map<ParameterId, active*>& params, const ColumnPosition& colPos, int comp, int parType, int bnd, double val) const = 0;
virtual double getValue(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd, double val) const = 0;

/**
* @brief Evaluates the parameter
* @details This function is called simultaneously from multiple threads.
*
* @param [in] parTypeIdx Index of the particle type this parameter dependence belongs to
* @param [in] params Parameters of the unit operation
* @param [in] model Model that owns this parameter dependence
* @param [in] colPos Position in normalized coordinates (column inlet = 0, column outlet = 1; outer shell = 1, inner center = 0)
* @param [in] comp Index of the component the parameter belongs to (or @c -1 if independent of components)
* @param [in] parType Index of the particle type the parameter belongs to (or @c -1 if independent of particle types)
* @param [in] bnd Index of the bound state the parameter belongs to (or @c -1 if independent of bound states)
* @param [in] val Additional parameter-dependent value
* @return Actual parameter value
*/
virtual active getValue(UnitOpIdx unitOpIdx, const std::unordered_map<ParameterId, active*>& params, const ColumnPosition& colPos, int comp, int parType, int bnd, const active& val) const = 0;
virtual active getValue(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd, const active& val) const = 0;

protected:
};
Expand Down
Loading

0 comments on commit f8f20c0

Please sign in to comment.