diff --git a/src/libcadet/CMakeLists.txt b/src/libcadet/CMakeLists.txt index aa8312a41..a0f5e212f 100644 --- a/src/libcadet/CMakeLists.txt +++ b/src/libcadet/CMakeLists.txt @@ -146,7 +146,6 @@ set(LIBCADET_SOURCES ${CMAKE_SOURCE_DIR}/src/libcadet/model/paramdep/ParameterDependenceBase.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/model/paramdep/LiquidSaltSolidParameterDependence.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/model/paramdep/DummyParameterDependence.cpp - ${CMAKE_SOURCE_DIR}/src/libcadet/model/paramdep/IdentityParameterDependence.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/model/paramdep/PowerLawParameterDependence.cpp ) diff --git a/src/libcadet/ParameterDependenceFactory.cpp b/src/libcadet/ParameterDependenceFactory.cpp index ad165dfe9..359cdd14d 100644 --- a/src/libcadet/ParameterDependenceFactory.cpp +++ b/src/libcadet/ParameterDependenceFactory.cpp @@ -22,8 +22,6 @@ namespace cadet void registerLiquidSaltSolidParamDependence(std::unordered_map>& paramDeps); void registerDummyParamDependence(std::unordered_map>& paramDeps); void registerDummyParamDependence(std::unordered_map>& paramDeps); - void registerIdentityParamDependence(std::unordered_map>& paramDeps); - void registerIdentityParamDependence(std::unordered_map>& paramDeps); void registerPowerLawParamDependence(std::unordered_map>& paramDeps); } } @@ -33,11 +31,9 @@ namespace cadet // Register all ParamState dependencies model::paramdep::registerLiquidSaltSolidParamDependence(_paramStateDeps); model::paramdep::registerDummyParamDependence(_paramStateDeps); - model::paramdep::registerIdentityParamDependence(_paramStateDeps); // Register all ParamParam dependencies model::paramdep::registerDummyParamDependence(_paramParamDeps); - model::paramdep::registerIdentityParamDependence(_paramParamDeps); model::paramdep::registerPowerLawParamDependence(_paramParamDeps); } diff --git a/src/libcadet/model/GeneralRateModel.cpp b/src/libcadet/model/GeneralRateModel.cpp index 5154f6b37..7c8dae0da 100644 --- a/src/libcadet/model/GeneralRateModel.cpp +++ b/src/libcadet/model/GeneralRateModel.cpp @@ -23,6 +23,7 @@ #include "model/ReactionModel.hpp" #include "model/ParameterDependence.hpp" #include "model/parts/BindingCellKernel.hpp" +#include "model/ParameterDependence.hpp" #include "SimulationTypes.hpp" #include "linalg/DenseMatrix.hpp" #include "linalg/BandMatrix.hpp" @@ -67,7 +68,7 @@ int schurComplementMultiplierGRM(void* userData, double const* x, double* z) template GeneralRateModel::GeneralRateModel(UnitOpIdx unitOpIdx) : UnitOperationBase(unitOpIdx), - _hasSurfaceDiffusion(0, false), _dynReactionBulk(nullptr), + _hasSurfaceDiffusion(0, false), _dynReactionBulk(nullptr), _filmDiffDep(nullptr), _jacP(nullptr), _jacPdisc(nullptr), _jacPF(nullptr), _jacFP(nullptr), _jacInlet(), _hasParDepSurfDiffusion(false), _analyticJac(true), _jacobianAdDirs(0), _factorizeJacobian(false), _tempState(nullptr), _initC(0), _initCp(0), _initQ(0), _initState(0), _initStateDot(0) @@ -86,6 +87,7 @@ GeneralRateModel::~GeneralRateModel() CADET_NOEXCEPT delete[] _jacPdisc; delete _dynReactionBulk; + delete _filmDiffDep; clearParDepSurfDiffusion(); @@ -447,6 +449,18 @@ bool GeneralRateModel::configureModelDiscretization(IParameter 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("IDENTITY"); + // ==== Construct and configure binding model clearBindingModels(); _binding = std::vector(_disc.nParType, nullptr); @@ -651,6 +665,12 @@ bool GeneralRateModel::configure(IParameterProvider& paramProv 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); @@ -1897,53 +1917,44 @@ int GeneralRateModel::residualFlux(double t, unsigned int secI const ParamType jacCF_val = invBetaC * surfaceToVolumeRatio; const ParamType jacPF_val = -outerAreaPerVolume / epsP; - // Discretized film diffusion kf for finite volumes - if (cadet_likely(_colParBoundaryOrder == 2)) + for (unsigned int colCell = 0; colCell < _disc.nCol * _disc.nComp; ++colCell) { + // Discretized film diffusion kf for finite volumes const ParamType absOuterShellHalfRadius = 0.5 * static_cast(_parCellSize[_disc.nParCellsBeforeType[type]]); - for (unsigned int comp = 0; comp < _disc.nComp; ++comp) - kf_FV[comp] = 1.0 / (absOuterShellHalfRadius / epsP / static_cast(_poreAccessFactor[type * _disc.nComp + comp]) / static_cast(parDiff[comp]) + 1.0 / static_cast(filmDiff[comp])); - } - else - { - for (unsigned int comp = 0; comp < _disc.nComp; ++comp) - kf_FV[comp] = static_cast(filmDiff[comp]); - } - // J_{0,f} block, adds flux to column void / bulk volume equations - for (unsigned int i = 0; i < _disc.nCol * _disc.nComp; ++i) - { - const unsigned int colCell = i / _disc.nComp; - resCol[i] += jacCF_val * static_cast(_parTypeVolFrac[type + colCell * _disc.nParType]) * yFluxType[i]; - } + for (unsigned int comp = 0; comp < _disc.nComp; ++comp) { - // J_{f,0} block, adds bulk volume state c_i to flux equation - for (unsigned int bnd = 0; bnd < _disc.nCol; ++bnd) - { + const double relPos = _convDispOp.relativeCoordinate(colCell); + const active curVelocity = _convDispOp.currentVelocity(relPos); + const ParamType filmDiffDep = static_cast(_filmDiffDep->getValue(*this, ColumnPosition{ relPos, 0.0, 0.0 }, comp, ParTypeIndep, BoundStateIndep, filmDiff[comp], curVelocity)); + if (cadet_likely(_colParBoundaryOrder == 2)) + kf_FV[comp] = 1.0 / (absOuterShellHalfRadius / epsP / static_cast(_poreAccessFactor[type * _disc.nComp + comp]) / static_cast(parDiff[comp]) + 1.0 / filmDiffDep); + else + kf_FV[comp] = filmDiffDep; + } + + // J_{0,f} block, adds flux to column void / bulk volume equations + resCol[colCell] += jacCF_val * static_cast(_parTypeVolFrac[type + colCell * _disc.nParType]) * yFluxType[colCell]; + + // J_{f,0} block, adds bulk volume state c_i to flux equation for (unsigned int comp = 0; comp < _disc.nComp; ++comp) { - const unsigned int eq = bnd * idxr.strideColCell() + comp * idxr.strideColComp(); + const unsigned int eq = colCell * idxr.strideColCell() + comp * idxr.strideColComp(); resFluxType[eq] -= kf_FV[comp] * yCol[eq]; } - } - // J_{p,f} block, implements bead boundary condition in outer bead shell equation - for (unsigned int pblk = 0; pblk < _disc.nCol; ++pblk) - { + // J_{p,f} block, implements bead boundary condition in outer bead shell equation for (unsigned int comp = 0; comp < _disc.nComp; ++comp) { - const unsigned int eq = pblk * idxr.strideColCell() + comp * idxr.strideColComp(); - resParType[pblk * idxr.strideParBlock(type) + comp] += jacPF_val / static_cast(_poreAccessFactor[type * _disc.nComp + comp]) * yFluxType[eq]; + const unsigned int eq = colCell * idxr.strideColCell() + comp * idxr.strideColComp(); + resParType[colCell * idxr.strideParBlock(type) + comp] += jacPF_val / static_cast(_poreAccessFactor[type * _disc.nComp + comp]) * yFluxType[eq]; } - } - // J_{f,p} block, adds outer bead shell state c_{p,i} to flux equation - for (unsigned int pblk = 0; pblk < _disc.nCol; ++pblk) - { + // J_{f,p} block, adds outer bead shell state c_{p,i} to flux equation for (unsigned int comp = 0; comp < _disc.nComp; ++comp) { - const unsigned int eq = pblk * idxr.strideColCell() + comp * idxr.strideColComp(); - resFluxType[eq] += kf_FV[comp] * yParType[comp + pblk * idxr.strideParBlock(type)]; + const unsigned int eq = colCell * idxr.strideColCell() + comp * idxr.strideColComp(); + resFluxType[eq] += kf_FV[comp] * yParType[comp + colCell * idxr.strideParBlock(type)]; } } @@ -1957,11 +1968,17 @@ int GeneralRateModel::residualFlux(double t, unsigned int secI active const* const parCenterRadius = _parCenterRadius.data() + _disc.nParCellsBeforeType[type]; const ParamType absOuterShellHalfRadius = 0.5 * static_cast(_parCellSize[_disc.nParCellsBeforeType[type]]); - for (unsigned int comp = 0; comp < _disc.nComp; ++comp) - kf_FV[comp] = (1.0 - static_cast(_parPorosity[type])) / (1.0 + epsP * static_cast(_poreAccessFactor[type * _disc.nComp + comp]) * static_cast(parDiff[comp]) / (absOuterShellHalfRadius * static_cast(filmDiff[comp]))); - for (unsigned int pblk = 0; pblk < _disc.nCol; ++pblk) { + for (unsigned int comp = 0; comp < _disc.nComp; ++comp) { + + const double relPos = _convDispOp.relativeCoordinate(pblk); + const active curVelocity = _convDispOp.currentVelocity(relPos); + const ParamType filmDiffDep = static_cast(_filmDiffDep->getValue(*this, ColumnPosition{ relPos, 0.0, 0.0 }, comp, ParTypeIndep, BoundStateIndep, filmDiff[comp], curVelocity)); + + kf_FV[comp] = (1.0 - static_cast(_parPorosity[type])) / (1.0 + epsP * static_cast(_poreAccessFactor[type * _disc.nComp + comp]) * static_cast(parDiff[comp]) / (absOuterShellHalfRadius * filmDiffDep)); + } + const ColumnPosition colPos{(0.5 + static_cast(pblk)) / static_cast(_disc.nCol), 0.0, static_cast(parCenterRadius[0]) / static_cast(_parRadius[type])}; const ParamType dr = static_cast(parCenterRadius[0]) - static_cast(parCenterRadius[1]); @@ -2063,54 +2080,49 @@ void GeneralRateModel::assembleOffdiagJac(double t, unsigned i const double jacCF_val = invBetaC * surfaceToVolumeRatio; const double jacPF_val = -outerAreaPerVolume / epsP; - // Discretized film diffusion kf for finite volumes - if (cadet_likely(_colParBoundaryOrder == 2)) - { + for (unsigned int pblk = 0; pblk < _disc.nCol; ++pblk) { + + // Discretized film diffusion kf for finite volumes const double absOuterShellHalfRadius = 0.5 * static_cast(_parCellSize[_disc.nParCellsBeforeType[type]]); - for (unsigned int comp = 0; comp < _disc.nComp; ++comp) - kf_FV[comp] = 1.0 / (absOuterShellHalfRadius / epsP / static_cast(_poreAccessFactor[type * _disc.nComp + comp]) / static_cast(parDiff[comp]) + 1.0 / static_cast(filmDiff[comp])); - } - else - { - for (unsigned int comp = 0; comp < _disc.nComp; ++comp) - kf_FV[comp] = static_cast(filmDiff[comp]); - } - // J_{0,f} block, adds flux to column void / bulk volume equations - for (unsigned int eq = 0; eq < _disc.nCol * _disc.nComp; ++eq) - { - const unsigned int colCell = eq / _disc.nComp; + for (unsigned int comp = 0; comp < _disc.nComp; ++comp) { - // Main diagonal corresponds to j_{f,i} (flux) state variable - _jacCF.addElement(eq, eq + typeOffset, jacCF_val * static_cast(_parTypeVolFrac[type + colCell * _disc.nParType])); - } + const double relPos = _convDispOp.relativeCoordinate(pblk); + const double curVelocity = static_cast(_convDispOp.currentVelocity(relPos)); + const double filmDiffDep = static_cast(_filmDiffDep->getValue(*this, ColumnPosition{ relPos, 0.0, 0.0 }, comp, ParTypeIndep, BoundStateIndep, filmDiff[comp], curVelocity)); - // J_{f,0} block, adds bulk volume state c_i to flux equation - for (unsigned int col = 0; col < _disc.nCol; ++col) - { + if (cadet_likely(_colParBoundaryOrder == 2)) + kf_FV[comp] = 1.0 / (absOuterShellHalfRadius / epsP / static_cast(_poreAccessFactor[type * _disc.nComp + comp]) / static_cast(parDiff[comp]) + 1.0 / filmDiffDep); + else + kf_FV[comp] = filmDiffDep; + } + + // J_{0,f} block, adds flux to column void / bulk volume equations + for (unsigned int comp = 0; comp < _disc.nComp; comp++) + { + const unsigned int eq = pblk * _disc.nComp + comp; + // Main diagonal corresponds to j_{f,i} (flux) state variable + _jacCF.addElement(eq, eq + typeOffset, jacCF_val * static_cast(_parTypeVolFrac[type + pblk * _disc.nParType])); + } + + // J_{f,0} block, adds bulk volume state c_i to flux equation for (unsigned int comp = 0; comp < _disc.nComp; ++comp) { // Main diagonal corresponds to c_i state variable in each column cell - const unsigned int eq = col * idxr.strideColCell() + comp * idxr.strideColComp(); + const unsigned int eq = pblk * idxr.strideColCell() + comp * idxr.strideColComp(); _jacFC.addElement(eq + typeOffset, eq, -kf_FV[comp]); } - } - // J_{p,f} block, implements bead boundary condition in outer bead shell equation - linalg::DoubleSparseMatrix* const jacPFtype = _jacPF + type * _disc.nCol; - for (unsigned int pblk = 0; pblk < _disc.nCol; ++pblk) - { + // J_{p,f} block, implements bead boundary condition in outer bead shell equation + linalg::DoubleSparseMatrix* const jacPFtype = _jacPF + type * _disc.nCol; for (unsigned int comp = 0; comp < _disc.nComp; ++comp) { const unsigned int eq = typeOffset + pblk * idxr.strideColCell() + comp * idxr.strideColComp(); jacPFtype[pblk].addElement(comp, eq, jacPF_val / static_cast(_poreAccessFactor[comp])); } - } - // J_{f,p} block, adds outer bead shell state c_{p,i} to flux equation - linalg::DoubleSparseMatrix* const jacFPtype = _jacFP + type * _disc.nCol; - for (unsigned int pblk = 0; pblk < _disc.nCol; ++pblk) - { + // J_{f,p} block, adds outer bead shell state c_{p,i} to flux equation + linalg::DoubleSparseMatrix* const jacFPtype = _jacFP + type * _disc.nCol; for (unsigned int comp = 0; comp < _disc.nComp; ++comp) { const unsigned int eq = typeOffset + pblk * idxr.strideColCell() + comp * idxr.strideColComp(); @@ -2122,21 +2134,29 @@ void GeneralRateModel::assembleOffdiagJac(double t, unsigned i { int const* const qsReaction = _binding[type]->reactionQuasiStationarity(); + linalg::DoubleSparseMatrix* const jacFPtype = _jacFP + type * _disc.nCol; + // Ordering of particle surface diffusion: // bnd0comp0, bnd0comp1, bnd0comp2, bnd1comp0, bnd1comp1, bnd1comp2 active const* const parSurfDiff = getSectionDependentSlice(_parSurfDiffusion, _disc.strideBound[_disc.nParType], secIdx) + _disc.nBoundBeforeType[type]; active const* const parCenterRadius = _parCenterRadius.data() + _disc.nParCellsBeforeType[type]; const double absOuterShellHalfRadius = 0.5 * static_cast(_parCellSize[_disc.nParCellsBeforeType[type]]); - for (unsigned int comp = 0; comp < _disc.nComp; ++comp) - kf_FV[comp] = (1.0 - static_cast(_parPorosity[type])) / (1.0 + epsP * static_cast(_poreAccessFactor[type * _disc.nComp + comp]) * static_cast(parDiff[comp]) / (absOuterShellHalfRadius * static_cast(filmDiff[comp]))); - for (unsigned int pblk = 0; pblk < _disc.nCol; ++pblk) { - const ColumnPosition colPos{(0.5 + static_cast(pblk)) / static_cast(_disc.nCol), 0.0, static_cast(parCenterRadius[0]) / static_cast(_parRadius[type])}; + for (unsigned int comp = 0; comp < _disc.nComp; ++comp) + { + const double relPos = _convDispOp.relativeCoordinate(pblk); + const double curVelocity = static_cast(_convDispOp.currentVelocity(relPos)); + const double filmDiffDep = static_cast(_filmDiffDep->getValue(*this, ColumnPosition{ relPos, 0.0, 0.0 }, comp, ParTypeIndep, BoundStateIndep, filmDiff[comp], curVelocity)); + + kf_FV[comp] = (1.0 - static_cast(_parPorosity[type])) / (1.0 + epsP * static_cast(_poreAccessFactor[type * _disc.nComp + comp]) * static_cast(parDiff[comp]) / (absOuterShellHalfRadius * filmDiffDep)); + } + + const ColumnPosition colPos{ (0.5 + static_cast(pblk)) / static_cast(_disc.nCol), 0.0, static_cast(parCenterRadius[0]) / static_cast(_parRadius[type]) }; const double dr = static_cast(parCenterRadius[0]) - static_cast(parCenterRadius[1]); - double const* const yCell = vecStateY + idxr.offsetCp(ParticleTypeIndex{type}, ParticleIndex{pblk}); + double const* const yCell = vecStateY + idxr.offsetCp(ParticleTypeIndex{ type }, ParticleIndex{ pblk }); for (unsigned int comp = 0; comp < _disc.nComp; ++comp) { @@ -2145,7 +2165,7 @@ void GeneralRateModel::assembleOffdiagJac(double t, unsigned i for (unsigned int i = 0; i < nBound; ++i) { - const int idxBnd = idxr.offsetBoundComp(ParticleTypeIndex{type}, ComponentIndex{comp}) + i; + const int idxBnd = idxr.offsetBoundComp(ParticleTypeIndex{ type }, ComponentIndex{ comp }) + i; // Skip quasi-stationary bound states if (!qsReaction[idxBnd]) @@ -2838,6 +2858,15 @@ bool GeneralRateModel::setParameter(const ParameterId& pId, do if (_convDispOp.setParameter(pId, value)) return true; + + if (_filmDiffDep) + { + if (_filmDiffDep->hasParameter(pId)) + { + _filmDiffDep->setParameter(pId, value); + return true; + } + } } const bool result = UnitOperationBase::setParameter(pId, value); @@ -2918,6 +2947,16 @@ void GeneralRateModel::setSensitiveParameterValue(const Parame 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); @@ -3012,6 +3051,17 @@ bool GeneralRateModel::setSensitiveParameter(const ParameterId 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; + } + } } const bool result = UnitOperationBase::setSensitiveParameter(pId, adDirection, adValue); diff --git a/src/libcadet/model/GeneralRateModel.hpp b/src/libcadet/model/GeneralRateModel.hpp index f0dee376f..fd817316b 100644 --- a/src/libcadet/model/GeneralRateModel.hpp +++ b/src/libcadet/model/GeneralRateModel.hpp @@ -314,6 +314,7 @@ class GeneralRateModel : 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 linalg::BandMatrix* _jacP; //!< Particle jacobian diagonal blocks (all of them) linalg::FactorizableBandMatrix* _jacPdisc; //!< Particle jacobian diagonal blocks (all of them) with time derivatives from BDF method diff --git a/src/libcadet/model/LumpedRateModelWithPores.cpp b/src/libcadet/model/LumpedRateModelWithPores.cpp index a92857508..00ef01665 100644 --- a/src/libcadet/model/LumpedRateModelWithPores.cpp +++ b/src/libcadet/model/LumpedRateModelWithPores.cpp @@ -275,7 +275,7 @@ bool LumpedRateModelWithPores::configureModelDiscretization(IP _filmDiffDep->configureModelDiscretization(paramProvider); } else - _filmDiffDep = helper.createParameterParameterDependence("CONSTANT_ONE"); + _filmDiffDep = helper.createParameterParameterDependence("IDENTITY"); // Allocate memory Indexer idxr(_disc); @@ -1091,9 +1091,9 @@ int LumpedRateModelWithPores::residualFlux(double t, unsigned const double relPos = _convDispOp.relativeCoordinate(colCell); const active curVelocity = _convDispOp.currentVelocity(relPos); - const active modifier = _filmDiffDep->getValue(*this, ColumnPosition{relPos, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, curVelocity); + const ParamType filmDiffDep = static_cast(_filmDiffDep->getValue(*this, ColumnPosition{relPos, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, filmDiff[comp], curVelocity)); - resCol[i] += jacCF_val * static_cast(filmDiff[comp]) * static_cast(modifier) * static_cast(_parTypeVolFrac[type + _disc.nParType * colCell]) * yFluxType[i]; + resCol[i] += jacCF_val * filmDiffDep * static_cast(_parTypeVolFrac[type + _disc.nParType * colCell]) * yFluxType[i]; } // J_{f,0} block, adds bulk volume state c_i to flux equation @@ -1114,10 +1114,10 @@ int LumpedRateModelWithPores::residualFlux(double t, unsigned for (unsigned int comp = 0; comp < _disc.nComp; ++comp) { - const active modifier = _filmDiffDep->getValue(*this, ColumnPosition{relPos, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, curVelocity); + const ParamType filmDiffDep = static_cast(_filmDiffDep->getValue(*this, ColumnPosition{relPos, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, filmDiff[comp], curVelocity)); const unsigned int eq = pblk * idxr.strideColCell() + comp * idxr.strideColComp(); - resParType[pblk * idxr.strideParBlock(type) + comp] += jacPF_val / static_cast(poreAccFactor[comp]) * static_cast(filmDiff[comp]) * static_cast(modifier) * yFluxType[eq]; + resParType[pblk * idxr.strideParBlock(type) + comp] += jacPF_val / static_cast(poreAccFactor[comp]) * filmDiffDep * yFluxType[eq]; } } @@ -1180,10 +1180,10 @@ void LumpedRateModelWithPores::assembleOffdiagJac(double t, un 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); + const double filmDiffDep = static_cast(_filmDiffDep->getValue(*this, ColumnPosition{relPos, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, filmDiff[comp], curVelocity)); // Main diagonal corresponds to j_{f,i} (flux) state variable - _jacCF.addElement(eq, eq + typeOffset, jacCF_val * static_cast(filmDiff[comp]) * modifier * static_cast(_parTypeVolFrac[type + _disc.nParType * colCell])); + _jacCF.addElement(eq, eq + typeOffset, jacCF_val * filmDiffDep * static_cast(_parTypeVolFrac[type + _disc.nParType * colCell])); } // J_{f,0} block, adds bulk volume state c_i to flux equation @@ -1203,8 +1203,8 @@ void LumpedRateModelWithPores::assembleOffdiagJac(double t, un const unsigned int eq = typeOffset + pblk * idxr.strideColCell() + comp * idxr.strideColComp(); const unsigned int col = pblk * idxr.strideParBlock(type) + comp; - const double modifier = _filmDiffDep->getValue(*this, ColumnPosition{relPos, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, curVelocity); - _jacPF[type].addElement(col, eq, jacPF_val / static_cast(poreAccFactor[comp]) * static_cast(filmDiff[comp]) * modifier); + const double filmDiffDep = static_cast(_filmDiffDep->getValue(*this, ColumnPosition{relPos, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, filmDiff[comp], curVelocity)); + _jacPF[type].addElement(col, eq, jacPF_val / static_cast(poreAccFactor[comp]) * filmDiffDep); } } diff --git a/src/libcadet/model/ParameterDependence.hpp b/src/libcadet/model/ParameterDependence.hpp index a463ca02a..1e981a2cc 100644 --- a/src/libcadet/model/ParameterDependence.hpp +++ b/src/libcadet/model/ParameterDependence.hpp @@ -460,9 +460,11 @@ class IParameterParameterDependence * @param [in] comp Index of the component the parameter belongs to (or @c -1 if independent of components) * @param [in] parType Index of the particle type the parameter belongs to (or @c -1 if independent of particle types) * @param [in] bnd Index of the bound state the parameter belongs to (or @c -1 if independent of bound states) + * @param [in] depVal parameter-dependent value + * @param [in] indepVal parameter depVal depends on * @return Actual parameter value */ - virtual double getValue(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd) const = 0; + virtual double getValue(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd, double depVal, double indepVal) const = 0; /** * @brief Evaluates the parameter @@ -473,37 +475,11 @@ class IParameterParameterDependence * @param [in] comp Index of the component the parameter belongs to (or @c -1 if independent of components) * @param [in] parType Index of the particle type the parameter belongs to (or @c -1 if independent of particle types) * @param [in] bnd Index of the bound state the parameter belongs to (or @c -1 if independent of bound states) + * @param [in] depVal parameter-dependent value + * @param [in] indepVal parameter depVal depends on * @return Actual parameter value */ - virtual active getValueActive(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd) const = 0; - - /** - * @brief Evaluates the parameter - * @details This function is called simultaneously from multiple threads. - * - * @param [in] model Model that owns this parameter dependence - * @param [in] colPos Position in normalized coordinates (column inlet = 0, column outlet = 1; outer shell = 1, inner center = 0) - * @param [in] comp Index of the component the parameter belongs to (or @c -1 if independent of components) - * @param [in] parType Index of the particle type the parameter belongs to (or @c -1 if independent of particle types) - * @param [in] bnd Index of the bound state the parameter belongs to (or @c -1 if independent of bound states) - * @param [in] val Additional parameter-dependent value - * @return Actual parameter value - */ - virtual double getValue(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd, double val) const = 0; - - /** - * @brief Evaluates the parameter - * @details This function is called simultaneously from multiple threads. - * - * @param [in] model Model that owns this parameter dependence - * @param [in] colPos Position in normalized coordinates (column inlet = 0, column outlet = 1; outer shell = 1, inner center = 0) - * @param [in] comp Index of the component the parameter belongs to (or @c -1 if independent of components) - * @param [in] parType Index of the particle type the parameter belongs to (or @c -1 if independent of particle types) - * @param [in] bnd Index of the bound state the parameter belongs to (or @c -1 if independent of bound states) - * @param [in] val Additional parameter-dependent value - * @return Actual parameter value - */ - virtual active getValue(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd, const active& val) const = 0; + virtual active getValue(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd, const active& depVal, const active& indepVal) const = 0; protected: }; diff --git a/src/libcadet/model/paramdep/DummyParameterDependence.cpp b/src/libcadet/model/paramdep/DummyParameterDependence.cpp index c73c92753..b406155b9 100644 --- a/src/libcadet/model/paramdep/DummyParameterDependence.cpp +++ b/src/libcadet/model/paramdep/DummyParameterDependence.cpp @@ -26,6 +26,36 @@ namespace cadet namespace model { +/** + * @brief Defines a dummy parameter dependence that outputs the (independent) parameter itself + */ +class IdentityParameterParameterDependence : public ParameterParameterDependenceBase +{ +public: + + IdentityParameterParameterDependence() { } + virtual ~IdentityParameterParameterDependence() CADET_NOEXCEPT { } + + static const char* identifier() { return "IDENTITY"; } + virtual const char* name() const CADET_NOEXCEPT { return IdentityParameterParameterDependence::identifier(); } + + CADET_PARAMETERPARAMETERDEPENDENCE_BOILERPLATE + +protected: + + virtual bool configureImpl(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, BoundStateIdx bndIdx, const std::string& name) + { + return true; + } + + template + ParamType getValueImpl(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd, const ParamType& depVal, const ParamType& indepVal) const + { + return depVal; + } + +}; + /** * @brief Defines a parameter dependence that outputs constant 0.0 */ @@ -83,7 +113,6 @@ class ConstantZeroParameterStateDependence : public ParameterStateDependenceBase void analyticJacobianCombinedAddSolidImpl(const ColumnPosition& colPos, double param, double const* yLiquid, double const* ySolid, int bnd, double factor, int offset, RowIterator jac) const { } }; - /** * @brief Defines a parameter dependence that outputs constant 1.0 */ @@ -141,81 +170,6 @@ class ConstantOneParameterStateDependence : public ParameterStateDependenceBase void analyticJacobianCombinedAddSolidImpl(const ColumnPosition& colPos, double param, double const* yLiquid, double const* ySolid, int bnd, double factor, int offset, RowIterator jac) const { } }; - -/** - * @brief Defines a parameter dependence that outputs constant 0.0 - */ -class ConstantZeroParameterParameterDependence : public ParameterParameterDependenceBase -{ -public: - - ConstantZeroParameterParameterDependence() { } - virtual ~ConstantZeroParameterParameterDependence() CADET_NOEXCEPT { } - - static const char* identifier() { return "CONSTANT_ZERO"; } - virtual const char* name() const CADET_NOEXCEPT { return ConstantZeroParameterParameterDependence::identifier(); } - - CADET_PARAMETERPARAMETERDEPENDENCE_BOILERPLATE - -protected: - - virtual bool configureImpl(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, BoundStateIdx bndIdx, const std::string& name) - { - return true; - } - - template - ParamType getValueImpl(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd) const - { - return 0.0; - } - - template - ParamType getValueImpl(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd, const ParamType& val) const - { - return 0.0; - } - -}; - - -/** - * @brief Defines a parameter dependence that outputs constant 1.0 - */ -class ConstantOneParameterParameterDependence : public ParameterParameterDependenceBase -{ -public: - - ConstantOneParameterParameterDependence() { } - virtual ~ConstantOneParameterParameterDependence() CADET_NOEXCEPT { } - - static const char* identifier() { return "CONSTANT_ONE"; } - virtual const char* name() const CADET_NOEXCEPT { return ConstantOneParameterParameterDependence::identifier(); } - - CADET_PARAMETERPARAMETERDEPENDENCE_BOILERPLATE - -protected: - - virtual bool configureImpl(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, BoundStateIdx bndIdx, const std::string& name) - { - return true; - } - - template - 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) @@ -227,9 +181,8 @@ namespace paramdep void registerDummyParamDependence(std::unordered_map>& paramDeps) { - paramDeps[ConstantOneParameterParameterDependence::identifier()] = []() { return new ConstantOneParameterParameterDependence(); }; - paramDeps[ConstantZeroParameterParameterDependence::identifier()] = []() { return new ConstantZeroParameterParameterDependence(); }; - paramDeps["NONE"] = []() { return new ConstantOneParameterParameterDependence(); }; + paramDeps[IdentityParameterParameterDependence::identifier()] = []() { return new IdentityParameterParameterDependence(); }; + paramDeps["NONE"] = []() { return new IdentityParameterParameterDependence(); }; } } // namespace paramdep diff --git a/src/libcadet/model/paramdep/IdentityParameterDependence.cpp b/src/libcadet/model/paramdep/IdentityParameterDependence.cpp deleted file mode 100644 index 60092ed6b..000000000 --- a/src/libcadet/model/paramdep/IdentityParameterDependence.cpp +++ /dev/null @@ -1,141 +0,0 @@ -// ============================================================================= -// CADET -// -// Copyright ┬® 2008-2022: The CADET Authors -// Please see the AUTHORS and CONTRIBUTORS file. -// -// All rights reserved. This program and the accompanying materials -// are made available under the terms of the GNU Public License v3.0 (or, at -// your option, any later version) which accompanies this distribution, and -// is available at http://www.gnu.org/licenses/gpl.html -// ============================================================================= - -#include "model/paramdep/ParameterDependenceBase.hpp" -#include "ParamReaderHelper.hpp" -#include "cadet/Exceptions.hpp" -#include "SimulationTypes.hpp" -#include "cadet/ParameterId.hpp" - -#include -#include -#include - -namespace cadet -{ - -namespace model -{ - -/** - * @brief Defines an identity parameter dependence - */ -class IdentityParameterStateDependence : public ParameterStateDependenceBase -{ -public: - - IdentityParameterStateDependence() { } - virtual ~IdentityParameterStateDependence() CADET_NOEXCEPT { } - - static const char* identifier() { return "IDENTITY"; } - virtual const char* name() const CADET_NOEXCEPT { return IdentityParameterStateDependence::identifier(); } - - virtual int jacobianElementsPerRowLiquid() const CADET_NOEXCEPT { return 0; } - virtual int jacobianElementsPerRowCombined() const CADET_NOEXCEPT { return 0; } - - virtual void analyticJacobianLiquidAdd(const ColumnPosition& colPos, double param, double const* y, int comp, double factor, int offset, int row, linalg::DoubleSparseMatrix& jac) const { } - virtual void analyticJacobianCombinedAddLiquid(const ColumnPosition& colPos, double param, double const* yLiquid, double const* ySolid, int comp, double factor, int offset, int row, linalg::DoubleSparseMatrix& jac) const { } - virtual void analyticJacobianCombinedAddSolid(const ColumnPosition& colPos, double param, double const* yLiquid, double const* ySolid, int bnd, double factor, int offset, int row, linalg::DoubleSparseMatrix& jac) const { } - - CADET_PARAMETERSTATEDEPENDENCE_BOILERPLATE - -protected: - - virtual bool configureImpl(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, const std::string& name) - { - return true; - } - - template - typename DoubleActivePromoter::type liquidParameterImpl(const ColumnPosition& colPos, const ParamType& param, StateType const* y, int comp) const - { - return param; - } - - 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 param; - } - - 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 param; - } - - 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 an identity parameter dependence - */ -class IdentityParameterParameterDependence : public ParameterParameterDependenceBase -{ -public: - - IdentityParameterParameterDependence() { } - virtual ~IdentityParameterParameterDependence() CADET_NOEXCEPT { } - - static const char* identifier() { return "IDENTITY"; } - virtual const char* name() const CADET_NOEXCEPT { return IdentityParameterParameterDependence::identifier(); } - - CADET_PARAMETERPARAMETERDEPENDENCE_BOILERPLATE - -protected: - - virtual bool configureImpl(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, BoundStateIdx bndIdx, const std::string& name) - { - return true; - } - - template - ParamType getValueImpl(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd) const - { - return 0.0; - } - - template - ParamType getValueImpl(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd, const ParamType& val) const - { - return val; - } - -}; - - -namespace paramdep -{ - void registerIdentityParamDependence(std::unordered_map>& paramDeps) - { - paramDeps[IdentityParameterStateDependence::identifier()] = []() { return new IdentityParameterStateDependence(); }; - paramDeps["NONE"] = []() { return new IdentityParameterStateDependence(); }; - } - - void registerIdentityParamDependence(std::unordered_map>& paramDeps) - { - paramDeps[IdentityParameterParameterDependence::identifier()] = []() { return new IdentityParameterParameterDependence(); }; - paramDeps["NONE"] = []() { return new IdentityParameterParameterDependence(); }; - } -} // namespace paramdep - -} // namespace model - -} // namespace cadet diff --git a/src/libcadet/model/paramdep/ParameterDependenceBase.hpp b/src/libcadet/model/paramdep/ParameterDependenceBase.hpp index 1f72db332..80da1a2c3 100644 --- a/src/libcadet/model/paramdep/ParameterDependenceBase.hpp +++ b/src/libcadet/model/paramdep/ParameterDependenceBase.hpp @@ -241,27 +241,16 @@ class ParameterParameterDependenceBase : public IParameterParameterDependence * * The implementation is inserted inline in the class declaration. */ -#define CADET_PARAMETERPARAMETERDEPENDENCE_BOILERPLATE \ - virtual double getValue(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd, double val) const \ - { \ - return getValueImpl(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); \ - } - +#define CADET_PARAMETERPARAMETERDEPENDENCE_BOILERPLATE \ + virtual double getValue(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd, double depVal, double indepVal) const \ + { \ + return getValueImpl(model, colPos, comp, parType, bnd, depVal, indepVal); \ + } \ + \ + virtual active getValue(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd, const active& depVal, const active& indepVal) const \ + { \ + return getValueImpl(model, colPos, comp, parType, bnd, depVal, indepVal); \ + } \ } // namespace model } // namespace cadet diff --git a/src/libcadet/model/paramdep/PowerLawParameterDependence.cpp b/src/libcadet/model/paramdep/PowerLawParameterDependence.cpp index 99d9e0d2e..71fc27b9d 100644 --- a/src/libcadet/model/paramdep/PowerLawParameterDependence.cpp +++ b/src/libcadet/model/paramdep/PowerLawParameterDependence.cpp @@ -71,21 +71,15 @@ class PowerLawParameterParameterDependence : public ParameterParameterDependence } template - ParamType getValueImpl(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd) const - { - return 0.0; - } - - template - ParamType getValueImpl(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd, ParamType val) const + ParamType getValueImpl(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd, const ParamType& depVal, const ParamType& indepVal) const { using std::pow; using std::abs; if (_useAbs) - return static_cast(_base) * pow(abs(val), static_cast(_exponent)); + return depVal * static_cast(_base) * pow(abs(indepVal), static_cast(_exponent)); else - return static_cast(_base) * pow(val, static_cast(_exponent)); + return depVal * static_cast(_base) * pow(indepVal, static_cast(_exponent)); } }; diff --git a/src/libcadet/model/parts/AxialConvectionDispersionKernel.hpp b/src/libcadet/model/parts/AxialConvectionDispersionKernel.hpp index 49d9fb31b..30a00aab3 100644 --- a/src/libcadet/model/parts/AxialConvectionDispersionKernel.hpp +++ b/src/libcadet/model/parts/AxialConvectionDispersionKernel.hpp @@ -124,7 +124,7 @@ namespace impl if (cadet_likely(col < p.nCol - 1)) { const double relCoord = static_cast(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)); + const ParamType d_ax_right = p.parDep->getValue(p.model, ColumnPosition{relCoord, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, d_ax, static_cast(p.u)); resBulkComp[col * p.strideCell] -= d_ax_right / h2 * (stencil[1] - stencil[0]); // Jacobian entries if (wantJac) @@ -138,7 +138,7 @@ namespace impl if (cadet_likely(col > 0)) { const double relCoord = static_cast(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)); + const ParamType d_ax_left = p.parDep->getValue(p.model, ColumnPosition{relCoord, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, d_ax, static_cast(p.u)); resBulkComp[col * p.strideCell] -= d_ax_left / h2 * (stencil[-1] - stencil[0]); // Jacobian entries if (wantJac) @@ -271,7 +271,7 @@ namespace impl if (cadet_likely(col < p.nCol - 1)) { const double relCoord = static_cast(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)); + const ParamType d_ax_right = p.parDep->getValue(p.model, ColumnPosition{relCoord, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, d_ax, static_cast(p.u)); resBulkComp[col * p.strideCell] -= d_ax_right / h2 * (stencil[-1] - stencil[0]); // Jacobian entries if (wantJac) @@ -285,7 +285,7 @@ namespace impl if (cadet_likely(col > 0)) { const double relCoord = static_cast(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)); + const ParamType d_ax_left = p.parDep->getValue(p.model, ColumnPosition{relCoord, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, d_ax, static_cast(p.u)); resBulkComp[col * p.strideCell] -= d_ax_left / h2 * (stencil[1] - stencil[0]); // Jacobian entries if (wantJac) diff --git a/src/libcadet/model/parts/ConvectionDispersionOperator.cpp b/src/libcadet/model/parts/ConvectionDispersionOperator.cpp index 0e508127d..82017c919 100644 --- a/src/libcadet/model/parts/ConvectionDispersionOperator.cpp +++ b/src/libcadet/model/parts/ConvectionDispersionOperator.cpp @@ -78,7 +78,7 @@ bool AxialConvectionDispersionOperatorBase::configureModelDiscretization(IParame _dispersionDep->configureModelDiscretization(paramProvider); } else - _dispersionDep = helper.createParameterParameterDependence("CONSTANT_ONE"); + _dispersionDep = helper.createParameterParameterDependence("IDENTITY"); paramProvider.pushScope("discretization"); @@ -577,7 +577,7 @@ bool RadialConvectionDispersionOperatorBase::configureModelDiscretization(IParam _dispersionDep->configureModelDiscretization(paramProvider); } else - _dispersionDep = helper.createParameterParameterDependence("CONSTANT_ONE"); + _dispersionDep = helper.createParameterParameterDependence("IDENTITY"); paramProvider.pushScope("discretization"); diff --git a/src/libcadet/model/parts/RadialConvectionDispersionKernel.hpp b/src/libcadet/model/parts/RadialConvectionDispersionKernel.hpp index 33e8b5f23..30e240798 100644 --- a/src/libcadet/model/parts/RadialConvectionDispersionKernel.hpp +++ b/src/libcadet/model/parts/RadialConvectionDispersionKernel.hpp @@ -126,7 +126,7 @@ namespace impl if (cadet_likely(col < p.nCol - 1)) { const double relCoord = (static_cast(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])); + const ParamType d_rad_right = p.parDep->getValue(p.model, ColumnPosition{relCoord, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, d_rad, 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) @@ -142,11 +142,11 @@ namespace impl { const double relCoord = (static_cast(p.cellBounds[col]) - static_cast(p.cellBounds[0])) / (static_cast(p.cellBounds[p.nCol - 1]) - static_cast(p.cellBounds[0])); const ParamType d_rad_left = d_rad * p.parDep->getValue(p.model, ColumnPosition{ relCoord, 0.0, 0.0 }, comp, ParTypeIndep, BoundStateIndep, static_cast(p.u) / static_cast(p.cellBounds[col])); - resBulkComp[col * p.strideCell] -= d_rad_left * static_cast(p.cellBounds[col]) / denom * (stencil[-1] - stencil[0]) / (static_cast(p.cellCenters[col]) - static_cast(p.cellCenters[col-1])); + 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_left) * static_cast(p.cellBounds[col]) / static_cast(denom) / (static_cast(p.cellCenters[col]) - static_cast(p.cellCenters[col-1])); + 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; } @@ -284,7 +284,7 @@ namespace impl if (cadet_likely(col < p.nCol - 1)) { 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])); + const ParamType d_rad_right = p.parDep->getValue(p.model, ColumnPosition{relCoord, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, d_rad, 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) @@ -299,7 +299,7 @@ namespace impl if (cadet_likely(col > 0)) { 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])); + const ParamType d_rad_left = p.parDep->getValue(p.model, ColumnPosition{relCoord, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, d_rad, 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)