Skip to content

Commit

Permalink
return parameter value (instead of modifying factor) in parameter par…
Browse files Browse the repository at this point in the history 
…ameter dependence
  • Loading branch information
@jbreue16 @schmoelder
jbreue16 authored and schmoelder committed Jul 21, 2023
1 parent 866db2c commit 68801ae
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 @@ -440,7 +440,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 @@ -1907,11 +1907,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 @@ -1955,9 +1955,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 @@ -2070,12 +2070,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 @@ -2129,9 +2129,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 @@ -249,7 +249,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 @@ -1065,9 +1065,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 @@ -1088,10 +1088,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 @@ -1154,10 +1154,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 @@ -1177,8 +1177,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 68801ae

Please sign in to comment.