Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "ReactionNetworkUtils.h"
      11             : #include "MooseTypes.h"
      12             : #include "MooseUtils.h"
      13             : #include "libmesh/utility.h"
      14             : #include "peglib.h"
      15             : 
      16             : // not needed
      17             : #include "Conversion.h"
      18             : 
      19             : namespace ReactionNetworkUtils
      20             : {
      21             : 
      22             : std::vector<Reaction>
      23          51 : parseReactionNetwork(const std::string & reaction_network_string, bool output_to_cout)
      24             : {
      25             :   using namespace peg;
      26             :   parser parser(R"(
      27             :     Reactions     <- (_ Reaction _ Newline?)* _
      28             :     Reaction      <- TermList _ '->' _ TermList _ Metadata?
      29             :     TermList      <- Term (_ Term)*
      30             :     Term          <- Coefficient? Species
      31             :     Coefficient   <- Float
      32             :     Float         <- [+-]? (_)? ([0-9]+)? ('.' [0-9]+)?
      33             :     Species       <- BaseSpecies State? Charge?
      34             :     BaseSpecies   <- [a-zA-Z][a-zA-Z0-9_]*
      35             :     State         <- '(' [a-z]{1,3} ')'
      36             :     Charge        <- ('^'? [0-9]* [+-])
      37             :     Metadata      <- '[' _ Pair (_ ',' _ Pair)* _ ']'
      38             :     Pair          <- Key _ '=' _ Value
      39             :     Key           <- [A-Za-z_][A-Za-z0-9_]*
      40             :     Value         <- [A-Za-z0-9_.-]+
      41             :     Newline       <- [\r\n]+
      42             :     ~_            <- [ \t]*
      43          51 :   )");
      44             :   mooseAssert(parser, "Failed to build parser");
      45             : 
      46             :   // Rules for how to parse each type
      47         153 :   parser["Reactions"] = [](const SemanticValues & vs)
      48             :   {
      49             :     std::vector<Reaction> reactions;
      50         162 :     for (const auto & v : vs)
      51             :     {
      52         111 :       if (v.type() == typeid(Reaction))
      53             :       {
      54         162 :         reactions.push_back(std::any_cast<Reaction>(v));
      55             :       }
      56             :     }
      57          51 :     return reactions;
      58           0 :   };
      59             : 
      60         514 :   parser["Float"] = [](const SemanticValues & vs)
      61             :   {
      62             :     // We need to remove whitespace in case the reaction was A - B -> C
      63         824 :     std::string vs_str = std::string(std::any_cast<std::string_view>(vs.token()));
      64         412 :     vs_str.erase(std::remove(vs_str.begin(), vs_str.end(), ' '), vs_str.end());
      65         412 :     if (vs_str == "-")
      66             :       return -1.;
      67         299 :     else if (vs_str == "+" || vs_str == "")
      68             :       return 1.;
      69             :     else
      70         100 :       return std::stod(std::string(vs_str));
      71             :   };
      72             : 
      73         514 :   parser["Coefficient"] = [](const SemanticValues & vs) { return std::any_cast<double>(vs[0]); };
      74             : 
      75         102 :   parser["BaseSpecies"] = [](const SemanticValues & vs) { return vs.token(); };
      76             : 
      77         102 :   parser["State"] = [](const SemanticValues & vs) { return vs.token(); };
      78             : 
      79         102 :   parser["Charge"] = [](const SemanticValues & vs) { return vs.token(); };
      80             : 
      81         312 :   parser["Species"] = [](const SemanticValues & vs)
      82             :   {
      83         210 :     std::string base = std::string(std::any_cast<std::string_view>(vs[0]));
      84         210 :     std::optional<std::string> state;
      85         210 :     std::optional<std::string> charge;
      86             : 
      87             :     // Parse the state and charge
      88         210 :     if (vs.size() >= 2)
      89             :     {
      90         158 :       for (size_t i = 1; i < vs.size(); ++i)
      91          81 :         if (vs[i].type() == typeid(std::string_view))
      92             :         {
      93         162 :           std::string val = std::string(std::any_cast<std::string_view>(vs[i]));
      94          81 :           if (val.front() == '(' && val.back() == ')')
      95           7 :             state = val;
      96             :           else
      97          74 :             charge = val;
      98             :         }
      99             :     }
     100             : 
     101         420 :     return Term{1.0, base, state, charge}; // default coefficient = 1.0
     102             :   };
     103             : 
     104         312 :   parser["Term"] = [](const SemanticValues & vs)
     105             :   {
     106             :     mooseAssert(vs.size() <= 2, "Should only have a coefficient and a species");
     107         210 :     if (vs.size() == 2)
     108             :     {
     109         210 :       double coeff = std::any_cast<double>(vs[0]);
     110         210 :       Term term = std::any_cast<Term>(vs[1]);
     111         210 :       term.coefficient = coeff;
     112         210 :       return term;
     113         210 :     }
     114             :     else
     115           0 :       return std::any_cast<Term>(vs[0]);
     116             :   };
     117             : 
     118         264 :   parser["TermList"] = [](const SemanticValues & vs)
     119             :   {
     120             :     TermList terms;
     121         372 :     for (const auto & v : vs)
     122             :     {
     123             :       // One term
     124         210 :       if (v.type() == typeid(Term))
     125         420 :         terms.push_back(std::any_cast<Term>(v));
     126             :       // Multiple terms
     127           0 :       else if (v.type() == typeid(std::vector<Term>))
     128             :       {
     129           0 :         const auto & more = std::any_cast<std::vector<Term>>(v);
     130           0 :         terms.insert(terms.end(), more.begin(), more.end());
     131           0 :       }
     132             :     }
     133         162 :     return terms;
     134           0 :   };
     135             : 
     136             :   // Leaf types can just use 'vs.token()'
     137         102 :   parser["Key"] = [](const SemanticValues & vs) { return vs.token(); };
     138             : 
     139         102 :   parser["Value"] = [](const SemanticValues & vs) { return vs.token(); };
     140             : 
     141         175 :   parser["Pair"] = [](const SemanticValues & vs)
     142             :   {
     143             :     mooseAssert(vs.size() == 2, "Metadata entry should have exactly one key and one value");
     144         146 :     return std::make_pair(std::string(std::any_cast<std::string_view>(vs[0])),
     145         146 :                           std::string(std::any_cast<std::string_view>(vs[1])));
     146             :   };
     147             : 
     148         170 :   parser["Metadata"] = [](const SemanticValues & vs)
     149             :   {
     150             :     Metadata meta;
     151         141 :     for (const auto & v : vs)
     152             :     {
     153          73 :       if (v.type() == typeid(std::pair<std::string, std::string>))
     154             :       {
     155          73 :         const auto & [k, val] = std::any_cast<std::pair<std::string, std::string>>(v);
     156          73 :         meta[k] = val;
     157             :       }
     158           0 :       else if (v.type() == typeid(std::vector<std::pair<std::string, std::string>>))
     159             :       {
     160           0 :         const auto & pairs = std::any_cast<std::vector<std::pair<std::string, std::string>>>(v);
     161           0 :         for (const auto & [k, val] : pairs)
     162           0 :           meta[k] = val;
     163           0 :       }
     164             :     }
     165          68 :     return meta;
     166             :   };
     167             : 
     168         132 :   parser["Reaction"] = [](const SemanticValues & vs)
     169             :   {
     170             :     Reaction r;
     171          81 :     if (vs.size() != 2 && vs.size() != 3)
     172           0 :       mooseError("Reaction failed to read as 'reactors -> products [metadata]");
     173          81 :     r.reactants = std::any_cast<TermList>(vs[0]);
     174          81 :     r.products = std::any_cast<TermList>(vs[1]);
     175          81 :     if (vs.size() == 3)
     176         136 :       r.metadata = std::any_cast<Metadata>(vs[2]);
     177          81 :     return r;
     178           0 :   };
     179             : 
     180          51 :   parser.enable_packrat_parsing();
     181             : 
     182             :   std::vector<Reaction> reactions;
     183          51 :   if (!parser.parse(reaction_network_string, reactions))
     184           0 :     mooseError("Failed to parse reaction.");
     185             : 
     186             :   // Output the reaction network
     187          51 :   if (output_to_cout)
     188             :   {
     189          18 :     for (size_t i = 0; i < reactions.size(); ++i)
     190             :     {
     191             :       const auto & r = reactions[i];
     192           9 :       Moose::out << "Reaction " << i + 1 << ":\n";
     193          18 :       Moose::out << Moose::stringify(r) << std::endl;
     194             :     }
     195             :   }
     196             : 
     197          51 :   return reactions;
     198          51 : }
     199             : 
     200             : std::vector<VariableName>
     201          57 : Reaction::getSpecies() const
     202             : {
     203             :   std::vector<VariableName> species;
     204         163 :   for (const auto & term : reactants)
     205         212 :     species.push_back(term.species + term.state.value_or("") + term.charge.value_or(""));
     206         148 :   for (const auto & term : products)
     207         182 :     species.push_back(term.species + term.state.value_or("") + term.charge.value_or(""));
     208          57 :   return species;
     209           0 : }
     210             : 
     211             : std::vector<VariableName>
     212           6 : Reaction::getUniqueSpecies() const
     213             : {
     214             :   std::vector<VariableName> species;
     215             :   std::set<Term> terms_gathered;
     216          20 :   for (const auto & term : reactants)
     217          14 :     if (terms_gathered.find(term) == terms_gathered.end())
     218             :     {
     219          20 :       species.push_back(term.species + term.state.value_or("") + term.charge.value_or(""));
     220          10 :       terms_gathered.insert(term);
     221             :     }
     222          20 :   for (const auto & term : products)
     223          14 :     if (terms_gathered.find(term) == terms_gathered.end())
     224             :     {
     225          20 :       species.push_back(term.species + term.state.value_or("") + term.charge.value_or(""));
     226          10 :       terms_gathered.insert(term);
     227             :     }
     228           6 :   return species;
     229           0 : }
     230             : 
     231             : std::vector<Real>
     232         119 : Reaction::getStoichiometricCoefficients() const
     233             : {
     234             :   std::vector<Real> coeffs;
     235         312 :   for (const auto & term : reactants)
     236         193 :     coeffs.push_back(term.coefficient);
     237         272 :   for (const auto & term : products)
     238         153 :     coeffs.push_back(term.coefficient);
     239         119 :   return coeffs;
     240           0 : }
     241             : 
     242             : std::map<VariableName, Real>
     243           6 : Reaction::getUniqueStoichiometricCoefficients() const
     244             : {
     245             :   std::map<VariableName, Real> terms_gathered;
     246          20 :   for (const auto & term : reactants)
     247             :   {
     248          28 :     const auto term_str = term.species + term.state.value_or("") + term.charge.value_or("");
     249          14 :     terms_gathered[term_str] -= term.coefficient;
     250             :   }
     251          20 :   for (const auto & term : products)
     252             :   {
     253          28 :     const auto term_str = term.species + term.state.value_or("") + term.charge.value_or("");
     254          14 :     terms_gathered[term_str] += term.coefficient;
     255             :   }
     256           6 :   return terms_gathered;
     257             : }
     258             : 
     259             : std::vector<VariableName>
     260          92 : Reaction::getReactantSpecies() const
     261             : {
     262             :   std::vector<VariableName> species;
     263         235 :   for (const auto & term : reactants)
     264         286 :     species.push_back(term.species + term.state.value_or("") + term.charge.value_or(""));
     265          92 :   return species;
     266           0 : }
     267             : 
     268             : std::vector<VariableName>
     269           2 : Reaction::getUniqueReactantSpecies() const
     270             : {
     271             :   std::vector<VariableName> species;
     272             :   std::set<Term> terms_gathered;
     273           8 :   for (const auto & term : reactants)
     274           6 :     if (terms_gathered.find(term) == terms_gathered.end())
     275             :     {
     276           4 :       species.push_back(term.species + term.state.value_or("") + term.charge.value_or(""));
     277           2 :       terms_gathered.insert(term);
     278             :     }
     279           2 :   return species;
     280           0 : }
     281             : 
     282             : std::vector<std::string>
     283          88 : Reaction::getProductSpecies() const
     284             : {
     285             :   std::vector<std::string> species;
     286         192 :   for (const auto & term : products)
     287         208 :     species.push_back(term.species + term.state.value_or("") + term.charge.value_or(""));
     288          88 :   return species;
     289           0 : }
     290             : 
     291             : std::vector<std::string>
     292           3 : Reaction::getUniqueProductSpecies() const
     293             : {
     294             :   std::vector<std::string> species;
     295             :   std::set<Term> terms_gathered;
     296          10 :   for (const auto & term : products)
     297           7 :     if (terms_gathered.find(term) == terms_gathered.end())
     298             :     {
     299          10 :       species.push_back(term.species + term.state.value_or("") + term.charge.value_or(""));
     300           5 :       terms_gathered.insert(term);
     301             :     }
     302           3 :   return species;
     303           0 : }
     304             : 
     305             : template <>
     306             : bool
     307          78 : Reaction::hasMetaData<Real>(const std::string & key) const
     308             : {
     309          69 :   return metadata.count(key) && MooseUtils::isFloat(libmesh_map_find(metadata, key));
     310             : }
     311             : 
     312             : template <typename T>
     313             : bool
     314             : Reaction::hasMetaData(const std::string & /*key*/) const
     315             : {
     316             :   mooseError("Has metadata not implemented for type " + MooseUtils::prettyCppType<T>());
     317             :   return false;
     318             : }
     319             : 
     320             : bool
     321          94 : Reaction::hasMetaData(const std::string & key) const
     322             : {
     323          94 :   return metadata.count(key);
     324             : }
     325             : 
     326             : const std::string &
     327          70 : Reaction::getMetaData(const std::string & key) const
     328             : {
     329             :   // Get a better error
     330          70 :   if (hasMetaData(key))
     331          69 :     return libmesh_map_find(metadata, key);
     332             :   else
     333           3 :     mooseError("MetaData item '" + key + "' was not found in reaction:\n" +
     334           1 :                Moose::stringify(*this));
     335             : }
     336             : 
     337             : // namespace ReactionNetworkUtils
     338             : }
     339             : 
     340             : namespace Moose
     341             : {
     342             : std::string
     343          10 : stringify(const ReactionNetworkUtils::Reaction & r)
     344             : {
     345          10 :   std::string str_out = "  Reactants:\n";
     346          20 :   for (const auto & t : r.reactants)
     347             :   {
     348          30 :     str_out += "    " + std::to_string(t.coefficient) + " " + t.species;
     349          10 :     if (t.charge)
     350             :       str_out += t.charge.value();
     351          10 :     if (t.state)
     352           0 :       str_out += " (" + t.state.value() + ")";
     353             :     str_out += "\n";
     354             :   }
     355             :   str_out += "  Products:\n";
     356          20 :   for (const auto & t : r.products)
     357             :   {
     358          30 :     str_out += "    " + std::to_string(t.coefficient) + " " + t.species;
     359          10 :     if (t.charge)
     360             :       str_out += t.charge.value();
     361          10 :     if (t.state)
     362           0 :       str_out += " (" + t.state.value() + ")";
     363             :     str_out += "\n";
     364             :   }
     365          10 :   if (!r.metadata.empty())
     366             :   {
     367             :     str_out += "  Metadata:\n";
     368          18 :     for (const auto & [k, v] : r.metadata)
     369          27 :       str_out += "    " + k + " = " + v + "\n";
     370             :   }
     371          10 :   return str_out;
     372             : }
     373             : }