www.mooseframework.org
ActuallyExplicitEuler.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 // MOOSE includes
11 #include "ActuallyExplicitEuler.h"
12 #include "NonlinearSystem.h"
13 #include "FEProblem.h"
14 #include "PetscSupport.h"
15 
16 // libMesh includes
17 #include "libmesh/sparse_matrix.h"
18 #include "libmesh/nonlinear_solver.h"
19 #include "libmesh/preconditioner.h"
20 #include "libmesh/enum_convergence_flags.h"
21 
23 
24 template <>
27 {
29 
30  MooseEnum solve_type("consistent lumped lump_preconditioned", "consistent");
31 
32  params.addParam<MooseEnum>(
33  "solve_type",
34  solve_type,
35  "The way to solve the system. A 'consistent' solve uses the full mass matrix and actually "
36  "needs to use a linear solver to solve the problem. 'lumped' uses a lumped mass matrix with "
37  "a simple inversion - incredibly fast but may be less accurate. 'lump_preconditioned' uses "
38  "the lumped mass matrix as a preconditioner for the 'consistent' solve");
39 
40  params.addClassDescription(
41  "Implementation of Explicit/Forward Euler without invoking any of the nonlinear solver");
42 
43  return params;
44 }
45 
49 class LumpedPreconditioner : public Preconditioner<Real>
50 {
51 public:
52  LumpedPreconditioner(const NumericVector<Real> & diag_inverse)
53  : Preconditioner(diag_inverse.comm()), _diag_inverse(diag_inverse)
54  {
55  }
56 
57  virtual void init() override
58  {
59  // No more initialization needed here
60  _is_initialized = true;
61  }
62 
63  virtual void apply(const NumericVector<Real> & x, NumericVector<Real> & y) override
64  {
65  y.pointwise_mult(_diag_inverse, x);
66  }
67 
68 protected:
70  const NumericVector<Real> & _diag_inverse;
71 };
72 
74  : TimeIntegrator(parameters),
75  MeshChangedInterface(parameters),
76  _solve_type(getParam<MooseEnum>("solve_type")),
77  _explicit_residual(_nl.addVector("explicit_residual", false, PARALLEL)),
78  _explicit_euler_update(_nl.addVector("explicit_euler_update", true, PARALLEL)),
79  _mass_matrix_diag(_nl.addVector("mass_matrix_diag", false, PARALLEL))
80 {
82 
83  // Try to keep MOOSE from doing any nonlinear stuff
85 
87  _ones = &_nl.addVector("ones", false, PARALLEL);
88 }
89 
90 void
92 {
93  meshChanged();
94 }
95 
96 void
98 {
99 }
100 
101 void
103 {
104 }
105 
106 void
108 {
109  if (!_sys.solutionUDot())
110  mooseError("ActuallyExplicitEuler: Time derivative of solution (`u_dot`) is not stored. Please "
111  "set uDotRequested() to true in FEProblemBase befor requesting `u_dot`.");
112 
113  NumericVector<Number> & u_dot = *_sys.solutionUDot();
114  u_dot = *_solution;
116  u_dot.close();
117 
118  _du_dot_du = 1.0 / _dt;
119 }
120 
121 void
122 ActuallyExplicitEuler::computeADTimeDerivatives(DualReal & ad_u_dot, const dof_id_type & dof) const
123 {
125 }
126 
127 void
129 {
130  auto & es = _fe_problem.es();
131 
132  auto & nonlinear_system = _fe_problem.getNonlinearSystemBase();
133 
134  auto & libmesh_system = dynamic_cast<NonlinearImplicitSystem &>(nonlinear_system.system());
135 
136  auto & mass_matrix = *libmesh_system.matrix;
137 
139 
140  // Set time back so that we're evaluating the interior residual at the old time
142 
143  libmesh_system.update();
144 
145  // Must compute the residual first
146  _explicit_residual.zero();
147  _fe_problem.computeResidual(*libmesh_system.current_local_solution, _explicit_residual);
148 
149  // The residual is on the RHS
150  _explicit_residual *= -1.;
151 
152  // Compute the mass matrix
153  _fe_problem.computeJacobianTag(*libmesh_system.current_local_solution, mass_matrix, _Ke_time_tag);
154 
155  // Still testing whether leaving the old update is a good idea or not
156  // _explicit_euler_update = 0;
157 
158  auto converged = false;
159 
160  switch (_solve_type)
161  {
162  case CONSISTENT:
163  {
164  const auto num_its_and_final_tol = _linear_solver->solve(
165  mass_matrix,
168  es.parameters.get<Real>("linear solver tolerance"),
169  es.parameters.get<unsigned int>("linear solver maximum iterations"));
170 
171  converged = checkLinearConvergence();
172 
173  _n_linear_iterations = num_its_and_final_tol.first;
174 
175  break;
176  }
177  case LUMPED:
178  {
179  // Computes the sum of each row (lumping)
180  // Note: This is actually how PETSc does it
181  // It's not "perfectly optimal" - but it will be fast (and universal)
182  mass_matrix.vector_mult(_mass_matrix_diag, *_ones);
183 
184  // "Invert" the diagonal mass matrix
185  _mass_matrix_diag.reciprocal();
186 
187  // Multiply the inversion by the RHS
189 
190  // Check for convergence by seeing if there is a nan or inf
191  auto sum = _explicit_euler_update.sum();
192  converged = std::isfinite(sum);
193 
195 
196  break;
197  }
198  case LUMP_PRECONDITIONED:
199  {
200  mass_matrix.vector_mult(_mass_matrix_diag, *_ones);
201  _mass_matrix_diag.reciprocal();
202 
203  const auto num_its_and_final_tol = _linear_solver->solve(
204  mass_matrix,
207  es.parameters.get<Real>("linear solver tolerance"),
208  es.parameters.get<unsigned int>("linear solver maximum iterations"));
209 
210  converged = checkLinearConvergence();
211 
212  _n_linear_iterations = num_its_and_final_tol.first;
213 
214  break;
215  }
216  default:
217  mooseError("Unknown solve_type in ActuallyExplicitEuler ");
218  }
219 
220  *libmesh_system.solution = nonlinear_system.solutionOld();
221  *libmesh_system.solution += _explicit_euler_update;
222 
223  // Enforce contraints on the solution
224  DofMap & dof_map = libmesh_system.get_dof_map();
225  dof_map.enforce_constraints_exactly(libmesh_system, libmesh_system.solution.get());
226 
227  libmesh_system.update();
228 
229  nonlinear_system.setSolution(*libmesh_system.current_local_solution);
230 
231  libmesh_system.nonlinear_solver->converged = converged;
232 }
233 
234 void
235 ActuallyExplicitEuler::postResidual(NumericVector<Number> & residual)
236 {
237  residual += _Re_time;
238  residual += _Re_non_time;
239  residual.close();
240 
241  // Reset time - the boundary conditions (which is what comes next) are applied at the final time
243 }
244 
245 void
247 {
248  // Can only be done after the system is inited
250  *_ones = 1.;
251 
253  _linear_solver = LinearSolver<Number>::build(comm());
254 
256  {
257  _preconditioner = libmesh_make_unique<LumpedPreconditioner>(_mass_matrix_diag);
258  _linear_solver->attach_preconditioner(_preconditioner.get());
259  _linear_solver->init();
260  }
261 
264 }
265 
266 bool
268 {
269  auto reason = _linear_solver->get_converged_reason();
270 
271  switch (reason)
272  {
273  case CONVERGED_RTOL_NORMAL:
274  case CONVERGED_ATOL_NORMAL:
275  case CONVERGED_RTOL:
276  case CONVERGED_ATOL:
277  case CONVERGED_ITS:
278  case CONVERGED_CG_NEG_CURVE:
279  case CONVERGED_CG_CONSTRAINED:
280  case CONVERGED_STEP_LENGTH:
281  case CONVERGED_HAPPY_BREAKDOWN:
282  return true;
283  case DIVERGED_NULL:
284  case DIVERGED_ITS:
285  case DIVERGED_DTOL:
286  case DIVERGED_BREAKDOWN:
287  case DIVERGED_BREAKDOWN_BICG:
288  case DIVERGED_NONSYMMETRIC:
289  case DIVERGED_INDEFINITE_PC:
290  case DIVERGED_NAN:
291  case DIVERGED_INDEFINITE_MAT:
292  case CONVERGED_ITERATING:
294  return false;
295  default:
296  mooseError("Unknown convergence flat in ActuallyExplicitEuler");
297  }
298 }
NonlinearSystemBase & _nl
virtual NumericVector< Number > * solutionUDot()=0
virtual void initialSetup() override
Called to setup datastructures.
virtual Real & time() const
SolverParams & solverParams()
Get the solver parameters.
virtual void postResidual(NumericVector< Number > &residual) override
Callback to the TimeIntegrator called immediately after the residuals are computed in NonlinearSystem...
NonlinearSystemBase & getNonlinearSystemBase()
NumericVector< Real > & _mass_matrix_diag
Diagonal of the lumped mass matrix (and its inversion)
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
FEProblemBase & _fe_problem
NumericVector< Real > * _ones
Just a vector of 1&#39;s to help with creating the lumped mass matrix.
InputParameters validParams< ActuallyExplicitEuler >()
NumericVector< Number > & _Re_non_time
residual vector for non-time contributions
SystemBase & _sys
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
DualNumber< Real, NumberArray< AD_MAX_DOFS_PER_ELEM, Real > > DualReal
Definition: DualReal.h:29
Solving a linear problem.
Definition: MooseTypes.h:595
void mooseError(Args &&... args) const
Definition: MooseObject.h:147
virtual void init() override
virtual void computeJacobianTag(const NumericVector< Number > &soln, SparseMatrix< Number > &jacobian, TagID tag)
Form a Jacobian matrix for a given tag.
static PetscErrorCode Vec x
bool checkLinearConvergence()
Check for the linear solver convergence.
virtual EquationSystems & es() override
void computeADTimeDerivatives(DualReal &ad_u_dot, const dof_id_type &dof) const override
method for computing local automatic differentiation time derivatives
Interface for notifications that the mesh has changed.
Real _current_time
Save off current time to reset it back and forth.
virtual void init() override
Called only before the very first timestep (t_step = 0) Never called again (not even during recover/r...
const NumericVector< Real > & _diag_inverse
The inverse of the diagonal of the lumped matrix.
DofMap & dof_map
InputParameters validParams< TimeIntegrator >()
unsigned int _n_linear_iterations
Total number of linear iterations over all stages of the time step.
This is a "smart" enum class intended to replace many of the shortcomings in the C++ enum type It sho...
Definition: MooseEnum.h:31
registerMooseObject("MooseApp", ActuallyExplicitEuler)
Moose::SolveType _type
Definition: SolverParams.h:19
NumericVector< Real > & _explicit_residual
Residual used for the RHS.
virtual void preSolve() override
std::unique_ptr< LumpedPreconditioner > _preconditioner
For solving with lumped preconditioning.
Real & _du_dot_du
solution vector for
const NumericVector< Number > *const & _solution
solution vectors
LumpedPreconditioner(const NumericVector< Real > &diag_inverse)
virtual void solve() override
Solves the time step and sets the number of nonlinear and linear iterations.
Base class for time integrators.
ActuallyExplicitEuler(const InputParameters &parameters)
virtual TagID getMatrixTagID(const TagName &tag_name)
Get a TagID from a TagName.
Definition: SubProblem.C:143
void computeTimeDerivativeHelper(T &u_dot, const T2 &u_old) const
Helper function that actually does the math for computing the time derivative.
MPI_Comm comm
NumericVector< Number > & _Re_time
residual vector for time contributions
void addClassDescription(const std::string &doc_string)
This method adds a description of the class that will be displayed in the input file syntax dump...
void setLinearSolverDefaults(FEProblemBase &problem, LinearSolver< T > &linear_solver)
Set the defaults for a libMesh LinearSolver.
Definition: PetscSupport.h:74
Implements a truly explicit (no nonlinear solve) first-order, forward Euler time integration scheme...
virtual Real & timeOld() const
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...
const NumericVector< Number > & _solution_old
std::unique_ptr< LinearSolver< Number > > _linear_solver
For solving with the consistent matrix.
virtual void apply(const NumericVector< Real > &x, NumericVector< Real > &y) override
void computeResidual(NonlinearImplicitSystem &sys, const NumericVector< Number > &soln, NumericVector< Number > &residual)
This function is called by Libmesh to form a residual.
virtual void computeTimeDerivatives() override
virtual void meshChanged() override
Called on this object when the mesh changes.
NumericVector< Real > & _explicit_euler_update
Solution vector for the linear solve.
TagID _Ke_time_tag
For computing the mass matrix.
Helper class to apply preconditioner.