Preconditioners for Saddle Point Problems on Truncated Domains in Phase Separation Modelling

09/28/2017 ∙ by Pawan Kumar, et al. ∙ ProtonMail 0

The discretization of Cahn-Hilliard equation with obstacle potential leads to a block 2 by 2 non-linear system, where the p1, 1q block has a non-linear and non-smooth term. Recently a globally convergent Newton Schur method was proposed for the non-linear Schur complement corresponding to this non-linear system. The solver may be seen as an inexact Uzawa method which has the falvour of an active set method in the sense that the active sets are first identified by solving a quadratic obstacle problem corresponding to the p1, 1q block of the block 2 by 2 nonlinear system, and a new decent direction is obtained after discarding the active set region. The problem becomes linear on nonactive set, and corresponds to solving a linear saddle point problem on truncated domains. For solving the quadratic obstacle problem, various optimal multigrid like methods have been proposed. In this paper solvers for the truncated saddle point problem is considered. Three preconditioners are considered, two of them have block diagonal structure, and the third one has block tridiagonal structure. One of the block diagonal preconditioners is obtained by adding certain scaling of stiffness and mass matrices, whereas, the remaining two involves Schur complement. Eigenvalue bound and condition number estimates are derived for the preconditioned untruncated problem. It is shown that the extreme eigenvalues of the preconditioned truncated system remain bounded by the extreme eigenvalues of the preconditioned untruncated system. Numerical experiments confirm the optimality of the solvers.



There are no comments yet.


page 23

page 24

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

The Cahn-Hilliard equation was first proposed in 1958 by Cahn and Hilliard [6] to study the phase separation process in a binary alloy. Here the term phase stands for the concentration of different components in an alloy. It has been empirically observed that the concentration changes from a given mixed state to a spatially separated two phase state when the alloy under preparation is subjected to a rapid cooling below a critical temperature. This rapid reduction in the temperature the so-called deep quench limit has been found to be modeled efficiently by obstacle potential proposed by Oono and Puri [19, Fig. 7, p. 439] in 1987, and analyzed by Blowey and Elliot [3, p. 237, (1.14)]. The phase separation has been noted to be highly non-linear, and the obstacle potential emulates the nonlinearity and non-smoothness that is empirically observed and much desired in numerical smulations. Consequently, handling the non-smoothness as well as designing robust iterative procedure have been a subject of much active research during last decades. Assuming semi-implicit time discretizations [4] to alleviate the time step restrictions, most of the proposed methods essentially differ in the way the nonlinearity and non-smoothness are handled. There seems to be two main approaches to handle the non-smoothness: regularization around the non-smooth region and subsequently using a variant of smooth solvers, for example, as in [5], or an active set like approach [11], i.e., where one identifies the active sets via a nonlinear solver, subsequently, after discarding the active set nodes, we obtain a reduced (or truncated) problem which is linear. Moreover, the global convergence of the nonlinear solver may be ensured by a proper damping parameter, for example, as done in [11].

The non-linear problem corresponding to Cahn-Hilliard equation with obstacle potential could be written as a non-linear system in block matrix form as follows:


where and are unknowns corresponding to order parameter and chemical potential respectively, where denotes the indicator functional for corresponding to the admissible set We note that is a set valued mapping due to the presence of set-valued operator hence, we have inclusion in (1) instead of equality. The matrix corresponds to Laplacian with Neumann boundary conditions perturbed by a rank one term, and is multiplied by a parameter corresponding to interface width. On the other hand, is also Laplacian with Neumann boundary condition, but multiplied by the time step parameter. Both nonlinearity and non-smoothness are due to in Various non-linear and nonsmooth solvers have been proposed for (1) [1, 5].

By nonlinear Gaussian elimination of the system above could be reduced to a nonlinear Schur complement system in variables [11], where the “negative” nonlinear Schur complement is given by Here is understood as inversion in the nonlinear sense. In [11], a globally convergent Newton method is proposed for this nonlinear Schur complement system, which is interpreted as a preconditioned Uzawa iteration. To solve the inclusion corresponding to the quadratic obstacle problem, many methods have been proposed such as block Gauss-Seidel [2, 7], monotone multigrid method [14, 15, 17], truncated monotone multigrid [12], and truncated Newton multigrid [12].

Once active sets are identified, the corresponding rows and columns are anhilated, we then obtain a reduced linear system as follows


Here a solution to (2) is a new descent direction in the Uzawa iteration. By a choice of appropriate step size along this descent direction, global convergence of the Uzawa method is ensured. As the active sets change during each iteration, the linear system, and hence the preconditioners need to be updated.

In this paper, our goal is to design effective preconditioner and hence an iterative solver for (2) such that the convergence rate is independent of problem parameters and mesh size. Three preconditioners are considered; two of them involves Schur complement. Two of these preconditioners have block diagonal structure and they correspond to non-standard norms proposed in [25]. To approximate the Schur complement, we consider an approximation proposed in [5]. It turns out that the building blocks of these preconditioners are same, their analysis is remarkably similar, even though, they may look structurally different from the outset. Eigenvalue bound and condition number estimates are derived for these preconditioners for the untruncated problem. The obtained eigenvalue bounds seem to be tight when compared to numerically computed extreme eigenvalues. Subsequently, it is shown that the extreme eigenvalues of the preconditioned truncated problem are bounded from above and below by the extreme eigenvalues of the corresponding preconditioned untruncated problem. We also verify the effectiveness of these preconditioners numerically for various evolutions.

The rest of this paper is organized as follows. In Section 3, we describe the Cahn-Hilliard model with obstacle potential, we discuss the time and space discretizations and variational formulations. In Section 4, we discuss briefly the solver for Cahn-Hilliard with obstacle problem. The preconditioners for the truncated linear saddle point problem (2), and their eigenvalue analysis are discussed in Section 5. Finally, in Section 6, we show numerical experiments with the proposed preconditioners.

2 Notations

Let SPD and SPSD denote symmetric positive definite and symmetric positive semi definite respectively. Let denote the condition number of SPD matrix For denotes the absolute value of whereas, for any set denotes the number of elements in Let

denote the identity matrix. Let

denote For a matrix with all real eigenvalues, the eigenvalues will be denoted and ordered as follows


3 Cahn-Hilliard Problem with Obstacle Potential

3.1 The Model

We will consider a model for phase separation of two components in a binary alloy mixture. Here phase stands for concentration of two components in the mixture. Let be the concentration of two components in the mixture, then we set The phase separation is modelled using Cahn-Hilliard equations, which is obtained by gradient flow of Ginzburg-Landau (GL) energy functional which is given as follows


Here the constant relates to interfacial thickness, and the obstacle potential which is used to model deep quench phenomena is given as follows


Here the subscript of indicator function above denotes the range of admissible values of Here is defined as follows


Moreover, is assumed to be conserved. The gradient flow of leads to the Cahn-Hilliard equation in PDE form


The unknowns and are called order parameter and chemical potential respectively. For a given final time and initial condition where


the equivalent initial value problem for Cahn-Hilliard equation with obstacle potential interpreted as variational inequality reads


where we use the notation to denote the duality pairing of and Note that we used the fact that in the second term on the left of inequality (13) above. The inequalities (12) and (13) are defined on constrained set the variational inequality of first kind is also equivalently represented on unconstrained set using the indicator functional [7, p. 2]. The existence and uniqueness of the solution of (12) and (13) above have been established in Blowey and Elliot [3]. We next consider an appropriate discretization in time and space for (12) and (13).

3.2 Time and space discretizations

We consider a fixed non-adaptive grid in time interval and in space defined in (4). The time step is kept uniform. We consider the semi-implicit Euler discretization in time and finite element discretization in space as in Barrett et. al. [2] with triangulation with the following spaces


which leads to the following discrete Cahn-Hilliard problem with obstacle potential:

Find s.t.


holds for each The initial solution is taken to be the discrete projection

Existence and uniqueness of the discrete Cahn-Hilliard equations has been established in [4]. The discrete Cahn-Hilliard equation is equivalent to the set valued saddle point block nonlinear system (1) with and


We write the above in more compact notations as follows


where and are usual notations for mass and stiffness matrices respectively.

4 Iterative solver for Cahn-Hilliard with obstacle potential

In [11], a nonsmooth Newton Schur method is proposed which is also interpreted as a preconditioned Uzawa iteration. For a given time step the Uzawa iteration reads:


for the saddle point problem (1). Here denotes the Uzawa step, and denotes the time step. Here and are defined as follows


The time loop starts with an initial value for which can be taken arbitrary as the method is globally convergent, and with the initial value obtained from (23). The Uzawa iteration requires three main computations that we describe below.

4.1 Computing

The first step (22) corresponds to solving a quadratic obstacle problem interpreted as a minimization problem as follows


As mentioned in the introduction, this problem has been extensively studied during last decades [2, 12, 14, 15].

4.1.1 Algebraic Monotone Multigrid for Obstacle Problem

To solve the quadratic obstacle problem (22), we use the monotone multigrid method proposed in [14]. In Algorithm 1, we describe an algebraic variant of the method. The algorithm performs one V-cycle of multigrid; it takes from the previous iteration, and outputs the improved solution

The initial set of interpolation operators are constructed using aggregation based coarsening


4.2 Computing

The quantity in (23) is obtained as a solution of the following reduced linear block system:




Here truncation matrices and are defined as follows:


where is the th component of and is the th diagonal entry of In words, is the matrix obtained from by replacing the th row and

th column by the unit vector

corresponding to the active sets identified by diagonal entries of Similarly, is the matrix obtained from by annihilating rows, and is the matrix obtained from by annihilating columns. Rewriting untruncated version of (26) in simpler notation as follows


where By a change of variable we obtain


where Furthermore, we modify the term of the system matrix above as follows


where Now the untruncated system may be rewritten as


where is a rank one term with proper extension by zero. Now we are in a position to use Sherman-Woodbury inversion for matrix plus rank-one term to invert . In this paper, we shall develop efficient solvers to solve the truncated system


where is defined next. We denote


Note that the notation appearing above has now been redefined. Thus, the truncated system corresponding to (32) reads


Thus, Sherman-Woodbury inversion formula may be used to invert and it is enough to find an efficient solver for (33).

4.3 Computing step length

The step length can be computed using a bisection method. We refer the interested reader to [10][p. 88].

0:  Let and let
0:   solution from previous cycle or a given initial solution
1:  Compute residual:
2:  Compute defect obstacles:
3:  for  do
4:     Projected Gauss-Seidel Solve using Algorithm 2:
5:     Update
6:     Restrict and compute new obstacle
7:  end for
8:  Solve
9:  for  do
10:     Add corrections
11:  end for
12:  Compute
12:  improved solution
Algorithm 1 Monotone Multigrid (MMG) V cycle
0:    current iterate
0:  new iterate
1:  Compute residual:
2:  Compute defect obstacles:
3:  for  do
4:     for  do
5:        Compute
6:     end for
7:  end for
Algorithm 2 PGS()

4.4 Mixed Finite Element Formulation of Reduced Linear System

We choose suitable Hilbert spaces for trial and test spaces as follows


where is the domain where truncation takes place. Indeed, if is empty, then and we set

The weak form of the partial differential equations corresponding to the truncated system (

33) reads




The mixed variational problem above can also be written as a variational form on product spaces


where and are defined as follows


for and The corresponding bilinear form for the untruncated system is given as follows


for and where The mixed variational problem corresponding to untruncated system now reads


In the rest of this paper, we shall consider norms proposed in [25] as follows


where and are inner products of Hilbert spaces and respectively. As will see shortly such norms lead to block diagonal preconditioners. The boundedness condition that we seek for the mixed problem for the untruncated problem reads


We have the following conjecture for the truncated problem


Similarly, for well-posedness of (53), following well known Babuska-Brezzi condition needs to be satisfied


Similarly, it is not evident whether the following inequality must hold.


We shall call the norms and optimal if the constants and remain independent of the problem parameters: and , moreover, in the discrete space, also remains independent of the mesh size The reason why we are interested in the inequalities (56) and (58) is that any optimal norm that is found for untruncated problem shall lead to optimal norm for truncated problem as well. Note that boundary of untruncated problem has certian regularity (for example Lipschitz continuity), but for the truncated problem no such regularity is to be assumed, because the truncations are assumed to be arbitrary. Our plan of attack is to use the approach of [25], which is readily applicable for our untruncated problem. Although, (56) and (58

) are left as conjecture for the moment, we shall try to answer this in the discrete case: we shall show a related result that the extreme eigenvalues of the truncated preconditioned operator are bounded by the extreme eigenvalues of the corresponding untruncated preconditioned operator. Hence, in the following, we first derive optimal preconditioners for the untruncated problem.

We shall provide equivalent conditions as in [25] for (55) and (57) that lead to deriving the optimal norms, i.e., optimal preconditioners. But first we introduce some notations for operators corresponding to bilinear forms. It is easy to see that is a Hilbert space itself as and are themselves Hilbert spaces. It is convenient to associate linear operators for the bilinear forms and as follows


Consequently, the operator corresponding to mixed bilinear form and the right hand side (reusing the notation) in operator notation are given as follows


The untruncated problem is denoted as follows


and the corresponding truncated problem reads


where is given as follows


where analogous to (59), we have following definitions for truncated operators


where is defined as in (59). In [25], starting from the abstract theory on Hilbert spaces that lead to representation of isometries, a preconditioner is proposed; it is based on non-standard norms, or isometries that correspond to block diagonal preconditioner of the following form


In the next section, our goal is to determine and

4.5 Choice of norm: a brief introduction to Zulehner’s idea

Before we move further, we introduce some notations. The duality pairing on is defined as follows

Let be an isometric isomorphism defined as follows

The inverse is Riesz-isomorphism, by which functionals in can be identified with elements in and we have

We already chose the type of norm in (54), we now look for explicit representation of isometries or norms in terms of operators defined in (59). The main ingredient is the following theorem.

Theorem 4.1.

[25][p. 543, Th. 2.6] If there are constants such that






is satisfied with constants that depend only on and on And, vice versa, if the estimates (68) are satisfied with constants then the estimates (66) and (67) are satisfied.

Equivalently, as conjectured for (56) and (58), and recalling that we may ask whether the following bounds hold for truncated system


However, we shall show a similar result in finite dimension using Fischer’s theorem in Lemma 5.2. In [25], the terms and in (66) and (67) respectively are defined using isometries and as follows:


Using (70), (66) and (67) are equivalently written as follows

In short, in new notation meaning “spectrally similar”, we obtain the following equivalent conditions for isometries

Let and be any SPD matrices, consequently, they define inner products and a Hilbert space structure in The intermediate Hilbert spaces between and are given as follows

Continuing from above, when and are non singular, the generic form of the norms are given by the following lemma.

Lemma 4.1.

Let , consequently, be nonsingular. Then


See [25][p. 547-548]. ∎

The isometries and above provide a general template for obtaining a variety of preconditioners. Obviously, our goal is to find those that are easier to compute with numerically. Before we propose preconditioners, we shall need some properties of the block of and that for the negative Schur complement Such properties will be useful in developing preconditioners using and

4.6 Properties of the system matrix and Schur complement

An important property that we shall need shortly when analyzing preconditioners is that the eigenvalues of the truncated matrix is bounded from above and below by the eigenvalues of the untruncated matrix.

Lemma 4.2.

The operator is symmetric and indefinite.


Symmetry is obvious. Indefiniteness follows from below:


For the choice of

Lemma 4.3.

is SPD.


From (34), we recall that Here being a stiffness matrix corresponding to natural boundary condition is SPD except on the span of vector which is in the kernel of but Also, being a mass matrix is SPD, is SPD. ∎

Fact 4.1 (Permutation preserves eigenvalues).

Let be a permutation matrix, then and are similar.


being a permutation matrix, hence the proof. ∎

Lemma 4.4 (Poincare separation theorem for eigenvalues).

Let be any symmetric matrix with eigenvalues and let be a semi-orthogonal matrix such that Then the eigenvalues of are separated by the eigenvalues of as follows