Add axial dispersion and film diffusion parameter dependence

Axial dispersion for all models, film diffusion parameter dependence
only for LRMP.

Author: sleweke-bayer

src/libcadet/model/GeneralRateModel2D.cpp

@@ -703,7 +703,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()];
@@ -1316,7 +1316,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;
@@ -237,6 +239,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);
 
@@ -376,6 +390,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);
@@ -1042,7 +1062,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
@@ -1058,10 +1083,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];
 			}
 		}
 
@@ -1122,8 +1152,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
@@ -1135,11 +1169,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);
 			}
 		}
 
@@ -1442,6 +1481,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);
@@ -1483,6 +1531,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);
@@ -1549,6 +1607,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

@@ -775,7 +775,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/paramdep/DummyParameterDependence.cpp

@@ -27,17 +27,17 @@ namespace model
 {
 
 /**
- * @brief Defines a dummy parameter dependence
+ * @brief Defines a parameter dependence that outputs constant 0.0
  */
-class DummyParameterStateDependence : public ParameterStateDependenceBase
+class ConstantZeroParameterStateDependence : public ParameterStateDependenceBase
 {
 public:
 
-	DummyParameterStateDependence() { }
-	virtual ~DummyParameterStateDependence() CADET_NOEXCEPT { }
+	ConstantZeroParameterStateDependence() { }
+	virtual ~ConstantZeroParameterStateDependence() CADET_NOEXCEPT { }
 
-	static const char* identifier() { return "DUMMY"; }
-	virtual const char* name() const CADET_NOEXCEPT { return DummyParameterStateDependence::identifier(); }
+	static const char* identifier() { return "CONSTANT_ZERO"; }
+	virtual const char* name() const CADET_NOEXCEPT { return ConstantZeroParameterStateDependence::identifier(); }
 
 	virtual int jacobianElementsPerRowLiquid() const CADET_NOEXCEPT { return 0; }
 	virtual int jacobianElementsPerRowCombined() const CADET_NOEXCEPT { return 0; }
@@ -85,17 +85,75 @@ class DummyParameterStateDependence : public ParameterStateDependenceBase
 
 
 /**
- * @brief Defines a dummy parameter dependence
+ * @brief Defines a parameter dependence that outputs constant 1.0
  */
-class DummyParameterParameterDependence : public ParameterParameterDependenceBase
+class ConstantOneParameterStateDependence : public ParameterStateDependenceBase
 {
 public:
 
-	DummyParameterParameterDependence() { }
-	virtual ~DummyParameterParameterDependence() CADET_NOEXCEPT { }
+	ConstantOneParameterStateDependence() { }
+	virtual ~ConstantOneParameterStateDependence() CADET_NOEXCEPT { }
 
-	static const char* identifier() { return "DUMMY"; }
-	virtual const char* name() const CADET_NOEXCEPT { return DummyParameterParameterDependence::identifier(); }
+	static const char* identifier() { return "CONSTANT_ONE"; }
+	virtual const char* name() const CADET_NOEXCEPT { return ConstantOneParameterStateDependence::identifier(); }
+
+	virtual int jacobianElementsPerRowLiquid() const CADET_NOEXCEPT { return 0; }
+	virtual int jacobianElementsPerRowCombined() const CADET_NOEXCEPT { return 0; }
+
+	virtual void analyticJacobianLiquidAdd(const ColumnPosition& colPos, double param, double const* y, int comp, double factor, int offset, int row, linalg::DoubleSparseMatrix& jac) const { }
+	virtual void analyticJacobianCombinedAddLiquid(const ColumnPosition& colPos, double param, double const* yLiquid, double const* ySolid, int comp, double factor, int offset, int row, linalg::DoubleSparseMatrix& jac) const { }
+	virtual void analyticJacobianCombinedAddSolid(const ColumnPosition& colPos, double param, double const* yLiquid, double const* ySolid, int bnd, double factor, int offset, int row, linalg::DoubleSparseMatrix& jac) const { }
+
+	CADET_PARAMETERSTATEDEPENDENCE_BOILERPLATE
+
+protected:
+
+	virtual bool configureImpl(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, const std::string& name)
+	{
+		return true;
+	}
+
+	template <typename StateType, typename ParamType>
+	typename DoubleActivePromoter<StateType, ParamType>::type liquidParameterImpl(const ColumnPosition& colPos, const ParamType& param, StateType const* y, int comp) const
+	{
+		return 1.0;
+	}
+
+	template <typename RowIterator>
+	void analyticJacobianLiquidAddImpl(const ColumnPosition& colPos, double param, double const* y, int comp, double factor, int offset, RowIterator jac) const { }
+
+	template <typename StateType, typename ParamType>
+	typename DoubleActivePromoter<StateType, ParamType>::type combinedParameterLiquidImpl(const ColumnPosition& colPos, const ParamType& param, StateType const* yLiquid, StateType const* ySolid, int comp) const
+	{
+		return 1.0;
+	}
+
+	template <typename RowIterator>
+	void analyticJacobianCombinedAddLiquidImpl(const ColumnPosition& colPos, double param, double const* yLiquid, double const* ySolid, int comp, double factor, int offset, RowIterator jac) const { }
+
+	template <typename StateType, typename ParamType>
+	typename DoubleActivePromoter<StateType, ParamType>::type combinedParameterSolidImpl(const ColumnPosition& colPos, const ParamType& param, StateType const* yLiquid, StateType const* ySolid, int bnd) const
+	{
+		return 1.0;
+	}
+
+	template <typename RowIterator>
+	void analyticJacobianCombinedAddSolidImpl(const ColumnPosition& colPos, double param, double const* yLiquid, double const* ySolid, int bnd, double factor, int offset, RowIterator jac) const { }
+};
+
+
+/**
+ * @brief Defines a parameter dependence that outputs constant 0.0
+ */
+class ConstantZeroParameterParameterDependence : public ParameterParameterDependenceBase
+{
+public:
+
+	ConstantZeroParameterParameterDependence() { }
+	virtual ~ConstantZeroParameterParameterDependence() CADET_NOEXCEPT { }
+
+	static const char* identifier() { return "CONSTANT_ZERO"; }
+	virtual const char* name() const CADET_NOEXCEPT { return ConstantZeroParameterParameterDependence::identifier(); }
 
 	CADET_PARAMETERPARAMETERDEPENDENCE_BOILERPLATE
 
@@ -121,18 +179,57 @@ class DummyParameterParameterDependence : public ParameterParameterDependenceBas
 };
 
 
+/**
+ * @brief Defines a parameter dependence that outputs constant 1.0
+ */
+class ConstantOneParameterParameterDependence : public ParameterParameterDependenceBase
+{
+public:
+
+	ConstantOneParameterParameterDependence() { }
+	virtual ~ConstantOneParameterParameterDependence() CADET_NOEXCEPT { }
+
+	static const char* identifier() { return "CONSTANT_ONE"; }
+	virtual const char* name() const CADET_NOEXCEPT { return ConstantOneParameterParameterDependence::identifier(); }
+
+	CADET_PARAMETERPARAMETERDEPENDENCE_BOILERPLATE
+
+protected:
+
+	virtual bool configureImpl(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, BoundStateIdx bndIdx, const std::string& name)
+	{
+		return true;
+	}
+
+	template <typename ParamType>
+	ParamType getValueImpl(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd) const
+	{
+		return 1.0;
+	}
+
+	template <typename ParamType>
+	ParamType getValueImpl(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd, const ParamType& val) const
+	{
+		return 1.0;
+	}
+
+};
+
+
 namespace paramdep
 {
 	void registerDummyParamDependence(std::unordered_map<std::string, std::function<model::IParameterStateDependence*()>>& paramDeps)
 	{
-		paramDeps[DummyParameterStateDependence::identifier()] = []() { return new DummyParameterStateDependence(); };
-		paramDeps["NONE"] = []() { return new DummyParameterStateDependence(); };
+		paramDeps[ConstantOneParameterStateDependence::identifier()] = []() { return new ConstantOneParameterStateDependence(); };
+		paramDeps[ConstantZeroParameterStateDependence::identifier()] = []() { return new ConstantZeroParameterStateDependence(); };
+		paramDeps["NONE"] = []() { return new ConstantOneParameterStateDependence(); };
 	}
 
 	void registerDummyParamDependence(std::unordered_map<std::string, std::function<model::IParameterParameterDependence*()>>& paramDeps)
 	{
-		paramDeps[DummyParameterParameterDependence::identifier()] = []() { return new DummyParameterParameterDependence(); };
-		paramDeps["NONE"] = []() { return new DummyParameterParameterDependence(); };
+		paramDeps[ConstantOneParameterParameterDependence::identifier()] = []() { return new ConstantOneParameterParameterDependence(); };
+		paramDeps[ConstantZeroParameterParameterDependence::identifier()] = []() { return new ConstantZeroParameterParameterDependence(); };
+		paramDeps["NONE"] = []() { return new ConstantOneParameterParameterDependence(); };
 	}
 }  // namespace paramdep
 

src/libcadet/model/paramdep/PowerLawParameterDependence.cpp

@@ -51,7 +51,12 @@ class PowerLawParameterParameterDependence : public ParameterParameterDependence
 		const std::string baseName = name + "_BASE";
 		const std::string expName = name + "_EXPONENT";
 		const std::string flagName = name + "_ABS";
-		_base = paramProvider.getDouble(baseName);
+
+		if (paramProvider.exists(baseName))
+			_base = paramProvider.getDouble(baseName);
+		else
+			_base = 1.0;
+
 		_exponent = paramProvider.getDouble(expName);
 
 		if (paramProvider.exists(flagName))