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 201 : AddCoupledEqSpeciesAction::validParams()
27 : {
28 201 : InputParameters params = Action::validParams();
29 402 : params.addRequiredParam<std::vector<NonlinearVariableName>>(
30 : "primary_species", "The list of primary variables to add");
31 402 : params.addParam<std::vector<AuxVariableName>>(
32 : "secondary_species", "The list of aqueous equilibrium species to be output as aux variables");
33 402 : params.addParam<std::string>("reactions", "The list of aqueous equilibrium reactions");
34 402 : params.addParam<std::vector<VariableName>>("pressure", {}, "Pressure variable");
35 : RealVectorValue g(0, 0, 0);
36 402 : params.addParam<RealVectorValue>("gravity", g, "Gravity vector (default is (0, 0, 0))");
37 201 : params.addClassDescription("Adds coupled equilibrium Kernels and AuxKernels for primary species");
38 201 : return params;
39 0 : }
40 :
41 201 : AddCoupledEqSpeciesAction::AddCoupledEqSpeciesAction(const InputParameters & params)
42 : : Action(params),
43 201 : _primary_species(getParam<std::vector<NonlinearVariableName>>("primary_species")),
44 603 : _secondary_species(getParam<std::vector<AuxVariableName>>("secondary_species")),
45 201 : _weights(_primary_species.size()),
46 201 : _primary_participation(_primary_species.size()),
47 201 : _sto_u(_primary_species.size()),
48 201 : _sto_v(_primary_species.size()),
49 201 : _coupled_v(_primary_species.size()),
50 603 : _input_reactions(getParam<std::string>("reactions")),
51 402 : _pressure_var(getParam<std::vector<VariableName>>("pressure")),
52 603 : _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 201 : pcrecpp::RE_Options().set_extended(true));
66 :
67 201 : 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 201 : 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 868 : while (re_reaction.FindAndConsume(&input, &single_reaction_str, &eq_coeff))
81 : {
82 667 : _reactions.push_back(single_reaction_str);
83 667 : _eq_const.push_back(eq_coeff);
84 : }
85 :
86 201 : _num_reactions = _reactions.size();
87 :
88 201 : if (_num_reactions == 0)
89 0 : mooseError("No equilibrium reaction provided!");
90 :
91 402 : if (_pars.isParamValid("secondary_species"))
92 868 : for (unsigned int k = 0; k < _secondary_species.size(); ++k)
93 667 : _aux_species.insert(_secondary_species[k]);
94 :
95 : // Start parsing each reaction
96 868 : for (unsigned int i = 0; i < _num_reactions; ++i)
97 : {
98 667 : 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 3851 : while (re_terms.FindAndConsume(&single_reaction, &term))
111 : {
112 : // Separating the stoichiometric coefficients from species
113 3184 : if (re_coeff_and_species.PartialMatch(term, &coeff_str, &species))
114 : {
115 1871 : if (coeff_str.length())
116 130 : coeff = std::stod(coeff_str);
117 : else
118 1741 : coeff = 1.0;
119 :
120 1871 : coeff *= sign;
121 :
122 1871 : if (secondary)
123 1334 : _eq_species.push_back(species);
124 : else
125 : {
126 1204 : local_stos.push_back(coeff);
127 1204 : local_species_list.push_back(species);
128 : }
129 : }
130 : // Finding the operators and assign value of -1.0 to "-" sign
131 1313 : else if (term == "+" || term == "=" || term == "-")
132 : {
133 1313 : if (term == "-")
134 : {
135 : sign = -1;
136 341 : term = "+";
137 : }
138 : else
139 : sign = 1;
140 :
141 1313 : if (term == "=")
142 : secondary = true;
143 : }
144 : else
145 0 : mooseError("Error parsing term: ", term.as_string());
146 : }
147 :
148 667 : _stos.push_back(local_stos);
149 667 : _primary_species_involved.push_back(local_species_list);
150 667 : }
151 :
152 : // Start picking out primary species and coupled primary species and assigning
153 : // corresponding stoichiomentric coefficients
154 636 : for (unsigned int i = 0; i < _primary_species.size(); ++i)
155 : {
156 435 : _sto_u[i].resize(_num_reactions, 0.0);
157 435 : _sto_v[i].resize(_num_reactions);
158 435 : _coupled_v[i].resize(_num_reactions);
159 435 : _weights[i].resize(_num_reactions, 0.0);
160 :
161 435 : _primary_participation[i].resize(_num_reactions, false);
162 2157 : for (unsigned int j = 0; j < _num_reactions; ++j)
163 : {
164 4944 : for (unsigned int k = 0; k < _primary_species_involved[j].size(); ++k)
165 3222 : if (_primary_species_involved[j][k] == _primary_species[i])
166 : _primary_participation[i][j] = true;
167 :
168 1722 : if (_primary_participation[i][j])
169 : {
170 3624 : for (unsigned int k = 0; k < _primary_species_involved[j].size(); ++k)
171 : {
172 2420 : if (_primary_species_involved[j][k] == _primary_species[i])
173 : {
174 1204 : _sto_u[i][j] = _stos[j][k];
175 1204 : _weights[i][j] = _stos[j][k];
176 : }
177 : else
178 : {
179 1216 : _sto_v[i][j].push_back(_stos[j][k]);
180 1216 : _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 201 : _console << "Aqueous equilibrium reactions:\n";
189 868 : for (unsigned int i = 0; i < _num_reactions; ++i)
190 1334 : _console << " Reaction " << i + 1 << ": " << _reactions[i] << ", log10(Keq) = " << _eq_const[i]
191 667 : << "\n";
192 201 : _console << std::endl;
193 201 : }
194 :
195 : void
196 402 : AddCoupledEqSpeciesAction::act()
197 : {
198 402 : if (_current_task == "add_kernel")
199 : {
200 : // Add Kernels for each primary species
201 636 : for (unsigned int i = 0; i < _primary_species.size(); ++i)
202 : {
203 2157 : for (unsigned int j = 0; j < _num_reactions; ++j)
204 : {
205 1722 : if (_primary_participation[i][j])
206 : {
207 2408 : InputParameters params_sub = _factory.getValidParams("CoupledBEEquilibriumSub");
208 2408 : params_sub.set<NonlinearVariableName>("variable") = _primary_species[i];
209 1204 : params_sub.set<Real>("weight") = _weights[i][j];
210 2408 : params_sub.defaultCoupledValue("log_k", _eq_const[j]);
211 1204 : params_sub.set<Real>("sto_u") = _sto_u[i][j];
212 2408 : params_sub.set<std::vector<Real>>("sto_v") = _sto_v[i][j];
213 1204 : params_sub.set<std::vector<VariableName>>("v") = _coupled_v[i][j];
214 3612 : _problem->addKernel("CoupledBEEquilibriumSub",
215 1204 : _primary_species[i] + "_" + _eq_species[j] + "_sub",
216 : params_sub);
217 :
218 2408 : InputParameters params_cd = _factory.getValidParams("CoupledDiffusionReactionSub");
219 2408 : params_cd.set<NonlinearVariableName>("variable") = _primary_species[i];
220 1204 : params_cd.set<Real>("weight") = _weights[i][j];
221 2408 : params_cd.defaultCoupledValue("log_k", _eq_const[j]);
222 1204 : params_cd.set<Real>("sto_u") = _sto_u[i][j];
223 2408 : params_cd.set<std::vector<Real>>("sto_v") = _sto_v[i][j];
224 1204 : params_cd.set<std::vector<VariableName>>("v") = _coupled_v[i][j];
225 3612 : _problem->addKernel("CoupledDiffusionReactionSub",
226 1204 : _primary_species[i] + "_" + _eq_species[j] + "_cd",
227 : params_cd);
228 :
229 : // If pressure is coupled, add a CoupledConvectionReactionSub Kernel as well
230 2408 : if (_pars.isParamValid("pressure"))
231 : {
232 2408 : InputParameters params_conv = _factory.getValidParams("CoupledConvectionReactionSub");
233 2408 : params_conv.set<NonlinearVariableName>("variable") = _primary_species[i];
234 1204 : params_conv.set<Real>("weight") = _weights[i][j];
235 2408 : params_conv.defaultCoupledValue("log_k", _eq_const[j]);
236 1204 : params_conv.set<Real>("sto_u") = _sto_u[i][j];
237 2408 : params_conv.set<std::vector<Real>>("sto_v") = _sto_v[i][j];
238 1204 : params_conv.set<std::vector<VariableName>>("v") = _coupled_v[i][j];
239 1204 : params_conv.set<std::vector<VariableName>>("p") = _pressure_var;
240 1204 : params_conv.set<RealVectorValue>("gravity") = _gravity;
241 3612 : _problem->addKernel("CoupledConvectionReactionSub",
242 1204 : _primary_species[i] + "_" + _eq_species[j] + "_conv",
243 : params_conv);
244 1204 : }
245 1204 : }
246 : }
247 : }
248 : }
249 :
250 402 : if (_current_task == "add_aux_kernel")
251 : {
252 : // Add AqueousEquilibriumRxnAux AuxKernels for equilibrium species
253 868 : for (unsigned int j = 0; j < _num_reactions; ++j)
254 : {
255 667 : if (_aux_species.find(_eq_species[j]) != _aux_species.end())
256 : {
257 1334 : InputParameters params_eq = _factory.getValidParams("AqueousEquilibriumRxnAux");
258 1334 : params_eq.set<AuxVariableName>("variable") = _eq_species[j];
259 1334 : params_eq.defaultCoupledValue("log_k", _eq_const[j]);
260 1334 : params_eq.set<std::vector<Real>>("sto_v") = _stos[j];
261 667 : params_eq.set<std::vector<VariableName>>("v") = _primary_species_involved[j];
262 1334 : _problem->addAuxKernel("AqueousEquilibriumRxnAux", "aux_" + _eq_species[j], params_eq);
263 667 : }
264 : }
265 : }
266 402 : }
|
