Linear and nonlinear substructured Restricted Additive Schwarz iterations and preconditioning

Iterative substructuring Domain Decomposition (DD) methods have been extensively studied, and they are usually associated with nonoverlapping decompositions. It is less known that classical overlapping DD methods can also be formulated in substructured form, i.e., as iterative methods acting on variables defined exclusively on the interfaces of the overlapping domain decomposition. We call such formulations substructured domain decomposition methods. We introduce here a substructured version of Restricted Additive Schwarz (RAS) which we call SRAS. We show that RAS and SRAS are equivalent when used as iterative solvers, as they produce the same iterates, while they are substantially different when used as preconditioners for GMRES. We link the volume and substructured Krylov spaces and show that the iterates are different by deriving the least squares problems solved at each GMRES iteration. When used as iterative solvers, SRAS presents computational advantages over RAS, as it avoids computations with matrices and vectors at the volume level. When used as preconditioners, SRAS has the further advantage of allowing GMRES to store smaller vectors and perform orthogonalization in a lower dimensional space. We then consider nonlinear problems, and we introduce SRASPEN (Substructured Restricted Additive Schwarz Preconditioned Exact Newton), where SRAS is used as a preconditioner for Newton’s method. In contrast to the linear case, we prove that Newton’s method applied to the preconditioned volume and substructured formulation produces the same iterates in the nonlinear case. Next, we introduce two-level versions of nonlinear SRAS and SRASPEN. Finally, we validate our theoretical results with numerical experiments.


Introduction
We consider a boundary value problem posed in a Lipschitz domain We assume that (1) admits a unique solution in some Hilbert space V.If the boundary value problem is linear, a discretization of (1) with N v degrees of freedom leads to a linear system where A ∈ R Nv×Nv , u ∈ V ( ∼ = R Nv ), and f ∈ V .If the boundary value problem is nonlinear, we obtain a nonlinear system where F : V → V is a nonlinear function and u ∈ V .Several numerical methods have been proposed in the last decades for the efficient solution of such boundary value problems, e.g., multigrid methods [21,37] and domain decomposition (DD) methods [36,34].We will focus on DD methods, which are usually divided into two distinct classes, that is overlapping methods, which include the AS (Additive Schwarz) and RAS (Restricted Additive Schwarz) methods [36,6], and nonoverlapping methods such as FETI (Finite Element Tearing and Interconnect) and Neumann-Neumann methods [14,32,26].Concerning nonlinear problems, DD methods can be applied either as nonlinear iterative methods, that is by just solving nonlinear problems in each subdomain and then exchanging information between subdomains as in the linear case [29,30,2], or as preconditioners to solve the Jacobian linear system inside Newton's iteration.In the latter case, the term Newton-Krylov-DD is employed, where DD is replaced by the domain decomposition preconditioner used [23].
An alternative is to use a DD method as a preconditioner for Newton's method.Preconditioning a nonlinear system F (u) = 0 means that we aim to replace the original nonlinear system with a new nonlinear system, still having the same solution, but for which the nonlinearities are more balanced and Newton's method converges faster [3,17].Seminal contributions in nonlinear preconditioning have been made by Cai and Keyes in [3,4], where they introduced ASPIN (Additive Schwarz Preconditioned Inexact Newton).The development of good preconditioners is not an easy task even in the linear case.One useful strategy is to study efficient iterative methods, and then to use the associated preconditioners in combination with Krylov methods [17].The same logical path paved the way to the development of RASPEN (Restricted Additive Schwarz Preconditioned Exact Newton) in [13], which in short applies Newton's method to the fixed point equation defined by the nonlinear RAS iteration at convergence.Extensions of this idea to Dirichlet-Neumann are presented in [7].In [22], the authors describe and analyze the scalability of the two-level variants of the aformentioned methods (ASPIN and RASPEN).In particular, they discuss several approaches of adding the coarse space correction, and a numerical comparison of all these methods is reported for different types of coarse spaces.All these methods are left preconditioners.
Right preconditioners are usually based on the concept of nonlinear elimination, presented in [27], and they are very efficient as shown in [20,19,5,31].
While left preconditioners aim to transform the original nonlinear function into a better behaved one, right preconditioners aim to provide a better initial guess for the next outer Newton iteration.
Nonoverlapping methods are sometimes called substructuring methods (a term borrowed from Przemieniecki's work [33]), as in these methods the unknowns in the interior of the nonoverlapping domains are eliminated through static condensation so that one needs to solve a smaller system involving only the degrees of freedom on the interfaces between the nonoverlapping subdomains [36].However, it is also possible to write an overlapping method, such as Lions' Parallel Schwarz Method (PSM) [28], which is equivalent to RAS [16], in substructured form, even though this approach is much less common in the literature.For a two subdomain decomposition, a substructuring procedure applied to the PSM is carried out in [15, Section 5], [18,Section 3.4] and [10].
In [10,11], the authors introduced a substructured formulation of the PSM at the continuous level for decompositions with many subdomains and crosspoints, and further studied ad-hoc spectral and geometric two-level methods.
In this particular framework, the substructured unknowns are now the degrees of freedom located on the portions of a subdomain boundary that lie in the interior of another subdomain; that is where the overlapping DD method takes the information to compute the new iterate.We emphasize that, at a given iteration n, any iterative DD method (overlapping or nonoverlapping) needs only a few values of u n to compute the new approximation u n+1 .The major part of u n is useless.
In this manuscript, we define a substructured version of RAS, that is we define an iterative scheme based on RAS which acts only over unknowns that are located on the portions of a subdomain boundary that lie in the interior of another subdomain.We study in detail the effects that such a substructuring procedure has on RAS when the latter is applied either as an iterative solver or as a preconditioner to solve linear and nonlinear boundary value problems.Does the substructured iterative version converge faster than the volume one?Is the convergence of GMRES affected by substructuring?What about nonlinear problems when instead of preconditioned GMRES we rely on preconditioned Newton?We prove that substructuring does not influence the convergence of the iterative methods both in the linear and nonlinear case, by showing that at each iteration, the restriction on the interfaces of the volume iterates coincides with the iterates of the substructured iterative method.
Nevertheless, we discuss in Section 3.1 and corroborate by numerical experiments that a substructured formulation presents computational advantages.
The equivalence of iterates does not hold anymore when considering preconditioned GMRES.Specifically, our study shows that GMRES should be applied to the substructured system, since it is computationally less expensive, requiring to perform orthogonalization on a much smaller space, and thus needs also less memory.In contrast to the linear case, we prove that, surprisingly, the nonlinear preconditioners RASPEN and SRASPEN (Substructured RASPEN) for Newton produce the same iterates once these are restricted to the interfaces.
However, SRASPEN has again more favorable properties when assembling and solving the Jacobian matrices at each Newton iteration.Finally, we also extend the work in [10,11] defining substructured two-level methods to the nonlinear case, where both smoother and coarse correction are defined directly on the interfaces between subdomains.
This paper is organized as follows: we introduce in Section 2 the mathematical setting with the domain, subdomains and operators defined on them.In Section 3, devoted to the linear case, we study the effects of substructuring on RAS and on GMRES applied to the preconditioned system.In Section 4, we extend our analysis to nonlinear boundary value problems.Section 5 contains two-level substructured methods for the nonlinear problems.Finally Section 6 presents numerical tests to corroborate the framework proposed.

Notation
Let us decompose the domain Ω into N nonoverlapping subdomains Ω j , that is Ω = j∈J Ω j with J := {1, 2, . . ., N }.The nonoverlapping subdomains Ω j are then enlarged to obtain subdomains Ω j which form an overlapping decomposition of Ω.For each subdomain Ω j , we define V j as the restriction of V to Ω j , that is V j collects the degrees of freedom on Ω j .Further, we introduce the classical restriction and prolongation operators R j : V → V j , P j : V j → V , and the restricted prolongation operators P j : V j → V .We assume that these operators satisfy R j P j = I Vj , and j∈J Fig. 1: The domain Ω is divided into nine nonoverlapping subdomains (left).
The center panel shows how the diagonal nonoverlapping subdomains are enlarged to form overlapping subdomains.On the right, we denote the unknowns represented in V (blue line) and the unknowns of a coarse space of V (red crosses).
where I Vj is the identity on V j and I is the identity on V .
We now define the substructured skeleton.In the following, we use the notation introduced in [9].For any j ∈ J , we define the set of neighboring indices N j := { ∈ J : Ω j ∩ ∂Ω = ∅}.Given a j ∈ J , we introduce the substructure of Ω j defined as S j := ∈Nj ∂Ω ∩ Ω j , that is the union of all the portions of ∂Ω with ∈ N j .The substructure of the whole domain Ω is defined as S := j∈J S j .A graphical representation of S is given in Fig. 1 for a decomposition of a square into nine subdomains.We now introduce the space V as the trace space of V onto the substructure S. Associated to V , we consider the restriction operator R : V → V and a prolongation operator The restriction operator R takes an element v ∈ V and restricts it to the skeleton S. The prolongation operator P extends an element v ∈ V to the global space V .In our numerical experiments, P extends an element v S ∈ V by zero in Ω\S.However, we can consider several different prolongation operators.How this extension is done is not crucial as we will use P inside a domain decomposition algorithm, and thus only the values on the skeleton S will play a role.Hence, as of now, we will need only one assumption on the restriction and prolongation operator, namely where I is the identity over V .

The linear case
In this section, we focus on the linear problem Au = f .After defining a substructured variant of RAS called SRAS, we prove the equivalence between RAS and SRAS.Then, we study in detail how GMRES performs if applied to the volume preconditioned system or the substructured system.

Linear iterative methods
To introduce our analysis, we recall the classical definition of RAS to solve the linear system (2).RAS starts from an approximation u 0 and computes for n = 1, 2, . . ., where A j := R j AP j , that is, we use exact local solvers.Let us now rewrite the iteration (6) in an equivalent form using the hypothesis in (4) and the definition of A j , We emphasize that (P j R j − I) u n−1 contains non-zero elements only outside subdomain Ω j , and in particular the terms A (P j R j − I) u n−1 represent precisely the boundary conditions for Ω j given the old approximation u n−1 .This observation suggests that RAS, like most domain decomposition methods, can be written in a substructured form.Indeed, despite iteration (7) being written in volume form, involving the entire vector u n−1 , only very few elements of u n−1 are needed to compute the new approximation u n .A substructured method iterates only on those values of u which are really needed at each iteration, avoiding thus superfluous operations on the whole volume vector u (e.g. the volume residual computation f − Au n−1 and the summation with the old iterate u n−1 in the RAS method (6)).For further details about a substructured formulation of the parallel Schwarz method at the continuous level, we refer to [18] for the two subdomain case, and [10,11] for a general decomposition into several subdomains with crosspoints.
In Section 2, we introduced the substructured space V geometrically, but we can also provide an algebraic characterization using the RAS operators R j and P j .We consider that is, K is the set of indices such that the canonical vectors e k represent a Dirichlet boundary condition at least for a subdomain, and its complement The cardinality of K is |K| =: N .We can thus introduce Finally R is the Boolean restriction operator, mapping a vector of R Nv onto a vector of R N , keeping only the indices in K. Hence, V := ImR( ∼ = R N ) and To define SRAS, we need one more assumption on the restriction and prolongation operators, namely where M −1 is the preconditioner for RAS, formally defined as Heuristically, this assumption means that the operator P R preserves all the information needed by G RAS (defined in (7)) to compute correctly the values of the new iterate on the skeleton S. Indeed a direct calculation shows that ( 8) is equivalent to the condition Given a substructured approximation v 0 ∈ V , for n = 1, 2, . . ., we define SRAS as RAS and SRAS are tightly linked, but when are they equivalent?Clearly, we must impose some conditions on P and R. The next theorem shows that assumption (8) is in fact sufficient for equivalence.
Theorem 1 (Equivalence between RAS and SRAS) Assume that the operators R and P satisfy (8).Given an initial guess u 0 ∈ V and its substructured restriction v 0 := Ru 0 ∈ V , define the sequences {u n } and {v n } such that Then, Ru n = v n for every iteration n ≥ 1.
Proof We prove this statement for n = 1 by a direct calculation.Taking the restriction of u 1 we have where we used assumption (8), and the definition of v 0 and G SRAS .For a general n, the proof is obtained by induction.
Remark 1 (Implementation of SRAS) In ( 10), we have introduced the substructured operator G SRAS directly through the volume operator G RAS .This definition is very useful from the theoretical point of view as it permits to link the volume and substructured methods and facilitates the theoretical analysis.
However, we stress that one should not implement G SRAS by directly calling the volume routine G RAS onto a vector P v. Doing so, one would lose all the computational advantages as computations on the volume vector P v would be performed.There are two strategies to implement a fully substructured SRAS method.The first one is to implement a routine that, for each j = 1, . . ., N , extracts from v those values which lie on the boundary of Ω j and rescale them appropriately as boundary conditions for the local subdomain solve.A second possibility arises from ( 10) and ( 7), by observing that where b := R N j=1 P j A −1 j R j f .One can then pre-assemble the matrices , where M j is the number of degrees of freedom in Ω j .The matrices R j are very sparse and their goal is just to extract the values needed for the j-th subdomain solve from v n−1 .Similarly weights a subdomain solution with the partition of unity and maps it to the substructured vector.We finally obtain the equivalent iterative method where no computations are performed with matrices or vectors at the volume level, except for the subdomain solves.

Linear preconditioners for GMRES
It is well known that any stationary iterative method should be used in practice as a preconditioner for a Krylov method, since the Krylov method finds in general a much better residual polynomial with certain optimality properties, compared to the residual polynomial of the stationary iteration, see e.g.[8].The preconditioner associated to RAS is M −1 and is defined in (9).The preconditioned volume system then reads To discover the preconditioner associated with SRAS, we consider the fixed point limit of (10), where in the second line we used the identity RP = I.We thus consider the preconditioned substructured system Observe that as A = M − N , ( 15) can be written as (I − G)v = b, where G := RM −1 N P and b := RM −1 f , thus recovering the classical form of the substructured PSM, see [18,10,11].
It is then natural to ask how a Krylov method like GMRES performs if applied to (13), compared to (15).Let us consider an initial guess in volume u 0 , its restriction v 0 := Ru 0 and the initial residuals Then GMRES applied to the preconditioned systems ( 13) and ( 15) looks for solutions in the affine Krylov spaces where k ≥ 1.The two Krylov spaces are tightly linked, as Theorem 2 below will show.To prove it, we need the following Lemma.
Lemma 1 If the restriction and prolongation operators R and P satisfy (8) Proof Multiplying equation ( 8) from the right by M −1 A we get where in the second equality we have used once more (8).Using induction, one gets for every k ≥ 1, and this completes the proof.

Theorem 2 (Relation between RAS and SRAS Krylov subspaces)
Let us consider operators R and P satisfying (8), an initial guess u 0 ∈ V , its restriction v 0 := Ru 0 ∈ V and the residuals r Then for every k ≥ 1, we have Proof First, due to (8) we have Let us now show the first inclusion.
) and we achieve the desired relation (18).
Theorem 2 shows that the restriction to the substructure of the affine volume Krylov space of RAS coincides with the affine substructured Krylov space of SRAS.One could then wonder if the restrictions of the iterates of GMRES applied to the preconditioned volume system (13) coincide with the iterates of GMRES applied to the preconditioned substructured system (15).
However, this does not turn out to be true.Nevertheless, we can further link the action of GMRES on these two preconditioned systems.
It is well known, (see e.g [35, Section 6.5.1]), that GMRES applied to (13) and ( 15) generates a sequence of iterates u k k and v k k such that and The iterates u k and v k can be characterized using orthogonal and Hessenberg matrices obtained with the Arnoldi iteration.In particular, the k-th iteration of Arnoldi provides orthogonal matrices Q k , Q k+1 and a Hessenberg matrix and the columns of Q k form an orthonormal basis for the Krylov subspace K k (M −1 A, r 0 ).Using these matrices, one writes and e 1 is the canonical vector of R k+1 .Similarly, one characterizes the vector where Q k , H k are the orthogonal and Hessenberg matrices obtained through the Arnoldi method applied to the matrix RM −1 AP .
The next theorem provides a link between the volume least square problem (19) and the substructured one (20).
Theorem 3 Under the hypothesis of Theorem 2, the k-th iterate of GMRES applied to (15) is equal to , where y satisfies (22 while Proof It is clear that Thus the columns of RQ k form an orthonormal basis of K k (RM −1 AP , r 0 ) and hence, v k can be expressed as a linear combination of the columns of RQ k with coefficients in the vector t ∈ R k plus v 0 .We are then left to show (23).We have min Using the relation Im(RQ k ) =Im(Q k ), Lemma (1), the Arnoldi relation H k and that r 0 coincides with the first column of Q k except for a normalization constant, we conclude min and this completes the proof.
Few comments are in order here.First, GMRES applied to (15) converges in maximum N iterations as the preconditioned matrix RM −1 AP has size N ×N .
Second, Theorem (2) states that R(u 0 + K N (M −1 A, r 0 )) already contains the exact substructured solution, that is the exact substructured solution lies in the restriction of the volume Krylov space after N iterations.Theoretically, if one could get the exact substructured solution from R(u then N iterations of GMRES applied to (13), plus an harmonic extension of the substructured data into the subdomains, would be sufficient to get the exact volume solution.
On the other hand, we can say a bit more analyzing the structure of M −1 A.

Using the splitting
that is the Krylov space generated by M −1 A is equal to the Krylov space generated by the RAS iteration matrix for error equation.We denote this linear operator with G RAS 0 which is defined as in (7) with f = 0. We now consider the orthogonal complement V ⊥ := (span {e k } k∈K ) ⊥ = span {e i } i∈K c , and Using the rank-nullity theorem, we obtain dim Im(G RAS 0 hence GMRES applied to the preconditioned volume system encounters a lucky Arnoldi breakdown after at most N + 1 iterations (in exact arithmetic).This rank argument can be used for the substructured preconditioned system as well.Indeed as RM −1 AP = I − RM −1 N P , the substructured Krylov space is generated by the matrix RM −1 N P , whose rank is equal to the rank of Heuristically, choosing a zero initial guess, r 0 := M −1 f corresponds to a solution of subdomains problem with the correct right hand side, but with zero Dirichlet boundary conditions along the interfaces of each subdomain.Thus, GMRES applied to (13) needs only to find the correct boundary conditions for each subdomain, and this can be achieved in at most N iterations as Theorem 2 shows.
Finally, we remark that each GMRES iteration on ( 15) is computationally less expensive than a GMRES iteration on (13) as the orthogonalization of the Arnoldi method is carried out in a much smaller space.From the memory point of view, this implies that GMRES needs to store shorter vectors.Thus, a saturation of the memory is less likely, and restarted versions of GMRES may be avoided.

The nonlinear case
In this section, we study iterative and preconditioned domain decomposition methods to solve the nonlinear system (3).

Nonlinear iterative methods
RAS can be generalized to solve the nonlinear equation (3).To show this, we introduce the solution operators G j which are defined through where the operators R j and P j are defined in Section 2. Nonlinear RAS for N subdomains then reads It is possible to show that (26) reduces to (7) if F (u) is a linear function: assuming that F (u) = Au − f , equation ( 25) becomes , and thus (26) reduces to (7).
Similarly to the linear case, we introduce the nonlinear SRAS.Defining we obtain the nonlinear substructured iteration which is the nonlinear counterpart of (10).
The same calculations of Theorem 1 allow one to obtain an equivalence result between nonlinear RAS and nonlinear SRAS.
Theorem 4 (Equivalence between nonlinear RAS and SRAS) Assume that the operators R and P satisfy R j∈J P j G j (u) = R j∈J P j G j (P Ru).
Let us consider an initial guess u 0 ∈ V and its substructured restriction v 0 := Ru 0 ∈ V , and define the sequences {u n }, {v n } such that Then for every n ≥ 1, Ru n = v n .

Nonlinear preconditioners for Newton's method
In [13], it was proposed to use the fixed point equation of nonlinear RAS as a preconditioner for Newton's method, in a spirit that goes back to [4,3].
This method has been called RASPEN (Restricted Additive Schwarz Preconditioned Exact Newton) and it consists in applying Newton's method to the fixed point equation of nonlinear RAS, that is, For a comprehensive discussion of this method, we refer to [13].As done in (14) for the linear case, we now introduce a substructured variant of RASPEN and we call it SRASPEN (Substructured Restricted Additive Schwarz Preconditioned Exact Newton).SRASPEN is obtained by applying Newton's method to the fixed point equation of nonlinear SRAS, that is, One can verify that the above equation F(v) = 0 can also be written as This formulation of SRASPEN provides its relation with RASPEN and simplifies the task of computing the Jacobian of SRASPEN.

Computation of the Jacobian and implementation details
To apply Newton's method, we need to compute the Jacobian of SRASPEN.
Let J F (w) and J F (w) denote the action of the Jacobian of RASPEN and SRASPEN on a vector w.Since these methods are closely related, indeed F(v) = RF(P v), we can immediately compute the Jacobian of F once we have the Jacobian of F, using the chain rule, J F (v) = RJ F (P v)P .The Jacobian of F has been derived in [13] and we report here the main steps for the sake of completeness.Differentiating equation ( 29) with respect to u leads to Recall that the local inverse operators G j : V → V j are defined in equation (25) as the solutions of R j F (P j G j (u) + (I − P j R j )u) = 0. Differentiating this relation yields where u (j) := P j G j (u) + (I − P j R j )u is the volume solution vector in subdomain j and J is the Jacobian of the original nonlinear function F .Combining the above equations ( 31)-( 32) and defining u (j) := P j G j (P v) + (I − P j R j )P v, we get and where we used the assumptions j∈J P j R j = I and RP = I.We remark that to assemble J F (u) or to compute its action on a given vector, one needs to calculate J u (j) , that is, evaluate the Jacobian of the original nonlinear function F on the subdomain solutions u (j) .The subdomain solutions u (j)     are obtained evaluating F(u), that is performing one step of RAS with initial guess equal to u.A smart implementation can use the local Jacobian matrices R j J u (j) P j that are already computed by the inner Newton solvers while solving the nonlinear problem on each subdomain, and hence no extra cost is required to assemble this term.Further, the matrices R j J u (j) are different from the local Jacobian matrices at very few columns corresponding to the degrees of freedom on the interfaces and thus it suffices to only modify those specific entries.In a non-optimized implementation, one can also directly evaluate the Jacobian of F on the subdomain solutions u (j) , without relying on already computed quantities.Concerning J F (v), we emphasize that u (j) is the volume subdomain solution obtained by substructured RAS starting from a substructured function v. Thus, like u (j) , u (j) is readily available in Newton iteration after evaluating the function F.
From the computational point of view, SRASPEN has several advantages over RASPEN.From ( 33) and (34), we note that J F is a matrix of dimension N × N where N is the number of unknowns on S, and thus is a much smaller matrix than J F , whose size is N v × N v , with N v the number of unknowns in volume.On the one hand, if one prefers to assemble the Jacobian matrix, either because one wants to use a direct solver or because one wants to recycle the Jacobian for several iterations, then SRASPEN dramatically reduces the cost of the assembly of the Jacobian matrix.On the other hand, we remark that (33) and ( 34) have the same structure of the volume and substructured preconditioned matrices ( 13) and (15), by just identifying M −1 = j∈J P j R j J u (j) P j −1 .Similarly to Remark 1 in the linear case, we can have a fully substructured formulation, by writing where P j := R P j and R j := R j J( u j )P .It follows that if one prefers to use a Krylov method such as GMRES, then according to the discussion in Section 3.2, SRASPEN better exploits the properties of the underlying domain decomposition method, and saves computational time by permitting to perform the orthogonalization in a much smaller space.Further implementation details and a more extensive comparison are available in the numerical section 6.

Convergence analysis of RASPEN and SRASPEN
Theorem 4 gives an equivalence between nonlinear RAS and nonlinear SRAS.
Are RASPEN and SRASPEN equivalent?Does Newton's method behave differently if applied to the volume or to the substructured fixed point equation, like it happens with GMRES (see Section 3.2)?In this section, we aim to answer these questions by discussing the convergence properties of the exact Newton's method applied to F and F.
Let us recall that, given two approximations u 0 and v 0 , the exact Newton's method computes for n ≥ 1, where J F u n−1 and J F v n−1 are the Jacobian matrices respectively of F and F evaluated at u n−1 and v n−1 .In this paragraph, we do not need a precise expression for J F and J F .However we recall that, the definition F(v) = RF(P v) and the chain rule derivation provides us the relation J F (v) = RJ F (P v)P .If the operators R and P were square matrices, we would immediately obtain that RASPEN and SRASPEN are equivalent, due to the affine invariance theory for Newton's method [12].However, in our case, R and P are rectangular matrices and they map between spaces of different dimensions.Nevertheless, in the following theorem, we show that RASPEN and SRASPEN provide the same iterates restricted to the interfaces under further assumptions on R and P , which is a direct generalization of (8) to the nonlinear case.
Theorem 5 (Equivalence between RASPEN and SRASPEN) Assume that the operators R and P satisfy Given an initial guess u 0 ∈ V and its substructured restriction v 0 := Ru 0 ∈ V , define the sequences {u n } and {v n } such that Then for every n ≥ 1, Proof We first prove the equality Ru 1 = v 1 by direct calculations.Taking the restriction of the RASPEN iteration, we obtain Now, due to the definition of F and of v 0 , and to the assumption (35), we have Further, taking the Jacobian of assumption (35), we have RJ F (u 0 ) = J F (Ru 0 )R, which simplifies by taking the inverse of the Jacobians to Finally substituting relations (37) and ( 38) into (36) leads to and the general case is obtained by induction.
5 Two-level nonlinear methods RAS and SRAS can be generalized to two-level iterative schemes.This has already been treated in detail for the linear case in [10,11].In this section, we introduce two-level variants for nonlinear RAS and SRAS, and also for the associated RASPEN and SRASPEN.

Two-Level iterative methods
To define a two-level method, we introduce a coarse space V 0 ⊂ V , a restriction operator R 0 : V → V 0 and an interpolation operator P 0 : V 0 → V .The nonlinear system F can be projected onto the coarse space V 0 , defining the coarse nonlinear function F 0 (u 0 ) := R 0 F (P 0 u 0 ), for every u 0 ∈ V 0 .Due to this definition, it follows immediately that J F0 (u 0 ) = R 0 J F (P 0 u 0 ) P 0 , ∀u 0 ∈ V 0 .To compute a coarse correction we rely on the FAS approach [1].Given a current approximation u, the coarse correction C 0 (u) is computed as the solution of Two-level nonlinear RAS is described by Algorithm 1 and it consists of a coarse correction followed by one iteration of nonlinear RAS (see [22] for different approaches).
We now focus on its substructured counterpart.We introduce a coarse substructured space V 0 ⊂ V , a restriction operator R 0 : V → V 0 and a prolongation operator P 0 : V 0 → V .We define the coarse substructured function as From the definition it follows that There is a profound difference between two-level nonlinear RAS and two-level Algorithm 1 Two-level nonlinear RAS 2: Add the coarse correction to the current iterate, u k+ 1 2 = u k + P 0 C 0 u k .3: Compute one step of nonlinear RAS, u k+1 = j∈J P j G j u k+ 1 2 .4: Repeat steps 1 to 3 until convergence.
Algorithm 2 Two-level iterative nonlinear SRAS 2: Add the coarse correction to the current iterate, 2 .4: Repeat steps 1 to 3 until convergence.
nonlinear SRAS: in the first one (Algorithm 1), the coarse function is obtained restricting the original nonlinear system F (u) = 0 onto a coarse mesh.In the substructured version, the coarse substructured function is defined restricting the fixed point equation of nonlinear SRAS to V 0 .That is, the coarse substructured function corresponds to a coarse version of SRASPEN.Hence, we remark that this algorithm is the nonlinear counterpart of the linear two-level algorithm described in [10,11].Two-level nonlinear SRAS is then defined in Algorithm 2. As in the linear case, numerical experiments will show that twolevel iterative nonlinear SRAS exhibits faster convergence in terms of iteration counts compared to two-level nonlinear RAS.However, we remark that evaluating F 0 is rather cheap, while evaluating F 0 could be quite expensive as it requires to perform subdomain solves on the fine mesh.One possible improvement is to approximate F 0 replacing F in its definition with another function which performs subdomain solves on a coarse mesh.Further, we emphasize that a prerequisite of any domain decomposition method is that the subdomain solves are cheap to compute in a high performance parallel implementation, so that in such a setting evaluating F 0 needs to be cheap as well.

Two-level preconditioners for Newton's method
Once we have defined the two-level iterative methods, we are ready to introduce the two-level versions of RASPEN and SRASPEN.The fixed point equation of two-level nonlinear RAS is where we have introduced the correction operators C j (u) := G j (u) − R j u.
Thus two-level RASPEN defined in [13] consists in applying Newton's method to the fixed point equation (41).
Similarly, the fixed point equation of two-level nonlinear SRAS is where the correction operators C j are defined as Two-level SRASPEN consists in applying Newton's method to the fixed point equation (42).

Numerical results
We discuss three different examples in this section to illustrate our theoretical results.In the first example, we consider a linear problem where we study the performance of GMRES when it is applied to the preconditioned volume system and the preconditioned substructured system.In the next two examples, we present numerical results in order to compare Newton's method, NKRAS [5], nonlinear RAS, nonlinear SRAS, RASPEN, and SRASPEN for the solution of a one-dimensional Forchheimer equation and for a two-dimensional nonlinear diffusion equation.

Linear example
We consider the diffusion equation −∆u = f , with source term f ≡ 1 and homogeneous boundary conditions inside the unit cube Ω := (0, 1) 3 decomposed into N equally-sized bricks with overlap, each discretized with 27000 degrees of freedom.The size of the overlap is δ := 4 × h.In Table 1, we study the computational effort and memory required by GMRES when applied to the preconditioned volume system (13) (GMRES-RAS) and to the preconditioned substructured system (15) (GMRES-SRAS).We let the number of subdomains grow, while keeping their sizes constant, that is the global problem becomes larger as N increases.We report the computational times to reach a relative residual smaller than 10 −8 , and the number of gigabytes required to store the orthogonal matrices of the Arnoldi iteration, both for the volume and substructured implementations.The subdomain solves are performed in a serial fashion, we precompute the Cholesky factorizations for the subdomain matrices A j , and for SRAS we use the fixed point equation related to formulation (12).
Table 1 shows that GMRES applied to the preconditioned substructured system is faster in terms of computational time compared to the volume implementation.This advantage becomes more evident as the global problem becomes larger.We emphasize that GMRES required the same number of iterations to reach the tolerance for both methods in all cases considered.Thus the faster time to solution of GMRES-SRAS is due to the smaller number of floating point operations that GMRES-SRAS has to perform, since the orthogonalization steps are performed in a much smaller space, and SRAS avoids unnecessary volume computations.Furthermore, GMRES-SRAS significantly outperforms GMRES-RAS in terms of memory requirements; in this particular case, GMRES-SRAS computes and stores orthogonal matrices which are about seven times smaller than the ones used by GMRES-RAS.

Forchheimer equation in 1D
Forchheimer equation is an extension of the Darcy equation for high flow rates, where the linear relation between the flow velocity and the gradient flow does not hold anymore.In a one dimensional domain Ω := (0, 1), the Forchheimer model is where u L , u R ∈ R, λ(x) is a positive and bounded permeability field and q(y) := sign(y) , with γ > 0. To discretize (43), we use the finite volume scheme described in detail in [13].In our numerical experiments, we set λ(x) = 2 + cos(5πx), f (x) = 50 sin(5πx)e x , γ = 1, u(0) = 1 and u(1) = e 1 .
The solution field u(x) and the force field f (x) are shown in Fig. 2. We then study the convergence behavior of our different methods.Fig. 3 shows how the relative error decays for the different methods and for a decomposition into 20 subdomains (left panel) and 50 subdomains (right panel).The initial guess is equal to zero for all these methods.Both plots in Fig. 3 show that the convergence rate of iterative nonlinear RAS and nonlinear SRAS is the same and very slow.As expected, NKRAS with line search converges better than Newton's method and further RASPEN and SRASPEN converge in the same number of outer Newton iterations as they produce the same iterates.
Moreover, it seems that the convergence of RASPEN and SRASPEN is not affected by the number of subdomains.However, these plots do not tell the whole story, as one should focus not only on the number of iterations but also  on the cost of each iteration.To compare the cost of an iteration of RASPEN and SRASPEN, we have to distinguish two cases, that is, if one solves the Jacobian system directly or with some Krylov methods, e.g., GMRES.First, suppose that we want to solve the Jacobian system with a direct method and thus we need to assemble and store the Jacobians.From the expressions in equation (34) we remark that the assembly of the Jacobian of RASPEN requires N × N v subdomain solves, where N is the number of subdomains and N v is the number of unknowns in volume.On the other hand, the assembly of the Jacobian of SRASPEN requires N × N solves, where N is the number of unknowns on the substructures and N N v .Thus, while the assembly of J F is prohibitive, it can still be affordable to assemble J F .Further, the direct solution of the Jacobian system is feasible as J F has size N × N .Suppose now N steps, which is the size of the substructured Jacobian.Under these circumstances, it can be convenient to actually assemble J F , as it requires N × N subdomain solves which is the total cost of GMRES.Furthermore, the N × N subdomain solves are embarrassingly parallel, while the N × N solves of GM-RES can be parallelized in the spatial direction, but not in the iterative one.
As future work, we believe it will be interesting to study the convergence of a Quasi-Newton method based on SRASPEN, where one assembles the Jacobian substructured matrix after every few outer Newton iterations, reducing the overall computational cost.
As a final remark, we specify that Fig. 4 has been obtained setting a zero initial guess for the nonlinear subdomain problems.However, at the iteration k of RASPEN one can use the subdomain restriction of the updated volume solution, that is R j u k−1 , which has been obtained by solving the volume Jacobian system at iteration k − 1, and is thus generally a better initial guess for the next iteration.On the other hand in SRASPEN, one could use the subdomain solutions computed at iteration k − 1, i.e. u k−1 i , as initial guess for the nonlinear subdomain problems, as the substructured Jacobian system corrects only the substructured values.Numerical experiments showed that with this particular choice of initial guess for the nonlinear subdomain problems, SRASPEN requires generally more Newton iterations to solve the local problems.In this setting, there is not a method that is constantly faster than the other as it depends on a delicate trade-off between the better GMRES performance and the need to perform more Newton iterations for the nonlinear local problems in SRASPEN.

Nonlinear Diffusion
In this subsection we consider the nonlinear diffusion problem on a square domain where the right hand side f is chosen such that u(x) = sin(πx) sin(πy) is the exact solution.We start all these methods with an initial guess u 0 (x) = 10 5 , so that we start far away from the exact solution, and hence Newton's method exhibits a long plateau before quadratic convergence begins.Fig. 5 shows the convergence behavior for the different methods as function of the number of iterations and the number of linear solves.The average number of GMRES iterations is 8.1667 for both RASPEN and SRASPEN for the four subdomain decomposition.For a decomposition into 25 subdomains, the average number of GMRES iterations is 19.14 for RASPEN and 19.57for SRASPEN.We remark that as the number of subdomains increases, GMRES needs more iterations to solve the Jacobian system.This is consistent with the interpretation of (34) as a Jacobian matrix J u (j) preconditioned by the additive operator j∈J R j J u (j) P j −1 ; We expect this preconditioner not to be scalable since it does not involve a coarse correction.In Table 2 we compare the computational time in seconds to reach a tolerance of 10 −8 by RASPEN and SRASPEN.SRASPEN is faster due to the less expensive GMRES iteration which is inherited by the linear analysis, see Table 1.We conclude this section by showing the convergence behavior for the twolevel variants of nonlinear RAS, nonlinear SRAS, RASPEN, and SRASPEN.
We use a coarse grid in volume taking half of the points in x and y, and a coarse substructured grid taking half of the unknowns as depicted in Fig. 1.
The interpolation and restriction operators P 0 , R 0 , P 0 and R 0 are the classical linear interpolation and fully weighting restriction operators defined in Section 5. From Fig. 6, we note that two-level nonlinear SRAS is much faster than twolevel nonlinear RAS, and this observation is in agreement with the linear case treated in [10,11].Since the two-level iterative methods are not equivalent, we also remark that two-level SRASPEN shows a better performance than two-level RASPEN in terms of iteration count.As the one-level smoother is the same in all methods, the better convergence of the substructured methods implies that the coarse equation involving F 0 provides a much better coarse correction than the classical volume one involving F 0 .
Even though the two-level substructured methods are faster in terms of iteration count, the solution of the FAS problem involving F 0 = R 0 F(P 0 (v 0 )) is rather expensive as it requires to evaluate twice the substructured function F (each evaluation requires subdomain solves) to compute the right hand side, to solve a Jacobian system involving J F 0 , and to evaluate F on the iterates, which again require the solution of subdomain problems.Unless one has a fully parallel implementation available, the coarse correction involving F 0 is doomed to represent a bottleneck.

Conclusions
We presented an analysis of the effects of substructuring on RAS when it is applied as an iterative solver and as a preconditioner.We proved that iterative RAS and iterative SRAS converge at the same rate, both in the linear and nonlinear case.For the nonlinear case, we showed that the preconditioned methods, namely RASPEN and SRASPEN, also have the same rate of convergence as they produce the same iterates once these are restricted to the interfaces.Surprisingly, the equivalence between volume and substructured RAS breaks down when they are considered as preconditioners for Krylov methods.We showed that the Krylov spaces are equivalent, once the volume one is restricted to the substructure, however we obtained that the iterates are different by carefully deriving the least squares problems solved by GMRES.Our analysis shows that GMRES should be applied to the substructured system as it converges similarly when applied to the volume formulation, but needs much less memory.This allows us to state that, while nonlinear RASPEN and SRASPEN produce the same iterates, SRASPEN has advantages when solving the Jacobian system, either because the use of a direct solve is feasible, or because the Krylov method can work at the substructured level.Finally, we introduced substructured two-level nonlinear SRAS and SRASPEN, and showed numerically that these methods have better convergence properties than their volume counterparts in terms of iteration count, although they are quite expensive in the present form per iteration.Future efforts will be in the direction of approximating F 0 , by replacing the function F, which is defined on a fine mesh, with an approximation on a very coarse mesh, thus reducing the overall cost of the substructured coarse correction, or by using spectral coarse spaces.

Fig. 3 :
Fig.3: Convergence behavior for Newton's method, NKRAS, nonlinear RAS, nonlinear SRAS, RASPEN and SRASPEN applied to Forchheimer equation.On the left, the simulation refers to a decomposition into 20 subdomains while on the right we consider 50 subdomains.The mesh size is h = 10 −3 and the overlap is 8h.

Fig. 4 :
Fig.4: Relative error decay for RASPEN and SRASPEN applied to Forchheimer equation with respect to the number of linear solves.On the left, the simulation refers to a decomposition into 20 subdomains while on the right we consider 50 subdomains.The mesh size is h = 10 −3 .

Fig. 5 :
Fig.5: Relative error decay versus the number of iterations (top row) and error decay versus the number of linear solves (bottom row).The left figures refer to a decomposition into four subdomains, the right figures to a decomposition into 25 subdomains.The mesh size is h = 0.012 and the overlap is 8h.

Fig. 6 :
Fig. 6: Relative error decay versus the number of iterations for Newton's method, iterative two-level nonlinear RAS and SRAS, and the two-level variants of RASPEN and SRASPEN.The left figure refers to a decomposition into 4 subdomains, while the right figure refers to a decomposition into 16 subdomains.The mesh size is h = 0.012 and the overlap is 4h.
as the first equality follows from standard GMRES literature (see e.g [35, Section 6.5.1]).The second equality follows from Theorem (2) as we have shown that v 0 +K k

Table 1 :
On the top, time in seconds required by GMRES-RAS and GMRES-SRAS to reach a relative error smaller than 10 −8 for increasingly larger problems.At the bottom, memory use expressed in gigabytes to store the Arnoldi orthogonal matrices in both GMRES implementations.

Table 2 :
Time in seconds required by RASPEN and SRASPEN to reach a relative error smaller than 10 −8 for the nonlinear diffusion equation with 4, 25 and 49 subdomains.