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 "GeochemistryConsoleOutput.h"
11 : #include "GeochemistryConstants.h"
12 : #include "GeochemistryFormattedOutput.h"
13 : #include "GeochemistrySortedIndices.h"
14 :
15 : registerMooseObject("GeochemistryApp", GeochemistryConsoleOutput);
16 :
17 : InputParameters
18 1469 : GeochemistryConsoleOutput::sharedParams()
19 : {
20 1469 : InputParameters params = emptyInputParameters();
21 2938 : params.addParam<unsigned int>("precision", 4, "Precision for printing values");
22 2938 : params.addParam<Real>(
23 : "mol_cutoff",
24 2938 : 1E-40,
25 : "Information regarding species with molalities less than this amount will not be outputted");
26 2938 : params.addParam<bool>(
27 : "solver_info",
28 2938 : false,
29 : "Print information (to the console) from the solver including residuals, swaps, etc");
30 1469 : return params;
31 0 : }
32 :
33 : InputParameters
34 846 : GeochemistryConsoleOutput::validParams()
35 : {
36 846 : InputParameters params = Output::validParams();
37 846 : params += GeochemistryConsoleOutput::sharedParams();
38 1692 : params.addRequiredParam<UserObjectName>("geochemistry_reactor",
39 : "The name of the GeochemistryReactor UserObject");
40 2538 : params.addRangeCheckedParam<Real>("stoichiometry_tolerance",
41 1692 : 1E-6,
42 : "stoichiometry_tolerance >= 0.0",
43 : "if abs(any stoichiometric coefficient) < stoi_tol then it is "
44 : "set to zero, and so will not appear in the output");
45 1692 : params.addRequiredParam<UserObjectName>(
46 : "nearest_node_number_UO",
47 : "The NearestNodeNumber UserObject that defines the physical point at which to query the "
48 : "GeochemistryReactor");
49 846 : params.addClassDescription("Outputs results from a GeochemistryReactor at a particular point");
50 846 : return params;
51 0 : }
52 :
53 423 : GeochemistryConsoleOutput::GeochemistryConsoleOutput(const InputParameters & parameters)
54 : : Output(parameters),
55 : UserObjectInterface(this),
56 423 : _reactor(getUserObject<GeochemistryReactorBase>("geochemistry_reactor")),
57 423 : _nnn(getUserObject<NearestNodeNumberUO>("nearest_node_number_UO")),
58 846 : _precision(getParam<unsigned int>("precision")),
59 846 : _stoi_tol(getParam<Real>("stoichiometry_tolerance")),
60 846 : _solver_info(getParam<bool>("solver_info")),
61 1269 : _mol_cutoff(getParam<Real>("mol_cutoff"))
62 : {
63 423 : }
64 :
65 : void
66 506 : GeochemistryConsoleOutput::output()
67 : {
68 506 : const Node * closest_node = _nnn.getClosestNode();
69 506 : if (!closest_node)
70 192 : return;
71 : const dof_id_type closest_id = closest_node->id();
72 :
73 314 : if (_solver_info)
74 435 : _console << _reactor.getSolverOutput(closest_id).str();
75 :
76 : // retrieve information
77 314 : const GeochemicalSystem & egs = _reactor.getGeochemicalSystem(closest_id);
78 314 : const unsigned num_basis = egs.getNumInBasis();
79 314 : const unsigned num_eqm = egs.getNumInEquilibrium();
80 314 : const unsigned num_kin = egs.getNumKinetic();
81 314 : const std::vector<Real> & basis_molality = egs.getSolventMassAndFreeMolalityAndMineralMoles();
82 314 : const std::vector<Real> & basis_activity = egs.getBasisActivity();
83 314 : const std::vector<Real> & basis_act_coef = egs.getBasisActivityCoefficient();
84 314 : const std::vector<Real> & bulk_moles = egs.getBulkMolesOld();
85 314 : const std::vector<Real> & eqm_molality = egs.getEquilibriumMolality();
86 314 : const std::vector<Real> & eqm_act_coef = egs.getEquilibriumActivityCoefficient();
87 314 : const std::vector<Real> & eqm_SI = egs.getSaturationIndices();
88 314 : const std::vector<Real> & kin_moles = egs.getKineticMoles();
89 314 : const ModelGeochemicalDatabase & mgd = egs.getModelGeochemicalDatabase();
90 :
91 314 : _console << std::setprecision(_precision);
92 :
93 314 : _console << "\nSummary:\n";
94 :
95 314 : _console << "Total number of iterations required = " << _reactor.getSolverIterations(closest_id)
96 314 : << "\n";
97 314 : _console << "Error in calculation = " << _reactor.getSolverResidual(closest_id) << "mol\n";
98 314 : _console << "Charge of solution = " << egs.getTotalChargeOld() << "mol";
99 314 : _console << " (charge-balance species = "
100 314 : << mgd.basis_species_name[egs.getChargeBalanceBasisIndex()] << ")\n";
101 :
102 314 : _console << "Mass of solvent water = " << basis_molality[0] << "kg\n";
103 :
104 314 : Real mass = bulk_moles[0] / GeochemistryConstants::MOLES_PER_KG_WATER;
105 1757 : for (unsigned i = 1; i < num_basis; ++i) // do not loop over water
106 1443 : mass += bulk_moles[i] * mgd.basis_species_molecular_weight[i] / 1000.0;
107 314 : _console << "Total mass = " << mass << "kg";
108 314 : if (num_kin == 0)
109 264 : _console << "\n";
110 : else
111 : {
112 50 : _console << " (including kinetic species and free minerals)\n";
113 100 : for (unsigned k = 0; k < num_kin; ++k)
114 50 : mass -= kin_moles[k] * mgd.kin_species_molecular_weight[k] / 1000.0;
115 50 : _console << "Mass without kinetic species but including free minerals = " << mass << "kg\n";
116 : }
117 : // remove the free minerals
118 1757 : for (unsigned i = 1; i < num_basis; ++i) // do not loop over water
119 1443 : if (mgd.basis_species_mineral[i])
120 90 : mass -= basis_molality[i] * mgd.basis_species_molecular_weight[i] / 1000.0;
121 : // remove surface complexes
122 314 : for (const auto & name_info :
123 369 : mgd.surface_complexation_info) // all minerals involved in surface complexation
124 55 : for (const auto & name_frac :
125 165 : name_info.second.sorption_sites) // all sorption sites on the given mineral
126 : {
127 : const unsigned i =
128 110 : mgd.basis_species_index.at(name_frac.first); // i = basis_index_of_sorption_site
129 110 : mass -= basis_molality[i] * mgd.basis_species_molecular_weight[i] / 1000.0;
130 : }
131 314 : _console << "Mass of aqueous solution = " << mass << "kg";
132 314 : if (num_kin == 0)
133 264 : _console << " (without free minerals)\n";
134 : else
135 50 : _console << " (without kinetic species and without free minerals)\n";
136 :
137 : // Output the aqueous solution pH, if relevant
138 628 : if (mgd.basis_species_index.count("H+"))
139 822 : _console << "pH = " << -std::log10(basis_activity[mgd.basis_species_index.at("H+")]) << "\n";
140 628 : if (mgd.eqm_species_index.count("H+"))
141 30 : _console << "pH = "
142 120 : << -std::log10(eqm_molality[mgd.eqm_species_index.at("H+")] *
143 90 : eqm_act_coef[mgd.eqm_species_index.at("H+")])
144 30 : << "\n";
145 :
146 : // Output the aqueous solution pe, if relevant
147 314 : if (mgd.redox_stoichiometry.m() > 0)
148 50 : _console << "pe = " << egs.getRedoxLog10K(0) - egs.log10RedoxActivityProduct(0) << "\n";
149 :
150 : // Output ionic strengths
151 314 : _console << "Ionic strength = " << egs.getIonicStrength() << "mol/kg(solvent water)\n";
152 314 : _console << "Stoichiometric ionic strength = " << egs.getStoichiometricIonicStrength()
153 314 : << "mol/kg(solvent water)\n";
154 :
155 : // Output activity of water
156 314 : _console << "Activity of water = " << basis_activity[0] << "\n";
157 :
158 : // Output temperature
159 314 : _console << "Temperature = " << egs.getTemperature() << "\n";
160 :
161 : // Output the basis species information, sorted by molality
162 : std::vector<unsigned> basis_order =
163 314 : GeochemistrySortedIndices::sortedIndices(basis_molality, false);
164 314 : _console << "\nBasis Species:\n";
165 2071 : for (const auto & i : basis_order)
166 1757 : if (i == 0 || mgd.basis_species_gas[i])
167 329 : continue;
168 : else
169 : {
170 1428 : _console << mgd.basis_species_name[i] << "; bulk_moles = " << bulk_moles[i]
171 1428 : << "mol; bulk_conc = "
172 1428 : << bulk_moles[i] * mgd.basis_species_molecular_weight[i] * 1000.0 / mass
173 1428 : << "mg/kg(soln);";
174 1428 : if (!mgd.basis_species_mineral[i])
175 1338 : _console << " molality = " << basis_molality[i] << "mol/kg(solvent water); free_conc = "
176 1338 : << basis_molality[i] * basis_molality[0] / mass *
177 1338 : mgd.basis_species_molecular_weight[i] * 1000.0
178 1338 : << "mg/kg(soln); act_coeff = " << basis_act_coef[i]
179 1338 : << "; log10(a) = " << std::log10(basis_activity[i]) << "\n";
180 90 : else if (mgd.basis_species_mineral[i])
181 90 : _console << " free_moles = " << basis_molality[i] << "mol; free_mass = "
182 90 : << basis_molality[i] * mgd.basis_species_molecular_weight[i] * 1000.0 << "mg\n";
183 : }
184 2071 : for (unsigned i = 0; i < num_basis; ++i)
185 1757 : if (mgd.basis_species_gas[i])
186 15 : _console << mgd.basis_species_name[i] << "; fugacity = " << basis_activity[i] << "\n";
187 :
188 : // Output the equilibrium species info, sorted by molality
189 314 : std::vector<unsigned> eqm_order = GeochemistrySortedIndices::sortedIndices(eqm_molality, false);
190 314 : _console << "\nEquilibrium Species:\n";
191 1922 : for (const auto & i : eqm_order)
192 1678 : if (eqm_molality[i] <= _mol_cutoff)
193 : break;
194 1608 : else if (mgd.eqm_species_gas[i])
195 0 : continue;
196 : else
197 1608 : _console << mgd.eqm_species_name[i] << "; molality = " << eqm_molality[i]
198 1608 : << "mol/kg(solvent water); free_conc = "
199 1608 : << eqm_molality[i] * basis_molality[0] / mass * mgd.eqm_species_molecular_weight[i] *
200 1608 : 1000.0
201 1608 : << "mg/kg(soln); act_coeff = " << eqm_act_coef[i]
202 1608 : << "; log10(a) = " << std::log10(eqm_molality[i] * eqm_act_coef[i]) << "; "
203 1608 : << mgd.eqm_species_name[i] << " = "
204 1608 : << GeochemistryFormattedOutput::reaction(
205 1608 : mgd.eqm_stoichiometry, i, mgd.basis_species_name, _stoi_tol, _precision)
206 1608 : << "; log10K = " << egs.getLog10K(i) << "\n";
207 6472 : for (unsigned i = 0; i < num_eqm; ++i)
208 6158 : if (mgd.eqm_species_gas[i])
209 : {
210 : Real log10f = 0;
211 0 : for (unsigned basis_i = 0; basis_i < num_basis; ++basis_i)
212 0 : log10f += mgd.eqm_stoichiometry(i, basis_i) * std::log10(basis_activity[basis_i]);
213 0 : log10f -= egs.getLog10K(i);
214 0 : _console << mgd.eqm_species_name[i]
215 0 : << "; act_coeff = " << egs.getEquilibriumActivityCoefficient(i)
216 0 : << "; fugacity = " << std::pow(10.0, log10f) << "; " << mgd.eqm_species_name[i]
217 0 : << " = "
218 0 : << GeochemistryFormattedOutput::reaction(
219 0 : mgd.eqm_stoichiometry, i, mgd.basis_species_name, _stoi_tol, _precision)
220 0 : << "; log10K = " << egs.getLog10K(i) << "\n";
221 : }
222 :
223 : // Output the kinetic species information, sorted by mole number
224 314 : std::vector<unsigned> kin_order = GeochemistrySortedIndices::sortedIndices(kin_moles, false);
225 314 : _console << "\nKinetic Species:\n";
226 364 : for (const auto & k : kin_order)
227 : {
228 50 : _console << mgd.kin_species_name[k] << "; moles = " << kin_moles[k]
229 50 : << "; mass = " << kin_moles[k] * mgd.kin_species_molecular_weight[k] * 1000.0
230 50 : << "mg; ";
231 50 : if (mgd.kin_species_mineral[k])
232 50 : _console << "volume = " << kin_moles[k] * mgd.kin_species_molecular_volume[k] << "cm^3; ";
233 50 : _console << mgd.kin_species_name[k] << " = "
234 50 : << GeochemistryFormattedOutput::reaction(
235 50 : mgd.kin_stoichiometry, k, mgd.basis_species_name, _stoi_tol, _precision)
236 50 : << "; log10(Q) = " << egs.log10KineticActivityProduct(k)
237 50 : << "; log10K = " << egs.getKineticLog10K(k)
238 100 : << "; dissolution_rate*dt = " << -_reactor.getMoleAdditions(closest_id)(num_basis + k)
239 50 : << "\n";
240 : }
241 :
242 : // Output the mineral info, sorted by saturation indices
243 314 : std::vector<unsigned> mineral_order = GeochemistrySortedIndices::sortedIndices(eqm_SI, false);
244 314 : _console << "\nMinerals:\n";
245 6472 : for (const auto & i : mineral_order)
246 6158 : if (mgd.eqm_species_mineral[i])
247 325 : _console << mgd.eqm_species_name[i] << " = "
248 325 : << GeochemistryFormattedOutput::reaction(
249 325 : mgd.eqm_stoichiometry, i, mgd.basis_species_name, _stoi_tol, _precision)
250 325 : << "; log10K = " << egs.getLog10K(i) << "; SI = " << eqm_SI[i] << "\n";
251 :
252 : // Output the Nernst potentials, if relevant
253 314 : _console << "\nNernst potentials:\n";
254 314 : if (mgd.redox_stoichiometry.m() > 0)
255 50 : outputNernstInfo(egs);
256 :
257 314 : const unsigned num_pot = egs.getNumSurfacePotentials();
258 314 : if (num_pot > 0)
259 : {
260 55 : _console << "\nSorbing surfaces:\n";
261 55 : const std::vector<Real> area = egs.getSorbingSurfaceArea();
262 110 : for (unsigned sp = 0; sp < num_pot; ++sp)
263 55 : _console << mgd.surface_sorption_name[sp] << "; area = " << area[sp]
264 55 : << "m^2; specific charge = " << egs.getSurfaceCharge(sp)
265 55 : << "C/m^2; surface potential = " << egs.getSurfacePotential(sp) << "V\n";
266 : }
267 :
268 314 : DenseVector<Real> bulk_in_original_basis = egs.getBulkOldInOriginalBasis();
269 314 : DenseVector<Real> transported_bulk_in_original_basis = egs.getTransportedBulkInOriginalBasis();
270 : std::vector<std::string> original_basis_names =
271 314 : _reactor.getPertinentGeochemicalSystem().originalBasisNames();
272 314 : _console << "\nIn original basis:\n";
273 2071 : for (unsigned i = 0; i < num_basis; ++i)
274 1757 : _console << original_basis_names[i] << "; total_bulk_moles = " << bulk_in_original_basis(i)
275 1757 : << "; transported_bulk_moles = " << transported_bulk_in_original_basis(i) << "\n";
276 :
277 314 : _console << std::flush;
278 314 : }
279 :
280 : void
281 50 : GeochemistryConsoleOutput::outputNernstInfo(const GeochemicalSystem & egs_ref) const
282 : {
283 : // Copy egs_ref so we can call non-const methods (viz, swaps). This does not copy the data in
284 : // egs.getModelGeochemicalDatabase(), only the reference (both references refer to the same
285 : // block of memory) and unfortunately the swaps below manipulate that memory
286 50 : GeochemicalSystem egs = egs_ref;
287 : // Since we only want to do the swaps to print out some Nernst info, but don't want to impact
288 : // the rest of the simulation, copy the mgd, so that we can copy it back into the aforementioned
289 : // block of memory at the end of this method
290 50 : ModelGeochemicalDatabase mgd_without_nernst_swaps = egs.getModelGeochemicalDatabase();
291 :
292 50 : const ModelGeochemicalDatabase & mgd_ref = egs.getModelGeochemicalDatabase();
293 : // attempt to undo the swaps that have been done to mgd_ref.
294 : // NOTE FOR FUTURE: In the current code (2020, June 1) swapping the redox stuff in mgd seems
295 : // redundant. While the current code does these swaps, here we just undo all the swaps again to
296 : // write the redox stuff in the original basis. Since this is the only place that the redox
297 : // stuff is used, doing any swaps on the redox stuff is currently just a waste of time! take a
298 : // copy of the following because they get modified during the swaps
299 50 : const std::vector<unsigned> have_swapped_out_of_basis = mgd_ref.have_swapped_out_of_basis;
300 50 : const std::vector<unsigned> have_swapped_into_basis = mgd_ref.have_swapped_into_basis;
301 150 : for (int sw = have_swapped_out_of_basis.size() - 1; sw >= 0; --sw)
302 : {
303 : // Don't check for gases being swapped in or out of the basis (which usually can't happen in
304 : // the middle of a Newton process) because we're going to trash all the swaps at the end of
305 : // this method, and we don't have to worry about bulk moles, etc: we just want the
306 : // stoichiometries and the activities
307 : try
308 : {
309 100 : egs.performSwapNoCheck(have_swapped_out_of_basis[sw], have_swapped_into_basis[sw]);
310 : }
311 0 : catch (const MooseException & e)
312 : {
313 0 : const std::string to_swap_in = mgd_ref.eqm_species_name[have_swapped_into_basis[sw]];
314 0 : const std::string to_swap_out = mgd_ref.basis_species_name[have_swapped_out_of_basis[sw]];
315 : // Don't crash the entire simulation, just because Nernst-related swapping does not work
316 0 : mooseWarning("Swapping ", to_swap_out, " and ", to_swap_in, ": ", e.what());
317 0 : }
318 : }
319 :
320 50 : const Real prefactor = -GeochemistryConstants::LOGTEN * GeochemistryConstants::GAS_CONSTANT *
321 50 : (egs.getTemperature() + GeochemistryConstants::CELSIUS_TO_KELVIN) /
322 50 : GeochemistryConstants::FARADAY;
323 :
324 50 : if (mgd_ref.redox_lhs == egs_ref.getOriginalRedoxLHS())
325 : {
326 50 : std::vector<Real> eh(mgd_ref.redox_stoichiometry.m());
327 140 : for (unsigned red = 0; red < mgd_ref.redox_stoichiometry.m(); ++red)
328 90 : eh[red] = prefactor * (egs.log10RedoxActivityProduct(red) - egs.getRedoxLog10K(red));
329 50 : const std::vector<unsigned> eh_order = GeochemistrySortedIndices::sortedIndices(eh, false);
330 140 : for (const auto & red : eh_order)
331 90 : _console << mgd_ref.redox_lhs << " = "
332 90 : << GeochemistryFormattedOutput::reaction(mgd_ref.redox_stoichiometry,
333 : red,
334 90 : mgd_ref.basis_species_name,
335 90 : _stoi_tol,
336 90 : _precision)
337 90 : << "; Eh = " << eh[red] << "V\n";
338 : }
339 :
340 50 : _console << std::flush;
341 :
342 : // restore the original mgd by copying the data in copy_of_mgd into the memory referenced by
343 : // egs.getModelGeochemicalDatabase()
344 50 : egs.setModelGeochemicalDatabase(mgd_without_nernst_swaps);
345 50 : }
|
