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>

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

src/libcadet/model/ParameterDependence.hpp

@@ -33,6 +33,7 @@ namespace cadet
 
 class IParameterProvider;
 struct ColumnPosition;
+class IModel;
 
 namespace model
 {
@@ -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:
 };