template<typename Metavariables, size_t Dim, typename FieldsTag, typename OptionsGroup, typename ResidualIsMassiveTag, typename SourceTag = db::add_tag_prefix<::Tags::FixedSource, FieldsTag>>
struct LinearSolver::multigrid::Multigrid< Metavariables, Dim, FieldsTag, OptionsGroup, ResidualIsMassiveTag, SourceTag >
A V-cycle geometric multgrid solver for linear equations \(Ax = b\).
This linear solver iteratively corrects an initial guess \(x_0\) by restricting the residual \(b - Ax\) to a series of coarser grids, solving for a correction on the coarser grids, and then prolongating (interpolating) the correction back to the finer grids. The solves on grids with different scales can very effectively solve large-scale modes in the solution, which Krylov-type linear solvers such as GMRES or Conjugate Gradients typically struggle with. Therefore, a multigrid solver can be an effective preconditioner for Krylov-type linear solvers (see LinearSolver::gmres::Gmres and LinearSolver::cg::ConjugateGradient). See [33] for an introduction to multigrid methods.
- Grid hierarchy
- This geometric multigrid solver relies on a strategy to coarsen the computational grid in a way that removes small-scale modes. We currently h-coarsen the initial domain, meaning that we create multigrid levels by successively combining two elements into one along every dimension of the grid. We only p-coarsen the grid in the sense that we choose the smaller of the two polynomial degrees when combining elements. This strategy follows [208]. See LinearSolver::multigrid::ElementsAllocator and LinearSolver::multigrid::coarsen for the code that creates the initial multigrid hierarchy. Once the initial grids are created, AMR steps can create incrementally finer grids. For this to work, AMR must be configured with the compile-time option keep_coarse_grids = true (see amr::protocols::AmrMetavariables).
- Inter-mesh operators
- The algorithm relies on operations that project data between grids. Residuals are projected from finer to coarser grids ("restriction") and solutions are projected from coarser to finer grids ("prolongation"). We use the standard \(L_2\)-projections implemented in Spectral::projection_matrix_child_to_parent and Spectral::projection_matrix_parent_to_child, where "child" means the finer grid and "parent" means the coarser grid. Note that the residuals \(b -
Ax\) may or may not already include a mass matrix and Jacobian factors, depending on the implementation of the linear operator \(A\). To account for this, provide a tag as the ResidualIsMassiveTag that holds a bool. See section 3 in [78] for details on the inter-mesh operators.
- Smoother and bottom solver
- On every level of the grid hierarchy the multigrid solver relies on a "smoother" that performs an approximate linear solve on that level. The smoother can be any linear solver that solves the smooth_fields_tag for the source in the smooth_source_tag. Useful smoothers are asynchronous linear solvers that parallelize well, such as LinearSolver::Schwarz::Schwarz. The multigrid algorithm doesn't assume anything about the smoother, but requires only that the PreSmootherActions and the PostSmootherActions passed to solve leave the smooth_fields_tag in a state that represents an approximate solution to the smooth_source_tag, and that db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo,
smooth_fields_tag> is left up-to-date as well. The smoother can assume that, on entry, these two tags represent the initial fields and the linear operator applied to the initial fields, respectively. Here's an example of setting up a smoother for the multigrid solver:
using multigrid = LinearSolver::multigrid::Multigrid<
Metavariables, volume_dim, helpers_mg::fields_tag, MultigridSolver,
helpers_mg::OperatorIsMassive, helpers_mg::sources_tag>;
using smoother = LinearSolver::Richardson::Richardson<
typename multigrid::smooth_fields_tag, RichardsonSmoother,
typename multigrid::smooth_source_tag, ::amr::Tags::GridIndex>;
On the coarsest level of the grid hierarchy, the multigrid algorithm can also use a separate "bottom solver". If enabled, the bottom solver replaces the pre-smoother on the coarsest grid. The bottom solver can also be any linear solver, but it is typically a direct solver that solves the linear problem exactly (see LinearSolver::Actions::BuildMatrix for a good bottom solver).
The smoother and bottom solver can be used to construct an action list like this:
template <typename Label>
using smooth_actions = typename smoother::template solve<
compute_operator_action<typename smoother::operand_tag>, Label>;
using solve_actions =
tmpl::list<compute_operator_action<typename multigrid::fields_tag>,
typename multigrid::template solve<
compute_operator_action<typename smoother::fields_tag>,
smooth_actions<LinearSolver::multigrid::VcycleDownLabel>,
smooth_actions<LinearSolver::multigrid::VcycleUpLabel>,
bottom_solver_actions>,
Parallel::Actions::TerminatePhase>;
- Algorithm overview
- Every iteration of the multigrid algorithm performs a V-cycle over the grid hierarchy. One V-cycle consists of first "going down" the grid hierarchy, from the finest to successively coarser grids, smoothing on every level ("pre-smoothing", can be controlled by the option LinearSolver::multigrid::Tags::EnablePreSmoothing), and then "going up" the grid hierarchy again, smoothing on every level again ("post-smoothing"). When going down, the algorithm projects the remaining residual of the smoother to the next-coarser grid, setting it as the source for the smoother on the coarser grid. When going up again, the algorithm projects the solution of the smoother to the next-finer grid, adding it to the solution on the finer grid as a correction. The bottom-most coarsest grid (the "tip" of the V-cycle) may run the bottom solver instead of the pre-smoother. It may also skip the post-smoother, so the result of the bottom solver is immediately projected up to the finer grid (controlled by the options LinearSolver::multigrid::Tags::UseBottomSolver and LinearSolver::multigrid::Tags::EnablePostSmoothingAtBottom). On the top-most finest grid (the "original" grid that represents the overall solution) the algorithm applies the smoothing and the corrections from the coarser grids directly to the solution fields.