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