Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <cstddef>
       7             : 
       8             : #include "DataStructures/DataBox/Prefixes.hpp"
       9             : #include "IO/Observer/Actions/RegisterWithObservers.hpp"
      10             : #include "IO/Observer/Helpers.hpp"
      11             : #include "ParallelAlgorithms/Actions/Goto.hpp"
      12             : #include "ParallelAlgorithms/LinearSolver/AsynchronousSolvers/ElementActions.hpp"
      13             : #include "ParallelAlgorithms/LinearSolver/Multigrid/ElementActions.hpp"
      14             : #include "ParallelAlgorithms/LinearSolver/Multigrid/ObserveVolumeData.hpp"
      15             : #include "Utilities/TMPL.hpp"
      16             : 
      17             : /// Items related to the multigrid linear solver
      18             : ///
      19             : /// \see `LinearSolver::multigrid::Multigrid`
      20             : namespace LinearSolver::multigrid {
      21             : 
      22             : /// A label indicating the pre-smoothing step in a V-cycle multigrid algorithm,
      23             : /// i.e. the smoothing step before sending the residual to the coarser (parent)
      24             : /// grid
      25           1 : struct VcycleDownLabel {};
      26             : 
      27             : /// A label indicating the post-smoothing step in a V-cycle multigrid algorithm,
      28             : /// i.e. the smoothing step before sending the correction to the finer (child)
      29             : /// grid
      30           1 : struct VcycleUpLabel {};
      31             : 
      32             : /*!
      33             :  * \brief A V-cycle geometric multgrid solver for linear equations \f$Ax = b\f$
      34             :  *
      35             :  * This linear solver iteratively corrects an initial guess \f$x_0\f$ by
      36             :  * restricting the residual \f$b - Ax\f$ to a series of coarser grids, solving
      37             :  * for a correction on the coarser grids, and then prolongating (interpolating)
      38             :  * the correction back to the finer grids. The solves on grids with different
      39             :  * scales can very effectively solve large-scale modes in the solution, which
      40             :  * Krylov-type linear solvers such as GMRES or Conjugate Gradients typically
      41             :  * struggle with. Therefore, a multigrid solver can be an effective
      42             :  * preconditioner for Krylov-type linear solvers (see
      43             :  * `LinearSolver::gmres::Gmres` and `LinearSolver::cg::ConjugateGradient`). See
      44             :  * \cite Briggs2000jp for an introduction to multigrid methods.
      45             :  *
      46             :  * \par Grid hierarchy
      47             :  * This geometric multigrid solver relies on a strategy to coarsen the
      48             :  * computational grid in a way that removes small-scale modes. We currently
      49             :  * h-coarsen the initial domain, meaning that we create multigrid levels by
      50             :  * successively combining two elements into one along every dimension of the
      51             :  * grid. We only p-coarsen the grid in the sense that we choose the smaller of
      52             :  * the two polynomial degrees when combining elements. This strategy follows
      53             :  * \cite Vincent2019qpd. See `LinearSolver::multigrid::ElementsAllocator` and
      54             :  * `LinearSolver::multigrid::coarsen` for the code that creates the initial
      55             :  * multigrid hierarchy.
      56             :  * Once the initial grids are created, AMR steps can create incrementally finer
      57             :  * grids. For this to work, AMR must be configured with the compile-time option
      58             :  * `keep_coarse_grids = true` (see `amr::protocols::AmrMetavariables`).
      59             :  *
      60             :  * \par Inter-mesh operators
      61             :  * The algorithm relies on operations that project data between grids. Residuals
      62             :  * are projected from finer to coarser grids ("restriction") and solutions are
      63             :  * projected from coarser to finer grids ("prolongation"). We use the standard
      64             :  * \f$L_2\f$-projections implemented in
      65             :  * `Spectral::projection_matrix_child_to_parent` and
      66             :  * `Spectral::projection_matrix_parent_to_child`, where "child" means the finer
      67             :  * grid and "parent" means the coarser grid. Note that the residuals \f$b -
      68             :  * Ax\f$ may or may not already include a mass matrix and Jacobian factors,
      69             :  * depending on the implementation of the linear operator \f$A\f$. To account
      70             :  * for this, provide a tag as the `ResidualIsMassiveTag` that holds a `bool`.
      71             :  * See section 3 in \cite Fortunato2019jl for details on the inter-mesh
      72             :  * operators.
      73             :  *
      74             :  * \par Smoother and bottom solver
      75             :  * On every level of the grid hierarchy the multigrid solver relies on a
      76             :  * "smoother" that performs an approximate linear solve on that level. The
      77             :  * smoother can be any linear solver that solves the `smooth_fields_tag` for the
      78             :  * source in the `smooth_source_tag`. Useful smoothers are asynchronous linear
      79             :  * solvers that parallelize well, such as `LinearSolver::Schwarz::Schwarz`. The
      80             :  * multigrid algorithm doesn't assume anything about the smoother, but requires
      81             :  * only that the `PreSmootherActions` and the `PostSmootherActions` passed to
      82             :  * `solve` leave the `smooth_fields_tag` in a state that represents an
      83             :  * approximate solution to the `smooth_source_tag`, and that
      84             :  * `db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo,
      85             :  * smooth_fields_tag>` is left up-to-date as well.
      86             :  * The smoother can assume that, on entry, these two tags represent the initial
      87             :  * fields and the linear operator applied to the initial fields, respectively.
      88             :  * Here's an example of setting up a smoother for the multigrid solver:
      89             :  *
      90             :  * \snippet Test_MultigridAlgorithm.cpp setup_smoother
      91             :  *
      92             :  * On the coarsest level of the grid hierarchy, the multigrid algorithm can also
      93             :  * use a separate "bottom solver". If enabled, the bottom solver replaces the
      94             :  * pre-smoother on the coarsest grid. The bottom solver can also be any linear
      95             :  * solver, but it is typically a direct solver that solves the linear problem
      96             :  * exactly (see `LinearSolver::Actions::BuildMatrix` for a good bottom solver).
      97             :  *
      98             :  * The smoother and bottom solver can be used to construct an action list like
      99             :  * this:
     100             :  *
     101             :  * \snippet Test_MultigridAlgorithm.cpp action_list
     102             :  *
     103             :  * \par Algorithm overview
     104             :  * Every iteration of the multigrid algorithm performs a V-cycle over the grid
     105             :  * hierarchy. One V-cycle consists of first "going down" the grid hierarchy,
     106             :  * from the finest to successively coarser grids, smoothing on every level
     107             :  * ("pre-smoothing", can be controlled by the option
     108             :  * `LinearSolver::multigrid::Tags::EnablePreSmoothing`),
     109             :  * and then "going up" the grid hierarchy again, smoothing on
     110             :  * every level again ("post-smoothing"). When going down, the algorithm projects
     111             :  * the remaining residual of the smoother to the next-coarser grid, setting it
     112             :  * as the source for the smoother on the coarser grid. When going up again, the
     113             :  * algorithm projects the solution of the smoother to the next-finer grid,
     114             :  * adding it to the solution on the finer grid as a correction. The bottom-most
     115             :  * coarsest grid (the "tip" of the V-cycle) may run the bottom solver instead of
     116             :  * the pre-smoother. It may also skip the post-smoother, so the result of the
     117             :  * bottom solver is immediately projected up to the finer grid (controlled by
     118             :  * the options `LinearSolver::multigrid::Tags::UseBottomSolver` and
     119             :  * `LinearSolver::multigrid::Tags::EnablePostSmoothingAtBottom`). On the
     120             :  * top-most finest grid (the "original" grid that represents the overall
     121             :  * solution) the algorithm applies the smoothing and the corrections from the
     122             :  * coarser grids directly to the solution fields.
     123             :  */
     124             : template <typename Metavariables, size_t Dim, typename FieldsTag,
     125             :           typename OptionsGroup, typename ResidualIsMassiveTag,
     126             :           typename SourceTag =
     127             :               db::add_tag_prefix<::Tags::FixedSource, FieldsTag>>
     128           1 : struct Multigrid {
     129           0 :   using fields_tag = FieldsTag;
     130           0 :   using options_group = OptionsGroup;
     131           0 :   using source_tag = SourceTag;
     132             : 
     133           0 :   using operand_tag = FieldsTag;
     134             : 
     135           0 :   using smooth_source_tag = source_tag;
     136           0 :   using smooth_fields_tag = fields_tag;
     137             : 
     138           0 :   using component_list = tmpl::list<>;
     139             : 
     140           0 :   using observed_reduction_data_tags = observers::make_reduction_data_tags<
     141             :       tmpl::list<async_solvers::reduction_data>>;
     142             : 
     143           0 :   using initialize_element = tmpl::list<
     144             :       async_solvers::InitializeElement<FieldsTag, OptionsGroup, SourceTag>,
     145             :       detail::InitializeElement<Dim, FieldsTag, OptionsGroup, SourceTag>>;
     146             : 
     147           0 :   using amr_projectors = initialize_element;
     148             : 
     149           0 :   using register_element = tmpl::list<
     150             :       async_solvers::RegisterElement<FieldsTag, OptionsGroup, SourceTag,
     151             :                                      ::amr::Tags::IsFinestGrid>,
     152             :       observers::Actions::RegisterWithObservers<
     153             :           detail::RegisterWithVolumeObserver<OptionsGroup>>>;
     154             : 
     155             :   template <typename ApplyOperatorActions, typename PreSmootherActions,
     156             :             typename PostSmootherActions,
     157             :             typename BottomSolverActions = tmpl::list<>,
     158             :             typename Label = OptionsGroup>
     159           0 :   using solve = tmpl::list<
     160             :       async_solvers::PrepareSolve<FieldsTag, OptionsGroup, SourceTag, Label,
     161             :                                   ::amr::Tags::IsFinestGrid, false>,
     162             :       detail::ReceiveResidualFromFinerGrid<Dim, FieldsTag, OptionsGroup,
     163             :                                            SourceTag>,
     164             :       detail::PreparePreSmoothing<
     165             :           FieldsTag, OptionsGroup, SourceTag,
     166             :           not std::is_same_v<BottomSolverActions, tmpl::list<>>>,
     167             :       // No need to apply the linear operator here:
     168             :       // - On the finest grid, the operator applied to the fields should have
     169             :       //   already been computed at this point, either applied to the initial
     170             :       //   fields before the multigrid solver is invoked, or by the smoother at
     171             :       //   the end of the previous V-cycle.
     172             :       // - On coarser grids, the initial fields are zero, so the operator
     173             :       //   applied to them is also zero.
     174             :       PreSmootherActions,
     175             :       detail::SkipBottomSolver<FieldsTag, OptionsGroup, SourceTag>,
     176             :       BottomSolverActions,
     177             :       detail::SkipPostSmoothingAtBottom<FieldsTag, OptionsGroup, SourceTag>,
     178             :       detail::SendResidualToCoarserGrid<FieldsTag, OptionsGroup,
     179             :                                         ResidualIsMassiveTag, SourceTag>,
     180             :       detail::ReceiveCorrectionFromCoarserGrid<Dim, FieldsTag, OptionsGroup,
     181             :                                                SourceTag>,
     182             :       ApplyOperatorActions, ::Actions::Label<detail::PostSmoothingBeginLabel>,
     183             :       PostSmootherActions,
     184             :       detail::SendCorrectionToFinerGrid<FieldsTag, OptionsGroup, SourceTag>,
     185             :       detail::ObserveVolumeData<FieldsTag, OptionsGroup, SourceTag>,
     186             :       async_solvers::CompleteStep<FieldsTag, OptionsGroup, SourceTag, Label,
     187             :                                   ::amr::Tags::IsFinestGrid, false>>;
     188             : };
     189             : 
     190             : }  // namespace LinearSolver::multigrid