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 : #pragma once
11 :
12 : #include "MooseTypes.h"
13 : #include "DenseMatrix.h"
14 : #include "GeochemistryConstants.h"
15 :
16 : /**
17 : * This controls the direction of a kinetic rate
18 : * BOTH: both dissolution and precipitation are allowed
19 : * PRECIPITATION: if log10(activity_product) < log10(equilibrium_constant) so that dissolution
20 : * should occur, then the rate is set to zero
21 : * DISSOLUTION: if log10(activity_product) > log10(equilibrium_constant) so that precipitation
22 : * should occur, then the rate is set to zero
23 : * RAW: rate is unimpacted by sign of log10(activity_product) / log10(equilibrium_constant)
24 : * DEATH: rate is unimpacted by sign of log10(activity_product) / log10(equilibrium_constant). This
25 : * is different from RAW because no reaction products are produced: the kinetic rate is *only* used
26 : * to modify the kinetic species' mass. This is useful to model biological mortality.
27 : */
28 : enum DirectionChoiceEnum
29 : {
30 : BOTH,
31 : DISSOLUTION,
32 : PRECIPITATION,
33 : RAW,
34 : DEATH
35 : };
36 :
37 : /**
38 : * Holds a user-specified description of a kinetic rate
39 : *
40 : * @param kinetic_species name of the kinetic species
41 : * @param intrinsic_rate_constant Note that intrinsic_rate_constant * area_quantity *
42 : * [kinetic_species mass] must have dimensions mol.s^-1
43 : * @param area_quantity Either 1, or the fixed surface area of the kinetic species, or a specific
44 : * surface area (m^2/g)
45 : * @param multiply_by_mass whether the rate should be multiplied by the kinetic_species mass
46 : * @param kinetic_molal_index, rate is multiplied by kinetic_species_molality^kinetic_molal_index
47 : * @param kinetic_monod_index, rate is multiplied by 1.0 /
48 : * (kinetic_species_molality^kinetic_molal_index +
49 : * kinetic_half_saturation^kinetic_molal_index)^kinetic_monod_index
50 : * @param kinetic_half_saturation, rate is multiplied by 1.0 /
51 : * (kinetic_species_molality^kinetic_molal_index +
52 : * kinetic_half_saturation^kinetic_molal_index)^kinetic_monod_index
53 : * @param promoting_species names of species (which must be primary or secondary species in the
54 : * system)
55 : * @param promoting_indices indices of mass, fugacity, activity or mobility (as appropriate)
56 : * @param promoting_monod_indices monod indices of mass, fugacity, activity or mobility (as
57 : * appropriate)
58 : * @param promoting_half_saturation half saturation values of all promoting species
59 : * @param theta exponent of (Q/K)
60 : * @param eta exponent of |1-(Q/K)^theta|
61 : * @param activation_energy in J.mol^-1
62 : * @param one_over_T0 measured in 1/Kelvin
63 : * @param direction: whether this kinetic rate is designed for: "both" precipitation and
64 : * dissolution; "precipitation" only; "dissolution" only; and "raw" or "death" mean the rate does
65 : * not depend on the sign of 1-(Q/K)
66 : * @param progeny: A non-kinetic species that catalyses the reaction, and potentially gets produced
67 : * or consumed by it
68 : * @param progeny_efficiency: When one mole of reaction is catalysed, progeny_efficiency moles of
69 : * the progeny is created
70 : * @param kinetic_bio_efficiency: the efficiency of a biologically-catalysed reaction, that is, when
71 : * one mole of reaction is catalysed, the biomass increases by kinetic_bio_efficiency moles
72 : * @param energy captured: energy captured by a biologically-catalysed reaction, essentially this
73 : * reduces the equilibrium constant of the reaction by exp(-energy_capture / R / Tk)
74 : */
75 : struct KineticRateUserDescription
76 : {
77 108 : KineticRateUserDescription(const std::string & kinetic_species_name,
78 : Real intrinsic_rate_constant,
79 : Real area_quantity,
80 : bool multiply_by_mass,
81 : Real kinetic_molal_index,
82 : Real kinetic_monod_index,
83 : Real kinetic_half_saturation,
84 : const std::vector<std::string> & promoting_species,
85 : const std::vector<Real> & promoting_indices,
86 : const std::vector<Real> & promoting_monod_indices,
87 : const std::vector<Real> & promoting_half_saturation,
88 : Real theta,
89 : Real eta,
90 : Real activation_energy,
91 : Real one_over_T0,
92 : DirectionChoiceEnum direction,
93 : std::string progeny,
94 : Real progeny_efficiency,
95 : Real kinetic_bio_efficiency,
96 : Real energy_captured)
97 108 : : kinetic_species_name(kinetic_species_name),
98 108 : intrinsic_rate_constant(intrinsic_rate_constant),
99 108 : area_quantity(area_quantity),
100 108 : multiply_by_mass(multiply_by_mass),
101 108 : kinetic_molal_index(kinetic_molal_index),
102 108 : kinetic_monod_index(kinetic_monod_index),
103 108 : kinetic_half_saturation(kinetic_half_saturation),
104 108 : promoting_species(promoting_species),
105 108 : promoting_indices(promoting_indices),
106 108 : promoting_monod_indices(promoting_monod_indices),
107 108 : promoting_half_saturation(promoting_half_saturation),
108 108 : theta(theta),
109 108 : eta(eta),
110 108 : activation_energy(activation_energy),
111 108 : one_over_T0(one_over_T0),
112 108 : direction(direction),
113 108 : progeny(progeny),
114 108 : progeny_efficiency(progeny_efficiency),
115 108 : kinetic_bio_efficiency(kinetic_bio_efficiency),
116 108 : energy_captured(energy_captured)
117 : {
118 108 : if (promoting_species.size() != promoting_indices.size())
119 1 : mooseError("The promoting_species and promoting_indices vectors must be the same size");
120 107 : if (promoting_species.size() != promoting_monod_indices.size())
121 1 : mooseError("The promoting_species and promoting_monod_indices vectors must be the same size");
122 106 : if (promoting_species.size() != promoting_half_saturation.size())
123 1 : mooseError(
124 : "The promoting_species and promoting_half_saturation vectors must be the same size");
125 : std::unordered_map<std::string, int> check_for_repeats;
126 230 : for (const std::string & name : promoting_species)
127 : if (check_for_repeats.count(name) == 1)
128 1 : mooseError("Promoting species ", name, " has already been provided with an exponent");
129 : else
130 125 : check_for_repeats[name] = 1;
131 108 : };
132 :
133 0 : bool operator==(const KineticRateUserDescription & rhs) const
134 : {
135 0 : return (kinetic_species_name == rhs.kinetic_species_name) &&
136 0 : (intrinsic_rate_constant == rhs.intrinsic_rate_constant) &&
137 0 : (area_quantity == rhs.area_quantity) && (multiply_by_mass == rhs.multiply_by_mass) &&
138 0 : (kinetic_molal_index == rhs.kinetic_molal_index) &&
139 0 : (kinetic_monod_index == rhs.kinetic_monod_index) &&
140 0 : (kinetic_half_saturation == rhs.kinetic_half_saturation) &&
141 0 : (promoting_species == rhs.promoting_species) &&
142 0 : (promoting_indices == rhs.promoting_indices) &&
143 0 : (promoting_monod_indices == rhs.promoting_monod_indices) &&
144 0 : (promoting_half_saturation == rhs.promoting_half_saturation) && (theta == rhs.theta) &&
145 0 : (eta == rhs.eta) && (activation_energy == rhs.activation_energy) &&
146 0 : (one_over_T0 == rhs.one_over_T0) && (direction == rhs.direction) &&
147 0 : (progeny == rhs.progeny) && (progeny_efficiency == rhs.progeny_efficiency) &&
148 0 : (kinetic_bio_efficiency == rhs.kinetic_bio_efficiency) &&
149 0 : (energy_captured == rhs.energy_captured);
150 : };
151 :
152 : std::string kinetic_species_name;
153 : Real intrinsic_rate_constant;
154 : Real area_quantity;
155 : bool multiply_by_mass;
156 : Real kinetic_molal_index;
157 : Real kinetic_monod_index;
158 : Real kinetic_half_saturation;
159 : std::vector<std::string> promoting_species;
160 : std::vector<Real> promoting_indices;
161 : std::vector<Real> promoting_monod_indices;
162 : std::vector<Real> promoting_half_saturation;
163 : Real theta;
164 : Real eta;
165 : Real activation_energy;
166 : Real one_over_T0;
167 : DirectionChoiceEnum direction;
168 : std::string progeny;
169 : Real progeny_efficiency;
170 : Real kinetic_bio_efficiency;
171 : Real energy_captured;
172 : };
173 :
174 : /**
175 : * Provides a parametric description of a general kinetic rate. This is a separate class with very
176 : * little functionality partly so that it can be used by multiple other classes
177 : * (PertinentGeochemicalSystem, GeochemistryKineticRate) but also so that it can be easily added to
178 : * if more general kinetic rates are ever required.
179 : *
180 : * The rate is a product of the following terms:
181 : *
182 : * intrinsic_rate_constant
183 : * area_quantity
184 : * mass of the kinetic_species, measured in grams, if multiply_by_mass is true
185 : * (molality of kinetic_species)^kinetic_molal_index
186 : * 1 / ((molality of kinetic_species)^kinetic_molal_index +
187 : * kinetic_half_saturation^kinetic_molal_index)^kinetic_monod_index
188 : * product over the promoting_species of m^(promoting_index) / (m^(promoting_index) +
189 : * half_saturation^(promoting_index))^(promoting_monod_index)
190 : * |1 - (Q/K)^theta|^eta, where K = eqm_constant_in_database * exp(- energy_captured / R / TK)
191 : * exp(activation_energy / R * (1/T0 - 1/T)) D(1 - (Q/K))
192 : *
193 : * Some explanation may be useful:
194 : *
195 : * intrinsic_rate_constant * area_quantity * [mass of kinetic_species] has units mol.s^-1. This
196 : * allows for the following common possibilities:
197 : * - intrinsic_rate_constant is measured in mol.s^-1. In this case, the user should set
198 : * area_quantity=1 and multiply_by_mass=false
199 : * - intrinsic_rate_constant is measured in mol.s^-1/kg(solvent_water). In this case the user
200 : * should set area_quantity=1 and multipliy_by_mass=false, but ensure that
201 : * promoting_species={"H2O",...} and promoting_index={1,...}.
202 : * - intrinsic_rate_constant is measured in mol.s^-1/area_of_kinetic_mineral, and the area of the
203 : * kinetic mineral is fixed. In this case the user should set
204 : * area_quantity=area_of_kinetic_mineral and multiply_by_mass=false.
205 : * - intrinsic_rate_constant is measured in mol.s^-1/area_of_kinetic_mineral, and the area of the
206 : * kinetic mineral depends on its mass. In this case, the user should set
207 : * area_quantity=specific_area_of_mineral (in m^2/g) and multiply_by_mass=true
208 : *
209 : * The "m" in the promoting species product is:
210 : * - mass of solvent water (in kg) if promoting_species="H2O"
211 : * - fugacity of a gas if promoting_species is a gas
212 : * - activity if promoting_species is either "H+" or "OH-"
213 : * - mobility, otherwise
214 : *
215 : * Q is the activity product, defined by the kinetic_species reaction (defined in the database
216 : * file).
217 : * K is the reaction's equilibrium constant (defined in the database file) multiplied by
218 : * exp(-energy_captured / (RT)).
219 : * R = 8.314472 m^2.kg.s^-2.K^-1.mol^-1 = 8.314472 J.K^-1.mol^-1 is
220 : * the gas constant. T is the temperature in Kelvin. T0 is a reference temperature, in Kelvin. It
221 : * is inputted as 1/T0 so that 1/T0 = 0 is possible.
222 : *
223 : * D(x) depends on direction. If direction == BOTH then D(x) = sgn(x). If direction == DISSOLUTION
224 : * then D(x) = (x>0)?1:0. If direction == PRECIPITATION then D(x) = (x<0)?-1:0. If direction ==
225 : * RAW then D=1. If direction == DEATH then D=1.
226 : *
227 : * The amount of kinetic species that is generated is kinetic_bio_efficiency * rate. Note that rate
228 : * > 0 for dissolution of the kinetic species, so kinetic_bio_efficiency defaults to -1. However,
229 : * for biogeochemistry, it is appropriate to set kinetic_bio_efficiency > 0, so that "dissolution"
230 : * of the biomass (with rate > 0) generates further biomass
231 : */
232 : namespace GeochemistryKineticRateCalculator
233 : {
234 : /**
235 : * Calclates a kinetic rate and its derivative. The rate is a product of the following terms:
236 : *
237 : * intrinsic_rate_constant
238 : * area_quantity
239 : * mass of the kinetic_species, measured in grams, if multiply_by_mass is true
240 : * product over the promoting_species of m^(promoting_index)
241 : * |1 - (Q/K)^theta|^eta
242 : * exp(activation_energy / R * (1/T0 - 1/T))
243 : * D(1 - (Q/K))
244 : *
245 : * @param promoting_indices of the basis species and equilibrium species. Note that this is
246 : * different from description.promoting_indices (which is paired with description.promoting_species)
247 : * because it is the promoting indices for the current basis and equilibrium. For computational
248 : * efficiency promoting_indices[0:num_basis-1] are the promoting indices for the current basis
249 : * species while promoting_indices[num_basis:] are the promoting indices for the current eqm species
250 : * @param description contains definition of intrinsic_rate_constant, area_quantity, etc
251 : * @param basis_species_name vector of basis species names
252 : * @param basis_species_gas i^th element is true if the i^th basis species is a gas
253 : * @param basis_molality vector of the basis molalities (zeroth element is kg of solvent water, and
254 : * this is mole number for mineral species)
255 : * @param basis_activity vector of basis activities
256 : * @param basis_activity_known i^th element is true if the activity for the i^th basis species has
257 : * been constrained by the user
258 : * @param eqm_species_name vector of eqm species names
259 : * @param eqm_species_gas i^th element is true if the i^th eqm species is a gas
260 : * @param eqm_molality vector of the eqm molalities (zeroth element is kg of solvent water, and this
261 : * is mole number for mineral species)
262 : * @param eqm_activity vector of eqm activities
263 : * @param eqm_stoichiometry matrix of stoichiometric coefficient
264 : * @param kin_moles moles of the kinetic species
265 : * @param kin_species_molecular_weight molecular weight (g/mol) of the kinetic species
266 : * @param kin_stoichiometry kinetic stoichiometry matrix (only row=kin is used)
267 : * @param kin the row of kin_stoichiometry that corresponds to the kinetic species
268 : * @param log10K log10(equilibrium constant) of the kinetic species
269 : * @param log10_activity_product log10(activity product) of the kinetic species
270 : * @param temp_deg temperature in degC
271 : * @param[out] rate the rate computed by this method
272 : * @param[out] drate_dkin d(rate)/d(mole number of the kinetic species)
273 : * @param[out] drate_dmol d(rate)/d(basis molality[i])
274 : */
275 : void calculateRate(const std::vector<Real> & promoting_indices,
276 : const std::vector<Real> & promoting_monod_indices,
277 : const std::vector<Real> & promoting_half_saturation,
278 : const KineticRateUserDescription & description,
279 : const std::vector<std::string> & basis_species_name,
280 : const std::vector<bool> & basis_species_gas,
281 : const std::vector<Real> & basis_molality,
282 : const std::vector<Real> & basis_activity,
283 : const std::vector<bool> & basis_activity_known,
284 : const std::vector<std::string> & eqm_species_name,
285 : const std::vector<bool> & eqm_species_gas,
286 : const std::vector<Real> & eqm_molality,
287 : const std::vector<Real> & eqm_activity,
288 : const DenseMatrix<Real> & eqm_stoichiometry,
289 : Real kin_moles,
290 : Real kin_species_molecular_weight,
291 : Real log10K,
292 : Real log10_activity_product,
293 : const DenseMatrix<Real> & kin_stoichiometry,
294 : unsigned kin,
295 : Real temp_degC,
296 : Real & rate,
297 : Real & drate_dkin,
298 : std::vector<Real> & drate_dmol);
299 : }
|
