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 "AqueousReactionsEquilibriumPhysics.h"
11 : #include "ReactionNetworkUtils.h"
12 :
13 : // Register the actions for the objects actually used
14 : registerMooseAction("ChemicalReactionsApp", AqueousReactionsEquilibriumPhysics, "add_kernel");
15 : registerMooseAction("ChemicalReactionsApp", AqueousReactionsEquilibriumPhysics, "add_aux_kernel");
16 : registerReactionNetworkPhysicsBaseTasks("ChemicalReactionsApp", AqueousReactionsEquilibriumPhysics);
17 :
18 : InputParameters
19 36 : AqueousReactionsEquilibriumPhysics::validParams()
20 : {
21 36 : InputParameters params = ReactionNetworkPhysicsBase::validParams();
22 36 : params.addClassDescription("Forms the equations for the chemical reaction networks in a fluid"
23 : " medium using a continuous Galerkin finite element discretization.");
24 :
25 : // Rename parameters to match the previously existing actions for AqueousEquilibrium
26 : // ReactionNetwork
27 72 : params.renameParam("solver_variables", "primary_species", "The list of primary species to add");
28 72 : params.renameParam(
29 : "auxiliary_variables", "secondary_species", "The list of secondary species to add");
30 72 : params.renameParam(
31 : "variable_order", "order", "Order of both the primary and secondary species variables");
32 72 : params.renameParam("equation_scaling", "scaling", "");
33 72 : params.setDocString("reactions", "The list of equilibrium reactions occuring in the fluid");
34 :
35 : // To preserve the legacy option to perform Darcy-advection with AqueousEquilibrium
36 : // ReactionNetwork
37 72 : params.addParam<bool>("add_darcy_advection_term",
38 72 : false,
39 : "Whether to add an advection term using the Darcy equation to compute the "
40 : "advecting velocity");
41 72 : params.addParam<std::vector<VariableName>>(
42 : "pressure", {}, "Pressure variable (for Darcy advection)");
43 36 : params.addParam<RealVectorValue>(
44 36 : "gravity", RealVectorValue(0, 0, -9.81), "Gravity vector (for Darcy advection)");
45 :
46 36 : return params;
47 0 : }
48 :
49 36 : AqueousReactionsEquilibriumPhysics::AqueousReactionsEquilibriumPhysics(
50 36 : const InputParameters & parameters)
51 : : ReactionNetworkPhysicsBase(parameters),
52 36 : _sto_u(_num_solver_species),
53 36 : _sto_v(_num_solver_species),
54 36 : _weights(_num_solver_species),
55 36 : _primary_participation(_num_solver_species),
56 36 : _coupled_v(_num_solver_species),
57 108 : _pressure_var(getParam<std::vector<VariableName>>("pressure"))
58 : {
59 : // Further parse the reactions
60 : // Keeping the data in these vectors of vectors as an intermediate processing step
61 99 : for (const auto i : index_range(_reactions))
62 : {
63 : const auto & reaction = _reactions[i];
64 126 : _solver_species_involved.push_back(reaction.getReactantSpecies());
65 126 : _stos.push_back(reaction.getStoichiometricCoefficients());
66 :
67 : // There are only one equilibrium species. We look at the output species
68 63 : const auto & products = reaction.getProductSpecies();
69 63 : if (products.size() != 1)
70 0 : mooseError("Reaction:\n" + _reactions_input[i] +
71 : "\n has more than one product species (on the RHS). This Physics only supports "
72 : "one product species per reaction");
73 63 : _eq_species.push_back(products[0]);
74 :
75 : // If user passed log_K use that, else use K
76 126 : if (reaction.hasMetaData<Real>("log10_K"))
77 162 : _log_eq_const.push_back(std::stod(reaction.getMetaData("log10_K")));
78 18 : else if (reaction.hasMetaData<Real>("K"))
79 27 : _log_eq_const.push_back(std::log10(std::stod(reaction.getMetaData("K"))));
80 0 : else if (reaction.hasMetaData("K") || reaction.hasMetaData("log10_K"))
81 0 : paramError("reactions",
82 : "Equilibrium constant species in square brackets must be specified as a float, "
83 : "not a name");
84 : else
85 0 : paramError("reactions",
86 0 : "Reaction: '" + _reactions_input[i] +
87 : "'\n is missing an equilibrium constant [K=] or its logarithm [log10_K=] in "
88 : "its metadata");
89 63 : }
90 :
91 : // For each primary/solver species, we examine the reactions to get the coefficients
92 : // and various constants
93 90 : for (unsigned int i = 0; i < _solver_species.size(); ++i)
94 : {
95 54 : _sto_u[i].resize(_num_reactions, 0.0);
96 54 : _sto_v[i].resize(_num_reactions);
97 54 : _coupled_v[i].resize(_num_reactions);
98 54 : _weights[i].resize(_num_reactions, 0.0);
99 :
100 54 : _primary_participation[i].resize(_num_reactions, false);
101 162 : for (unsigned int j = 0; j < _num_reactions; ++j)
102 : {
103 : // Set the booleans for involvement in a reaction
104 270 : for (unsigned int k = 0; k < _solver_species_involved[j].size(); ++k)
105 162 : if (_solver_species_involved[j][k] == _solver_species[i])
106 : _primary_participation[i][j] = true;
107 :
108 108 : if (_primary_participation[i][j])
109 : {
110 234 : for (unsigned int k = 0; k < _solver_species_involved[j].size(); ++k)
111 : {
112 : // Re-sort the coefficients into a more convenient format
113 144 : if (_solver_species_involved[j][k] == _solver_species[i])
114 : {
115 90 : _sto_u[i][j] = _stos[j][k];
116 90 : _weights[i][j] = _stos[j][k];
117 : }
118 : else
119 : {
120 54 : _sto_v[i][j].push_back(_stos[j][k]);
121 54 : _coupled_v[i][j].push_back(_solver_species_involved[j][k]);
122 : }
123 : }
124 : }
125 : }
126 : }
127 :
128 : // Parameter checks
129 72 : checkSecondParamSetOnlyIfFirstOneTrue("add_darcy_advection_term", "pressure");
130 72 : checkSecondParamSetOnlyIfFirstOneTrue("add_darcy_advection_term", "gravity");
131 36 : }
132 :
133 : void
134 36 : AqueousReactionsEquilibriumPhysics::addFEKernels()
135 : {
136 : // Add Kernels for each primary species
137 : // Note that the equations are on a per-species basis!
138 : // So equations are organized differently than reactions
139 90 : for (const auto i : index_range(_solver_species))
140 : {
141 162 : for (const auto j : make_range(_num_reactions))
142 : {
143 108 : if (_primary_participation[i][j])
144 : {
145 : {
146 90 : InputParameters params_sub = _factory.getValidParams("CoupledBEEquilibriumSub");
147 90 : assignBlocks(params_sub, _blocks);
148 180 : params_sub.set<NonlinearVariableName>("variable") = _solver_species[i];
149 90 : params_sub.set<Real>("weight") = _weights[i][j];
150 180 : params_sub.defaultCoupledValue("log_k", _log_eq_const[j]);
151 90 : params_sub.set<Real>("sto_u") = _sto_u[i][j];
152 180 : params_sub.set<std::vector<Real>>("sto_v") = _sto_v[i][j];
153 90 : params_sub.set<std::vector<VariableName>>("v") = _coupled_v[i][j];
154 270 : _problem->addKernel("CoupledBEEquilibriumSub",
155 90 : _solver_species[i] + "_" + _eq_species[j] + "_sub",
156 : params_sub);
157 90 : }
158 :
159 : {
160 90 : InputParameters params_cd = _factory.getValidParams("CoupledDiffusionReactionSub");
161 90 : assignBlocks(params_cd, _blocks);
162 180 : params_cd.set<NonlinearVariableName>("variable") = _solver_species[i];
163 90 : params_cd.set<Real>("weight") = _weights[i][j];
164 180 : params_cd.defaultCoupledValue("log_k", _log_eq_const[j]);
165 90 : params_cd.set<Real>("sto_u") = _sto_u[i][j];
166 180 : params_cd.set<std::vector<Real>>("sto_v") = _sto_v[i][j];
167 90 : params_cd.set<std::vector<VariableName>>("v") = _coupled_v[i][j];
168 270 : _problem->addKernel("CoupledDiffusionReactionSub",
169 90 : _solver_species[i] + "_" + _eq_species[j] + "_cd",
170 : params_cd);
171 90 : }
172 :
173 : // If pressure is coupled, add a CoupledConvectionReactionSub Kernel as well
174 360 : if (getParam<bool>("add_darcy_advection_term") && isParamValid("pressure"))
175 : {
176 90 : InputParameters params_conv = _factory.getValidParams("CoupledConvectionReactionSub");
177 90 : assignBlocks(params_conv, _blocks);
178 180 : params_conv.set<NonlinearVariableName>("variable") = _solver_species[i];
179 90 : params_conv.set<Real>("weight") = _weights[i][j];
180 180 : params_conv.defaultCoupledValue("log_k", _log_eq_const[j]);
181 90 : params_conv.set<Real>("sto_u") = _sto_u[i][j];
182 180 : params_conv.set<std::vector<Real>>("sto_v") = _sto_v[i][j];
183 90 : params_conv.set<std::vector<VariableName>>("v") = _coupled_v[i][j];
184 90 : params_conv.set<std::vector<VariableName>>("p") = _pressure_var;
185 180 : params_conv.set<RealVectorValue>("gravity") = getParam<RealVectorValue>("gravity");
186 270 : _problem->addKernel("CoupledConvectionReactionSub",
187 90 : _solver_species[i] + "_" + _eq_species[j] + "_conv",
188 : params_conv);
189 90 : }
190 : }
191 : }
192 : }
193 36 : }
194 :
195 : void
196 36 : AqueousReactionsEquilibriumPhysics::addAuxiliaryKernels()
197 : {
198 : // Add AqueousEquilibriumRxnAux AuxKernels for equilibrium species
199 99 : for (const auto j : make_range(_num_reactions))
200 : {
201 : // Add these aux-kernels only for the aux species involved in the reaction
202 : // we should not be adding them twice, since only 1 eq_species per reaction
203 63 : if (std::find(_aux_species.begin(), _aux_species.end(), _eq_species[j]) != _aux_species.end())
204 : {
205 63 : InputParameters params_eq = _factory.getValidParams("AqueousEquilibriumRxnAux");
206 63 : assignBlocks(params_eq, _blocks);
207 126 : params_eq.set<AuxVariableName>("variable") = _eq_species[j];
208 126 : params_eq.defaultCoupledValue("log_k", _log_eq_const[j]);
209 : mooseAssert(_stos[j].size() >= _solver_species_involved[j].size(),
210 : "Coefs are for solver + auxiliary");
211 : std::vector<Real> stos_primary(_stos[j].begin(),
212 63 : _stos[j].begin() + _solver_species_involved[j].size());
213 126 : params_eq.set<std::vector<Real>>("sto_v") = stos_primary;
214 63 : params_eq.set<std::vector<VariableName>>("v") = _solver_species_involved[j];
215 126 : getProblem().addAuxKernel("AqueousEquilibriumRxnAux", "aux_" + _eq_species[j], params_eq);
216 63 : }
217 : }
218 36 : }
|
