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 "GeochemistryKineticRateCalculator.h"
11 :
12 : namespace GeochemistryKineticRateCalculator
13 : {
14 : void
15 5963 : calculateRate(const std::vector<Real> & promoting_indices,
16 : const std::vector<Real> & promoting_monod_indices,
17 : const std::vector<Real> & promoting_half_saturation,
18 : const KineticRateUserDescription & description,
19 : const std::vector<std::string> & basis_species_name,
20 : const std::vector<bool> & basis_species_gas,
21 : const std::vector<Real> & basis_molality,
22 : const std::vector<Real> & basis_activity,
23 : const std::vector<bool> & basis_activity_known,
24 : const std::vector<std::string> & eqm_species_name,
25 : const std::vector<bool> & eqm_species_gas,
26 : const std::vector<Real> & eqm_molality,
27 : const std::vector<Real> & eqm_activity,
28 : const DenseMatrix<Real> & eqm_stoichiometry,
29 : Real kin_moles,
30 : Real kin_species_molecular_weight,
31 : Real log10K,
32 : Real log10_activity_product,
33 : const DenseMatrix<Real> & kin_stoichiometry,
34 : unsigned kin,
35 : Real temp_degC,
36 : Real & rate,
37 : Real & drate_dkin,
38 : std::vector<Real> & drate_dmol)
39 : {
40 5963 : const unsigned num_basis = basis_species_name.size();
41 5963 : const unsigned num_eqm = eqm_species_name.size();
42 :
43 5963 : if (num_basis + num_eqm != promoting_indices.size())
44 1 : mooseError("kinetic_rate: promoting_indices incorrectly sized ", promoting_indices.size());
45 5962 : if (num_basis + num_eqm != promoting_monod_indices.size())
46 1 : mooseError("kinetic_rate: promoting_monod_indices incorrectly sized ",
47 1 : promoting_monod_indices.size());
48 5961 : if (num_basis + num_eqm != promoting_half_saturation.size())
49 1 : mooseError("kinetic_rate: promoting_half_saturation incorrectly sized ",
50 1 : promoting_half_saturation.size());
51 11916 : if (!(num_basis == basis_species_gas.size() && num_basis == basis_molality.size() &&
52 5957 : num_basis == basis_activity.size() && num_basis == basis_activity_known.size() &&
53 : num_basis == drate_dmol.size()))
54 5 : mooseError("kinetic_rate: incorrectly sized basis-species vectors ",
55 : num_basis,
56 : " ",
57 5 : basis_species_gas.size(),
58 : " ",
59 5 : basis_molality.size(),
60 : " ",
61 5 : basis_activity.size(),
62 : " ",
63 5 : basis_activity_known.size(),
64 : " ",
65 5 : drate_dmol.size());
66 5955 : if (!(num_eqm == eqm_species_gas.size() && num_eqm == eqm_molality.size() &&
67 : num_eqm == eqm_activity.size()))
68 3 : mooseError("kinetic_rate: incorrectly sized equilibrium-species vectors ",
69 : num_eqm,
70 : " ",
71 3 : eqm_species_gas.size(),
72 : " ",
73 3 : eqm_molality.size(),
74 : " ",
75 3 : eqm_activity.size());
76 5952 : if (!(eqm_stoichiometry.m() == num_eqm && eqm_stoichiometry.n() == num_basis))
77 2 : mooseError("kinetic_rate: incorrectly sized eqm stoichiometry matrix ",
78 2 : eqm_stoichiometry.m(),
79 : " ",
80 2 : eqm_stoichiometry.n());
81 5950 : if (!(kin_stoichiometry.m() > kin && kin_stoichiometry.n() == num_basis))
82 2 : mooseError("kinetic_rate: incorrectly sized kinetic stoichiometry matrix ",
83 2 : kin_stoichiometry.m(),
84 : " ",
85 2 : kin_stoichiometry.n());
86 :
87 5948 : rate = description.intrinsic_rate_constant;
88 5948 : rate *= description.area_quantity;
89 5948 : if (description.multiply_by_mass)
90 5600 : rate *= kin_moles * kin_species_molecular_weight;
91 :
92 5948 : const Real kin_molality = kin_moles / basis_molality[0];
93 5948 : const Real dkin_molality_dkin_moles = 1.0 / basis_molality[0];
94 5948 : const Real dkin_molality_dnw = -kin_molality / basis_molality[0];
95 5948 : if (description.kinetic_molal_index != 0.0)
96 185 : rate *= std::pow(kin_molality, description.kinetic_molal_index);
97 5948 : if (description.kinetic_monod_index != 0.0)
98 185 : rate /=
99 185 : std::pow(std::pow(kin_molality, description.kinetic_molal_index) +
100 185 : std::pow(description.kinetic_half_saturation, description.kinetic_molal_index),
101 : description.kinetic_monod_index);
102 :
103 : // promoting species numerators
104 52391 : for (unsigned i = 0; i < num_basis; ++i)
105 : {
106 46443 : if (promoting_indices[i] == 0.0)
107 45301 : continue;
108 2805 : if (basis_species_gas[i] || basis_species_name[i] == "H+" || basis_species_name[i] == "OH-")
109 591 : rate *= std::pow(basis_activity[i], promoting_indices[i]);
110 : else
111 551 : rate *= std::pow(basis_molality[i], promoting_indices[i]);
112 : }
113 184871 : for (unsigned j = 0; j < num_eqm; ++j)
114 : {
115 178923 : const unsigned index = num_basis + j;
116 178923 : if (promoting_indices[index] == 0.0)
117 175645 : continue;
118 7078 : if (eqm_species_gas[j] || eqm_species_name[j] == "H+" || eqm_species_name[j] == "OH-")
119 2756 : rate *= std::pow(eqm_activity[j], promoting_indices[index]);
120 : else
121 522 : rate *= std::pow(eqm_molality[j], promoting_indices[index]);
122 : }
123 :
124 : // promoting species denominators: the monod terms
125 52391 : for (unsigned i = 0; i < num_basis; ++i)
126 : {
127 46443 : if (promoting_monod_indices[i] == 0.0)
128 46211 : continue;
129 618 : if (basis_species_gas[i] || basis_species_name[i] == "H+" || basis_species_name[i] == "OH-")
130 59 : rate /= std::pow(std::pow(basis_activity[i], promoting_indices[i]) +
131 59 : std::pow(promoting_half_saturation[i], promoting_indices[i]),
132 : promoting_monod_indices[i]);
133 : else
134 173 : rate /= std::pow(std::pow(basis_molality[i], promoting_indices[i]) +
135 173 : std::pow(promoting_half_saturation[i], promoting_indices[i]),
136 : promoting_monod_indices[i]);
137 : }
138 184871 : for (unsigned j = 0; j < num_eqm; ++j)
139 : {
140 178923 : const unsigned index = num_basis + j;
141 178923 : if (promoting_monod_indices[index] == 0.0)
142 178390 : continue;
143 1563 : if (eqm_species_gas[j] || eqm_species_name[j] == "H+" || eqm_species_name[j] == "OH-")
144 36 : rate /= std::pow(std::pow(eqm_activity[j], promoting_indices[index]) +
145 36 : std::pow(promoting_half_saturation[index], promoting_indices[index]),
146 : promoting_monod_indices[index]);
147 : else
148 497 : rate /= std::pow(std::pow(eqm_molality[j], promoting_indices[index]) +
149 497 : std::pow(promoting_half_saturation[index], promoting_indices[index]),
150 : promoting_monod_indices[index]);
151 : }
152 :
153 : // temperature dependence
154 5948 : rate *= std::exp(
155 5948 : description.activation_energy / GeochemistryConstants::GAS_CONSTANT *
156 5948 : (description.one_over_T0 - 1.0 / (temp_degC + GeochemistryConstants::CELSIUS_TO_KELVIN)));
157 :
158 : // dependence on activity and equilibrium constant
159 5948 : const Real log10K_bio = log10K - description.energy_captured /
160 5948 : GeochemistryConstants::GAS_CONSTANT /
161 5948 : (temp_degC + GeochemistryConstants::CELSIUS_TO_KELVIN) /
162 : GeochemistryConstants::LOGTEN;
163 5948 : const Real ap_over_k = std::pow(10.0, log10_activity_product - log10K_bio);
164 5948 : const Real theta_term = std::pow(ap_over_k, description.theta);
165 5948 : switch (description.direction)
166 : {
167 5139 : case DirectionChoiceEnum::BOTH:
168 5139 : if (ap_over_k > 1.0) // precipitation
169 645 : rate = -rate;
170 : // leave the sign of rate unchanged for dissolution
171 : break;
172 572 : case DirectionChoiceEnum::DISSOLUTION:
173 572 : if (ap_over_k > 1.0) // precipitation
174 1 : rate = 0.0;
175 : // leave the sign of rate unchanged for dissolution
176 : break;
177 2 : case DirectionChoiceEnum::PRECIPITATION:
178 2 : if (ap_over_k > 1.0) // precipitation
179 1 : rate = -rate;
180 : else // dissolution
181 1 : rate = 0.0;
182 : break;
183 : case DirectionChoiceEnum::RAW:
184 : break; // no dependence on 1 - ap_over_k
185 : case DirectionChoiceEnum::DEATH:
186 : break; // no dependence on 1 - ap_over_k
187 : }
188 : // const Real rate_no_theta_term = rate; // needed for derivative calcs when theta_term == 1
189 11896 : rate *= (description.eta == 0.0)
190 5948 : ? 1.0
191 5560 : : ((theta_term == 1.0) ? 0.0 : std::pow(std::abs(1.0 - theta_term), description.eta));
192 :
193 52391 : for (unsigned i = 0; i < num_basis; ++i)
194 46443 : drate_dmol[i] = 0.0;
195 5948 : drate_dkin = 0.0;
196 :
197 : // derivatives with respect to the kinetic species moles
198 5948 : if (description.multiply_by_mass)
199 5600 : drate_dkin += rate / kin_moles;
200 :
201 5948 : if (description.kinetic_molal_index != 0.0)
202 : {
203 185 : const Real d_by_dkin_molality = description.kinetic_molal_index * rate / kin_molality;
204 185 : drate_dkin += d_by_dkin_molality * dkin_molality_dkin_moles;
205 185 : drate_dmol[0] += d_by_dkin_molality * dkin_molality_dnw;
206 : }
207 5948 : if (!(description.kinetic_molal_index == 0 || description.kinetic_monod_index == 0.0))
208 : {
209 : const Real d_by_dkin_molality =
210 185 : -description.kinetic_monod_index * description.kinetic_molal_index *
211 185 : std::pow(kin_molality, description.kinetic_molal_index - 1) * rate /
212 185 : (std::pow(kin_molality, description.kinetic_molal_index) +
213 185 : std::pow(description.kinetic_half_saturation, description.kinetic_molal_index));
214 185 : drate_dkin += d_by_dkin_molality * dkin_molality_dkin_moles;
215 185 : drate_dmol[0] += d_by_dkin_molality * dkin_molality_dnw;
216 : }
217 :
218 : // Derivatives of the promoting-species numerators (ignore all derivatives of activity
219 : // coefficients), so
220 : // d(molality^P) / dmolality = P * molality^P / molality , and
221 : // d(activity^P)/d(molality) = P activity^(P-1) d(activity)/d(molality) = P activity^(P-1)
222 : // activity_product = P * activity^P / molality
223 52391 : for (unsigned i = 0; i < num_basis; ++i)
224 : {
225 46443 : if (promoting_indices[i] == 0.0)
226 45301 : continue;
227 1142 : else if (basis_species_gas[i]) // molality is undefined
228 30 : continue;
229 : else
230 1112 : drate_dmol[i] += promoting_indices[i] * rate / basis_molality[i];
231 : }
232 184871 : for (unsigned j = 0; j < num_eqm; ++j)
233 : {
234 : // d(eqm^P)/d(basis_i) = P eqm^P / eqm * d(eqm)/d(basis_i) = P eqm^P * stoi / basis_i
235 178923 : const unsigned index = num_basis + j;
236 178923 : if (promoting_indices[index] == 0.0)
237 175645 : continue;
238 31568 : for (unsigned i = 1; i < num_basis; ++i) // deriv of water activity is ignored
239 28290 : if (!(basis_species_gas[i] || eqm_stoichiometry(j, i) == 0.0))
240 19587 : drate_dmol[i] +=
241 19587 : promoting_indices[index] * rate * eqm_stoichiometry(j, i) / basis_molality[i];
242 : }
243 : // Derivatives of the promoting-species denominators (ignore all derivatives of activity
244 : // coefficients) so
245 : // d((molality^P + K)^(-u))/d(molality) = -u (molality^P + K)^(-u) / (molality^P + K) * P *
246 : // molality^P / molality
247 : // d((activity^P + K)^(-u))/d(molality) = -u (activity^P + K)^(-u) /
248 : // (activity^P + K) * P * activity^P / molality
249 52391 : for (unsigned i = 0; i < num_basis; ++i)
250 : {
251 46443 : if (promoting_indices[i] == 0.0 || promoting_monod_indices[i] == 0.0)
252 46211 : continue;
253 232 : else if (basis_species_gas[i]) // molality is undefined
254 19 : continue;
255 386 : else if (basis_species_name[i] == "H+" || basis_species_name[i] == "OH-")
256 40 : drate_dmol[i] -= promoting_monod_indices[i] * promoting_indices[i] *
257 40 : std::pow(basis_activity[i], promoting_indices[i]) * rate /
258 40 : basis_molality[i] /
259 40 : (std::pow(basis_activity[i], promoting_indices[i]) +
260 40 : std::pow(promoting_half_saturation[i], promoting_indices[i]));
261 : else
262 173 : drate_dmol[i] -= promoting_monod_indices[i] * promoting_indices[i] *
263 173 : std::pow(basis_molality[i], promoting_indices[i] - 1) * rate /
264 173 : (std::pow(basis_molality[i], promoting_indices[i]) +
265 173 : std::pow(promoting_half_saturation[i], promoting_indices[i]));
266 : }
267 184871 : for (unsigned j = 0; j < num_eqm; ++j)
268 : {
269 178923 : const unsigned index = num_basis + j;
270 178923 : if (promoting_indices[index] == 0.0 || promoting_monod_indices[index] == 0.0)
271 178390 : continue;
272 4293 : for (unsigned i = 1; i < num_basis; ++i) // deriv of water activity is ignored
273 : {
274 3760 : if (basis_species_gas[i] || eqm_stoichiometry(j, i) == 0.0)
275 2892 : continue;
276 2568 : else if (eqm_species_gas[j] || eqm_species_name[j] == "H+" || eqm_species_name[j] == "OH-")
277 36 : drate_dmol[i] -= promoting_monod_indices[index] * promoting_indices[index] *
278 36 : std::pow(eqm_activity[j], promoting_indices[index]) * rate *
279 36 : eqm_stoichiometry(j, i) /
280 36 : (std::pow(eqm_activity[j], promoting_indices[index]) +
281 36 : std::pow(promoting_half_saturation[index], promoting_indices[index])) /
282 : basis_molality[i];
283 : else
284 832 : drate_dmol[i] -= promoting_monod_indices[index] * promoting_indices[index] *
285 832 : std::pow(eqm_molality[j], promoting_indices[index]) * rate *
286 832 : eqm_stoichiometry(j, i) /
287 832 : (std::pow(eqm_molality[j], promoting_indices[index]) +
288 832 : std::pow(promoting_half_saturation[index], promoting_indices[index])) /
289 : basis_molality[i];
290 : }
291 : }
292 : // derivative of the activity-product term
293 : Real deriv_ap_term = 0.0;
294 5948 : if (theta_term != 1.0)
295 5943 : deriv_ap_term =
296 5943 : description.eta * rate / std::abs(1 - theta_term) * (-description.theta) * theta_term;
297 : else // theta_term = 1, so rate = 0 (unless eta = 0 too)
298 : {
299 5 : if (description.eta > 1)
300 : deriv_ap_term = 0.0;
301 4 : else if (description.eta == 1.0)
302 : deriv_ap_term = 0.0; // rate_no_theta_term * description.theta * theta_term;
303 2 : else if (description.eta == 0.0)
304 : deriv_ap_term = 0.0;
305 : else
306 : deriv_ap_term = std::numeric_limits<Real>::max();
307 : }
308 5948 : if (theta_term > 1.0)
309 687 : deriv_ap_term = -deriv_ap_term;
310 46443 : for (unsigned i = 1; i < num_basis; ++i)
311 40495 : if (!(basis_species_gas[i] || kin_stoichiometry(kin, i) == 0.0 || deriv_ap_term == 0.0))
312 21420 : drate_dmol[i] += deriv_ap_term * kin_stoichiometry(kin, i) / basis_molality[i];
313 5948 : }
314 : }
|
