ReactionNetworkUtils
commentnote
For the [ReactionNetwork] syntax parser, see this page instead. This utility is currently used by the [Physics] syntax.
Reaction network parser using a Parsed Expression Grammar
The chemical reactions module includes a reaction network parser this utility that enables chemical reactions to be specified in a natural form in the input file. This parser is notably used by the AqueousReactionsEquilibriumPhysics.
The input file syntax for reactions has to be written in the following form:
Individual reactions are provided with the reactant species on the left hand side, while the product species follow the -> sign, followed by metadata in between square brackets. A linebreak is used to delimit reactions, so that multiple reactions can be entered.
commentnote
The AqueousReactionsEquilibriumPhysics currently only supports one product species per reaction, on the right hand side, as shown in the second equation in the example.
are coefficients which have to be floating point number in the regular non-scientific notation. Metadata keys and values are parsed as strings. The values may be floating point numbers or names of properties, though AqueousReactionsEquilibriumPhysics only supports numbers at this time. The parsing of the metadata is implemented in the classes using this utility.
Troubleshooting
If the reaction network parsing fails, it will error with:
mooseError("Failed to parse reaction.");
In that case, you may examine the expected syntax below. The grammar is used by the parser to understand the input. *, +, ? and other symbols use the same rules as for regular expression parsing.
parser parser(R"(
Reactions <- (_ Reaction _ Newline?)* _
Reaction <- TermList _ '->' _ TermList _ Metadata?
TermList <- Term (_ Term)*
Term <- Coefficient? Species
Coefficient <- Float
Float <- [+-]? (_)? ([0-9]+)? ('.' [0-9]+)?
Species <- BaseSpecies State? Charge?
BaseSpecies <- [a-zA-Z][a-zA-Z0-9_]*
State <- '(' [a-z]{1,3} ')'
Charge <- ('^'? [0-9]* [+-])
Metadata <- '[' _ Pair (_ ',' _ Pair)* _ ']'
Pair <- Key _ '=' _ Value
Key <- [A-Za-z_][A-Za-z0-9_]*
Value <- [A-Za-z0-9_.-]+
Newline <- [\r\n]+
~_ <- [ \t]*
)");
(modules/chemical_reactions/src/utils/ReactionNetworkUtils.C)
// This file is part of the MOOSE framework
// https://mooseframework.inl.gov
//
// All rights reserved, see COPYRIGHT for full restrictions
// https://github.com/idaholab/moose/blob/master/COPYRIGHT
//
// Licensed under LGPL 2.1, please see LICENSE for details
// https://www.gnu.org/licenses/lgpl-2.1.html
#include "ReactionNetworkUtils.h"
#include "MooseTypes.h"
#include "MooseUtils.h"
#include "libmesh/utility.h"
#include "peglib.h"
// not needed
#include "Conversion.h"
namespace ReactionNetworkUtils
{
std::vector<Reaction>
parseReactionNetwork(const std::string & reaction_network_string, bool output_to_cout)
{
using namespace peg;
parser parser(R"(
Reactions <- (_ Reaction _ Newline?)* _
Reaction <- TermList _ '->' _ TermList _ Metadata?
TermList <- Term (_ Term)*
Term <- Coefficient? Species
Coefficient <- Float
Float <- [+-]? (_)? ([0-9]+)? ('.' [0-9]+)?
Species <- BaseSpecies State? Charge?
BaseSpecies <- [a-zA-Z][a-zA-Z0-9_]*
State <- '(' [a-z]{1,3} ')'
Charge <- ('^'? [0-9]* [+-])
Metadata <- '[' _ Pair (_ ',' _ Pair)* _ ']'
Pair <- Key _ '=' _ Value
Key <- [A-Za-z_][A-Za-z0-9_]*
Value <- [A-Za-z0-9_.-]+
Newline <- [\r\n]+
~_ <- [ \t]*
)");
mooseAssert(parser, "Failed to build parser");
// Rules for how to parse each type
parser["Reactions"] = [](const SemanticValues & vs)
{
std::vector<Reaction> reactions;
for (const auto & v : vs)
{
if (v.type() == typeid(Reaction))
{
reactions.push_back(std::any_cast<Reaction>(v));
}
}
return reactions;
};
parser["Float"] = [](const SemanticValues & vs)
{
// We need to remove whitespace in case the reaction was A - B -> C
std::string vs_str = std::string(std::any_cast<std::string_view>(vs.token()));
vs_str.erase(std::remove(vs_str.begin(), vs_str.end(), ' '), vs_str.end());
if (vs_str == "-")
return -1.;
else if (vs_str == "+" || vs_str == "")
return 1.;
else
return std::stod(std::string(vs_str));
};
parser["Coefficient"] = [](const SemanticValues & vs) { return std::any_cast<double>(vs[0]); };
parser["BaseSpecies"] = [](const SemanticValues & vs) { return vs.token(); };
parser["State"] = [](const SemanticValues & vs) { return vs.token(); };
parser["Charge"] = [](const SemanticValues & vs) { return vs.token(); };
parser["Species"] = [](const SemanticValues & vs)
{
std::string base = std::string(std::any_cast<std::string_view>(vs[0]));
std::optional<std::string> state;
std::optional<std::string> charge;
// Parse the state and charge
if (vs.size() >= 2)
{
for (size_t i = 1; i < vs.size(); ++i)
if (vs[i].type() == typeid(std::string_view))
{
std::string val = std::string(std::any_cast<std::string_view>(vs[i]));
if (val.front() == '(' && val.back() == ')')
state = val;
else
charge = val;
}
}
return Term{1.0, base, state, charge}; // default coefficient = 1.0
};
parser["Term"] = [](const SemanticValues & vs)
{
mooseAssert(vs.size() <= 2, "Should only have a coefficient and a species");
if (vs.size() == 2)
{
double coeff = std::any_cast<double>(vs[0]);
Term term = std::any_cast<Term>(vs[1]);
term.coefficient = coeff;
return term;
}
else
return std::any_cast<Term>(vs[0]);
};
parser["TermList"] = [](const SemanticValues & vs)
{
TermList terms;
for (const auto & v : vs)
{
// One term
if (v.type() == typeid(Term))
terms.push_back(std::any_cast<Term>(v));
// Multiple terms
else if (v.type() == typeid(std::vector<Term>))
{
const auto & more = std::any_cast<std::vector<Term>>(v);
terms.insert(terms.end(), more.begin(), more.end());
}
}
return terms;
};
// Leaf types can just use 'vs.token()'
parser["Key"] = [](const SemanticValues & vs) { return vs.token(); };
parser["Value"] = [](const SemanticValues & vs) { return vs.token(); };
parser["Pair"] = [](const SemanticValues & vs)
{
mooseAssert(vs.size() == 2, "Metadata entry should have exactly one key and one value");
return std::make_pair(std::string(std::any_cast<std::string_view>(vs[0])),
std::string(std::any_cast<std::string_view>(vs[1])));
};
parser["Metadata"] = [](const SemanticValues & vs)
{
Metadata meta;
for (const auto & v : vs)
{
if (v.type() == typeid(std::pair<std::string, std::string>))
{
const auto & [k, val] = std::any_cast<std::pair<std::string, std::string>>(v);
meta[k] = val;
}
else if (v.type() == typeid(std::vector<std::pair<std::string, std::string>>))
{
const auto & pairs = std::any_cast<std::vector<std::pair<std::string, std::string>>>(v);
for (const auto & [k, val] : pairs)
meta[k] = val;
}
}
return meta;
};
parser["Reaction"] = [](const SemanticValues & vs)
{
Reaction r;
if (vs.size() != 2 && vs.size() != 3)
mooseError("Reaction failed to read as 'reactors -> products [metadata]");
r.reactants = std::any_cast<TermList>(vs[0]);
r.products = std::any_cast<TermList>(vs[1]);
if (vs.size() == 3)
r.metadata = std::any_cast<Metadata>(vs[2]);
return r;
};
parser.enable_packrat_parsing();
std::vector<Reaction> reactions;
if (!parser.parse(reaction_network_string, reactions))
mooseError("Failed to parse reaction.");
// Output the reaction network
if (output_to_cout)
{
for (size_t i = 0; i < reactions.size(); ++i)
{
const auto & r = reactions[i];
Moose::out << "Reaction " << i + 1 << ":\n";
Moose::out << Moose::stringify(r) << std::endl;
}
}
return reactions;
}
std::vector<VariableName>
Reaction::getSpecies() const
{
std::vector<VariableName> species;
for (const auto & term : reactants)
species.push_back(term.species + term.state.value_or("") + term.charge.value_or(""));
for (const auto & term : products)
species.push_back(term.species + term.state.value_or("") + term.charge.value_or(""));
return species;
}
std::vector<VariableName>
Reaction::getUniqueSpecies() const
{
std::vector<VariableName> species;
std::set<Term> terms_gathered;
for (const auto & term : reactants)
if (terms_gathered.find(term) == terms_gathered.end())
{
species.push_back(term.species + term.state.value_or("") + term.charge.value_or(""));
terms_gathered.insert(term);
}
for (const auto & term : products)
if (terms_gathered.find(term) == terms_gathered.end())
{
species.push_back(term.species + term.state.value_or("") + term.charge.value_or(""));
terms_gathered.insert(term);
}
return species;
}
std::vector<Real>
Reaction::getStoichiometricCoefficients() const
{
std::vector<Real> coeffs;
for (const auto & term : reactants)
coeffs.push_back(term.coefficient);
for (const auto & term : products)
coeffs.push_back(term.coefficient);
return coeffs;
}
std::map<VariableName, Real>
Reaction::getUniqueStoichiometricCoefficients() const
{
std::map<VariableName, Real> terms_gathered;
for (const auto & term : reactants)
{
const auto term_str = term.species + term.state.value_or("") + term.charge.value_or("");
terms_gathered[term_str] -= term.coefficient;
}
for (const auto & term : products)
{
const auto term_str = term.species + term.state.value_or("") + term.charge.value_or("");
terms_gathered[term_str] += term.coefficient;
}
return terms_gathered;
}
std::vector<VariableName>
Reaction::getReactantSpecies() const
{
std::vector<VariableName> species;
for (const auto & term : reactants)
species.push_back(term.species + term.state.value_or("") + term.charge.value_or(""));
return species;
}
std::vector<VariableName>
Reaction::getUniqueReactantSpecies() const
{
std::vector<VariableName> species;
std::set<Term> terms_gathered;
for (const auto & term : reactants)
if (terms_gathered.find(term) == terms_gathered.end())
{
species.push_back(term.species + term.state.value_or("") + term.charge.value_or(""));
terms_gathered.insert(term);
}
return species;
}
std::vector<std::string>
Reaction::getProductSpecies() const
{
std::vector<std::string> species;
for (const auto & term : products)
species.push_back(term.species + term.state.value_or("") + term.charge.value_or(""));
return species;
}
std::vector<std::string>
Reaction::getUniqueProductSpecies() const
{
std::vector<std::string> species;
std::set<Term> terms_gathered;
for (const auto & term : products)
if (terms_gathered.find(term) == terms_gathered.end())
{
species.push_back(term.species + term.state.value_or("") + term.charge.value_or(""));
terms_gathered.insert(term);
}
return species;
}
template <>
bool
Reaction::hasMetaData<Real>(const std::string & key) const
{
return metadata.count(key) && MooseUtils::isFloat(libmesh_map_find(metadata, key));
}
template <typename T>
bool
Reaction::hasMetaData(const std::string & /*key*/) const
{
mooseError("Has metadata not implemented for type " + MooseUtils::prettyCppType<T>());
return false;
}
bool
Reaction::hasMetaData(const std::string & key) const
{
return metadata.count(key);
}
const std::string &
Reaction::getMetaData(const std::string & key) const
{
// Get a better error
if (hasMetaData(key))
return libmesh_map_find(metadata, key);
else
mooseError("MetaData item '" + key + "' was not found in reaction:\n" +
Moose::stringify(*this));
}
// namespace ReactionNetworkUtils
}
namespace Moose
{
std::string
stringify(const ReactionNetworkUtils::Reaction & r)
{
std::string str_out = " Reactants:\n";
for (const auto & t : r.reactants)
{
str_out += " " + std::to_string(t.coefficient) + " " + t.species;
if (t.charge)
str_out += t.charge.value();
if (t.state)
str_out += " (" + t.state.value() + ")";
str_out += "\n";
}
str_out += " Products:\n";
for (const auto & t : r.products)
{
str_out += " " + std::to_string(t.coefficient) + " " + t.species;
if (t.charge)
str_out += t.charge.value();
if (t.state)
str_out += " (" + t.state.value() + ")";
str_out += "\n";
}
if (!r.metadata.empty())
{
str_out += " Metadata:\n";
for (const auto & [k, v] : r.metadata)
str_out += " " + k + " = " + v + "\n";
}
return str_out;
}
}
(modules/chemical_reactions/src/utils/ReactionNetworkUtils.C)
// This file is part of the MOOSE framework
// https://mooseframework.inl.gov
//
// All rights reserved, see COPYRIGHT for full restrictions
// https://github.com/idaholab/moose/blob/master/COPYRIGHT
//
// Licensed under LGPL 2.1, please see LICENSE for details
// https://www.gnu.org/licenses/lgpl-2.1.html
#include "ReactionNetworkUtils.h"
#include "MooseTypes.h"
#include "MooseUtils.h"
#include "libmesh/utility.h"
#include "peglib.h"
// not needed
#include "Conversion.h"
namespace ReactionNetworkUtils
{
std::vector<Reaction>
parseReactionNetwork(const std::string & reaction_network_string, bool output_to_cout)
{
using namespace peg;
parser parser(R"(
Reactions <- (_ Reaction _ Newline?)* _
Reaction <- TermList _ '->' _ TermList _ Metadata?
TermList <- Term (_ Term)*
Term <- Coefficient? Species
Coefficient <- Float
Float <- [+-]? (_)? ([0-9]+)? ('.' [0-9]+)?
Species <- BaseSpecies State? Charge?
BaseSpecies <- [a-zA-Z][a-zA-Z0-9_]*
State <- '(' [a-z]{1,3} ')'
Charge <- ('^'? [0-9]* [+-])
Metadata <- '[' _ Pair (_ ',' _ Pair)* _ ']'
Pair <- Key _ '=' _ Value
Key <- [A-Za-z_][A-Za-z0-9_]*
Value <- [A-Za-z0-9_.-]+
Newline <- [\r\n]+
~_ <- [ \t]*
)");
mooseAssert(parser, "Failed to build parser");
// Rules for how to parse each type
parser["Reactions"] = [](const SemanticValues & vs)
{
std::vector<Reaction> reactions;
for (const auto & v : vs)
{
if (v.type() == typeid(Reaction))
{
reactions.push_back(std::any_cast<Reaction>(v));
}
}
return reactions;
};
parser["Float"] = [](const SemanticValues & vs)
{
// We need to remove whitespace in case the reaction was A - B -> C
std::string vs_str = std::string(std::any_cast<std::string_view>(vs.token()));
vs_str.erase(std::remove(vs_str.begin(), vs_str.end(), ' '), vs_str.end());
if (vs_str == "-")
return -1.;
else if (vs_str == "+" || vs_str == "")
return 1.;
else
return std::stod(std::string(vs_str));
};
parser["Coefficient"] = [](const SemanticValues & vs) { return std::any_cast<double>(vs[0]); };
parser["BaseSpecies"] = [](const SemanticValues & vs) { return vs.token(); };
parser["State"] = [](const SemanticValues & vs) { return vs.token(); };
parser["Charge"] = [](const SemanticValues & vs) { return vs.token(); };
parser["Species"] = [](const SemanticValues & vs)
{
std::string base = std::string(std::any_cast<std::string_view>(vs[0]));
std::optional<std::string> state;
std::optional<std::string> charge;
// Parse the state and charge
if (vs.size() >= 2)
{
for (size_t i = 1; i < vs.size(); ++i)
if (vs[i].type() == typeid(std::string_view))
{
std::string val = std::string(std::any_cast<std::string_view>(vs[i]));
if (val.front() == '(' && val.back() == ')')
state = val;
else
charge = val;
}
}
return Term{1.0, base, state, charge}; // default coefficient = 1.0
};
parser["Term"] = [](const SemanticValues & vs)
{
mooseAssert(vs.size() <= 2, "Should only have a coefficient and a species");
if (vs.size() == 2)
{
double coeff = std::any_cast<double>(vs[0]);
Term term = std::any_cast<Term>(vs[1]);
term.coefficient = coeff;
return term;
}
else
return std::any_cast<Term>(vs[0]);
};
parser["TermList"] = [](const SemanticValues & vs)
{
TermList terms;
for (const auto & v : vs)
{
// One term
if (v.type() == typeid(Term))
terms.push_back(std::any_cast<Term>(v));
// Multiple terms
else if (v.type() == typeid(std::vector<Term>))
{
const auto & more = std::any_cast<std::vector<Term>>(v);
terms.insert(terms.end(), more.begin(), more.end());
}
}
return terms;
};
// Leaf types can just use 'vs.token()'
parser["Key"] = [](const SemanticValues & vs) { return vs.token(); };
parser["Value"] = [](const SemanticValues & vs) { return vs.token(); };
parser["Pair"] = [](const SemanticValues & vs)
{
mooseAssert(vs.size() == 2, "Metadata entry should have exactly one key and one value");
return std::make_pair(std::string(std::any_cast<std::string_view>(vs[0])),
std::string(std::any_cast<std::string_view>(vs[1])));
};
parser["Metadata"] = [](const SemanticValues & vs)
{
Metadata meta;
for (const auto & v : vs)
{
if (v.type() == typeid(std::pair<std::string, std::string>))
{
const auto & [k, val] = std::any_cast<std::pair<std::string, std::string>>(v);
meta[k] = val;
}
else if (v.type() == typeid(std::vector<std::pair<std::string, std::string>>))
{
const auto & pairs = std::any_cast<std::vector<std::pair<std::string, std::string>>>(v);
for (const auto & [k, val] : pairs)
meta[k] = val;
}
}
return meta;
};
parser["Reaction"] = [](const SemanticValues & vs)
{
Reaction r;
if (vs.size() != 2 && vs.size() != 3)
mooseError("Reaction failed to read as 'reactors -> products [metadata]");
r.reactants = std::any_cast<TermList>(vs[0]);
r.products = std::any_cast<TermList>(vs[1]);
if (vs.size() == 3)
r.metadata = std::any_cast<Metadata>(vs[2]);
return r;
};
parser.enable_packrat_parsing();
std::vector<Reaction> reactions;
if (!parser.parse(reaction_network_string, reactions))
mooseError("Failed to parse reaction.");
// Output the reaction network
if (output_to_cout)
{
for (size_t i = 0; i < reactions.size(); ++i)
{
const auto & r = reactions[i];
Moose::out << "Reaction " << i + 1 << ":\n";
Moose::out << Moose::stringify(r) << std::endl;
}
}
return reactions;
}
std::vector<VariableName>
Reaction::getSpecies() const
{
std::vector<VariableName> species;
for (const auto & term : reactants)
species.push_back(term.species + term.state.value_or("") + term.charge.value_or(""));
for (const auto & term : products)
species.push_back(term.species + term.state.value_or("") + term.charge.value_or(""));
return species;
}
std::vector<VariableName>
Reaction::getUniqueSpecies() const
{
std::vector<VariableName> species;
std::set<Term> terms_gathered;
for (const auto & term : reactants)
if (terms_gathered.find(term) == terms_gathered.end())
{
species.push_back(term.species + term.state.value_or("") + term.charge.value_or(""));
terms_gathered.insert(term);
}
for (const auto & term : products)
if (terms_gathered.find(term) == terms_gathered.end())
{
species.push_back(term.species + term.state.value_or("") + term.charge.value_or(""));
terms_gathered.insert(term);
}
return species;
}
std::vector<Real>
Reaction::getStoichiometricCoefficients() const
{
std::vector<Real> coeffs;
for (const auto & term : reactants)
coeffs.push_back(term.coefficient);
for (const auto & term : products)
coeffs.push_back(term.coefficient);
return coeffs;
}
std::map<VariableName, Real>
Reaction::getUniqueStoichiometricCoefficients() const
{
std::map<VariableName, Real> terms_gathered;
for (const auto & term : reactants)
{
const auto term_str = term.species + term.state.value_or("") + term.charge.value_or("");
terms_gathered[term_str] -= term.coefficient;
}
for (const auto & term : products)
{
const auto term_str = term.species + term.state.value_or("") + term.charge.value_or("");
terms_gathered[term_str] += term.coefficient;
}
return terms_gathered;
}
std::vector<VariableName>
Reaction::getReactantSpecies() const
{
std::vector<VariableName> species;
for (const auto & term : reactants)
species.push_back(term.species + term.state.value_or("") + term.charge.value_or(""));
return species;
}
std::vector<VariableName>
Reaction::getUniqueReactantSpecies() const
{
std::vector<VariableName> species;
std::set<Term> terms_gathered;
for (const auto & term : reactants)
if (terms_gathered.find(term) == terms_gathered.end())
{
species.push_back(term.species + term.state.value_or("") + term.charge.value_or(""));
terms_gathered.insert(term);
}
return species;
}
std::vector<std::string>
Reaction::getProductSpecies() const
{
std::vector<std::string> species;
for (const auto & term : products)
species.push_back(term.species + term.state.value_or("") + term.charge.value_or(""));
return species;
}
std::vector<std::string>
Reaction::getUniqueProductSpecies() const
{
std::vector<std::string> species;
std::set<Term> terms_gathered;
for (const auto & term : products)
if (terms_gathered.find(term) == terms_gathered.end())
{
species.push_back(term.species + term.state.value_or("") + term.charge.value_or(""));
terms_gathered.insert(term);
}
return species;
}
template <>
bool
Reaction::hasMetaData<Real>(const std::string & key) const
{
return metadata.count(key) && MooseUtils::isFloat(libmesh_map_find(metadata, key));
}
template <typename T>
bool
Reaction::hasMetaData(const std::string & /*key*/) const
{
mooseError("Has metadata not implemented for type " + MooseUtils::prettyCppType<T>());
return false;
}
bool
Reaction::hasMetaData(const std::string & key) const
{
return metadata.count(key);
}
const std::string &
Reaction::getMetaData(const std::string & key) const
{
// Get a better error
if (hasMetaData(key))
return libmesh_map_find(metadata, key);
else
mooseError("MetaData item '" + key + "' was not found in reaction:\n" +
Moose::stringify(*this));
}
// namespace ReactionNetworkUtils
}
namespace Moose
{
std::string
stringify(const ReactionNetworkUtils::Reaction & r)
{
std::string str_out = " Reactants:\n";
for (const auto & t : r.reactants)
{
str_out += " " + std::to_string(t.coefficient) + " " + t.species;
if (t.charge)
str_out += t.charge.value();
if (t.state)
str_out += " (" + t.state.value() + ")";
str_out += "\n";
}
str_out += " Products:\n";
for (const auto & t : r.products)
{
str_out += " " + std::to_string(t.coefficient) + " " + t.species;
if (t.charge)
str_out += t.charge.value();
if (t.state)
str_out += " (" + t.state.value() + ")";
str_out += "\n";
}
if (!r.metadata.empty())
{
str_out += " Metadata:\n";
for (const auto & [k, v] : r.metadata)
str_out += " " + k + " = " + v + "\n";
}
return str_out;
}
}