From f8f20c018295926f884c8fd8b5613f3bb05e7107 Mon Sep 17 00:00:00 2001 From: Samuel Leweke Date: Fri, 17 Mar 2023 14:31:48 +0100 Subject: [PATCH] WIP add parameter dependencies to chromatography models Support radial flow and param deps in createLWE Add support for radial flow models and parameter dependence (i.e., dispersion and film diffusion coefficients) to createLWE Add axial dispersion and film diffusion parameter dependence Axial dispersion for all models, film diffusion parameter dependence for LRMP and GRM Add tests for radial flow convection dispersion operators Enable radial dispersion coeff dependency in operators Co-authored-by: jbreue16 --- src/libcadet/model/GeneralRateModel.cpp | 2 +- src/libcadet/model/GeneralRateModel2D.cpp | 4 +- .../model/LumpedRateModelWithPores.cpp | 81 +++++- .../model/LumpedRateModelWithPores.hpp | 1 + .../model/LumpedRateModelWithoutPores.cpp | 18 +- src/libcadet/model/ParameterDependence.hpp | 21 +- .../paramdep/DummyParameterDependence.cpp | 133 +++++++-- .../paramdep/ParameterDependenceBase.hpp | 38 +-- .../paramdep/PowerLawParameterDependence.cpp | 25 +- .../parts/AxialConvectionDispersionKernel.hpp | 36 ++- .../parts/ConvectionDispersionOperator.cpp | 202 ++++++++++--- .../parts/ConvectionDispersionOperator.hpp | 61 ++-- .../RadialConvectionDispersionKernel.hpp | 30 +- ...imensionalConvectionDispersionOperator.cpp | 41 +-- ...imensionalConvectionDispersionOperator.hpp | 18 +- src/tools/ToolsHelper.hpp | 8 + src/tools/createLWE.cpp | 35 ++- test/ConvectionDispersionOperator.cpp | 266 +++++++++++++++--- test/Dummies.hpp | 65 +++++ test/JsonTestModels.cpp | 13 +- test/ModelSystem.cpp | 18 +- test/TwoDimConvectionDispersionOperator.cpp | 22 +- test/testRadialKernel.cpp | 2 +- 23 files changed, 874 insertions(+), 266 deletions(-) create mode 100644 test/Dummies.hpp 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 323928119..868069a4c 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 has to be set"); @@ -683,8 +760,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 @@ -693,44 +777,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); @@ -746,7 +830,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); @@ -832,6 +918,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")) { @@ -878,6 +974,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")) { @@ -930,6 +1037,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")) { @@ -1082,8 +1200,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 @@ -1140,6 +1257,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 @@ -1149,33 +1267,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 46d43da5a..4010d4ac0 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 6127785fe..3ee40e366 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-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; } @@ -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..8bf5e35a7 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,28 @@ 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("PAR_RADIUS", 4.5e-5); + writer.scalar("PAR_CORERADIUS", 0.0); writer.scalar("COL_POROSITY", 0.37); const double par_radius = 4.5e-5; const double par_coreradius = 0.0; @@ -447,12 +470,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 fd63b36e0..83d92ca9d 100644 --- a/test/JsonTestModels.cpp +++ b/test/JsonTestModels.cpp @@ -31,6 +31,9 @@ json createColumnWithSMAJson(const std::string& uoType) // Geometry config["COL_LENGTH"] = 0.014; config["COL_RADIUS"] = 0.01; + config["COL_RADIUS_INNER"] = 0.001; + config["COL_RADIUS_OUTER"] = 0.004; + config["CROSS_SECTION_AREA"] = 0.0003141592653589793; config["PAR_RADIUS"] = 4.5e-5; config["COL_POROSITY"] = 0.37; config["PAR_POROSITY"] = 0.75; @@ -159,6 +162,9 @@ json createColumnWithTwoCompLinearJson(const std::string& uoType) // Geometry config["COL_LENGTH"] = 0.014; config["COL_RADIUS"] = 0.01; + config["COL_RADIUS_INNER"] = 0.001; + config["COL_RADIUS_OUTER"] = 0.004; + config["CROSS_SECTION_AREA"] = 0.0003141592653589793; config["PAR_RADIUS"] = 4.5e-5; config["COL_POROSITY"] = 0.37; config["PAR_POROSITY"] = 0.75; @@ -304,12 +310,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; @@ -444,6 +450,9 @@ cadet::JsonParameterProvider createPulseInjectionColumn(const std::string& uoTyp // Geometry grm["COL_LENGTH"] = 0.014; grm["COL_RADIUS"] = 0.01; + grm["COL_RADIUS_INNER"] = 0.001; + grm["COL_RADIUS_OUTER"] = 0.004; + grm["CROSS_SECTION_AREA"] = 0.0003141592653589793; grm["PAR_RADIUS"] = 4.5e-5; grm["PAR_CORERADIUS"] = 0.0; grm["COL_POROSITY"] = 0.37; 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