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