diff --git a/doc/interface/parameter_dependencies.rst b/doc/interface/parameter_dependencies.rst new file mode 100644 index 000000000..c632de17a --- /dev/null +++ b/doc/interface/parameter_dependencies.rst @@ -0,0 +1,79 @@ +.. _parameter_dependencies: + +Parameter Dependencies +====================== + +Some parameters depend on other parameters (parameter-parameter dependency) or the solution variables (parameter-state dependency). +Parameter dependencies are defined in the unit operation scope. + +Parameter-Parameter Dependencies +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Group /input/model/unit_XXX +--------------------------- + +``COL_DISPERSION_DEP`` + + Parameter dependence of column dispersion on the interstitial velocity. Available for the LRM, LRMP and GRM units (with FV discretization only at the moment) + + ================ ===================================== ============= + **Type:** string **Range:** :math:`\texttt{POWER_LAW}` **Length:** 1 + ================ ===================================== ============= + +``FILM_DIFFUSION_DEP`` + + Parameter dependence of film diffusion on the interstitial velocity. Available for the LRMP unit (with FV discretization only at the moment) + + ================ ===================================== ============= + **Type:** string **Range:** :math:`\texttt{POWER_LAW}` **Length:** 1 + ================ ===================================== ============= + + +**Correlations** +"""""""""""""""" + +Different types of parameter correlations are can be applied. +The following correlations can be used for all parameter-parameter dependencies, but we specify the required input fields only for ``COL_DISPERSION_DEP``, for the sake of conciseness. + +**Power Law** + +.. math:: + + \begin{aligned} + p_{dep} &= p_{dep} \cdot b \ |p_{on}^x| + \end{aligned} + +Here, :math:`p_{dep}` is the dependent parameter and :math:`p_{on}` is the parameter it depends on. + +``COL_DISPERSION_DEP_BASE`` + + Base :math:`b` of the power law parameter dependence. Optional, defaults to :math:`1.0` + + ================ ============================= ============= + **Type:** double **Range:** :math:`\mathbb{R}` **Length:** 1 + ================ ============================= ============= + +``COL_DISPERSION_DEP_EXPONENT`` + + Exponent :math:`x` of the power law parameter dependence + + ================ ============================= ============= + **Type:** double **Range:** :math:`\mathbb{R}` **Length:** 1 + ================ ============================= ============= + +``COL_DISPERSION_DEP_ABS`` + + Specifies whether or not the absolute value should be computed. Optional, defaults to :math:`1` + + ============= =========================== ============= + **Type:** int **Range:** :math:`\{0, 1\}` **Length:** 1 + ============= =========================== ============= + + +Parameter-State Dependencies +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Group /input/model/unit_XXX +--------------------------- + +Parameter-State Dependencies are not fully implemented yet. diff --git a/doc/interface/unit_operations/radial_flow_models.rst b/doc/interface/unit_operations/radial_flow_models.rst index b40e74c17..3349785ad 100644 --- a/doc/interface/unit_operations/radial_flow_models.rst +++ b/doc/interface/unit_operations/radial_flow_models.rst @@ -34,7 +34,9 @@ For information on model equations, refer to :ref:`lumped_rate_model_without_por **Type:** double **Range:** :math:`\geq 0` **Length:** see :math:`\texttt{COL_DISPERSION_MULTIPLEX}` ================ ========================= ========================================================= - For multiplex specifications, e.g. for component dependency, see :ref:`general_rate_model_model`. + In addition to the multiplex specification (e.g. component dependency, see :ref:`general_rate_model_model`), the dispersion coefficient for radial flow model usually depends on other parameters. + Parameter dependencies are described here :ref:`parameter_dependencies`. + ``COL_RADIUS_INNER`` diff --git a/src/libcadet/model/GeneralRateModel.cpp b/src/libcadet/model/GeneralRateModel.cpp index 6d8f6e8fb..5154f6b37 100644 --- a/src/libcadet/model/GeneralRateModel.cpp +++ b/src/libcadet/model/GeneralRateModel.cpp @@ -1263,7 +1263,7 @@ template template int GeneralRateModel::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::enabled()); + _convDispOp.residual(*this, t, secIdx, yBase, yDotBase, resBase, wantJac, typename ParamSens::enabled()); if (!_dynReactionBulk || (_dynReactionBulk->numReactionsLiquid() == 0)) return 0; diff --git a/src/libcadet/model/GeneralRateModel2D.cpp b/src/libcadet/model/GeneralRateModel2D.cpp index 05f401179..7bd46dccf 100644 --- a/src/libcadet/model/GeneralRateModel2D.cpp +++ b/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 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::enabled()); + _convDispOp.residual(*this, t, secIdx, yBase, yDotBase, resBase, wantJac, typename ParamSens::enabled()); if (!_dynReactionBulk || (_dynReactionBulk->numReactionsLiquid() == 0)) return 0; diff --git a/src/libcadet/model/LumpedRateModelWithPores.cpp b/src/libcadet/model/LumpedRateModelWithPores.cpp index e39190825..a92857508 100644 --- a/src/libcadet/model/LumpedRateModelWithPores.cpp +++ b/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 LumpedRateModelWithPores::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::~LumpedRateModelWithPores() CADET_NO delete[] _tempState; delete _dynReactionBulk; + delete _filmDiffDep; delete[] _disc.parTypeOffset; delete[] _disc.nBound; @@ -263,6 +265,18 @@ bool LumpedRateModelWithPores::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::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 template int LumpedRateModelWithPores::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::enabled()); + _convDispOp.residual(*this, t, secIdx, yBase, yDotBase, resBase, wantJac, typename ParamSens::enabled()); if (!_dynReactionBulk || (_dynReactionBulk->numReactionsLiquid() == 0)) return 0; @@ -1068,7 +1088,12 @@ int LumpedRateModelWithPores::residualFlux(double t, unsigned { const unsigned int colCell = i / _disc.nComp; const unsigned int comp = i % _disc.nComp; - resCol[i] += jacCF_val * static_cast(filmDiff[comp]) * static_cast(_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(filmDiff[comp]) * static_cast(modifier) * static_cast(_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::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(poreAccFactor[comp]) * static_cast(filmDiff[comp]) * yFluxType[eq]; + resParType[pblk * idxr.strideParBlock(type) + comp] += jacPF_val / static_cast(poreAccFactor[comp]) * static_cast(filmDiff[comp]) * static_cast(modifier) * yFluxType[eq]; } } @@ -1148,8 +1178,12 @@ void LumpedRateModelWithPores::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(_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(filmDiff[comp]) * static_cast(_parTypeVolFrac[type + _disc.nParType * colCell])); + _jacCF.addElement(eq, eq + typeOffset, jacCF_val * static_cast(filmDiff[comp]) * modifier * static_cast(_parTypeVolFrac[type + _disc.nParType * colCell])); } // J_{f,0} block, adds bulk volume state c_i to flux equation @@ -1161,11 +1195,16 @@ void LumpedRateModelWithPores::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(_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(poreAccFactor[comp]) * static_cast(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(poreAccFactor[comp]) * static_cast(filmDiff[comp]) * modifier); } } @@ -1468,6 +1507,15 @@ bool LumpedRateModelWithPores::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::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::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); diff --git a/src/libcadet/model/LumpedRateModelWithPores.hpp b/src/libcadet/model/LumpedRateModelWithPores.hpp index f21e2c898..1fcc6fbd6 100644 --- a/src/libcadet/model/LumpedRateModelWithPores.hpp +++ b/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 _jacP; //!< Particle jacobian diagonal blocks (all of them for each particle type) std::vector _jacPdisc; //!< Particle jacobian diagonal blocks (all of them for each particle type) with time derivatives from BDF method diff --git a/src/libcadet/model/LumpedRateModelWithoutPores.cpp b/src/libcadet/model/LumpedRateModelWithoutPores.cpp index 7b044419d..4e04bf1a0 100644 --- a/src/libcadet/model/LumpedRateModelWithoutPores.cpp +++ b/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,9 +57,9 @@ namespace class ConvOpResidual { 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); } }; @@ -67,9 +67,9 @@ namespace class ConvOpResidual { 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::enabled()); + op.residual(*model, t, secIdx, y, yDot, res, typename cadet::ParamSens::enabled()); } }; @@ -77,9 +77,9 @@ namespace class ConvOpResidual { 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::enabled()); + op.residual(*model, t, secIdx, y, yDot, res, typename cadet::ParamSens::enabled()); } }; } @@ -609,7 +609,7 @@ template template int LumpedRateModelWithoutPores::residualImpl(double t, unsigned int secIdx, StateType const* const y, double const* const yDot, ResidualType* const res, util::ThreadLocalStorage& threadLocalMem) { - ConvOpResidual::call(_convDispOp, t, secIdx, y, yDot, res, _jac); + ConvOpResidual::call(this, _convDispOp, t, secIdx, y, yDot, res, _jac); Indexer idxr(_disc); @@ -785,7 +785,7 @@ template unsigned int LumpedRateModelWithoutPores::localOutletComponentIndex(unsigned int port) const CADET_NOEXCEPT { // Inlets are duplicated so need to be accounted for - if (static_cast(_convDispOp.currentVelocity()) >= 0.0) + if (_convDispOp.forwardFlow()) // Forward Flow: outlet is last cell return _disc.nComp + (_disc.nCol - 1) * (_disc.nComp + _disc.strideBound); else diff --git a/src/libcadet/model/ParameterDependence.hpp b/src/libcadet/model/ParameterDependence.hpp index 9363cf2f7..a463ca02a 100644 --- a/src/libcadet/model/ParameterDependence.hpp +++ b/src/libcadet/model/ParameterDependence.hpp @@ -33,6 +33,7 @@ namespace cadet class IParameterProvider; struct ColumnPosition; +class IModel; namespace model { @@ -454,36 +455,33 @@ 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& 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& 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) @@ -491,14 +489,13 @@ class IParameterParameterDependence * @param [in] val Additional parameter-dependent value * @return Actual parameter value */ - virtual double getValue(UnitOpIdx unitOpIdx, const std::unordered_map& 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) @@ -506,7 +503,7 @@ class IParameterParameterDependence * @param [in] val Additional parameter-dependent value * @return Actual parameter value */ - virtual active getValue(UnitOpIdx unitOpIdx, const std::unordered_map& 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: }; diff --git a/src/libcadet/model/paramdep/DummyParameterDependence.cpp b/src/libcadet/model/paramdep/DummyParameterDependence.cpp index 36c331411..c73c92753 100644 --- a/src/libcadet/model/paramdep/DummyParameterDependence.cpp +++ b/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 DoubleActivePromoter::type liquidParameterImpl(const ColumnPosition& colPos, const ParamType& param, StateType const* y, int comp) const + { + return 1.0; + } + + template + void analyticJacobianLiquidAddImpl(const ColumnPosition& colPos, double param, double const* y, int comp, double factor, int offset, RowIterator jac) const { } + + template + typename DoubleActivePromoter::type combinedParameterLiquidImpl(const ColumnPosition& colPos, const ParamType& param, StateType const* yLiquid, StateType const* ySolid, int comp) const + { + return 1.0; + } + + template + void analyticJacobianCombinedAddLiquidImpl(const ColumnPosition& colPos, double param, double const* yLiquid, double const* ySolid, int comp, double factor, int offset, RowIterator jac) const { } + + template + typename DoubleActivePromoter::type combinedParameterSolidImpl(const ColumnPosition& colPos, const ParamType& param, StateType const* yLiquid, StateType const* ySolid, int bnd) const + { + return 1.0; + } + + template + 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 @@ -107,13 +165,13 @@ class DummyParameterParameterDependence : public ParameterParameterDependenceBas } template - ParamType getValueImpl(UnitOpIdx unitOpIdx, const std::unordered_map& params, const ColumnPosition& colPos, int comp, int parType, int bnd) const + ParamType getValueImpl(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd) const { return 0.0; } template - ParamType getValueImpl(UnitOpIdx unitOpIdx, const std::unordered_map& params, const ColumnPosition& colPos, int comp, int parType, int bnd, const ParamType& val) const + ParamType getValueImpl(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd, const ParamType& val) const { return 0.0; } @@ -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 + ParamType getValueImpl(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd) const + { + return 1.0; + } + + template + 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>& 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>& 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 diff --git a/src/libcadet/model/paramdep/ParameterDependenceBase.hpp b/src/libcadet/model/paramdep/ParameterDependenceBase.hpp index 1adf69c04..1f72db332 100644 --- a/src/libcadet/model/paramdep/ParameterDependenceBase.hpp +++ b/src/libcadet/model/paramdep/ParameterDependenceBase.hpp @@ -241,25 +241,25 @@ class ParameterParameterDependenceBase : public IParameterParameterDependence * * The implementation is inserted inline in the class declaration. */ -#define CADET_PARAMETERPARAMETERDEPENDENCE_BOILERPLATE \ - virtual double getValue(UnitOpIdx unitOpIdx, const std::unordered_map& params, const ColumnPosition& colPos, int comp, int parType, int bnd, double val) const \ - { \ - return getValueImpl(unitOpIdx, params, colPos, comp, parType, bnd, val); \ - } \ - \ - virtual active getValue(UnitOpIdx unitOpIdx, const std::unordered_map& params, const ColumnPosition& colPos, int comp, int parType, int bnd, const active& val) const \ - { \ - return getValueImpl(unitOpIdx, params, colPos, comp, parType, bnd, val); \ - } \ - \ - virtual double getValue(UnitOpIdx unitOpIdx, const std::unordered_map& params, const ColumnPosition& colPos, int comp, int parType, int bnd) const \ - { \ - return getValueImpl(unitOpIdx, params, colPos, comp, parType, bnd); \ - } \ - \ - virtual active getValueActive(UnitOpIdx unitOpIdx, const std::unordered_map& params, const ColumnPosition& colPos, int comp, int parType, int bnd) const \ - { \ - return getValueImpl(unitOpIdx, params, colPos, comp, parType, bnd); \ +#define CADET_PARAMETERPARAMETERDEPENDENCE_BOILERPLATE \ + virtual double getValue(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd, double val) const \ + { \ + return getValueImpl(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(model, colPos, comp, parType, bnd, val); \ + } \ + \ + virtual double getValue(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd) const \ + { \ + return getValueImpl(model, colPos, comp, parType, bnd); \ + } \ + \ + virtual active getValueActive(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd) const \ + { \ + return getValueImpl(model, colPos, comp, parType, bnd); \ } diff --git a/src/libcadet/model/paramdep/PowerLawParameterDependence.cpp b/src/libcadet/model/paramdep/PowerLawParameterDependence.cpp index 91bb36c05..99d9e0d2e 100644 --- a/src/libcadet/model/paramdep/PowerLawParameterDependence.cpp +++ b/src/libcadet/model/paramdep/PowerLawParameterDependence.cpp @@ -44,14 +44,26 @@ class PowerLawParameterParameterDependence : public ParameterParameterDependence protected: active _base; active _exponent; + bool _useAbs; virtual bool configureImpl(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, BoundStateIdx bndIdx, const std::string& name) { const std::string baseName = name + "_BASE"; const std::string expName = name + "_EXPONENT"; - _base = paramProvider.getDouble(baseName); + const std::string flagName = name + "_ABS"; + + if (paramProvider.exists(baseName)) + _base = paramProvider.getDouble(baseName); + else + _base = 1.0; + _exponent = paramProvider.getDouble(expName); + if (paramProvider.exists(flagName)) + _useAbs = paramProvider.getBool(flagName); + else + _useAbs = true; + _parameters[makeParamId(hashStringRuntime(baseName), unitOpIdx, CompIndep, parTypeIdx, bndIdx, ReactionIndep, SectionIndep)] = &_base; _parameters[makeParamId(hashStringRuntime(expName), unitOpIdx, CompIndep, parTypeIdx, bndIdx, ReactionIndep, SectionIndep)] = &_exponent; @@ -59,16 +71,21 @@ class PowerLawParameterParameterDependence : public ParameterParameterDependence } template - ParamType getValueImpl(UnitOpIdx unitOpIdx, const std::unordered_map& params, const ColumnPosition& colPos, int comp, int parType, int bnd) const + ParamType getValueImpl(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd) const { return 0.0; } template - ParamType getValueImpl(UnitOpIdx unitOpIdx, const std::unordered_map& params, 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, ParamType val) const { using std::pow; - return static_cast(_base) * pow(val, static_cast(_exponent)); + using std::abs; + + if (_useAbs) + return static_cast(_base) * pow(abs(val), static_cast(_exponent)); + else + return static_cast(_base) * pow(val, static_cast(_exponent)); } }; diff --git a/src/libcadet/model/parts/AxialConvectionDispersionKernel.hpp b/src/libcadet/model/parts/AxialConvectionDispersionKernel.hpp index 42e1142f0..49d9fb31b 100644 --- a/src/libcadet/model/parts/AxialConvectionDispersionKernel.hpp +++ b/src/libcadet/model/parts/AxialConvectionDispersionKernel.hpp @@ -24,6 +24,8 @@ #include "Stencil.hpp" #include "linalg/CompressedSparseMatrix.hpp" #include "SimulationTypes.hpp" +#include "model/ParameterDependence.hpp" +#include "model/UnitOperation.hpp" namespace cadet { @@ -52,6 +54,8 @@ struct AxialFlowParameters unsigned int nCol; unsigned int offsetToInlet; //!< Offset to the first component of the inlet DOFs in the local state vector unsigned int offsetToBulk; //!< Offset to the first component of the first bulk cell in the local state vector + IParameterParameterDependence* parDep; + const IModel& model; }; @@ -119,24 +123,28 @@ namespace impl // Right side, leave out if we're in the last cell (boundary condition) if (cadet_likely(col < p.nCol - 1)) { - resBulkComp[col * p.strideCell] -= d_ax / h2 * (stencil[1] - stencil[0]); + const double relCoord = static_cast(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(p.u)); + resBulkComp[col * p.strideCell] -= d_ax_right / h2 * (stencil[1] - stencil[0]); // Jacobian entries if (wantJac) { - jac[0] += static_cast(d_ax) / static_cast(h2); - jac[p.strideCell] -= static_cast(d_ax) / static_cast(h2); + jac[0] += static_cast(d_ax_right) / static_cast(h2); + jac[p.strideCell] -= static_cast(d_ax_right) / static_cast(h2); } } // Left side, leave out if we're in the first cell (boundary condition) if (cadet_likely(col > 0)) { - resBulkComp[col * p.strideCell] -= d_ax / h2 * (stencil[-1] - stencil[0]); + const double relCoord = static_cast(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(p.u)); + resBulkComp[col * p.strideCell] -= d_ax_left / h2 * (stencil[-1] - stencil[0]); // Jacobian entries if (wantJac) { - jac[0] += static_cast(d_ax) / static_cast(h2); - jac[-p.strideCell] -= static_cast(d_ax) / static_cast(h2); + jac[0] += static_cast(d_ax_left) / static_cast(h2); + jac[-p.strideCell] -= static_cast(d_ax_left) / static_cast(h2); } } @@ -262,24 +270,28 @@ namespace impl // Right side, leave out if we're in the first cell (boundary condition) if (cadet_likely(col < p.nCol - 1)) { - resBulkComp[col * p.strideCell] -= d_ax / h2 * (stencil[-1] - stencil[0]); + const double relCoord = static_cast(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(p.u)); + resBulkComp[col * p.strideCell] -= d_ax_right / h2 * (stencil[-1] - stencil[0]); // Jacobian entries if (wantJac) { - jac[0] += static_cast(d_ax) / static_cast(h2); - jac[p.strideCell] -= static_cast(d_ax) / static_cast(h2); + jac[0] += static_cast(d_ax_right) / static_cast(h2); + jac[p.strideCell] -= static_cast(d_ax_right) / static_cast(h2); } } // Left side, leave out if we're in the last cell (boundary condition) if (cadet_likely(col > 0)) { - resBulkComp[col * p.strideCell] -= d_ax / h2 * (stencil[1] - stencil[0]); + const double relCoord = static_cast(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(p.u)); + resBulkComp[col * p.strideCell] -= d_ax_left / h2 * (stencil[1] - stencil[0]); // Jacobian entries if (wantJac) { - jac[0] += static_cast(d_ax) / static_cast(h2); - jac[-p.strideCell] -= static_cast(d_ax) / static_cast(h2); + jac[0] += static_cast(d_ax_left) / static_cast(h2); + jac[-p.strideCell] -= static_cast(d_ax_left) / static_cast(h2); } } diff --git a/src/libcadet/model/parts/ConvectionDispersionOperator.cpp b/src/libcadet/model/parts/ConvectionDispersionOperator.cpp index c93cef318..0e508127d 100644 --- a/src/libcadet/model/parts/ConvectionDispersionOperator.cpp +++ b/src/libcadet/model/parts/ConvectionDispersionOperator.cpp @@ -19,6 +19,7 @@ #include "SimulationTypes.hpp" #include "model/parts/AxialConvectionDispersionKernel.hpp" #include "model/parts/RadialConvectionDispersionKernel.hpp" +#include "model/ParameterDependence.hpp" #include "SensParamUtil.hpp" #include "ConfigurationHelper.hpp" @@ -41,12 +42,15 @@ namespace parts * @brief Creates an AxialConvectionDispersionOperatorBase */ AxialConvectionDispersionOperatorBase::AxialConvectionDispersionOperatorBase() : _stencilMemory(sizeof(active) * Weno::maxStencilSize()), - _wenoDerivatives(new double[Weno::maxStencilSize()]), _weno() + _wenoDerivatives(new double[Weno::maxStencilSize()]), _weno(), _dispersionDep(nullptr) { } AxialConvectionDispersionOperatorBase::~AxialConvectionDispersionOperatorBase() CADET_NOEXCEPT { + if (_dispersionDep) + delete _dispersionDep; + delete[] _wenoDerivatives; } @@ -64,6 +68,18 @@ bool AxialConvectionDispersionOperatorBase::configureModelDiscretization(IParame _nCol = nCol; _strideCell = strideCell; + if (paramProvider.exists("COL_DISPERSION_DEP")) + { + const std::string paramDepName = paramProvider.getString("COL_DISPERSION_DEP"); + _dispersionDep = helper.createParameterParameterDependence(paramDepName); + if (!_dispersionDep) + throw InvalidParameterException("Unknown parameter dependence " + paramDepName + " in COL_DISPERSION_DEP"); + + _dispersionDep->configureModelDiscretization(paramProvider); + } + else + _dispersionDep = helper.createParameterParameterDependence("CONSTANT_ONE"); + paramProvider.pushScope("discretization"); // Read WENO settings and apply them @@ -149,6 +165,12 @@ bool AxialConvectionDispersionOperatorBase::configure(UnitOpIdx unitOpIdx, IPara _colDispersion = std::move(expanded); } + if (_dispersionDep) + { + if (!_dispersionDep->configure(paramProvider, unitOpIdx, ParTypeIndep, BoundStateIndep, "COL_DISPERSION_DEP")) + throw InvalidParameterException("Failed to configure dispersion parameter dependency (COL_DISPERSION_DEP)"); + } + if (_velocity.empty() && (_crossSection <= 0.0)) { throw InvalidParameterException("At least one of CROSS_SECTION_AREA and VELOCITY has to be set"); @@ -233,6 +255,7 @@ void AxialConvectionDispersionOperatorBase::setFlowRates(const active& in, const /** * @brief Computes the residual of the transport equations + * @param [in] model Model that owns the operator * @param [in] t Current time point * @param [in] secIdx Index of the current section * @param [in] y Pointer to unit operation's state vector @@ -241,44 +264,44 @@ void AxialConvectionDispersionOperatorBase::setFlowRates(const active& in, const * @param [in] jac Matrix that holds the Jacobian * @return @c 0 on success, @c -1 on non-recoverable error, and @c +1 on recoverable error */ -int AxialConvectionDispersionOperatorBase::residual(double t, unsigned int secIdx, double const* y, double const* yDot, double* res, linalg::BandMatrix& jac) +int AxialConvectionDispersionOperatorBase::residual(const IModel& model, double t, unsigned int secIdx, double const* y, double const* yDot, double* res, linalg::BandMatrix& jac) { // Reset Jacobian jac.setAll(0.0); - return residualImpl(t, secIdx, y, yDot, res, jac.row(0)); + return residualImpl(model, t, secIdx, y, yDot, res, jac.row(0)); } -int AxialConvectionDispersionOperatorBase::residual(double t, unsigned int secIdx, double const* y, double const* yDot, double* res, WithoutParamSensitivity) +int AxialConvectionDispersionOperatorBase::residual(const IModel& model, double t, unsigned int secIdx, double const* y, double const* yDot, double* res, WithoutParamSensitivity) { - return residualImpl(t, secIdx, y, yDot, res, linalg::BandMatrix::RowIterator()); + return residualImpl(model, t, secIdx, y, yDot, res, linalg::BandMatrix::RowIterator()); } -int AxialConvectionDispersionOperatorBase::residual(double t, unsigned int secIdx, active const* y, double const* yDot, active* res, WithoutParamSensitivity) +int AxialConvectionDispersionOperatorBase::residual(const IModel& model, double t, unsigned int secIdx, active const* y, double const* yDot, active* res, WithoutParamSensitivity) { - return residualImpl(t, secIdx, y, yDot, res, linalg::BandMatrix::RowIterator()); + return residualImpl(model, t, secIdx, y, yDot, res, linalg::BandMatrix::RowIterator()); } -int AxialConvectionDispersionOperatorBase::residual(double t, unsigned int secIdx, double const* y, double const* yDot, active* res, linalg::BandMatrix& jac) +int AxialConvectionDispersionOperatorBase::residual(const IModel& model, double t, unsigned int secIdx, double const* y, double const* yDot, active* res, linalg::BandMatrix& jac) { // Reset Jacobian jac.setAll(0.0); - return residualImpl(t, secIdx, y, yDot, res, jac.row(0)); + return residualImpl(model, t, secIdx, y, yDot, res, jac.row(0)); } -int AxialConvectionDispersionOperatorBase::residual(double t, unsigned int secIdx, double const* y, double const* yDot, active* res, WithParamSensitivity) +int AxialConvectionDispersionOperatorBase::residual(const IModel& model, double t, unsigned int secIdx, double const* y, double const* yDot, active* res, WithParamSensitivity) { - return residualImpl(t, secIdx, y, yDot, res, linalg::BandMatrix::RowIterator()); + return residualImpl(model, t, secIdx, y, yDot, res, linalg::BandMatrix::RowIterator()); } -int AxialConvectionDispersionOperatorBase::residual(double t, unsigned int secIdx, active const* y, double const* yDot, active* res, WithParamSensitivity) +int AxialConvectionDispersionOperatorBase::residual(const IModel& model, double t, unsigned int secIdx, active const* y, double const* yDot, active* res, WithParamSensitivity) { - return residualImpl(t, secIdx, y, yDot, res, linalg::BandMatrix::RowIterator()); + return residualImpl(model, t, secIdx, y, yDot, res, linalg::BandMatrix::RowIterator()); } template -int AxialConvectionDispersionOperatorBase::residualImpl(double t, unsigned int secIdx, StateType const* y, double const* yDot, ResidualType* res, RowIteratorType jacBegin) +int AxialConvectionDispersionOperatorBase::residualImpl(const IModel& model, double t, unsigned int secIdx, StateType const* y, double const* yDot, ResidualType* res, RowIteratorType jacBegin) { const ParamType u = static_cast(_curVelocity); active const* const d_c = getSectionDependentSlice(_colDispersion, _nComp, secIdx); @@ -297,7 +320,9 @@ int AxialConvectionDispersionOperatorBase::residualImpl(double t, unsigned int s _nComp, _nCol, 0u, - _nComp + _nComp, + _dispersionDep, + model }; return convdisp::residualKernelAxial(SimulationTime{t, secIdx}, y, yDot, res, jacBegin, fp); @@ -381,6 +406,16 @@ double AxialConvectionDispersionOperatorBase::inletJacobianFactor() const CADET_ bool AxialConvectionDispersionOperatorBase::setParameter(const ParameterId& pId, double value) { + // Check if parameter is in parameter dependence of column dispersion coefficient + if (_dispersionDep) + { + if (_dispersionDep->hasParameter(pId)) + { + _dispersionDep->setParameter(pId, value); + return true; + } + } + // We only need to do something if COL_DISPERSION is component independent if (!_dispersionCompIndep) return false; @@ -412,6 +447,17 @@ bool AxialConvectionDispersionOperatorBase::setParameter(const ParameterId& pId, bool AxialConvectionDispersionOperatorBase::setSensitiveParameterValue(const std::unordered_set& sensParams, const ParameterId& pId, double value) { + // Check if parameter is in parameter dependence of column dispersion coefficient + if (_dispersionDep) + { + active* const param = _dispersionDep->getParameter(pId); + if (param) + { + param->setValue(value); + return true; + } + } + // We only need to do something if COL_DISPERSION is component independent if (!_dispersionCompIndep) return false; @@ -449,6 +495,17 @@ bool AxialConvectionDispersionOperatorBase::setSensitiveParameterValue(const std bool AxialConvectionDispersionOperatorBase::setSensitiveParameter(std::unordered_set& sensParams, const ParameterId& pId, unsigned int adDirection, double adValue) { + // Check if parameter is in parameter dependence of column dispersion coefficient + if (_dispersionDep) + { + active* const param = _dispersionDep->getParameter(pId); + if (param) + { + param->setADValue(adDirection, adValue); + return true; + } + } + // We only need to do something if COL_DISPERSION is component independent if (!_dispersionCompIndep) return false; @@ -486,12 +543,14 @@ bool AxialConvectionDispersionOperatorBase::setSensitiveParameter(std::unordered /** * @brief Creates a RadialConvectionDispersionOperatorBase */ -RadialConvectionDispersionOperatorBase::RadialConvectionDispersionOperatorBase() : _stencilMemory(sizeof(active) * Weno::maxStencilSize()) +RadialConvectionDispersionOperatorBase::RadialConvectionDispersionOperatorBase() : _stencilMemory(sizeof(active) * Weno::maxStencilSize()), _dispersionDep(nullptr) { } RadialConvectionDispersionOperatorBase::~RadialConvectionDispersionOperatorBase() CADET_NOEXCEPT { + if (_dispersionDep) + delete _dispersionDep; } /** @@ -508,6 +567,18 @@ bool RadialConvectionDispersionOperatorBase::configureModelDiscretization(IParam _nCol = nCol; _strideCell = strideCell; + if (paramProvider.exists("COL_DISPERSION_DEP")) + { + const std::string paramDepName = paramProvider.getString("COL_DISPERSION_DEP"); + _dispersionDep = helper.createParameterParameterDependence(paramDepName); + if (!_dispersionDep) + throw InvalidParameterException("Unknown parameter dependence " + paramDepName + " in COL_DISPERSION_DEP"); + + _dispersionDep->configureModelDiscretization(paramProvider); + } + else + _dispersionDep = helper.createParameterParameterDependence("CONSTANT_ONE"); + paramProvider.pushScope("discretization"); // Read WENO settings and apply them @@ -596,6 +667,12 @@ bool RadialConvectionDispersionOperatorBase::configure(UnitOpIdx unitOpIdx, IPar _colDispersion = std::move(expanded); } + if (_dispersionDep) + { + if (!_dispersionDep->configure(paramProvider, unitOpIdx, ParTypeIndep, BoundStateIndep, "COL_DISPERSION_DEP")) + throw InvalidParameterException("Failed to configure dispersion parameter dependency (COL_DISPERSION_DEP)"); + } + if (_velocity.empty() && (_colLength <= 0.0)) { throw InvalidParameterException("At least one of COL_LENGTH and VELOCITY_COEFF has to be set"); @@ -680,8 +757,15 @@ void RadialConvectionDispersionOperatorBase::setFlowRates(const active& in, cons _curVelocity = _dir * in / (2.0 * pi * _colLength * colPorosity); } +active RadialConvectionDispersionOperatorBase::currentVelocity(double pos) const CADET_NOEXCEPT +{ + const active radius = pos * (_outerRadius - _innerRadius) + _innerRadius; + return _curVelocity / radius; +} + /** * @brief Computes the residual of the transport equations + * @param [in] model Model that owns the operator * @param [in] t Current time point * @param [in] secIdx Index of the current section * @param [in] y Pointer to unit operation's state vector @@ -690,44 +774,44 @@ void RadialConvectionDispersionOperatorBase::setFlowRates(const active& in, cons * @param [in] jac Matrix that holds the Jacobian * @return @c 0 on success, @c -1 on non-recoverable error, and @c +1 on recoverable error */ -int RadialConvectionDispersionOperatorBase::residual(double t, unsigned int secIdx, double const* y, double const* yDot, double* res, linalg::BandMatrix& jac) +int RadialConvectionDispersionOperatorBase::residual(const IModel& model, double t, unsigned int secIdx, double const* y, double const* yDot, double* res, linalg::BandMatrix& jac) { // Reset Jacobian jac.setAll(0.0); - return residualImpl(t, secIdx, y, yDot, res, jac.row(0)); + return residualImpl(model, t, secIdx, y, yDot, res, jac.row(0)); } -int RadialConvectionDispersionOperatorBase::residual(double t, unsigned int secIdx, double const* y, double const* yDot, double* res, WithoutParamSensitivity) +int RadialConvectionDispersionOperatorBase::residual(const IModel& model, double t, unsigned int secIdx, double const* y, double const* yDot, double* res, WithoutParamSensitivity) { - return residualImpl(t, secIdx, y, yDot, res, linalg::BandMatrix::RowIterator()); + return residualImpl(model, t, secIdx, y, yDot, res, linalg::BandMatrix::RowIterator()); } -int RadialConvectionDispersionOperatorBase::residual(double t, unsigned int secIdx, active const* y, double const* yDot, active* res, WithoutParamSensitivity) +int RadialConvectionDispersionOperatorBase::residual(const IModel& model, double t, unsigned int secIdx, active const* y, double const* yDot, active* res, WithoutParamSensitivity) { - return residualImpl(t, secIdx, y, yDot, res, linalg::BandMatrix::RowIterator()); + return residualImpl(model, t, secIdx, y, yDot, res, linalg::BandMatrix::RowIterator()); } -int RadialConvectionDispersionOperatorBase::residual(double t, unsigned int secIdx, double const* y, double const* yDot, active* res, linalg::BandMatrix& jac) +int RadialConvectionDispersionOperatorBase::residual(const IModel& model, double t, unsigned int secIdx, double const* y, double const* yDot, active* res, linalg::BandMatrix& jac) { // Reset Jacobian jac.setAll(0.0); - return residualImpl(t, secIdx, y, yDot, res, jac.row(0)); + return residualImpl(model, t, secIdx, y, yDot, res, jac.row(0)); } -int RadialConvectionDispersionOperatorBase::residual(double t, unsigned int secIdx, double const* y, double const* yDot, active* res, WithParamSensitivity) +int RadialConvectionDispersionOperatorBase::residual(const IModel& model, double t, unsigned int secIdx, double const* y, double const* yDot, active* res, WithParamSensitivity) { - return residualImpl(t, secIdx, y, yDot, res, linalg::BandMatrix::RowIterator()); + return residualImpl(model, t, secIdx, y, yDot, res, linalg::BandMatrix::RowIterator()); } -int RadialConvectionDispersionOperatorBase::residual(double t, unsigned int secIdx, active const* y, double const* yDot, active* res, WithParamSensitivity) +int RadialConvectionDispersionOperatorBase::residual(const IModel& model, double t, unsigned int secIdx, active const* y, double const* yDot, active* res, WithParamSensitivity) { - return residualImpl(t, secIdx, y, yDot, res, linalg::BandMatrix::RowIterator()); + return residualImpl(model, t, secIdx, y, yDot, res, linalg::BandMatrix::RowIterator()); } template -int RadialConvectionDispersionOperatorBase::residualImpl(double t, unsigned int secIdx, StateType const* y, double const* yDot, ResidualType* res, RowIteratorType jacBegin) +int RadialConvectionDispersionOperatorBase::residualImpl(const IModel& model, double t, unsigned int secIdx, StateType const* y, double const* yDot, ResidualType* res, RowIteratorType jacBegin) { const ParamType u = static_cast(_curVelocity); active const* const d_rad = getSectionDependentSlice(_colDispersion, _nComp, secIdx); @@ -743,7 +827,9 @@ int RadialConvectionDispersionOperatorBase::residualImpl(double t, unsigned int _nComp, _nCol, 0u, - _nComp + _nComp, + _dispersionDep, + model }; return convdisp::residualKernelRadial(SimulationTime{t, secIdx}, y, yDot, res, jacBegin, fp); @@ -829,6 +915,16 @@ double RadialConvectionDispersionOperatorBase::inletJacobianFactor() const CADET bool RadialConvectionDispersionOperatorBase::setParameter(const ParameterId& pId, double value) { + // Check if parameter is in parameter dependence of column dispersion coefficient + if (_dispersionDep) + { + if (_dispersionDep->hasParameter(pId)) + { + _dispersionDep->setParameter(pId, value); + return true; + } + } + // Check whether column radius has changed and update discretization if necessary if (pId.name == hashString("COL_RADIUS_INNER")) { @@ -875,6 +971,17 @@ bool RadialConvectionDispersionOperatorBase::setParameter(const ParameterId& pId bool RadialConvectionDispersionOperatorBase::setSensitiveParameterValue(const std::unordered_set& sensParams, const ParameterId& pId, double value) { + // Check if parameter is in parameter dependence of column dispersion coefficient + if (_dispersionDep) + { + active* const param = _dispersionDep->getParameter(pId); + if (param) + { + param->setValue(value); + return true; + } + } + // Check whether column radius has changed and update discretization if necessary if (pId.name == hashString("COL_RADIUS_INNER")) { @@ -927,6 +1034,17 @@ bool RadialConvectionDispersionOperatorBase::setSensitiveParameterValue(const st bool RadialConvectionDispersionOperatorBase::setSensitiveParameter(std::unordered_set& sensParams, const ParameterId& pId, unsigned int adDirection, double adValue) { + // Check if parameter is in parameter dependence of column dispersion coefficient + if (_dispersionDep) + { + active* const param = _dispersionDep->getParameter(pId); + if (param) + { + param->setADValue(adDirection, adValue); + return true; + } + } + // Check whether column radius has changed and update discretization if necessary if (pId.name == hashString("COL_RADIUS_INNER")) { @@ -1079,8 +1197,7 @@ bool ConvectionDispersionOperator::notifyDiscontinuousSectionTransitio const unsigned int lb = _baseOp.jacobianLowerBandwidth(); const unsigned int ub = _baseOp.jacobianUpperBandwidth(); - const double u = static_cast(_baseOp.currentVelocity()); - if (u >= 0.0) + if (_baseOp.forwardFlow()) { // Forwards flow @@ -1137,6 +1254,7 @@ void ConvectionDispersionOperator::setFlowRates(const active& in, cons /** * @brief Computes the residual of the transport equations + * @param [in] model Model that owns the operator * @param [in] t Current time point * @param [in] secIdx Index of the current section * @param [in] y Pointer to unit operation's state vector @@ -1146,33 +1264,33 @@ void ConvectionDispersionOperator::setFlowRates(const active& in, cons * @return @c 0 on success, @c -1 on non-recoverable error, and @c +1 on recoverable error */ template -int ConvectionDispersionOperator::residual(double t, unsigned int secIdx, double const* y, double const* yDot, double* res, bool wantJac, WithoutParamSensitivity) +int ConvectionDispersionOperator::residual(const IModel& model, double t, unsigned int secIdx, double const* y, double const* yDot, double* res, bool wantJac, WithoutParamSensitivity) { if (wantJac) - return _baseOp.residual(t, secIdx, y, yDot, res, _jacC); + return _baseOp.residual(model, t, secIdx, y, yDot, res, _jacC); else - return _baseOp.residual(t, secIdx, y, yDot, res, WithoutParamSensitivity()); + return _baseOp.residual(model, t, secIdx, y, yDot, res, WithoutParamSensitivity()); } template -int ConvectionDispersionOperator::residual(double t, unsigned int secIdx, active const* y, double const* yDot, active* res, bool wantJac, WithoutParamSensitivity) +int ConvectionDispersionOperator::residual(const IModel& model, double t, unsigned int secIdx, active const* y, double const* yDot, active* res, bool wantJac, WithoutParamSensitivity) { - return _baseOp.residual(t, secIdx, y, yDot, res, WithoutParamSensitivity()); + return _baseOp.residual(model, t, secIdx, y, yDot, res, WithoutParamSensitivity()); } template -int ConvectionDispersionOperator::residual(double t, unsigned int secIdx, active const* y, double const* yDot, active* res, bool wantJac, WithParamSensitivity) +int ConvectionDispersionOperator::residual(const IModel& model, double t, unsigned int secIdx, active const* y, double const* yDot, active* res, bool wantJac, WithParamSensitivity) { - return _baseOp.residual(t, secIdx, y, yDot, res, WithParamSensitivity()); + return _baseOp.residual(model, t, secIdx, y, yDot, res, WithParamSensitivity()); } template -int ConvectionDispersionOperator::residual(double t, unsigned int secIdx, double const* y, double const* yDot, active* res, bool wantJac, WithParamSensitivity) +int ConvectionDispersionOperator::residual(const IModel& model, double t, unsigned int secIdx, double const* y, double const* yDot, active* res, bool wantJac, WithParamSensitivity) { if (wantJac) - return _baseOp.residual(t, secIdx, y, yDot, res, _jacC); + return _baseOp.residual(model, t, secIdx, y, yDot, res, _jacC); else - return _baseOp.residual(t, secIdx, y, yDot, res, WithParamSensitivity()); + return _baseOp.residual(model, t, secIdx, y, yDot, res, WithParamSensitivity()); } /** diff --git a/src/libcadet/model/parts/ConvectionDispersionOperator.hpp b/src/libcadet/model/parts/ConvectionDispersionOperator.hpp index 4e3b778b7..41bd14860 100644 --- a/src/libcadet/model/parts/ConvectionDispersionOperator.hpp +++ b/src/libcadet/model/parts/ConvectionDispersionOperator.hpp @@ -37,10 +37,13 @@ class IParameterProvider; class IConfigHelper; struct AdJacobianParams; struct SimulationTime; +class IModel; namespace model { +class IParameterParameterDependence; + namespace parts { @@ -75,12 +78,12 @@ class AxialConvectionDispersionOperatorBase bool configure(UnitOpIdx unitOpIdx, IParameterProvider& paramProvider, std::unordered_map& parameters); bool notifyDiscontinuousSectionTransition(double t, unsigned int secIdx); - int residual(double t, unsigned int secIdx, double const* y, double const* yDot, double* res, linalg::BandMatrix& jac); - int residual(double t, unsigned int secIdx, double const* y, double const* yDot, active* res, linalg::BandMatrix& jac); - int residual(double t, unsigned int secIdx, double const* y, double const* yDot, double* res, WithoutParamSensitivity); - int residual(double t, unsigned int secIdx, double const* y, double const* yDot, active* res, WithParamSensitivity); - int residual(double t, unsigned int secIdx, active const* y, double const* yDot, active* res, WithParamSensitivity); - int residual(double t, unsigned int secIdx, active const* y, double const* yDot, active* res, WithoutParamSensitivity); + int residual(const IModel& model, double t, unsigned int secIdx, double const* y, double const* yDot, double* res, linalg::BandMatrix& jac); + int residual(const IModel& model, double t, unsigned int secIdx, double const* y, double const* yDot, active* res, linalg::BandMatrix& jac); + int residual(const IModel& model, double t, unsigned int secIdx, double const* y, double const* yDot, double* res, WithoutParamSensitivity); + int residual(const IModel& model, double t, unsigned int secIdx, double const* y, double const* yDot, active* res, WithParamSensitivity); + int residual(const IModel& model, double t, unsigned int secIdx, active const* y, double const* yDot, active* res, WithParamSensitivity); + int residual(const IModel& model, double t, unsigned int secIdx, active const* y, double const* yDot, active* res, WithoutParamSensitivity); void multiplyWithDerivativeJacobian(const SimulationTime& simTime, double const* sDot, double* ret) const; void addTimeDerivativeToJacobian(double alpha, linalg::FactorizableBandMatrix& jacDisc); @@ -112,13 +115,7 @@ class AxialConvectionDispersionOperatorBase protected: template - int residualImpl(double t, unsigned int secIdx, StateType const* y, double const* yDot, ResidualType* res, RowIteratorType jacBegin); - - template - int residualForwardsFlow(double t, unsigned int secIdx, StateType const* y, double const* yDot, ResidualType* res, RowIteratorType jacBegin); - - template - int residualBackwardsFlow(double t, unsigned int secIdx, StateType const* y, double const* yDot, ResidualType* res, RowIteratorType jacBegin); + int residualImpl(const IModel& model, double t, unsigned int secIdx, StateType const* y, double const* yDot, ResidualType* res, RowIteratorType jacBegin); unsigned int _nComp; //!< Number of components unsigned int _nCol; //!< Number of axial cells @@ -140,6 +137,8 @@ class AxialConvectionDispersionOperatorBase bool _dispersionCompIndep; //!< Determines whether dispersion is component independent + IParameterParameterDependence* _dispersionDep; + // Indexer functionality // Strides @@ -182,12 +181,12 @@ class RadialConvectionDispersionOperatorBase bool configure(UnitOpIdx unitOpIdx, IParameterProvider& paramProvider, std::unordered_map& parameters); bool notifyDiscontinuousSectionTransition(double t, unsigned int secIdx); - int residual(double t, unsigned int secIdx, double const* y, double const* yDot, double* res, linalg::BandMatrix& jac); - int residual(double t, unsigned int secIdx, double const* y, double const* yDot, active* res, linalg::BandMatrix& jac); - int residual(double t, unsigned int secIdx, double const* y, double const* yDot, double* res, WithoutParamSensitivity); - int residual(double t, unsigned int secIdx, double const* y, double const* yDot, active* res, WithParamSensitivity); - int residual(double t, unsigned int secIdx, active const* y, double const* yDot, active* res, WithParamSensitivity); - int residual(double t, unsigned int secIdx, active const* y, double const* yDot, active* res, WithoutParamSensitivity); + int residual(const IModel& model, double t, unsigned int secIdx, double const* y, double const* yDot, double* res, linalg::BandMatrix& jac); + int residual(const IModel& model, double t, unsigned int secIdx, double const* y, double const* yDot, active* res, linalg::BandMatrix& jac); + int residual(const IModel& model, double t, unsigned int secIdx, double const* y, double const* yDot, double* res, WithoutParamSensitivity); + int residual(const IModel& model, double t, unsigned int secIdx, double const* y, double const* yDot, active* res, WithParamSensitivity); + int residual(const IModel& model, double t, unsigned int secIdx, active const* y, double const* yDot, active* res, WithParamSensitivity); + int residual(const IModel& model, double t, unsigned int secIdx, active const* y, double const* yDot, active* res, WithoutParamSensitivity); void multiplyWithDerivativeJacobian(const SimulationTime& simTime, double const* sDot, double* ret) const; void addTimeDerivativeToJacobian(double alpha, linalg::FactorizableBandMatrix& jacDisc); @@ -195,7 +194,7 @@ class RadialConvectionDispersionOperatorBase inline const active& columnLength() const CADET_NOEXCEPT { return _colLength; } inline const active& innerRadius() const CADET_NOEXCEPT { return _innerRadius; } inline const active& outerRadius() const CADET_NOEXCEPT { return _outerRadius; } - inline const active& currentVelocity() const CADET_NOEXCEPT { return _curVelocity; } + active currentVelocity(double pos) const CADET_NOEXCEPT; inline bool forwardFlow() const CADET_NOEXCEPT { return _curVelocity >= 0.0; } inline double cellCenter(unsigned int idx) const CADET_NOEXCEPT { return static_cast(_cellCenters[idx]); } @@ -217,13 +216,7 @@ class RadialConvectionDispersionOperatorBase protected: template - int residualImpl(double t, unsigned int secIdx, StateType const* y, double const* yDot, ResidualType* res, RowIteratorType jacBegin); - - template - int residualForwardsFlow(double t, unsigned int secIdx, StateType const* y, double const* yDot, ResidualType* res, RowIteratorType jacBegin); - - template - int residualBackwardsFlow(double t, unsigned int secIdx, StateType const* y, double const* yDot, ResidualType* res, RowIteratorType jacBegin); + int residualImpl(const IModel& model, double t, unsigned int secIdx, StateType const* y, double const* yDot, ResidualType* res, RowIteratorType jacBegin); void equidistantCells(); @@ -248,6 +241,8 @@ class RadialConvectionDispersionOperatorBase bool _dispersionCompIndep; //!< Determines whether dispersion is component independent + IParameterParameterDependence* _dispersionDep; + // Grid info std::vector _cellCenters; std::vector _cellSizes; @@ -288,10 +283,10 @@ class ConvectionDispersionOperator bool configure(UnitOpIdx unitOpIdx, IParameterProvider& paramProvider, std::unordered_map& parameters); bool notifyDiscontinuousSectionTransition(double t, unsigned int secIdx, const AdJacobianParams& adJac); - int residual(double t, unsigned int secIdx, double const* y, double const* yDot, double* res, bool wantJac, WithoutParamSensitivity); - int residual(double t, unsigned int secIdx, active const* y, double const* yDot, active* res, bool wantJac, WithParamSensitivity); - int residual(double t, unsigned int secIdx, active const* y, double const* yDot, active* res, bool wantJac, WithoutParamSensitivity); - int residual(double t, unsigned int secIdx, double const* y, double const* yDot, active* res, bool wantJac, WithParamSensitivity); + int residual(const IModel& model, double t, unsigned int secIdx, double const* y, double const* yDot, double* res, bool wantJac, WithoutParamSensitivity); + int residual(const IModel& model, double t, unsigned int secIdx, active const* y, double const* yDot, active* res, bool wantJac, WithParamSensitivity); + int residual(const IModel& model, double t, unsigned int secIdx, active const* y, double const* yDot, active* res, bool wantJac, WithoutParamSensitivity); + int residual(const IModel& model, double t, unsigned int secIdx, double const* y, double const* yDot, active* res, bool wantJac, WithParamSensitivity); void prepareADvectors(const AdJacobianParams& adJac) const; void extractJacobianFromAD(active const* const adRes, unsigned int adDirOffset); @@ -307,8 +302,8 @@ class ConvectionDispersionOperator #endif inline const active& columnLength() const CADET_NOEXCEPT { return _baseOp.columnLength(); } - inline const active& currentVelocity() const CADET_NOEXCEPT { return _baseOp.currentVelocity(); } - inline bool forwardFlow() const CADET_NOEXCEPT { return _baseOp.currentVelocity() >= 0.0; } + inline active currentVelocity(double pos) const CADET_NOEXCEPT { return _baseOp.currentVelocity(pos); } + inline bool forwardFlow() const CADET_NOEXCEPT { return _baseOp.forwardFlow(); } inline double inletJacobianFactor() const CADET_NOEXCEPT { return _baseOp.inletJacobianFactor(); } inline double cellCenter(unsigned int idx) const CADET_NOEXCEPT { return _baseOp.cellCenter(idx); } diff --git a/src/libcadet/model/parts/RadialConvectionDispersionKernel.hpp b/src/libcadet/model/parts/RadialConvectionDispersionKernel.hpp index 0a2c19b93..33e8b5f23 100644 --- a/src/libcadet/model/parts/RadialConvectionDispersionKernel.hpp +++ b/src/libcadet/model/parts/RadialConvectionDispersionKernel.hpp @@ -24,10 +24,14 @@ #include "Stencil.hpp" #include "linalg/CompressedSparseMatrix.hpp" #include "SimulationTypes.hpp" +#include "model/ParameterDependence.hpp" +#include "model/UnitOperation.hpp" namespace cadet { +class IModel; + namespace model { @@ -51,6 +55,8 @@ struct RadialFlowParameters unsigned int nCol; unsigned int offsetToInlet; //!< Offset to the first component of the inlet DOFs in the local state vector unsigned int offsetToBulk; //!< Offset to the first component of the first bulk cell in the local state vector + IParameterParameterDependence* parDep; + const IModel& model; }; @@ -119,11 +125,13 @@ namespace impl // Right side, leave out if we're in the last cell (boundary condition) if (cadet_likely(col < p.nCol - 1)) { - resBulkComp[col * p.strideCell] -= d_rad * static_cast(p.cellBounds[col + 1]) / denom * (stencil[1] - stencil[0]) / (static_cast(p.cellCenters[col + 1]) - static_cast(p.cellCenters[col])); + const double relCoord = (static_cast(p.cellBounds[col+1]) - static_cast(p.cellBounds[0])) / (static_cast(p.cellBounds[p.nCol - 1]) - static_cast(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(p.u) / static_cast(p.cellBounds[col+1])); + resBulkComp[col * p.strideCell] -= d_rad_right * static_cast(p.cellBounds[col+1]) / denom * (stencil[1] - stencil[0]) / (static_cast(p.cellCenters[col+1]) - static_cast(p.cellCenters[col])); // Jacobian entries if (wantJac) { - const double val = static_cast(d_rad) * static_cast(p.cellBounds[col+1]) / static_cast(denom) / (static_cast(p.cellCenters[col+1]) - static_cast(p.cellCenters[col])); + const double val = static_cast(d_rad_right) * static_cast(p.cellBounds[col+1]) / static_cast(denom) / (static_cast(p.cellCenters[col+1]) - static_cast(p.cellCenters[col])); jac[0] += val; jac[p.strideCell] -= val; } @@ -132,11 +140,13 @@ namespace impl // Left side, leave out if we're in the first cell (boundary condition) if (cadet_likely(col > 0)) { - resBulkComp[col * p.strideCell] -= d_rad * static_cast(p.cellBounds[col]) / denom * (stencil[-1] - stencil[0]) / (static_cast(p.cellCenters[col - 1]) - static_cast(p.cellCenters[col])); + const double relCoord = (static_cast(p.cellBounds[col]) - static_cast(p.cellBounds[0])) / (static_cast(p.cellBounds[p.nCol - 1]) - static_cast(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(p.u) / static_cast(p.cellBounds[col])); + resBulkComp[col * p.strideCell] -= d_rad_left * static_cast(p.cellBounds[col]) / denom * (stencil[-1] - stencil[0]) / (static_cast(p.cellCenters[col]) - static_cast(p.cellCenters[col-1])); // Jacobian entries if (wantJac) { - const double val = static_cast(d_rad) * static_cast(p.cellBounds[col]) / static_cast(denom) / (static_cast(p.cellCenters[col - 1]) - static_cast(p.cellCenters[col])); + const double val = static_cast(d_rad_left) * static_cast(p.cellBounds[col]) / static_cast(denom) / (static_cast(p.cellCenters[col]) - static_cast(p.cellCenters[col-1])); jac[0] += val; jac[-p.strideCell] -= val; } @@ -273,11 +283,13 @@ namespace impl // Right side, leave out if we're in the first cell (boundary condition) if (cadet_likely(col < p.nCol - 1)) { - resBulkComp[col * p.strideCell] -= d_rad * static_cast(p.cellBounds[col+1]) / denom * (stencil[-1] - stencil[0]) / (static_cast(p.cellCenters[col+1]) - static_cast(p.cellCenters[col])); + const double relCoord = (static_cast(p.cellBounds[col+1]) - static_cast(p.cellBounds[0])) / (static_cast(p.cellBounds[p.nCol - 1]) - static_cast(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(p.u) / static_cast(p.cellBounds[col+1])); + resBulkComp[col * p.strideCell] -= d_rad_right * static_cast(p.cellBounds[col+1]) / denom * (stencil[-1] - stencil[0]) / (static_cast(p.cellCenters[col+1]) - static_cast(p.cellCenters[col])); // Jacobian entries if (wantJac) { - const double val = static_cast(d_rad) * static_cast(p.cellBounds[col+1]) / static_cast(denom) / (static_cast(p.cellCenters[col+1]) - static_cast(p.cellCenters[col])); + const double val = static_cast(d_rad_right) * static_cast(p.cellBounds[col+1]) / static_cast(denom) / (static_cast(p.cellCenters[col+1]) - static_cast(p.cellCenters[col])); jac[0] += val; jac[p.strideCell] -= val; } @@ -286,11 +298,13 @@ namespace impl // Left side, leave out if we're in the last cell (boundary condition) if (cadet_likely(col > 0)) { - resBulkComp[col * p.strideCell] -= d_rad * static_cast(p.cellBounds[col]) / denom * (stencil[1] - stencil[0]) / (static_cast(p.cellCenters[col-1]) - static_cast(p.cellCenters[col])); + const double relCoord = (static_cast(p.cellBounds[col]) - static_cast(p.cellBounds[0])) / (static_cast(p.cellBounds[p.nCol - 1]) - static_cast(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(p.u) / static_cast(p.cellBounds[col])); + resBulkComp[col * p.strideCell] -= d_rad_left * static_cast(p.cellBounds[col]) / denom * (stencil[1] - stencil[0]) / (static_cast(p.cellCenters[col-1]) - static_cast(p.cellCenters[col])); // Jacobian entries if (wantJac) { - const double val = static_cast(d_rad) * static_cast(p.cellBounds[col]) / static_cast(denom) / (static_cast(p.cellCenters[col-1]) - static_cast(p.cellCenters[col])); + const double val = static_cast(d_rad_left) * static_cast(p.cellBounds[col]) / static_cast(denom) / (static_cast(p.cellCenters[col-1]) - static_cast(p.cellCenters[col])); jac[0] += val; jac[-p.strideCell] -= val; } diff --git a/src/libcadet/model/parts/TwoDimensionalConvectionDispersionOperator.cpp b/src/libcadet/model/parts/TwoDimensionalConvectionDispersionOperator.cpp index b2abc3c1b..adefb8726 100644 --- a/src/libcadet/model/parts/TwoDimensionalConvectionDispersionOperator.cpp +++ b/src/libcadet/model/parts/TwoDimensionalConvectionDispersionOperator.cpp @@ -20,6 +20,9 @@ #include "SimulationTypes.hpp" #include "linalg/CompressedSparseMatrix.hpp" #include "model/parts/AxialConvectionDispersionKernel.hpp" +#include "model/ParameterDependence.hpp" +#include "ConfigurationHelper.hpp" + #ifdef SUPERLU_FOUND #include "linalg/SuperLUSparseMatrix.hpp" @@ -654,7 +657,7 @@ class TwoDimensionalConvectionDispersionOperator::DenseDirectSolver : public Two * @brief Creates a TwoDimensionalConvectionDispersionOperator */ TwoDimensionalConvectionDispersionOperator::TwoDimensionalConvectionDispersionOperator() : _colPorosities(0), _dir(0), _stencilMemory(sizeof(active) * Weno::maxStencilSize()), - _wenoDerivatives(new double[Weno::maxStencilSize()]), _weno(), _linearSolver(nullptr) + _wenoDerivatives(new double[Weno::maxStencilSize()]), _weno(), _linearSolver(nullptr), _dispersionDep(nullptr) { } @@ -662,6 +665,7 @@ TwoDimensionalConvectionDispersionOperator::~TwoDimensionalConvectionDispersionO { delete[] _wenoDerivatives; delete _linearSolver; + delete _dispersionDep; } /** @@ -673,13 +677,16 @@ TwoDimensionalConvectionDispersionOperator::~TwoDimensionalConvectionDispersionO * @param [in] dynamicReactions Determines whether the sparsity pattern accounts for dynamic reactions * @return @c true if configuration went fine, @c false otherwise */ -bool TwoDimensionalConvectionDispersionOperator::configureModelDiscretization(IParameterProvider& paramProvider, unsigned int nComp, unsigned int nCol, unsigned int nRad, bool dynamicReactions) +bool TwoDimensionalConvectionDispersionOperator::configureModelDiscretization(IParameterProvider& paramProvider, const IConfigHelper& helper, unsigned int nComp, unsigned int nCol, unsigned int nRad, bool dynamicReactions) { _nComp = nComp; _nCol = nCol; _nRad = nRad; _hasDynamicReactions = dynamicReactions; + // TODO: Add support for parameter dependent dispersion + _dispersionDep = helper.createParameterParameterDependence("CONSTANT_ONE"); + paramProvider.pushScope("discretization"); // Read WENO settings and apply them @@ -944,40 +951,40 @@ const active& TwoDimensionalConvectionDispersionOperator::radialDispersion(unsig * @param [in] wantJac Determines whether the Jacobian is computed or not * @return @c 0 on success, @c -1 on non-recoverable error, and @c +1 on recoverable error */ -int TwoDimensionalConvectionDispersionOperator::residual(double t, unsigned int secIdx, double const* y, double const* yDot, double* res, bool wantJac, WithoutParamSensitivity) +int TwoDimensionalConvectionDispersionOperator::residual(const IModel& model, double t, unsigned int secIdx, double const* y, double const* yDot, double* res, bool wantJac, WithoutParamSensitivity) { if (wantJac) - return residualImpl(t, secIdx, y, yDot, res); + return residualImpl(model, t, secIdx, y, yDot, res); else - return residualImpl(t, secIdx, y, yDot, res); + return residualImpl(model, t, secIdx, y, yDot, res); } -int TwoDimensionalConvectionDispersionOperator::residual(double t, unsigned int secIdx, active const* y, double const* yDot, active* res, bool wantJac, WithoutParamSensitivity) +int TwoDimensionalConvectionDispersionOperator::residual(const IModel& model, double t, unsigned int secIdx, active const* y, double const* yDot, active* res, bool wantJac, WithoutParamSensitivity) { if (wantJac) - return residualImpl(t, secIdx, y, yDot, res); + return residualImpl(model, t, secIdx, y, yDot, res); else - return residualImpl(t, secIdx, y, yDot, res); + return residualImpl(model, t, secIdx, y, yDot, res); } -int TwoDimensionalConvectionDispersionOperator::residual(double t, unsigned int secIdx, double const* y, double const* yDot, active* res, bool wantJac, WithParamSensitivity) +int TwoDimensionalConvectionDispersionOperator::residual(const IModel& model, double t, unsigned int secIdx, double const* y, double const* yDot, active* res, bool wantJac, WithParamSensitivity) { if (wantJac) - return residualImpl(t, secIdx, y, yDot, res); + return residualImpl(model, t, secIdx, y, yDot, res); else - return residualImpl(t, secIdx, y, yDot, res); + return residualImpl(model, t, secIdx, y, yDot, res); } -int TwoDimensionalConvectionDispersionOperator::residual(double t, unsigned int secIdx, active const* y, double const* yDot, active* res, bool wantJac, WithParamSensitivity) +int TwoDimensionalConvectionDispersionOperator::residual(const IModel& model, double t, unsigned int secIdx, active const* y, double const* yDot, active* res, bool wantJac, WithParamSensitivity) { if (wantJac) - return residualImpl(t, secIdx, y, yDot, res); + return residualImpl(model, t, secIdx, y, yDot, res); else - return residualImpl(t, secIdx, y, yDot, res); + return residualImpl(model, t, secIdx, y, yDot, res); } template -int TwoDimensionalConvectionDispersionOperator::residualImpl(double t, unsigned int secIdx, StateType const* y, double const* yDot, ResidualType* res) +int TwoDimensionalConvectionDispersionOperator::residualImpl(const IModel& model, double t, unsigned int secIdx, StateType const* y, double const* yDot, ResidualType* res) { if (wantJac) { @@ -1003,7 +1010,9 @@ int TwoDimensionalConvectionDispersionOperator::residualImpl(double t, unsigned _nComp, _nCol, _nComp * i, // Offset to the first component of the inlet DOFs in the local state vector - _nComp * (_nRad + i) // Offset to the first component of the first bulk cell in the local state vector + _nComp * (_nRad + i), // Offset to the first component of the first bulk cell in the local state vector + _dispersionDep, + model }; if (wantJac) diff --git a/src/libcadet/model/parts/TwoDimensionalConvectionDispersionOperator.hpp b/src/libcadet/model/parts/TwoDimensionalConvectionDispersionOperator.hpp index b90db84d3..960cfa206 100644 --- a/src/libcadet/model/parts/TwoDimensionalConvectionDispersionOperator.hpp +++ b/src/libcadet/model/parts/TwoDimensionalConvectionDispersionOperator.hpp @@ -34,10 +34,14 @@ namespace cadet { class IParameterProvider; +class IModel; +class IConfigHelper; namespace model { +class IParameterParameterDependence; + namespace parts { @@ -70,14 +74,14 @@ class TwoDimensionalConvectionDispersionOperator void setFlowRates(int compartment, const active& in, const active& out) CADET_NOEXCEPT; void setFlowRates(active const* in, active const* out) CADET_NOEXCEPT; - bool configureModelDiscretization(IParameterProvider& paramProvider, unsigned int nComp, unsigned int nCol, unsigned int nRad, bool dynamicReactions); + bool configureModelDiscretization(IParameterProvider& paramProvider, const IConfigHelper& helper, unsigned int nComp, unsigned int nCol, unsigned int nRad, bool dynamicReactions); bool configure(UnitOpIdx unitOpIdx, IParameterProvider& paramProvider, std::unordered_map& parameters); bool notifyDiscontinuousSectionTransition(double t, unsigned int secIdx); - int residual(double t, unsigned int secIdx, double const* y, double const* yDot, double* res, bool wantJac, WithoutParamSensitivity); - int residual(double t, unsigned int secIdx, active const* y, double const* yDot, active* res, bool wantJac, WithoutParamSensitivity); - int residual(double t, unsigned int secIdx, active const* y, double const* yDot, active* res, bool wantJac, WithParamSensitivity); - int residual(double t, unsigned int secIdx, double const* y, double const* yDot, active* res, bool wantJac, WithParamSensitivity); + int residual(const IModel& model, double t, unsigned int secIdx, double const* y, double const* yDot, double* res, bool wantJac, WithoutParamSensitivity); + int residual(const IModel& model, double t, unsigned int secIdx, active const* y, double const* yDot, active* res, bool wantJac, WithoutParamSensitivity); + int residual(const IModel& model, double t, unsigned int secIdx, active const* y, double const* yDot, active* res, bool wantJac, WithParamSensitivity); + int residual(const IModel& model, double t, unsigned int secIdx, double const* y, double const* yDot, active* res, bool wantJac, WithParamSensitivity); bool solveTimeDerivativeSystem(const SimulationTime& simTime, double* const rhs); void multiplyWithDerivativeJacobian(const SimulationTime& simTime, double const* sDot, double* ret) const; @@ -118,7 +122,7 @@ class TwoDimensionalConvectionDispersionOperator void assembleDiscretizedJacobian(double alpha); template - int residualImpl(double t, unsigned int secIdx, StateType const* y, double const* yDot, ResidualType* res); + int residualImpl(const IModel& model, double t, unsigned int secIdx, StateType const* y, double const* yDot, ResidualType* res); void setSparsityPattern(); @@ -177,6 +181,8 @@ class TwoDimensionalConvectionDispersionOperator linalg::CompressedSparseMatrix _jacC; //!< Jacobian LinearSolver* _linearSolver; //!< Solves linear system with time discretized Jacobian + + IParameterParameterDependence* _dispersionDep; }; } // namespace parts diff --git a/src/tools/ToolsHelper.hpp b/src/tools/ToolsHelper.hpp index 153a95ac1..f8d83e38f 100644 --- a/src/tools/ToolsHelper.hpp +++ b/src/tools/ToolsHelper.hpp @@ -75,6 +75,14 @@ inline void parseUnitType(std::string& unitType) unitType = "GENERAL_RATE_MODEL_2D"; } +inline void parseUnitType(std::string& unitType, bool radialFlow) +{ + parseUnitType(unitType); + + if (radialFlow) + unitType = "RADIAL_" + unitType; +} + inline void addSensitivitiyParserToCmdLine(TCLAP::CmdLine& cmd, std::vector& sensitivities) { cmd >> (new TCLAP::MultiArg("S", "sens", diff --git a/src/tools/createLWE.cpp b/src/tools/createLWE.cpp index adb541efb..43baf037f 100644 --- a/src/tools/createLWE.cpp +++ b/src/tools/createLWE.cpp @@ -33,6 +33,8 @@ struct ProgramOptions double endTime; double constAlg; double stddevAlg; + bool radialFlow; + bool velocityDependence; bool reverseFlow; bool adJacobian; int nPar; @@ -64,6 +66,8 @@ int main(int argc, char** argv) cmd >> (new TCLAP::ValueArg("c", "constAlg", "Set all algebraic variables to constant value", false, nanVal, "Value"))->storeIn(&opts.constAlg); cmd >> (new TCLAP::ValueArg("s", "stddevAlg", "Perturb algebraic variables with normal variates", false, nanVal, "Value"))->storeIn(&opts.stddevAlg); cmd >> (new TCLAP::SwitchArg("", "reverseFlow", "Reverse the flow for column"))->storeIn(&opts.reverseFlow); + cmd >> (new TCLAP::SwitchArg("", "radialFlow", "Use radial flow column"))->storeIn(&opts.radialFlow); + cmd >> (new TCLAP::SwitchArg("", "velDep", "Use velocity dependent dispersion and film diffusion"))->storeIn(&opts.velocityDependence); cmd >> (new TCLAP::ValueArg("", "rad", "Number of radial cells (default: 3)", false, 3, "Value"))->storeIn(&opts.nRad); cmd >> (new TCLAP::ValueArg("", "parTypes", "Number of particle types (default: 1)", false, 1, "Value"))->storeIn(&opts.nParType); addMiscToCmdLine(cmd, opts); @@ -83,7 +87,7 @@ int main(int argc, char** argv) writer.openFile(opts.fileName, "co"); writer.pushGroup("input"); - parseUnitType(opts.unitType); + parseUnitType(opts.unitType, opts.radialFlow); const bool isGRM2D = (opts.unitType == "GENERAL_RATE_MODEL_2D"); // Model @@ -100,9 +104,9 @@ int main(int argc, char** argv) // Transport if (!opts.reverseFlow) - writer.scalar("VELOCITY", 5.75e-4); + writer.scalar("VELOCITY", 1.0); else - writer.scalar("VELOCITY", -5.75e-4); + writer.scalar("VELOCITY", -1.0); writer.scalar("COL_DISPERSION", 5.75e-8); writer.scalar("COL_DISPERSION_RADIAL", 1e-6); @@ -135,9 +139,26 @@ int main(int argc, char** argv) writer.vector("PAR_SURFDIFFUSION", 4, parSurfDiff); } + if (opts.velocityDependence) + { + writer.scalar("COL_DISPERSION_DEP", "POWER_LAW"); + writer.scalar("COL_DISPERSION_DEP_BASE", 1.25); + writer.scalar("COL_DISPERSION_DEP_EXPONENT", 1.0); + + writer.scalar("FILM_DIFFUSION_DEP", "POWER_LAW"); + writer.scalar("FILM_DIFFUSION_DEP_BASE", 1.25); + writer.scalar("FILM_DIFFUSION_DEP_EXPONENT", 1.0); + } + // Geometry - writer.scalar("COL_LENGTH", 0.014); + if (opts.radialFlow) + writer.scalar("COL_LENGTH", 0.0014); + else + writer.scalar("COL_LENGTH", 0.014); writer.scalar("COL_RADIUS", 0.01); + writer.scalar("COL_RADIUS_INNER", 0.01); + writer.scalar("COL_RADIUS_OUTER", 0.04); + writer.scalar("CROSS_SECTION_AREA", 0.0003141592653589793); writer.scalar("COL_POROSITY", 0.37); const double par_radius = 4.5e-5; const double par_coreradius = 0.0; @@ -168,7 +189,7 @@ int main(int argc, char** argv) } else { - writer.scalar("PAR_GEOM", std::string("SPHERE")); + writer.scalar("PAR_GEOM", std::string("SPHERE")); writer.scalar("PAR_RADIUS", par_radius); writer.scalar("PAR_CORERADIUS", par_coreradius); writer.scalar("PAR_POROSITY", par_porosity); @@ -447,12 +468,12 @@ int main(int argc, char** argv) { // Connection list is 1x7 since we have 1 connection between // the two unit operations (and we need to have 7 columns) - const double connMatrix[] = {1, 0, -1, -1, -1, -1, 1.0}; + const double connMatrix[] = {1, 0, -1, -1, -1, -1, 6.683738370512285e-8}; // Connections: From unit operation 1 port -1 (i.e., all ports) // to unit operation 0 port -1 (i.e., all ports), // connect component -1 (i.e., all components) // to component -1 (i.e., all components) with - // a flow rate of 1.0 + // a flow rate of 6.683738370512285e-8 m^3/s writer.vector("CONNECTIONS", 7, connMatrix); } diff --git a/test/ConvectionDispersionOperator.cpp b/test/ConvectionDispersionOperator.cpp index 0907603a4..475c243a1 100644 --- a/test/ConvectionDispersionOperator.cpp +++ b/test/ConvectionDispersionOperator.cpp @@ -21,14 +21,17 @@ #include "AdUtils.hpp" #include "SimulationTypes.hpp" #include "ModelBuilderImpl.hpp" +#include "ParameterDependenceFactory.hpp" #include "JsonTestModels.hpp" #include "ColumnTests.hpp" #include "JacobianHelper.hpp" #include "Utils.hpp" +#include "Dummies.hpp" #include #include +#include namespace { @@ -114,7 +117,8 @@ namespace } } - inline cadet::active* createAndConfigureOperator(cadet::model::parts::AxialConvectionDispersionOperator& convDispOp, int& nComp, int& nCol, int wenoOrder) + template + inline cadet::active* createAndConfigureOperator(OperatorType& convDispOp, int& nComp, int& nCol, int wenoOrder) { // Obtain parameters from some test case cadet::JsonParameterProvider jpp = createColumnWithSMA("GENERAL_RATE_MODEL"); @@ -142,13 +146,14 @@ namespace } // namespace +template void testResidualBulkWenoForwardBackward(int wenoOrder) { SECTION("Forward vs backward flow residual bulk (WENO=" + std::to_string(wenoOrder) + ")") { int nComp = 0; int nCol = 0; - cadet::model::parts::AxialConvectionDispersionOperator convDispOp; + OperatorType convDispOp; cadet::active* const velocity = createAndConfigureOperator(convDispOp, nComp, nCol, wenoOrder); const double origVelocity = velocity->getValue(); @@ -168,7 +173,7 @@ void testResidualBulkWenoForwardBackward(int wenoOrder) SECTION("Forward flow yields backwards flow residual (zero state)") { // Forward flow residual - convDispOp.residual(0.0, 0u, y.data(), nullptr, res.data(), false, cadet::WithoutParamSensitivity()); + convDispOp.residual(DummyModel(), 0.0, 0u, y.data(), nullptr, res.data(), false, cadet::WithoutParamSensitivity()); // Reverse flow velocity->setValue(-origVelocity); @@ -176,7 +181,7 @@ void testResidualBulkWenoForwardBackward(int wenoOrder) std::vector resRev(nDof, 0.0); // Backward flow residual - convDispOp.residual(0.0, 0u, y.data(), nullptr, resRev.data(), false, cadet::WithoutParamSensitivity()); + convDispOp.residual(DummyModel(), 0.0, 0u, y.data(), nullptr, resRev.data(), false, cadet::WithoutParamSensitivity()); // Compare compareResidualBulkFwdBwd(res.data(), resRev.data(), nComp, nCol); @@ -187,7 +192,7 @@ void testResidualBulkWenoForwardBackward(int wenoOrder) std::vector resFwd2(nDof, 0.0); // Forward flow residual - convDispOp.residual(0.0, 0u, y.data(), nullptr, resFwd2.data(), false, cadet::WithoutParamSensitivity()); + convDispOp.residual(DummyModel(), 0.0, 0u, y.data(), nullptr, resFwd2.data(), false, cadet::WithoutParamSensitivity()); // Compare against first forward flow residual compareResidualBulkFwdFwd(res.data(), resFwd2.data(), nComp, nCol); @@ -199,7 +204,7 @@ void testResidualBulkWenoForwardBackward(int wenoOrder) // Fill state vector with some values fillStateBulkFwd(y.data(), [](unsigned int comp, unsigned int col, unsigned int idx) { return std::abs(std::sin(idx * 0.13)); }, nComp, nCol); - convDispOp.residual(0.0, 0u, y.data(), nullptr, res.data(), false, cadet::WithoutParamSensitivity()); + convDispOp.residual(DummyModel(), 0.0, 0u, y.data(), nullptr, res.data(), false, cadet::WithoutParamSensitivity()); // Reverse state for backwards flow fillStateBulkBwd(y.data(), [](unsigned int comp, unsigned int col, unsigned int idx) { return std::abs(std::sin(idx * 0.13)); }, nComp, nCol); @@ -210,7 +215,7 @@ void testResidualBulkWenoForwardBackward(int wenoOrder) std::vector resRev(nDof, 0.0); // Backward flow residual - convDispOp.residual(0.0, 0u, y.data(), nullptr, resRev.data(), false, cadet::WithoutParamSensitivity()); + convDispOp.residual(DummyModel(), 0.0, 0u, y.data(), nullptr, resRev.data(), false, cadet::WithoutParamSensitivity()); // Compare compareResidualBulkFwdBwd(res.data(), resRev.data(), nComp, nCol); @@ -224,7 +229,7 @@ void testResidualBulkWenoForwardBackward(int wenoOrder) fillStateBulkFwd(y.data(), [](unsigned int comp, unsigned int col, unsigned int idx) { return std::abs(std::sin(idx * 0.13)); }, nComp, nCol); // Forward flow residual - convDispOp.residual(0.0, 0u, y.data(), nullptr, resFwd2.data(), false, cadet::WithoutParamSensitivity()); + convDispOp.residual(DummyModel(), 0.0, 0u, y.data(), nullptr, resFwd2.data(), false, cadet::WithoutParamSensitivity()); // Compare against first forward flow residual compareResidualBulkFwdFwd(res.data(), resFwd2.data(), nComp, nCol); @@ -232,11 +237,12 @@ void testResidualBulkWenoForwardBackward(int wenoOrder) } } +template void testTimeDerivativeBulkJacobianFD(double h, double absTol, double relTol) { int nComp = 0; int nCol = 0; - cadet::model::parts::AxialConvectionDispersionOperator convDispOp; + OperatorType convDispOp; createAndConfigureOperator(convDispOp, nComp, nCol, cadet::Weno::maxOrder()); // Setup matrices @@ -256,19 +262,20 @@ void testTimeDerivativeBulkJacobianFD(double h, double absTol, double relTol) // Compare Jacobians cadet::test::compareJacobianFD( - [&](double const* dir, double* res) -> void { convDispOp.residual(0.0, 0u, y.data(), dir, res, false, cadet::WithoutParamSensitivity()); }, + [&](double const* dir, double* res) -> void { convDispOp.residual(DummyModel(), 0.0, 0u, y.data(), dir, res, false, cadet::WithoutParamSensitivity()); }, [&](double const* dir, double* res) -> void { convDispOp.multiplyWithDerivativeJacobian(cadet::SimulationTime{0.0, 0u}, dir, res); }, yDot.data(), jacDir.data(), jacCol1.data(), jacCol2.data(), nDof, h, absTol, relTol); } +template void testBulkJacobianWenoForwardBackward(int wenoOrder) { SECTION("Forward vs backward flow Jacobian (WENO=" + std::to_string(wenoOrder) + ")") { int nComp = 0; int nCol = 0; - cadet::model::parts::AxialConvectionDispersionOperator opAna; - cadet::model::parts::AxialConvectionDispersionOperator opAD; + OperatorType opAna; + OperatorType opAD; cadet::active* const anaVelocity = createAndConfigureOperator(opAna, nComp, nCol, wenoOrder); cadet::active* const adVelocity = createAndConfigureOperator(opAD, nComp, nCol, wenoOrder); @@ -296,17 +303,17 @@ void testBulkJacobianWenoForwardBackward(int wenoOrder) SECTION("Forward then backward flow (nonzero state)") { // Compute state Jacobian - opAna.residual(0.0, 0u, y.data(), nullptr, jacDir.data(), true, cadet::WithoutParamSensitivity()); + opAna.residual(DummyModel(), 0.0, 0u, y.data(), nullptr, jacDir.data(), true, cadet::WithoutParamSensitivity()); std::fill(jacDir.begin(), jacDir.end(), 0.0); cadet::ad::copyToAd(y.data(), adY, nDof); cadet::ad::resetAd(adRes, nDof); - opAD.residual(0.0, 0u, adY, nullptr, adRes, false, cadet::WithoutParamSensitivity()); + opAD.residual(DummyModel(), 0.0, 0u, adY, nullptr, adRes, false, cadet::WithoutParamSensitivity()); opAD.extractJacobianFromAD(adRes, 0); const std::function anaResidual = [&](double const* lDir, double* res) -> void { - opAna.residual(0.0, 0u, lDir - nComp, nullptr, res - nComp, false, cadet::WithoutParamSensitivity()); + opAna.residual(DummyModel(), 0.0, 0u, lDir - nComp, nullptr, res - nComp, false, cadet::WithoutParamSensitivity()); }; const std::function anaMultJac = [&](double const* lDir, double* res) -> void @@ -338,12 +345,12 @@ void testBulkJacobianWenoForwardBackward(int wenoOrder) opAD.notifyDiscontinuousSectionTransition(0.0, 0u, cadet::AdJacobianParams{adRes, adY, 0u}); // Compute state Jacobian - opAna.residual(0.0, 0u, y.data(), nullptr, jacDir.data(), true, cadet::WithoutParamSensitivity()); + opAna.residual(DummyModel(), 0.0, 0u, y.data(), nullptr, jacDir.data(), true, cadet::WithoutParamSensitivity()); std::fill(jacDir.begin(), jacDir.end(), 0.0); cadet::ad::copyToAd(y.data(), adY, nDof); cadet::ad::resetAd(adRes, nDof); - opAD.residual(0.0, 0u, adY, nullptr, adRes, false, cadet::WithoutParamSensitivity()); + opAD.residual(DummyModel(), 0.0, 0u, adY, nullptr, adRes, false, cadet::WithoutParamSensitivity()); opAD.extractJacobianFromAD(adRes, 0); // Compare Jacobians @@ -362,6 +369,121 @@ void testBulkJacobianWenoForwardBackward(int wenoOrder) } } + +struct AxialFlow +{ + typedef cadet::model::parts::convdisp::AxialFlowParameters Params; + + static void sparsityPattern(cadet::linalg::SparsityPatternRowIterator itBegin, unsigned int nComp, unsigned int nCol, int strideCell, double u, cadet::Weno& weno) + { + cadet::model::parts::convdisp::sparsityPatternAxial(itBegin, nComp, nCol, strideCell, u, weno); + } + + static void residual(double const* y, double const* yDot, double* res, const Params& fp) + { + cadet::model::parts::convdisp::residualKernelAxial(cadet::SimulationTime{0.0, 0u}, y, yDot, res, cadet::linalg::BandedSparseRowIterator(), fp); + } + + template + static void residualWithJacobian(double const* y, double const* yDot, double* res, IteratorType jacBegin, const Params& fp) + { + cadet::model::parts::convdisp::residualKernelAxial(cadet::SimulationTime{0.0, 0u}, y, yDot, res, jacBegin, fp); + } + + std::unique_ptr parDep; + + AxialFlow() + { + cadet::ParameterDependenceFactory paramDepFactory; + parDep.reset(paramDepFactory.createParameterDependence("CONSTANT_ONE")); + } + + Params makeParams(double u, cadet::active const* d_ax, double h, double* wenoDerivatives, cadet::Weno* weno, cadet::ArrayPool* stencilMemory, int strideCell, int nComp, int nCol) + { + return Params { + u, + d_ax, + h, + wenoDerivatives, + weno, + stencilMemory, + 1e-12, + strideCell, + static_cast(nComp), + static_cast(nCol), + 0u, + static_cast(nComp), + parDep.get(), + DummyModel() + }; + } +}; + +struct RadialFlow +{ + typedef cadet::model::parts::convdisp::RadialFlowParameters Params; + + static void sparsityPattern(cadet::linalg::SparsityPatternRowIterator itBegin, unsigned int nComp, unsigned int nCol, int strideCell, double u, cadet::Weno& weno) + { + cadet::model::parts::convdisp::sparsityPatternRadial(itBegin, nComp, nCol, strideCell, u, weno); + } + + static void residual(double const* y, double const* yDot, double* res, const Params& fp) + { + cadet::model::parts::convdisp::residualKernelRadial(cadet::SimulationTime{0.0, 0u}, y, yDot, res, cadet::linalg::BandedSparseRowIterator(), fp); + } + + template + static void residualWithJacobian(double const* y, double const* yDot, double* res, IteratorType jacBegin, const Params& fp) + { + cadet::model::parts::convdisp::residualKernelRadial(cadet::SimulationTime{0.0, 0u}, y, yDot, res, jacBegin, fp); + } + + std::unique_ptr parDep; + std::vector centers; + std::vector sizes; + std::vector bounds; + + RadialFlow() + { + cadet::ParameterDependenceFactory paramDepFactory; + parDep.reset(paramDepFactory.createParameterDependence("CONSTANT_ONE")); + } + + Params makeParams(double u, cadet::active const* d_rad, double h, double* wenoDerivatives, cadet::Weno* weno, cadet::ArrayPool* stencilMemory, int strideCell, int nComp, int nCol) + { + centers.resize(10); + sizes.resize(10); + bounds.resize(11); + + for (int i = 0; i < 10; ++i) + { + centers[i] = (i + 1.5) * h; + sizes[i] = h; + bounds[i] = (i + 1) * h; + } + bounds.back() = 11 * h; + + return Params { + u, + d_rad, + centers.data(), + sizes.data(), + bounds.data(), + stencilMemory, + strideCell, + static_cast(nComp), + static_cast(nCol), + 0u, + static_cast(nComp), + parDep.get(), + DummyModel() + }; + } +}; + + +template void testBulkJacobianSparsityWeno(int wenoOrder, bool forwardFlow) { SECTION("WENO=" + std::to_string(wenoOrder)) @@ -381,24 +503,25 @@ void testBulkJacobianSparsityWeno(int wenoOrder, bool forwardFlow) weno.order(wenoOrder); weno.boundaryTreatment(cadet::Weno::BoundaryTreatment::ReduceOrder); - cadet::model::parts::convdisp::AxialFlowParameters fp{ + cadet::ParameterDependenceFactory paramDepFactory; + std::unique_ptr parDep(paramDepFactory.createParameterDependence("CONSTANT_ONE")); + + FlowType ft; + typename FlowType::Params fp = ft.makeParams( u, d_c.data(), h, wenoDerivatives.data(), &weno, &stencilMemory, - 1e-12, strideCell, - static_cast(nComp), - static_cast(nCol), - 0u, - static_cast(nComp) - }; + nComp, + nCol + ); // Obtain sparsity pattern cadet::linalg::SparsityPattern pattern(nComp * nCol, std::max(weno.lowerBandwidth() + 1u, 1u) + 1u + std::max(weno.upperBandwidth(), 1u)); - cadet::model::parts::convdisp::sparsityPatternAxial(pattern.row(0), nComp, nCol, strideCell, u, weno); + FlowType::sparsityPattern(pattern.row(0), nComp, nCol, strideCell, u, weno); // Obtain memory for state, Jacobian columns const int nDof = nComp + nComp * nCol; @@ -415,10 +538,10 @@ void testBulkJacobianSparsityWeno(int wenoOrder, bool forwardFlow) // Central finite differences y[nComp + col] = ref * (1.0 + 1e-6); - cadet::model::parts::convdisp::residualKernelAxial(cadet::SimulationTime{0.0, 0u}, y.data(), nullptr, jacCol1.data(), cadet::linalg::BandedSparseRowIterator(), fp); + FlowType::residual(y.data(), nullptr, jacCol1.data(), fp); y[nComp + col] = ref * (1.0 - 1e-6); - cadet::model::parts::convdisp::residualKernelAxial(cadet::SimulationTime{0.0, 0u}, y.data(), nullptr, jacCol2.data(), cadet::linalg::BandedSparseRowIterator(), fp); + FlowType::residual(y.data(), nullptr, jacCol2.data(), fp); y[nComp + col] = ref; @@ -436,6 +559,7 @@ void testBulkJacobianSparsityWeno(int wenoOrder, bool forwardFlow) } } +template void testBulkJacobianSparseBandedWeno(int wenoOrder, bool forwardFlow) { SECTION("WENO=" + std::to_string(wenoOrder)) @@ -455,26 +579,27 @@ void testBulkJacobianSparseBandedWeno(int wenoOrder, bool forwardFlow) weno.order(wenoOrder); weno.boundaryTreatment(cadet::Weno::BoundaryTreatment::ReduceOrder); - cadet::model::parts::convdisp::AxialFlowParameters fp{ + cadet::ParameterDependenceFactory paramDepFactory; + std::unique_ptr parDep(paramDepFactory.createParameterDependence("CONSTANT_ONE")); + + FlowType ft; + typename FlowType::Params fp = ft.makeParams( u, d_c.data(), h, wenoDerivatives.data(), &weno, &stencilMemory, - 1e-12, strideCell, - static_cast(nComp), - static_cast(nCol), - 0u, - static_cast(nComp) - }; + nComp, + nCol + ); // Obtain sparsity pattern const unsigned int lowerBandwidth = std::max(weno.lowerBandwidth() + 1u, 1u); const unsigned int upperBandwidth = std::max(weno.upperBandwidth(), 1u); cadet::linalg::SparsityPattern pattern(nComp * nCol, lowerBandwidth + 1u + upperBandwidth); - cadet::model::parts::convdisp::sparsityPatternAxial(pattern.row(0), nComp, nCol, strideCell, u, weno); + FlowType::sparsityPattern(pattern.row(0), nComp, nCol, strideCell, u, weno); // Obtain memory for state const int nDof = nComp + nComp * nCol; @@ -486,7 +611,7 @@ void testBulkJacobianSparseBandedWeno(int wenoOrder, bool forwardFlow) // Populate sparse matrix cadet::linalg::CompressedSparseMatrix sparseMat(pattern); - cadet::model::parts::convdisp::residualKernelAxial(cadet::SimulationTime{0.0, 0u}, y.data(), nullptr, res.data(), sparseMat.row(0), fp); + FlowType::template residualWithJacobian(y.data(), nullptr, res.data(), sparseMat.row(0), fp); // Populate dense matrix cadet::linalg::BandMatrix bandMat; @@ -494,7 +619,8 @@ void testBulkJacobianSparseBandedWeno(int wenoOrder, bool forwardFlow) bandMat.resize(nComp * nCol, lowerBandwidth * strideCell, upperBandwidth * strideCell); else bandMat.resize(nComp * nCol, upperBandwidth * strideCell, lowerBandwidth * strideCell); - cadet::model::parts::convdisp::residualKernelAxial(cadet::SimulationTime{0.0, 0u}, y.data(), nullptr, res.data(), bandMat.row(0), fp); + + FlowType::template residualWithJacobian(y.data(), nullptr, res.data(), bandMat.row(0), fp); for (int col = 0; col < bandMat.rows(); ++col) { @@ -517,19 +643,19 @@ TEST_CASE("AxialConvectionDispersionOperator residual forward vs backward flow", { // Test all WENO orders for (unsigned int i = 1; i <= cadet::Weno::maxOrder(); ++i) - testResidualBulkWenoForwardBackward(i); + testResidualBulkWenoForwardBackward(i); } TEST_CASE("AxialConvectionDispersionOperator time derivative Jacobian vs FD", "[Operator],[AxialFlow],[Residual],[Jacobian]") { - testTimeDerivativeBulkJacobianFD(1e-6, 0.0, 1e-5); + testTimeDerivativeBulkJacobianFD(1e-6, 0.0, 1e-5); } TEST_CASE("AxialConvectionDispersionOperator Jacobian forward vs backward flow", "[Operator],[AxialFlow],[Residual],[Jacobian],[AD]") { // Test all WENO orders for (unsigned int i = 1; i <= cadet::Weno::maxOrder(); ++i) - testBulkJacobianWenoForwardBackward(i); + testBulkJacobianWenoForwardBackward(i); } TEST_CASE("AxialConvectionDispersionKernel Jacobian sparsity pattern vs FD", "[Operator],[AxialFlow],[Residual],[Jacobian],[SparseMatrix]") @@ -538,13 +664,13 @@ TEST_CASE("AxialConvectionDispersionKernel Jacobian sparsity pattern vs FD", "[O { // Test all WENO orders for (unsigned int i = 1; i <= cadet::Weno::maxOrder(); ++i) - testBulkJacobianSparsityWeno(i, true); + testBulkJacobianSparsityWeno(i, true); } SECTION("Backward flow") { // Test all WENO orders for (unsigned int i = 1; i <= cadet::Weno::maxOrder(); ++i) - testBulkJacobianSparsityWeno(i, false); + testBulkJacobianSparsityWeno(i, false); } } @@ -554,12 +680,64 @@ TEST_CASE("AxialConvectionDispersionKernel Jacobian sparse vs banded", "[Operato { // Test all WENO orders for (unsigned int i = 1; i <= cadet::Weno::maxOrder(); ++i) - testBulkJacobianSparseBandedWeno(i, true); + testBulkJacobianSparseBandedWeno(i, true); } SECTION("Backward flow") { // Test all WENO orders for (unsigned int i = 1; i <= cadet::Weno::maxOrder(); ++i) - testBulkJacobianSparseBandedWeno(i, false); + testBulkJacobianSparseBandedWeno(i, false); + } +} + + +TEST_CASE("RadialConvectionDispersionOperator residual forward vs backward flow", "[Operator],[RadialFlow],[Residual]") +{ + // Test all WENO orders + for (unsigned int i = 1; i <= 1; ++i) + testResidualBulkWenoForwardBackward(i); +} + +TEST_CASE("RadialConvectionDispersionOperator time derivative Jacobian vs FD", "[Operator],[RadialFlow],[Residual],[Jacobian]") +{ + testTimeDerivativeBulkJacobianFD(1e-6, 0.0, 1e-5); +} + +TEST_CASE("RadialConvectionDispersionOperator Jacobian forward vs backward flow", "[Operator],[RadialFlow],[Residual],[Jacobian],[AD]") +{ + // Test all WENO orders + for (unsigned int i = 1; i <= 1; ++i) + testBulkJacobianWenoForwardBackward(i); +} + +TEST_CASE("RadialConvectionDispersionKernel Jacobian sparsity pattern vs FD", "[Operator],[RadialFlow],[Residual],[Jacobian],[SparseMatrix]") +{ + SECTION("Forward flow") + { + // Test all WENO orders + for (unsigned int i = 1; i <= 1; ++i) + testBulkJacobianSparsityWeno(i, true); + } + SECTION("Backward flow") + { + // Test all WENO orders + for (unsigned int i = 1; i <= 1; ++i) + testBulkJacobianSparsityWeno(i, false); + } +} + +TEST_CASE("RadialConvectionDispersionKernel Jacobian sparse vs banded", "[Operator],[RadialFlow],[Jacobian],[SparseMatrix]") +{ + SECTION("Forward flow") + { + // Test all WENO orders + for (unsigned int i = 1; i <= 1; ++i) + testBulkJacobianSparseBandedWeno(i, true); + } + SECTION("Backward flow") + { + // Test all WENO orders + for (unsigned int i = 1; i <= 1; ++i) + testBulkJacobianSparseBandedWeno(i, false); } } diff --git a/test/Dummies.hpp b/test/Dummies.hpp new file mode 100644 index 000000000..b839779b9 --- /dev/null +++ b/test/Dummies.hpp @@ -0,0 +1,65 @@ +// ============================================================================= +// 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 "cadet/Model.hpp" +#include "cadet/ParameterId.hpp" + +#include "ConfigurationHelper.hpp" + +namespace +{ + class DummyConfigHelper : public cadet::IConfigHelper + { + public: + + DummyConfigHelper() { } + + virtual cadet::IInletProfile* createInletProfile(const std::string& type) const { return nullptr; } + virtual cadet::model::IBindingModel* createBindingModel(const std::string& name) const { return nullptr; } + virtual bool isValidBindingModel(const std::string& name) const { return false; } + virtual cadet::IExternalFunction* createExternalFunction(const std::string& type) const { return nullptr; } + virtual cadet::model::IDynamicReactionModel* createDynamicReactionModel(const std::string& name) const { return nullptr; } + virtual bool isValidDynamicReactionModel(const std::string& name) const { return false; } + virtual cadet::model::IParameterStateDependence* createParameterStateDependence(const std::string& name) const { return nullptr; } + virtual bool isValidParameterStateDependence(const std::string& name) const { return false; } + virtual cadet::model::IParameterParameterDependence* createParameterParameterDependence(const std::string& name) const { return nullptr; } + virtual bool isValidParameterParameterDependence(const std::string& name) const { return false; } + }; + + class DummyModel : public cadet::IModel + { + public: + DummyModel() { } + virtual ~DummyModel() CADET_NOEXCEPT { } + + virtual cadet::UnitOpIdx unitOperationId() const CADET_NOEXCEPT { return 0; }; + virtual const char* unitOperationName() const CADET_NOEXCEPT { return "DUMMY"; } + + virtual bool setParameter(const cadet::ParameterId& pId, int value) { return false; } + virtual bool setParameter(const cadet::ParameterId& pId, double value) { return false; } + virtual bool setParameter(const cadet::ParameterId& pId, bool value) { return false; } + + virtual bool hasParameter(const cadet::ParameterId& pId) const { return false; } + + virtual std::unordered_map getAllParameterValues() const { return std::unordered_map(); } + + virtual double getParameterDouble(const cadet::ParameterId& pId) const { return 0.0; } + + virtual void useAnalyticJacobian(const bool analyticJac) { } + +#ifdef CADET_BENCHMARK_MODE + virtual std::vector benchmarkTimings() const { return std::vector(0); } + virtual char const* const* benchmarkDescriptions() const { return nullptr; } +#endif + }; + +} diff --git a/test/JsonTestModels.cpp b/test/JsonTestModels.cpp index ccb570239..45cc69de4 100644 --- a/test/JsonTestModels.cpp +++ b/test/JsonTestModels.cpp @@ -30,7 +30,6 @@ json createColumnWithSMAJson(const std::string& uoType) config["PAR_SURFDIFFUSION"] = {0.0, 0.0, 0.0, 0.0}; // Geometry - config["COL_RADIUS"] = 0.01; if (uoType.substr(0, 6) == "RADIAL") { config["CROSS_SECTION_AREA"] = 0.0003141592653589793; @@ -38,7 +37,10 @@ json createColumnWithSMAJson(const std::string& uoType) config["COL_RADIUS_OUTER"] = 0.004; } else + { config["COL_LENGTH"] = 0.014; + config["COL_RADIUS"] = 0.01; + } config["PAR_RADIUS"] = 4.5e-5; config["COL_POROSITY"] = 0.37; config["PAR_POROSITY"] = 0.75; @@ -321,12 +323,12 @@ json createLWEJson(const std::string& uoType) { // Connection list is 1x7 since we have 1 connection between // the two unit operations (and we need to have 7 columns) - sw["CONNECTIONS"] = {1.0, 0.0, -1.0, -1.0, -1.0, -1.0, 1.0}; + sw["CONNECTIONS"] = {1.0, 0.0, -1.0, -1.0, -1.0, -1.0, 6.683738370512285e-8}; // Connections: From unit operation 1 port -1 (i.e., all ports) // to unit operation 0 port -1 (i.e., all ports), // connect component -1 (i.e., all components) // to component -1 (i.e., all components) with - // volumetric flow rate 1.0 m^3/s + // volumetric flow rate 6.683738370512285e-8 m^3/s } con["switch_000"] = sw; diff --git a/test/ModelSystem.cpp b/test/ModelSystem.cpp index 342a7c516..a14e4ba54 100644 --- a/test/ModelSystem.cpp +++ b/test/ModelSystem.cpp @@ -24,6 +24,7 @@ #include "JacobianHelper.hpp" #include "ColumnTests.hpp" #include "Utils.hpp" +#include "Dummies.hpp" #include "model/UnitOperation.hpp" #include @@ -32,23 +33,6 @@ namespace { - class DummyConfigHelper : public cadet::IConfigHelper - { - public: - - DummyConfigHelper() { } - - virtual cadet::IInletProfile* createInletProfile(const std::string& type) const { return nullptr; } - virtual cadet::model::IBindingModel* createBindingModel(const std::string& name) const { return nullptr; } - virtual bool isValidBindingModel(const std::string& name) const { return false; } - virtual cadet::IExternalFunction* createExternalFunction(const std::string& type) const { return nullptr; } - virtual cadet::model::IDynamicReactionModel* createDynamicReactionModel(const std::string& name) const { return nullptr; } - virtual bool isValidDynamicReactionModel(const std::string& name) const { return false; } - virtual cadet::model::IParameterStateDependence* createParameterStateDependence(const std::string& name) const { return nullptr; } - virtual bool isValidParameterStateDependence(const std::string& name) const { return false; } - virtual cadet::model::IParameterParameterDependence* createParameterParameterDependence(const std::string& name) const { return nullptr; } - virtual bool isValidParameterParameterDependence(const std::string& name) const { return false; } - }; class DummyUnitOperation : public cadet::IUnitOperation { diff --git a/test/TwoDimConvectionDispersionOperator.cpp b/test/TwoDimConvectionDispersionOperator.cpp index da7d1bf3f..8cd5df32e 100644 --- a/test/TwoDimConvectionDispersionOperator.cpp +++ b/test/TwoDimConvectionDispersionOperator.cpp @@ -15,9 +15,11 @@ #include "model/parts/TwoDimensionalConvectionDispersionOperator.hpp" #include "Weno.hpp" +#include "ModelBuilderImpl.hpp" #include "ColumnTests.hpp" #include "Utils.hpp" +#include "Dummies.hpp" #include "common/JsonParameterProvider.hpp" @@ -30,6 +32,9 @@ namespace cadet::JsonParameterProvider jpp(R"json({ "COL_LENGTH": 10, "COL_RADIUS": 1, + "COL_RADIUS_INNER": 0.001, + "COL_RADIUS_OUTER": 0.004, + "CROSS_SECTION_AREA": 0.0003141592653589793, "COL_POROSITY": 0.37, "COL_DISPERSION": 1e-6, "COL_DISPERSION_RADIAL": [1e-4, 1e-4, 1e-4, 1e-4, 1e-4], @@ -61,7 +66,8 @@ namespace // Configure the operator typedef std::unordered_map ParameterMap; ParameterMap parameters; - REQUIRE(convDispOp.configureModelDiscretization(jpp, nComp, nCol, nRad, false)); + cadet::ModelBuilder builder; + REQUIRE(convDispOp.configureModelDiscretization(jpp, builder, nComp, nCol, nRad, false)); REQUIRE(convDispOp.configure(0, jpp, parameters)); } @@ -73,10 +79,10 @@ namespace // Central finite differences y[nInletDof + col] = ref * (1.0 + h); - convDispOp.residual(0.0, 0u, y, nullptr, jacCol1, false, cadet::WithoutParamSensitivity()); + convDispOp.residual(DummyModel(), 0.0, 0u, y, nullptr, jacCol1, false, cadet::WithoutParamSensitivity()); y[nInletDof + col] = ref * (1.0 - h); - convDispOp.residual(0.0, 0u, y, nullptr, jacCol2, false, cadet::WithoutParamSensitivity()); + convDispOp.residual(DummyModel(), 0.0, 0u, y, nullptr, jacCol2, false, cadet::WithoutParamSensitivity()); y[nInletDof + col] = ref; @@ -128,14 +134,14 @@ void testBulk2DJacobianWenoForwardBackward(int wenoOrder) convDispOp.setFlowRates(i, 1e-2 * convDispOp.crossSection(i) * convDispOp.columnPorosity(i), 0.0); convDispOp.notifyDiscontinuousSectionTransition(0.0, 0u); - convDispOp.residual(0.0, 0u, y.data(), nullptr, jacCol1.data(), true, cadet::WithoutParamSensitivity()); + convDispOp.residual(DummyModel(), 0.0, 0u, y.data(), nullptr, jacCol1.data(), true, cadet::WithoutParamSensitivity()); // Compare Jacobian pattern against FD compareSparseJacobianAgainstFD(convDispOp, nInletDof, nPureDof, y.data(), jacCol1.data(), jacCol2.data(), h, relTol, absTol); // Reverse flow convDispOp.notifyDiscontinuousSectionTransition(0.0, 1u); - convDispOp.residual(0.0, 1u, y.data(), nullptr, jacCol1.data(), true, cadet::WithoutParamSensitivity()); + convDispOp.residual(DummyModel(), 0.0, 1u, y.data(), nullptr, jacCol1.data(), true, cadet::WithoutParamSensitivity()); // Compare Jacobian pattern against FD compareSparseJacobianAgainstFD(convDispOp, nInletDof, nPureDof, y.data(), jacCol1.data(), jacCol2.data(), h, relTol, absTol); @@ -170,7 +176,7 @@ void testBulk2DJacobianSparsityWeno(int wenoOrder, bool forwardFlow) convDispOp.setFlowRates(i, 1e-2 * convDispOp.crossSection(i) * convDispOp.columnPorosity(i), 0.0); convDispOp.notifyDiscontinuousSectionTransition(0.0, 0u); - convDispOp.residual(0.0, 0u, y.data(), nullptr, jacCol1.data(), true, cadet::WithoutParamSensitivity()); + convDispOp.residual(DummyModel(), 0.0, 0u, y.data(), nullptr, jacCol1.data(), true, cadet::WithoutParamSensitivity()); // Compare Jacobian pattern with FD for (int col = 0; col < nPureDof; ++col) @@ -179,10 +185,10 @@ void testBulk2DJacobianSparsityWeno(int wenoOrder, bool forwardFlow) // Central finite differences y[nInletDof + col] = ref * (1.0 + h); - convDispOp.residual(0.0, 0u, y.data(), nullptr, jacCol1.data(), false, cadet::WithoutParamSensitivity()); + convDispOp.residual(DummyModel(), 0.0, 0u, y.data(), nullptr, jacCol1.data(), false, cadet::WithoutParamSensitivity()); y[nInletDof + col] = ref * (1.0 - h); - convDispOp.residual(0.0, 0u, y.data(), nullptr, jacCol2.data(), false, cadet::WithoutParamSensitivity()); + convDispOp.residual(DummyModel(), 0.0, 0u, y.data(), nullptr, jacCol2.data(), false, cadet::WithoutParamSensitivity()); y[nInletDof + col] = ref; diff --git a/test/testRadialKernel.cpp b/test/testRadialKernel.cpp index 9985b4527..7a4287889 100644 --- a/test/testRadialKernel.cpp +++ b/test/testRadialKernel.cpp @@ -79,7 +79,7 @@ class RadialFlowModel : public cadet::test::IDiffEqModel _params.nComp = _nComp; _params.offsetToInlet = 0; _params.strideCell = _nComp; - _params.parDep = new cadet::model::DummyParameterParameterDependence(); + _params.parDep = new cadet::model::ConstantOneParameterParameterDependence(); } virtual ~RadialFlowModel() CADET_NOEXCEPT