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 "GeochemistryTimeDependentReactor.h"
11 :
12 : registerMooseObject("GeochemistryApp", GeochemistryTimeDependentReactor);
13 :
14 : InputParameters
15 827 : GeochemistryTimeDependentReactor::sharedParams()
16 : {
17 827 : InputParameters params = emptyInputParameters();
18 1654 : params.addParam<unsigned>(
19 : "ramp_max_ionic_strength_subsequent",
20 1654 : 0,
21 : "The number of iterations over which to progressively increase the maximum ionic strength "
22 : "(from zero to max_ionic_strength) during time-stepping. Unless a great deal occurs in each "
23 : "time step, this parameter can be set quite small");
24 1654 : params.addCoupledVar(
25 : "mode",
26 : 0.0,
27 : "This may vary temporally. If mode=1 then 'dump' mode is used, which means all non-kinetic "
28 : "mineral masses are removed from the system before the equilibrium solution is sought (ie, "
29 : "removal occurs at the beginning of the time step). If mode=2 then 'flow-through' mode is "
30 : "used, which means all mineral masses are removed from the system after it the "
31 : "equilbrium solution has been found (ie, at the end of a time step). If mode=3 then 'flush' "
32 : "mode is used, then before the equilibrium solution is sought (ie, at the start of a time "
33 : "step) water+species is removed from the system at the same rate as pure water + non-mineral "
34 : "solutes are entering the system (specified in source_species_rates). If mode=4 then "
35 : "'heat-exchanger' mode is used, which means the entire current aqueous solution is removed, "
36 : "then the source_species are added, then the temperature is set to 'cold_temperature', the "
37 : "system is solved and any precipitated minerals are removed, then the temperature is set to "
38 : "'temperature', the system re-solved and any precipitated minerals are removed. If mode is "
39 : "any other number, no special mode is active (the system simply responds to the "
40 : "source_species_rates, controlled_activity_value, etc).");
41 1654 : params.addCoupledVar(
42 : "temperature",
43 : 25,
44 : "Temperature. This has two different meanings if mode!=4. (1) If no species are being "
45 : "added to the solution (no source_species_rates are positive) then this is the temperature "
46 : "of the aqueous solution. (2) If species are being added, this is the temperature of the "
47 : "species being added. In case (2), the final aqueous-solution temperature is computed "
48 : "assuming the species are added, temperature is equilibrated and then, if species are also "
49 : "being removed, they are removed. If you wish to add species and simultaneously alter the "
50 : "temperature, you will have to use a sequence of heat-add-heat-add, etc steps. In the case "
51 : "that mode=4, temperature is the final temperature of the aqueous solution");
52 1654 : params.addCoupledVar(
53 : "cold_temperature",
54 : 25,
55 : "This is only used if mode=4, where it is the cold temperature of the heat exchanger.");
56 2481 : params.addRangeCheckedParam<unsigned>(
57 : "heating_increments",
58 1654 : 1,
59 : "heating_increments > 0",
60 : "This is only used if mode=4. Internal to this object, the temperature is ramped from "
61 : "cold_temperature to temperature in heating_increments increments. This helps difficult "
62 : "problems converge");
63 1654 : params.addParam<Real>("initial_temperature",
64 1654 : 25.0,
65 : "The initial aqueous solution is equilibrated at this system before adding "
66 : "reactants, changing temperature, etc.");
67 1654 : params.addParam<Real>("close_system_at_time",
68 1654 : 0.0,
69 : "Time at which to 'close' the system, that is, change a kg_solvent_water "
70 : "constraint to moles_bulk_water, and all free_molality and "
71 : "free_moles_mineral_species to moles_bulk_species");
72 1654 : params.addParam<std::vector<std::string>>(
73 : "remove_fixed_activity_name",
74 : {},
75 : "The name of the species that should have their activity or fugacity constraint removed at "
76 : "time given in remove_fixed_activity_time. There should be an equal number of these names "
77 : "as times given in remove_fixed_activity_time. Each of these must be in the basis and have "
78 : "an activity or fugacity constraint");
79 1654 : params.addParam<std::vector<Real>>("remove_fixed_activity_time",
80 : {},
81 : "The times at which the species in remove_fixed_activity_name "
82 : "should have their activity or fugacity constraint removed.");
83 1654 : params.addParam<std::vector<std::string>>(
84 : "source_species_names",
85 : {},
86 : "The name of the species that are added at rates given in source_species_rates. There must "
87 : "be an equal number of these as source_species_rates.");
88 1654 : params.addCoupledVar("source_species_rates",
89 : "Rates, in mols/time_unit, of addition of the species with names given in "
90 : "source_species_names. A negative value corresponds to removing a species: "
91 : "be careful that you don't cause negative mass problems!");
92 1654 : params.addParam<std::vector<std::string>>(
93 : "controlled_activity_name",
94 : {},
95 : "The names of the species that have their activity or fugacity constrained. There should be "
96 : "an equal number of these names as values given in controlled_activity_value. NOTE: if "
97 : "these species are not in the basis, or they do not have an activity (or fugacity) "
98 : "constraint then their activity cannot be controlled: in this case MOOSE will ignore the "
99 : "value you prescribe in controlled_activity_value.");
100 1654 : params.addCoupledVar("controlled_activity_value",
101 : "Values of the activity or fugacity of the species in "
102 : "controlled_activity_name list. These should always be positive");
103 1654 : params.addParam<bool>(
104 : "evaluate_kinetic_rates_always",
105 1654 : true,
106 : "If true, then, evaluate the kinetic rates at every Newton step during the solve using the "
107 : "current values of molality, activity, etc (ie, implement an implicit solve). If false, "
108 : "then evaluate the kinetic rates using the values of molality, activity, etc, at the start "
109 : "of the current time step (ie, implement an explicit solve)");
110 1654 : params.addParam<std::vector<std::string>>(
111 : "kinetic_species_name",
112 : {},
113 : "Names of the kinetic species given initial values in kinetic_species_initial_value");
114 1654 : params.addParam<std::vector<Real>>(
115 : "kinetic_species_initial_value",
116 : {},
117 : "Initial number of moles, mass or volume (depending on kinetic_species_unit) for each of the "
118 : "species named in kinetic_species_name");
119 : MultiMooseEnum kin_species_unit("dimensionless moles molal kg g mg ug kg_per_kg_solvent "
120 827 : "g_per_kg_solvent mg_per_kg_solvent ug_per_kg_solvent cm3");
121 1654 : params.addParam<MultiMooseEnum>(
122 : "kinetic_species_unit",
123 : kin_species_unit,
124 : "Units of the numerical values given in kinetic_species_initial_value. Moles: mole number. "
125 : "kg: kilograms. g: grams. mg: milligrams. ug: micrograms. cm3: cubic centimeters");
126 827 : return params;
127 827 : }
128 :
129 : InputParameters
130 581 : GeochemistryTimeDependentReactor::validParams()
131 : {
132 581 : InputParameters params = GeochemistryReactorBase::validParams();
133 581 : params += GeochemistryTimeDependentReactor::sharedParams();
134 581 : params.addClassDescription("UserObject that controls the time-dependent geochemistry reaction "
135 : "processes. Spatial dependence is not possible using this class");
136 581 : return params;
137 0 : }
138 :
139 338 : GeochemistryTimeDependentReactor::GeochemistryTimeDependentReactor(
140 338 : const InputParameters & parameters)
141 : : GeochemistryReactorBase(parameters),
142 338 : _temperature(coupledValue("temperature")),
143 338 : _cold_temperature(coupledValue("cold_temperature")),
144 676 : _heating_increments(getParam<unsigned>("heating_increments")),
145 676 : _new_temperature(getParam<Real>("initial_temperature")),
146 676 : _previous_temperature(getParam<Real>("initial_temperature")),
147 4056 : _egs(_mgd,
148 : _gac,
149 338 : _is,
150 338 : _swapper,
151 : getParam<std::vector<std::string>>("swap_out_of_basis"),
152 : getParam<std::vector<std::string>>("swap_into_basis"),
153 676 : getParam<std::string>("charge_balance_species"),
154 : getParam<std::vector<std::string>>("constraint_species"),
155 : getParam<std::vector<Real>>("constraint_value"),
156 : getParam<MultiMooseEnum>("constraint_unit"),
157 : getParam<MultiMooseEnum>("constraint_meaning"),
158 : _previous_temperature,
159 676 : getParam<unsigned>("extra_iterations_to_make_consistent"),
160 676 : getParam<Real>("min_initial_molality"),
161 : getParam<std::vector<std::string>>("kinetic_species_name"),
162 : getParam<std::vector<Real>>("kinetic_species_initial_value"),
163 : getParam<MultiMooseEnum>("kinetic_species_unit")),
164 1014 : _solver(_mgd.basis_species_name.size(),
165 : _mgd.kin_species_name.size(),
166 : _is,
167 676 : getParam<Real>("abs_tol"),
168 676 : getParam<Real>("rel_tol"),
169 676 : getParam<unsigned>("max_iter"),
170 338 : getParam<Real>("max_initial_residual"),
171 338 : _small_molality,
172 338 : _max_swaps_allowed,
173 : getParam<std::vector<std::string>>("prevent_precipitation"),
174 676 : getParam<Real>("max_ionic_strength"),
175 338 : getParam<unsigned>("ramp_max_ionic_strength_initial"),
176 676 : getParam<bool>("evaluate_kinetic_rates_always")),
177 338 : _num_kin(_egs.getNumKinetic()),
178 676 : _close_system_at_time(getParam<Real>("close_system_at_time")),
179 338 : _closed_system(false),
180 1014 : _source_species_names(getParam<std::vector<std::string>>("source_species_names")),
181 338 : _num_source_species(_source_species_names.size()),
182 338 : _source_species_rates(coupledValues("source_species_rates")),
183 676 : _remove_fixed_activity_name(getParam<std::vector<std::string>>("remove_fixed_activity_name")),
184 1014 : _remove_fixed_activity_time(getParam<std::vector<Real>>("remove_fixed_activity_time")),
185 338 : _num_removed_fixed(_remove_fixed_activity_name.size()),
186 338 : _removed_fixed_activity(_num_removed_fixed, false),
187 1014 : _controlled_activity_species_names(
188 : getParam<std::vector<std::string>>("controlled_activity_name")),
189 338 : _num_controlled_activity(_controlled_activity_species_names.size()),
190 338 : _controlled_activity_species_values(coupledValues("controlled_activity_value")),
191 338 : _mole_additions(_num_basis + _num_kin),
192 338 : _dmole_additions(_num_basis + _num_kin, _num_basis + _num_kin),
193 338 : _mode(coupledValue("mode")),
194 338 : _minerals_dumped(),
195 1014 : _ramp_subsequent(getParam<unsigned>("ramp_max_ionic_strength_subsequent"))
196 : {
197 : // check sources and set the rates
198 338 : if (coupledComponents("source_species_rates") != _num_source_species)
199 2 : paramError("source_species_names", "must have the same size as source_species_rates");
200 492 : for (unsigned i = 0; i < _num_source_species; ++i)
201 158 : if (!(_mgd.basis_species_index.count(_source_species_names[i]) == 1 ||
202 : _mgd.eqm_species_index.count(_source_species_names[i]) == 1 ||
203 : _mgd.kin_species_index.count(_source_species_names[i]) == 1))
204 4 : paramError("source_species_names",
205 2 : "The name " + _source_species_names[i] +
206 : " does not appear in the basis, equilibrium or kinetic species list");
207 :
208 : // check and set activity controllers
209 444 : for (unsigned i = 0; i < _num_removed_fixed; ++i)
210 : {
211 114 : if (_mgd.basis_species_index.count(_remove_fixed_activity_name[i]) == 0)
212 2 : paramError(
213 : "remove_fixed_activity_name",
214 : "The species ",
215 : _remove_fixed_activity_name[i],
216 : " is not in the basis, so cannot have its activity or fugacity constraint removed");
217 : else
218 : {
219 112 : const unsigned basis_ind = _mgd.basis_species_index.at(_remove_fixed_activity_name[i]);
220 112 : const GeochemicalSystem::ConstraintMeaningEnum cm = _egs.getConstraintMeaning()[basis_ind];
221 112 : if (!(cm == GeochemicalSystem::ConstraintMeaningEnum::ACTIVITY ||
222 : cm == GeochemicalSystem::ConstraintMeaningEnum::FUGACITY))
223 2 : paramError("remove_fixed_activity_name",
224 : "The species ",
225 : _remove_fixed_activity_name[i],
226 : " is does not have an activity or fugacity constraint so cannot have such a "
227 : "constraint removed");
228 : }
229 : }
230 330 : if (_num_removed_fixed != _remove_fixed_activity_time.size())
231 2 : paramError("remove_fixed_activity_name",
232 : "must be of the same size as remove_fixed_activity_time");
233 328 : if (coupledComponents("controlled_activity_value") != _num_controlled_activity)
234 2 : paramError("controlled_activity_name", "must have the same size as controlled_activity_value");
235 :
236 : // record coupled variables
237 : const std::vector<MooseVariableFEBase *> & coupled_vars = getCoupledMooseVars();
238 422 : for (unsigned int i = 0; i < coupled_vars.size(); i++)
239 192 : addMooseVariableDependency(coupled_vars[i]);
240 :
241 : // setup minerals_dumped
242 2548 : for (unsigned i = 0; i < _mgd.basis_species_name.size(); ++i)
243 2222 : if (_mgd.basis_species_mineral[i])
244 180 : _minerals_dumped[_mgd.basis_species_name[i]] = 0.0;
245 7422 : for (unsigned j = 0; j < _mgd.eqm_species_name.size(); ++j)
246 7096 : if (_mgd.eqm_species_mineral[j])
247 336 : _minerals_dumped[_mgd.eqm_species_name[j]] = 0.0;
248 496 : for (unsigned k = 0; k < _mgd.kin_species_name.size(); ++k)
249 170 : if (_mgd.kin_species_mineral[k])
250 146 : _minerals_dumped[_mgd.kin_species_name[k]] = 0.0;
251 326 : }
252 :
253 : void
254 5543 : GeochemistryTimeDependentReactor::initialize()
255 : {
256 : GeochemistryReactorBase::initialize();
257 5543 : }
258 : void
259 4157 : GeochemistryTimeDependentReactor::finalize()
260 : {
261 : GeochemistryReactorBase::finalize();
262 4157 : }
263 :
264 : void
265 311 : GeochemistryTimeDependentReactor::initialSetup()
266 : {
267 : // solve the geochemical system with its initial composition and with dt=0 so no kinetic additions
268 311 : if (_num_my_nodes == 0)
269 : return; // rather peculiar case where user has used many processors
270 : _mole_additions.zero();
271 207 : _dmole_additions.zero();
272 207 : _solver.solveSystem(_egs,
273 : _solver_output[0],
274 : _tot_iter[0],
275 : _abs_residual[0],
276 : 0.0,
277 207 : _mole_additions,
278 : _dmole_additions);
279 : }
280 :
281 : void
282 5542 : GeochemistryTimeDependentReactor::execute()
283 : {
284 5542 : if (_current_node->id() != 0)
285 : return;
286 :
287 2771 : _solver.setRampMaxIonicStrength(_ramp_subsequent);
288 :
289 : _mole_additions.zero();
290 2771 : _dmole_additions.zero();
291 :
292 : // remove appropriate constraints
293 2771 : if (!_closed_system && _t >= _close_system_at_time)
294 : {
295 155 : _egs.closeSystem();
296 155 : _closed_system = true;
297 : }
298 3911 : for (unsigned i = 0; i < _num_removed_fixed; ++i)
299 : {
300 1140 : if (!_removed_fixed_activity[i] && _t >= _remove_fixed_activity_time[i])
301 : {
302 : if (_mgd.basis_species_index.count(_remove_fixed_activity_name[i]))
303 54 : _egs.changeConstraintToBulk(_mgd.basis_species_index.at(_remove_fixed_activity_name[i]));
304 : _removed_fixed_activity[i] = true;
305 : }
306 : }
307 :
308 : // control activity
309 3221 : for (unsigned ca = 0; ca < _num_controlled_activity; ++ca)
310 : {
311 450 : const std::vector<GeochemicalSystem::ConstraintMeaningEnum> & cm = _egs.getConstraintMeaning();
312 450 : if (_mgd.basis_species_index.count(_controlled_activity_species_names[ca]))
313 : {
314 : const unsigned basis_ind =
315 450 : _mgd.basis_species_index.at(_controlled_activity_species_names[ca]);
316 450 : if (cm[basis_ind] == GeochemicalSystem::ConstraintMeaningEnum::ACTIVITY ||
317 : cm[basis_ind] == GeochemicalSystem::ConstraintMeaningEnum::FUGACITY)
318 378 : _egs.setConstraintValue(basis_ind, (*_controlled_activity_species_values[ca])[0]);
319 : }
320 : }
321 :
322 : // compute moles added
323 4421 : for (unsigned i = 0; i < _num_source_species; ++i)
324 : {
325 1650 : const Real this_rate = (*_source_species_rates[i])[0];
326 : if (_mgd.basis_species_index.count(_source_species_names[i]))
327 : {
328 456 : const unsigned basis_ind = _mgd.basis_species_index.at(_source_species_names[i]);
329 456 : _mole_additions(basis_ind) += this_rate;
330 : }
331 : else if (_mgd.eqm_species_index.count(_source_species_names[i]))
332 : {
333 1194 : const unsigned eqm_j = _mgd.eqm_species_index.at(_source_species_names[i]);
334 11346 : for (unsigned basis_ind = 0; basis_ind < _num_basis; ++basis_ind)
335 10152 : _mole_additions(basis_ind) += _mgd.eqm_stoichiometry(eqm_j, basis_ind) * this_rate;
336 : }
337 : else
338 : {
339 0 : const unsigned kin_ind = _mgd.kin_species_index.at(_source_species_names[i]);
340 0 : _mole_additions(_num_basis + kin_ind) += this_rate;
341 : }
342 : }
343 23280 : for (unsigned i = 0; i < _num_basis + _num_kin; ++i)
344 20509 : _mole_additions(i) *= _dt;
345 :
346 : // activate special modes
347 2771 : if (_mode[0] == 1.0)
348 60 : preSolveDump();
349 2711 : else if (_mode[0] == 3.0)
350 120 : preSolveFlush();
351 2591 : else if (_mode[0] == 4.0)
352 : {
353 0 : removeCurrentSpecies();
354 0 : _new_temperature = _cold_temperature[0];
355 : }
356 : else // nothing special: simply compute the desired temperature
357 2591 : _new_temperature = newTemperature(_mole_additions);
358 :
359 : // set temperature, if necessary
360 2771 : if (_new_temperature != _previous_temperature)
361 : {
362 468 : _egs.setTemperature(_new_temperature);
363 468 : _egs.computeConsistentConfiguration();
364 : }
365 2771 : _previous_temperature = _new_temperature;
366 :
367 : // solve the geochemical system
368 2771 : _solver.solveSystem(_egs,
369 : _solver_output[0],
370 : _tot_iter[0],
371 : _abs_residual[0],
372 2771 : _dt,
373 2771 : _mole_additions,
374 : _dmole_additions);
375 :
376 : // activate special modes
377 2771 : if (_mode[0] == 2.0)
378 66 : postSolveFlowThrough();
379 2705 : else if (_mode[0] == 4.0)
380 0 : postSolveExchanger();
381 : }
382 :
383 : const GeochemicalSystem &
384 682626 : GeochemistryTimeDependentReactor::getGeochemicalSystem(dof_id_type /*node_id*/) const
385 : {
386 682626 : return _egs;
387 : }
388 :
389 : const DenseVector<Real> &
390 2168 : GeochemistryTimeDependentReactor::getMoleAdditions(dof_id_type /*node_id*/) const
391 : {
392 2168 : return _mole_additions;
393 : }
394 :
395 : const std::stringstream &
396 11 : GeochemistryTimeDependentReactor::getSolverOutput(dof_id_type /*node_id*/) const
397 : {
398 11 : return _solver_output[0];
399 : }
400 :
401 64 : unsigned GeochemistryTimeDependentReactor::getSolverIterations(dof_id_type /*node_id*/) const
402 : {
403 64 : return _tot_iter[0];
404 : }
405 :
406 64 : Real GeochemistryTimeDependentReactor::getSolverResidual(dof_id_type /*node_id*/) const
407 : {
408 64 : return _abs_residual[0];
409 : }
410 :
411 : Real
412 528 : GeochemistryTimeDependentReactor::getMolesDumped(dof_id_type /*node_id*/,
413 : const std::string & species) const
414 : {
415 : if (_minerals_dumped.count(species) == 1)
416 528 : return _minerals_dumped.at(species);
417 : return 0.0;
418 : }
419 :
420 : Real
421 2771 : GeochemistryTimeDependentReactor::newTemperature(const DenseVector<Real> & mole_additions) const
422 : {
423 2771 : if (_temperature[0] == _previous_temperature)
424 : return _temperature[0];
425 :
426 : // if no reactants are being added, the system temperature will be _temperature[0]
427 : bool any_additions = false;
428 3144 : for (unsigned i = 0; i < _num_basis + _num_kin; ++i)
429 2676 : if (mole_additions(i) > 0)
430 : {
431 : any_additions = true;
432 : break;
433 : }
434 468 : if (!any_additions)
435 : return _temperature[0];
436 :
437 : // assume heat capacities of inputs and outputs are the same, so final temperature is dictated
438 : // by masses, also assume that the input happens first, then temperature equilibration, then
439 : // the outputs occur
440 : Real new_temperature = _temperature[0];
441 0 : const std::vector<Real> & current_bulk = _egs.getBulkMolesOld();
442 0 : Real current_kg = current_bulk[0] / GeochemistryConstants::MOLES_PER_KG_WATER;
443 0 : Real input_kg = std::max(mole_additions(0), 0.0) / GeochemistryConstants::MOLES_PER_KG_WATER;
444 0 : for (unsigned i = 1; i < _num_basis; ++i)
445 : {
446 0 : current_kg += current_bulk[i] * _mgd.basis_species_molecular_weight[i] / 1000.0;
447 0 : input_kg += std::max(mole_additions(i), 0.0) * _mgd.basis_species_molecular_weight[i] / 1000.0;
448 : }
449 0 : for (unsigned k = 0; k < _num_kin; ++k)
450 0 : input_kg += std::max(mole_additions(k + _num_basis), 0.0) *
451 0 : _mgd.kin_species_molecular_weight[k] / 1000.0;
452 0 : new_temperature =
453 0 : (_previous_temperature * current_kg + _temperature[0] * input_kg) / (current_kg + input_kg);
454 0 : return new_temperature;
455 : }
456 :
457 : void
458 60 : GeochemistryTimeDependentReactor::preSolveDump()
459 : {
460 : // remove basis mineral moles
461 60 : const std::vector<Real> & current_molal = _egs.getSolventMassAndFreeMolalityAndMineralMoles();
462 360 : for (unsigned i = 1; i < _num_basis; ++i)
463 300 : if (_mgd.basis_species_mineral[i])
464 : {
465 6 : _mole_additions(i) = -current_molal[i]; // might overwrite the rates set above, which is good
466 6 : _minerals_dumped[_mgd.basis_species_name[i]] += current_molal[i];
467 : }
468 :
469 60 : _new_temperature = newTemperature(_mole_additions);
470 :
471 : // add the chemicals immediately instead of during the solve (as occurs for other modes)
472 420 : for (unsigned basis_ind = 0; basis_ind < _num_basis; ++basis_ind)
473 : {
474 360 : _egs.addToBulkMoles(basis_ind, _mole_additions(basis_ind));
475 360 : _mole_additions(basis_ind) = 0.0;
476 : }
477 : // dump needs free mineral moles to be exactly zero and the above addToBulkMoles will have set
478 : // this for standard minerals, but not things related to sorption or kinetic minerals, so:
479 60 : _egs.setMineralRelatedFreeMoles(0.0);
480 :
481 : // Now need to swap all minerals out of the basis
482 60 : const std::vector<Real> & eqm_molality = _egs.getEquilibriumMolality();
483 60 : unsigned swap_into_basis = 0;
484 420 : for (unsigned i = 0; i < _num_basis; ++i)
485 360 : if (_mgd.basis_species_mineral[i])
486 : {
487 6 : const bool legitimate_swap_found = _egs.getSwapper().findBestEqmSwap(
488 6 : i, _mgd, eqm_molality, false, false, false, swap_into_basis);
489 6 : if (legitimate_swap_found)
490 : {
491 : try
492 : {
493 6 : _egs.performSwap(i, swap_into_basis);
494 : }
495 0 : catch (const MooseException & e)
496 : {
497 0 : mooseError(e.what());
498 0 : }
499 : }
500 : }
501 60 : }
502 :
503 : void
504 120 : GeochemistryTimeDependentReactor::preSolveFlush()
505 : {
506 120 : _new_temperature = newTemperature(_mole_additions);
507 :
508 : // Here we conserve mass, so compute the mass of the solution, without the free mineral moles.
509 : // We don't include the free mineral moles because users of GeochemistWorkbench will want
510 : // "flush" to operate like Bethke Eqn(13.14)
511 : // I assume we also don't include kinetic-mineral moles
512 120 : Real kg_in = _mole_additions(0) / GeochemistryConstants::MOLES_PER_KG_WATER;
513 1200 : for (unsigned i = 1; i < _num_basis; ++i)
514 1080 : if (!_mgd.basis_species_mineral[i])
515 294 : kg_in += _mole_additions(i) * _mgd.basis_species_molecular_weight[i] / 1000.0;
516 240 : for (unsigned kin = 0; kin < _num_kin; ++kin)
517 120 : if (!_mgd.kin_species_mineral[kin])
518 0 : kg_in += _mole_additions(kin + _num_basis) * _mgd.kin_species_molecular_weight[kin] / 1000.0;
519 :
520 120 : const std::vector<Real> & current_bulk = _egs.getBulkMolesOld();
521 120 : const std::vector<Real> & current_molal = _egs.getSolventMassAndFreeMolalityAndMineralMoles();
522 120 : const std::vector<Real> & kin_moles = _egs.getKineticMoles();
523 :
524 : // compute the current mass, without moles from free minerals and without kinetic minerals
525 120 : Real current_kg = current_bulk[0] / GeochemistryConstants::MOLES_PER_KG_WATER;
526 1200 : for (unsigned i = 1; i < _num_basis; ++i)
527 : {
528 : Real kinetic_contribution = 0.0;
529 2160 : for (unsigned k = 0; k < _num_kin; ++k)
530 1080 : if (_mgd.kin_species_mineral[k])
531 1080 : kinetic_contribution += kin_moles[k] * _mgd.kin_stoichiometry(k, i);
532 1080 : if (_mgd.basis_species_mineral[i])
533 786 : current_kg += (current_bulk[i] - current_molal[i] - kinetic_contribution) *
534 786 : _mgd.basis_species_molecular_weight[i] / 1000.0;
535 : else
536 294 : current_kg += (current_bulk[i] - kinetic_contribution) *
537 294 : _mgd.basis_species_molecular_weight[i] / 1000.0;
538 : }
539 :
540 120 : const Real fraction_to_remove = kg_in / current_kg;
541 1320 : for (unsigned i = 0; i < _num_basis; ++i)
542 : {
543 : Real all_kinetic_contribution = 0.0;
544 2400 : for (unsigned k = 0; k < _num_kin; ++k)
545 1200 : all_kinetic_contribution += kin_moles[k] * _mgd.kin_stoichiometry(k, i);
546 1200 : if (_mgd.basis_species_mineral[i])
547 786 : _mole_additions(i) -=
548 786 : fraction_to_remove * (current_bulk[i] - current_molal[i] - all_kinetic_contribution);
549 : else
550 414 : _mole_additions(i) -= fraction_to_remove * (current_bulk[i] - all_kinetic_contribution);
551 : }
552 240 : for (unsigned k = 0; k < _num_kin; ++k)
553 120 : if (!_mgd.kin_species_mineral[k])
554 0 : _mole_additions(k + _num_basis) -= fraction_to_remove * kin_moles[k];
555 120 : }
556 :
557 : void
558 66 : GeochemistryTimeDependentReactor::postSolveFlowThrough()
559 : {
560 : // copy the current_molal values into a new vector
561 66 : const std::vector<Real> current_molal = _egs.getSolventMassAndFreeMolalityAndMineralMoles();
562 : // remove minerals
563 594 : for (unsigned i = 1; i < _num_basis; ++i)
564 528 : if (_mgd.basis_species_mineral[i])
565 : {
566 84 : const Real to_remove = current_molal[i] - _small_molality * 10.0;
567 84 : _egs.addToBulkMoles(i, -to_remove);
568 84 : _minerals_dumped[_mgd.basis_species_name[i]] += to_remove;
569 : }
570 :
571 : // copy the current kinetic moles into a new vector
572 66 : const std::vector<Real> kin_moles = _egs.getKineticMoles();
573 : // remove minerals
574 66 : for (unsigned k = 0; k < _num_kin; ++k)
575 0 : if (_mgd.kin_species_mineral[k])
576 : {
577 0 : const Real to_remove = kin_moles[k] - _small_molality;
578 0 : for (unsigned i = 0; i < _num_basis; ++i)
579 0 : if (_mgd.kin_stoichiometry(k, i) != 0)
580 0 : _egs.addToBulkMoles(i,
581 0 : -_mgd.kin_stoichiometry(k, i) *
582 : to_remove); // remember bulk moles contains kinetic contributions
583 0 : _egs.setKineticMoles(k, _small_molality);
584 0 : _minerals_dumped[_mgd.kin_species_name[k]] += to_remove;
585 : }
586 :
587 66 : _egs.setMineralRelatedFreeMoles(_small_molality * 10.0);
588 66 : }
589 :
590 : void
591 0 : GeochemistryTimeDependentReactor::removeCurrentSpecies()
592 : {
593 0 : const std::vector<Real> & current_bulk = _egs.getBulkMolesOld();
594 0 : for (unsigned i = 0; i < _num_basis; ++i)
595 0 : _mole_additions(i) -= current_bulk[i];
596 0 : const std::vector<Real> & kin_moles = _egs.getKineticMoles();
597 0 : for (unsigned k = 0; k < _num_kin; ++k)
598 0 : _mole_additions(k + _num_basis) -= kin_moles[k];
599 0 : }
600 :
601 : void
602 0 : GeochemistryTimeDependentReactor::postSolveExchanger()
603 : {
604 : // remove precipitates
605 0 : postSolveFlowThrough();
606 :
607 0 : DenseVector<Real> zero_mole_additions(_num_basis + _num_kin);
608 0 : DenseMatrix<Real> zero_dmole_additions(_num_basis + _num_kin, _num_basis + _num_kin);
609 :
610 0 : const Real del_temperature = (_temperature[0] - _previous_temperature) / _heating_increments;
611 :
612 0 : for (unsigned incr = 0; incr < _heating_increments; ++incr)
613 : {
614 : // set temperature
615 0 : _new_temperature = _previous_temperature + del_temperature;
616 0 : _egs.setTemperature(_new_temperature);
617 0 : _egs.computeConsistentConfiguration();
618 0 : _previous_temperature = _new_temperature;
619 :
620 : zero_mole_additions.zero();
621 : zero_dmole_additions.zero();
622 :
623 : // solve the geochemical system
624 0 : _solver.solveSystem(_egs,
625 : _solver_output[0],
626 : _tot_iter[0],
627 : _abs_residual[0],
628 0 : _dt,
629 : zero_mole_additions,
630 : zero_dmole_additions);
631 :
632 : // remove precipitates
633 0 : postSolveFlowThrough();
634 : }
635 0 : }
|
