Skip to content

Commit

Permalink
Change return logic of parameter parameter dependence
Browse files Browse the repository at this point in the history 
Logic changed such that the parameter value itself is returned, instead of a modifying factor

Co-authored-by: schmoelder <j.schmoelder@fz-juelich.de>
  • Loading branch information
@jbreue16 @schmoelder
jbreue16 and schmoelder committed Apr 19, 2024
1 parent 484dd54 commit 2c63233
Show file tree
Hide file tree
Showing 12 changed files with 81 additions and 315 deletions.
1 change: 0 additions & 1 deletion src/libcadet/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
)

Expand Down
4 changes: 0 additions & 4 deletions src/libcadet/ParameterDependenceFactory.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,8 +22,6 @@ namespace cadet
void registerLiquidSaltSolidParamDependence(std::unordered_map<std::string, std::function<model::IParameterStateDependence*()>>& paramDeps);
void registerDummyParamDependence(std::unordered_map<std::string, std::function<model::IParameterStateDependence*()>>& paramDeps);
void registerDummyParamDependence(std::unordered_map<std::string, std::function<model::IParameterParameterDependence*()>>& paramDeps);
void registerIdentityParamDependence(std::unordered_map<std::string, std::function<model::IParameterStateDependence*()>>& paramDeps);
void registerIdentityParamDependence(std::unordered_map<std::string, std::function<model::IParameterParameterDependence*()>>& paramDeps);
void registerPowerLawParamDependence(std::unordered_map<std::string, std::function<model::IParameterParameterDependence*()>>& paramDeps);
}
}
Expand All @@ -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);
}

Expand Down
22 changes: 11 additions & 11 deletions src/libcadet/model/GeneralRateModel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -459,7 +459,7 @@ bool GeneralRateModel<ConvDispOperator>::configureModelDiscretization(IParameter
_filmDiffDep->configureModelDiscretization(paramProvider);
}
else
_filmDiffDep = helper.createParameterParameterDependence("CONSTANT_ONE");
_filmDiffDep = helper.createParameterParameterDependence("IDENTITY");

// ==== Construct and configure binding model
clearBindingModels();
Expand Down Expand Up @@ -1926,11 +1926,11 @@ int GeneralRateModel<ConvDispOperator>::residualFlux(double t, unsigned int secI

const double relPos = _convDispOp.relativeCoordinate(colCell);
const active curVelocity = _convDispOp.currentVelocity(relPos);
const active modifier = _filmDiffDep->getValue(*this, ColumnPosition{ relPos, 0.0, 0.0 }, comp, ParTypeIndep, BoundStateIndep, curVelocity);
const ParamType filmDiffDep = static_cast<ParamType>(_filmDiffDep->getValue(*this, ColumnPosition{ relPos, 0.0, 0.0 }, comp, ParTypeIndep, BoundStateIndep, filmDiff[comp], curVelocity));
if (cadet_likely(_colParBoundaryOrder == 2))
kf_FV[comp] = 1.0 / (absOuterShellHalfRadius / epsP / static_cast<ParamType>(_poreAccessFactor[type * _disc.nComp + comp]) / static_cast<ParamType>(parDiff[comp]) + 1.0 / (static_cast<ParamType>(filmDiff[comp]) * static_cast<ParamType>(modifier)));
kf_FV[comp] = 1.0 / (absOuterShellHalfRadius / epsP / static_cast<ParamType>(_poreAccessFactor[type * _disc.nComp + comp]) / static_cast<ParamType>(parDiff[comp]) + 1.0 / filmDiffDep);
else
kf_FV[comp] = static_cast<ParamType>(filmDiff[comp]) * static_cast<ParamType>(modifier);
kf_FV[comp] = filmDiffDep;
}

// J_{0,f} block, adds flux to column void / bulk volume equations
Expand Down Expand Up @@ -1974,9 +1974,9 @@ int GeneralRateModel<ConvDispOperator>::residualFlux(double t, unsigned int secI

const double relPos = _convDispOp.relativeCoordinate(pblk);
const active curVelocity = _convDispOp.currentVelocity(relPos);
const active modifier = _filmDiffDep->getValue(*this, ColumnPosition{ relPos, 0.0, 0.0 }, comp, ParTypeIndep, BoundStateIndep, curVelocity);
const ParamType filmDiffDep = static_cast<ParamType>(_filmDiffDep->getValue(*this, ColumnPosition{ relPos, 0.0, 0.0 }, comp, ParTypeIndep, BoundStateIndep, filmDiff[comp], curVelocity));

kf_FV[comp] = (1.0 - static_cast<ParamType>(_parPorosity[type])) / (1.0 + epsP * static_cast<ParamType>(_poreAccessFactor[type * _disc.nComp + comp]) * static_cast<ParamType>(parDiff[comp]) / (absOuterShellHalfRadius * static_cast<ParamType>(filmDiff[comp]) * static_cast<ParamType>(modifier)));
kf_FV[comp] = (1.0 - static_cast<ParamType>(_parPorosity[type])) / (1.0 + epsP * static_cast<ParamType>(_poreAccessFactor[type * _disc.nComp + comp]) * static_cast<ParamType>(parDiff[comp]) / (absOuterShellHalfRadius * filmDiffDep));
}

const ColumnPosition colPos{(0.5 + static_cast<double>(pblk)) / static_cast<double>(_disc.nCol), 0.0, static_cast<double>(parCenterRadius[0]) / static_cast<double>(_parRadius[type])};
Expand Down Expand Up @@ -2089,12 +2089,12 @@ void GeneralRateModel<ConvDispOperator>::assembleOffdiagJac(double t, unsigned i

const double relPos = _convDispOp.relativeCoordinate(pblk);
const double curVelocity = static_cast<double>(_convDispOp.currentVelocity(relPos));
const double modifier = _filmDiffDep->getValue(*this, ColumnPosition{ relPos, 0.0, 0.0 }, comp, ParTypeIndep, BoundStateIndep, curVelocity);
const double filmDiffDep = static_cast<double>(_filmDiffDep->getValue(*this, ColumnPosition{ relPos, 0.0, 0.0 }, comp, ParTypeIndep, BoundStateIndep, filmDiff[comp], curVelocity));

if (cadet_likely(_colParBoundaryOrder == 2))
kf_FV[comp] = 1.0 / (absOuterShellHalfRadius / epsP / static_cast<double>(_poreAccessFactor[type * _disc.nComp + comp]) / static_cast<double>(parDiff[comp]) + 1.0 / (static_cast<double>(filmDiff[comp]) * modifier));
kf_FV[comp] = 1.0 / (absOuterShellHalfRadius / epsP / static_cast<double>(_poreAccessFactor[type * _disc.nComp + comp]) / static_cast<double>(parDiff[comp]) + 1.0 / filmDiffDep);
else
kf_FV[comp] = static_cast<double>(filmDiff[comp]) * modifier;
kf_FV[comp] = filmDiffDep;
}

// J_{0,f} block, adds flux to column void / bulk volume equations
Expand Down Expand Up @@ -2148,9 +2148,9 @@ void GeneralRateModel<ConvDispOperator>::assembleOffdiagJac(double t, unsigned i
{
const double relPos = _convDispOp.relativeCoordinate(pblk);
const double curVelocity = static_cast<double>(_convDispOp.currentVelocity(relPos));
const double modifier = _filmDiffDep->getValue(*this, ColumnPosition{ relPos, 0.0, 0.0 }, comp, ParTypeIndep, BoundStateIndep, curVelocity);
const double filmDiffDep = static_cast<double>(_filmDiffDep->getValue(*this, ColumnPosition{ relPos, 0.0, 0.0 }, comp, ParTypeIndep, BoundStateIndep, filmDiff[comp], curVelocity));

kf_FV[comp] = (1.0 - static_cast<double>(_parPorosity[type])) / (1.0 + epsP * static_cast<double>(_poreAccessFactor[type * _disc.nComp + comp]) * static_cast<double>(parDiff[comp]) / (absOuterShellHalfRadius * static_cast<double>(filmDiff[comp]) * modifier));
kf_FV[comp] = (1.0 - static_cast<double>(_parPorosity[type])) / (1.0 + epsP * static_cast<double>(_poreAccessFactor[type * _disc.nComp + comp]) * static_cast<double>(parDiff[comp]) / (absOuterShellHalfRadius * filmDiffDep));
}

const ColumnPosition colPos{ (0.5 + static_cast<double>(pblk)) / static_cast<double>(_disc.nCol), 0.0, static_cast<double>(parCenterRadius[0]) / static_cast<double>(_parRadius[type]) };
Expand Down
18 changes: 9 additions & 9 deletions src/libcadet/model/LumpedRateModelWithPores.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -275,7 +275,7 @@ bool LumpedRateModelWithPores<ConvDispOperator>::configureModelDiscretization(IP
_filmDiffDep->configureModelDiscretization(paramProvider);
}
else
_filmDiffDep = helper.createParameterParameterDependence("CONSTANT_ONE");
_filmDiffDep = helper.createParameterParameterDependence("IDENTITY");

// Allocate memory
Indexer idxr(_disc);
Expand Down Expand Up @@ -1091,9 +1091,9 @@ int LumpedRateModelWithPores<ConvDispOperator>::residualFlux(double t, unsigned

const double relPos = _convDispOp.relativeCoordinate(colCell);
const active curVelocity = _convDispOp.currentVelocity(relPos);
const active modifier = _filmDiffDep->getValue(*this, ColumnPosition{relPos, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, curVelocity);
const ParamType filmDiffDep = static_cast<ParamType>(_filmDiffDep->getValue(*this, ColumnPosition{relPos, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, filmDiff[comp], curVelocity));

resCol[i] += jacCF_val * static_cast<ParamType>(filmDiff[comp]) * static_cast<ParamType>(modifier) * static_cast<ParamType>(_parTypeVolFrac[type + _disc.nParType * colCell]) * yFluxType[i];
resCol[i] += jacCF_val * filmDiffDep * static_cast<ParamType>(_parTypeVolFrac[type + _disc.nParType * colCell]) * yFluxType[i];
}

// J_{f,0} block, adds bulk volume state c_i to flux equation
Expand All @@ -1114,10 +1114,10 @@ int LumpedRateModelWithPores<ConvDispOperator>::residualFlux(double t, unsigned

for (unsigned int comp = 0; comp < _disc.nComp; ++comp)
{
const active modifier = _filmDiffDep->getValue(*this, ColumnPosition{relPos, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, curVelocity);
const ParamType filmDiffDep = static_cast<ParamType>(_filmDiffDep->getValue(*this, ColumnPosition{relPos, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, filmDiff[comp], curVelocity));

const unsigned int eq = pblk * idxr.strideColCell() + comp * idxr.strideColComp();
resParType[pblk * idxr.strideParBlock(type) + comp] += jacPF_val / static_cast<ParamType>(poreAccFactor[comp]) * static_cast<ParamType>(filmDiff[comp]) * static_cast<ParamType>(modifier) * yFluxType[eq];
resParType[pblk * idxr.strideParBlock(type) + comp] += jacPF_val / static_cast<ParamType>(poreAccFactor[comp]) * filmDiffDep * yFluxType[eq];
}
}

Expand Down Expand Up @@ -1180,10 +1180,10 @@ void LumpedRateModelWithPores<ConvDispOperator>::assembleOffdiagJac(double t, un

const double relPos = _convDispOp.relativeCoordinate(colCell);
const double curVelocity = static_cast<double>(_convDispOp.currentVelocity(relPos));
const double modifier = _filmDiffDep->getValue(*this, ColumnPosition{relPos, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, curVelocity);
const double filmDiffDep = static_cast<double>(_filmDiffDep->getValue(*this, ColumnPosition{relPos, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, filmDiff[comp], curVelocity));

// Main diagonal corresponds to j_{f,i} (flux) state variable
_jacCF.addElement(eq, eq + typeOffset, jacCF_val * static_cast<double>(filmDiff[comp]) * modifier * static_cast<double>(_parTypeVolFrac[type + _disc.nParType * colCell]));
_jacCF.addElement(eq, eq + typeOffset, jacCF_val * filmDiffDep * static_cast<double>(_parTypeVolFrac[type + _disc.nParType * colCell]));
}

// J_{f,0} block, adds bulk volume state c_i to flux equation
Expand All @@ -1203,8 +1203,8 @@ void LumpedRateModelWithPores<ConvDispOperator>::assembleOffdiagJac(double t, un
const unsigned int eq = typeOffset + pblk * idxr.strideColCell() + comp * idxr.strideColComp();
const unsigned int col = pblk * idxr.strideParBlock(type) + comp;

const double modifier = _filmDiffDep->getValue(*this, ColumnPosition{relPos, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, curVelocity);
_jacPF[type].addElement(col, eq, jacPF_val / static_cast<double>(poreAccFactor[comp]) * static_cast<double>(filmDiff[comp]) * modifier);
const double filmDiffDep = static_cast<double>(_filmDiffDep->getValue(*this, ColumnPosition{relPos, 0.0, 0.0}, comp, ParTypeIndep, BoundStateIndep, filmDiff[comp], curVelocity));
_jacPF[type].addElement(col, eq, jacPF_val / static_cast<double>(poreAccFactor[comp]) * filmDiffDep);
}
}

Expand Down
36 changes: 6 additions & 30 deletions src/libcadet/model/ParameterDependence.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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:
};
Expand Down
111 changes: 32 additions & 79 deletions src/libcadet/model/paramdep/DummyParameterDependence.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,36 @@ namespace cadet
namespace model
{

/**
* @brief Defines a dummy parameter dependence that outputs the (independent) parameter itself
*/
class IdentityParameterParameterDependence : public ParameterParameterDependenceBase
{
public:

IdentityParameterParameterDependence() { }
virtual ~IdentityParameterParameterDependence() CADET_NOEXCEPT { }

static const char* identifier() { return "IDENTITY"; }
virtual const char* name() const CADET_NOEXCEPT { return IdentityParameterParameterDependence::identifier(); }

CADET_PARAMETERPARAMETERDEPENDENCE_BOILERPLATE

protected:

virtual bool configureImpl(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, BoundStateIdx bndIdx, const std::string& name)
{
return true;
}

template <typename ParamType>
ParamType getValueImpl(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd, const ParamType& depVal, const ParamType& indepVal) const
{
return depVal;
}

};

/**
* @brief Defines a parameter dependence that outputs constant 0.0
*/
Expand Down Expand Up @@ -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
*/
Expand Down Expand Up @@ -141,81 +170,6 @@ class ConstantOneParameterStateDependence : public ParameterStateDependenceBase
void analyticJacobianCombinedAddSolidImpl(const ColumnPosition& colPos, double param, double const* yLiquid, double const* ySolid, int bnd, double factor, int offset, RowIterator jac) const { }
};


/**
* @brief Defines a parameter dependence that outputs constant 0.0
*/
class ConstantZeroParameterParameterDependence : public ParameterParameterDependenceBase
{
public:

ConstantZeroParameterParameterDependence() { }
virtual ~ConstantZeroParameterParameterDependence() CADET_NOEXCEPT { }

static const char* identifier() { return "CONSTANT_ZERO"; }
virtual const char* name() const CADET_NOEXCEPT { return ConstantZeroParameterParameterDependence::identifier(); }

CADET_PARAMETERPARAMETERDEPENDENCE_BOILERPLATE

protected:

virtual bool configureImpl(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, BoundStateIdx bndIdx, const std::string& name)
{
return true;
}

template <typename ParamType>
ParamType getValueImpl(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd) const
{
return 0.0;
}

template <typename ParamType>
ParamType getValueImpl(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd, const ParamType& val) const
{
return 0.0;
}

};


/**
* @brief Defines a parameter dependence that outputs constant 1.0
*/
class ConstantOneParameterParameterDependence : public ParameterParameterDependenceBase
{
public:

ConstantOneParameterParameterDependence() { }
virtual ~ConstantOneParameterParameterDependence() CADET_NOEXCEPT { }

static const char* identifier() { return "CONSTANT_ONE"; }
virtual const char* name() const CADET_NOEXCEPT { return ConstantOneParameterParameterDependence::identifier(); }

CADET_PARAMETERPARAMETERDEPENDENCE_BOILERPLATE

protected:

virtual bool configureImpl(IParameterProvider& paramProvider, UnitOpIdx unitOpIdx, ParticleTypeIdx parTypeIdx, BoundStateIdx bndIdx, const std::string& name)
{
return true;
}

template <typename ParamType>
ParamType getValueImpl(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd) const
{
return 1.0;
}

template <typename ParamType>
ParamType getValueImpl(const IModel& model, const ColumnPosition& colPos, int comp, int parType, int bnd, const ParamType& val) const
{
return 1.0;
}

};


namespace paramdep
{
void registerDummyParamDependence(std::unordered_map<std::string, std::function<model::IParameterStateDependence*()>>& paramDeps)
Expand All @@ -227,9 +181,8 @@ namespace paramdep

void registerDummyParamDependence(std::unordered_map<std::string, std::function<model::IParameterParameterDependence*()>>& paramDeps)
{
paramDeps[ConstantOneParameterParameterDependence::identifier()] = []() { return new ConstantOneParameterParameterDependence(); };
paramDeps[ConstantZeroParameterParameterDependence::identifier()] = []() { return new ConstantZeroParameterParameterDependence(); };
paramDeps["NONE"] = []() { return new ConstantOneParameterParameterDependence(); };
paramDeps[IdentityParameterParameterDependence::identifier()] = []() { return new IdentityParameterParameterDependence(); };
paramDeps["NONE"] = []() { return new IdentityParameterParameterDependence(); };
}
} // namespace paramdep

Expand Down

0 comments on commit 2c63233

Please sign in to comment.