www.mooseframework.org
PicardSolve.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
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 "PicardSolve.h"
11 
12 #include "FEProblem.h"
13 #include "Executioner.h"
14 #include "MooseMesh.h"
15 #include "NonlinearSystem.h"
17 #include "Console.h"
18 #include "EigenExecutionerBase.h"
19 
20 template <>
23 {
25 
26  params.addParam<unsigned int>(
27  "picard_max_its",
28  1,
29  "Specifies the maximum number of Picard iterations. Mainly used when "
30  "wanting to do Picard iterations with MultiApps that are set to "
31  "execute_on timestep_end or timestep_begin. Setting this parameter to 1 turns off the Picard "
32  "iterations.");
33  params.addParam<bool>(
34  "accept_on_max_picard_iteration",
35  false,
36  "True to treat reaching the maximum number of Picard iterations as converged.");
37  params.addParam<bool>("disable_picard_residual_norm_check",
38  false,
39  "Disable the Picard residual norm evaluation thus the three parameters "
40  "picard_rel_tol, picard_abs_tol and picard_force_norms.");
41  params.addParam<Real>("picard_rel_tol",
42  1e-8,
43  "The relative nonlinear residual drop to shoot for "
44  "during Picard iterations. This check is "
45  "performed based on the Master app's nonlinear "
46  "residual.");
47  params.addParam<Real>("picard_abs_tol",
48  1e-50,
49  "The absolute nonlinear residual to shoot for "
50  "during Picard iterations. This check is "
51  "performed based on the Master app's nonlinear "
52  "residual.");
53  params.addParam<bool>(
54  "picard_force_norms",
55  false,
56  "Force the evaluation of both the TIMESTEP_BEGIN and TIMESTEP_END norms regardless of the "
57  "existance of active MultiApps with those execute_on flags, default: false.");
58 
59  params.addRangeCheckedParam<Real>("relaxation_factor",
60  1.0,
61  "relaxation_factor>0 & relaxation_factor<2",
62  "Fraction of newly computed value to keep."
63  "Set between 0 and 2.");
64  params.addParam<std::vector<std::string>>("relaxed_variables",
65  std::vector<std::string>(),
66  "List of variables to relax during Picard Iteration");
67 
68  params.addParamNamesToGroup("picard_max_its accept_on_max_picard_iteration "
69  "disable_picard_residual_norm_check picard_rel_tol "
70  "picard_abs_tol picard_force_norms "
71  "relaxation_factor relaxed_variables",
72  "Picard");
73 
74  params.addParam<unsigned int>(
75  "max_xfem_update",
76  std::numeric_limits<unsigned int>::max(),
77  "Maximum number of times to update XFEM crack topology in a step due to evolving cracks");
78  params.addParam<bool>("update_xfem_at_timestep_begin",
79  false,
80  "Should XFEM update the mesh at the beginning of the timestep");
81 
82  return params;
83 }
84 
86  : SolveObject(ex),
87  _picard_max_its(getParam<unsigned int>("picard_max_its")),
88  _has_picard_its(_picard_max_its > 1),
89  _accept_max_it(getParam<bool>("accept_on_max_picard_iteration")),
90  _has_picard_norm(!getParam<bool>("disable_picard_residual_norm_check")),
91  _picard_rel_tol(getParam<Real>("picard_rel_tol")),
92  _picard_abs_tol(getParam<Real>("picard_abs_tol")),
93  _picard_force_norms(getParam<bool>("picard_force_norms")),
94  _relax_factor(getParam<Real>("relaxation_factor")),
95  _relaxed_vars(getParam<std::vector<std::string>>("relaxed_variables")),
96  // this value will be set by MultiApp
97  _picard_self_relaxation_factor(1.0),
98  _max_xfem_update(getParam<unsigned int>("max_xfem_update")),
99  _update_xfem_at_timestep_begin(getParam<bool>("update_xfem_at_timestep_begin")),
100  _picard_timer(registerTimedSection("PicardSolve", 1)),
101  _picard_it(0),
102  _picard_status(MoosePicardConvergenceReason::UNSOLVED),
103  _xfem_update_count(0),
104  _xfem_repeat_step(false),
105  _previous_entering_time(_problem.time() - 1),
106  _solve_message(_problem.shouldSolve() ? "Solve Converged!" : "Solve Skipped!")
107 {
108  if (_relax_factor != 1.0)
109  // Store a copy of the previous solution here
110  _problem.getNonlinearSystemBase().addVector("relax_previous", false, PARALLEL);
111 }
112 
113 bool
115 {
116  TIME_SECTION(_picard_timer);
117 
118  Real current_dt = _problem.dt();
119 
124 
125  bool converged = true;
126 
127  // need to back up multi-apps even when not doing Picard iteration for recovering from failed
128  // multiapp solve
131 
132  // Prepare to relax variables as a master
133  std::set<dof_id_type> relaxed_dofs;
134  if (_relax_factor != 1.0)
135  {
136  // Snag all of the local dof indices for all of these variables
137  System & libmesh_nl_system = _nl.system();
138  AllLocalDofIndicesThread aldit(libmesh_nl_system, _relaxed_vars);
139  ConstElemRange & elem_range = *_problem.mesh().getActiveLocalElementRange();
140  Threads::parallel_reduce(elem_range, aldit);
141 
142  relaxed_dofs = aldit._all_dof_indices;
143  }
144 
145  // Prepare to relax variables as a subapp
146  std::set<dof_id_type> self_relaxed_dofs;
148  {
149  // Snag all of the local dof indices for all of these variables
150  System & libmesh_nl_system = _nl.system();
152  ConstElemRange & elem_range = *_problem.mesh().getActiveLocalElementRange();
153  Threads::parallel_reduce(elem_range, aldit);
154 
155  self_relaxed_dofs = aldit._all_dof_indices;
156 
158  {
159  NumericVector<Number> & solution = _nl.solution();
160  NumericVector<Number> & relax_previous = _nl.getVector("self_relax_previous");
161  relax_previous = solution;
162  }
163  }
164 
166  {
167  if (_has_picard_its)
168  {
169  if (_picard_it == 0)
170  {
171  if (_has_picard_norm)
172  {
173  // First Picard iteration - need to save off the initial nonlinear residual
175  _console << COLOR_MAGENTA << "Initial Picard Norm: " << COLOR_DEFAULT;
176  if (_picard_initial_norm == std::numeric_limits<Real>::max())
177  _console << " MAX ";
178  else
179  _console << std::scientific << _picard_initial_norm;
180  _console << COLOR_DEFAULT << "\n\n";
181  }
182  }
183  else
184  {
185  // For every iteration other than the first, we need to restore the state of the MultiApps
188  }
189 
190  _console << COLOR_MAGENTA << "Beginning Picard Iteration " << _picard_it << COLOR_DEFAULT
191  << '\n';
192  }
193 
194  Real begin_norm_old = (_picard_it > 0 ? _picard_timestep_begin_norm[_picard_it - 1]
195  : std::numeric_limits<Real>::max());
196  Real end_norm_old = (_picard_it > 0 ? _picard_timestep_end_norm[_picard_it - 1]
197  : std::numeric_limits<Real>::max());
198  bool relax = (_relax_factor != 1) && (_picard_it > 0);
199  bool solve_converged = solveStep(begin_norm_old,
201  end_norm_old,
203  relax,
204  relaxed_dofs);
205 
206  if (solve_converged)
207  {
208  if (_has_picard_its)
209  {
210  if (_has_picard_norm)
211  {
212  _console << "\n 0 Picard |R| = "
213  << Console::outputNorm(std::numeric_limits<Real>::max(), _picard_initial_norm)
214  << '\n';
215 
216  for (unsigned int i = 0; i <= _picard_it; ++i)
217  {
218  Real max_norm = std::max(_picard_timestep_begin_norm[i], _picard_timestep_end_norm[i]);
219  _console << std::setw(2) << i + 1
220  << " Picard |R| = " << Console::outputNorm(_picard_initial_norm, max_norm)
221  << '\n';
222  }
223 
224  Real max_norm = std::max(_picard_timestep_begin_norm[_picard_it],
226 
227  Real max_relative_drop = max_norm / _picard_initial_norm;
228 
229  if (max_norm < _picard_abs_tol)
230  {
232  break;
233  }
234  if (max_relative_drop < _picard_rel_tol)
235  {
237  break;
238  }
239  }
241  {
243  break;
244  }
245  if (_picard_it + 1 == _picard_max_its)
246  {
247  if (_accept_max_it)
248  {
250  converged = true;
251  }
252  else
253  {
255  converged = false;
256  }
257  break;
258  }
259  }
260  }
261  else
262  {
263  // If the last solve didn't converge then we need to exit this step completely (even in the
264  // case of Picard). So we can retry...
265  converged = false;
266  break;
267  }
268 
269  _problem.dt() =
270  current_dt; // _dt might be smaller than this at this point for multistep methods
271  }
272 
273  if (converged && _picard_self_relaxation_factor != 1.0)
274  {
276  {
277  NumericVector<Number> & solution = _nl.solution();
278  NumericVector<Number> & relax_previous = _nl.getVector("self_relax_previous");
279  Real factor = _picard_self_relaxation_factor;
280  for (const auto & dof : self_relaxed_dofs)
281  solution.set(dof, (relax_previous(dof) * (1.0 - factor)) + (solution(dof) * factor));
282  solution.close();
283  _nl.update();
284  }
286  }
287 
288  if (_has_picard_its)
289  {
290  _console << "Picard converged reason: ";
291  switch (_picard_status)
292  {
294  _console << "CONVERGED_ABS";
295  break;
297  _console << "CONVERGED_RELATIVE";
298  break;
300  _console << "CONVERGED_CUSTOM";
301  break;
303  _console << "REACH_MAX_ITS";
304  break;
306  _console << "DIVERGED_MAX_ITS";
307  break;
309  _console << "DIVERGED_NONLINEAR";
310  break;
312  _console << "DIVERGED_FAILED_MULTIAPP";
313  break;
314  default:
315  // UNSOLVED and CONVERGED_NONLINEAR should not be hit when Picard
316  // iteration is not on here
317  mooseError("Internal error: wrong Picard status!");
318  break;
319  }
320  _console << std::endl;
321  }
322  return converged;
323 }
324 
325 bool
326 PicardSolve::solveStep(Real begin_norm_old,
327  Real & begin_norm,
328  Real end_norm_old,
329  Real & end_norm,
330  bool relax,
331  const std::set<dof_id_type> & relaxed_dofs)
332 {
333  bool auto_advance = !(_has_picard_its && _problem.isTransient());
334 
335  if (dynamic_cast<EigenExecutionerBase *>(&_executioner) && _has_picard_its)
336  auto_advance = true;
337 
339 
341  if (!_problem.execMultiApps(EXEC_TIMESTEP_BEGIN, auto_advance))
342  {
344  return false;
345  }
346 
349 
351 
354  {
355  begin_norm = _problem.computeResidualL2Norm();
356 
357  _console << COLOR_MAGENTA << "Picard Norm after TIMESTEP_BEGIN MultiApps: "
358  << Console::outputNorm(begin_norm_old, begin_norm) << '\n';
359  }
360 
361  // Perform output for timestep begin
363 
364  // Update warehouse active objects
366 
367  if (relax)
368  {
369  NumericVector<Number> & solution = _nl.solution();
370  NumericVector<Number> & relax_previous = _nl.getVector("relax_previous");
371 
372  // Save off the current solution
373  relax_previous = solution;
374  }
375 
376  if (_has_picard_its)
377  _console << COLOR_MAGENTA << "\nMaster solve:\n" << COLOR_DEFAULT;
378  if (!_inner_solve->solve())
379  {
381 
382  _console << COLOR_RED << " Solve Did NOT Converge!" << COLOR_DEFAULT << std::endl;
383  // Perform the output of the current, failed time step (this only occurs if desired)
385  return false;
386  }
387  else
389 
390  _console << COLOR_GREEN << ' ' << _solve_message << COLOR_DEFAULT << std::endl;
391 
392  // Relax the "relaxed_variables"
393  if (relax)
394  {
395  NumericVector<Number> & solution = _nl.solution();
396  NumericVector<Number> & relax_previous = _nl.getVector("relax_previous");
397  Real factor = _relax_factor;
398  for (const auto & dof : relaxed_dofs)
399  solution.set(dof, (relax_previous(dof) * (1.0 - factor)) + (solution(dof) * factor));
400  solution.close();
401  _nl.update();
402  }
403 
405  {
406  _console << "XFEM modifying mesh, repeating step" << std::endl;
407  _xfem_repeat_step = true;
409  }
410  else
411  {
412  if (_problem.haveXFEM())
413  {
414  _xfem_repeat_step = false;
415  _xfem_update_count = 0;
416  _console << "XFEM not modifying mesh, continuing" << std::endl;
417  }
418 
421 
423  if (!_problem.execMultiApps(EXEC_TIMESTEP_END, auto_advance))
424  {
426  return false;
427  }
428  }
429 
431 
434  {
435  end_norm = _problem.computeResidualL2Norm();
436 
437  _console << COLOR_MAGENTA << "Picard Norm after TIMESTEP_END MultiApps: "
438  << Console::outputNorm(end_norm_old, end_norm) << '\n';
439  }
440 
441  return true;
442 }
const std::string _solve_message
Definition: PicardSolve.h:156
FEProblemBase & _problem
Reference to FEProblem.
Definition: SolveObject.h:41
Real _picard_self_relaxation_factor
Relaxation factor outside of Picard iteration (used as a subapp)
Definition: PicardSolve.h:122
const ExecFlagType EXEC_FAILED
ConstElemRange * getActiveLocalElementRange()
Return pointers to range objects for various types of ranges (local nodes, boundary elems...
Definition: MooseMesh.C:737
unsigned int _picard_max_its
Maximum Picard iterations.
Definition: PicardSolve.h:103
unsigned int _max_xfem_update
Maximum number of xfem updates per step.
Definition: PicardSolve.h:127
bool solveStep(Real begin_norm_old, Real &begin_norm, Real end_norm_old, Real &end_norm, bool relax, const std::set< dof_id_type > &relaxed_dofs)
Perform one Picard iteration or a full solve.
Definition: PicardSolve.C:326
virtual Real & time() const
std::vector< std::string > _relaxed_vars
The transferred variables that are going to be relaxed.
Definition: PicardSolve.h:119
NonlinearSystemBase & getNonlinearSystemBase()
PicardSolve(Executioner *ex)
Definition: PicardSolve.C:85
const PerfID _picard_timer
Timer for Picard iteration.
Definition: PicardSolve.h:133
std::vector< Real > _picard_timestep_end_norm
Full history of residual norm after evaluation of timestep_end.
Definition: PicardSolve.h:143
virtual NumericVector< Number > & addVector(const std::string &vector_name, const bool project, const ParallelType type)
Adds a solution length vector to the system.
Definition: SystemBase.C:542
virtual void preSolve()
Override this for actions that should take place before execution, called by PicardSolve.
Definition: Executioner.h:72
unsigned int _picard_it
Definition: PicardSolve.h:137
virtual void updateActiveObjects()
Update the active objects in the warehouses.
bool _xfem_repeat_step
Whether step should be repeated due to xfem modifying the mesh.
Definition: PicardSolve.h:151
virtual bool augmentedPicardConvergenceCheck() const
Augmented Picard convergence check that to be called by PicardSolve and can be overridden by derived ...
Definition: Executioner.h:114
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
bool _update_xfem_at_timestep_begin
Controls whether xfem should update the mesh at the beginning of the time step.
Definition: PicardSolve.h:129
bool _accept_max_it
Whether or not to treat reaching maximum number of Picard iteration as converged. ...
Definition: PicardSolve.h:107
virtual void onTimestepEnd() override
bool haveXFEM()
Find out whether the current analysis is using XFEM.
const ExecFlagType EXEC_TIMESTEP_END
void mooseError(Args &&... args) const
Definition: MooseObject.h:147
MoosePicardConvergenceReason _picard_status
Status of Picard solve.
Definition: PicardSolve.h:145
InputParameters emptyInputParameters()
virtual void update()
Update the system (doing libMesh magic)
Definition: SystemBase.C:1020
std::vector< std::string > _picard_self_relaxed_variables
Variables to be relaxed outside of Picard iteration (used as a subapp)
Definition: PicardSolve.h:124
Grab all the local dof indices for the variables passed in, in the system passed in.
virtual void execute(const ExecFlagType &exec_type)
Convenience function for performing execution of MOOSE systems.
InputParameters validParams< PicardSolve >()
Definition: PicardSolve.C:22
Executioner & _executioner
Executioner used to construct this.
Definition: SolveObject.h:39
bool _has_picard_its
Whether or not we activate Picard iteration.
Definition: PicardSolve.h:105
virtual bool solve() override
Picard solve the FEProblem.
Definition: PicardSolve.C:114
std::set< dof_id_type > _all_dof_indices
const ExecFlagType EXEC_TIMESTEP_BEGIN
bool _has_picard_norm
Whether or not to use residual norm to check the Picard convergence.
Definition: PicardSolve.h:109
SolveObject * _inner_solve
SolveObject wrapped by this solve object.
Definition: SolveObject.h:53
Real _picard_initial_norm
Initial residual norm.
Definition: PicardSolve.h:139
virtual bool updateMeshXFEM()
Update the mesh due to changing XFEM cuts.
Real _relax_factor
Relaxation factor for Picard Iteration.
Definition: PicardSolve.h:117
Executioners are objects that do the actual work of solving your problem.
Definition: Executioner.h:32
Real _picard_abs_tol
Absolute tolerance on residual norm.
Definition: PicardSolve.h:113
void backupMultiApps(ExecFlagType type)
Backup the MultiApps associated with the ExecFlagType.
NumericVector< Number > & solution() override
MoosePicardConvergenceReason
Enumeration for Picard convergence reasons.
Definition: PicardSolve.h:34
unsigned int _xfem_update_count
Counter for number of xfem updates that have been performed in the current step.
Definition: PicardSolve.h:149
virtual void postSolve()
Override this for actions that should take place after execution, called by PicardSolve.
Definition: Executioner.h:77
virtual System & system() override
Get the reference to the libMesh system.
std::vector< Real > _picard_timestep_begin_norm
Full history of residual norm after evaluation of timestep_begin.
Definition: PicardSolve.h:141
void restoreMultiApps(ExecFlagType type, bool force=false)
Restore the MultiApps associated with the ExecFlagType.
virtual Real computeResidualL2Norm()
Computes the residual using whatever is sitting in the current solution vector then returns the L2 no...
virtual MooseMesh & mesh() override
void addParam(const std::string &name, const S &value, const std::string &doc_string)
These methods add an option parameter and a documentation string to the InputParameters object...
void addRangeCheckedParam(const std::string &name, const T &value, const std::string &parsed_function, const std::string &doc_string)
const ConsoleStream _console
An instance of helper class to write streams to the Console objects.
bool execMultiApps(ExecFlagType type, bool auto_advance=true)
Execute the MultiApps associated with the ExecFlagType.
Real _picard_rel_tol
Relative tolerance on residual norm.
Definition: PicardSolve.h:111
virtual bool isTransient() const override
bool hasMultiApps() const
Returns whether or not the current simulation has any multiapps.
bool _picard_force_norms
Whether or not we force evaluation of residual norms even without multiapps.
Definition: PicardSolve.h:115
NonlinearSystemBase & _nl
Reference to nonlinear system base for faster access.
Definition: SolveObject.h:49
virtual bool solve()=0
Solve routine provided by this object.
virtual NumericVector< Number > & getVector(const std::string &name)
Get a raw NumericVector.
Definition: SystemBase.C:751
virtual Real & dt() const
virtual void outputStep(ExecFlagType type)
Output the current step.
Real _previous_entering_time
Time of previous Picard solve as a subapp.
Definition: PicardSolve.h:154
void execTransfers(ExecFlagType type)
Execute the Transfers associated with the ExecFlagType.
static std::string outputNorm(const Real &old_norm, const Real &norm)
A helper function for outputting norms in color.
Definition: Console.C:526
void addParamNamesToGroup(const std::string &space_delim_names, const std::string group_name)
This method takes a space delimited list of parameter names and adds them to the specified group name...