Skip to content

add more parameter dependencies #189

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open

add more parameter dependencies #189

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
b6c98d7
Add film diffusion parameter dependence to GRM
jbreue16 May 9, 2023
f33ab8d
Change return logic of parameter parameter dependence
jbreue16 May 26, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 0 additions & 1 deletion src/libcadet/CMakeLists.txt
View file
Open in desktop
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
Open in desktop
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
198 changes: 124 additions & 74 deletions src/libcadet/model/GeneralRateModel.cpp
View file
Open in desktop

Large diffs are not rendered by default.

1 change: 1 addition & 0 deletions src/libcadet/model/GeneralRateModel.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
18 changes: 9 additions & 9 deletions src/libcadet/model/LumpedRateModelWithPores.cpp
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Loading
Loading