Skip to content
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.

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

Merge GasKinetics into BulkKinetics and simplify rate updates #1483

Merged
merged 8 commits into from
Apr 21, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
2 changes: 2 additions & 0 deletions include/cantera/base/ExtensionManager.h
Original file line number Diff line number Diff line change
Expand Up @@ -92,6 +92,8 @@ class ExtensionManager
//! @param rateName The name of the reaction rate type
//! @param wrapperName The name used for Solution wrappers to be used with this
//! object, corresponding to a type registered with registerSolutionLinker().
//! @param link Function that creates ReactionData wrapper and links it to the
//! provided C++ object
static void registerReactionDataLinker(const string& rateName,
const string& wrapperName, function<void(ReactionDataDelegator&)> link);

Expand Down
1 change: 1 addition & 0 deletions include/cantera/base/SolutionArray.h
Original file line number Diff line number Diff line change
Expand Up @@ -285,6 +285,7 @@ class SolutionArray
*
* @param fname Name of container file (YAML or HDF)
* @param id Identifier of SolutionArray within the container file
* @param sub Name of the subgroup holding actual data
*/
AnyMap restore(const string& fname, const string& id, const string& sub);

Expand Down
1 change: 1 addition & 0 deletions include/cantera/clib/ct.h
Original file line number Diff line number Diff line change
Expand Up @@ -122,6 +122,7 @@ extern "C" {
CANTERA_CAPI int thermo_setState_Psat(int n, double p, double x);
CANTERA_CAPI int thermo_setState_Tsat(int n, double t, double x);

//! @since Starting in Cantera 3.0, the "phasename" argument should be blank
CANTERA_CAPI int kin_newFromFile(const char* filename, const char* phasename,
int reactingPhase, int neighbor1, int neighbor2,
int neighbor3, int neighbor4);
Expand Down
9 changes: 7 additions & 2 deletions include/cantera/kinetics/Arrhenius.h
Original file line number Diff line number Diff line change
Expand Up @@ -71,15 +71,20 @@ class ArrheniusBase : public ReactionRate

//! Perform object setup based on AnyValue node information
/*!
* @param rate AnyValue containing rate information
* Used to set parameters from a child of the reaction node, which may have
* different names for different rate parameterizations, such as falloff rates.
*
* @param rate Child of the reaction node containing Arrhenius rate parameters.
* For example, the `rate-coefficient` node for a standard Arrhenius reaction.
* @param units Unit system
* @param rate_units Unit definitions specific to rate information
*/
void setRateParameters(const AnyValue& rate,
const UnitSystem& units,
const UnitStack& rate_units);

//! Return parameters
//! Get Arrhenius parameters used to populate the `rate-coefficient` or
//! equivalent field
void getRateParameters(AnyMap& node) const;

virtual void setParameters(
Expand Down
170 changes: 141 additions & 29 deletions include/cantera/kinetics/BulkKinetics.h
Original file line number Diff line number Diff line change
Expand Up @@ -16,47 +16,148 @@
namespace Cantera
{

//! Partial specialization of Kinetics for chemistry in a single bulk phase
//! Specialization of Kinetics for chemistry in a single bulk phase
//! @ingroup kinetics
class BulkKinetics : public Kinetics
{
public:
BulkKinetics() = default;
//! @name Constructors and General Information
//! @{
BulkKinetics();

//! @deprecated To be removed after Cantera 3.0; code base only uses default.
BulkKinetics(ThermoPhase* thermo);

virtual void resizeReactions();
string kineticsType() const override {
return "bulk";
}

virtual bool isReversible(size_t i);

virtual void getDeltaGibbs(doublereal* deltaG);
virtual void getDeltaEnthalpy(doublereal* deltaH);
virtual void getDeltaEntropy(doublereal* deltaS);

virtual void getDeltaSSGibbs(doublereal* deltaG);
virtual void getDeltaSSEnthalpy(doublereal* deltaH);
virtual void getDeltaSSEntropy(doublereal* deltaS);

virtual void getRevRateConstants(double* krev,
bool doIrreversible = false);

virtual bool addReaction(shared_ptr<Reaction> r, bool resize=true);
virtual void modifyReaction(size_t i, shared_ptr<Reaction> rNew);

virtual void resizeSpecies();

virtual void setMultiplier(size_t i, double f);
virtual void invalidateCache();
bool isReversible(size_t i) override;
//! @}

//! @name Reaction Mechanism Setup Routines
//! @{
bool addReaction(shared_ptr<Reaction> r, bool resize=true) override;
void addThirdBody(shared_ptr<Reaction> r);
void modifyReaction(size_t i, shared_ptr<Reaction> rNew) override;
void resizeSpecies() override;
void resizeReactions() override;
void setMultiplier(size_t i, double f) override;
void invalidateCache() override;
//! @}

//! @name Reaction rate constants, rates of progress, and thermodynamic properties
//! @{
void getFwdRateConstants(double* kfwd) override;
void getEquilibriumConstants(double* kc) override;
void getRevRateConstants(double* krev, bool doIrreversible=false) override;

void getDeltaGibbs(double* deltaG) override;
void getDeltaEnthalpy(double* deltaH) override;
void getDeltaEntropy(double* deltaS) override;

void getDeltaSSGibbs(double* deltaG) override;
void getDeltaSSEnthalpy(double* deltaH) override;
void getDeltaSSEntropy(double* deltaS) override;
//! @}

//! @name Derivatives of rate constants and rates of progress
//! @{
void getDerivativeSettings(AnyMap& settings) const override;
void setDerivativeSettings(const AnyMap& settings) override;
void getFwdRateConstants_ddT(double* dkfwd) override;
void getFwdRatesOfProgress_ddT(double* drop) override;
void getRevRatesOfProgress_ddT(double* drop) override;
void getNetRatesOfProgress_ddT(double* drop) override;
void getFwdRateConstants_ddP(double* dkfwd) override;
void getFwdRatesOfProgress_ddP(double* drop) override;
void getRevRatesOfProgress_ddP(double* drop) override;
void getNetRatesOfProgress_ddP(double* drop) override;
void getFwdRateConstants_ddC(double* dkfwd) override;
void getFwdRatesOfProgress_ddC(double* drop) override;
void getRevRatesOfProgress_ddC(double* drop) override;
void getNetRatesOfProgress_ddC(double* drop) override;
Eigen::SparseMatrix<double> fwdRatesOfProgress_ddX() override;
Eigen::SparseMatrix<double> revRatesOfProgress_ddX() override;
Eigen::SparseMatrix<double> netRatesOfProgress_ddX() override;
Eigen::SparseMatrix<double> fwdRatesOfProgress_ddCi() override;
Eigen::SparseMatrix<double> revRatesOfProgress_ddCi() override;
Eigen::SparseMatrix<double> netRatesOfProgress_ddCi() override;
//! @}

//! @name Rate calculation intermediate methods
//! @{

void updateROP() override;

void getThirdBodyConcentrations(double* concm) override;
const vector_fp& thirdBodyConcentrations() const override {
return m_concm;
}

//! @}

protected:
//! @name Internal service methods
//!
//! @note These methods are for internal use, and seek to avoid code duplication
//! while evaluating terms used for rate constants, rates of progress, and
//! their derivatives.
//! @{

//! Multiply rate with third-body collider concentrations
void processThirdBodies(double* rop);

//! Multiply rate with inverse equilibrium constant
void applyEquilibriumConstants(double* rop);

//! Multiply rate with scaled temperature derivatives of the inverse
//! equilibrium constant
/*!
* This (scaled) derivative is handled by a finite difference.
*/
void applyEquilibriumConstants_ddT(double* drkcn);

//! Process temperature derivative
//! @param in rate expression used for the derivative calculation
//! @param drop pointer to output buffer
void process_ddT(const vector_fp& in, double* drop);

//! Process pressure derivative
//! @param in rate expression used for the derivative calculation
//! @param drop pointer to output buffer
void process_ddP(const vector_fp& in, double* drop);

//! Process concentration (molar density) derivative
//! @param stoich stoichiometry manager
//! @param in rate expression used for the derivative calculation
//! @param drop pointer to output buffer
//! @param mass_action boolean indicating whether law of mass action applies
void process_ddC(StoichManagerN& stoich, const vector_fp& in,
double* drop, bool mass_action=true);

//! Process derivatives
//! @param stoich stoichiometry manager
//! @param in rate expression used for the derivative calculation
//! @param ddX true: w.r.t mole fractions false: w.r.t species concentrations
//! @return a sparse matrix of derivative contributions for each reaction of
//! dimensions nTotalReactions by nTotalSpecies
Eigen::SparseMatrix<double> calculateCompositionDerivatives(
StoichManagerN& stoich, const vector_fp& in, bool ddX=true);

//! Helper function ensuring that all rate derivatives can be calculated
//! @param name method name used for error output
//! @throw CanteraError if ideal gas assumption does not hold
void assertDerivativesValid(const string& name);

//! @}

//! Vector of rate handlers
std::vector<unique_ptr<MultiRateBase>> m_bulk_rates;
std::map<std::string, size_t> m_bulk_types; //!< Mapping of rate handlers
vector<unique_ptr<MultiRateBase>> m_bulk_rates;
map<string, size_t> m_bulk_types; //!< Mapping of rate handlers

std::vector<size_t> m_revindex; //!< Indices of reversible reactions
std::vector<size_t> m_irrev; //!< Indices of irreversible reactions
vector<size_t> m_revindex; //!< Indices of reversible reactions
vector<size_t> m_irrev; //!< Indices of irreversible reactions

//! Difference between the global reactants order and the global products
//! order. Of type "double" to account for the fact that we can have real-
Expand All @@ -74,10 +175,21 @@ class BulkKinetics : public Kinetics
//! Physical concentrations, as calculated by ThermoPhase::getConcentrations
vector_fp m_phys_conc;

vector_fp m_grt;
//! Derivative settings
bool m_jac_skip_third_bodies;
bool m_jac_skip_falloff;
double m_jac_rtol_delta;

bool m_ROP_ok = false;
double m_temp = 0.0;

//! Buffers for partial rop results with length nReactions()
vector_fp m_rbuf0;
vector_fp m_rbuf1;
vector_fp m_rbuf2;
vector_fp m_kf0; //!< Forward rate constants without perturbation
vector_fp m_sbuf0;
vector_fp m_state;
vector_fp m_grt; //!< Standard chemical potentials for each species
};

}
Expand Down
Loading