WIP add parameter dependencies to chromatography models

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 <jan.breuer@fz-juelich.de>

Author: sleweke-bayer

src/libcadet/model/GeneralRateModel.cpp

@@ -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;
 

src/libcadet/model/GeneralRateModel2D.cpp

@@ -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()];
@@ -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;
 

src/libcadet/model/GeneralRateModelDGFV.cpp

@@ -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");

src/libcadet/model/LumpedRateModelWithPores.cpp

@@ -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"
@@ -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)
 {
@@ -75,6 +76,7 @@ LumpedRateModelWithPores<ConvDispOperator>::~LumpedRateModelWithPores() CADET_NO
 	delete[] _tempState;
 
 	delete _dynReactionBulk;
+	delete _filmDiffDep;
 
 	delete[] _disc.parTypeOffset;
 	delete[] _disc.nBound;
@@ -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);
 
@@ -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);
@@ -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;
 
@@ -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
@@ -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];
 			}
 		}
 
@@ -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
@@ -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);
 			}
 		}
 
@@ -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);
@@ -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);
@@ -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);

src/libcadet/model/LumpedRateModelWithPores.hpp

@@ -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

src/libcadet/model/LumpedRateModelWithoutPores.cpp

@@ -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);
@@ -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());
 		}
 	};
 }
@@ -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);
 
@@ -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