deal.II version GIT relicensing-2250-g88cb8ba3cb 2024-12-13 12:20:00+00:00
|
#include <deal.II/lac/petsc_ts.h>
Public Types | |
using | real_type = PetscReal |
Public Member Functions | |
TimeStepper (const TimeStepperData &data=TimeStepperData(), const MPI_Comm mpi_comm=PETSC_COMM_WORLD) | |
~TimeStepper () | |
operator TS () const | |
TS | petsc_ts () |
MPI_Comm | get_mpi_communicator () const |
void | reinit () |
void | reinit (const TimeStepperData &data) |
void | set_matrix (PMatrixType &P) |
void | set_matrices (AMatrixType &A, PMatrixType &P) |
real_type | get_time () |
real_type | get_time_step () |
unsigned int | get_step_number () |
unsigned int | solve (VectorType &y) |
unsigned int | solve (VectorType &y, PMatrixType &P) |
unsigned int | solve (VectorType &y, AMatrixType &A, PMatrixType &P) |
Public Attributes | |
std::function< void(const real_type t, const VectorType &y, const VectorType &y_dot, VectorType &res)> | implicit_function |
std::function< void(const real_type t, const VectorType &y, const VectorType &y_dot, const real_type alpha, AMatrixType &A, PMatrixType &P)> | implicit_jacobian |
std::function< void(const real_type t, const VectorType &y, VectorType &res)> | explicit_function |
std::function< void(const real_type t, const VectorType &y, AMatrixType &A, PMatrixType &P)> | explicit_jacobian |
std::function< void(const real_type t, const VectorType &y, const unsigned int step_number)> | monitor |
std::function< void(const real_type t, const VectorType &y, const VectorType &ydot, const real_type alpha)> | setup_jacobian |
std::function< void(const VectorType &src, VectorType &dst)> | solve_with_jacobian |
std::function< IndexSet()> | algebraic_components |
std::function< void(const real_type t, VectorType &y)> | distribute |
std::function< void(const real_type t, VectorType &y)> | update_constrained_components |
std::function< void(const real_type t, const unsigned int step, const VectorType &y, bool &resize)> | decide_for_coarsening_and_refinement |
std::function< bool(const real_type t, const unsigned int step, const VectorType &y)> | decide_and_prepare_for_remeshing |
std::function< void(const std::vector< VectorType > &all_in, std::vector< VectorType > &all_out)> | interpolate |
std::function< void(const real_type t, const std::vector< VectorType > &all_in, std::vector< VectorType > &all_out)> | transfer_solution_vectors_to_new_mesh |
Private Member Functions | |
void | setup_callbacks () |
void | setup_algebraic_constraints (const VectorType &y) |
Private Attributes | |
TS | ts |
ObserverPointer< AMatrixType, TimeStepper > | A |
ObserverPointer< PMatrixType, TimeStepper > | P |
PreconditionShell | solve_with_jacobian_pc |
bool | restart_if_remesh |
bool | need_dae_tolerances |
bool | need_dummy_assemble |
std::exception_ptr | pending_exception |
bool | error_in_function |
Interface to the PETSc TS solver for Ordinary Differential Equations and Differential-Algebraic Equations. The TS solver is described in the PETSc manual. This class is used and extensively discussed in step-86.
This class supports two kinds of formulations. The explicit formulation:
\[ \begin{cases} \dot y = G(t,y)\, , \\ y(t_0) = y_0\, , \\ \end{cases} \]
and the implicit formulation:
\[ \begin{cases} F(t,y,\dot y) = 0\, , \\ y(t_0) = y_0\, . \\ \end{cases} \]
The interface to PETSc is realized by means of std::function callbacks like in the SUNDIALS::IDA (which also solves implicit ODES) and SUNDIALS::ARKode classes (which solves a slightly generalized form of the explicit formulation above that also allows for a mass matrix on the left hand side).
TimeStepper supports any vector and matrix type having constructors and methods:
In particular, the supported types are the ones that can wrap PETSc's native vector and matrix classes, that are able to modify them in place, and that can return PETSc native types when requested.
To use explicit solvers (like for example explicit Runge-Kutta methods), the user only needs to provide the implementation of \(G\) via the TimeStepper::explicit_function. For implicit solvers, users have also the alternative of providing the \(F\) function via TimeStepper::implicit_function. IMEX methods are also supported by providing both callbacks.
The default linearization procedure of an implicit solver instantiated with this class consists in using Jacobian-Free-Newton-Krylov; the action of tangent matrices inside a linear solver process are approximated via matrix-free finite-differencing of the nonlinear residual equations that are ODE-solver specific. For details, consult the PETSc manual.
Users can also provide the implementations of the Jacobians. This can be accomplished in two ways:
TimeStepper::set_matrices() must be used in case the user wants to provide the iteration matrix of the tangent system in the deal.II style approach, thus replacing the matrix-free linearization.
The correctness of the constructed Jacobians passed using TimeStepper::set_matrix() can be checked using
See TimeStepper::set_matrix() and TimeStepper::set_matrices() for additional details.
The deal.II style approach still allows command line customization, like for example,
in case the user wants to change the default nonlinear solver to a trust region solver and iterate on the tangent system with CG, still using TimeStepper::solve_with_jacobian as a preconditioner.
The PETSc style approach has instead the advantage that only the matrix assembly procedure has to be implemented, thus allowing quicker implementations and faster turnaround for experimenting with linear solver preconditioning configurations via command line customizations, like for example,
Definition at line 338 of file petsc_ts.h.
using PETScWrappers::TimeStepper< VectorType, PMatrixType, AMatrixType >::real_type = PetscReal |
Type that holds real-valued numbers.
Used to represent time and norms tolerances.
Definition at line 346 of file petsc_ts.h.
PETScWrappers::TimeStepper< VectorType, PMatrixType, AMatrixType >::TimeStepper | ( | const TimeStepperData & | data = TimeStepperData() , |
const MPI_Comm | mpi_comm = PETSC_COMM_WORLD |
||
) |
Constructor.
PETScWrappers::TimeStepper< VectorType, PMatrixType, AMatrixType >::~TimeStepper | ( | ) |
Destructor.
PETScWrappers::TimeStepper< VectorType, PMatrixType, AMatrixType >::operator TS | ( | ) | const |
Conversion operator to gain access to the underlying PETSc type. If you do this, you cut this class off some information it may need, so this conversion operator should only be used if you know what you do.
TS PETScWrappers::TimeStepper< VectorType, PMatrixType, AMatrixType >::petsc_ts | ( | ) |
Return the PETSc TS object.
MPI_Comm PETScWrappers::TimeStepper< VectorType, PMatrixType, AMatrixType >::get_mpi_communicator | ( | ) | const |
Return the underlying MPI communicator.
void PETScWrappers::TimeStepper< VectorType, PMatrixType, AMatrixType >::reinit | ( | ) |
Reset the solver, it does not change the customization.
void PETScWrappers::TimeStepper< VectorType, PMatrixType, AMatrixType >::reinit | ( | const TimeStepperData & | data | ) |
Reset solver. Change customization according to data
.
void PETScWrappers::TimeStepper< VectorType, PMatrixType, AMatrixType >::set_matrix | ( | PMatrixType & | P | ) |
Set the preconditioning matrix only.
When used with TimeStepper::setup_jacobian and TimeStepper::solve_with_jacobian, PETSc will approximate the linear system matrix-vector product using an internal matrix-free representation.
When used with TimeStepper::implicit_jacobian or TimeStepper::explicit_jacobian, PETSc will use the same matrix for both preconditioning and matrix-vector products.
void PETScWrappers::TimeStepper< VectorType, PMatrixType, AMatrixType >::set_matrices | ( | AMatrixType & | A, |
PMatrixType & | P | ||
) |
Set both the linear system matrix and the preconditioning matrix that PETSc will use (can be the same matrix). In this case, the Jacobian-Free-Newton-Krylov approach will not be used.
real_type PETScWrappers::TimeStepper< VectorType, PMatrixType, AMatrixType >::get_time | ( | ) |
Return current time.
real_type PETScWrappers::TimeStepper< VectorType, PMatrixType, AMatrixType >::get_time_step | ( | ) |
Return current time step.
unsigned int PETScWrappers::TimeStepper< VectorType, PMatrixType, AMatrixType >::get_step_number | ( | ) |
Return current step number.
unsigned int PETScWrappers::TimeStepper< VectorType, PMatrixType, AMatrixType >::solve | ( | VectorType & | y | ) |
Integrate the differential-algebraic equations starting from y
.
This function returns the final number of computed steps. Upon returning, the y
vector contains the solution of the DAE at the end time.
unsigned int PETScWrappers::TimeStepper< VectorType, PMatrixType, AMatrixType >::solve | ( | VectorType & | y, |
PMatrixType & | P | ||
) |
Integrate the differential-algebraic equations starting from y
.
This function returns the final number of computed steps. Upon returning, the y
vector contains the solution of the DAE at the end time.
Here we also set the matrix to precondition the tangent system.
unsigned int PETScWrappers::TimeStepper< VectorType, PMatrixType, AMatrixType >::solve | ( | VectorType & | y, |
AMatrixType & | A, | ||
PMatrixType & | P | ||
) |
Integrate the differential-algebraic equations starting from y
.
This function returns the final number of computed steps. Upon returning, the y
vector contains the solution of the DAE at the end time.
Here we also set the matrices to describe and precondition the tangent system.
|
private |
Setup callbacks.
This function is called inside TimeStepper::solve routines and does not need to be called by the user. It is used to reinitialize the solver if mesh adaption has been performed.
|
private |
Setup algebraic constraints.
This function is called inside TimeStepper::solve routines and does not need to be called by the user. It is used to reinitialize the solver if mesh adaption has been performed.
std::function<void(const real_type t, const VectorType &y, const VectorType &y_dot, VectorType &res)> PETScWrappers::TimeStepper< VectorType, PMatrixType, AMatrixType >::implicit_function |
Callback for the computation of the implicit residual \(F(t, y, \dot y)\).
Definition at line 479 of file petsc_ts.h.
std::function<void(const real_type t, const VectorType &y, const VectorType &y_dot, const real_type alpha, AMatrixType &A, PMatrixType &P)> PETScWrappers::TimeStepper< VectorType, PMatrixType, AMatrixType >::implicit_jacobian |
Callback for the computation of the implicit Jacobian \(\dfrac{\partial F}{\partial y} + \alpha \dfrac{\partial F}{\partial \dot y}\).
All implicit solvers implementations are recast to use the above linearization. The \(\alpha\) parameter is time-step and solver-type specific.
Definition at line 501 of file petsc_ts.h.
std::function<void(const real_type t, const VectorType &y, VectorType &res)> PETScWrappers::TimeStepper< VectorType, PMatrixType, AMatrixType >::explicit_function |
Callback for the computation of the explicit residual \(G(t, y)\).
Definition at line 512 of file petsc_ts.h.
std::function<void(const real_type t, const VectorType &y, AMatrixType &A, PMatrixType &P)> PETScWrappers::TimeStepper< VectorType, PMatrixType, AMatrixType >::explicit_jacobian |
Callback for the computation of the explicit Jacobian \(\dfrac{\partial G}{\partial y}\).
Definition at line 527 of file petsc_ts.h.
std::function<void(const real_type t, const VectorType &y, const unsigned int step_number)> PETScWrappers::TimeStepper< VectorType, PMatrixType, AMatrixType >::monitor |
Callback for monitoring the solution process.
This function is called by TimeStepper at the beginning of each time step.
Definition at line 543 of file petsc_ts.h.
std::function<void(const real_type t, const VectorType &y, const VectorType &ydot, const real_type alpha)> PETScWrappers::TimeStepper< VectorType, PMatrixType, AMatrixType >::setup_jacobian |
Callback for the set up of the Jacobian system.
This callback gives full control to users to set up the linearized equations \(\dfrac{\partial F}{\partial y} + \alpha \dfrac{\partial F}{\partial \dot y}\).
All implicit solvers implementations are recast to use the above linearization. The \(\alpha\) parameter is time-step and solver-type specific.
Solvers must be provided via TimeStepper::solve_with_jacobian.
Definition at line 568 of file petsc_ts.h.
std::function<void(const VectorType &src, VectorType &dst)> PETScWrappers::TimeStepper< VectorType, PMatrixType, AMatrixType >::solve_with_jacobian |
Callback for the solution of the tangent system set up with TimeStepper::setup_jacobian.
This is used as a preconditioner inside the Krylov process.
Definition at line 582 of file petsc_ts.h.
std::function<IndexSet()> PETScWrappers::TimeStepper< VectorType, PMatrixType, AMatrixType >::algebraic_components |
Callback to return an index set containing the algebraic components. Algebraic components are degrees of freedom for which the differential equation does not provide a time derivative. This can either be because the degree of freedom is constrained (say, because it is a hanging node, or because it is part of Dirichlet boundary values), or because the differential equation simply does not contain time derivatives for a specific solution variable. An example for the latter case is the pressure in the time dependent Stokes equations,
\begin{align*} \frac{\partial \mathbf u(\mathbf x,t)}{\partial t} - \nu \Delta \mathbf u(\mathbf x,t) + \nabla p(\mathbf x,t) &= \mathbf f(\mathbf x,t), \\ \nabla \cdot \mathbf u(\mathbf x,t) &= 0. \end{align*}
The documentation of the SUNDIALS::IDA class has an extensive documentation of algebraic variables as part of differential-algebraic equations.
Implementation of this function is optional. If your equation is also algebraic (i.e., it contains algebraic constraints, or Lagrange multipliers), you should implement this function in order to return only these components of your system.
Definition at line 614 of file petsc_ts.h.
std::function<void(const real_type t, VectorType &y)> PETScWrappers::TimeStepper< VectorType, PMatrixType, AMatrixType >::distribute |
update_constrained_components
, but is deprecated. Use update_constrained_components
instead. Definition at line 621 of file petsc_ts.h.
std::function<void(const real_type t, VectorType &y)> PETScWrappers::TimeStepper< VectorType, PMatrixType, AMatrixType >::update_constrained_components |
Callback to set the values of constrained components to their correct values. Constrained components are a subset of the algebraic components discussed in the documentation of the TimeStepper::algebraic_components callback. In practice, the constrained components are typically hanging nodes and degrees of freedom constrained by Dirichlet boundary conditions.
Implementation of this function is optional. It is called at the end of each successful stage. The same functionality can be equivalently implemented in TimeStepper::solve_with_jacobian.
Definition at line 642 of file petsc_ts.h.
std::function<void(const real_type t, const unsigned int step, const VectorType &y, bool &resize)> PETScWrappers::TimeStepper< VectorType, PMatrixType, AMatrixType >::decide_for_coarsening_and_refinement |
decide_and_prepare_for_remeshing
except that it returns the decision whether or not to stop operations via the last reference argument of the function object instead of a plain return value. This callback is deprecated. Use decide_and_prepare_for_remeshing
instead. Definition at line 656 of file petsc_ts.h.
std::function< bool(const real_type t, const unsigned int step, const VectorType &y)> PETScWrappers::TimeStepper< VectorType, PMatrixType, AMatrixType >::decide_and_prepare_for_remeshing |
A callback that returns whether or not to stop time stepping at the current moment for mesh adaptation. If the callback returns true
, then the time stepper stops time integration for now, saves some state, calls the interpolate
callback, and then resumes time integration. Either in the current callback or the interpolate
callback, the user code needs to perform the mesh refinement.
Implementation of this function is optional. The callback must return true
if mesh adaption is to be performed, false
otherwise. The y
vector contains the current solution, t
the current time, @ step the step number.
Definition at line 679 of file petsc_ts.h.
std::function<void(const std::vector<VectorType> &all_in, std::vector<VectorType> &all_out)> PETScWrappers::TimeStepper< VectorType, PMatrixType, AMatrixType >::interpolate |
transfer_solution_vectors_to_new_mesh
, but is deprecated. Use transfer_solution_vectors_to_new_mesh
instead. Definition at line 688 of file petsc_ts.h.
std::function<void(const real_type t, const std::vector<VectorType> &all_in, std::vector<VectorType> &all_out)> PETScWrappers::TimeStepper< VectorType, PMatrixType, AMatrixType >::transfer_solution_vectors_to_new_mesh |
Callback to perform mesh adaptation and transfer solution vectors from the old to the new mesh.
Implementation of this function is mandatory if TimeStepper::decide_and_prepare_for_remeshing is used. This function must perform mesh adaption and interpolate the discrete functions that are stored in all_in
onto the refined and/or coarsened grid. The output vectors must be sized correctly within this callback.
Definition at line 708 of file petsc_ts.h.
|
private |
The PETSc object.
Definition at line 714 of file petsc_ts.h.
|
private |
Pointers to the internal PETSc matrix objects.
Definition at line 719 of file petsc_ts.h.
|
private |
Definition at line 720 of file petsc_ts.h.
|
private |
Object to apply solve_with_jacobian.
Definition at line 725 of file petsc_ts.h.
|
private |
This flag is used to decide whether or not to restart the step if remeshing has been performed
Definition at line 731 of file petsc_ts.h.
|
private |
This flag is set when changing the customization and used within solve.
Definition at line 736 of file petsc_ts.h.
|
private |
This flag is used to support versions of PETSc older than 3.13.
Definition at line 741 of file petsc_ts.h.
|
mutableprivate |
A pointer to any exception that may have been thrown in user-defined call-backs and that we have to deal after the KINSOL function we call has returned.
Definition at line 748 of file petsc_ts.h.
|
private |
Internal data to handle recoverable errors.
Definition at line 753 of file petsc_ts.h.