Nonlinear Monolithic Two-Level Schwarz Methods for the Navier-Stokes Equations
Axel Klawonn Martin Lanser



Nonlinear Monolithic Two-Level Schwarz Methods for the Navier-Stokes Equations
Nonlinear Two-Level Schwarz Methods for Navier-Stokes



Axel Klawonn^1,2, Martin Lanser^1,2 ^1Department of Mathematics and Computer Science, Division of Mathematics, University of Cologne, Weyertal 86-90, 50931 Cologne, Germany, axel.klawonn@uni-koeln.de, martin.lanser@uni-koeln.de, url: <https://www.numerik.uni-koeln.de> ^2Center for Data and Simulation Science, University of Cologne, Germany, url: <https://www.cds.uni-koeln.de>





*
    Axel Klawonn0000-0003-4765-7387 
 Martin Lanser0000-0002-4232-9395
    
======================================================================



Nonlinear domain decomposition methods became popular in recent years since they can improve the nonlinear convergence behavior of Newton's method significantly for many complex problems. In this article, a nonlinear two-level Schwarz approach is considered and, for the first time, equipped with monolithic GDSW (Generalized Dryja-Smith-Widlund) coarse basis functions for the Navier-Stokes equations. Results for lid-driven cavity problems with high Reynolds numbers are presented and compared with classical global Newton's method equipped with a linear Schwarz preconditioner. Different options, for example, local pressure corrections on the subdomain and recycling of coarse basis functions are discussed in the nonlinear Schwarz approach for the first time.



§ INTRODUCTION AND NONLINEAR PROBLEMS



In recent years, nonlinear preconditioning techniques became popular to improve the convergence speed of nonlinear solvers for complex partial differential equations. Many of these nonlinear preconditioners are based on the ideas and basic principles of linear domain decomposition methods (DDMs) and consequently they are denoted as nonlinear DDMs. In nonlinear DDMs, e.g., <cit.>, the discretized nonlinear problem is restricted to smaller local ones on subdomains and the local solutions of those are recombined to the original one in a Newton-like iteration, which hopefully converges faster than classical Newton's method applied to the original problem. As in the linear case it often is beneficial to include a (nonlinear) coarse space for robustness and convergence improvement, which leads to two-level nonlinear DDMs. In this article, we consider some variants of the nonlinear two-level Schwarz approaches first introduced in <cit.>; older two-level nonlinear Schwarz approaches for less general coarse spaces can be found in <cit.>. Here, for the first time, we combine two-level nonlinear Schwarz methods with a monolithic RGDSW (Reduced Generalized Dryja-Smith-Widlund) coarse space for Navier-Stokes problems (see <cit.>) and apply it to the well-known lid-driven cavity flow test problem considering high Reynolds numbers. 

We compare our solvers with classical Newton-Krylov-Schwarz approaches equipped with a simple globalization, namely a simple backtracking approach.   



  
Lid-driven Cavity Problem
As stated above, we consider the classical lid-driven cavity flow problem modeled by the stationary, dimensionless, and incompressible Navier-Stokes equations in two-dimensions

    [ -1/ReΔ v + (v ·∇)v + ∇ p                      = 0                     in Ω;                   div(v)                      = 0                     in Ω;                        v                    = v_0                    on ∂Ω;                        p                       =0         in  (x,y)=(0,0). ]

Here, v denotes the velocity, p the pressure, and Re the Reynolds number. In case of the lid-driven cavity problem Ω is the unity square and v_0=[1,0] for the Dirichlet boundary constraints on the lid, that is, the upper boundary of Ω where y is equal to one. On all other parts of the boundary ∂Ω we have v_0=[0,0]. We usually use v^(0)|_∂Ω=v_0 and v^(0)=0 within Ω and p^(0)=0 as initial value for Newton's method, regardless if we use a classical Newton-Krylov-Schwarz approach or a modern nonlinear Schwarz method. Therefore, since v^(0) already fulfills all boundary constraints, we can enforce a zero Dirichlet boundary condition on ∂Ω for the velocity of the Newton update in each Newton iteration. Note that for the pressure we enforce a Dirichlet boundary condition in a single corner, that is, in (x,y) = (0,0) to obtain a unique solution. 



§ NONLINEAR TWO-LEVEL SCHWARZ METHODS



Let us briefly describe the two-level nonlinear Schwarz methods considered in our experiments, which were first introduced in <cit.> and extended in <cit.> for more general coarse spaces. We assume to have a nonlinear problem 

    F(u)=0

which originates from a finite element discretization of a given nonlinear partial differential equation. In the following, we define nonlinearly left-preconditioned systems of the form

    ℱ(u) := G(F(u)) = 0,

where the left preconditioner G is given implicitly by a domain decomposition approach. The nonlinearly preconditioned problem is then solved with Newton's method instead of solving the original formulation from <ref> directly. The goal of this approach is to accelerate the nonlinear convergence of Newton's method by a proper choice of G. 

We first consider a decomposition of the computational domain Ω⊂ℝ^d,  d=2,3, into N nonoverlapping subdomains Ω_i',  i=1,...,N, such that Ω = ⋃_i=1^N Ω_i'.
By adding m rows of finite elements around the boundary of each subdomain, we obtain an overlapping domain decomposition Ω = ⋃_i=1^N Ω_i with overlap δ =m · h, where h is the typical diameter of a finite element. 

For the first level of the nonlinear Schwarz preconditioner, we define local nonlinear corrections T_i(u) by solving 

    R_i F(u- P_i T_i(u))=0,   i=1,...,N,

where R_i: V → V_i is the restriction from the global finite element space V to the local space V_i belonging to the overlapping subdomain Ω_i,   i=1,...,N. Here, P_i: V_i → V,   i=1,...,N, is a corresponding prolongation operator; we always use the symmetric choice of P_i := R_i^T throughout this article. 

For the definition of the second and coarse level, we assume to have a coarse space V_0 and a prolongation operator P_0: V_0 → V. In practice, we usually compute the operator P_0, which is a column-wise collection of coarse basis functions discretized in V and afterwards define 

    V_0 := { v_0 = P_0^T v   |   v ∈ V }.

Let us note that we are usually interested in coarse spaces which can be built without having a coarse discretization of Ω. Nonetheless, if such coarse finite elements defining V_0, the corresponding operators R_0 and P_0 can be built by simply interpolating the finite element shape functions from V_0 to the fine space V. Details on the construction of P_0 for the Navier-Stokes equations without having a coarse discretization are given in <ref>. With R_0 := P_0^T, we can define a nonlinear coarse correction T_0(u) by

    R_0 F(u- P_0 T_0(u))=0.


Now we have all ingredients to define new left-preconditioned nonlinear operators and we start with an additive variant. To obtain an additive nonlinear two-level Schwarz method, the nonlinear problem

    ℱ_a(u) := ∑_i=1^N P_i T_i(u) + P_0 T_0(u)

is solved using Newton's method. We refer to this method as two-level ASPEN (Additive Schwarz Preconditioned Exact Newton) based on the one-level approaches from <cit.>, or, even simpler, as additive nonlinear two-level Schwarz method. Similarly, a hybrid nonlinear two-level Schwarz method can be obtained by solving

    ℱ_h(u) := ∑_i=1^N P_i T_i(u-P_0 T_0(u)) + P_0 T_0(u)

with Newton's method. Here, the coupling between the levels is multiplicative while the subdomains are still coupled additively against each other. For more details and further nonlinear two-level variants, we refer to <cit.>.

Let us remark that in each Newton iteration a linear system with the Jacobian Dℱ_X(u^(k)) has to be solved iteratively with, for example, GMRES (Generalized Minimal Residual), where X ∈{a,h }. Fortunately, Dℱ_X(u^(k)) has already the favorable structure of a linear two-level Schwarz preconditioner, either additive or hybrid, times the Jacobian DF(·). Therefore, no additional linear preconditioner is necessary to solve the linearized systems efficiently. For details on the exact shape and the different linearization points of the local and the coarse part of the preconditioner, we refer to <cit.>. For experts, let us remark that we use always the exakt Jacobian and not the approximation used, for example, in the well-known ASPIN approach <cit.>. All in all, nonlinear two-level Schwarz algorithms have a special structure. First, there is an outer Newton loop where in each iteration a linearized system with the matrix D ℱ_X(u^(k)) is solved iteratively using GMRES. Second, to compute the residual ℱ_X(u^(k)) and the Jacobian D ℱ_X(u^(k)), all corrections T_i(u^(k)),  i=0,...,N, have to be computed in each outer Newton iteration, which is done by solving <ref> and <ref> with Newton's method. This immediately leads to inner Newton iterations which are denoted by local inner loops in case of the local corrections and the coarse inner loop in case of the coarse correction. Let us note that all local inner loops can be easily run in parallel but the coarse inner loop can only be carried out in parallel in the case of an additive coupling between levels. In the hybrid variant, the coarse loop is carried out before the subdomain corrections are computed. We provide a brief algorithmic overview of two-level nonlinear Schwarz methods in <ref>. 




To get a good estimate of the performance of our nonlinear Schwarz methods, we compare against classical Newton-Krylov-Schwarz, where Newton's method is applied to solve F(u)=0 directly and all linear systems are solved using GMRES and a linear two-level Schwarz preconditioner, either hybrid or additive. 
In our comparison, we always use the same coarse spaces for Newton-Krylov-Schwarz and nonlinear Schwarz. 




§ COARSE SPACES




For our Navier-Stokes problems we use the monolithic basis functions introduced in <cit.>, which are based on the GDSW coarse spaces; see <cit.> for scalar elliptic problems. 

We only provide a brief description here. Let us first assume that we have a symmetric and positive definite matrix K which is obtained from the discretization of a scalar linear partial differential equation, for example, a linear diffusion equation with sufficient boundary constraints of Dirichlet type. Thus, we have a single degree of freedom in each node of the finite element mesh and, for the moment, do not have to distinguish between nodes and degrees of freedom. We partition K into interface and interior nodes

    K = [  K_II K_I Γ; K_Γ I  K_ΓΓ ],

where the interface Γ contains all nodes on the boundaries of the nonoverlapping subdomains Ω_i',  i=1,...,N, except the nodes belonging to a part of ∂Ω where a Dirichlet boundary condition is imposed. The set I contains all remaining nodes. Furthermore, we partition the interface Γ into n_c overlapping or nonoverlapping patches Γ_j,  j=1,...,n_c, and define associated functions or, more precisely, vectors Φ_Γ^j,   j=1,...,n_c, which are defined on Γ but are only nonzero on the corresponding patch Γ_j,  j=1,...,n_c. Additionally, we construct Φ_Γ^j such that we obtain a partition of unity ∑_j=1^n_cΦ_Γ^j = 1 on the complete interface. 
Using the discrete energy minimizing extension

    Φ_I^j := - K_II^-1 K_I ΓΦ_Γ^j

we define the coarse basis function belonging to the patch Γ_j,   j=1,...,n_c, by 
    Φ^j := [ Φ_I^j; Φ_Γ^j ]   and   by  Φ := [ Φ^1,...,Φ^n_c ]

the matrix which contains all coarse basis functions. Let us remark that P_0 := Φ gives us the prolongation from the coarse space V_0 to the fine space V and that we have a partition of unity

    ∑_j=1^n_cΦ^j = 1

on Ω with 1 being a vector of ones, since K is symmetric positive definite. 

Using this concept, different partitions into patches Γ_j and different definitions of Φ_Γ^j are possible. Partitioning Γ into n_e edges and n_v vertices in two dimensions and defining Φ_Γ^j,   j=1,...,n_c, with n_c=n_v+n_e to be one on the corresponding edge or vertex leads to the classical GDSW coarse space. Alternatively, we can define one patch Γ_j,   j=1,...,n_v, for each vertex, consisting of the vertex itself and all adjacent edges. Then, Φ_Γ^j can be defined to be one in the vertex and to linearly decrease to zero along the adjacent edges. This gives a smaller coarse space of size n_v but also fulfills all necessary properties. This coarse space is called MsFEM-D coarse space in <cit.> and can also be interpreted as one instance of the class of reduced GDSW (RGDSW) coarse spaces; see <cit.> for more details on RGDSW. In the present article, we denote this coarse space as RGDSW coarse space of type A. An alternative approach to define Φ_Γ^j is to set it to 1 in the vertex and to 0.5 on the adjacent edges. Let us note that if an edge ends at the boundary ∂Ω we set it to 1 instead of 0.5. This coarse space can be implemented in an algebraic way without information about the geometry of the subdomains. We refer to this as RGDSW type B coarse space.

In general, for a nonlinear problem, one can use the Jacobian K := DF(u^(k)) in the k-th linearization u^(k) of Newton's method to compute the coarse basis functions. We consider two approaches: 1) compute Φ once using K:= DF(u^(0)) in the initial value u^(0) and then recycle Φ in all further steps or 2) recompute Φ with K := DF(u^(k)) at the beginning of each Newton step. Both versions can be used either for Newton-Krylov-Schwarz or nonlinear two-level Schwarz, where Φ can be recomputed at the beginning of each outer iteration. 

Considering the Navier-Stokes equations in two dimensions, we have three degrees of freedom for each node, that is, two velocities v=(v_x,v_y) and the pressure p. To build a partition of unity on the interface, we now define three basis functions Φ_Γ^j,i, i=1,2,3,  j=1,...,n_c, for each patch Γ_j,   j=1,...,n_c; one for each degree of freedom (velocities and pressure). Sorting the degrees of freedom appropriately, that is, first all x-velocities, then all y-velocities, and finally all pressures, we simply have

    Φ_Γ^j,1=[ Φ_Γ^j;     0;     0 ],   Φ_Γ^j,2=[     0; Φ_Γ^j;     0 ],    and  Φ_Γ^j,3=[     0;     0; Φ_Γ^j ]

and, of course,

    ∑_i=1^3 ∑_j=1^n_cΦ_Γ^j,i =1

on the interface. For the extension to the interior degrees of freedom, the complete block matrix

    K = DF(u^(k)) = [ A(u^(k))      B^T;        B        0 ]
 
from the linearization of the discretized Navier-Stokes equation is used, which is clearly a monolithic approach to compute the coarse basis. Here, as usual, only the velocity part A depends on the current velocities and pressures u^(k) = (v^(k),p^(k)). We note that, since K is no longer a symmetric positive definite matrix, we have a lack of theory and cannot show that the coarse basis functions Φ build a partition of unity on the interior degrees of freedom. Nonetheless, in case of Newton-Krylov-Schwarz, this approach has already shown good performance; see <cit.>. Let us remark that the resulting Φ can also be written as

    Φ = [ Φ_vv Φ_vp; Φ_pv Φ_pp ],

where the first column contains all basis functions belonging to the velocities v and the second column contains all basis functions belonging to the pressure variables. In our computations, we use a modification by deleting the coupling blocks, that is,

    Φ = [ Φ_vv    0;    0 Φ_pp ].

These modified coarse basis functions yield better results in the case of nonlinear Schwarz methods and have already been used in <cit.>. For the definition of the patches Γ_j and corresponding scalar basis functions Φ_Γ^j we use both RGDSW type approaches suggested above.

To further improve the linear convergence, we additionally enforce a zero pressure average for the subdomain corrections in the linear solves. Following <cit.>, we define the local projection 
    P_i := I_p - a_i^T (a_i   a_i^T)^-1 a_i,  i=1,...,N,
 on the pressure variables of the subdomain Ω_i,  i=1,...,N, with a_i=1 is a vector of ones; one for each local pressure variable. We then define the local projection

    𝒫_i := [ I_v   0;   0 P_i ]

for each subdomain Ω_i, i=1,...,N, which is the identity for all local velocities but enforces a zero average in the pressure.

Usually, we prefer to use restricted versions of nonlinear Schwarz. In the restricted variants, the local projections P_i,  i=1,...,N, which are used to recombine the nonlinear local corrections are replaced by P_i,  i=1,...,N, where the P_i fulfill a partition of unity property on Ω, that is, ∑_i=1^N P_i R_i = I.
 In our numerical results, due to space limitatons, we only consider the hybrid two-level restricted nonlinear Schwarz approach with or without local pressure projections, that is, we either solve

    ℱ_h(u) := ∑_i=1^N P_i 𝒫_i T_i(u-P_0 T_0(u)) + P_0 T_0(u) = 0

or 

    ℱ_h(u) := ∑_i=1^N P_i T_i(u-P_0 T_0(u)) + P_0 T_0(u) = 0

with Newton's method. Let us remark that the local projections have successfully been used to reduce the number of GMRES iterations in the case of Newton-Krylov-Schwarz but now, in the case of nonlinear Schwarz, 𝒫_i actually appears in the nonlinear formulation and might influence the convergence of Newton's method.  

Let us give some details on the initial value of Newton's method. We decided to use an initial value which fulfills all Dirichlet boundary constraints such that we can use homogeneous boundary conditions in all linear solves. We thus choose v^(0)|_∂Ω=v_0 and v^(0)=0 elsewhere. We use p^(0)=0 for the pressure. Since we have zero Dirichlet boundary conditions for the velocity now, we only choose coarse basis functions in the vertices which are not part of the boundary for both, velocity as well as for the pressure.


To summarize, we vary three different aspects in the numerical tests: 1) the coarse space using two different types of RGDSW, 2) using or not using local pressure projections P_i, and 3) recycling the coarse basis functions Φ or recomputing them in each Newton iteration.   




 


§ NUMERICAL RESULTS FOR THE LID-DRIVEN CAVITY FLOW PROBLEM

 
As already stated, we consider the classical lid-driven cavity problem; see <ref> for details. As a classical approach, we use Newton-Krylov-Schwarz with a linear hybrid restricted two-level Schwarz preconditioner. We always use the monolithic RGDSW basis functions described above to build the second level using either type A or B to define them on the interface. The linear preconditioner is not varied since it will not affect the nonlinear convergence in case of Newton-Krylov-Schwarz. We choose only the best possible option which is, in this specific case,  recycling the coarse basis functions Φ and using the local pressure projections 𝒫_i,  i=1,...,N. As a nonlinear Schwarz approach, we use the hybrid restricted two-level Schwarz approach with or without local pressure projection, that is, we  consider one of the residual formulations <ref> or <ref> and solve with Newton's method. We can optionally also use recycling of the coarse basis functions in nonlinear Schwarz methods and also compare both initial values.

We stop the nonlinear iteration when the relative residual ||F(u^(k))||_2/||F(u^(0))||_2 is smaller than 10^-6. All inner local and coarse nonlinear iterations in nonlinear Schwarz are stopped early when a reduction of the relative residual of at least three orders of magnitude is reached. Furthermore, all iterative solves with GMRES are stopped when the relative residual is below 10^-8. We test for different Reynolds numbers Re up to 2500.

To increase the robustness of Newton's method, we use a simple backtracking approach for inexact Newton methods suggested in <cit.>. We use this approach for the local and coarse inner loops in nonlinear Schwarz as well as for all classical Newton-Krylov-Schwarz approaches. The step size cannot get smaller then 0.01 in the line search, which gave the best results for out specific problem. We always use 256 subdomains, since using less subdomains leads to very small coarse spaces which have no positive effect in the case of nonlinear two-level Schwarz methods. Using more subdomains is too time consuming using our MATLAB implementation. We use 128 finite elements (Taylor-Hood elements) to discretize each nonoverlapping subdomain and an overlap of δ=3h. As a consequence, the original global problem has 148 739 degrees of freedom in total. In <ref> and <ref> we present the obtained results for the different combinations.
 







In general, we can observe that nonlinear Schwarz methods can improve the nonlinear convergence, especially for high Reynolds numbers where the classical Newton approach does not converge within 20 iterations. For small Reynolds numbers the global Newton method is, as expected, the best choice, since it provides comparably fast convergence without the additional work introduced by the inner loops. The best setup in our tests was to always recycle the coarse basis functions and to use the local pressure projection 𝒫_i. The latter ones improve the linear convergence without deteriorating the fast convergence of the nonlinear solver. Comparing both RGDSW coarse spaces, the nonlinear convergence is very similar but slightly better for type B coarse spaces up to a Reynolds number of 2000. Additionally, the inner coarse loop converges a bit faster. On the other hand, the type A coarse space leads to a better linear convergence. For the case with Re=2500, nonlinear Schwarz with type A coarse basis functions divergence caused by a divergent inner course iteration. To summarize, type B RGDSW coarse basis functions proved to be more robust in our experiments. Additionally, the type B coarse spaces can be built without geometry information on the subdomains while for the type A coarse basis functions the coordinates of the nodes have to be known. For a more detailed picture of the convergence behavior, we provide <ref> including the most promising variants for Reynolds numbers of 1000 and 1500. 

Furthermore, we provide a visualization of the coarse corrections T_0 for a Reynolds number of 500 and a decomposition into 64 subdomains in <ref>. Although the convergence speed is nearly identical for both coarse spaces, it can be seen that, in general, the convergence strongly depends on the choice of the coarse basis and different coarse corrections are computed. 

Finally, we show the Newton iterates for the case with Re=1 500 for Newton-Krylov-Schwarz and nonlinear Schwarz to compare the convergence behavior; see <ref>. Already after one iteration, nonlinear Schwarz visibly converges against the solution while Newton-Krylov-Schwarz tends to form more vortices and clearly diverges, although a backtracking approach is used to control the step lengths.







The authors acknowledge the financial support by the German Federal Ministry of Education and Research (BMBF) for the project Stroemungsraum within the exascale computing programme SCALEXA.	


spmpsci     



