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 : // MOOSE includes
11 : #include "ReferenceResidualConvergence.h"
12 : #include "ReferenceResidualProblem.h"
13 : #include "FEProblemBase.h"
14 : #include "PetscSupport.h"
15 : #include "Executioner.h"
16 : #include "NonlinearSystemBase.h"
17 : #include "TaggingInterface.h"
18 : #include "AuxiliarySystem.h"
19 : #include "MooseVariableScalar.h"
20 : #include "NonlinearSystem.h"
21 :
22 : // PETSc includes
23 : #include <petscsnes.h>
24 :
25 : registerMooseObject("MooseApp", ReferenceResidualConvergence);
26 :
27 : InputParameters
28 15080 : ReferenceResidualConvergence::validParams()
29 : {
30 15080 : InputParameters params = DefaultNonlinearConvergence::validParams();
31 15080 : params += ReferenceResidualInterface::validParams();
32 :
33 15080 : params.addClassDescription(
34 : "Check the convergence of a problem with respect to a user-supplied reference solution."
35 : " Replaces ReferenceResidualProblem, currently still used in conjunction with it.");
36 :
37 15080 : return params;
38 0 : }
39 :
40 435 : ReferenceResidualConvergence::ReferenceResidualConvergence(const InputParameters & parameters)
41 : : DefaultNonlinearConvergence(parameters),
42 : ReferenceResidualInterface(this),
43 435 : _reference_vector(nullptr),
44 435 : _converge_on(getParam<std::vector<NonlinearVariableName>>("converge_on")),
45 435 : _zero_ref_type(
46 435 : getParam<MooseEnum>("zero_reference_residual_treatment").getEnum<ZeroReferenceType>()),
47 435 : _reference_vector_tag_id(Moose::INVALID_TAG_ID),
48 870 : _initialized(false)
49 : {
50 435 : if (parameters.isParamValid("solution_variables"))
51 : {
52 13 : if (parameters.isParamValid("reference_vector"))
53 13 : mooseDeprecated("The `solution_variables` parameter is deprecated, has no effect when "
54 : "the tagging system is used, and will be removed on January 1, 2020. "
55 : "Please simply delete this parameter from your input file.");
56 13 : _soln_var_names = parameters.get<std::vector<NonlinearVariableName>>("solution_variables");
57 : }
58 :
59 870 : if (parameters.isParamValid("reference_residual_variables") &&
60 435 : parameters.isParamValid("reference_vector"))
61 0 : mooseError(
62 : "You may specify either the `reference_residual_variables` "
63 : "or `reference_vector` parameter, not both. `reference_residual_variables` is deprecated "
64 : "so we recommend using `reference_vector`");
65 :
66 435 : if (parameters.isParamValid("reference_residual_variables"))
67 : {
68 0 : mooseDeprecated(
69 : "The save-in method for composing reference residual quantities is deprecated "
70 : "and will be removed on January 1, 2020. Please use the tagging system instead; "
71 : "specifically, please assign a TagName to the `reference_vector` parameter");
72 :
73 : _ref_resid_var_names =
74 0 : parameters.get<std::vector<AuxVariableName>>("reference_residual_variables");
75 :
76 0 : if (_soln_var_names.size() != _ref_resid_var_names.size())
77 0 : mooseError("Size of solution_variables (",
78 0 : _soln_var_names.size(),
79 : ") != size of reference_residual_variables (",
80 0 : _ref_resid_var_names.size(),
81 : ")");
82 : }
83 435 : else if (parameters.isParamValid("reference_vector"))
84 : {
85 391 : if (_fe_problem.numNonlinearSystems() > 1)
86 0 : paramError(
87 : "nl_sys_names",
88 : "reference residual problem does not currently support multiple nonlinear systems");
89 391 : _reference_vector_tag_id = _fe_problem.getVectorTagID(getParam<TagName>("reference_vector"));
90 383 : _reference_vector = &_fe_problem.getNonlinearSystemBase(0).getVector(_reference_vector_tag_id);
91 : }
92 : else
93 44 : mooseInfo("Neither the `reference_residual_variables` nor `reference_vector` parameter is "
94 : "specified, which means that no reference "
95 : "quantites are set. Because of this, the standard technique of comparing the "
96 : "norm of the full residual vector with its initial value will be used.");
97 :
98 427 : _accept_mult = parameters.get<Real>("acceptable_multiplier");
99 427 : _accept_iters = parameters.get<unsigned int>("acceptable_iterations");
100 :
101 : const auto norm_type_enum =
102 427 : parameters.get<MooseEnum>("normalization_type").getEnum<NormalizationType>();
103 427 : if (norm_type_enum == NormalizationType::LOCAL_L2)
104 : {
105 34 : _norm_type = libMesh::DISCRETE_L2;
106 34 : _local_norm = true;
107 : }
108 393 : else if (norm_type_enum == NormalizationType::GLOBAL_L2)
109 : {
110 341 : _norm_type = libMesh::DISCRETE_L2;
111 341 : _local_norm = false;
112 : }
113 52 : else if (norm_type_enum == NormalizationType::LOCAL_LINF)
114 : {
115 26 : _norm_type = libMesh::DISCRETE_L_INF;
116 26 : _local_norm = true;
117 : }
118 26 : else if (norm_type_enum == NormalizationType::GLOBAL_LINF)
119 : {
120 26 : _norm_type = libMesh::DISCRETE_L_INF;
121 26 : _local_norm = false;
122 : }
123 : else
124 : {
125 : mooseAssert(false, "This point should not be reached.");
126 : }
127 :
128 427 : if (_local_norm && !parameters.isParamValid("reference_vector"))
129 8 : paramError("reference_vector", "If local norm is used, a reference_vector must be provided.");
130 419 : }
131 :
132 : void
133 397 : ReferenceResidualConvergence::initialSetup()
134 : {
135 397 : DefaultNonlinearConvergence::initialSetup();
136 :
137 397 : NonlinearSystemBase & nonlinear_sys = _fe_problem.getNonlinearSystemBase(/*nl_sys=*/0);
138 397 : AuxiliarySystem & aux_sys = _fe_problem.getAuxiliarySystem();
139 397 : System & s = nonlinear_sys.system();
140 397 : auto & as = aux_sys.sys();
141 :
142 397 : if (_soln_var_names.empty())
143 : {
144 : // If the user provides reference_vector, that implies that they want the
145 : // individual variables compared against their reference quantities in the
146 : // tag vector. The code depends on having _soln_var_names populated,
147 : // so fill that out if they didn't specify solution_variables.
148 384 : if (_reference_vector)
149 1052 : for (unsigned int var_num = 0; var_num < s.n_vars(); var_num++)
150 694 : _soln_var_names.push_back(s.variable_name(var_num));
151 :
152 : // If they didn't provide reference_vector, that implies that they
153 : // want to skip the individual variable comparison, so leave it alone.
154 : }
155 13 : else if (_soln_var_names.size() != s.n_vars())
156 0 : mooseError("Size of solution_variables (",
157 0 : _soln_var_names.size(),
158 : ") != number of variables in system (",
159 0 : s.n_vars(),
160 : ")");
161 :
162 397 : const auto n_soln_vars = _soln_var_names.size();
163 397 : _variable_group_num_index.resize(n_soln_vars);
164 :
165 397 : if (!_converge_on.empty())
166 : {
167 52 : _converge_on_var.assign(n_soln_vars, false);
168 156 : for (std::size_t i = 0; i < n_soln_vars; ++i)
169 134 : for (const auto & c : _converge_on)
170 104 : if (MooseUtils::globCompare(_soln_var_names[i], c))
171 : {
172 74 : _converge_on_var[i] = true;
173 74 : break;
174 : }
175 : }
176 : else
177 345 : _converge_on_var.assign(n_soln_vars, true);
178 :
179 397 : unsigned int group_variable_num = 0;
180 397 : if (_use_group_variables)
181 : {
182 16 : for (unsigned int i = 0; i < _group_variables.size(); ++i)
183 : {
184 8 : group_variable_num += _group_variables[i].size();
185 8 : if (_group_variables[i].size() == 1)
186 0 : mooseError(" In the 'group_variables' parameter, variable ",
187 0 : _group_variables[i][0],
188 : " is not grouped with other variables.");
189 : }
190 :
191 8 : unsigned int size = n_soln_vars - group_variable_num + _group_variables.size();
192 8 : _group_ref_resid.resize(size);
193 8 : _group_resid.resize(size);
194 8 : _group_output_resid.resize(size);
195 8 : _group_soln_var_names.resize(size);
196 8 : _group_ref_resid_var_names.resize(size);
197 : }
198 : // If not using groups, use one group for each variable
199 : else
200 : {
201 389 : _group_ref_resid.resize(n_soln_vars);
202 389 : _group_resid.resize(n_soln_vars);
203 389 : _group_output_resid.resize(n_soln_vars);
204 389 : _group_soln_var_names.resize(n_soln_vars);
205 389 : _group_ref_resid_var_names.resize(n_soln_vars);
206 : }
207 :
208 397 : std::set<std::string> check_duplicate;
209 397 : if (_use_group_variables)
210 : {
211 16 : for (unsigned int i = 0; i < _group_variables.size(); ++i)
212 24 : for (unsigned int j = 0; j < _group_variables[i].size(); ++j)
213 16 : check_duplicate.insert(_group_variables[i][j]);
214 :
215 8 : if (check_duplicate.size() != group_variable_num)
216 0 : mooseError(
217 : "A variable cannot be included in multiple groups in the 'group_variables' parameter.");
218 : }
219 :
220 397 : _soln_vars.clear();
221 1104 : for (unsigned int i = 0; i < n_soln_vars; ++i)
222 : {
223 707 : bool found_match = false;
224 1043 : for (unsigned int var_num = 0; var_num < s.n_vars(); var_num++)
225 1043 : if (_soln_var_names[i] == s.variable_name(var_num))
226 : {
227 707 : _soln_vars.push_back(var_num);
228 707 : found_match = true;
229 707 : break;
230 : }
231 :
232 707 : if (!found_match)
233 0 : mooseError("Could not find solution variable '",
234 0 : _soln_var_names[i],
235 : "' in system '",
236 0 : s.name(),
237 : "'.");
238 : }
239 :
240 397 : if (!_reference_vector)
241 : {
242 26 : _ref_resid_vars.clear();
243 26 : for (unsigned int i = 0; i < _ref_resid_var_names.size(); ++i)
244 : {
245 0 : bool foundMatch = false;
246 0 : for (unsigned int var_num = 0; var_num < as.n_vars(); var_num++)
247 0 : if (_ref_resid_var_names[i] == as.variable_name(var_num))
248 : {
249 0 : _ref_resid_vars.push_back(var_num);
250 0 : foundMatch = true;
251 0 : break;
252 : }
253 :
254 0 : if (!foundMatch)
255 0 : mooseError("Could not find variable '", _ref_resid_var_names[i], "' in auxiliary system");
256 : }
257 : }
258 :
259 397 : unsigned int ungroup_index = 0;
260 397 : if (_use_group_variables)
261 8 : ungroup_index = _group_variables.size();
262 :
263 1088 : for (const auto i : index_range(_soln_vars))
264 : {
265 699 : bool find_group = false;
266 699 : if (_use_group_variables)
267 : {
268 8 : for (unsigned int j = 0; j < _group_variables.size(); ++j)
269 8 : if (std::find(_group_variables[j].begin(),
270 8 : _group_variables[j].end(),
271 16 : s.variable_name(_soln_vars[i])) != _group_variables[j].end())
272 : {
273 8 : if (!_converge_on_var[i])
274 8 : paramError("converge_on",
275 : "You added variable '",
276 8 : _soln_var_names[i],
277 : "' to a group but excluded it from the convergence check. This is not "
278 : "permitted.");
279 :
280 0 : _variable_group_num_index[i] = j;
281 0 : find_group = true;
282 0 : break;
283 : }
284 :
285 0 : if (!find_group)
286 : {
287 0 : _variable_group_num_index[i] = ungroup_index;
288 0 : ungroup_index++;
289 : }
290 : }
291 : else
292 691 : _variable_group_num_index[i] = i;
293 : }
294 :
295 389 : if (_use_group_variables)
296 : {
297 : // Check for variable groups containing both field and scalar variables
298 0 : for (unsigned int i = 0; i < _group_variables.size(); ++i)
299 : {
300 0 : unsigned int num_scalar_vars = 0;
301 0 : unsigned int num_field_vars = 0;
302 0 : if (_group_variables[i].size() > 1)
303 : {
304 0 : for (unsigned int j = 0; j < _group_variables[i].size(); ++j)
305 0 : for (unsigned int var_num = 0; var_num < s.n_vars(); var_num++)
306 0 : if (_group_variables[i][j] == s.variable_name(var_num))
307 : {
308 0 : if (nonlinear_sys.isScalarVariable(_soln_vars[var_num]))
309 0 : ++num_scalar_vars;
310 : else
311 0 : ++num_field_vars;
312 0 : break;
313 : }
314 : }
315 0 : if (num_scalar_vars > 0 && num_field_vars > 0)
316 0 : mooseWarning("In the 'group_variables' parameter, standard variables and scalar variables "
317 : "are grouped together in group ",
318 : i);
319 : }
320 : }
321 :
322 : // Keep track of the names of the variables in each group (of 1 variable if not using groups)
323 1080 : for (unsigned int i = 0; i < n_soln_vars; ++i)
324 : {
325 691 : if (_group_soln_var_names[_variable_group_num_index[i]].empty())
326 : {
327 691 : _group_soln_var_names[_variable_group_num_index[i]] = _soln_var_names[i];
328 691 : if (_use_group_variables && _variable_group_num_index[i] < _group_variables.size())
329 0 : _group_soln_var_names[_variable_group_num_index[i]] += " (grouped) ";
330 : }
331 :
332 691 : if (!_reference_vector && _group_ref_resid_var_names[_variable_group_num_index[i]].empty())
333 : {
334 0 : _group_ref_resid_var_names[_variable_group_num_index[i]] = _ref_resid_var_names[i];
335 0 : if (_use_group_variables && _variable_group_num_index[i] < _group_variables.size())
336 0 : _group_ref_resid_var_names[_variable_group_num_index[i]] += " (grouped) ";
337 : }
338 : }
339 :
340 389 : if (!_reference_vector)
341 : {
342 26 : const unsigned int size_soln_vars = _soln_vars.size();
343 26 : _scaling_factors.resize(size_soln_vars);
344 26 : for (unsigned int i = 0; i < size_soln_vars; ++i)
345 0 : if (nonlinear_sys.isScalarVariable(_soln_vars[i]))
346 0 : _scaling_factors[i] = nonlinear_sys.getScalarVariable(0, _soln_vars[i]).scalingFactor();
347 : else
348 0 : _scaling_factors[i] = nonlinear_sys.getVariable(/*tid*/ 0, _soln_vars[i]).scalingFactor();
349 : }
350 389 : _initialized = true;
351 389 : }
352 :
353 : void
354 22201 : ReferenceResidualConvergence::updateReferenceResidual()
355 : {
356 22201 : NonlinearSystemBase & _current_nl_sys = _fe_problem.currentNonlinearSystem();
357 22201 : AuxiliarySystem & aux_sys = _fe_problem.getAuxiliarySystem();
358 22201 : System & s = _current_nl_sys.system();
359 22201 : auto & as = aux_sys.sys();
360 :
361 22201 : std::fill(_group_resid.begin(), _group_resid.end(), 0.0);
362 22201 : std::fill(_group_output_resid.begin(), _group_output_resid.end(), 0.0);
363 22201 : if (_local_norm)
364 1144 : std::fill(_group_ref_resid.begin(), _group_ref_resid.end(), 1.0);
365 : else
366 21057 : std::fill(_group_ref_resid.begin(), _group_ref_resid.end(), 0.0);
367 :
368 65302 : for (const auto i : index_range(_soln_vars))
369 : {
370 43101 : Real resid = 0.0;
371 43101 : const auto group = _variable_group_num_index[i];
372 43101 : if (_local_norm)
373 : {
374 : mooseAssert(_current_nl_sys.RHS().size() == (*_reference_vector).size(),
375 : "Sizes of nonlinear RHS and reference vector should be the same.");
376 : mooseAssert((*_reference_vector).size(), "Reference vector must be provided.");
377 : // Add a tiny number to the reference to prevent a divide by zero.
378 2288 : auto ref = _reference_vector->clone();
379 2288 : ref->add(std::numeric_limits<Number>::min());
380 2288 : auto div = _current_nl_sys.RHS().clone();
381 2288 : *div /= *ref;
382 2288 : resid = Utility::pow<2>(s.calculate_norm(*div, _soln_vars[i], _norm_type));
383 2288 : }
384 : else
385 : {
386 40813 : resid = Utility::pow<2>(s.calculate_norm(_current_nl_sys.RHS(), _soln_vars[i], _norm_type));
387 40813 : if (_reference_vector)
388 : {
389 40813 : const auto ref_resid = s.calculate_norm(*_reference_vector, _soln_vars[i], _norm_type);
390 40813 : _group_ref_resid[group] += Utility::pow<2>(ref_resid);
391 : }
392 : }
393 :
394 43101 : _group_resid[group] += _converge_on_var[i] ? resid : 0;
395 43101 : _group_output_resid[group] += resid;
396 : }
397 :
398 22201 : if (!_reference_vector)
399 : {
400 484 : for (unsigned int i = 0; i < _ref_resid_vars.size(); ++i)
401 : {
402 : const auto ref_resid =
403 0 : as.calculate_norm(*as.current_local_solution, _ref_resid_vars[i], _norm_type) *
404 0 : _scaling_factors[i];
405 0 : _group_ref_resid[_variable_group_num_index[i]] += Utility::pow<2>(ref_resid);
406 : }
407 : }
408 :
409 65302 : for (unsigned int i = 0; i < _group_resid.size(); ++i)
410 : {
411 43101 : _group_resid[i] = std::sqrt(_group_resid[i]);
412 43101 : _group_output_resid[i] = std::sqrt(_group_output_resid[i]);
413 43101 : _group_ref_resid[i] = std::sqrt(_group_ref_resid[i]);
414 : }
415 22201 : }
416 :
417 : void
418 22201 : ReferenceResidualConvergence::nonlinearConvergenceSetup()
419 : {
420 22201 : if (!_initialized)
421 0 : initialSetup();
422 :
423 22201 : updateReferenceResidual();
424 :
425 22201 : std::ostringstream out;
426 :
427 22201 : if (_group_soln_var_names.size() > 0)
428 : {
429 21717 : out << std::setprecision(2) << std::scientific
430 21717 : << " Solution, reference convergence variable norms:\n";
431 21717 : unsigned int maxwsv = 0;
432 21717 : unsigned int maxwrv = 0;
433 64818 : for (unsigned int i = 0; i < _group_soln_var_names.size(); ++i)
434 : {
435 43101 : if (_group_soln_var_names[i].size() > maxwsv)
436 21717 : maxwsv = _group_soln_var_names[i].size();
437 43101 : if (!_reference_vector && _group_ref_resid_var_names[i].size() > maxwrv)
438 0 : maxwrv = _group_ref_resid_var_names[i].size();
439 : }
440 21717 : if (_reference_vector)
441 : // maxwrv is the width of maxwsv plus the length of "_ref" (e.g. 4)
442 21717 : maxwrv = maxwsv + 4;
443 :
444 64818 : for (unsigned int i = 0; i < _group_soln_var_names.size(); ++i)
445 : {
446 43101 : out << " " << std::setw(maxwsv + (_local_norm ? 5 : 2)) << std::left
447 43101 : << (_local_norm ? "norm " : "") + _group_soln_var_names[i] + ": ";
448 :
449 43101 : if (_group_output_resid[i] == _group_resid[i])
450 43001 : out << std::setw(8) << _group_output_resid[i];
451 : else
452 100 : out << std::setw(8) << _group_resid[i] << " (" << _group_output_resid[i] << ')';
453 :
454 43101 : if (!_local_norm)
455 : {
456 : const auto ref_var_name =
457 40813 : _reference_vector ? _group_soln_var_names[i] + "_ref" : _group_ref_resid_var_names[i];
458 81626 : out << " " << std::setw(maxwrv + 2) << ref_var_name + ":" << std::setw(8)
459 40813 : << _group_ref_resid[i] << " (" << std::setw(8)
460 40813 : << (_group_ref_resid[i] ? _group_resid[i] / _group_ref_resid[i] : _group_resid[i])
461 40813 : << ")";
462 40813 : }
463 43101 : out << '\n';
464 : }
465 21717 : _console << out.str() << std::flush;
466 : }
467 22201 : }
468 :
469 : bool
470 38108 : ReferenceResidualConvergence::checkConvergenceIndividVars(
471 : const Real fnorm,
472 : const Real abstol,
473 : const Real rtol,
474 : const Real initial_residual_before_preset_bcs)
475 : {
476 : // Convergence is checked via:
477 : // 1) if group residual is less than group reference residual by relative tolerance
478 : // 2) if group residual is less than absolute tolerance
479 : // 3) if group reference residual is zero and:
480 : // 3.1) Convergence type is ZERO_TOLERANCE and group residual is zero (rare, but possible, and
481 : // historically implemented that way)
482 : // 3.2) Convergence type is RELATIVE_TOLERANCE and group residual
483 : // is less than relative tolerance. (i.e., using the relative tolerance to check group
484 : // convergence in an absolute way)
485 :
486 38108 : bool convergedRelative = true;
487 38108 : if (_group_resid.size() > 0)
488 : {
489 113212 : for (unsigned int i = 0; i < _group_resid.size(); ++i)
490 75412 : convergedRelative &=
491 144138 : (_group_resid[i] < _group_ref_resid[i] * rtol || _group_resid[i] < abstol ||
492 68726 : (_group_ref_resid[i] == 0.0 &&
493 65368 : ((_zero_ref_type == ZeroReferenceType::ZERO_TOLERANCE && _group_resid[i] == 0.0) ||
494 65368 : (_zero_ref_type == ZeroReferenceType::RELATIVE_TOLERANCE &&
495 968 : _group_resid[i] <= rtol))));
496 : }
497 :
498 308 : else if (fnorm > initial_residual_before_preset_bcs * rtol)
499 88 : convergedRelative = false;
500 :
501 38108 : return convergedRelative;
502 : }
503 :
504 : bool
505 19929 : ReferenceResidualConvergence::checkRelativeConvergence(const unsigned int it,
506 : const Real fnorm,
507 : const Real the_residual,
508 : const Real rtol,
509 : const Real abstol,
510 : std::ostringstream & oss)
511 : {
512 19929 : if (checkConvergenceIndividVars(fnorm, abstol, rtol, the_residual))
513 : {
514 1750 : oss << "Converged due to function norm " << fnorm << " < relative tolerance (" << rtol
515 1750 : << ") or absolute tolerance (" << abstol << ") for all solution variables\n";
516 1750 : return true;
517 : }
518 36358 : else if (it >= _accept_iters &&
519 18179 : checkConvergenceIndividVars(
520 18179 : fnorm, abstol * _accept_mult, rtol * _accept_mult, the_residual))
521 : {
522 40 : oss << "Converged due to function norm " << fnorm << " < acceptable relative tolerance ("
523 40 : << rtol * _accept_mult << ") or acceptable absolute tolerance (" << abstol * _accept_mult
524 40 : << ") for all solution variables\n";
525 40 : _console << "Converged due to ACCEPTABLE tolerances" << std::endl;
526 40 : return true;
527 : }
528 :
529 18139 : return false;
530 : }
|
