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