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 453 : GeochemistrySpatialReactor::sharedParams()
16 : {
17 453 : InputParameters params = emptyInputParameters();
18 906 : params.addParam<unsigned>(
19 : "ramp_max_ionic_strength_subsequent",
20 906 : 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 906 : params.addParam<Real>("initial_temperature",
25 906 : 25,
26 : "Temperature at which the initial system is equilibrated. This is uniform "
27 : "over the entire mesh.");
28 906 : params.addCoupledVar("temperature", 25, "Temperature");
29 906 : params.addParam<Real>("close_system_at_time",
30 906 : 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 906 : 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 906 : 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 906 : 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 906 : 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 906 : 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 906 : 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 906 : params.addParam<bool>(
66 : "evaluate_kinetic_rates_always",
67 906 : 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 906 : 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 906 : 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 453 : "g_per_kg_solvent mg_per_kg_solvent ug_per_kg_solvent cm3");
83 906 : 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 906 : params.addParam<bool>("adaptive_timestepping",
89 906 : 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 906 : params.addParam<Real>(
94 : "dt_min",
95 906 : 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 1359 : params.addRangeCheckedParam<Real>(
103 : "dt_dec",
104 906 : 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 1359 : params.addRangeCheckedParam<Real>(
110 : "dt_inc",
111 906 : 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 453 : return params;
116 453 : }
117 :
118 : InputParameters
119 317 : GeochemistrySpatialReactor::validParams()
120 : {
121 317 : InputParameters params = GeochemistryReactorBase::validParams();
122 317 : params += GeochemistrySpatialReactor::sharedParams();
123 317 : params.addClassDescription("UserObject that controls the space-dependent and time-dependent "
124 : "geochemistry reaction processes");
125 317 : return params;
126 0 : }
127 :
128 193 : GeochemistrySpatialReactor::GeochemistrySpatialReactor(const InputParameters & parameters)
129 : : GeochemistryReactorBase(parameters),
130 193 : _initial_temperature(getParam<Real>("initial_temperature")),
131 193 : _temperature(coupledValue("temperature")),
132 193 : _num_kin(_mgd.kin_species_name.size()),
133 : // NOTE: initialize _mgd_at_node before the swaps are performed
134 193 : _mgd_at_node(_num_my_nodes, _mgd),
135 193 : _egs_at_node(),
136 : // NOTE: the following implements the swaps in _mgd
137 2123 : _egs_copy(_mgd,
138 : _gac,
139 193 : _is,
140 193 : _swapper,
141 : getParam<std::vector<std::string>>("swap_out_of_basis"),
142 : getParam<std::vector<std::string>>("swap_into_basis"),
143 386 : 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 193 : _initial_temperature,
149 386 : getParam<unsigned>("extra_iterations_to_make_consistent"),
150 386 : 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 579 : _solver(_mgd.basis_species_name.size(),
155 : _mgd.kin_species_name.size(),
156 : _is,
157 386 : getParam<Real>("abs_tol"),
158 386 : getParam<Real>("rel_tol"),
159 386 : getParam<unsigned>("max_iter"),
160 193 : getParam<Real>("max_initial_residual"),
161 193 : _small_molality,
162 193 : _max_swaps_allowed,
163 : getParam<std::vector<std::string>>("prevent_precipitation"),
164 386 : getParam<Real>("max_ionic_strength"),
165 193 : getParam<unsigned>("ramp_max_ionic_strength_initial"),
166 386 : getParam<bool>("evaluate_kinetic_rates_always")),
167 386 : _close_system_at_time(getParam<Real>("close_system_at_time")),
168 193 : _closed_system(false),
169 579 : _source_species_names(getParam<std::vector<std::string>>("source_species_names")),
170 193 : _num_source_species(_source_species_names.size()),
171 193 : _source_species_rates(0),
172 386 : _remove_fixed_activity_name(getParam<std::vector<std::string>>("remove_fixed_activity_name")),
173 579 : _remove_fixed_activity_time(getParam<std::vector<Real>>("remove_fixed_activity_time")),
174 193 : _num_removed_fixed(_remove_fixed_activity_name.size()),
175 193 : _removed_fixed_activity(_num_my_nodes, std::vector<bool>(_num_removed_fixed, false)),
176 579 : _controlled_activity_species_names(
177 : getParam<std::vector<std::string>>("controlled_activity_name")),
178 193 : _num_controlled_activity(_controlled_activity_species_names.size()),
179 193 : _controlled_activity_species_values(0),
180 193 : _mole_rates(_num_basis + _num_kin),
181 193 : _mole_additions(_num_my_nodes, DenseVector<Real>(_num_basis + _num_kin)),
182 193 : _dmole_additions(_num_my_nodes,
183 386 : DenseMatrix<Real>(_num_basis + _num_kin, _num_basis + _num_kin)),
184 386 : _ramp_subsequent(getParam<unsigned>("ramp_max_ionic_strength_subsequent")),
185 193 : _my_node_number(),
186 193 : _execute_done(_num_my_nodes, false),
187 386 : _adaptive_timestepping(getParam<bool>("adaptive_timestepping")),
188 220 : _dt_min(_adaptive_timestepping ? getParam<Real>("dt_min") : std::numeric_limits<Real>::max()),
189 386 : _dt_dec(getParam<Real>("dt_dec")),
190 386 : _dt_inc(getParam<Real>("dt_inc")),
191 193 : _nthreads(1)
192 : {
193 : // build _egs_at_node
194 1223 : for (unsigned i = 0; i < _num_my_nodes; ++i)
195 : _egs_at_node.push_back(
196 13390 : 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 2060 : 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 1030 : _initial_temperature,
208 2060 : getParam<unsigned>("extra_iterations_to_make_consistent"),
209 2060 : 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 193 : if (coupledComponents("source_species_rates") != _num_source_species)
216 2 : paramError("source_species_names", "must have the same size as source_species_rates");
217 311 : for (unsigned i = 0; i < _num_source_species; ++i)
218 240 : _source_species_rates.push_back(&coupledValue("source_species_rates", i));
219 309 : for (unsigned i = 0; i < _num_source_species; ++i)
220 120 : 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 203 : for (unsigned i = 0; i < _num_removed_fixed; ++i)
229 : {
230 18 : 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 16 : else if (_num_my_nodes >
237 : 0) // don't consider the silly (but possible) case that there are no nodes
238 : {
239 16 : const unsigned basis_ind = _mgd.basis_species_index.at(_remove_fixed_activity_name[i]);
240 : const GeochemicalSystem::ConstraintMeaningEnum cm =
241 16 : _egs_at_node[0].getConstraintMeaning()[basis_ind];
242 16 : 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 185 : 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 183 : if (coupledComponents("controlled_activity_value") != _num_controlled_activity)
255 2 : paramError("controlled_activity_name", "must have the same size as controlled_activity_value");
256 226 : for (unsigned i = 0; i < _num_controlled_activity; ++i)
257 90 : _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 350 : for (unsigned int i = 0; i < coupled_vars.size(); i++)
262 338 : addMooseVariableDependency(coupled_vars[i]);
263 :
264 181 : buildMyNodeNumber();
265 181 : }
266 :
267 : void
268 181 : GeochemistrySpatialReactor::buildMyNodeNumber()
269 : {
270 181 : 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 1549 : for (const auto & node : as_range(msh.local_nodes_begin(), msh.local_nodes_end()))
277 : {
278 1006 : if (_my_node_number.count(node->id()) == 0)
279 1006 : _my_node_number[node->id()] = num_nodes_inserted;
280 : else
281 0 : mooseError(
282 : "GeochemistrySpatialReactor: something wrong with node numbering in buildMyNodeNumber");
283 1006 : num_nodes_inserted += 1;
284 181 : }
285 : mooseAssert(
286 : _my_node_number.size() == _num_my_nodes,
287 : "GeochemistrySpatialReactor: something wrong with node numbering in buildMyNodeNumber");
288 181 : }
289 :
290 : void
291 181 : GeochemistrySpatialReactor::initialSetup()
292 : {
293 : // Solve the geochemical system with its initial composition and with dt=0 so no kinetic additions
294 1187 : for (unsigned i = 0; i < _num_my_nodes; ++i)
295 : {
296 1006 : _mole_additions[i].zero();
297 : _dmole_additions[i].zero();
298 1006 : _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 181 : }
307 :
308 : void
309 305 : GeochemistrySpatialReactor::initialize()
310 : {
311 : GeochemistryReactorBase::initialize();
312 305 : _execute_done.assign(_num_my_nodes, false);
313 305 : _nthreads = 1;
314 305 : }
315 :
316 : void
317 1110 : GeochemistrySpatialReactor::execute()
318 : {
319 1110 : if (_my_node_number.count(_current_node->id()) == 0)
320 0 : mooseError(
321 : "GeochemistrySpatialReactor: something wrong with node numbering in buildMyNodeNumber");
322 1110 : 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 1110 : _egs_at_node[my_node_number].getModelGeochemicalDatabase();
327 :
328 : // close system
329 1110 : if (!_closed_system && _t >= _close_system_at_time)
330 670 : _egs_at_node[my_node_number].closeSystem();
331 :
332 : // remove fixed-activity constraints.
333 1242 : for (unsigned i = 0; i < _num_removed_fixed; ++i)
334 : {
335 132 : 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 66 : _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 1352 : for (unsigned ca = 0; ca < _num_controlled_activity; ++ca)
346 : {
347 : const std::vector<GeochemicalSystem::ConstraintMeaningEnum> & cm =
348 242 : _egs_at_node[my_node_number].getConstraintMeaning();
349 242 : if (mgd_ref.basis_species_index.count(_controlled_activity_species_names[ca]))
350 : {
351 : const unsigned basis_ind =
352 242 : mgd_ref.basis_species_index.at(_controlled_activity_species_names[ca]);
353 242 : if (cm[basis_ind] == GeochemicalSystem::ConstraintMeaningEnum::ACTIVITY ||
354 : cm[basis_ind] == GeochemicalSystem::ConstraintMeaningEnum::FUGACITY)
355 242 : _egs_at_node[my_node_number].setConstraintValue(
356 242 : basis_ind, (*_controlled_activity_species_values[ca])[aux_comp_number]);
357 : }
358 : }
359 :
360 1110 : _solver.setRampMaxIonicStrength(_ramp_subsequent);
361 :
362 1110 : Real temperature0 = _egs_at_node[my_node_number].getTemperature();
363 1110 : 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 1110 : if (_adaptive_timestepping)
368 44 : _egs_copy = _egs_at_node[my_node_number];
369 :
370 : Real done_dt = 0.0;
371 1110 : Real my_dt = _dt;
372 :
373 : // the following loop implements adaptive timestepping at the node
374 2252 : 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 1920 : for (unsigned i = 0; i < _num_source_species; ++i)
380 : {
381 774 : 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 66 : const unsigned basis_ind = mgd_ref.basis_species_index.at(_source_species_names[i]);
385 66 : _mole_rates(basis_ind) += this_rate;
386 : }
387 : else if (mgd_ref.eqm_species_index.count(_source_species_names[i]))
388 : {
389 708 : const unsigned eqm_j = mgd_ref.eqm_species_index.at(_source_species_names[i]);
390 2832 : for (unsigned basis_ind = 0; basis_ind < _num_basis; ++basis_ind)
391 2124 : _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 1146 : Real temperature = temperature0 + my_dt * temperature_rate;
401 4584 : for (unsigned i = 0; i < _num_basis + _num_kin; ++i)
402 3438 : _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 1146 : if (temperature != _egs_at_node[my_node_number].getTemperature())
407 : {
408 94 : _egs_at_node[my_node_number].setTemperature(temperature);
409 94 : _egs_at_node[my_node_number].computeConsistentConfiguration();
410 : }
411 :
412 : // solve the geochemical system
413 : try
414 : {
415 1146 : _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 1124 : done_dt += my_dt;
423 1124 : if (done_dt < _dt)
424 : {
425 18 : temperature0 = _egs_at_node[my_node_number].getTemperature();
426 : _egs_copy =
427 18 : _egs_at_node[my_node_number]; // use the copy-assignment operator of GeochemicalSystem
428 : }
429 1124 : my_dt *= _dt_inc;
430 : }
431 22 : catch (const MooseException & e)
432 : {
433 : // use the copy-assignment operator of GeochemicalSystem to restore to the original
434 22 : if (_adaptive_timestepping)
435 20 : _egs_at_node[my_node_number] = _egs_copy;
436 22 : my_dt *= _dt_dec;
437 22 : if (my_dt < _dt_min)
438 8 : mooseException(
439 : "Geochemistry solve failed with dt = ", my_dt, " at node: ", _current_node->get_info());
440 22 : }
441 :
442 1142 : if (done_dt + my_dt > _dt)
443 1124 : my_dt = _dt - done_dt; // avoid overstepping
444 : }
445 :
446 : _execute_done[my_node_number] = true;
447 1106 : }
448 :
449 : void
450 99 : GeochemistrySpatialReactor::threadJoin(const UserObject & uo)
451 : {
452 99 : _nthreads += 1;
453 : const auto & gsr = static_cast<const GeochemistrySpatialReactor &>(uo);
454 655 : for (unsigned i = 0; i < _num_my_nodes; ++i)
455 : {
456 556 : if (!_execute_done[i] && gsr._execute_done[i])
457 : {
458 600 : _solver_output[i].str("");
459 600 : _solver_output[i] << gsr._solver_output[i].str();
460 300 : _tot_iter[i] = gsr._tot_iter[i];
461 300 : _abs_residual[i] = gsr._abs_residual[i];
462 : _mole_additions[i] = gsr._mole_additions[i];
463 300 : _egs_at_node[i] = gsr._egs_at_node[i];
464 300 : _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 99 : }
471 :
472 : void
473 206 : GeochemistrySpatialReactor::finalize()
474 : {
475 : GeochemistryReactorBase::finalize();
476 : // if relevant, record that system is closed
477 206 : if (!_closed_system && _t >= _close_system_at_time)
478 115 : _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 305 : for (unsigned thrd = 1; thrd < _nthreads; ++thrd)
483 : {
484 : std::vector<GeochemistrySpatialReactor *> objects;
485 99 : _fe_problem.theWarehouse()
486 99 : .query()
487 99 : .condition<AttribSystem>("UserObject")
488 99 : .condition<AttribThread>(thrd)
489 99 : .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 99 : objects[0]->_removed_fixed_activity = _removed_fixed_activity;
495 99 : objects[0]->_egs_at_node = _egs_at_node;
496 99 : objects[0]->_closed_system = _closed_system;
497 : }
498 206 : }
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 23286 : 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 23286 : 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 144 : 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 144 : return _tot_iter[_my_node_number.at(node_id)];
552 : }
553 :
554 : Real
555 144 : 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 144 : 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 : }
|
