www.mooseframework.org
AddCoupledEqSpeciesAction.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
8 #include "Parser.h"
9 #include "FEProblem.h"
10 #include "Factory.h"
11 #include "MooseError.h"
12 
13 #include "libmesh/string_to_enum.h"
14 
15 // Regular expression includes
16 #include "pcrecpp.h"
17 
18 template <>
19 InputParameters
21 {
22  InputParameters params = validParams<Action>();
23  params.addRequiredParam<std::vector<NonlinearVariableName>>(
24  "primary_species", "The list of primary variables to add");
25  params.addParam<std::vector<AuxVariableName>>(
26  "secondary_species", "The list of aqueous equilibrium species to be output as aux variables");
27  params.addParam<std::string>("reactions", "The list of aqueous equilibrium reactions");
28  params.addParam<std::vector<VariableName>>("pressure", "Pressure variable");
29  RealVectorValue g(0, 0, 0);
30  params.addParam<RealVectorValue>("gravity", g, "Gravity vector (default is (0, 0, 0))");
31  params.addClassDescription("Adds coupled equilibrium Kernels and AuxKernels for primary species");
32  return params;
33 }
34 
36  : Action(params),
37  _primary_species(getParam<std::vector<NonlinearVariableName>>("primary_species")),
38  _secondary_species(getParam<std::vector<AuxVariableName>>("secondary_species")),
39  _weights(_primary_species.size()),
40  _primary_participation(_primary_species.size()),
41  _sto_u(_primary_species.size()),
42  _sto_v(_primary_species.size()),
43  _coupled_v(_primary_species.size()),
44  _input_reactions(getParam<std::string>("reactions")),
45  _pressure_var(getParam<std::vector<VariableName>>("pressure")),
46  _gravity(getParam<RealVectorValue>("gravity"))
47 {
48  // Parse the aqueous equilibrium reactions
49  pcrecpp::RE re_reaction(
50  "(.+?)" // single reaction (any character until the equilibrium coefficient appears)
51  "\\s" // word boundary
52  "(" // start capture
53  "-?" // optional minus sign
54  "\\d+(?:\\.\\d*)?" // digits followed by optional decimal and more digits
55  ")" // end capture
56  "\\b" // word boundary
57  "(?:,?\\s*|$)" // optional comma or end of string
58  ,
59  pcrecpp::RE_Options().set_extended(true));
60 
61  pcrecpp::RE re_terms("(\\S+)");
62  pcrecpp::RE re_coeff_and_species("(?: \\(? (.*?) \\)? )" // match the leading coefficent
63  "([A-Za-z].*)" // match the species
64  ,
65  pcrecpp::RE_Options().set_extended(true));
66 
67  pcrecpp::StringPiece input(_input_reactions);
68 
69  pcrecpp::StringPiece single_reaction, term;
70  std::string single_reaction_str;
71  Real eq_coeff;
72 
73  // Parse reaction network to extract each individual reaction
74  while (re_reaction.FindAndConsume(&input, &single_reaction_str, &eq_coeff))
75  {
76  _reactions.push_back(single_reaction_str);
77  _eq_const.push_back(eq_coeff);
78  }
79 
80  _num_reactions = _reactions.size();
81 
82  if (_num_reactions == 0)
83  mooseError("No equilibrium reaction provided!");
84 
85  if (_pars.isParamValid("secondary_species"))
86  for (unsigned int k = 0; k < _secondary_species.size(); ++k)
88 
89  // Start parsing each reaction
90  for (unsigned int i = 0; i < _num_reactions; ++i)
91  {
92  single_reaction = _reactions[i];
93 
94  // Capture all of the terms
95  std::string species, coeff_str;
96  Real coeff;
97  int sign = 1;
98  bool secondary = false;
99 
100  std::vector<Real> local_stos;
101  std::vector<VariableName> local_species_list;
102 
103  // Find every single term in this reaction (species and operators)
104  while (re_terms.FindAndConsume(&single_reaction, &term))
105  {
106  // Separating the stoichiometric coefficients from species
107  if (re_coeff_and_species.PartialMatch(term, &coeff_str, &species))
108  {
109  if (coeff_str.length())
110  coeff = std::stod(coeff_str);
111  else
112  coeff = 1.0;
113 
114  coeff *= sign;
115 
116  if (secondary)
117  _eq_species.push_back(species);
118  else
119  {
120  local_stos.push_back(coeff);
121  local_species_list.push_back(species);
122  }
123  }
124  // Finding the operators and assign value of -1.0 to "-" sign
125  else if (term == "+" || term == "=" || term == "-")
126  {
127  if (term == "-")
128  {
129  sign = -1;
130  term = "+";
131  }
132  else
133  sign = 1;
134 
135  if (term == "=")
136  secondary = true;
137  }
138  else
139  mooseError("Error parsing term: ", term.as_string());
140  }
141 
142  _stos.push_back(local_stos);
143  _primary_species_involved.push_back(local_species_list);
144  }
145 
146  // Start picking out primary species and coupled primary species and assigning
147  // corresponding stoichiomentric coefficients
148  for (unsigned int i = 0; i < _primary_species.size(); ++i)
149  {
150  _sto_u[i].resize(_num_reactions, 0.0);
151  _sto_v[i].resize(_num_reactions);
152  _coupled_v[i].resize(_num_reactions);
153  _weights[i].resize(_num_reactions, 0.0);
154 
155  _primary_participation[i].resize(_num_reactions, false);
156  for (unsigned int j = 0; j < _num_reactions; ++j)
157  {
158  for (unsigned int k = 0; k < _primary_species_involved[j].size(); ++k)
160  _primary_participation[i][j] = true;
161 
162  if (_primary_participation[i][j])
163  {
164  for (unsigned int k = 0; k < _primary_species_involved[j].size(); ++k)
165  {
167  {
168  _sto_u[i][j] = _stos[j][k];
169  _weights[i][j] = _stos[j][k];
170  }
171  else
172  {
173  _sto_v[i][j].push_back(_stos[j][k]);
174  _coupled_v[i][j].push_back(_primary_species_involved[j][k]);
175  }
176  }
177  }
178  }
179  }
180 
181  // Print out details of the equilibrium reactions to the console
182  _console << "Aqueous equilibrium reactions:\n";
183  for (unsigned int i = 0; i < _num_reactions; ++i)
184  _console << " Reaction " << i + 1 << ": " << _reactions[i] << ", log10(Keq) = " << _eq_const[i]
185  << "\n";
186  _console << "\n";
187 }
188 
189 void
191 {
192  if (_current_task == "add_kernel")
193  {
194  // Add Kernels for each primary species
195  for (unsigned int i = 0; i < _primary_species.size(); ++i)
196  {
197  for (unsigned int j = 0; j < _num_reactions; ++j)
198  {
199  if (_primary_participation[i][j])
200  {
201  InputParameters params_sub = _factory.getValidParams("CoupledBEEquilibriumSub");
202  params_sub.set<NonlinearVariableName>("variable") = _primary_species[i];
203  params_sub.set<Real>("weight") = _weights[i][j];
204  params_sub.set<Real>("log_k") = _eq_const[j];
205  params_sub.set<Real>("sto_u") = _sto_u[i][j];
206  params_sub.set<std::vector<Real>>("sto_v") = _sto_v[i][j];
207  params_sub.set<std::vector<VariableName>>("v") = _coupled_v[i][j];
208  _problem->addKernel("CoupledBEEquilibriumSub",
209  _primary_species[i] + "_" + _eq_species[j] + "_sub",
210  params_sub);
211 
212  InputParameters params_cd = _factory.getValidParams("CoupledDiffusionReactionSub");
213  params_cd.set<NonlinearVariableName>("variable") = _primary_species[i];
214  params_cd.set<Real>("weight") = _weights[i][j];
215  params_cd.set<Real>("log_k") = _eq_const[j];
216  params_cd.set<Real>("sto_u") = _sto_u[i][j];
217  params_cd.set<std::vector<Real>>("sto_v") = _sto_v[i][j];
218  params_cd.set<std::vector<VariableName>>("v") = _coupled_v[i][j];
219  _problem->addKernel("CoupledDiffusionReactionSub",
220  _primary_species[i] + "_" + _eq_species[j] + "_cd",
221  params_cd);
222 
223  // If pressure is coupled, add a CoupledConvectionReactionSub Kernel as well
224  if (_pars.isParamValid("pressure"))
225  {
226  InputParameters params_conv = _factory.getValidParams("CoupledConvectionReactionSub");
227  params_conv.set<NonlinearVariableName>("variable") = _primary_species[i];
228  params_conv.set<Real>("weight") = _weights[i][j];
229  params_conv.set<Real>("log_k") = _eq_const[j];
230  params_conv.set<Real>("sto_u") = _sto_u[i][j];
231  params_conv.set<std::vector<Real>>("sto_v") = _sto_v[i][j];
232  params_conv.set<std::vector<VariableName>>("v") = _coupled_v[i][j];
233  params_conv.set<std::vector<VariableName>>("p") = _pressure_var;
234  params_conv.set<RealVectorValue>("gravity") = _gravity;
235  _problem->addKernel("CoupledConvectionReactionSub",
236  _primary_species[i] + "_" + _eq_species[j] + "_conv",
237  params_conv);
238  }
239  }
240  }
241  }
242  }
243 
244  if (_current_task == "add_aux_kernel")
245  {
246  // Add AqueousEquilibriumRxnAux AuxKernels for equilibrium species
247  for (unsigned int j = 0; j < _num_reactions; ++j)
248  {
249  if (_aux_species.find(_eq_species[j]) != _aux_species.end())
250  {
251  InputParameters params_eq = _factory.getValidParams("AqueousEquilibriumRxnAux");
252  params_eq.set<AuxVariableName>("variable") = _eq_species[j];
253  params_eq.set<Real>("log_k") = _eq_const[j];
254  params_eq.set<std::vector<Real>>("sto_v") = _stos[j];
255  params_eq.set<std::vector<VariableName>>("v") = _primary_species_involved[j];
256  _problem->addAuxKernel("AqueousEquilibriumRxnAux", "aux_" + _eq_species[j], params_eq);
257  }
258  }
259  }
260 }
AddCoupledEqSpeciesAction(const InputParameters &params)
std::vector< Real > _eq_const
Equilibrium constants for each reaction.
const std::vector< AuxVariableName > _secondary_species
Secondary species added as AuxVariables.
Real sign(Real x)
Definition: MathUtils.h:24
std::vector< VariableName > _eq_species
Equilibrium species.
std::vector< std::vector< bool > > _primary_participation
Participation of primary species in each reaction.
std::vector< std::string > _reactions
Vector of parsed reactions.
std::vector< std::vector< Real > > _sto_u
Stoichiometric coefficients of primary variables in each reaction.
const RealVectorValue _gravity
Gravity (default is (0, 0, 0))
std::set< std::string > _aux_species
Set of auxillary species.
std::vector< std::vector< Real > > _weights
Weight of each primary species in each reaction.
std::vector< std::vector< std::vector< Real > > > _sto_v
Stoichiometric coefficients of coupled primary variables in each reaction.
std::string _input_reactions
Reaction network read from input file.
const std::vector< NonlinearVariableName > _primary_species
Basis set of primary species.
std::vector< std::vector< Real > > _stos
Stoichiometric coefficients for each primary species in each reaction.
unsigned int _num_reactions
Number of reactions.
std::vector< std::vector< std::vector< VariableName > > > _coupled_v
Coupled primary species for each reaction.
InputParameters validParams< AddCoupledEqSpeciesAction >()
std::vector< std::vector< VariableName > > _primary_species_involved
Primary species involved in the ith equilibrium reaction.
const std::vector< VariableName > _pressure_var
Pressure variable.