diff --git a/src/libcadet/CMakeLists.txt b/src/libcadet/CMakeLists.txt
index 64f24d48d..31a5af600 100644
--- a/src/libcadet/CMakeLists.txt
+++ b/src/libcadet/CMakeLists.txt
@@ -146,7 +146,6 @@ set(LIBCADET_SOURCES
${CMAKE_SOURCE_DIR}/src/libcadet/model/paramdep/ParameterDependenceBase.cpp
${CMAKE_SOURCE_DIR}/src/libcadet/model/paramdep/LiquidSaltSolidParameterDependence.cpp
${CMAKE_SOURCE_DIR}/src/libcadet/model/paramdep/DummyParameterDependence.cpp
- ${CMAKE_SOURCE_DIR}/src/libcadet/model/paramdep/IdentityParameterDependence.cpp
${CMAKE_SOURCE_DIR}/src/libcadet/model/paramdep/PowerLawParameterDependence.cpp
)
diff --git a/src/libcadet/ParameterDependenceFactory.cpp b/src/libcadet/ParameterDependenceFactory.cpp
index ad165dfe9..359cdd14d 100644
--- a/src/libcadet/ParameterDependenceFactory.cpp
+++ b/src/libcadet/ParameterDependenceFactory.cpp
@@ -22,8 +22,6 @@ namespace cadet
void registerLiquidSaltSolidParamDependence(std::unordered_map<std::string, std::function<model::IParameterStateDependence*()>>& paramDeps);
void registerDummyParamDependence(std::unordered_map<std::string, std::function<model::IParameterStateDependence*()>>& paramDeps);
void registerDummyParamDependence(std::unordered_map<std::string, std::function<model::IParameterParameterDependence*()>>& paramDeps);
- void registerIdentityParamDependence(std::unordered_map<std::string, std::function<model::IParameterStateDependence*()>>& paramDeps);
- void registerIdentityParamDependence(std::unordered_map<std::string, std::function<model::IParameterParameterDependence*()>>& paramDeps);
void registerPowerLawParamDependence(std::unordered_map<std::string, std::function<model::IParameterParameterDependence*()>>& paramDeps);
}
}
@@ -33,11 +31,9 @@ namespace cadet
// Register all ParamState dependencies
model::paramdep::registerLiquidSaltSolidParamDependence(_paramStateDeps);
model::paramdep::registerDummyParamDependence(_paramStateDeps);
- model::paramdep::registerIdentityParamDependence(_paramStateDeps);
// Register all ParamParam dependencies
model::paramdep::registerDummyParamDependence(_paramParamDeps);
- model::paramdep::registerIdentityParamDependence(_paramParamDeps);
model::paramdep::registerPowerLawParamDependence(_paramParamDeps);
}
diff --git a/src/libcadet/model/GeneralRateModel.cpp b/src/libcadet/model/GeneralRateModel.cpp
index fef6a1c76..7c8dae0da 100644
--- a/src/libcadet/model/GeneralRateModel.cpp
+++ b/src/libcadet/model/GeneralRateModel.cpp
@@ -459,7 +459,7 @@ bool GeneralRateModel<ConvDispOperator>::configureModelDiscretization(IParameter
_filmDiffDep->configureModelDiscretization(paramProvider);
}
else
- _filmDiffDep = helper.createParameterParameterDependence("CONSTANT_ONE");
+ _filmDiffDep = helper.createParameterParameterDependence("IDENTITY");
// ==== Construct and configure binding model
clearBindingModels();
@@ -1926,11 +1926,11 @@ int GeneralRateModel<ConvDispOperator>::residualFlux(double t, unsigned int secI
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);
+ const ParamType filmDiffDep = static_cast<ParamType>(_filmDiffDep->getValue(*this, ColumnPosition{ relPos, 0.0, 0.0 }, comp, ParTypeIndep, BoundStateIndep, filmDiff[comp], curVelocity));
if (cadet_likely(_colParBoundaryOrder == 2))
- kf_FV[comp] = 1.0 / (absOuterShellHalfRadius / epsP / static_cast<ParamType>(_poreAccessFactor[type * _disc.nComp + comp]) / static_cast<ParamType>(parDiff[comp]) + 1.0 / (static_cast<ParamType>(filmDiff[comp]) * static_cast<ParamType>(modifier)));
+ kf_FV[comp] = 1.0 / (absOuterShellHalfRadius / epsP / static_cast<ParamType>(_poreAccessFactor[type * _disc.nComp + comp]) / static_cast<ParamType>(parDiff[comp]) + 1.0 / filmDiffDep);
else
- kf_FV[comp] = static_cast<ParamType>(filmDiff[comp]) * static_cast<ParamType>(modifier);
+ kf_FV[comp] = filmDiffDep;
}
// J_{0,f} block, adds flux to column void / bulk volume equations
@@ -1974,9 +1974,9 @@ int GeneralRateModel<ConvDispOperator>::residualFlux(double t, unsigned int secI
const double relPos = _convDispOp.relativeCoordinate(pblk);
const active curVelocity = _convDispOp.currentVelocity(relPos);
- const active modifier = _filmDiffDep->getValue(*this, ColumnPosition{ relPos, 0.0, 0.0 }, comp, ParTypeIndep, BoundStateIndep, curVelocity);
+ const ParamType filmDiffDep = static_cast<ParamType>(_filmDiffDep->getValue(*this, ColumnPosition{ relPos, 0.0, 0.0 }, comp, ParTypeIndep, BoundStateIndep, filmDiff[comp], curVelocity));
- kf_FV[comp] = (1.0 - static_cast<ParamType>(_parPorosity[type])) / (1.0 + epsP * static_cast<ParamType>(_poreAccessFactor[type * _disc.nComp + comp]) * static_cast<ParamType>(parDiff[comp]) / (absOuterShellHalfRadius * static_cast<ParamType>(filmDiff[comp]) * static_cast<ParamType>(modifier)));
+ kf_FV[comp] = (1.0 - static_cast<ParamType>(_parPorosity[type])) / (1.0 + epsP * static_cast<ParamType>(_poreAccessFactor[type * _disc.nComp + comp]) * static_cast<ParamType>(parDiff[comp]) / (absOuterShellHalfRadius * filmDiffDep));
}
const ColumnPosition colPos{(0.5 + static_cast<double>(pblk)) / static_cast<double>(_disc.nCol), 0.0, static_cast<double>(parCenterRadius[0]) / static_cast<double>(_parRadius[type])};
@@ -2089,12 +2089,12 @@ void GeneralRateModel<ConvDispOperator>::assembleOffdiagJac(double t, unsigned i
const double relPos = _convDispOp.relativeCoordinate(pblk);
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);
+ const double filmDiffDep = static_cast<double>(_filmDiffDep->getValue(*this, ColumnPosition{ relPos, 0.0, 0.0 }, comp, ParTypeIndep, BoundStateIndep, filmDiff[comp], curVelocity));
if (cadet_likely(_colParBoundaryOrder == 2))
- kf_FV[comp] = 1.0 / (absOuterShellHalfRadius / epsP / static_cast<double>(_poreAccessFactor[type * _disc.nComp + comp]) / static_cast<double>(parDiff[comp]) + 1.0 / (static_cast<double>(filmDiff[comp]) * modifier));
+ kf_FV[comp] = 1.0 / (absOuterShellHalfRadius / epsP / static_cast<double>(_poreAccessFactor[type * _disc.nComp + comp]) / static_cast<double>(parDiff[comp]) + 1.0 / filmDiffDep);
else
- kf_FV[comp] = static_cast<double>(filmDiff[comp]) * modifier;
+ kf_FV[comp] = filmDiffDep;
}
// J_{0,f} block, adds flux to column void / bulk volume equations
@@ -2148,9 +2148,9 @@ void GeneralRateModel<ConvDispOperator>::assembleOffdiagJac(double t, unsigned i
{
const double relPos = _convDispOp.relativeCoordinate(pblk);
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);
+ const double filmDiffDep = static_cast<double>(_filmDiffDep->getValue(*this, ColumnPosition{ relPos, 0.0, 0.0 }, comp, ParTypeIndep, BoundStateIndep, filmDiff[comp], curVelocity));
- kf_FV[comp] = (1.0 - static_cast<double>(_parPorosity[type])) / (1.0 + epsP * static_cast<double>(_poreAccessFactor[type * _disc.nComp + comp]) * static_cast<double>(parDiff[comp]) / (absOuterShellHalfRadius * static_cast<double>(filmDiff[comp]) * modifier));
+ kf_FV[comp] = (1.0 - static_cast<double>(_parPorosity[type])) / (1.0 + epsP * static_cast<double>(_poreAccessFactor[type * _disc.nComp + comp]) * static_cast<double>(parDiff[comp]) / (absOuterShellHalfRadius * filmDiffDep));
}
const ColumnPosition colPos{ (0.5 + static_cast<double>(pblk)) / static_cast<double>(_disc.nCol), 0.0, static_cast<double>(parCenterRadius[0]) / static_cast<double>(_parRadius[type]) };
diff --git a/src/libcadet/model/LumpedRateModelWithPores.cpp b/src/libcadet/model/LumpedRateModelWithPores.cpp
index a92857508..00ef01665 100644
--- a/src/libcadet/model/LumpedRateModelWithPores.cpp
+++ b/src/libcadet/model/LumpedRateModelWithPores.cpp
@@ -275,7 +275,7 @@ bool LumpedRateModelWithPores<ConvDispOperator>::configureModelDiscretization(IP
_filmDiffDep->configureModelDiscretization(paramProvider);
}
else
- _filmDiffDep = helper.createParameterParameterDependence("CONSTANT_ONE");
+ _filmDiffDep = helper.createParameterParameterDependence("IDENTITY");
// Allocate memory
Indexer idxr(_disc);
@@ -1091,9 +1091,9 @@ int LumpedRateModelWithPores<ConvDispOperator>::residualFlux(double t, unsigned
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);
+ const ParamType filmDiffDep = static_cast<ParamType>(_filmDiffDep->getValue(*this, ColumnPosition{relPos, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, filmDiff[comp], curVelocity));
- resCol[i] += jacCF_val * static_cast<ParamType>(filmDiff[comp]) * static_cast<ParamType>(modifier) * static_cast<ParamType>(_parTypeVolFrac[type + _disc.nParType * colCell]) * yFluxType[i];
+ resCol[i] += jacCF_val * filmDiffDep * static_cast<ParamType>(_parTypeVolFrac[type + _disc.nParType * colCell]) * yFluxType[i];
}
// J_{f,0} block, adds bulk volume state c_i to flux equation
@@ -1114,10 +1114,10 @@ int LumpedRateModelWithPores<ConvDispOperator>::residualFlux(double t, unsigned
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 ParamType filmDiffDep = static_cast<ParamType>(_filmDiffDep->getValue(*this, ColumnPosition{relPos, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, filmDiff[comp], 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]) * static_cast<ParamType>(modifier) * yFluxType[eq];
+ resParType[pblk * idxr.strideParBlock(type) + comp] += jacPF_val / static_cast<ParamType>(poreAccFactor[comp]) * filmDiffDep * yFluxType[eq];
}
}
@@ -1180,10 +1180,10 @@ void LumpedRateModelWithPores<ConvDispOperator>::assembleOffdiagJac(double t, un
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);
+ const double filmDiffDep = static_cast<double>(_filmDiffDep->getValue(*this, ColumnPosition{relPos, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, filmDiff[comp], curVelocity));
// Main diagonal corresponds to j_{f,i} (flux) state variable
- _jacCF.addElement(eq, eq + typeOffset, jacCF_val * static_cast<double>(filmDiff[comp]) * modifier * static_cast<double>(_parTypeVolFrac[type + _disc.nParType * colCell]));
+ _jacCF.addElement(eq, eq + typeOffset, jacCF_val * filmDiffDep * static_cast<double>(_parTypeVolFrac[type + _disc.nParType * colCell]));
}
// J_{f,0} block, adds bulk volume state c_i to flux equation
@@ -1203,8 +1203,8 @@ void LumpedRateModelWithPores<ConvDispOperator>::assembleOffdiagJac(double t, un
const unsigned int eq = typeOffset + pblk * idxr.strideColCell() + comp * idxr.strideColComp();
const unsigned int col = pblk * idxr.strideParBlock(type) + 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);
+ const double filmDiffDep = static_cast<double>(_filmDiffDep->getValue(*this, ColumnPosition{relPos, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, filmDiff[comp], curVelocity));
+ _jacPF[type].addElement(col, eq, jacPF_val / static_cast<double>(poreAccFactor[comp]) * filmDiffDep);
}
}
diff --git a/src/libcadet/model/ParameterDependence.hpp b/src/libcadet/model/ParameterDependence.hpp
index a463ca02a..1e981a2cc 100644
--- a/src/libcadet/model/ParameterDependence.hpp
+++ b/src/libcadet/model/ParameterDependence.hpp
@@ -460,9 +460,11 @@ class IParameterParameterDependence
* @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] depVal parameter-dependent value
+ * @param [in] indepVal parameter depVal depends on
* @return Actual parameter value
*/
- virtual double getValue(const IModel& model, 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, double depVal, double indepVal) const = 0;
/**
* @brief Evaluates the parameter
@@ -473,37 +475,11 @@ class IParameterParameterDependence
* @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] depVal parameter-dependent value
+ * @param [in] indepVal parameter depVal depends on
* @return Actual parameter value
*/
- 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] 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(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] 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(const IModel& model, 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& depVal, const active& indepVal) const = 0;
protected:
};
diff --git a/src/libcadet/model/paramdep/DummyParameterDependence.cpp b/src/libcadet/model/paramdep/DummyParameterDependence.cpp
index c73c92753..b406155b9 100644
--- a/src/libcadet/model/paramdep/DummyParameterDependence.cpp
+++ b/src/libcadet/model/paramdep/DummyParameterDependence.cpp
@@ -26,6 +26,36 @@ namespace cadet
namespace model
{
+/**
+ * @brief Defines a dummy parameter dependence that outputs the (independent) parameter itself
+ */
+class IdentityParameterParameterDependence : public ParameterParameterDependenceBase
+{
+public:
+
+ IdentityParameterParameterDependence() { }
+ virtual ~IdentityParameterParameterDependence() CADET_NOEXCEPT { }
+
+ static const char* identifier() { return "IDENTITY"; }
+ virtual const char* name() const CADET_NOEXCEPT { return IdentityParameterParameterDependence::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 ParamType& depVal, const ParamType& indepVal) const
+ {
+ return depVal;
+ }
+
+};
+
/**
* @brief Defines a parameter dependence that outputs constant 0.0
*/
@@ -83,7 +113,6 @@ class ConstantZeroParameterStateDependence : public ParameterStateDependenceBase
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 1.0
*/
@@ -141,81 +170,6 @@ class ConstantOneParameterStateDependence : public ParameterStateDependenceBase
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
-
-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 0.0;
- }
-
- template <typename ParamType>
- ParamType getValueImpl(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd, const ParamType& val) const
- {
- return 0.0;
- }
-
-};
-
-
-/**
- * @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)
@@ -227,9 +181,8 @@ namespace paramdep
void registerDummyParamDependence(std::unordered_map<std::string, std::function<model::IParameterParameterDependence*()>>& paramDeps)
{
- paramDeps[ConstantOneParameterParameterDependence::identifier()] = []() { return new ConstantOneParameterParameterDependence(); };
- paramDeps[ConstantZeroParameterParameterDependence::identifier()] = []() { return new ConstantZeroParameterParameterDependence(); };
- paramDeps["NONE"] = []() { return new ConstantOneParameterParameterDependence(); };
+ paramDeps[IdentityParameterParameterDependence::identifier()] = []() { return new IdentityParameterParameterDependence(); };
+ paramDeps["NONE"] = []() { return new IdentityParameterParameterDependence(); };
}
} // namespace paramdep
diff --git a/src/libcadet/model/paramdep/IdentityParameterDependence.cpp b/src/libcadet/model/paramdep/IdentityParameterDependence.cpp
deleted file mode 100644
index 60092ed6b..000000000
--- a/src/libcadet/model/paramdep/IdentityParameterDependence.cpp
+++ /dev/null
@@ -1,141 +0,0 @@
-// =============================================================================
-// CADET
-//
-// Copyright ┬® 2008-2022: The CADET Authors
-// Please see the AUTHORS and CONTRIBUTORS file.
-//
-// All rights reserved. This program and the accompanying materials
-// are made available under the terms of the GNU Public License v3.0 (or, at
-// your option, any later version) which accompanies this distribution, and
-// is available at http://www.gnu.org/licenses/gpl.html
-// =============================================================================
-
-#include "model/paramdep/ParameterDependenceBase.hpp"
-#include "ParamReaderHelper.hpp"
-#include "cadet/Exceptions.hpp"
-#include "SimulationTypes.hpp"
-#include "cadet/ParameterId.hpp"
-
-#include <sstream>
-#include <vector>
-#include <iomanip>
-
-namespace cadet
-{
-
-namespace model
-{
-
-/**
- * @brief Defines an identity parameter dependence
- */
-class IdentityParameterStateDependence : public ParameterStateDependenceBase
-{
-public:
-
- IdentityParameterStateDependence() { }
- virtual ~IdentityParameterStateDependence() CADET_NOEXCEPT { }
-
- static const char* identifier() { return "IDENTITY"; }
- virtual const char* name() const CADET_NOEXCEPT { return IdentityParameterStateDependence::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 param;
- }
-
- 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 param;
- }
-
- 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 param;
- }
-
- 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 an identity parameter dependence
- */
-class IdentityParameterParameterDependence : public ParameterParameterDependenceBase
-{
-public:
-
- IdentityParameterParameterDependence() { }
- virtual ~IdentityParameterParameterDependence() CADET_NOEXCEPT { }
-
- static const char* identifier() { return "IDENTITY"; }
- virtual const char* name() const CADET_NOEXCEPT { return IdentityParameterParameterDependence::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 0.0;
- }
-
- template <typename ParamType>
- ParamType getValueImpl(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd, const ParamType& val) const
- {
- return val;
- }
-
-};
-
-
-namespace paramdep
-{
- void registerIdentityParamDependence(std::unordered_map<std::string, std::function<model::IParameterStateDependence*()>>& paramDeps)
- {
- paramDeps[IdentityParameterStateDependence::identifier()] = []() { return new IdentityParameterStateDependence(); };
- paramDeps["NONE"] = []() { return new IdentityParameterStateDependence(); };
- }
-
- void registerIdentityParamDependence(std::unordered_map<std::string, std::function<model::IParameterParameterDependence*()>>& paramDeps)
- {
- paramDeps[IdentityParameterParameterDependence::identifier()] = []() { return new IdentityParameterParameterDependence(); };
- paramDeps["NONE"] = []() { return new IdentityParameterParameterDependence(); };
- }
-} // namespace paramdep
-
-} // namespace model
-
-} // namespace cadet
diff --git a/src/libcadet/model/paramdep/ParameterDependenceBase.hpp b/src/libcadet/model/paramdep/ParameterDependenceBase.hpp
index 1f72db332..80da1a2c3 100644
--- a/src/libcadet/model/paramdep/ParameterDependenceBase.hpp
+++ b/src/libcadet/model/paramdep/ParameterDependenceBase.hpp
@@ -241,27 +241,16 @@ class ParameterParameterDependenceBase : public IParameterParameterDependence
*
* The implementation is inserted inline in the class declaration.
*/
-#define CADET_PARAMETERPARAMETERDEPENDENCE_BOILERPLATE \
- virtual double getValue(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd, double val) const \
- { \
- return getValueImpl<double>(model, colPos, comp, parType, bnd, val); \
- } \
- \
- virtual active getValue(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd, const active& val) const \
- { \
- return getValueImpl<active>(model, colPos, comp, parType, bnd, val); \
- } \
- \
- virtual double getValue(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd) const \
- { \
- return getValueImpl<double>(model, colPos, comp, parType, bnd); \
- } \
- \
- virtual active getValueActive(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd) const \
- { \
- return getValueImpl<active>(model, colPos, comp, parType, bnd); \
- }
-
+#define CADET_PARAMETERPARAMETERDEPENDENCE_BOILERPLATE \
+ virtual double getValue(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd, double depVal, double indepVal) const \
+ { \
+ return getValueImpl<double>(model, colPos, comp, parType, bnd, depVal, indepVal); \
+ } \
+ \
+ virtual active getValue(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd, const active& depVal, const active& indepVal) const \
+ { \
+ return getValueImpl<active>(model, colPos, comp, parType, bnd, depVal, indepVal); \
+ } \
} // namespace model
} // namespace cadet
diff --git a/src/libcadet/model/paramdep/PowerLawParameterDependence.cpp b/src/libcadet/model/paramdep/PowerLawParameterDependence.cpp
index 99d9e0d2e..71fc27b9d 100644
--- a/src/libcadet/model/paramdep/PowerLawParameterDependence.cpp
+++ b/src/libcadet/model/paramdep/PowerLawParameterDependence.cpp
@@ -71,21 +71,15 @@ class PowerLawParameterParameterDependence : public ParameterParameterDependence
}
template <typename ParamType>
- ParamType getValueImpl(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd) const
- {
- return 0.0;
- }
-
- template <typename ParamType>
- ParamType getValueImpl(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd, ParamType val) const
+ ParamType getValueImpl(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd, const ParamType& depVal, const ParamType& indepVal) const
{
using std::pow;
using std::abs;
if (_useAbs)
- return static_cast<ParamType>(_base) * pow(abs(val), static_cast<ParamType>(_exponent));
+ return depVal * static_cast<ParamType>(_base) * pow(abs(indepVal), static_cast<ParamType>(_exponent));
else
- return static_cast<ParamType>(_base) * pow(val, static_cast<ParamType>(_exponent));
+ return depVal * static_cast<ParamType>(_base) * pow(indepVal, static_cast<ParamType>(_exponent));
}
};
diff --git a/src/libcadet/model/parts/AxialConvectionDispersionKernel.hpp b/src/libcadet/model/parts/AxialConvectionDispersionKernel.hpp
index 49d9fb31b..30a00aab3 100644
--- a/src/libcadet/model/parts/AxialConvectionDispersionKernel.hpp
+++ b/src/libcadet/model/parts/AxialConvectionDispersionKernel.hpp
@@ -124,7 +124,7 @@ namespace impl
if (cadet_likely(col < p.nCol - 1))
{
const double relCoord = static_cast<double>(col+1) / p.nCol;
- const ParamType d_ax_right = d_ax * p.parDep->getValue(p.model, ColumnPosition{relCoord, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, static_cast<ParamType>(p.u));
+ const ParamType d_ax_right = p.parDep->getValue(p.model, ColumnPosition{relCoord, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, d_ax, static_cast<ParamType>(p.u));
resBulkComp[col * p.strideCell] -= d_ax_right / h2 * (stencil[1] - stencil[0]);
// Jacobian entries
if (wantJac)
@@ -138,7 +138,7 @@ namespace impl
if (cadet_likely(col > 0))
{
const double relCoord = static_cast<double>(col) / p.nCol;
- const ParamType d_ax_left = d_ax * p.parDep->getValue(p.model, ColumnPosition{relCoord, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, static_cast<ParamType>(p.u));
+ const ParamType d_ax_left = p.parDep->getValue(p.model, ColumnPosition{relCoord, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, d_ax, static_cast<ParamType>(p.u));
resBulkComp[col * p.strideCell] -= d_ax_left / h2 * (stencil[-1] - stencil[0]);
// Jacobian entries
if (wantJac)
@@ -271,7 +271,7 @@ namespace impl
if (cadet_likely(col < p.nCol - 1))
{
const double relCoord = static_cast<double>(col+1) / p.nCol;
- const ParamType d_ax_right = d_ax * p.parDep->getValue(p.model, ColumnPosition{relCoord, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, static_cast<ParamType>(p.u));
+ const ParamType d_ax_right = p.parDep->getValue(p.model, ColumnPosition{relCoord, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, d_ax, static_cast<ParamType>(p.u));
resBulkComp[col * p.strideCell] -= d_ax_right / h2 * (stencil[-1] - stencil[0]);
// Jacobian entries
if (wantJac)
@@ -285,7 +285,7 @@ namespace impl
if (cadet_likely(col > 0))
{
const double relCoord = static_cast<double>(col) / p.nCol;
- const ParamType d_ax_left = d_ax * p.parDep->getValue(p.model, ColumnPosition{relCoord, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, static_cast<ParamType>(p.u));
+ const ParamType d_ax_left = p.parDep->getValue(p.model, ColumnPosition{relCoord, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, d_ax, static_cast<ParamType>(p.u));
resBulkComp[col * p.strideCell] -= d_ax_left / h2 * (stencil[1] - stencil[0]);
// Jacobian entries
if (wantJac)
diff --git a/src/libcadet/model/parts/ConvectionDispersionOperator.cpp b/src/libcadet/model/parts/ConvectionDispersionOperator.cpp
index 868069a4c..cdf707a56 100644
--- a/src/libcadet/model/parts/ConvectionDispersionOperator.cpp
+++ b/src/libcadet/model/parts/ConvectionDispersionOperator.cpp
@@ -78,7 +78,7 @@ bool AxialConvectionDispersionOperatorBase::configureModelDiscretization(IParame
_dispersionDep->configureModelDiscretization(paramProvider);
}
else
- _dispersionDep = helper.createParameterParameterDependence("CONSTANT_ONE");
+ _dispersionDep = helper.createParameterParameterDependence("IDENTITY");
paramProvider.pushScope("discretization");
@@ -577,7 +577,7 @@ bool RadialConvectionDispersionOperatorBase::configureModelDiscretization(IParam
_dispersionDep->configureModelDiscretization(paramProvider);
}
else
- _dispersionDep = helper.createParameterParameterDependence("CONSTANT_ONE");
+ _dispersionDep = helper.createParameterParameterDependence("IDENTITY");
paramProvider.pushScope("discretization");
diff --git a/src/libcadet/model/parts/RadialConvectionDispersionKernel.hpp b/src/libcadet/model/parts/RadialConvectionDispersionKernel.hpp
index 3ee40e366..686fe525e 100644
--- a/src/libcadet/model/parts/RadialConvectionDispersionKernel.hpp
+++ b/src/libcadet/model/parts/RadialConvectionDispersionKernel.hpp
@@ -126,7 +126,7 @@ namespace impl
if (cadet_likely(col < p.nCol - 1))
{
const double relCoord = (static_cast<double>(p.cellBounds[col+1]) - static_cast<double>(p.cellBounds[0])) / (static_cast<double>(p.cellBounds[p.nCol - 1]) - static_cast<double>(p.cellBounds[0]));
- const ParamType d_rad_right = d_rad * p.parDep->getValue(p.model, ColumnPosition{relCoord, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, static_cast<ParamType>(p.u) / static_cast<ParamType>(p.cellBounds[col+1]));
+ const ParamType d_rad_right = p.parDep->getValue(p.model, ColumnPosition{relCoord, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, d_rad, static_cast<ParamType>(p.u) / static_cast<ParamType>(p.cellBounds[col+1]));
resBulkComp[col * p.strideCell] -= d_rad_right * static_cast<ParamType>(p.cellBounds[col+1]) / denom * (stencil[1] - stencil[0]) / (static_cast<ParamType>(p.cellCenters[col+1]) - static_cast<ParamType>(p.cellCenters[col]));
// Jacobian entries
if (wantJac)
@@ -141,7 +141,7 @@ namespace impl
if (cadet_likely(col > 0))
{
const double relCoord = (static_cast<double>(p.cellBounds[col]) - static_cast<double>(p.cellBounds[0])) / (static_cast<double>(p.cellBounds[p.nCol - 1]) - static_cast<double>(p.cellBounds[0]));
- const ParamType d_rad_left = d_rad * p.parDep->getValue(p.model, ColumnPosition{relCoord, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, static_cast<ParamType>(p.u) / static_cast<ParamType>(p.cellBounds[col]));
+ const ParamType d_rad_left = p.parDep->getValue(p.model, ColumnPosition{relCoord, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, d_rad, static_cast<ParamType>(p.u) / static_cast<ParamType>(p.cellBounds[col]));
resBulkComp[col * p.strideCell] -= d_rad_left * static_cast<ParamType>(p.cellBounds[col]) / denom * (stencil[-1] - stencil[0]) / (static_cast<ParamType>(p.cellCenters[col-1]) - static_cast<ParamType>(p.cellCenters[col]));
// Jacobian entries
if (wantJac)
@@ -284,7 +284,7 @@ namespace impl
if (cadet_likely(col < p.nCol - 1))
{
const double relCoord = (static_cast<double>(p.cellBounds[col+1]) - static_cast<double>(p.cellBounds[0])) / (static_cast<double>(p.cellBounds[p.nCol - 1]) - static_cast<double>(p.cellBounds[0]));
- const ParamType d_rad_right = d_rad * p.parDep->getValue(p.model, ColumnPosition{relCoord, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, static_cast<ParamType>(p.u) / static_cast<ParamType>(p.cellBounds[col+1]));
+ const ParamType d_rad_right = p.parDep->getValue(p.model, ColumnPosition{relCoord, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, d_rad, static_cast<ParamType>(p.u) / static_cast<ParamType>(p.cellBounds[col+1]));
resBulkComp[col * p.strideCell] -= d_rad_right * static_cast<ParamType>(p.cellBounds[col+1]) / denom * (stencil[-1] - stencil[0]) / (static_cast<ParamType>(p.cellCenters[col+1]) - static_cast<ParamType>(p.cellCenters[col]));
// Jacobian entries
if (wantJac)
@@ -299,7 +299,7 @@ namespace impl
if (cadet_likely(col > 0))
{
const double relCoord = (static_cast<double>(p.cellBounds[col]) - static_cast<double>(p.cellBounds[0])) / (static_cast<double>(p.cellBounds[p.nCol - 1]) - static_cast<double>(p.cellBounds[0]));
- const ParamType d_rad_left = d_rad * p.parDep->getValue(p.model, ColumnPosition{relCoord, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, static_cast<ParamType>(p.u) / static_cast<ParamType>(p.cellBounds[col]));
+ const ParamType d_rad_left = p.parDep->getValue(p.model, ColumnPosition{relCoord, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, d_rad, static_cast<ParamType>(p.u) / static_cast<ParamType>(p.cellBounds[col]));
resBulkComp[col * p.strideCell] -= d_rad_left * static_cast<ParamType>(p.cellBounds[col]) / denom * (stencil[1] - stencil[0]) / (static_cast<ParamType>(p.cellCenters[col-1]) - static_cast<ParamType>(p.cellCenters[col]));
// Jacobian entries
if (wantJac)