An additive Schwarz subdomain solver for linear systems of equations \(Ax=b\). More...

#include <Schwarz.hpp>

Public Types
using	operand_tag = FieldsTag

using	fields_tag = FieldsTag

using	source_tag = SourceTag

using	options_group = OptionsGroup

using	component_list = tmpl::list<>

using	observed_reduction_data_tags = observers::make_reduction_data_tags< tmpl::list< async_solvers::reduction_data, detail::reduction_data > >

using	subdomain_solver = detail::subdomain_solver< FieldsTag, SubdomainOperator, SubdomainPreconditioners >

using	initialize_element = tmpl::list< async_solvers::InitializeElement< FieldsTag, OptionsGroup, SourceTag >, detail::InitializeElement< FieldsTag, OptionsGroup, SubdomainOperator, SubdomainPreconditioners > >

using	amr_projectors = initialize_element

using	register_element = tmpl::list< async_solvers::RegisterElement< FieldsTag, OptionsGroup, SourceTag, ArraySectionIdTag >, detail::RegisterElement< FieldsTag, OptionsGroup, SourceTag, ArraySectionIdTag > >

template<typename ApplyOperatorActions , typename Label = OptionsGroup>
using	solve = tmpl::list< async_solvers::PrepareSolve< FieldsTag, OptionsGroup, SourceTag, Label, ArraySectionIdTag >, detail::SendOverlapData< FieldsTag, OptionsGroup, SubdomainOperator >, detail::SolveSubdomain< FieldsTag, OptionsGroup, SubdomainOperator, ArraySectionIdTag >, detail::ReceiveOverlapSolution< FieldsTag, OptionsGroup, SubdomainOperator >, ApplyOperatorActions, async_solvers::CompleteStep< FieldsTag, OptionsGroup, SourceTag, Label, ArraySectionIdTag > >

Detailed Description

template<typename FieldsTag, typename OptionsGroup, typename SubdomainOperator, typename SubdomainPreconditioners = tmpl::list<>, typename SourceTag = db::add_tag_prefix<::Tags::FixedSource, FieldsTag>, typename ArraySectionIdTag = void>
struct LinearSolver::Schwarz::Schwarz< FieldsTag, OptionsGroup, SubdomainOperator, SubdomainPreconditioners, SourceTag, ArraySectionIdTag >

An additive Schwarz subdomain solver for linear systems of equations \(Ax=b\).

This Schwarz-type linear solver works by solving many sub-problems in parallel and combining their solutions as a weighted sum to converge towards the global solution. Each sub-problem is the restriction of the global problem to an element-centered subdomain, which consists of a central element and an overlap region with its neighbors. The decomposition into independent sub-problems makes this linear solver very parallelizable, avoiding global synchronization points altogether. It is commonly used as a preconditioner to Krylov-type solvers such as LinearSolver::gmres::Gmres or LinearSolver::cg::ConjugateGradient, or as the "smoother" for a Multigrid solver (which may in turn precondition a Krylov-type solver).

This linear solver relies on an implementation of the global linear operator \(A(x)\) and its restriction to a subdomain \(A_{ss}(x)\). Each step of the algorithm expects that \(A(x)\) is computed and stored in the DataBox as db::add_tag_prefix<LinearSolver::Tags::OperatorAppliedTo, operand_tag>. To perform a solve, add the solve action list to an array parallel component. Pass the actions that compute \(A(x)\), as well as any further actions you wish to run in each step of the algorithm, as the first template parameter to solve. If you add the solve action list multiple times, use the second template parameter to label each solve with a different type. Pass the subdomain operator as the SubdomainOperator template parameter to this class (not to the solve action list template). See the paragraph below for information on implementing a subdomain operator.

Subdomain geometry:: The image below gives an overview of the structure of an element-centered subdomain:

Fig. 1: An element-centered subdomain

Fig. 1 shows part of a 2-dimensional computational domain. The domain is composed of elements (light gray) that each carry a Legendre-Gauss-Lobatto mesh of grid points (black dots). The diagonally shaded region to the left illustrates an external domain boundary. The lines between neighboring element faces illustrate mortar meshes, which are relevant when implementing a subdomain operator in a DG context but play no role in the Schwarz algorithm. The dashed line gives an example of an element-centered subdomain with a maximum of 2 grid points of overlap into neighboring elements. For details on how the number of overlap points is chosen see LinearSolver::Schwarz::overlap_extent. Note that the subdomain does not extend into corner- or edge-neighbors, which is different to both [173] and [186]. The reason we don't include such diagonal couplings is that in a DG context information only propagates across faces, as noted in [173]. By eliminating the corner- and edge-neighbors we significantly reduce the complexity of the subdomain geometry and potentially also the communication costs. However, the advantages and disadvantages of this choice have yet to be carefully evaluated.

Subdomain operator:: The Schwarz subdomain solver relies on a restriction of the global linear problem \(Ax=b\) to the subdomains. The subdomain operator \(A_{ss}=R_s A R_s^T\), where \(R_s\) is the restriction operator, essentially assumes that its operand is zero everywhere but within the subdomain (see Section 3.1 in [173]). Therefore it can be evaluated entirely on an element-centered subdomain with no need to communicate with neighbors within each subdomain-operator application, as opposed to the global linear operator that typically requires nearest-neighbor communications. See LinearSolver::Schwarz::SubdomainOperator for details on how to implement a subdomain operator for your problem.

Algorithm overview:: In each iteration, the Schwarz algorithm computes the residual \(r=b-Ax\), restricts it to all subdomains as \(r_s=R_s r\) and communicates data on overlap regions with neighboring elements. Once an element has received all data on its subdomain it solves the sub-problem \(A_{ss}\delta x_s=r_s\) for the correction \(\delta x_s\), where \(\delta x_s\) and \(r_s\) have data only on the points of the subdomain and \(A_{ss}\) is the subdomain operator. Since all elements perform such a subdomain-solve, we end up with a subdomain solution \(\delta x_s\) on every subdomain, and the solutions overlap. Therefore, the algorithm communicates subdomain solutions on overlap regions with neighboring elements and adds them to the solution field \(x\) as a weighted sum.

Applying the global linear operator:: In order to compute the residual \(r=b-Ax\) that will be restricted to the subdomains to serve as source for the subdomain solves, we must apply the global linear operator \(A\) to the solution field \(x\) once per Schwarz iteration. The global linear operator typically introduces nearest-neighbor couplings between elements but no global synchronization point (as opposed to Krylov-type solvers such as LinearSolver::gmres::Gmres that require global inner products in each iteration). The algorithm assumes that the action list passed to solve applies the global linear operator, as described above.

Subdomain solves:: Each subdomain solve is local to an element, since the data on the subdomain has been made available through communication with the neighboring elements. We can now use any means to solve the subdomain problem that's appropriate for the subdomain operator \(A_{ss}\). For example, if the subdomain operator was available as an explicit matrix we could invert it directly. Since it is typically a matrix-free operator a common approach to solve the subdomain problem is by means of an iterative Krylov-type method, such as GMRES or Conjugate Gradient, ideally with an appropriate preconditioner (yes, this would be preconditioned Krylov-methods solving the subdomains of the Schwarz solver, which might in turn precondition a global Krylov-solver - it's preconditioners all the way down). Currently, the linear solvers listed in LinearSolver::Serial are available to solve subdomains. They include Krylov-type methods that support nested preconditioning with any other linear solver. Additional problem-specific subdomain preconditioners can be added with the SubdomainPreconditioners template parameter. Note that the choice of subdomain solver (and, by extension, the choice of subdomain preconditioner) affects only the performance of the Schwarz solver, not its convergence or parallelization properties (assuming the subdomain solutions it produces are sufficiently precise).

Weighting:: Once the subdomain solutions \(\delta x_s\) have been found they must be combined where they have multiple values, i.e. on overlap regions of the subdomains. We use an additive approach where we combine the subdomain solutions as a weighted sum, which has the advantage over "multiplicative" Schwarz methods that the subdomains decouple and can be solved in parallel. See LinearSolver::Schwarz::element_weight and Section 3.1 of [173] for details. Note that we must account for the missing corner- and edge-neighbors when constructing the weights. See LinearSolver::Schwarz::intruding_weight for a discussion.

Array sections: This linear solver requires no synchronization between elements, so it runs on all elements in the array parallel component. Partitioning of the elements in sections is only relevant for observing residual norms. Pass the section ID tag for the ArraySectionIdTag template parameter if residual norms should be computed over a section. Pass void (default) to compute residual norms over all elements in the array.

Possible improvements:

Specify the number of overlap points as a fraction of the element width instead of a fixed number. This was shown in [172] to achieve scale-independent convergence at the cost of increasing subdomain sizes.

The documentation for this struct was generated from the following file:

src/ParallelAlgorithms/LinearSolver/Schwarz/Schwarz.hpp

Public Types

Detailed Description