Skip to content

Commit

Permalink
Extracted base class ExplicitTimeIntegrator from ActuallyExplicitEuler
Browse files Browse the repository at this point in the history 
Refs idaholab#11855
  • Loading branch information
@joshuahansel
joshuahansel committed Sep 27, 2019
1 parent e6e9216 commit 78704e2
Show file tree
Hide file tree
Showing 6 changed files with 455 additions and 346 deletions.
95 changes: 2 additions & 93 deletions framework/doc/content/source/timeintegrators/ActuallyExplicitEuler.md
View file
Original file line number Diff line number Diff line change
@@ -1,101 +1,10 @@
# ActuallyExplicitEuler

!syntax description /Executioner/TimeIntegrator/ActuallyExplicitEuler

## Description

`ActuallyExplicitEuler` implements true, first-order, forward Euler. It differentiates itself from [/ExplicitEuler.md] by implementing solution methods that don't involve the "nonlinear" solver in MOOSE at all and, in the case of `lumped`, actually don't involve a linear solve at all. This makes it incredibly fast.

## `solve_type`

This object has three different ways of solving. You can switch between them by using the `solve_type` parameter. The three methods are:

1. `consistent`: (The Default) A full mass matrix is built and used in a linear solve for the update
2. `lumped`: A "lumped" mass matrix is built, hand-inverted and applied to the RHS - fast, but possibly less accurate
3. `lump_preconditioned`: The inversion of the "lumped" mass matrix is used to precondition the `consistent` solve. Right now this is still experimental.

All three methods are solved similarly: by solving for the update $\delta u$ then adding that to the existing solution.

Below is some more explanation of each of these `solve_type` options:

### `consistent`

The `consistent` option builds a full ("consistent") "mass matrix" and uses it in a linear solve to get the update. This is done by calling `FEProblem::computeJacobianTag()` and specifying the `TIME` tag which includes all of the `TimeKernel` derived Kernels and `NodalBC` derived BoundaryConditions to compute $\mathbf{M}$:

!listing framework/src/timeintegrators/ActuallyExplicitEuler.C line=computeJacobianTag

A residual computation is also completed to use as the RHS ($R$):

!listing framework/src/timeintegrators/ActuallyExplicitEuler.C line=computeResidual

Creating the equation:

\begin{equation}
\mathbf{M} \delta u = -R
\end{equation}

This is solved using the default linear solver from libMesh (usually PETSc... including the application of any command-line parameters specified).

### `lumped`

The `lumped` option creates a "lumped mass matrix" to use in the solve. A lumped mass matrix is a diagonal matrix where the diagonal is the sum of all elements on the row from the original matrix.

To achieve the lumping here a matrix-vector product is performed with a vector of all ones. That creates a vector where each entry is what would be on the diagonal of the lumped mass matrix.

The inverse of a diagonal matrix is simply the reciprocal of each diagonal entry - easily applied to our vector. Then the matrix-vector product of the "inverse" lumped diagonal matrix is applied by simply doing a pointwise multiplication with the RHS.

This means that the `lumped` option actually doesn't need to solve a system of linear equations at all... making it incredibly fast. However, the use of a lumped mass matrix may lead to unacceptable phase errors.

### `lump_preconditioned`

This option is the combination of the above two. The consistent mass matrix is built and used to solve... but the preconditioner is applied as simply the inverse of the lumped mass matrix. This means that solving the true (consistent) system can be done with simply using point-wise multiplications. This makes it incredibly fast and memory efficient while still accurate.

## Advanced Details

A few notes on some of the implementation details of this object:

### Update Form

Note that even though we're doing an explicit solve we are currently doing it in "update form" similar to a single step Newton solve. This gives us good parity with the rest of MOOSE. We may want to change this in the future to make better use of the fact that the mass-matrix can be constant for a wider class of problems if we remove `dt` from it.

### `_ones`

To get the sum of each row of the mass matrix for "lumping" purposes a vector consisting of all `1`s is used in a matrix-vector product:

!listing framework/src/timeintegrators/ActuallyExplicitEuler.C line=mass_matrix.vector_mult

This is actually the very same way `MatGetRowSum` is implemented in PETSc. Doing it ourselves though cuts down on vector creation/destruction and a few other bookkeeping bits.

In the future we might change this to use `MatGetRowSum` if a specialization for `MPI_Aij` format is created

### Time

Time in an explicit solve must be handled carefully. When evaluating the weak form (the integral parts) of the residual time needs to actually be the "old" time (the time we are solving "from"):

!listing framework/src/timeintegrators/ActuallyExplicitEuler.C line=timeOld

However, `DirichletBC` derived boundary conditions need to use the **final** time to evaluate themselves. Think of it this way: you're integrating forward the "forces" as if evaluated from the beginning of the step... but ultimately the value on the boundary must end up being what it is supposed to be at the final time... no matter what. To achieve that we reset time to the `_current_time` in-between weak form evalution and `NodalBC` boundary condition application in `postResidual()`. `postResidual()` gets called at exactly this time to allow us to combine the `time` and `nontime` residuals into a single residual. So it's convenient to simply do:

!listing framework/src/timeintegrators/ActuallyExplicitEuler.C line=_fe_problem.time() = _current_time;

After `postResidual()` the `NodalBC` BCs will get applied with the time at the final time for the step.

### `meshChanged()`

When the mesh changes the linear solver needs to be destroyed and recreated. This is done by simply building a new one and setting it up again. This happens automatically just by "overwriting" the `std::unique_ptr` to the LinearSolver.

### `lump_preconditioned`

The `lump_preconditioned` option invokes a `LumpedPreconditioner` helper object:

!listing framework/src/timeintegrators/ActuallyExplicitEuler.C line=class LumpedPreconditioner

This helper object simply applies the inverse of the diagonal, lumped mass-matrix as the preconditioner for the linear solve. This is extremely efficient. Note that when this option is applied you shouldn't specify any other preconditioners using command-line syntax or they will override this option. In my testing this worked well.
This time integrator is an implementation of forward/explicit Euler based on
[/ExplicitTimeIntegrator.md], which avoids a nonlinear solve.

!syntax parameters /Executioner/TimeIntegrator/ActuallyExplicitEuler

!syntax inputs /Executioner/TimeIntegrator/ActuallyExplicitEuler

!syntax children /Executioner/TimeIntegrator/ActuallyExplicitEuler

!bibtex bibliography
124 changes: 124 additions & 0 deletions framework/doc/content/source/timeintegrators/ExplicitTimeIntegrator.md
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,124 @@
# ExplicitTimeIntegrator

## Description

`ExplicitTimeIntegrator` is a base class for explicit time integrators that are
implemented without performing any nonlinear solve, which reduces runtime. Unlike
explicit time integrators that are not derived from this class, it is not
necessary to set `implicit` to `false` for all of the non-time residual objects.

## Methods of Solution

Time integrators deriving from this class have three solve options, provided via
the `solve_type` parameter:

- `consistent`: (the default) A full mass matrix is built and used in a linear solve for the update
- `lumped`: A "lumped" mass matrix is built, inverted, and applied to the RHS, which
is faster but can possibly be less accurate.
- `lump_preconditioned`: The inversion of the "lumped" mass matrix is used to
precondition the `consistent` solve.

All three methods are solved similarly: a linear solve is performed to obtain a
solution update $\delta u$ that is added to the existing solution.

Below is some more explanation of each of these `solve_type` options:

### `consistent`

The `consistent` option builds a full ("consistent") "mass matrix" and uses it
in a linear solve to get the update. This is done by calling
`FEProblem::computeJacobianTag()` and specifying the `TIME` tag which includes
all of the `TimeKernel` derived Kernels and `NodalBC` derived BoundaryConditions
to compute $\mathbf{M}$:

!listing framework/src/timeintegrators/ActuallyExplicitEuler.C line=computeJacobianTag

A residual computation is also completed to use as the RHS ($R$):

!listing framework/src/timeintegrators/ActuallyExplicitEuler.C line=computeResidual

Finally, the following linear system is solved to obtain the solution update
$\delta u$ using the default linear solver from libMesh (usually PETSc, including
the application of any command-line parameters specified):

\begin{equation}
\mathbf{M} \delta u = -R
\end{equation}

### `lumped`

The `lumped` option creates a "lumped" mass matrix to use in the solve. A lumped mass matrix is a diagonal matrix where the diagonal is the sum of all elements on the row from the original matrix.
Here, to achieve the lumping, a matrix-vector product is performed between the
consistent mass matrix and a vector of all ones.

The inverse of a diagonal matrix is simply the reciprocal of each diagonal entry.
Then the matrix-vector product of the "inverse" lumped diagonal matrix is applied by simply doing a pointwise multiplication with the RHS.

Thus the `lumped` option does not actually solve a system of linear equations,
allowing it be much faster. However, the lumping of the mass matrix may lead to
unacceptable phase errors.

### `lump_preconditioned`

This option is the combination of the other two options: the consistent mass matrix
is used in the linear system, but the linear solve is preconditioned using the
lumped mass matrix. This compromise retains the accuracy of the `consistent`
option while benefiting from some of the speedup offered by the `lumped` option.

The lumped mass matrix preconditioner is applied with a helper class, `LumpedPreconditioner`:

!listing framework/include/timeintegrators/ExplicitTimeIntegrator.h line=class LumpedPreconditioner

This helper object simply applies the inverse of the diagonal, lumped
mass-matrix as the preconditioner for the linear solve, which is very
efficient. Note that when this option is applied you shouldn't specify any other
preconditioners using command-line syntax or they will override this option.

## Additional Details

Some notes on some of the implementation details of this class follow.

### Update Form

Note that even though we're doing an explicit solve we are currently doing it in
"update form" similar to a single step Newton solve. This gives us good parity
with the rest of MOOSE. We may want to change this in the future to make better
use of the fact that the mass-matrix can be constant for a wider class of
problems if we remove `dt` from it.

### `_ones`

To get the sum of each row of the mass matrix for "lumping" purposes a vector consisting of all `1`s is used in a matrix-vector product:

!listing framework/src/timeintegrators/ActuallyExplicitEuler.C line=mass_matrix.vector_mult

This is actually the very same way `MatGetRowSum` is implemented in PETSc; however,
doing it manually cuts down on vector creation/destruction and a few other book-keeping operations.

In the future this could be changed to use `MatGetRowSum` if a specialization for `MPI_Aij` format is created.

### Time

Time in an explicit solve must be handled carefully. When evaluating the weak
form (the integral parts) of the residual, time needs to be set to be the "old"
time (the time we are solving "from"):

!listing framework/src/timeintegrators/ActuallyExplicitEuler.C line=timeOld

However, `DirichletBC` derived boundary conditions need to use the **final** time,
since the strong constraints they represent use the final time and are not affected
by the time integrator. To achieve this, time is reset to the `_current_time` after
the weak form residual evaluation and before `NodalBC` boundary condition application,
which makes `postResidual()` the correct place to reset time for this purpose:

!listing framework/src/timeintegrators/ActuallyExplicitEuler.C line=_fe_problem.time() = _current_time;

After `postResidual()` the `NodalBC` BCs are applied with the time at the final time for the step.

### `meshChanged()`

When the mesh changes the linear solver needs to be destroyed and recreated.
This is done by simply building a new one and setting it up again. This happens
automatically just by "overwriting" the `std::unique_ptr` to the LinearSolver.

!syntax children /Executioner/TimeIntegrator/ExplicitTimeIntegrator
53 changes: 2 additions & 51 deletions framework/include/timeintegrators/ActuallyExplicitEuler.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,14 +9,9 @@

#pragma once

#include "TimeIntegrator.h"
#include "MeshChangedInterface.h"
#include "ExplicitTimeIntegrator.h"

#include "libmesh/linear_solver.h"

// Forward declarations
class ActuallyExplicitEuler;
class LumpedPreconditioner;

template <>
InputParameters validParams<ActuallyExplicitEuler>();
Expand All @@ -25,66 +20,23 @@ InputParameters validParams<ActuallyExplicitEuler>();
* Implements a truly explicit (no nonlinear solve) first-order, forward Euler
* time integration scheme.
*/
class ActuallyExplicitEuler : public TimeIntegrator, public MeshChangedInterface
class ActuallyExplicitEuler : public ExplicitTimeIntegrator
{
public:
ActuallyExplicitEuler(const InputParameters & parameters);

virtual void initialSetup() override;
virtual void init() override;
virtual void preSolve() override;
virtual int order() override { return 1; }
virtual void computeTimeDerivatives() override;
void computeADTimeDerivatives(DualReal & ad_u_dot, const dof_id_type & dof) const override;
virtual void solve() override;
virtual void postResidual(NumericVector<Number> & residual) override;

virtual void meshChanged() override;

protected:
enum SolveType
{
CONSISTENT,
LUMPED,
LUMP_PRECONDITIONED
};

/**
* Check for the linear solver convergence
*/
bool checkLinearConvergence();

/**
* Helper function that actually does the math for computing the time derivative
*/
template <typename T, typename T2>
void computeTimeDerivativeHelper(T & u_dot, const T2 & u_old) const;

MooseEnum _solve_type;

/// Residual used for the RHS
NumericVector<Real> & _explicit_residual;

/// Solution vector for the linear solve
NumericVector<Real> & _explicit_euler_update;

/// Diagonal of the lumped mass matrix (and its inversion)
NumericVector<Real> & _mass_matrix_diag;

/// Just a vector of 1's to help with creating the lumped mass matrix
NumericVector<Real> * _ones;

/// For computing the mass matrix
TagID _Ke_time_tag;

/// For solving with the consistent matrix
std::unique_ptr<LinearSolver<Number>> _linear_solver;

/// For solving with lumped preconditioning
std::unique_ptr<LumpedPreconditioner> _preconditioner;

/// Save off current time to reset it back and forth
Real _current_time;
};

template <typename T, typename T2>
Expand All @@ -94,4 +46,3 @@ ActuallyExplicitEuler::computeTimeDerivativeHelper(T & u_dot, const T2 & u_old)
u_dot -= u_old;
u_dot *= 1. / _dt;
}

0 comments on commit 78704e2

Please sign in to comment.