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 "PorousFlowAqueousPreDisChemistry.h"
11 :
12 : registerMooseObject("PorousFlowApp", PorousFlowAqueousPreDisChemistry);
13 :
14 : InputParameters
15 2260 : PorousFlowAqueousPreDisChemistry::validParams()
16 : {
17 2260 : InputParameters params = PorousFlowMaterialVectorBase::validParams();
18 4520 : params.addRequiredCoupledVar(
19 : "primary_concentrations",
20 : "List of MOOSE Variables that represent the concentrations of the primary species");
21 4520 : params.addRequiredParam<unsigned>("num_reactions",
22 : "Number of equations in the system of chemical reactions");
23 4520 : params.addParam<bool>("equilibrium_constants_as_log10",
24 4520 : false,
25 : "If true, the equilibrium constants are written in their log10 form, eg, "
26 : "-2. If false, the equilibrium constants are written in absolute terms, "
27 : "eg, 0.01");
28 4520 : params.addRequiredCoupledVar("equilibrium_constants",
29 : "Equilibrium constant for each equation (dimensionless). If these "
30 : "are temperature dependent AuxVariables, the Jacobian will not be "
31 : "exact");
32 4520 : params.addRequiredParam<std::vector<Real>>(
33 : "primary_activity_coefficients",
34 : "Activity coefficients for the primary species (dimensionless) (one for each)");
35 4520 : params.addRequiredParam<std::vector<Real>>(
36 : "reactions",
37 : "A matrix defining the aqueous reactions. The matrix is entered as a long vector: the first "
38 : "row is "
39 : "entered first, followed by the second row, etc. There should be num_reactions rows. All "
40 : "primary species should appear only on the LHS of each reaction (and there should be just "
41 : "one secondary species on the RHS, by definition) so they may have negative coefficients. "
42 : "Each row should have number of primary_concentrations entries, which are the stoichiometric "
43 : "coefficients. The first coefficient must always correspond to the first primary species, "
44 : "etc");
45 4520 : params.addRequiredParam<std::vector<Real>>("specific_reactive_surface_area",
46 : "Specific reactive surface area in m^2/(L solution).");
47 4520 : params.addRequiredParam<std::vector<Real>>(
48 : "kinetic_rate_constant",
49 : "Kinetic rate constant in mol/(m^2 s), at the reference temperature (one for each reaction)");
50 4520 : params.addRequiredParam<std::vector<Real>>("molar_volume",
51 : "Volume occupied by one mole of the secondary species "
52 : "(L(solution)/mol) (one for each reaction)");
53 4520 : params.addRequiredParam<std::vector<Real>>("activation_energy",
54 : "Activation energy, J/mol (one for each reaction)");
55 4520 : params.addParam<Real>("gas_constant", 8.31434, "Gas constant, in J/(mol K)");
56 4520 : params.addParam<Real>("reference_temperature", 298.15, "Reference temperature, K");
57 4520 : params.addParam<std::vector<Real>>("theta_exponent",
58 : "Theta exponent. Defaults to 1. (one for each reaction)");
59 4520 : params.addParam<std::vector<Real>>("eta_exponent",
60 : "Eta exponent. Defaults to 1. (one for each reaction)");
61 4520 : params.addPrivateParam<std::string>("pf_material_type", "chemistry");
62 2260 : params.addClassDescription("This Material forms a std::vector of mineralisation reaction rates "
63 : "(L(precipitate)/L(solution)/s) appropriate to the aqueous "
64 : "precipitation-dissolution system provided. Note: the "
65 : "PorousFlowTemperature must be measured in Kelvin.");
66 2260 : return params;
67 0 : }
68 :
69 1777 : PorousFlowAqueousPreDisChemistry::PorousFlowAqueousPreDisChemistry(
70 1777 : const InputParameters & parameters)
71 : : PorousFlowMaterialVectorBase(parameters),
72 1074 : _porosity_old(_nodal_material ? getMaterialPropertyOld<Real>("PorousFlow_porosity_nodal")
73 3183 : : getMaterialPropertyOld<Real>("PorousFlow_porosity_qp")),
74 1777 : _aq_ph(_dictator.aqueousPhaseNumber()),
75 3554 : _saturation(_nodal_material
76 1777 : ? getMaterialProperty<std::vector<Real>>("PorousFlow_saturation_nodal")
77 3183 : : getMaterialProperty<std::vector<Real>>("PorousFlow_saturation_qp")),
78 :
79 1777 : _temperature(_nodal_material ? getMaterialProperty<Real>("PorousFlow_temperature_nodal")
80 3183 : : getMaterialProperty<Real>("PorousFlow_temperature_qp")),
81 1777 : _dtemperature_dvar(
82 1777 : _nodal_material
83 1777 : ? getMaterialProperty<std::vector<Real>>("dPorousFlow_temperature_nodal_dvar")
84 3183 : : getMaterialProperty<std::vector<Real>>("dPorousFlow_temperature_qp_dvar")),
85 :
86 1777 : _num_primary(coupledComponents("primary_concentrations")),
87 3554 : _num_reactions(getParam<unsigned>("num_reactions")),
88 3554 : _equilibrium_constants_as_log10(getParam<bool>("equilibrium_constants_as_log10")),
89 1777 : _num_equilibrium_constants(coupledComponents("equilibrium_constants")),
90 1777 : _equilibrium_constants(_num_equilibrium_constants),
91 3554 : _primary_activity_coefficients(getParam<std::vector<Real>>("primary_activity_coefficients")),
92 5331 : _reactions(getParam<std::vector<Real>>("reactions")),
93 :
94 1777 : _sec_conc_old(
95 1777 : _nodal_material
96 1777 : ? getMaterialPropertyOld<std::vector<Real>>("PorousFlow_mineral_concentration_nodal")
97 3183 : : getMaterialPropertyOld<std::vector<Real>>("PorousFlow_mineral_concentration_qp")),
98 :
99 1777 : _mineral_sat(_num_reactions),
100 1777 : _bounded_rate(_num_reactions),
101 1777 : _reaction_rate(
102 1777 : _nodal_material
103 1777 : ? declareProperty<std::vector<Real>>("PorousFlow_mineral_reaction_rate_nodal")
104 2480 : : declareProperty<std::vector<Real>>("PorousFlow_mineral_reaction_rate_qp")),
105 1777 : _dreaction_rate_dvar(_nodal_material ? declareProperty<std::vector<std::vector<Real>>>(
106 : "dPorousFlow_mineral_reaction_rate_nodal_dvar")
107 2480 : : declareProperty<std::vector<std::vector<Real>>>(
108 : "dPorousFlow_mineral_reaction_rate_qp_dvar")),
109 :
110 3554 : _r_area(getParam<std::vector<Real>>("specific_reactive_surface_area")),
111 3554 : _molar_volume(getParam<std::vector<Real>>("molar_volume")),
112 3554 : _ref_kconst(getParam<std::vector<Real>>("kinetic_rate_constant")),
113 3554 : _e_act(getParam<std::vector<Real>>("activation_energy")),
114 3554 : _gas_const(getParam<Real>("gas_constant")),
115 3554 : _one_over_ref_temp(1.0 / getParam<Real>("reference_temperature")),
116 5762 : _theta_exponent(isParamValid("theta_exponent") ? getParam<std::vector<Real>>("theta_exponent")
117 1777 : : std::vector<Real>(_num_reactions, 1.0)),
118 5762 : _eta_exponent(isParamValid("eta_exponent") ? getParam<std::vector<Real>>("eta_exponent")
119 3554 : : std::vector<Real>(_num_reactions, 1.0))
120 : {
121 1777 : if (_dictator.numPhases() < 1)
122 2 : mooseError("PorousFlowAqueousPreDisChemistry: The number of fluid phases must not be zero");
123 :
124 1775 : if (_num_primary != _num_components - 1)
125 0 : mooseError("PorousFlowAqueousPreDisChemistry: The number of mass_fraction_vars is ",
126 0 : _num_components,
127 : " which must be one greater than the number of primary concentrations (which is ",
128 0 : _num_primary,
129 : ")");
130 :
131 : // correct number of equilibrium constants
132 1775 : if (_num_equilibrium_constants != _num_reactions)
133 0 : mooseError("PorousFlowAqueousPreDisChemistry: The number of equilibrium constants is ",
134 0 : _num_equilibrium_constants,
135 : " which must be equal to the number of reactions (",
136 0 : _num_reactions,
137 : ")");
138 :
139 : // correct number of activity coefficients
140 1775 : if (_primary_activity_coefficients.size() != _num_primary)
141 0 : mooseError("PorousFlowAqueousPreDisChemistry: The number of primary activity "
142 : "coefficients is ",
143 0 : _primary_activity_coefficients.size(),
144 : " which must be equal to the number of primary species (",
145 0 : _num_primary,
146 : ")");
147 :
148 : // correct number of stoichiometry coefficients
149 1775 : if (_reactions.size() != _num_reactions * _num_primary)
150 0 : mooseError("PorousFlowAqueousPreDisChemistry: The number of stoichiometric "
151 : "coefficients specified in 'reactions' (",
152 0 : _reactions.size(),
153 : ") must be equal to the number of reactions (",
154 0 : _num_reactions,
155 : ") multiplied by the number of primary species (",
156 0 : _num_primary,
157 : ")");
158 :
159 1775 : if (_r_area.size() != _num_reactions)
160 2 : mooseError("PorousFlowAqueousPreDisChemistry: The number of specific reactive "
161 : "surface areas provided is ",
162 2 : _r_area.size(),
163 : " which must be equal to the number of reactions (",
164 2 : _num_reactions,
165 : ")");
166 :
167 1773 : if (_ref_kconst.size() != _num_reactions)
168 2 : mooseError("PorousFlowAqueousPreDisChemistry: The number of kinetic rate constants is ",
169 2 : _ref_kconst.size(),
170 : " which must be equal to the number of reactions (",
171 2 : _num_reactions,
172 : ")");
173 :
174 1771 : if (_e_act.size() != _num_reactions)
175 2 : mooseError("PorousFlowAqueousPreDisChemistry: The number of activation energies is ",
176 2 : _e_act.size(),
177 : " which must be equal to the number of reactions (",
178 2 : _num_reactions,
179 : ")");
180 :
181 1769 : if (_molar_volume.size() != _num_reactions)
182 2 : mooseError("PorousFlowAqueousPreDisChemistry: The number of molar volumes is ",
183 2 : _molar_volume.size(),
184 : " which must be equal to the number of reactions (",
185 2 : _num_reactions,
186 : ")");
187 :
188 1767 : if (_theta_exponent.size() != _num_reactions)
189 2 : mooseError("PorousFlowAqueousPreDisChemistry: The number of theta exponents is ",
190 2 : _theta_exponent.size(),
191 : " which must be equal to the number of reactions (",
192 2 : _num_reactions,
193 : ")");
194 :
195 1765 : if (_eta_exponent.size() != _num_reactions)
196 2 : mooseError("PorousFlowAqueousPreDisChemistry: The number of eta exponents is ",
197 2 : _eta_exponent.size(),
198 : " which must be equal to the number of reactions (",
199 2 : _num_reactions,
200 : ")");
201 :
202 1763 : if (_num_reactions != _dictator.numAqueousKinetic())
203 4 : mooseError("PorousFlowAqueousPreDisChemistry: You have specified the number of "
204 : "reactions to be ",
205 2 : _num_reactions,
206 : " but the Dictator knows that the number of aqueous kinetic "
207 : "(precipitation-dissolution) reactions is ",
208 2 : _dictator.numAqueousKinetic());
209 :
210 1761 : _primary_var_num.resize(_num_primary);
211 1761 : _primary.resize(_num_primary);
212 4869 : for (unsigned i = 0; i < _num_primary; ++i)
213 : {
214 3108 : _primary_var_num[i] = coupled("primary_concentrations", i);
215 6216 : _primary[i] = (_nodal_material ? &coupledDofValues("primary_concentrations", i)
216 4200 : : &coupledValue("primary_concentrations", i));
217 : }
218 :
219 3852 : for (unsigned i = 0; i < _num_equilibrium_constants; ++i)
220 : {
221 : // If equilibrium_constants are elemental AuxVariables (or constants), we want to use
222 : // coupledGenericValue() rather than coupledGenericDofValue()
223 2091 : const bool is_nodal = isCoupled("equilibrium_constants")
224 4182 : ? getFieldVar("equilibrium_constants", i)->isNodal()
225 : : false;
226 :
227 2091 : _equilibrium_constants[i] =
228 4182 : (_nodal_material && is_nodal ? &coupledDofValues("equilibrium_constants", i)
229 2844 : : &coupledValue("equilibrium_constants", i));
230 : }
231 1761 : }
232 :
233 : void
234 22666 : PorousFlowAqueousPreDisChemistry::initQpStatefulProperties()
235 : {
236 22666 : _reaction_rate[_qp].resize(_num_reactions);
237 22666 : _dreaction_rate_dvar[_qp].resize(_num_reactions);
238 45532 : for (unsigned r = 0; r < _num_reactions; ++r)
239 22866 : _dreaction_rate_dvar[_qp][r].assign(_num_var, 0.0);
240 22666 : }
241 :
242 : void
243 951658 : PorousFlowAqueousPreDisChemistry::computeQpProperties()
244 : {
245 951658 : _reaction_rate[_qp].resize(_num_reactions);
246 951658 : _dreaction_rate_dvar[_qp].resize(_num_reactions);
247 1904556 : for (unsigned r = 0; r < _num_reactions; ++r)
248 952898 : _dreaction_rate_dvar[_qp][r].assign(_num_var, 0.0);
249 :
250 : // Compute the reaction rates
251 951658 : computeQpReactionRates();
252 :
253 : // Compute the derivatives of the reaction rates
254 951658 : std::vector<std::vector<Real>> drr(_num_reactions);
255 951658 : std::vector<Real> drr_dT(_num_reactions);
256 1904556 : for (unsigned r = 0; r < _num_reactions; ++r)
257 : {
258 952898 : dQpReactionRate_dprimary(r, drr[r]);
259 952898 : drr_dT[r] = dQpReactionRate_dT(r);
260 : }
261 :
262 : // compute _dreaction_rate_dvar[_qp]
263 4037576 : for (unsigned wrt = 0; wrt < _num_primary; ++wrt)
264 : {
265 : // derivative with respect to the "wrt"^th primary species concentration
266 3085918 : if (!_dictator.isPorousFlowVariable(_primary_var_num[wrt]))
267 222 : continue;
268 3085696 : const unsigned pf_wrt = _dictator.porousFlowVariableNum(_primary_var_num[wrt]);
269 :
270 : // run through the reactions, using drr in the appropriate places
271 6174972 : for (unsigned r = 0; r < _num_reactions; ++r)
272 3089276 : _dreaction_rate_dvar[_qp][r][pf_wrt] = drr[r][wrt];
273 : }
274 :
275 : // use the derivative wrt temperature
276 1904556 : for (unsigned r = 0; r < _num_reactions; ++r)
277 4044506 : for (unsigned v = 0; v < _num_var; ++v)
278 3091608 : _dreaction_rate_dvar[_qp][r][v] += drr_dT[r] * _dtemperature_dvar[_qp][v];
279 951658 : }
280 :
281 : Real
282 5519564 : PorousFlowAqueousPreDisChemistry::stoichiometry(unsigned reaction_num, unsigned primary_num) const
283 : {
284 5519564 : const unsigned index = reaction_num * _num_primary + primary_num;
285 5519564 : return _reactions[index];
286 : }
287 :
288 : void
289 615810 : PorousFlowAqueousPreDisChemistry::findZeroConcentration(unsigned & zero_conc_index,
290 : unsigned & zero_count) const
291 : {
292 615810 : zero_count = 0;
293 3043912 : for (unsigned i = 0; i < _num_primary; ++i)
294 : {
295 2428592 : if (_primary_activity_coefficients[i] * (*_primary[i])[_qp] <= 0.0)
296 : {
297 1350 : zero_count += 1;
298 1350 : zero_conc_index = i;
299 1350 : if (zero_count > 1)
300 : return;
301 : }
302 : }
303 : return;
304 : }
305 :
306 : void
307 951658 : PorousFlowAqueousPreDisChemistry::computeQpReactionRates()
308 : {
309 1904556 : for (unsigned r = 0; r < _num_reactions; ++r)
310 : {
311 952898 : _mineral_sat[r] =
312 952898 : (_equilibrium_constants_as_log10 ? std::pow(10.0, -(*_equilibrium_constants[r])[_qp])
313 952378 : : 1.0 / (*_equilibrium_constants[r])[_qp]);
314 3998490 : for (unsigned j = 0; j < _num_primary; ++j)
315 : {
316 3068812 : const Real gamp = _primary_activity_coefficients[j] * (*_primary[j])[_qp];
317 3068812 : if (gamp <= 0.0)
318 : {
319 23400 : if (stoichiometry(r, j) < 0.0)
320 20 : _mineral_sat[r] = std::numeric_limits<Real>::max();
321 23380 : else if (stoichiometry(r, j) == 0.0)
322 : _mineral_sat[r] *= 1.0;
323 : else
324 : {
325 23220 : _mineral_sat[r] = 0.0;
326 23220 : break;
327 : }
328 : }
329 : else
330 3045412 : _mineral_sat[r] *= std::pow(gamp, stoichiometry(r, j));
331 : }
332 952898 : const Real fac = 1.0 - std::pow(_mineral_sat[r], _theta_exponent[r]);
333 : // if fac > 0 then dissolution occurs; if fac < 0 then precipitation occurs.
334 952898 : const Real sgn = (fac < 0 ? -1.0 : 1.0);
335 952898 : const Real unbounded_rr = -sgn * rateConstantQp(r) * _r_area[r] * _molar_volume[r] *
336 952898 : std::pow(std::abs(fac), _eta_exponent[r]);
337 :
338 : /*
339 : *
340 : * Note the use of the OLD value of porosity here.
341 : * This strategy, which breaks the cyclic dependency between porosity
342 : * and mineral concentration, is used in
343 : * Kernel: PorousFlowPreDis
344 : * Material: PorousFlowPorosity
345 : * Material: PorousFlowAqueousPreDisChemistry
346 : * Material: PorousFlowAqueousPreDisMineral
347 : *
348 : */
349 :
350 : // bound the reaction so _sec_conc lies between zero and unity
351 952898 : const Real por_times_rr_dt = _porosity_old[_qp] * _saturation[_qp][_aq_ph] * unbounded_rr * _dt;
352 952898 : if (_sec_conc_old[_qp][r] + por_times_rr_dt > 1.0)
353 : {
354 : _bounded_rate[r] = true;
355 250 : _reaction_rate[_qp][r] =
356 250 : (1.0 - _sec_conc_old[_qp][r]) / _porosity_old[_qp] / _saturation[_qp][_aq_ph] / _dt;
357 : }
358 952648 : else if (_sec_conc_old[_qp][r] + por_times_rr_dt < 0.0)
359 : {
360 : _bounded_rate[r] = true;
361 336838 : _reaction_rate[_qp][r] =
362 336838 : -_sec_conc_old[_qp][r] / _porosity_old[_qp] / _saturation[_qp][_aq_ph] / _dt;
363 : }
364 : else
365 : {
366 : _bounded_rate[r] = false;
367 615810 : _reaction_rate[_qp][r] = unbounded_rr;
368 : }
369 : }
370 951658 : }
371 :
372 : void
373 952898 : PorousFlowAqueousPreDisChemistry::dQpReactionRate_dprimary(unsigned reaction_num,
374 : std::vector<Real> & drr) const
375 : {
376 952898 : drr.assign(_num_primary, 0.0);
377 :
378 : // handle corner case
379 952898 : if (_bounded_rate[reaction_num])
380 337758 : return;
381 :
382 : /*
383 : * Form the derivative of _mineral_sat, and store it in drr for now.
384 : * The derivatives are straightforward if all primary > 0.
385 : *
386 : * If more than one primary = 0 then I set the derivatives to zero, even though it could be
387 : * argued that with certain stoichiometric coefficients you might have derivative = 0/0 and it
388 : * might be appropriate to set this to a non-zero finite value.
389 : *
390 : * If exactly one primary = 0 and its stoichiometry = 1 then the derivative wrt this one is
391 : * nonzero.
392 : * If exactly one primary = 0 and its stoichiometry > 1 then all derivatives are zero.
393 : * If exactly one primary = 0 and its stoichiometry < 1 then the derivative wrt this one is
394 : * infinity
395 : */
396 :
397 615810 : unsigned zero_count = 0;
398 615810 : unsigned zero_conc_index = 0;
399 615810 : findZeroConcentration(zero_conc_index, zero_count);
400 615810 : if (zero_count == 0)
401 : {
402 3041752 : for (unsigned i = 0; i < _num_primary; ++i)
403 2426802 : drr[i] = stoichiometry(reaction_num, i) * _mineral_sat[reaction_num] /
404 2426802 : std::max((*_primary[i])[_qp], 0.0);
405 : }
406 : else
407 : {
408 860 : if (_theta_exponent[reaction_num] < 1.0)
409 : // dfac = infinity (see below) so the derivative may be inf, inf * 0, or inf/inf. I simply
410 : // return with drr = 0
411 : return;
412 :
413 : // count the number of primary <= 0, and record the one that's zero
414 720 : if (zero_count == 1 && stoichiometry(reaction_num, zero_conc_index) == 1.0)
415 : {
416 190 : Real conc_without_zero = (_equilibrium_constants_as_log10
417 190 : ? std::pow(10.0, -(*_equilibrium_constants[reaction_num])[_qp])
418 190 : : 1.0 / (*_equilibrium_constants[reaction_num])[_qp]);
419 580 : for (unsigned i = 0; i < _num_primary; ++i)
420 : {
421 390 : if (i == zero_conc_index)
422 190 : conc_without_zero *= _primary_activity_coefficients[i];
423 : else
424 200 : conc_without_zero *=
425 200 : std::pow(_primary_activity_coefficients[i] * std::max((*_primary[i])[_qp], 0.0),
426 : stoichiometry(reaction_num, i));
427 : }
428 190 : drr[zero_conc_index] = conc_without_zero;
429 : }
430 530 : else if (zero_count == 0 and stoichiometry(reaction_num, zero_conc_index) < 1.0)
431 0 : drr[zero_conc_index] = std::numeric_limits<Real>::max();
432 : else
433 : // all other cases have drr = 0, so return without performing any other calculations
434 530 : return;
435 : }
436 :
437 : // At the moment _drr = d(mineral_sat)/d(primary)
438 : // Multiply by appropriate term to give _drr = d(reaction_rate)/d(primary)
439 615140 : const Real fac = 1.0 - std::pow(_mineral_sat[reaction_num], _theta_exponent[reaction_num]);
440 615140 : const Real dfac = -_theta_exponent[reaction_num] *
441 615140 : std::pow(_mineral_sat[reaction_num], _theta_exponent[reaction_num] - 1.0);
442 615140 : const Real multiplier = -rateConstantQp(reaction_num) * _r_area[reaction_num] *
443 615140 : _molar_volume[reaction_num] *
444 615140 : std::pow(std::abs(fac), _eta_exponent[reaction_num] - 1.0) *
445 615140 : _eta_exponent[reaction_num] * dfac;
446 3042332 : for (unsigned i = 0; i < _num_primary; ++i)
447 2427192 : drr[i] *= multiplier;
448 : }
449 :
450 : Real
451 2183848 : PorousFlowAqueousPreDisChemistry::rateConstantQp(unsigned reaction_num) const
452 : {
453 2183848 : return _ref_kconst[reaction_num] * std::exp(_e_act[reaction_num] / _gas_const *
454 2183848 : (_one_over_ref_temp - 1.0 / _temperature[_qp]));
455 : }
456 :
457 : Real
458 952898 : PorousFlowAqueousPreDisChemistry::dQpReactionRate_dT(unsigned reaction_num) const
459 : {
460 : // handle corner case
461 952898 : if (_bounded_rate[reaction_num])
462 : return 0.0;
463 :
464 615810 : const Real drate_const = rateConstantQp(reaction_num) * _e_act[reaction_num] / _gas_const *
465 615810 : std::pow(_temperature[_qp], -2.0);
466 615810 : const Real fac = 1.0 - std::pow(_mineral_sat[reaction_num], _theta_exponent[reaction_num]);
467 615810 : const Real sgn = (fac < 0 ? -1.0 : 1.0);
468 615810 : const Real dkinetic_rate = -sgn * drate_const * _r_area[reaction_num] *
469 : _molar_volume[reaction_num] *
470 615810 : std::pow(std::abs(fac), _eta_exponent[reaction_num]);
471 :
472 615810 : return dkinetic_rate;
473 : }
|
