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