SpECTRE Documentation Coverage Report
 Current view: top level - ParallelAlgorithms/LinearSolver/Schwarz - Schwarz.hpp Hit Total Coverage Commit: edb1b5199a4a86c269aedbb831767801169f3e8a Lines: 1 12 8.3 % Date: 2021-04-19 16:23:01 Legend: Lines: hit not hit
  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. 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, we use 119 : * LinearSolver::Serial::Gmres to solve subdomains. It supports 120 : * preconditioning, and adding useful subdomain preconditioners will be the 121 : * subject of future work. Note that the choice of subdomain solver (and, by 122 : * extension, the choice of subdomain preconditioner) affects only the 123 : * _performance_ of the Schwarz solver, not its convergence or parallelization 124 : * properties (assuming the subdomain solutions it produces are sufficiently 125 : * precise). 126 : * 127 : * \par Weighting: 128 : * Once the subdomain solutions \f$\delta x_s\f$ have been found they must be 129 : * combined where they have multiple values, i.e. on overlap regions of the 130 : * subdomains. We use an additive approach where we combine the subdomain 131 : * solutions as a weighted sum, which has the advantage over "multiplicative" 132 : * Schwarz methods that the subdomains decouple and can be solved in parallel. 133 : * See LinearSolver::Schwarz::element_weight and Section 3.1 of 134 : * \cite Stiller2016a for details. Note that we must account for the missing 135 : * corner- and edge-neighbors when constructing the weights. See 136 : * LinearSolver::Schwarz::intruding_weight for a discussion. 137 : * 138 : * \par Possible improvements: 139 : * - Specify the number of overlap points as a fraction of the element width 140 : * instead of a fixed number. This was shown in \cite Stiller2016b to achieve 141 : * scale-independent convergence at the cost of increasing subdomain sizes. 142 : * - Subdomain preconditioning (see paragraph on subdomain solves above) 143 : */ 144 : template > 147 1 : struct Schwarz { 148 0 : using operand_tag = FieldsTag; 149 0 : using fields_tag = FieldsTag; 150 0 : using source_tag = SourceTag; 151 0 : using options_group = OptionsGroup; 152 : 153 0 : using component_list = tmpl::list<>; 154 0 : using observed_reduction_data_tags = observers::make_reduction_data_tags< 155 : tmpl::list>; 156 0 : using subdomain_solver = 157 : detail::subdomain_solver; 158 : 159 0 : using initialize_element = tmpl::list< 160 : async_solvers::InitializeElement, 161 : detail::InitializeElement>; 162 : 163 0 : using register_element = tmpl::list< 164 : async_solvers::RegisterElement, 165 : detail::RegisterElement>; 166 : 167 : template 168 0 : using solve = tmpl::list< 169 : async_solvers::PrepareSolve, 170 : detail::SendOverlapData, 171 : detail::ReceiveOverlapData, 172 : detail::SolveSubdomain, 173 : detail::ReceiveOverlapSolution, 175 : ApplyOperatorActions, 176 : async_solvers::CompleteStep>; 177 : }; 178 : 179 : } // namespace LinearSolver::Schwarz 

 Generated by: LCOV version 1.14