Preconditioned nonlinear iterations for overlapping Chebyshev discretizations with independent grids

TL;DR

This paper introduces a nonlinear preconditioning approach using overlapping Chebyshev discretizations with independent grids, reducing communication and enabling parallel solutions for PDEs.

Contribution

It presents a novel nonlinear Schwarz preconditioning method that operates on independent subdomain discretizations, avoiding restrictive updates and minimizing communication.

Findings

Reduces interprocess communication in nonlinear PDE solvers

Enables parallel independent subdomain problem solving

Offers a flexible preconditioning framework for nonlinear discretizations

Abstract

The additive Schwarz method is usually presented as a preconditioner for a PDE linearization based on overlapping subsets of nodes from a global discretization. It has previously been shown how to apply Schwarz preconditioning to a nonlinear problem. By first replacing the original global PDE with the Schwarz overlapping problem, the global discretization becomes a simple union of subdomain discretizations, and unknowns do not need to be shared. In this way restrictive-type updates can be avoided, and subdomains need to communicate only via interface interpolations. The resulting preconditioner can be applied linearly or nonlinearly. In the latter case nonlinear subdomain problems are solved independently in parallel, and the frequency and amount of interprocess communication can be greatly reduced compared to linearized preconditioning.

Tables

Table 1. Table 1: Parallel timing results for the Burgers equation experiment. “Jacobian time” is the total time spent within applications of the Jacobian to a given vector, and “Residual time” is the total amount of time spent evaluating the nonlinear SNK residual.

Number	Total time	Speedup	Jacobian	Residual	Residual
of cores	(sec.)		time	time	speedup
1	62.3	—	11.6	42.0	—
2	43.5	1.43	12.7	24.1	1.74
4	33.4	1.86	11.3	15.5	2.71
6	30.0	2.08	11.5	12.1	3.47
8	28.3	2.21	11.3	10.5	4.00
12	27.6	2.26	11.9	9.3	4.51
16	28.6	2.18	13.2	8.8	4.78
20	28.2	2.21	13.1	8.3	5.04

Equations

$$\varphi (\mathbf{x},u)$$

$$\beta (\mathbf{x},u)$$

$$\varphi (\mathbf{x},u_i)$$

$$\beta (\mathbf{x},u_i)$$

$$u_i$$

$${\mathbf{f}}_i({\mathbf{u}}_i)=\{\begin{array}{ll}\varphi (\mathbf{x},{\tilde{u}}_i)& \text{ for all }\mathbf{x}\in X_i,\\ \beta (\mathbf{x},{\tilde{u}}_i)& \text{ for all }\mathbf{x}\in G_{i0},\\ {\tilde{u}}_i(\mathbf{x})& \text{ for all }\mathbf{x}\in G_{ij},j=1,\dots ,i-1,i+1,\dots ,N.\end{array}$$

$${\mathbf{T}}_i\mathbf{u}=\{\begin{array}{ll}0& \text{ for all }\mathbf{x}\in X_i,\\ 0& \text{ for all }\mathbf{x}\in G_{i0},\\ {\tilde{u}}_j(\mathbf{x})& \text{ for all }\mathbf{x}\in G_{ij},j=1,\dots ,i-1,i+1,\dots ,N.\end{array}$$

$$\mathbf{f}(\mathbf{u})-\mathbf{T}\mathbf{u}=0.$$

$$\hat{\mathbf{f}}(\hat{\mathbf{e}}+\mathbf{R}\mathbf{u})-\hat{\mathbf{T}}\hat{\mathbf{e}}-\hat{\mathbf{f}}(\mathbf{R}\mathbf{u})+\mathbf{R}(\mathbf{f}(\mathbf{u})-\mathbf{T}\mathbf{u})=0$$

$$\bigl{[}\hat{\bm{f}}^{\prime}(\hat{\bm{e}}+\bm{R}\bm{u})-\hat{\bm{T}}\bigr{]}\frac{\partial\hat{\bm{e}}}{\partial\bm{u}}=-\bigl{(}\hat{\bm{f}}^{\prime}(\hat{\bm{e}}+\bm{R}\bm{u})-\hat{\bm{f}}^{\prime}(\bm{R}\bm{u})\bigr{)}\bm{R}-\bm{R}\bigl{(}\bm{f}^{\prime}(\bm{u})-\bm{T}\bigr{)}.$$

$${\mathbf{u}}^{†}$$

$${\mathbf{u}}^{\text{new}}$$

$$\mathbf{h}(\mathbf{u}):=\mathbf{c}(\mathbf{u})+(\mathbf{f}-\mathbf{T})(\mathbf{u}+\mathbf{c}(\mathbf{u})).$$

$${\mathbf{h}}^{\prime }(\mathbf{u})={\mathbf{c}}^{\prime }(\mathbf{u})+({\mathbf{f}}^{\prime }-\mathbf{T})(\mathbf{u}+\mathbf{c}(\mathbf{u}))\cdot (\mathbf{I}+{\mathbf{c}}^{\prime }(\mathbf{u})).$$

$$\mathbf{g}(\mathbf{u}):=\mathbf{u}-{\mathbf{f}}^{-1}(\mathbf{T}\mathbf{u}).$$

$${\mathbf{f}}_i^{\prime }({\mathbf{u}}_i-{\mathbf{z}}_i)\left[I-\frac{\partial {\mathbf{z}}_i}{\partial \mathbf{u}}\right]={\mathbf{T}}_i,$$

$$\Vert F(x_k)+F^{\prime }(x_k)s_k\Vert \le {\eta }_k\Vert F(x_k)\Vert$$

$${\eta }_k=10^{-4}{\left(\frac{\Vert F(x_k)\Vert }{\Vert F(x_{k-1})\Vert }\right)}^2.$$

$$-\Delta u-\frac{\partial \omega }{\partial y}$$

$$-\Delta v+\frac{\partial \omega }{\partial x}$$

$$-\frac{1}{\text{Re}}\Delta \omega +u\frac{\partial \omega }{\partial x}+v\frac{\partial \omega }{\partial y}$$

$$u$$

$$v$$

$$\omega$$

$$\nu \Delta u-u\cdot \nabla =0,$$

$$u=\arctan \left(\cos \left(\frac{3\pi }{16}\right)x+\sin \left(\frac{3\pi }{16}\right)y\right).$$

$$\exp \left(\frac{1}{1-x^{-20}}+\frac{1}{1-y^{-20}}\right)$$

Full text

{AMS}

65N55, 33F05, 97N40

1 Introduction

Overlapping domain decomposition has been recognized as a valuable aid in solving partial differential equations since Schwarz first described his alternating method in 1870. (For straightforward introductions to the topic, see [7, 12]; for a more historical perspective, see [9].) Overlapping decomposition provides a way to solve a problem on a global domain by exploiting its reduction to smaller subdomains. This creates geometric flexibility and allows special effort to be focused on small parts of the domain when appropriate. Domain decomposition also has a natural parallelism that is particularly attractive in the increasingly multicore context of scientific computing.

For a linear PDE, one typically seeks to apply a preconditioner for a Krylov iteration such as GMRES, in the form of solving problems on overlapping subdomains whose boundary data is in part determined by values of the solution in other subdomains. In the parallel context this is achieved by an additive Schwarz (AS) scheme. When one partitions the unknowns of a global discretization into overlapping subsets, the best form of AS are restricted AS (RAS) methods [4], which do not allow multiple domains to update shared unknowns independently and thus over-correct. Typically, then, the subdomain problems are solved on overlapping sets, but the results are distributed in a nonoverlapping fashion.

For nonlinear problems, the obvious extension of AS preconditioning is to apply it as described above on the linearized equations that are determined by a quasi-Newton iteration. We refer to this process as a Newton–Krylov–Schwarz (NKS) procedure, reflecting the nesting order of the different elements of linearization, linear solver, and preconditioning.

Cai and Keyes [3] proposed instead modifying the nonlinear problem using the Schwarz ansatz. In addition to yielding a preconditioned linearization for the Krylov solver, the preconditioned nonlinear problem exhibited more robust convergence for the Newton iteration than did the original nonlinear problem. They called their method ASPIN, short for additive Schwarz preconditioned inexact Newton. As a technical matter, they did not recommend applying the true Jacobian of the system, preferring an approximation that required less effort. Subsequently, Dolean et al. [6] pointed out that Cai and Keyes did not use the RAS form of AS preconditioning, and they proposed an improved variant called RASPEN that does. We refer to this type of nonlinear preconditioning as Schwarz–Newton–Krylov (SNK), because the Schwarz ansatz is applied before the linearization begins.

Our interest is in applying the nonlinear preconditioning technique to spectral collocation discretizations in overlapping rectangles or cuboids, leading to globally smooth approximations constructed from a partition of unity [1]. In this context, there is not naturally a single global discretization whose degrees of freedom are partitioned into overlapping sets, because the Chebyshev (or Legendre, or other classical) nodes will not generally coincide within the overlapping regions. In principle one could link the degrees of freedom within overlap regions by interpolating between subdomains, but this process adds complication, computational time, and (in the parallel context) communication of data.

Here we present an alternative strategy that begins by replacing the original PDE problem with the Schwarz problems on the union of the subdomains. That is, rather than regarding the subdomains as solving the global PDE on a region that includes portions shared with other subdomains, each subdomain has a “private copy” of its entire region and is free to have its own solution values throughout. Of course, the new global problem is not solved until the interface values of every subdomain agree with values interpolated from other subdomains that contain the interface. As a Schwarz starting point, our technique has both NKS and SNK variants.

One advantage of this new formulation is that interpolations need to be done only on lower-dimensional interfaces, rather than throughout the overlap regions. Another is that plain AS is preferred to RAS, because each subdomain has to update its own values separately. We show that it is straightforward to implement exact Jacobians for SNK with nothing more than the ability to do fully local PDE nonlinear and linearized solves, plus the ability to transfer values between subdomains through interface interpolations. We also derive a two-level method to prevent convergence degradation as the number of subdomains increases. The performance of the NKS and SNK methods is validated and compared through several numerical experiments.

2 PDE problem and multidomain formulation

The main goal of this work is to solve the PDE

[TABLE]

where $u(\bm{x})$ is the unknown solution and $\phi$ and $\beta$ are nonlinear differential operators (with $\phi$ being of higher order). (We can easily extend to the case where $u$, $\phi$, and $\beta$ are vector-valued, but we use scalars to calm the notation.) Many Schwarz-based algorithms for (1) begin with a global discretization whose solution is accelerated by an overlapping domain decomposition. In this situation, some of the numerical degrees of freedom are shared by multiple subdomains—either directly or through interpolation—and proper use of additive Schwarz (AS) calls for the restricted-AS (RAS) implementation, which essentially insures that updates of shared values are done only once from the global perspective, not independently by the subdomains.

We take a different approach, replacing the original problem (1) with

[TABLE]

where now $u_{1},\ldots,u_{N}$ are unknown functions on overlapping subdomains $\Omega_{i}$ that cover $\Omega$, and the $Z_{i}$ are nonoverlapping zones lying within the respective subdomains. Clearly any strong solution of (1) is also a solution of (2), and while the converse is not necessarily so in principle, we regard the possibility of finding a solution of (2) that is not also a solution of (1) as remote in practice.

The key consequence of starting from (2) as the global problem is that each overlapping region is covered separately by the involved subdomains; each is free to update its representation independently in order to converge to a solution. From one point of view, our discretizations of the overlap regions are redundant and somewhat wasteful. However, the fraction of redundant discrete unknowns is very modest. In return, we only need to interpolate on the interfaces, there is no need to use the RAS formulation, and the coarsening needed for a two-level variant is trivial (see section 2.3).

2.1 Discretization

We now describe a collocation discretization of (2) for concreteness. Each subfunction $u_{i}(\bm{x})$ is discretized by a vector $\bm{u}_{i}$ of length $n_{i}$. By \bm{u}=\left\llbracket\bm{u}_{i}\right\rrbracket we mean a concatenation of all the discrete unknowns over subdomains $i=1,\ldots,N$ into a single vector. Subdomain $\Omega_{i}$ is discretized by a node set $X_{i}\subset\overline{\Omega}_{i}$ and a boundary node set $B_{i}\subset\partial\Omega_{i}$. The total cardinality of $X_{i}$ and $B_{i}$ together is also $n_{i}$. The boundary nodes are subdivided into nonintersecting sets $G_{ij}=B_{i}\cap Z_{j}$ for $j\neq i$, and $G_{i0}=B_{i}\cap\partial\Omega$.

For each $i$, the vector $\bm{u}_{i}$ defines a function $\tilde{u}_{i}(\bm{x})$ on $\Omega_{i}$. These can be used to evaluate $\phi$ and $\beta$ from (2) anywhere in $\Omega_{i}$. We define an $n_{i}$-dimensional vector function $\bm{f}_{i}$ as the concatenation of three vectors:

[TABLE]

In addition, we have the linear transfer operator $\bm{T}_{i}$ defined by

[TABLE]

Note that while $\bm{f}_{i}$ is purely local to subdomain $i$, the transfer operator $\bm{T}_{i}$ operates on the complete discretization $\bm{u}$, as it interpolates from “foreign” subdomains onto the parts of $B_{i}$ lying inside $\Omega$. Finally, we are able to express the complete discretization of (2) through concatenations over the subdomains. Let \bm{u}=\left\llbracket\bm{u}_{i}\right\rrbracket, \bm{f}(\bm{u})=\left\llbracket\bm{f}_{i}(\bm{u}_{i})\right\rrbracket, and \bm{T}\bm{u}=\left\llbracket\bm{T}_{i}\bm{u}\right\rrbracket. Then the discrete form of (2) is the nonlinear equation

[TABLE]

For a square discretization, the goal is to solve (5), while in the least-squares case, the goal is to minimize $\bm{f}(\bm{u})-\bm{T}\bm{u}$ in the (possibly weighted) 2-norm.

2.2 Newton–Krylov–Schwarz

The standard approach to (5) for a large discretization is to apply an inexact Newton iteration with a Krylov subspace solver for finding correcting steps from the linearization. Within the Krylov solver we have a natural setting for applying an AS preconditioner. Specifically, if we have a proposed approximate solution $\bm{u}$, we evaluate the nonlinear residual $\bm{r}=\bm{f}(\bm{u})-\bm{T}\bm{u}$. We then (inexactly, perhaps) solve the linearization \bigl{[}\bm{f}^{\prime}(\bm{u})-\bm{T}\bm{u}\bigr{]}\bm{s}=-\bm{r} for the Newton correction $\bm{s}$, using a Krylov solver such as GMRES. These iterations are preconditioned by the block diagonal matrix $\bm{f}^{\prime}(\bm{u})$, which is simply the block diagonal of the subdomain Jacobians $\bm{f}_{i}^{\prime}(\bm{u}_{i})$. We refer to this method as Newton–Krylov–Schwarz, or NKS.

Implementation of NKS requires three major elements: the evaluations of $\bm{f}(\bm{u})$ and $\bm{T}\bm{u}$ for given $\bm{u}$, the application of the Jacobian $\bm{f}^{\prime}(\bm{u})$ to a given vector $\bm{v}$, and the inversion of $\bm{f}^{\prime}(\bm{u})$ for given data. All of the processes involving $\bm{f}$ are embarrassingly parallel and correspond to standard steps in solving the PDE on the local subdomains. Each application of the transfer operator $\bm{T}$, however, requires a communication from each subdomain to its overlapping neighbors, as outlined in Algorithm 1. This step occurs once in evaluating the nonlinear residual and in every GMRES iteration to apply the Jacobian. In a parallel code, the communication steps could be expected to be a major factor in the performance of the method.

2.3 Two-level scheme

As is well known [7], AS schemes should incorporate a coarse solution step in order to maintain convergence rates as the number of subdomains increases. The methods described above depend on the subdomain discretization sizes $n_{i}$ of the collocation nodes and solution representation, respectively. Now suppose we decrease the discretization sizes to $\hat{n}_{i}$, and denote the corresponding discretizations of (5) by $\hat{\bm{f}}(\hat{\bm{u}})-\hat{\bm{T}}\hat{\bm{u}}=\bm{0}$. We can define a restriction operator $\bm{R}$ that maps fine-scale vectors to their coarse counterparts. This operator is block diagonal, i.e., it can be applied independently within the subdomains. We can also construct a block diagonal prolongation operator $\bm{P}$ for mapping the solution representation from coarse to fine scales.

We are then able to apply the standard Full Approximation Scheme (FAS) using the coarsened problem [2]. Specifically, we solve the coarse problem

[TABLE]

for the coarse correction $\hat{\bm{e}}$, and define $\bm{c}(\bm{u})=\bm{P}\hat{\bm{e}}$ as the FAS corrector at the fine level. The procedure for calculating $\bm{c}$ is outlined in Algorithm 2.

We also require the action of the Jacobian $\frac{\partial\bm{c}}{\partial\bm{u}}=\bm{P}\frac{\partial\hat{\bm{e}}}{\partial\bm{u}}$ on a given vector $\bm{v}$. It is straightforward to derive from (6) that

[TABLE]

Note that the matrix $\hat{\bm{f}}^{\prime}(\hat{\bm{e}}+\bm{R}\bm{u})$ should be available at no extra cost from the end of the Newton solution of (6). Algorithm 3 describes the corresponding algorithm for computing the application of $\bm{c}^{\prime}(\bm{u})$ to any vector $\bm{v}$. Even though $\bm{c}^{\prime}$ is of the size of the fine discretization, the computation requires only coarse-dimension dense linear algebra.

Finally, we describe how to combine coarsening with the preconditioned fine scale into a two-level algorithm. If we were to alternate coarse and fine corrections in the classical fixed-point form,

[TABLE]

then we are effectively seeking a root of

[TABLE]

Finally, the Jacobian of the combined map is straightforwardly

[TABLE]

Thus the action of $\bm{h^{\prime}}$ on a vector can be calculated using the algorithms for $\bm{c}^{\prime}$, $\bm{f}^{\prime}$, and $\bm{T}$.

3 Preconditioned nonlinear iterations

As shown in section 2.2, the inner Krylov iterations of the NKS method are governed by the preconditioned Jacobian $\bm{I}-[\bm{f}^{\prime}(\bm{v})]^{-1}\bm{T}$. Following the observation of Cai and Keyes [3], we next derive a method that applies Krylov iterations to the same matrix, but arising as the natural result of preconditioning the nonlinear problem. Specifically, we precondition (5) by finding a root of the nonlinear operator

[TABLE]

Evaluation of $\bm{g}$ is feasible because of the block diagonal (that is, fully subdomain-local) action of the nonlinear $\bm{f}$. Since we are therefore applying the Schwarz ansatz even before linearizing the problem, we refer to the resulting method as Schwarz–Newton–Krylov (SNK).

We have several motivations for a method based on (10). First, one hopes that the nonlinear problem, being a (low-rank) perturbation of the identity operator, is somehow easier to solve by Newton stepping than the original form is. Second, the inversion of $\bm{f}$ means solving independent nonlinear problems on the subdomains with no communication, which well exploits parallelism. Finally, the same structure means that problems with relatively small highly active regions could be isolate the need to solve a nonlinear problem to that region, rather than having it be part of a fully coupled global nonlinear problem.

An algorithm for evaluating $\bm{g}$ is given in Algorithm 4. It requires one communication between subdomains to transfer interface data, followed by solving (in parallel if desired) the nonlinear subdomain problems $\bm{f}_{i}$ defined in (3). Note that the local problem in $\Omega_{i}$ is a discretization of the PDE with zero boundary data on the true boundary $\Gamma_{i0}$ and values transferred from the foreign subdomains on the interfaces.

Using the notation of Algorithm 4, we have that

[TABLE]

which implies that applying $[\partial\bm{z}_{i}/\partial\bm{u}]$ to a vector requires a single linear solve on a subdomain, with a matrix that is presumably already available at the end of the local Newton iteration used to compute $\bm{g}$. The process for applying $\bm{g}^{\prime}(\bm{u})$ to a vector is outlined in Algorithm 5.

3.1 Two-level scheme

The SNK method can be expected to require coarse correction steps to cope with a growing number of subdomains. An obvious approach to incorporating a coarse-grid correction step is to apply FAS directly, i.e., using the analog of (6) with fine $\bm{g}$ and coarse $\hat{\bm{g}}$. However, doing so means inverting $\hat{\bm{g}}$, which introduces another layer of iteration in the overall process.

We have found that it is simpler and successful to apply the FAS correction in the form of the original NKS method, as given in section 2.3. All we need to do is replace $\bm{f}(\bm{u})-\bm{T}\bm{u}$ and $\bm{f}^{\prime}-\bm{T}$ in (8) and (9) by $\bm{g}$ and $\bm{g}^{\prime}$, respectively.

4 Numerical experiments

For all experiments we compared three methods: NKS, SNK, and the two-level SNK2. The local nonlinear problems for SNK, and the coarse global problems in SNK2, were solved using fsolve from the Optimization Toolbox. Each solver used an inexact Newton method as the outer iteration, continued until the residual was less than $10^{-10}$ relative to the initial residual. For the inner iterations we used MATLAB’s gmres to solve for the Newton step $s_{k}$ such that

[TABLE]

where $\eta_{0}=10^{-4}$ and

[TABLE]

Given certain conditions on $F(x)$, if the intial solution is close enough to the true solution then this set of tolerances will yield a sequence with near q-2 convergence [8].

4.1 Regularized driven cavity flow

The first example is a regularized form of the lid-driven cavity flow problem [10], where we replace the boundary conditions with infinitely smooth ones. Using the velocity-vorticity formulation, in terms of the velocity $u,v$ and vorticity $\omega$ on $\Omega=[0,1]^{2}$ we have the nondimensionalized equations

[TABLE]

where Re is the Reynolds number. On the boundary $\partial\Omega$ we apply

[TABLE]

similar to the boundary conditions in [11].

We divided $\Omega$ into 16 overlapping patches of equal size (i.e. a 4 by 4 patch structure). Each subdomain was discretized by a second-kind Chebyshev grid of length 33 in each dimension. For the initial guess to the outer solver iterations we extended the boundary conditions (14) to all values of $x$ in the square $\Omega$.

The convergence of the three solvers is shown for $\text{Re}=100$, and SNK and NKS for $\text{Re}=1000$ in Figure 1. All three methods converge for $\text{Re}=100$. The number of GMRES iterations per nonlinear iteration are similar for the NKS and SNK methods; this is to be expected since the linear system used to solve the Newton step is similar in both methods. We do however see a dramatic reduction in the number of GMRES iterations with the two-level SNK method.

It is worth noting well that the computational time for each outer iteration varies greatly between solvers. The NKS residual requires only evaluating the discretized PDE and is thus is a good deal faster per iteration than the SNK solvers, which require solving the local nonlinear problems. In addition SNK2 must solve a global coarse problem in each outer iteration, but this added relatively little computing time.

For the higher Reynolds number $\text{Re}=1000$ we find that while SNK still converges, NKS does not, similar to what was reported in [3]. We also found that the coarse-level solver in SNK2 had trouble converging.

4.2 Burgers equation

The second test problem is Burgers’ equation,

[TABLE]

on $\Omega=[-1,1]^{2}$, with Dirichlet boundary condition

[TABLE]

This PDE was solved using a subdomain structure adapted to the function

[TABLE]

using the methods in [1], in order to help capture the boundary layers, as shown in Figure 2. For the initial guess of the outer iterations, the boundary condition (16) was extended throughout $\Omega$.

Convergence histories for (15)–(16) for $1/\nu=400,800,1000,1500$ are given in Figure 3. We observe again that the SNK and SNK2 solvers seem quite insensitive to the diffusion strength, while the number of outer iterations in NKS increases mildly as diffusion wanes. Furthermore, SNK2 converges in about half as many outer iterations as SNK.

4.3 Parallel efficiency

A fully parallel-aware implementation of the methods would presumably distribute all the data and solving steps across cores, which would handle communication of interface values with neighbors when necessary. A simpler step was to modify our serial MATLAB implementation to use the parfor capability of the Parallel Computing Toolbox for the most compute-intensive loop in the SNK methods, that for the independent local nonlinear solves in the SNK residual evaluation, followed by factoring the final Jacobian matrices of these solutions. While we also tried parallelizing the loop for applying the inverses of the local Jacobians as part of the full Jacobian application in an inner Krylov iteration, the effect on timings was minimal or even detrimental due to the greater importance of communication relative to that computation.

The parallel SNK2 solver was applied to the Burgers experiment as described in section 4.2, but with the domain $\Omega$ split into an 8-by-8 array of uniformly sized subdomains, resulting in 3 outer iterations and a total of 18 inner iterations. The code was run on a compute node equipped with two 18C Intel E5-2695 v4 (for 36 total cores), 32 GB of DDR4 memory, and a 100 Gbps Intel OmniPath cluster network. The timing results for different numbers of parallel computing cores are given in Table 1. There is a good amount of speedup in evaluations of the nonlinear SNK residuals, consisting mainly of the solution of local nonlinear problems, which dominate the computing time for a small number of cores. However, the parallel efficiency is limited by the other parts of the implementation, most notably the Jacobian evaluations.

5 Discussion

We have described a framework for overlapping domain decomposition in which overlap regions are discretized independently by the local subdomains, even in the formulation of the global problem. Communication between subdomains occurs only via interpolation of values to interface points. This formulation makes it straightforward to apply high-order or spectral discretization methods in the subdomains and to adaptively refine them.

The technique may be applied to precondition a linearized PDE, but it may also be used to precondition the nonlinear problem before linearization, to get what we call the Schwarz–Newton–Krylov (SNK) technique. In doing so, one gets the same benefit of faster Krylov inner convergence, but the resulting nonlinear problem is demonstrably easier to solve in terms of outer iterations and robustness. Although we have not given the derivation here, the Jacobian of the preconditioned nonlinear problem is readily shown to be a low-rank perturbation of the identity. Thus Kantorovich or other standard convergence theory for Newton’s method [5] may therefore suggest improved local convergence rates and larger basin of attraction. We have not yet pursued this analysis.

We have demonstrated that the SNK method can easily be part of a two-level Full Approximation Scheme in order to keep iteration counts from growing as the number of subdomains grows. The coarse level is simply a coarsening on each subdomain, so that restriction and prolongation steps can be done simply and in parallel. Indeed, the situation should make a fully multilevel implementation straightforward, as the multilevel coarsenings and refinements can all be done within subdomains.

The most time-consuming part of the SNK algorithm is expected to be typically in the solution of nonlinear PDE problems within each subdomain using given boundary data. These compute-intensive tasks require no communication and are therefore efficient to parallelize. By contrast, each inner Krylov iteration (i.e., Jacobian application) of both SNK and linearly preconditioned NKS requires a communication of interface data between overlapping subdomains, which appears to generate a more communication-bound form of parallelism. An additional feature of the SNK approach, mentioned also in [3], is that subdomains of low solution activity can be expected to be found relatively quickly. We observed this to be the case in the cavity flow problem of section 4.1, where local solutions in regions of low activity were sometimes 3-4 times faster as those in regions with steep solution gradients. This presents a natural way to limit the spatial scope of difficult nonlinear problems, though it also raises questions for load balancing in a parallel environment.

Finally, we remark that an important extension in [1] is to use least-squares approximation rather than interpolation to incorporate nonrectangular (sub)domains. We have been able to write a least-squares (as opposed to collocation) generalization of SNK and test it in one dimension. We hope to make it the subject of future work.