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 4100 : ReferenceResidualConvergence::validParams()
29 : {
30 4100 : InputParameters params = DefaultNonlinearConvergence::validParams();
31 4100 : params += ReferenceResidualInterface::validParams();
32 :
33 4100 : 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 4100 : return params;
38 0 : }
39 :
40 537 : ReferenceResidualConvergence::ReferenceResidualConvergence(const InputParameters & parameters)
41 : : DefaultNonlinearConvergence(parameters),
42 : ReferenceResidualInterface(this),
43 537 : _norm_type_enum(getParam<MooseEnum>("normalization_type")),
44 1074 : _accept_mult(getParam<Real>("acceptable_multiplier")),
45 1074 : _accept_iters(getParam<unsigned int>("acceptable_iterations")),
46 1074 : _nl_sys_num(_fe_problem.solverSysNum(getParam<SolverSystemName>("solver_sys"))),
47 537 : _residual_vector(nullptr),
48 537 : _reference_vector(nullptr),
49 537 : _zero_ref_type(
50 1074 : getParam<MooseEnum>("zero_reference_residual_treatment").getEnum<ZeroReferenceType>()),
51 1074 : _unscale_the_residual(getParam<bool>("unscale_the_residual")),
52 1611 : _reference_vector_tag_id(Moose::INVALID_TAG_ID)
53 : {
54 657 : if (_fe_problem.numNonlinearSystems() > 1 && !isParamSetByUser("solver_sys"))
55 0 : paramError("solver_sys",
56 : "Reference residual problem does not currently support multiple nonlinear systems "
57 : "in a single Convergence object. Multiple Convergence objects can be used, one for "
58 : "each nonlinear system, via the 'solver_sys' parameter.");
59 :
60 1074 : if (parameters.isParamValid("residual_vector"))
61 : {
62 : const auto residual_vector_tag_id =
63 52 : _fe_problem.getVectorTagID(getParam<TagName>("residual_vector"));
64 26 : _residual_vector =
65 26 : &_fe_problem.getNonlinearSystemBase(_nl_sys_num).getVector(residual_vector_tag_id);
66 : }
67 : else
68 511 : _residual_vector = &_fe_problem.getNonlinearSystemBase(_nl_sys_num).RHS();
69 :
70 1074 : if (parameters.isParamValid("reference_vector"))
71 : {
72 994 : _reference_vector_tag_id = _fe_problem.getVectorTagID(getParam<TagName>("reference_vector"));
73 491 : _reference_vector =
74 491 : &_fe_problem.getNonlinearSystemBase(_nl_sys_num).getVector(_reference_vector_tag_id);
75 : }
76 : else
77 40 : mooseDeprecated(
78 : "No `reference_vector` is provided, thus the Reference Residual convergence method will "
79 : "revert to default tolerance checking. `reference_vector` will become a required parameter "
80 : "on June 1st, 2027. If you are using `ReferenceResidualProblem`, either provide a "
81 : "reference_vector or use a standard problem type (e.g., remove "
82 : "Problem/type=ReferenceResidualProblem from your input file). If you are using "
83 : "`ReferenceResidualConvergence`, either provide a reference_vector or utilize "
84 : "`DefaultNonlinearConvergence` instead.");
85 :
86 531 : if (_norm_type_enum == "LOCAL_L2")
87 : {
88 82 : _norm_type = libMesh::DISCRETE_L2;
89 82 : _local_norm = true;
90 : }
91 449 : else if (_norm_type_enum == "GLOBAL_L2")
92 : {
93 397 : _norm_type = libMesh::DISCRETE_L2;
94 397 : _local_norm = false;
95 : }
96 52 : else if (_norm_type_enum == "LOCAL_LINF")
97 : {
98 26 : _norm_type = libMesh::DISCRETE_L_INF;
99 26 : _local_norm = true;
100 : }
101 26 : else if (_norm_type_enum == "GLOBAL_LINF")
102 : {
103 26 : _norm_type = libMesh::DISCRETE_L_INF;
104 26 : _local_norm = false;
105 : }
106 : else
107 : mooseAssert(false, "This point should not be reached.");
108 :
109 747 : if (_local_norm && !parameters.isParamValid("reference_vector"))
110 12 : paramError("reference_vector", "If local norm is used, a reference_vector must be provided.");
111 525 : }
112 :
113 : NonlinearSystemBase &
114 36599 : ReferenceResidualConvergence::nonlinearSystem()
115 : {
116 36599 : return _fe_problem.getNonlinearSystemBase(_nl_sys_num);
117 : }
118 :
119 : void
120 507 : ReferenceResidualConvergence::initialSetup()
121 : {
122 507 : DefaultNonlinearConvergence::initialSetup();
123 : // If no refernce_vector is provided, just revert to DefaultNonlinearConvergence behavior
124 507 : if (!_reference_vector)
125 26 : return;
126 :
127 481 : auto & nonlinear_sys = nonlinearSystem();
128 481 : auto & s = nonlinear_sys.system();
129 :
130 : // If the user provides reference_vector, that implies that they want the
131 : // individual variables compared against their reference quantities in the
132 : // tag vector. The code depends on having _soln_var_names populated,
133 : // so fill that out if they didn't specify solution_variables.
134 1700 : for (const auto var_num : make_range(s.n_vars()))
135 1219 : _soln_var_names.push_back(s.variable_name(var_num));
136 481 : const auto n_soln_vars = nonlinear_sys.nVariables();
137 :
138 962 : const auto converge_on = getParam<std::vector<NonlinearVariableName>>("converge_on");
139 481 : if (!converge_on.empty())
140 : {
141 146 : _converge_on_var.assign(n_soln_vars, false);
142 728 : for (std::size_t i = 0; i < n_soln_vars; ++i)
143 1578 : for (const auto & c : converge_on)
144 1432 : if (MooseUtils::globCompare(_soln_var_names[i], c))
145 : {
146 436 : _converge_on_var[i] = true;
147 436 : break;
148 : }
149 : }
150 : else
151 335 : _converge_on_var.assign(n_soln_vars, true);
152 :
153 481 : unsigned int num_variables_in_groups = 0;
154 587 : for (const auto i : index_range(_group_variables))
155 : {
156 106 : num_variables_in_groups += _group_variables[i].size();
157 106 : if (_group_variables[i].size() == 1)
158 0 : paramError("group_variables",
159 : "variable ",
160 0 : _group_variables[i][0],
161 : " is not grouped with other variables.");
162 : }
163 :
164 : // If no groups, size = n_soln_vars
165 481 : unsigned int n_groups = n_soln_vars - num_variables_in_groups + _group_variables.size();
166 481 : _group_ref_resid.resize(n_groups);
167 481 : _group_resid.resize(n_groups);
168 481 : _group_names.resize(n_groups);
169 481 : _converge_on_group.assign(n_groups, true);
170 481 : _scaling_factors.resize(n_soln_vars);
171 :
172 : // Check to make sure variables aren't in multiple groups
173 481 : if (_use_group_variables)
174 : {
175 56 : std::set<std::string> check_duplicate;
176 162 : for (const auto i : index_range(_group_variables))
177 368 : for (const auto j : index_range(_group_variables[i]))
178 262 : check_duplicate.insert(_group_variables[i][j]);
179 :
180 56 : if (check_duplicate.size() != num_variables_in_groups)
181 0 : paramError("group_variables", "A variable cannot be included in multiple groups.");
182 56 : }
183 :
184 481 : _soln_vars.clear();
185 1700 : for (const auto i : make_range(n_soln_vars))
186 : {
187 1219 : bool found_match = false;
188 2737 : for (const auto var_num : make_range(s.n_vars()))
189 2737 : if (_soln_var_names[i] == s.variable_name(var_num))
190 : {
191 1219 : _soln_vars.push_back(var_num);
192 1219 : found_match = true;
193 1219 : break;
194 : }
195 :
196 1219 : if (!found_match)
197 0 : mooseError("Could not find solution variable '",
198 0 : _soln_var_names[i],
199 : "' in system '",
200 0 : s.name(),
201 : "'.");
202 : }
203 :
204 481 : unsigned int ungroup_index = 0;
205 481 : if (_use_group_variables)
206 56 : ungroup_index = _group_variables.size();
207 :
208 : // Determine which group each variable belongs to
209 481 : _group_index.resize(n_soln_vars);
210 481 : _is_var_grouped.assign(n_soln_vars, false);
211 1688 : for (const auto i : index_range(_soln_vars))
212 : {
213 1213 : if (_use_group_variables)
214 : {
215 476 : for (const auto j : index_range(_group_variables))
216 436 : if (std::find(_group_variables[j].begin(),
217 436 : _group_variables[j].end(),
218 872 : s.variable_name(_soln_vars[i])) != _group_variables[j].end())
219 : {
220 256 : if (!_converge_on_var[i])
221 12 : paramError("converge_on",
222 : "You added variable '",
223 6 : _soln_var_names[i],
224 : "' to a group but excluded it from the convergence check. This is not "
225 : "permitted.");
226 :
227 250 : _group_index[i] = j;
228 250 : _is_var_grouped[i] = true;
229 250 : break;
230 : }
231 :
232 290 : if (!_is_var_grouped[i])
233 : {
234 40 : _group_index[i] = ungroup_index;
235 40 : ungroup_index++;
236 : }
237 : }
238 : else
239 917 : _group_index[i] = i;
240 : }
241 :
242 : // Check for variable groups containing both field and scalar variables
243 575 : for (const auto i : index_range(_group_variables))
244 : {
245 100 : unsigned int num_scalar_vars = 0;
246 100 : unsigned int num_field_vars = 0;
247 100 : if (_group_variables[i].size() > 1)
248 : {
249 350 : for (const auto j : index_range(_group_variables[i]))
250 770 : for (const auto var_num : make_range(s.n_vars()))
251 770 : if (_group_variables[i][j] == s.variable_name(var_num))
252 : {
253 250 : if (nonlinear_sys.isScalarVariable(_soln_vars[var_num]))
254 0 : ++num_scalar_vars;
255 : else
256 250 : ++num_field_vars;
257 250 : break;
258 : }
259 : }
260 100 : if (num_scalar_vars > 0 && num_field_vars > 0)
261 0 : paramWarning("group_variables",
262 : "standard variables and scalar variables are grouped together in group ",
263 : i);
264 : }
265 :
266 1532 : for (const auto i : index_range(_group_names))
267 : {
268 : // Accumulate names for a given group
269 1057 : std::vector<NonlinearVariableName> names;
270 4418 : for (const auto j : index_range(_group_index))
271 3361 : if (_group_index[j] == i)
272 : {
273 1207 : names.push_back(_soln_var_names[j]);
274 1207 : _converge_on_group[i] = _converge_on_group[i] && _converge_on_var[j];
275 : }
276 1057 : if (names.size() == 0)
277 0 : mooseError("Internal error, something is wrong with variable grouping");
278 1057 : else if (names.size() == 1)
279 957 : _group_names[i] = names[0];
280 : else
281 : {
282 100 : _group_names[i] = "(";
283 350 : for (const auto j : index_range(names))
284 : {
285 250 : _group_names[i] += names[j];
286 250 : if (j != names.size() - 1)
287 150 : _group_names[i] += ", ";
288 : }
289 100 : _group_names[i] += ")";
290 : }
291 1057 : }
292 475 : }
293 :
294 : void
295 17815 : ReferenceResidualConvergence::updateReferenceResidual()
296 : {
297 : // If no reference_vector is provided, this method is completely skipped
298 :
299 17815 : auto & nonlinear_sys = nonlinearSystem();
300 17815 : auto & s = nonlinear_sys.system();
301 :
302 54312 : for (const auto i : index_range(_scaling_factors))
303 36497 : if (nonlinear_sys.isScalarVariable(_soln_vars[i]))
304 0 : _scaling_factors[i] = nonlinear_sys.getScalarVariable(0, _soln_vars[i]).scalingFactor();
305 : else
306 36497 : _scaling_factors[i] = nonlinear_sys.getVariable(/*tid*/ 0, _soln_vars[i]).scalingFactor();
307 :
308 17815 : std::fill(_group_resid.begin(), _group_resid.end(), 0.0);
309 17815 : std::fill(_group_ref_resid.begin(), _group_ref_resid.end(), 0.0);
310 :
311 54312 : for (const auto i : index_range(_soln_vars))
312 : {
313 36497 : if (_converge_on_var[i])
314 : {
315 35921 : const auto group = _group_index[i];
316 :
317 : // Prepare residual
318 35921 : auto resid = Utility::pow<2>(s.calculate_norm(*_residual_vector, _soln_vars[i], _norm_type));
319 35921 : if (_unscale_the_residual)
320 : {
321 : mooseAssert(_scaling_factors[i], "Scaling factor must not be zero");
322 288 : resid /= Utility::pow<2>(_scaling_factors[i]);
323 : }
324 35921 : _group_resid[group] += resid;
325 :
326 : // Prepare reference residual. If local norm, this is actually the ratio of the residual
327 : // dividied by the reference at all DOF
328 : Real ref_resid;
329 35921 : if (_local_norm)
330 : {
331 : mooseAssert((*_residual_vector).size() == (*_reference_vector).size(),
332 : "Sizes of nonlinear RHS and reference vector should be the same.");
333 : mooseAssert((*_reference_vector).size(), "Reference vector must be provided.");
334 2754 : auto ref = _reference_vector->clone();
335 : // Add a tiny number to the reference to prevent a divide by zero.
336 2754 : ref->add(std::numeric_limits<Number>::min());
337 2754 : auto div = (*_residual_vector).clone();
338 2754 : *div /= *ref;
339 2754 : ref_resid = Utility::pow<2>(s.calculate_norm(*div, _soln_vars[i], _norm_type));
340 2754 : }
341 : else
342 : {
343 : ref_resid =
344 33167 : Utility::pow<2>(s.calculate_norm(*_reference_vector, _soln_vars[i], _norm_type));
345 33167 : if (_unscale_the_residual)
346 216 : ref_resid /= Utility::pow<2>(_scaling_factors[i]);
347 : }
348 35921 : _group_ref_resid[group] += ref_resid;
349 : }
350 : }
351 :
352 53691 : for (const auto i : index_range(_group_resid))
353 : {
354 35876 : _group_resid[i] = std::sqrt(_group_resid[i]);
355 35876 : _group_ref_resid[i] = std::sqrt(_group_ref_resid[i]);
356 : }
357 17815 : }
358 :
359 : void
360 18303 : ReferenceResidualConvergence::nonlinearConvergenceSetup()
361 : {
362 : // If no refernce_vector is provided, just revert to DefaultNonlinearConvergence behavior
363 18303 : if (!_reference_vector)
364 488 : return;
365 :
366 17815 : updateReferenceResidual();
367 :
368 17815 : std::ostringstream out;
369 17815 : out << _name << ": " << _norm_type_enum << " Reference Residual check\n";
370 :
371 17815 : if (_group_names.size() > 0)
372 : {
373 : // Set residual and references so that they always have a spacing of 8
374 17815 : out << std::setprecision(2) << std::scientific;
375 17815 : unsigned int var_space = 0;
376 53691 : for (const auto i : index_range(_group_names))
377 35876 : if (_group_names[i].size() > var_space)
378 17896 : var_space = _group_names[i].size();
379 :
380 53691 : for (const auto i : index_range(_group_names))
381 : {
382 35876 : if (_converge_on_group[i])
383 : {
384 : // Print residual
385 70600 : out << " " << std::setw(var_space + 8) << std::right << _group_names[i] + "-> res: "
386 35300 : << (_group_resid[i] < _abs_tol ? COLOR_YELLOW : COLOR_DEFAULT) << std::setw(8)
387 35300 : << _group_resid[i] << COLOR_DEFAULT;
388 :
389 : // Print res/ref ratio
390 35300 : if (_local_norm)
391 : out << " local res/ref: "
392 5508 : << (_group_resid[i] / _group_ref_resid[i] < _rel_tol ? COLOR_GREEN : COLOR_DEFAULT)
393 8262 : << std::setw(8) << _group_ref_resid[i] << COLOR_DEFAULT << "\n";
394 : else
395 : {
396 : // Print reference first if not local norm
397 32546 : out << " ref: " << std::setw(8) << _group_ref_resid[i] << " res/ref: ";
398 :
399 32546 : if (!_group_ref_resid[i])
400 28103 : out << _group_resid[i] << "\n";
401 : else
402 8886 : out << (_group_resid[i] / _group_ref_resid[i] < _rel_tol ? COLOR_GREEN : COLOR_DEFAULT)
403 8886 : << std::setw(8) << _group_resid[i] / _group_ref_resid[i] << COLOR_DEFAULT << "\n";
404 : }
405 : }
406 : }
407 17815 : _console << out.str() << std::flush;
408 : }
409 17815 : }
410 :
411 : bool
412 33563 : ReferenceResidualConvergence::checkConvergenceIndividVars(
413 : const Real /*fnorm*/,
414 : const Real abs_tol,
415 : const Real rel_tol,
416 : const Real /*initial_residual_before_preset_bcs*/)
417 : {
418 : // Convergence is checked via:
419 : // 1) Ratio of group residual to reference is less than relative tolerance
420 : // 2) if group residual is less than absolute tolerance
421 : // 3) if group reference residual is zero and:
422 : // 3.1) Convergence type is ZERO_TOLERANCE and group residual is zero (rare, but possible, and
423 : // historically implemented that way)
424 : // 3.2) Convergence type is RELATIVE_TOLERANCE and group residual
425 : // is less than relative tolerance. (i.e., using the relative tolerance to check group
426 : // convergence in an absolute way)
427 :
428 33563 : bool convergedRelative = true;
429 101466 : for (const auto i : index_range(_group_resid))
430 67903 : convergedRelative &=
431 62689 : ((!_local_norm && _group_resid[i] < _group_ref_resid[i] * rel_tol) ||
432 190240 : (_local_norm && _group_ref_resid[i] < rel_tol) || _group_resid[i] < abs_tol ||
433 59648 : (!_group_ref_resid[i] && !_local_norm &&
434 51918 : ((_zero_ref_type == ZeroReferenceType::ZERO_TOLERANCE && !_group_resid[i]) ||
435 51918 : (_zero_ref_type == ZeroReferenceType::RELATIVE_TOLERANCE &&
436 1970 : _group_resid[i] <= rel_tol))));
437 33563 : return convergedRelative;
438 : }
439 :
440 : bool
441 18145 : ReferenceResidualConvergence::checkResidualConvergence(const unsigned int n_iter,
442 : const Real fnorm,
443 : const Real ref_norm,
444 : const Real rel_tol,
445 : const Real abs_tol,
446 : std::ostringstream & oss)
447 : {
448 : // If no refernce_vector is provided, just revert to DefaultNonlinearConvergence behavior
449 18145 : if (!_reference_vector)
450 488 : return DefaultNonlinearConvergence::checkResidualConvergence(
451 488 : n_iter, fnorm, ref_norm, rel_tol, abs_tol, oss);
452 :
453 17657 : if (checkConvergenceIndividVars(fnorm, abs_tol, rel_tol, ref_norm))
454 : {
455 1697 : oss << "Converged normally";
456 1697 : return true;
457 : }
458 31866 : else if (n_iter >= _accept_iters &&
459 15906 : checkConvergenceIndividVars(
460 15906 : fnorm, abs_tol * _accept_mult, rel_tol * _accept_mult, ref_norm))
461 : {
462 : oss << " Converged due a larger acceptable tolerance due to `acceptible_multiplier` after "
463 36 : "`acceptible_iterations`.";
464 36 : _console << " Converged due to ACCEPTABLE tolerances" << std::endl;
465 36 : return true;
466 : }
467 :
468 15924 : return false;
469 : }
|
