Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include "IO/Observer/Helpers.hpp"
       7             : #include "ParallelAlgorithms/LinearSolver/AsynchronousSolvers/ElementActions.hpp"
       8             : #include "ParallelAlgorithms/LinearSolver/Schwarz/ElementActions.hpp"
       9             : #include "Utilities/ProtocolHelpers.hpp"
      10             : #include "Utilities/TMPL.hpp"
      11             : 
      12             : /// Items related to the %Schwarz linear solver
      13             : ///
      14             : /// \see `LinearSolver::Schwarz::Schwarz`
      15             : namespace LinearSolver::Schwarz {
      16             : 
      17             : /*!
      18             :  * \ingroup LinearSolverGroup
      19             :  * \brief An additive Schwarz subdomain solver for linear systems of equations
      20             :  * \f$Ax=b\f$.
      21             :  *
      22             :  * This Schwarz-type linear solver works by solving many sub-problems in
      23             :  * parallel and combining their solutions as a weighted sum to converge towards
      24             :  * the global solution. Each sub-problem is the restriction of the global
      25             :  * problem to an element-centered subdomain, which consists of a central element
      26             :  * and an overlap region with its neighbors. The decomposition into independent
      27             :  * sub-problems makes this linear solver very parallelizable, avoiding global
      28             :  * synchronization points altogether. It is commonly used as a preconditioner to
      29             :  * Krylov-type solvers such as `LinearSolver::gmres::Gmres` or
      30             :  * `LinearSolver::cg::ConjugateGradient`, or as the "smoother" for a Multigrid
      31             :  * solver (which may in turn precondition a Krylov-type solver).
      32             :  *
      33             :  * This linear solver relies on an implementation of the global linear operator
      34             :  * \f$A(x)\f$ and its restriction to a subdomain \f$A_{ss}(x)\f$. Each step of
      35             :  * the algorithm expects that \f$A(x)\f$ is computed and stored in the DataBox
      36             :  * as `db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo, operand_tag>`.
      37             :  * To perform a solve, add the `solve` action list to an array parallel
      38             :  * component. Pass the actions that compute \f$A(x)\f$, as well as any further
      39             :  * actions you wish to run in each step of the algorithm, as the first template
      40             :  * parameter to `solve`. If you add the `solve` action list multiple times, use
      41             :  * the second template parameter to label each solve with a different type.
      42             :  * Pass the subdomain operator as the `SubdomainOperator` template parameter to
      43             :  * this class (not to the `solve` action list template). See the paragraph below
      44             :  * for information on implementing a subdomain operator.
      45             :  *
      46             :  * \par Subdomain geometry:
      47             :  * The image below gives an overview of the structure of an element-centered
      48             :  * subdomain:
      49             :  *
      50             :  * \image html subdomain_structure.svg "Fig. 1: An element-centered subdomain" width=600px
      51             :  *
      52             :  * Fig. 1 shows part of a 2-dimensional computational domain. The domain is
      53             :  * composed of elements (light gray) that each carry a Legendre-Gauss-Lobatto
      54             :  * mesh of grid points (black dots). The diagonally shaded region to the left
      55             :  * illustrates an external domain boundary. The lines between neighboring
      56             :  * element faces illustrate mortar meshes, which are relevant when implementing
      57             :  * a subdomain operator in a DG context but play no role in the Schwarz
      58             :  * algorithm. The dashed line gives an example of an element-centered subdomain
      59             :  * with a maximum of 2 grid points of overlap into neighboring elements. For
      60             :  * details on how the number of overlap points is chosen see
      61             :  * `LinearSolver::Schwarz::overlap_extent`. Note that the subdomain does not
      62             :  * extend into corner- or edge-neighbors, which is different to both
      63             :  * \cite Stiller2016a and \cite Vincent2019qpd. The reason we don't include
      64             :  * such diagonal couplings is that in a DG context information only propagates
      65             :  * across faces, as noted in \cite Stiller2016a. By eliminating the corner-
      66             :  * and edge-neighbors we significantly reduce the complexity of the subdomain
      67             :  * geometry and potentially also the communication costs. However, the
      68             :  * advantages and disadvantages of this choice have yet to be carefully
      69             :  * evaluated.
      70             :  *
      71             :  * \par Subdomain operator:
      72             :  * The Schwarz subdomain solver relies on a restriction of the global linear
      73             :  * problem \f$Ax=b\f$ to the subdomains. The subdomain operator
      74             :  * \f$A_{ss}=R_s A R_s^T\f$, where \f$R_s\f$ is the restriction operator,
      75             :  * essentially assumes that its operand is zero everywhere but within the
      76             :  * subdomain (see Section 3.1 in \cite Stiller2016a). Therefore it can be
      77             :  * evaluated entirely on an element-centered subdomain with no need to
      78             :  * communicate with neighbors within each subdomain-operator application, as
      79             :  * opposed to the global linear operator that typically requires
      80             :  * nearest-neighbor communications. See
      81             :  * `LinearSolver::Schwarz::SubdomainOperator` for details on how to
      82             :  * implement a subdomain operator for your problem.
      83             :  *
      84             :  * \par Algorithm overview:
      85             :  * In each iteration, the Schwarz algorithm computes the residual \f$r=b-Ax\f$,
      86             :  * restricts it to all subdomains as \f$r_s=R_s r\f$ and communicates data on
      87             :  * overlap regions with neighboring elements. Once an element has received all
      88             :  * data on its subdomain it solves the sub-problem \f$A_{ss}\delta x_s=r_s\f$
      89             :  * for the correction \f$\delta x_s\f$, where \f$\delta x_s\f$ and \f$r_s\f$
      90             :  * have data only on the points of the subdomain and \f$A_{ss}\f$ is the
      91             :  * subdomain operator. Since all elements perform such a subdomain-solve, we end
      92             :  * up with a subdomain solution \f$\delta x_s\f$ on every subdomain, and the
      93             :  * solutions overlap. Therefore, the algorithm communicates subdomain solutions
      94             :  * on overlap regions with neighboring elements and adds them to the solution
      95             :  * field \f$x\f$ as a weighted sum.
      96             :  *
      97             :  * \par Applying the global linear operator:
      98             :  * In order to compute the residual \f$r=b-Ax\f$ that will be restricted to the
      99             :  * subdomains to serve as source for the subdomain solves, we must apply the
     100             :  * global linear operator \f$A\f$ to the solution field \f$x\f$ once per Schwarz
     101             :  * iteration. The global linear operator typically introduces nearest-neighbor
     102             :  * couplings between elements but no global synchronization point (as opposed to
     103             :  * Krylov-type solvers such as `LinearSolver::gmres::Gmres` that require global
     104             :  * inner products in each iteration). The algorithm assumes that the action list
     105             :  * passed to `solve` applies the global linear operator, as described above.
     106             :  *
     107             :  * \par Subdomain solves:
     108             :  * Each subdomain solve is local to an element, since the data on the subdomain
     109             :  * has been made available through communication with the neighboring elements.
     110             :  * We can now use any means to solve the subdomain problem that's appropriate
     111             :  * for the subdomain operator \f$A_{ss}\f$. For example, if the subdomain
     112             :  * operator was available as an explicit matrix we could invert it directly.
     113             :  * Since it is typically a matrix-free operator a common approach to solve the
     114             :  * subdomain problem is by means of an iterative Krylov-type method, such as
     115             :  * GMRES or Conjugate Gradient, ideally with an appropriate preconditioner (yes,
     116             :  * this would be preconditioned Krylov-methods solving the subdomains of the
     117             :  * Schwarz solver, which might in turn precondition a global Krylov-solver -
     118             :  * it's preconditioners all the way down). Currently, the linear solvers listed
     119             :  * in `LinearSolver::Serial` are available to solve subdomains. They include
     120             :  * Krylov-type methods that support nested preconditioning with any other linear
     121             :  * solver. Additional problem-specific subdomain preconditioners can be added
     122             :  * with the `SubdomainPreconditioners` template parameter. Note that the choice
     123             :  * of subdomain solver (and, by extension, the choice of subdomain
     124             :  * preconditioner) affects only the _performance_ of the Schwarz solver, not its
     125             :  * convergence or parallelization properties (assuming the subdomain solutions
     126             :  * it produces are sufficiently precise).
     127             :  *
     128             :  * \par Weighting:
     129             :  * Once the subdomain solutions \f$\delta x_s\f$ have been found they must be
     130             :  * combined where they have multiple values, i.e. on overlap regions of the
     131             :  * subdomains. We use an additive approach where we combine the subdomain
     132             :  * solutions as a weighted sum, which has the advantage over "multiplicative"
     133             :  * Schwarz methods that the subdomains decouple and can be solved in parallel.
     134             :  * See `LinearSolver::Schwarz::element_weight` and Section 3.1 of
     135             :  * \cite Stiller2016a for details. Note that we must account for the missing
     136             :  * corner- and edge-neighbors when constructing the weights. See
     137             :  * `LinearSolver::Schwarz::intruding_weight` for a discussion.
     138             :  *
     139             :  * \par Array sections
     140             :  * This linear solver requires no synchronization between elements, so it runs
     141             :  * on all elements in the array parallel component. Partitioning of the elements
     142             :  * in sections is only relevant for observing residual norms. Pass the section
     143             :  * ID tag for the `ArraySectionIdTag` template parameter if residual norms
     144             :  * should be computed over a section. Pass `void` (default) to compute residual
     145             :  * norms over all elements in the array.
     146             :  *
     147             :  * \par Possible improvements:
     148             :  * - Specify the number of overlap points as a fraction of the element width
     149             :  * instead of a fixed number. This was shown in \cite Stiller2016b to achieve
     150             :  * scale-independent convergence at the cost of increasing subdomain sizes.
     151             :  */
     152             : template <typename FieldsTag, typename OptionsGroup, typename SubdomainOperator,
     153             :           typename SubdomainPreconditioners = tmpl::list<>,
     154             :           typename SourceTag =
     155             :               db::add_tag_prefix<::Tags::FixedSource, FieldsTag>,
     156             :           typename ArraySectionIdTag = void>
     157           1 : struct Schwarz {
     158           0 :   using operand_tag = FieldsTag;
     159           0 :   using fields_tag = FieldsTag;
     160           0 :   using source_tag = SourceTag;
     161           0 :   using options_group = OptionsGroup;
     162             : 
     163           0 :   using component_list = tmpl::list<>;
     164           0 :   using observed_reduction_data_tags = observers::make_reduction_data_tags<
     165             :       tmpl::list<async_solvers::reduction_data, detail::reduction_data>>;
     166           0 :   using subdomain_solver =
     167             :       detail::subdomain_solver<FieldsTag, SubdomainOperator,
     168             :                                SubdomainPreconditioners>;
     169             : 
     170           0 :   using initialize_element = tmpl::list<
     171             :       async_solvers::InitializeElement<FieldsTag, OptionsGroup, SourceTag>,
     172             :       detail::InitializeElement<FieldsTag, OptionsGroup, SubdomainOperator,
     173             :                                 SubdomainPreconditioners>>;
     174             : 
     175           0 :   using amr_projectors = initialize_element;
     176             : 
     177           0 :   using register_element =
     178             :       tmpl::list<async_solvers::RegisterElement<FieldsTag, OptionsGroup,
     179             :                                                 SourceTag, ArraySectionIdTag>,
     180             :                  detail::RegisterElement<FieldsTag, OptionsGroup, SourceTag,
     181             :                                          ArraySectionIdTag>>;
     182             : 
     183             :   template <typename ApplyOperatorActions, typename Label = OptionsGroup>
     184           0 :   using solve = tmpl::list<
     185             :       async_solvers::PrepareSolve<FieldsTag, OptionsGroup, SourceTag, Label,
     186             :                                   ArraySectionIdTag>,
     187             :       detail::SendOverlapData<FieldsTag, OptionsGroup, SubdomainOperator>,
     188             :       detail::SolveSubdomain<FieldsTag, OptionsGroup, SubdomainOperator,
     189             :                              ArraySectionIdTag>,
     190             :       detail::ReceiveOverlapSolution<FieldsTag, OptionsGroup,
     191             :                                      SubdomainOperator>,
     192             :       ApplyOperatorActions,
     193             :       async_solvers::CompleteStep<FieldsTag, OptionsGroup, SourceTag, Label,
     194             :                                   ArraySectionIdTag>>;
     195             : };
     196             : 
     197             : }  // namespace LinearSolver::Schwarz