{'fraction_non_alphanumeric': 0.040795559666975026, 'fraction_numerical': 0.016975023126734504, 'mean_word_length': 4.052816078523019, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Preconditioning nonlocal multi-phase flow

We propose an efficient solver for saddle point problems arising from finite element approximations of nonlocal multi-phase Allen--Cahn variational inequalities. The solver is seen to behave mesh independently and to have only a very mild dependence on the number of phase field variables. In addition we prove convergence, in three GMRES iterations, of the approximation of the two phase problem, regardless of mesh size or interfacial width. Numerical results are presented that illustrate the competitiveness of this approach.

1. Introduction. The aim of this paper is to combine preconditioning methods for indefinite problems and multigrid preconditioning developed for elliptic systems to provide an efficient preconditioner for the solution of systems of multiphase Allen-Cahn variational inequalities of the form: P Q : For given u(·, 0) = u 0 ∈ G Q , find u ∈ L 2 (0, T ; G Q ) ∩ H 1 (0, T ; L 2 (Ω)) such that Here Ω ∈ R d , d = 1, 2 or 3 and u : Ω × (0, T ) → R N denotes the vector-valued phase field function which describes the fractions of the N phases, i.e. each component of u describes one phase, A ∈ R N ×N is a symmetric constant matrix that has at least one positive eigenvalue, We denote by L 2 (Ω) and H 1 (Ω) the spaces of vector-valued functions, (., .) is the standard L 2 inner product for scalar functions, (v, w) = |Ω| Ω f (x) dx. The system P Q arises from steepest descent dynamics with respect to the L 2 -norm of the Ginzburg-Landau energy, under the constraint − Ω u dx = Q. Here Ψ is the multi-obstacle potential Ψ(θ) := − 1 2 θ T Aθ + I G = − 1 2 θ T Aθ for θ ∈ G ∞ otherwise with I G denoting the indicator function for the Gibbs Simplex, G := {ξ ∈ R N | ξ ≥ 0, ξ ·1 = 1}, and the symmetric constant matrix A has at least one positive eigenvalue to allow for minima of Ψ to exist, see [13].
Remark 1.1. Steepest descent dynamics with respect to the L 2 -norm of the Ginzburg-Landau energy E(u), without the constraint − Ω u dx = Q, yields the system: P: For given u(·, 0) = u 0 ∈ G, find u ∈ L 2 (0, T ; G) ∩ H 1 (0, T ; L 2 (Ω)) such that where G := {v ∈ H 1 (Ω) | v ∈ G a.e.}. Since P is a simplified version of P Q the solver we propose in this paper can be applied to the corresponding finite element approximation of P.
P is a generalisation of the scalar Allen-Cahn equation that was introduced by Allen and Cahn [1] to describe the capillarity driven evolution of an interface separating two bulk phases. The parameter ε, with 0 < ε 1, is associated with the thickness of the diffuse interfacial layer in which the phase field variables rapidly change their value. The N phase extension of the scalar Allen-Cahn model was introduced in [7,16]. The nonlocal problem P Q models interface evolution with mass conservation.
Multiphase Allen-Cahn models have a variety of applications, including image segmentation, see for example [23], and identification of coefficients in elliptic equations [10]. Applications arising from identification of coefficients in elliptic equations include electric impedance tomography, [12,21,20], and flow in porous media with unknown permeabilities [9,11,6]. Applications of mass conserving multiphase Allen-Cahn models include structural topology optimisation [30,3].
Efficient and reliable, i.e., fast and globally converging, multigrid methods for solving implicit in time finite element approximations of P are presented in [24,25], while an explicit in time finite element approximation of P Q was introduced in [15]. In [4] (semi-)implicit in time finite element approximations of P Q are considered in which a primal-dual active set method, see [2,19], is used to solve the finite element approximations. By using Krylov-subspace solvers and suitable preconditioners the authors in [5] develop efficient, mesh independent, solvers for the (semi-)implicit approximations of P and P Q that were derived in [4]. In this work we introduce an alternative preconditioner to the ones in [5] resulting in a solver that is not only mesh independent, but also is only mildly dependent of the number of phases N .
We note that in [17] globally convergent nonsmooth Schur-Newton methods are introduced for the solution of discrete multicomponent Cahn-Hilliard systems with logarithmic and obstacle potentials. These methods could also be used to solve the multicomponent Allen-Cahn systems P and P Q .
When using iterative techniques to solve the linear system that arises when the primal-dual active set method is used to solve a finite element approximation of P Q , the majority of the work to be undertaken within each iteration is in the solving of the linear systems Kx = b where K is symmetric positive definite. Similar saddle point structures are common place within fluid dynamics, leading to much development of numerical solvers for Navier-Stokes equations, e.g. see [26,22,28]. In these papers, it is the choice of preconditioning matrix, P , that leads to improved convergence of the chosen iterative Krylov subspace scheme, e.g. see [27]. The preconditioning of the linear systems arising from P and P Q was initially considered in [5], where a preconditioning technique building upon Stokes type systems is proposed. In this work we will use a similar structural approach to that of [5], however, we provide an improved approximation to the Schur complement, BK −1 B T , leading to a solver that is almost independent of the number of phases. Moreover, this improvement does not lead to any significant increase in computational effort per iteration, ultimately leading to a more effective solver. In addition, when considering the two phase problem, with N = 2, the minimal polynomial of the resulting preconditioned system is of degree three and hence GMRES will converge within three iterations, see Theorem 3.1.
The paper is organised as follows. In Section 2 we reformulate P Q with the help of Lagrange multipliers, yielding the associated system R Q . We then introduce a finite element approximation of an implicit Euler-discretisation of R Q and we apply a primal-dual active set algorithm to this discretisation. In Section 3 a preconditioner for the primal-dual active set algorithm is developed and the implementation of the numerical solver is presented. In Section 4 we present numerical computations that illustrate the effectiveness of our approach, in particular they show how the iteration number is independent of the mesh size and only mildly dependent on the number of phases N .

Alternative Formulation and Finite Element Discretisation.
In this section we follow the authors in [4] in applying a primal-dual active set method, [2,19], to a finite element approximation of P Q , this method is well known in the context of optimisation with partial differential equations as constraints. To this end we first reformulate P Q with the help of Lagrange multipliers, yielding the associated system R Q , then we apply a Primal Dual Active Set algorithm to a finite element approximation of R Q .
R Q is complemented with the the boundary condition ∂u ∂ν = 0, were ν is the outer unit normal to ∂Ω.
Remark 2.1. The scaling 1 ε µ in R Q is introduced in order that µ is of order one, if we were to replace 1 ε µ by µ we would observe a severe ε-dependence of µ which in practice often results in oscillations in the discretised primal-dual active set method.

Finite Element Discretisation.
For simplicity we assume that Ω is a polyhedral domain. Let T h be a regular triangulation of Ω into disjoint open simplices, i.e. in particular Ω = ∪ T ∈T h T . Furthermore, we define h := max T ∈T h {diam T } the maximal element size of T h and we set J to be the set of nodes of T h and {p j } j∈J to be the coordinates of these nodes. Associated with T h is the piecewise linear finite element space where we denote by P 1 (T ) the set of all affine linear functions on T . Furthermore we denote the standard nodal basis functions of S h by {χ j } j∈J and we set For simplicity of presentation we denote by u h ∈ S h the discrete solution at time t n , while the solution at the previous time step will be denoted by u n−1 h ∈ S h , and similarly for µ h , λ h and Λ h . We consider the following fully discrete approximation of R Q : Remark 2.2. Due to the Au term in P Q , the problem is non-convex, in the above discretisation we have chosen to treat this term fully explicitly. Alternative choices would be to treat this term fully implicitly, or in a semi-implicit manner, neither of which would affect the performance of our proposed solver.
2.3. The Primal Dual Active Set Method. We use the nonlinear primal dual active set (PDAS) algorithm derived in [4] to solve R h,τ Q . The algorithm is obtained by reformulating the complementarity conditions (2.5) using active sets based on the primal variable u and the dual variable µ. Here we use the notation (u k i ) j and (u n−1 i ) j where k denotes the k-th iteration in the PDAS algorithm and n − 1 is the (n − 1)-st time step. This is of course a misuse of notation for k = n − 1. In addition for c > 0, we set and we define the lumped mass diagonal matrix M := (m ij ) i,j∈J with m ij = (χ i , χ j ) h and the stiffness matrix L := (l ij ) i,j∈J with l ij = (∇χ j , ∇χ i ). We define the mass lumped vector m := (M j ) j∈J , and the entries of the matrix A by a ij , .., N } stop, otherwise set k = k + 1 and goto 1.

Preconditioning.
3.1. Schur Complement. The main computational cost in the above algorithm is the solving of the system of equations (2.6), (2.8) and (2.7). To do this we firstly introduce some matrix notation. At the k-th iteration step we define K k , to be N diagonal blocks with the i-th block, K k i , being associated with the active set, D k i , of the i-th phase equation of (2.6). More precisely, At the k-th iteration, we may write this linear system in the form; find For convenience, from here on we will drop the superscript k.
To develop a preconditioner for K, we write it in the factored form Hence, if we choose the preconditioner the eigenvalues of the preconditioned system have value one, and it can be shown that only two GMRES iterations would be needed in this case, see [26].
Since each block of K consists of the standard finite element matrix for a reactiondiffusion type equation, there exists numerous practical preconditioners and solvers for this block. In particular standard algebraic, or geometric, multigrid can effectively be applied, see [18,29]. Hence, we are left to find a fully practical approximation to the action of the block matrix S exact := S 11 S 12 S 21 S 22 .

Approximate Schur Preconditioners.
In [5] the authors consider a block upper triangular preconditioner of the form: where the diagonal blocksS ii , i = 1, 2 are given bỹ where I is the identity matrix. This choice was shown to lead to the preconditioned system where E ii , i = 1, 2, are such that all eigenvalues are close to zero. Note that with this choice of preconditioner,S ij = 0, i = j, and hence in [5] the effect of the off diagonal blocks is not considered. In this work we propose an alternative approximation by building on the ideas presented in [14]. We firstly define, T , leading to S exact = BK −1 B T . Following [14] we approximate this with Note that, where α i = m T i m i , i = 1, 2, . . . , N and Remark 3.1. How well this preconditioner performs is closely related to how well B commutes with K. We note that the matrix D 22 is diagonal. Moreover, it is a relatively large proportion of the matrix B, and will always be close to a constant diagonal matrix in the regions D n,k i . This will even be the case when an adaptive mesh refinement strategy is used to accurately capture the interfacial region, since in this active region the elements are of similar size and shape. Hence, the commutator of [D 22 , C], with any square matrix C is close to zero. This is a major factor in the quality of the approximation used in (3.3).
We consider two preconditioners developed from the above methodology. Namely, the block diagonal choice and the full approximation 3.3. Practical Preconditioning. The implementation of any of the preconditioners P i , i = 1, 2, 3, requires a practical and scalable method to calculate the action of their inverses. All three preconditioners have inverses that may be written in the form, The main work in calculating the action of these inverses is in calculating the action of the inverses of K and S i , i = 1, 2, 3. As mentioned earlier, for the matrix K numerous efficient iterative solvers exist. In the following numerical results section we chose to apply three Algebraic Multigrid (AMG) V-cycles with simple Gauss-Seidel smoothing, see [29]. For S 1 and S 2 we are only required to invert a small dense matrixS 11 and in finding the action of the inverseS 22 we only require the inversion of the diagonal lumped mass matrix. Finally, for S 3 we write, We note from earlier remarks that the construction of the matrix D 11 −D 12 D −1 22 D 21 ∈ R N −1×N −1 and its inverse is inexpensive.

The Two Phase Problem.
Theorem 3.1. When considering a two phase problem, the right preconditioners P 2 and P 3 are identical. Moreover, the resulting preconditioned system, KP −1 i , i = 2, 3, is a lower triangular matrix with the diagonal consisting of 1's and a solitary Furthermore, the minimal polynomial of the resulting system is of degree three and hence GMRES will converge within three iterations.
Proof. In the two phase case the system to be solved is of the form This leads to the exact Schur complement, Hence, where a is a scalar. This leads to the full preconditioned system Clearly, this system has only two distinct eigenvalues, σ(KP −1 2 ) = {1, a} and the minimum polynomial for this preconditioned system is of order three. Hence, when using GMRES we would expect to obtain the exact solution in no more than three iterations, see [8].
4. Numerical Results. In this section we present numerical results that show the efficiency of our proposed preconditioner P 3 , (3.5). We begin by using exact solves for each matrix in the preconditioning system, then in later results we apply three Algebraic Multigrid (AMG) V-cycles with simple Gauss-Seidel smoothing, see [29], for calculating the action of the inverse of K. We denote the approximate preconditioners, in which we use the inexact AMG solver, by P i,AM G , i = 1, 2, 3.
We note that the fully explicit discretisation of A in R h,τ Q leads to an unconditionally well posed discrete problem, allowing the use of large time steps when slow dynamics are encountered, see [24]. Regarding time stepping, throughout we will use a simple adaptive time stepping strategy whereby: 1. If the number of PDAS steps required to obtain u n+1 are fewer than 5, we set τ n+1 = 1.1τ n . 2. If the number of PDAS steps required to obtain u n+1 is between 5 and 10, τ n remains unchanged.
3. If the number of PDAS steps required to obtain u n+1 exceeds 10, we recalculate u n+1 with a time step reduced by a half, τ n = 0.5τ n . At time step n with initial time step τ n and previous solution u n using the PDAS iteration scheme calculate u n+1 . This solution is taken when the active set size does not change between iterations. We found that this led to a practical stable method. In two space dimensions we set Ω = [0, 1] 2 and in three space dimensions we set Ω = [−0.5, 0.5] 3 .
Remark 4.1. During initial calculations, it was seen that the proposed PDAS scheme was only reliable when a high tolerance was enforced on linear solve. Hence, throughout the following results we will apply a tolerance on the relative GMRES error of 1e-10. Given this constraint on the PDAS scheme, it is critical that a robust and efficient solver is used.

Two Space Dimensions.
4.1.1. Grain Coarsening. The first problem we consider is that of grain coarsening in which we start with an initially well mixed mixture of N phases. The mixture rapidly separates into bulk regions of each phase, with typically each phase having multiple bulk regions. Once this fast dynamical process has taken place, the bulk regions then slowly diffuse, see Figure 4.1 where the motion of eight phases from T = 1 to T = 100 is presented. We compare the performance of the three preconditioners P 1 , (3.2), P 2 , (3.4), and P 3 , (3.5), with respect to: the number of phases, the interface width parameter ε, and the mesh size. In each case we use exact solves for each matrix in the preconditioning system. The initial mesh, Mesh 1, has mesh size h ≈ 1/32 and all other meshes are uniform refinements of this mesh, the number of nodes of each mesh is given by |J |. In Table 4.1 we display the maximum number of GMRES iteration counts together with the average number, in the form a/b, where a is the maximum number and b is the average number. We consider three meshes, three values of ε, and three values of N . From this table, for each of the three meshes and each of the three values of ε, we see the dependence of P 1 on the number of phases, N . A similar dependence can be seen for the choice P 2 , albeit a milder one. It is P 3 that outperforms the other two choices in this regard, as it shows almost no dependence on phase number for each of the three meshes and each of the three values of ε. In addition, for the two phase problem, N = 2, we observe the three iteration convergence, stated in Theorem 3.1.  In Figure 4.2 we further investigate how the number of phases, N , affects the performance of preconditioners P 1 and P 3 . We ignore P 2 since P 3 is computationally similar whilst having superior convergence rates. We display the number of GMRES iterations throughout a simulation together with the percentage of the total number of DOFs that the active sets make up. In addition we show the effect of the active time stepping by displaying the time step size throughout the simulation. We set ε = 0.01 and use Mesh 4, which has |J | = 236868. Similar to Table 4.1, we see a strong dependence for P 1 but a very weak dependence for P 3 . The results in Figure  4.3 are displayed in the same format as those in Figure 4.2, but here we investigate the effect that the mesh size has on the performance of preconditioners P 1 and P 3 . In particular we set ε = 0.02 and N = 4 and we show results for Mesh 3, for which |J | = 59217, and Mesh 4. We conclude with Figure 4.4 in which we set N = 4 and use Mesh 4, with ε = 0.02 and ε = 0.01 to see the effect that ε has on the performance of preconditioners P 1 and P 3 . In both Figures 4.3 and 4.4 we again see a strong dependence for P 1 but a very weak dependence for P 3 .

Quadruple Junction to Triple Junction.
We now turn to an initial condition of fully formed bulk regions. We set ε = 0.005 and N = 5. For the initial data we consider a square, consisting of four phases of bulk square regions, that is surrounded by a fifth phase. This unstable initial geometry rapidly evolves so that the quadruple junction is replaced by two triple junctions with 120 • angles, see Figure  4.5.
In Figure 4.6 we present the iteration counts and active set size when using P 3,AM G . Given the maximum number of iteration counts for a given time step is only 11, we conclude that the use of the inexact AM G solver on the K block has little effect on proposed solver.   Table 4.2 we present the GMRES iteration counts when using the preconditioner P 3 . As in the two dimensional case, we see that there is little dependence on any of the parameters, mesh size, ε, or number of phases, and again we observe the three iteration convergence of the two phase problem. Turning to the fully practical preconditiponer, P 3,AM G , in Table 4.3 we present CPU timings for this problem. It is not clear how to measure how these timings scale, since as the mesh is refined and more phases are added, the active set size changes considerably. However, we feel that these non-optimized CPU timings are an excellent indicator of the scalability of the proposed approach.
used has over a half a million nodes, this leads to a system size of more than five million degrees of freedom. The evolution of these phases can be seen in Figure 4.7.   Table 4.3: CPU timings using P 3,AM G in three space dimensions. Initial condition is well mixed and T = 2.
initial mesh has over half a million nodes and ε = 0.04. We see the evolution into bulk regions with spherical like minimal surfaces in contact with the fifth phase. Moreover, the central region shifts, in a similar way to the two dimensional problem, to remove any quadruple junctions. The final plot in Figure 4.9 shows the desired iteration counts for the solver.

Conclusions.
In this work we have presented a robust practical preconditioner for systems of multiphase Allen-Cahn variational inequalities. As mentioned earlier, see Remark 4.1, the need for a reliable and efficient solver is crucial when using iterative methods to solve the linear systems arising in the PDAS algorithm, where solve tolerances have to be small. Firstly, when exactly solving the matrices in the preconditioning system, in the case of two phases, it was shown, in Theorem 3.1, that GMRES will converge within three iterations. Secondly, in the case of multiple phases, it was shown experimentally, that the use of the precondioner P 3,AM G leads to low GMRES iteration counts on fine meshes. Finally, given the standard blocks used   in this solver, i.e., Multigrid, GMRES and simple smoothers, the proposed approach may immediately be applied in most of the software packages used to solve multiphase variational inequalities.