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 "GeochemistrySpatialReactor.h"
11 :
12 : registerMooseObject("GeochemistryApp", GeochemistrySpatialReactor);
13 :
14 : InputParameters
15 344 : GeochemistrySpatialReactor::sharedParams()
16 : {
17 344 : InputParameters params = emptyInputParameters();
18 688 : params.addParam<unsigned>(
19 : "ramp_max_ionic_strength_subsequent",
20 688 : 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 688 : params.addParam<Real>("initial_temperature",
25 688 : 25,
26 : "Temperature at which the initial system is equilibrated. This is uniform "
27 : "over the entire mesh.");
28 688 : params.addCoupledVar("temperature", 25, "Temperature");
29 688 : params.addParam<Real>("close_system_at_time",
30 688 : 0.0,
31 : "Time at which to 'close' the entire spatial system, that is, change a "
32 : "kg_solvent_water constraint to moles_bulk_water, and all free_molality "
33 : "and free_moles_mineral_species to moles_bulk_species");
34 688 : params.addParam<std::vector<std::string>>(
35 : "remove_fixed_activity_name",
36 : {},
37 : "The name of the species that should have their activity or fugacity constraint removed at "
38 : "time given in remove_fixed_activity_time. There should be an equal number of these names "
39 : "as times given in remove_fixed_activity_time. Each of these must be in the basis and have "
40 : "an activity or fugacity constraint");
41 688 : params.addParam<std::vector<Real>>("remove_fixed_activity_time",
42 : {},
43 : "The times at which the species in remove_fixed_activity_name "
44 : "should have their activity or fugacity constraint removed.");
45 688 : params.addParam<std::vector<std::string>>(
46 : "source_species_names",
47 : {},
48 : "The name of the species that are added at rates given in source_species_rates. There must "
49 : "be an equal number of these as source_species_rates.");
50 688 : params.addCoupledVar("source_species_rates",
51 : "Rates, in mols/time_unit, of addition of the species with names given in "
52 : "source_species_names. A negative value corresponds to removing a species: "
53 : "be careful that you don't cause negative mass problems!");
54 688 : params.addParam<std::vector<std::string>>(
55 : "controlled_activity_name",
56 : {},
57 : "The names of the species that have their activity or fugacity constrained. There should be "
58 : "an equal number of these names as values given in controlled_activity_value. NOTE: if "
59 : "these species are not in the basis, or they do not have an activity (or fugacity) "
60 : "constraint then their activity cannot be controlled: in this case MOOSE will ignore the "
61 : "value you prescribe in controlled_activity_value.");
62 688 : params.addCoupledVar("controlled_activity_value",
63 : "Values of the activity or fugacity of the species in "
64 : "controlled_activity_name list. These should always be positive");
65 688 : params.addParam<bool>(
66 : "evaluate_kinetic_rates_always",
67 688 : true,
68 : "If true, then, evaluate the kinetic rates at every Newton step during the solve using the "
69 : "current values of molality, activity, etc (ie, implement an implicit solve). If false, "
70 : "then evaluate the kinetic rates using the values of molality, activity, etc, at the start "
71 : "of the current time step (ie, implement an explicit solve)");
72 688 : params.addParam<std::vector<std::string>>(
73 : "kinetic_species_name",
74 : {},
75 : "Names of the kinetic species given initial values in kinetic_species_initial_value");
76 688 : params.addParam<std::vector<Real>>(
77 : "kinetic_species_initial_value",
78 : {},
79 : "Initial number of moles, mass or volume (depending on kinetic_species_unit) for each of the "
80 : "species named in kinetic_species_name");
81 : MultiMooseEnum kin_species_unit("dimensionless moles molal kg g mg ug kg_per_kg_solvent "
82 344 : "g_per_kg_solvent mg_per_kg_solvent ug_per_kg_solvent cm3");
83 688 : params.addParam<MultiMooseEnum>(
84 : "kinetic_species_unit",
85 : kin_species_unit,
86 : "Units of the numerical values given in kinetic_species_initial_value. Moles: mole number. "
87 : "kg: kilograms. g: grams. mg: milligrams. ug: micrograms. cm3: cubic centimeters");
88 688 : params.addParam<bool>("adaptive_timestepping",
89 688 : false,
90 : "Use adaptive timestepping at each node in an attempt to ensure "
91 : "convergence of the solver. Setting this parameter to false saves some "
92 : "compute time because copying of datastructures is avoided");
93 688 : params.addParam<Real>(
94 : "dt_min",
95 688 : 1E-10,
96 : "If, during adaptive timestepping at a node, the time-step fails below this value, "
97 : "MOOSE will give up trying to solve the geochemical system. Optimally, you should set this "
98 : "value bearing abs_tol in mind because as dt changes, the initial value of the residual will "
99 : "also typically change. For example, if you set dt_min very small relative to abs_tol MOOSE "
100 : "may think the system has converged just because dt is small. Also, bear in mind your "
101 : "typical timestep size: if dt_min < 1E-16*typical_dt then you will run out of precision");
102 1032 : params.addRangeCheckedParam<Real>(
103 : "dt_dec",
104 688 : 0.5,
105 : "dt_dec >= 0 & dt_dec < 1.0",
106 : "If a geochemical solve fails, then 'adpative timestepping' at the node is initiated "
107 : "(assuming adaptive_timestepping = true): the time-step at the node is multiplied by this "
108 : "amount, and the solve process re-tried");
109 1032 : params.addRangeCheckedParam<Real>(
110 : "dt_inc",
111 688 : 1.1,
112 : "dt_inc >= 1.0",
113 : "If a geochemical solve suceeds during adpative timestepping at a node, then the time-step "
114 : "at the node is multiplied by this amount before performing the next adaptive timestep");
115 344 : return params;
116 344 : }
117 :
118 : InputParameters
119 236 : GeochemistrySpatialReactor::validParams()
120 : {
121 236 : InputParameters params = GeochemistryReactorBase::validParams();
122 236 : params += GeochemistrySpatialReactor::sharedParams();
123 236 : params.addClassDescription("UserObject that controls the space-dependent and time-dependent "
124 : "geochemistry reaction processes");
125 236 : return params;
126 0 : }
127 :
128 140 : GeochemistrySpatialReactor::GeochemistrySpatialReactor(const InputParameters & parameters)
129 : : GeochemistryReactorBase(parameters),
130 140 : _initial_temperature(getParam<Real>("initial_temperature")),
131 140 : _temperature(coupledValue("temperature")),
132 140 : _num_kin(_mgd.kin_species_name.size()),
133 : // NOTE: initialize _mgd_at_node before the swaps are performed
134 140 : _mgd_at_node(_num_my_nodes, _mgd),
135 : _egs_at_node(),
136 : // NOTE: the following implements the swaps in _mgd
137 1540 : _egs_copy(_mgd,
138 : _gac,
139 140 : _is,
140 140 : _swapper,
141 : getParam<std::vector<std::string>>("swap_out_of_basis"),
142 : getParam<std::vector<std::string>>("swap_into_basis"),
143 280 : getParam<std::string>("charge_balance_species"),
144 : getParam<std::vector<std::string>>("constraint_species"),
145 : getParam<std::vector<Real>>("constraint_value"),
146 : getParam<MultiMooseEnum>("constraint_unit"),
147 : getParam<MultiMooseEnum>("constraint_meaning"),
148 140 : _initial_temperature,
149 280 : getParam<unsigned>("extra_iterations_to_make_consistent"),
150 280 : getParam<Real>("min_initial_molality"),
151 : getParam<std::vector<std::string>>("kinetic_species_name"),
152 : getParam<std::vector<Real>>("kinetic_species_initial_value"),
153 : getParam<MultiMooseEnum>("kinetic_species_unit")),
154 420 : _solver(_mgd.basis_species_name.size(),
155 : _mgd.kin_species_name.size(),
156 : _is,
157 280 : getParam<Real>("abs_tol"),
158 280 : getParam<Real>("rel_tol"),
159 280 : getParam<unsigned>("max_iter"),
160 140 : getParam<Real>("max_initial_residual"),
161 140 : _small_molality,
162 140 : _max_swaps_allowed,
163 : getParam<std::vector<std::string>>("prevent_precipitation"),
164 280 : getParam<Real>("max_ionic_strength"),
165 140 : getParam<unsigned>("ramp_max_ionic_strength_initial"),
166 280 : getParam<bool>("evaluate_kinetic_rates_always")),
167 280 : _close_system_at_time(getParam<Real>("close_system_at_time")),
168 140 : _closed_system(false),
169 420 : _source_species_names(getParam<std::vector<std::string>>("source_species_names")),
170 140 : _num_source_species(_source_species_names.size()),
171 140 : _source_species_rates(0),
172 280 : _remove_fixed_activity_name(getParam<std::vector<std::string>>("remove_fixed_activity_name")),
173 420 : _remove_fixed_activity_time(getParam<std::vector<Real>>("remove_fixed_activity_time")),
174 140 : _num_removed_fixed(_remove_fixed_activity_name.size()),
175 140 : _removed_fixed_activity(_num_my_nodes, std::vector<bool>(_num_removed_fixed, false)),
176 420 : _controlled_activity_species_names(
177 : getParam<std::vector<std::string>>("controlled_activity_name")),
178 140 : _num_controlled_activity(_controlled_activity_species_names.size()),
179 140 : _controlled_activity_species_values(0),
180 140 : _mole_rates(_num_basis + _num_kin),
181 140 : _mole_additions(_num_my_nodes, DenseVector<Real>(_num_basis + _num_kin)),
182 140 : _dmole_additions(_num_my_nodes,
183 280 : DenseMatrix<Real>(_num_basis + _num_kin, _num_basis + _num_kin)),
184 280 : _ramp_subsequent(getParam<unsigned>("ramp_max_ionic_strength_subsequent")),
185 140 : _my_node_number(),
186 140 : _execute_done(_num_my_nodes, false),
187 280 : _adaptive_timestepping(getParam<bool>("adaptive_timestepping")),
188 159 : _dt_min(_adaptive_timestepping ? getParam<Real>("dt_min") : std::numeric_limits<Real>::max()),
189 280 : _dt_dec(getParam<Real>("dt_dec")),
190 280 : _dt_inc(getParam<Real>("dt_inc")),
191 140 : _nthreads(1)
192 : {
193 : // build _egs_at_node
194 916 : for (unsigned i = 0; i < _num_my_nodes; ++i)
195 776 : _egs_at_node.push_back(
196 10088 : GeochemicalSystem(_mgd_at_node[i],
197 : _gac,
198 : _is,
199 : _swapper,
200 : getParam<std::vector<std::string>>("swap_out_of_basis"),
201 : getParam<std::vector<std::string>>("swap_into_basis"),
202 1552 : getParam<std::string>("charge_balance_species"),
203 : getParam<std::vector<std::string>>("constraint_species"),
204 : getParam<std::vector<Real>>("constraint_value"),
205 : getParam<MultiMooseEnum>("constraint_unit"),
206 : getParam<MultiMooseEnum>("constraint_meaning"),
207 776 : _initial_temperature,
208 1552 : getParam<unsigned>("extra_iterations_to_make_consistent"),
209 1552 : getParam<Real>("min_initial_molality"),
210 : getParam<std::vector<std::string>>("kinetic_species_name"),
211 : getParam<std::vector<Real>>("kinetic_species_initial_value"),
212 : getParam<MultiMooseEnum>("kinetic_species_unit")));
213 :
214 : // check sources and set the rates
215 140 : if (coupledComponents("source_species_rates") != _num_source_species)
216 2 : paramError("source_species_names", "must have the same size as source_species_rates");
217 226 : for (unsigned i = 0; i < _num_source_species; ++i)
218 176 : _source_species_rates.push_back(&coupledValue("source_species_rates", i));
219 224 : for (unsigned i = 0; i < _num_source_species; ++i)
220 88 : if (!(_mgd.basis_species_index.count(_source_species_names[i]) == 1 ||
221 : _mgd.eqm_species_index.count(_source_species_names[i]) == 1 ||
222 : _mgd.kin_species_index.count(_source_species_names[i]) == 1))
223 4 : paramError("source_species_names",
224 2 : "The name " + _source_species_names[i] +
225 : " does not appear in the basis, equilibrium or kinetic species list");
226 :
227 : // check and set activity controllers
228 146 : for (unsigned i = 0; i < _num_removed_fixed; ++i)
229 : {
230 14 : if (_mgd.basis_species_index.count(_remove_fixed_activity_name[i]) == 0)
231 2 : paramError(
232 : "remove_fixed_activity_name",
233 : "The species ",
234 : _remove_fixed_activity_name[i],
235 : " is not in the basis, so cannot have its activity or fugacity constraint removed");
236 12 : else if (_num_my_nodes >
237 : 0) // don't consider the silly (but possible) case that there are no nodes
238 : {
239 12 : const unsigned basis_ind = _mgd.basis_species_index.at(_remove_fixed_activity_name[i]);
240 : const GeochemicalSystem::ConstraintMeaningEnum cm =
241 12 : _egs_at_node[0].getConstraintMeaning()[basis_ind];
242 12 : if (!(cm == GeochemicalSystem::ConstraintMeaningEnum::ACTIVITY ||
243 : cm == GeochemicalSystem::ConstraintMeaningEnum::FUGACITY))
244 2 : paramError("remove_fixed_activity_name",
245 : "The species ",
246 : _remove_fixed_activity_name[i],
247 : " is does not have an activity or fugacity constraint so cannot have such a "
248 : "constraint removed");
249 : }
250 : }
251 132 : if (_num_removed_fixed != _remove_fixed_activity_time.size())
252 2 : paramError("remove_fixed_activity_name",
253 : "must be of the same size as remove_fixed_activity_time");
254 130 : if (coupledComponents("controlled_activity_value") != _num_controlled_activity)
255 2 : paramError("controlled_activity_name", "must have the same size as controlled_activity_value");
256 156 : for (unsigned i = 0; i < _num_controlled_activity; ++i)
257 56 : _controlled_activity_species_values.push_back(&coupledValue("controlled_activity_value", i));
258 :
259 : // record coupled variables
260 : const std::vector<MooseVariableFEBase *> & coupled_vars = getCoupledMooseVars();
261 248 : for (unsigned int i = 0; i < coupled_vars.size(); i++)
262 240 : addMooseVariableDependency(coupled_vars[i]);
263 :
264 128 : buildMyNodeNumber();
265 128 : }
266 :
267 : void
268 128 : GeochemistrySpatialReactor::buildMyNodeNumber()
269 : {
270 128 : const MeshBase & msh = _subproblem.mesh().getMesh();
271 :
272 : // create the map from MOOSE's node numbering to the local node numbering (_my_node_number) used
273 : // by this object
274 : _my_node_number.clear();
275 : unsigned num_nodes_inserted = 0;
276 1136 : for (const auto & node : as_range(msh.local_nodes_begin(), msh.local_nodes_end()))
277 : {
278 752 : if (_my_node_number.count(node->id()) == 0)
279 752 : _my_node_number[node->id()] = num_nodes_inserted;
280 : else
281 0 : mooseError(
282 : "GeochemistrySpatialReactor: something wrong with node numbering in buildMyNodeNumber");
283 752 : num_nodes_inserted += 1;
284 128 : }
285 : mooseAssert(
286 : _my_node_number.size() == _num_my_nodes,
287 : "GeochemistrySpatialReactor: something wrong with node numbering in buildMyNodeNumber");
288 128 : }
289 :
290 : void
291 128 : GeochemistrySpatialReactor::initialSetup()
292 : {
293 : // Solve the geochemical system with its initial composition and with dt=0 so no kinetic additions
294 880 : for (unsigned i = 0; i < _num_my_nodes; ++i)
295 : {
296 752 : _mole_additions[i].zero();
297 : _dmole_additions[i].zero();
298 752 : _solver.solveSystem(_egs_at_node[i],
299 : _solver_output[i],
300 : _tot_iter[i],
301 : _abs_residual[i],
302 : 0.0,
303 : _mole_additions[i],
304 : _dmole_additions[i]);
305 : }
306 128 : }
307 :
308 : void
309 220 : GeochemistrySpatialReactor::initialize()
310 : {
311 : GeochemistryReactorBase::initialize();
312 220 : _execute_done.assign(_num_my_nodes, false);
313 220 : _nthreads = 1;
314 220 : }
315 :
316 : void
317 902 : GeochemistrySpatialReactor::execute()
318 : {
319 902 : if (_my_node_number.count(_current_node->id()) == 0)
320 0 : mooseError(
321 : "GeochemistrySpatialReactor: something wrong with node numbering in buildMyNodeNumber");
322 902 : const unsigned my_node_number = _my_node_number.at(_current_node->id());
323 :
324 : const unsigned aux_comp_number = 0; // component number to use for AuxVariables
325 : const ModelGeochemicalDatabase & mgd_ref =
326 902 : _egs_at_node[my_node_number].getModelGeochemicalDatabase();
327 :
328 : // close system
329 902 : if (!_closed_system && _t >= _close_system_at_time)
330 534 : _egs_at_node[my_node_number].closeSystem();
331 :
332 : // remove fixed-activity constraints.
333 1012 : for (unsigned i = 0; i < _num_removed_fixed; ++i)
334 : {
335 110 : if (!_removed_fixed_activity[my_node_number][i] && _t >= _remove_fixed_activity_time[i])
336 : {
337 : if (mgd_ref.basis_species_index.count(_remove_fixed_activity_name[i]))
338 55 : _egs_at_node[my_node_number].changeConstraintToBulk(
339 : mgd_ref.basis_species_index.at(_remove_fixed_activity_name[i]));
340 : _removed_fixed_activity[my_node_number][i] = true;
341 : }
342 : }
343 :
344 : // control activity
345 1078 : for (unsigned ca = 0; ca < _num_controlled_activity; ++ca)
346 : {
347 : const std::vector<GeochemicalSystem::ConstraintMeaningEnum> & cm =
348 176 : _egs_at_node[my_node_number].getConstraintMeaning();
349 176 : if (mgd_ref.basis_species_index.count(_controlled_activity_species_names[ca]))
350 : {
351 : const unsigned basis_ind =
352 176 : mgd_ref.basis_species_index.at(_controlled_activity_species_names[ca]);
353 176 : if (cm[basis_ind] == GeochemicalSystem::ConstraintMeaningEnum::ACTIVITY ||
354 : cm[basis_ind] == GeochemicalSystem::ConstraintMeaningEnum::FUGACITY)
355 176 : _egs_at_node[my_node_number].setConstraintValue(
356 176 : basis_ind, (*_controlled_activity_species_values[ca])[aux_comp_number]);
357 : }
358 : }
359 :
360 902 : _solver.setRampMaxIonicStrength(_ramp_subsequent);
361 :
362 902 : Real temperature0 = _egs_at_node[my_node_number].getTemperature();
363 902 : const Real temperature_rate = (_temperature[aux_comp_number] - temperature0) / _dt;
364 :
365 : // record the system in case of solve failures using the copy-assignment operator of
366 : // GeochemicalSystem
367 902 : if (_adaptive_timestepping)
368 38 : _egs_copy = _egs_at_node[my_node_number];
369 :
370 : Real done_dt = 0.0;
371 902 : Real my_dt = _dt;
372 :
373 : // the following loop implements adaptive timestepping at the node
374 1830 : while (done_dt < _dt)
375 : {
376 : // compute moles added in the current basis (the basis might change during adaptive
377 : // timestepping)
378 : _mole_rates.zero();
379 1577 : for (unsigned i = 0; i < _num_source_species; ++i)
380 : {
381 645 : const Real this_rate = (*_source_species_rates[i])[aux_comp_number];
382 : if (mgd_ref.basis_species_index.count(_source_species_names[i]))
383 : {
384 55 : const unsigned basis_ind = mgd_ref.basis_species_index.at(_source_species_names[i]);
385 55 : _mole_rates(basis_ind) += this_rate;
386 : }
387 : else if (mgd_ref.eqm_species_index.count(_source_species_names[i]))
388 : {
389 590 : const unsigned eqm_j = mgd_ref.eqm_species_index.at(_source_species_names[i]);
390 2360 : for (unsigned basis_ind = 0; basis_ind < _num_basis; ++basis_ind)
391 1770 : _mole_rates(basis_ind) += mgd_ref.eqm_stoichiometry(eqm_j, basis_ind) * this_rate;
392 : }
393 : else
394 : {
395 0 : const unsigned kin_ind = mgd_ref.kin_species_index.at(_source_species_names[i]);
396 0 : _mole_rates(_num_basis + kin_ind) += this_rate;
397 : }
398 : }
399 :
400 932 : Real temperature = temperature0 + my_dt * temperature_rate;
401 3728 : for (unsigned i = 0; i < _num_basis + _num_kin; ++i)
402 2796 : _mole_additions[my_node_number](i) = my_dt * _mole_rates(i);
403 : _dmole_additions[my_node_number].zero();
404 :
405 : // set temperature, if needed
406 932 : if (temperature != _egs_at_node[my_node_number].getTemperature())
407 : {
408 79 : _egs_at_node[my_node_number].setTemperature(temperature);
409 79 : _egs_at_node[my_node_number].computeConsistentConfiguration();
410 : }
411 :
412 : // solve the geochemical system
413 : try
414 : {
415 932 : _solver.solveSystem(_egs_at_node[my_node_number],
416 : _solver_output[my_node_number],
417 : _tot_iter[my_node_number],
418 : _abs_residual[my_node_number],
419 : my_dt,
420 : _mole_additions[my_node_number],
421 : _dmole_additions[my_node_number]);
422 913 : done_dt += my_dt;
423 913 : if (done_dt < _dt)
424 : {
425 15 : temperature0 = _egs_at_node[my_node_number].getTemperature();
426 : _egs_copy =
427 15 : _egs_at_node[my_node_number]; // use the copy-assignment operator of GeochemicalSystem
428 : }
429 913 : my_dt *= _dt_inc;
430 : }
431 19 : catch (const MooseException & e)
432 : {
433 : // use the copy-assignment operator of GeochemicalSystem to restore to the original
434 19 : if (_adaptive_timestepping)
435 17 : _egs_at_node[my_node_number] = _egs_copy;
436 19 : my_dt *= _dt_dec;
437 19 : if (my_dt < _dt_min)
438 8 : mooseException(
439 : "Geochemistry solve failed with dt = ", my_dt, " at node: ", _current_node->get_info());
440 19 : }
441 :
442 928 : if (done_dt + my_dt > _dt)
443 913 : my_dt = _dt - done_dt; // avoid overstepping
444 : }
445 :
446 : _execute_done[my_node_number] = true;
447 898 : }
448 :
449 : void
450 58 : GeochemistrySpatialReactor::threadJoin(const UserObject & uo)
451 : {
452 58 : _nthreads += 1;
453 : const auto & gsr = static_cast<const GeochemistrySpatialReactor &>(uo);
454 428 : for (unsigned i = 0; i < _num_my_nodes; ++i)
455 : {
456 370 : if (!_execute_done[i] && gsr._execute_done[i])
457 : {
458 398 : _solver_output[i].str("");
459 398 : _solver_output[i] << gsr._solver_output[i].str();
460 199 : _tot_iter[i] = gsr._tot_iter[i];
461 199 : _abs_residual[i] = gsr._abs_residual[i];
462 : _mole_additions[i] = gsr._mole_additions[i];
463 199 : _egs_at_node[i] = gsr._egs_at_node[i];
464 199 : _removed_fixed_activity[i] = gsr._removed_fixed_activity[i];
465 : // _mgd_at_node does not need to be threadJoined, because _egs_at_node[i] =
466 : // gsr._egs_at_node[i] uses the copy-assignment operator to copy the data in
467 : // _egs_at_node[i]._mgd
468 : }
469 : }
470 58 : }
471 :
472 : void
473 162 : GeochemistrySpatialReactor::finalize()
474 : {
475 : GeochemistryReactorBase::finalize();
476 : // if relevant, record that system is closed
477 162 : if (!_closed_system && _t >= _close_system_at_time)
478 89 : _closed_system = true;
479 : // ensure that the non-main threads have the main-thread's copy of _egs_at_node (and hence
480 : // _mgd_at_node) and _removed_fixed_activity, since the main-thread's copy has correctly gathered
481 : // all the information during threadJoin
482 220 : for (unsigned thrd = 1; thrd < _nthreads; ++thrd)
483 : {
484 : std::vector<GeochemistrySpatialReactor *> objects;
485 58 : _fe_problem.theWarehouse()
486 58 : .query()
487 58 : .condition<AttribSystem>("UserObject")
488 58 : .condition<AttribThread>(thrd)
489 58 : .condition<AttribName>(name())
490 : .queryInto(objects);
491 : mooseAssert(objects.size() == 1,
492 : "GeochemistrySpatialReactor::finalize() failed to obtain a single thread copy of "
493 : "the GeochemistrySpatialReactor");
494 58 : objects[0]->_removed_fixed_activity = _removed_fixed_activity;
495 58 : objects[0]->_egs_at_node = _egs_at_node;
496 58 : objects[0]->_closed_system = _closed_system;
497 58 : }
498 162 : }
499 :
500 : void
501 0 : GeochemistrySpatialReactor::meshChanged()
502 : {
503 0 : mooseError("GeochemistrySpatialReactor cannot yet handle adaptive meshing");
504 : /*
505 : Note to future coders:
506 : - have to rebuild _my_node_number. _num_my_nodes has to be changed (so its not const anymore).
507 : The Action must not just execute_on = EXEC_INITIAL for NearestNodeNumberUO
508 : - have to populate the new nodes correctly. This might be easiest if, at the start of execute,
509 : _egs_at_node was always populated using AuxVariables (probably a RequiredCoupledVar that is
510 : actually a ArrayVariableValue & (constructed with coupledArrayValue instead of just coupledValue)
511 : that record the kg_solvent_water, free molality, surface_pot_expr, etc. That is use
512 : _egs_at_node[i].setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles with the
513 : values from the AuxVariables. The reason for this design is that then MOOSE handles the
514 : interpolation of the AuxVariables to the newly-created nodes during mesh adaptivity. The
515 : difficult thing is to figure out the variables at the new nodes. Probably should copy the basis
516 : species, swap stuff, etc (_mgd_at_node) from the nearest original node to the newly-created node.
517 : And then solve the new system to allow basis swaps if appropriate. This is a bit tricky, hence it
518 : hasn't yet been implemented.
519 : */
520 : }
521 :
522 : const GeochemicalSystem &
523 18890 : GeochemistrySpatialReactor::getGeochemicalSystem(dof_id_type node_id) const
524 : {
525 : if (_my_node_number.count(node_id) != 1)
526 0 : mooseError("GeochemistrySpatialReactor does not know about node ", node_id);
527 18890 : return _egs_at_node[_my_node_number.at(node_id)];
528 : }
529 :
530 : const DenseVector<Real> &
531 0 : GeochemistrySpatialReactor::getMoleAdditions(dof_id_type node_id) const
532 : {
533 : if (_my_node_number.count(node_id) != 1)
534 0 : mooseError("GeochemistrySpatialReactor does not know about node ", node_id);
535 0 : return _mole_additions[_my_node_number.at(node_id)];
536 : }
537 :
538 : const std::stringstream &
539 0 : GeochemistrySpatialReactor::getSolverOutput(dof_id_type node_id) const
540 : {
541 : if (_my_node_number.count(node_id) != 1)
542 0 : mooseError("GeochemistrySpatialReactor does not know about node ", node_id);
543 0 : return _solver_output[_my_node_number.at(node_id)];
544 : }
545 :
546 : unsigned
547 116 : GeochemistrySpatialReactor::getSolverIterations(dof_id_type node_id) const
548 : {
549 : if (_my_node_number.count(node_id) != 1)
550 0 : mooseError("GeochemistrySpatialReactor does not know about node ", node_id);
551 116 : return _tot_iter[_my_node_number.at(node_id)];
552 : }
553 :
554 : Real
555 116 : GeochemistrySpatialReactor::getSolverResidual(dof_id_type node_id) const
556 : {
557 : if (_my_node_number.count(node_id) != 1)
558 0 : mooseError("GeochemistrySpatialReactor does not know about node ", node_id);
559 116 : return _abs_residual[_my_node_number.at(node_id)];
560 : }
561 :
562 : Real
563 0 : GeochemistrySpatialReactor::getMolesDumped(dof_id_type /*node_id*/,
564 : const std::string & /*species*/) const
565 : {
566 0 : return 0.0;
567 : }
|
