Non-local Allen-Cahn systems: Analysis and a primal dual active set method

Universit¨ at Regensburg Mathematik

Non-local Allen-Cahn systems: Analysis and a primal dual active set method

Luise Blank, Harald Garcke, Lavinia Sarbu and Vanessa Styles

Preprint Nr. 02/2011

Non-local Allen-Cahn systems: Analysis and a primal dual active set method

Luise Blank

Harald Garcke

Lavinia Sarbu

Vanessa Styles

Abstract

We show existence and uniqueness of a solution for the non-local vector-valued Allen-Cahn variational inequality in a formulation involving Lagrange multi- pliers for local and non-local constraints. Furthermore, we propose and an- alyze a primal-dual active set method for local and non-local vector-valued Allen-Cahn variational inequalities. Convergence of the primal-dual active set algorithm is shown by interpreting the approach as a semi-smooth Newton method and numerical simulations are presented demonstrating its efficiency.

Key words: Allen-Cahn systems, non-local constraints, variational inequality, vector-valued obstacle problems, primal-dual active set method, semi-smooth New- ton method.

AMS subject classification. 35K55, 65K10, 90C33, 90C53, 82C24, 65M60

1 Introduction

The Allen-Cahn equation was introduced by Allen and Cahn [1] and describes the capillarity driven evolution of an interface separating two bulk phases. In the Allen- Cahn model interfaces are modelled to have a thickness of orderε where 0< ε1, and in the interfacial layer a phase field or order parameter rapidly changes its value.

The Allen-Cahn model (or phase field model) has a variety of applications, e.g. in materials science, image processing, biology and geology, see [6, 11, 14, 23, 25, 30, 39].

In many of these applications more than two phases occur. Therefore, the model has been extended to deal with N phases [13, 24]. The phase field takes now the form of a vector-valued function u : Ω×(0, T)→ R^N which describes the fractions of the phases, i.e. each component ofu describes one phase.

The underlying non-convex energy functional is based on the Ginzburg-Landau en- ergy for the vector-valued phase field u ∈R^N

E(u) :=R

Ω γε

2|∇u|²+ ^γ_εψ(u) dx

†Fakult¨at f¨ur Mathematik, Universit¨at Regensburg, 93040 Regensburg, Germany

‡Department of Mathematics, University of Sussex, Brighton BN1 9RF, UK

where Ω⊂R^d is a bounded domain, γ >0 is a parameter related to the interfacial energy andψ is a bulk potential. Since each component of ustands for the fraction of one phase, the phase space for the order parameter u is the Gibbs simplex

G:={ξ∈R^N |ξ ≥0,ξ·1 = 1}.

Hereξ ≥0 means ξ_i ≥0 for all i∈ {1, ..., N}, 1= (1, ...,1)^T and ξ·1=

N

P

i=1

ξ_i. For the bulk potential ψ :R^N →R⁺0 ∪ {∞}we consider the multi-obstacle potential

ψ(ξ) = ψ₀(ξ) +I_G =

ψ₀(ξ) := −¹₂ξ·Aξ for ξ∈G,

∞ otherwise, (1)

where I_G is the indicator function of the Gibbs simplex and A is a symmetric constant N ×N matrix [19]. Let σ_max(A) and σ_min(A) be the largest and lowest eigenvalues and k|Ak| the spectral norm of A. If all eigenvalues of A are negative ψ would be a convex potential. Different phases which correspond to minima of ψ only occur if Ahas at least one positive eigenvalue. We hence assume thatAhas at least one positive eigenvalue; the analysis in this paper would simplify if this were not the case.

Given an initial phase distribution u(.,0) = u₀ : Ω→Gat time t= 0 the interface motion can be modelled by the steepest descent dynamics of E with respect to the L²-norm which results, after suitable rescaling of time, in the following vector-valued Allen-Cahn equation

ε^∂u_∂t =−grad_L2E(u) =γε∆u+ ^γ_ε(Au−µ^∗) (2) where µ^∗ ∈ ∂I_G and ∂I_G denotes the subdifferential of I_G. As for the scalar case, see e.g. [9, 12], this equation leads to the following variational inequality

ε(^∂u_∂t,χ−u) +γε(∇u,∇(χ−u))− ^γ_ε(Au,χ−u)≥0 (3) which has to hold for almost all t and all χ ∈ H¹(Ω) with χ ∈ Ga.e.. Here, we denote by L²(Ω) and H¹(Ω) the spaces of vector-valued functions, (., .) is the standard L² inner product for scalar functions, (v,w) =

N

P

i=1

(v_i, w_i) forv,w∈L²(Ω) and (A,B) =

N

P

i=1 d

P

j=1

(a_ij, b_ij) for matrix-valued functions.

Often one considers systems in which the total spatial amount of the phases are conserved. In this case one studies the steepest descent of E under the constraint R

Ω

−udx = m = (mⁱ) where mⁱ ∈ (0,1) for i ∈ {1, ..., N} is a fixed number and we use the notation R

Ω

−f(x)dx := _|Ω|¹ R

Ω

f(x)dx with |Ω| being the Lebesgue measure of Ω. Here we use the notation m and mⁱ in order to avoid confusion with the mass

vectormand its componentsm_i which is introduced in Section 4. To ensure that all phases are present we require 0 < mⁱ < 1 and

N

P

i=1

mⁱ = 1, where the last condition makes sure that

N

P

i=1

u_i = 1 can be true. We define

G:={v∈H¹(Ω)|v∈G a.e.} and G^m :={v∈ G |R

Ω

−v=m}. Then the interface evolution with mass conservation can be formulated as:

(P_m)For given u(.,0) = u₀ ∈ G^m find u∈L²(0, T;G^m)∩H¹(0, T;L²(Ω)) such that ε(^∂u_∂t,χ−u) +γε(∇u,∇(χ−u))− ^γ_ε(Au,χ−u)≥0 (4) which has to hold for almost all t∈(0, T) and all χ∈ G^m.

This paper is organized as follows. In Section 2 we reformulate (P_m) with the help of Lagrange multipliersµ,λand Λ corresponding to the constraintsu ≥0,R

Ω

−u =m and PN

i=1u_i = 1 respectively. We show existence and uniqueness of a solution u of (P_m) and of the Lagrange multipliers,µ,λ and Λ.

In Sections 3 and 4 we introduce the main ideas of a primal-dual active set strat- egy. We apply the algorithm to a finite element discretization of an implicit Euler- discretization of (P_m). Using that the primal-dual active set method can be refor- mulated as a semi-smooth Newton method [27] we show local convergence of our algorithm.

Finally, in Section 5 we present numerical simulations for the non-local as well as for the local Allen-Cahn variational inequality with three and more phases. Using two model problems for which the explicit solution is known we show the efficiency and accuracy of the proposed method. We also discuss numerically how the primal-dual active set method depends on mesh parameters as well as on the number of phases.

2 Existence theory

In this section we show existence and uniqueness to the vector-valued Allen-Cahn variational inequality with integral constraints. As a first step we reformulate the problem (P_m) with the help of (scaled) Lagrange multipliers µ corresponding to the inequality constraint u≥0, Λ corresponding to the constraint

N

P

i=1

u_i = 1 andλ corresponding to the constraint R

Ω

−u=m. However, to be more precise, we have to reformulate the integral constraints in order to guarantee that the constraints are linearly independent, which is required to obtain uniqueness of the Lagrange multi- pliers. Let us therefore consider the sum constraint

N

P

i=1

ui = 1, which is equivalent

to (

N

P

i=1

u_i−1, v) = 0 for all v ∈L²(Ω), and the integral constraint R

Ω

−u =m, which is equivalent to (u−m,e_i) = 0 for all i = 1, ..., N where e_i is the function which is identical to 1 in the i-th component and 0 otherwise. Then for v ≡ 1 we have (

N

P

i=1

u_i−1,1) =

N

P

i=1

(u_i−mⁱ,1) =

N

P

i=1

(R

Ω

−u−m,e_i) which reveals the linear dependence of the constraints. For the reformulation we observe that if

N

P

i=1

u_i = 1 it follows that

N

P

i=1

(R

Ω

−u_i−mⁱ) = 0, i.e.

R

Ω

−u−m∈S :={v∈R^N |PN

i=1v_i = 0}. Noting the above it follows that if PN

i=1ui = 1 andPS(R

Ω

−u−m) =0 it follows that R

Ω

−u=m.

Since we have PN

i=1u_i = 1 we substitute the integral constraint R

Ω

− u = m by P_S(R

Ω

−u−m) =0 and we take λ(t)∈S, for almost every t, to be the corresponding Lagrange multiplier.

Lemma 2.1 Let T be a positive time and let Ω⊂R^d be a bounded domain which is either convex or fulfills∂Ω∈C^1,1. If there existu ∈L²(0, T;H²(Ω))∩H¹(0, T;L²(Ω))∩

L²(0, T;G^m),µ∈L²(0, T;L²(Ω)),λ∈L²(0, T;S)andΛ∈L²(0, T;L²(Ω))such that ε^∂u_∂t −γε∆u− ^γ_εAu− ¹_εµ− ¹_εΛ1− ¹_ελ=0 a.e. in Ω_T := Ω×(0, T), (5) u(0) =u₀, ^∂u_∂ν =0 a.e. on ∂Ω×(0, T), (6)

N

X

i=1

u_i = 1, u≥0, µ≥0 a.e. in Ω_T, (7) P_S(R

Ω

−u−m) =0, (µ,u) = 0 for almost all t∈(0, T), (8) then u solves (Pm).

Proof: We have u(·, t)∈ G^m for almost all tand now chooseχ∈ G^m. Multiplying (5) by (χ−u) and integrating by parts and noting (6) gives

Z

Ω

(ε^∂u_∂t − ^γ_εAu− ¹_εµ−¹_εΛ1− ¹_ελ)·(χ−u) + Z

Ω

γε∇u· ∇(χ−u) = 0

for a.e. t∈ (0, T). Using the property χ≥0 and (7) - (8) gives R

Ω

µ·(χ−u)≥0.

Using (8) and R

Ω

−χ=mwe have R

Ω(χ−u) =0. The fact that λ is independent of

x∈Ω yields

Z

Ω

λ·(χ−u) = 0.

Since 1·(χ−u) = 0 also R

ΩΛ1·(χ−u) = 0. Hence we obtain for all χ∈ G^m and almost all t∈(0, T)

Z

Ω

(ε^∂u_∂t − ^γ_εAu)·(χ−u) + Z

Ω

γε∇u· ∇(χ−u)≥0

and hence u solves (P_m).

To show the existence of (u,µ,λ,Λ) we now regularize ψ, i.e. I_G, and therefore substitute µ^∗ ∈ ∂I_G in (2). First, we introduce the following regularization of ψ₀(ξ) +I_{ξ≥0}, whereI_{ξ≥0} is the indicator function of the set{ξ ∈R^N |ξ ≥0},

ψ_δ(ξ) =ψ₀(ξ) + ¹_δψ(ξ)ˆ (9) where

ψ(ξ) =ˆ

N

X

i=1

(min(ξi,0))². (10)

Similar regularizations were used in [3, 4, 19].

In order to deal with the constraints

N

P

i=1

u_i = 1 and R

Ω

−u = m we project Dψ_δ or- thogonally onto the corresponding tangent space. This projection Pcan be realized by successive orthogonal projections P^R and P^P where

P^Rv:=v− Z

Ω

−v, P^Pv:=v−(P−v)1 and P−v:= _N¹

N

P

i=1

v_i for all v∈R^N. We note that P=P^RP^P =P^PP^R. This results in the following regularized version of the Allen-Cahn equation

ε^∂u_∂t^δ −γε∆u_δ+^γ_εP^RP^PDψ_δ(u_δ) = 0. (11) Equivalently we have to solve the following problem:

(P^δ_m)Givenu₀ ∈ G^m findu_δ ∈H¹(0, T;L²(Ω))∩L²(0, T;H¹(Ω))such thatu_δ(.,0) = u₀ and

ε(^∂u_∂t^δ,χ) +γε(∇u_δ,∇χ) + ^γ_ε(P^P(Dψ_δ(u_δ)),P^Rχ) = 0 (12) for all χ∈H¹(Ω) and almost all t∈(0, T).

Choosingχ=χ1, for anyχ∈H¹(Ω), in (12) gives that a solutionuδ of (P^δ_m) fulfills

∂

∂t(

N

P

i=1

(uδ)i)−∆(

N

P

i=1

(uδ)i) = 0 in a weak sense. Using the facts that u0 ∈ G^m and

that the initial value problem to the parabolic problem∂_tv−∆v = 0 with Neumann boundary conditions is uniquely solvable now gives

N

X

i=1

(uδ(x, t))i = 1 for almost all (x, t)∈ΩT. (13) Hence P^P(u_δ) =u_δ. Choosing constant test functions in (12) gives

d dt

R

Ω

uδ =0 (14)

and hence the total masses of the components of u_δ are preserved, i.e. R

Ω

−u_δ = m.

Before we show that the regularized problem (P^δ_m) has a unique solution, let us state some properties of ψ_δ and Dψ_δ.

We have Dψ_δ(ξ) = ¹_δφ(ξ)ˆ −Aξ with φ(ξ) :=ˆ Dψ(ξ) = ( ˆφ(ξ_i))^N_i=1 where ˆφ(r) = 2 min(r,0) = 2[r]− and [·]−:= min(·,0).

Furthermore:

• For all r, s∈R

0≤( ˆφ(r)−φ(s))(rˆ −s). (15)

• For all ξ,η ∈R^N

(ξ−η)·Dψδ(η) ≤ ¹_δψ(ξ)ˆ − ¹_δψ(η)ˆ −(ξ−η)·Aη

=ψ_δ(ξ)−ψ_δ(η) + ¹₂(ξ−η)·A(ξ−η)

≤ψ_δ(ξ)−ψ_δ(η) + ¹₂σ_max(A)kξ−ηk²,

(16)

where we have used that ˆψ is convex and the identity

−2(ξ−η)·Aη=η·Aη−ξ·Aξ+ (ξ−η)·A(ξ−η).

• For all ξ∈ M:={ξ∈R^N :

N

P

i=1

ξ_i = 1} and δ≤δ₀ := _4N(N−1)¹2σmax(A) we have that

ψ_δ(ξ)≥ _2δ¹

N

X

i=1

[ξ_i]²₋−C(N, σ_max(A)). (17) This follows from

ψ_δ(ξ) ≥ ¹_δ

N

X

i=1

[ξ_i]²₋−

N

X

i=1

σ_max(A)ξ_i²

≥ _2δ¹

N

X

i=1

[ξ_i]²₋+ _2δ¹[ξ_m]²₋−N σ_max(A)ξ_M² (18)

where ξ_m := min

i=1,...,Nξ_i and ξ_M is chosen such that|ξ_M|= max

i=1,...,N|ξ_i|.

Since ξ∈ M it follows that

[ξ_m]− ≤ξ_M ≤1−(N −1)[ξ_m]−

and

ξ_M² ≤2(1 + (N −1)²[ξ_m]²₋) and hence (17) follows from (18) for allδ ≤δ0.

Theorem 2.1 LetΩ⊂R^d be a bounded domain and assume that eitherΩis convex or ∂Ω ∈ C^1,1. Let u₀(x) ∈ H¹(Ω) with u₀ ≥ 0, R

−u₀ = m > 0 and

N

P

i=1

(u₀)_i = 1 a.e. in Ω. Then there exists a unique solution u_δ to (P^δ_m) for all δ ∈ (0, δ₀] and a constant C >0 which does not depend on δ such that

kuδk_L^∞_(0,T;H¹_(Ω))+kuδk_H¹_(0,T;L²_(Ω)) ≤C, (19) k[u_δ]−k_L^∞_(0,T;L²_(Ω)) ≤Cδ^1/2, (20)

1

δkφ(uˆ δ)k_L²_(0,T;L²_(Ω))≤C, (21)

and

ku_δk_L²_(0,T;H²_(Ω)) ≤C. (22) Proof: The existence of a solution to (P^δ_m) follows by using a Galerkin approxi- mation, a priori estimates and compactness arguments. Using the assumptions on Ω and the growth property ofDψ_δ regularity theory gives u_δ ∈L²(0, T;H²(Ω)), using methods from [21] and [26].

Now we show that the solution is unique. Therefore assume that (P^δ_m) has two solutions u¹_δ,u²_δ, subtracting them and choosing χ≡d :=u¹_δ−u²_δ gives

ε(^∂d_∂t,d) +γεk∇dk²+^γ_ε(P^P(¹_δφ(uˆ ¹_δ)−¹_δφ(uˆ ²_δ)−Ad),P^Rd) = 0.

Since P−d = 0 andR

Ω

−d=0 it follows that

ε 2

d

dtkdk²+γεk∇dk²+_δε^γ(φ(uˆ ¹_δ)−φ(uˆ ²_δ),d) = ^γ_ε(Ad,d)≤ ^γ_εσ_max(A)kdk² where k · kdenotes the L² -norm. With (15) we obtain

ε 2

d

dtkdk² ≤ ^γ_εσ_max(A)kdk².

Now Gr¨onwall’s inequality gives that kdk² = 0 and thus u¹_δ ≡u²_δ.

Choosing χ ≡ ∂u_δ/∂t in (12) and using P^PP^R∂u_∂t^δ = ^∂u_∂t^δ which is due to (13) and (14) we obtain

εk^∂u_∂t^δk²+^γε_{2 dt}^dk∇u_δk²+ ^γ_ε(Dψ_δ(u_δ),^∂u_∂t^δ) = 0.

Integrating this equation over (0, t) and rearranging gives for almost allt Z t

0

εk^∂u_∂t^δ(s)k²ds+ ^γε₂ k∇uδ(t)k²+^γ_ε(ψδ(uδ(t)),1) = ^γε₂ k∇u0k²+^γ_ε(ψδ(u0),1)≤C (23) whereC does not depend onδ which follows sinceu₀ ∈ G^m impliesψ_δ(u₀) =ψ(u₀).

In particular, it follows from (17) and (23) that

|(ψ_δ(u_δ)(t),1)| ≤C. (24)

Using (17) gives

N

X

n=1

k[(uδ(t))n]−k² = (

N

X

n=1

[(uδ(t))n]²₋,1)≤Cδ

for almost all t. In conclusion we have that k[u_δ(t)]−k ≤ Cδ^1/2 for almost all t∈(0, T) and thus (20) follows.

Furthermore, from (23) and (24) it follows that k∇u_δ(t)k ≤ C for almost all t and using the Poincar´e inequality

kηk ≤C_P(k∇ηk+|(η,1)|) for all η∈H¹(Ω) (25) gives that ku_δ(t)k²_H1(Ω) ≤ C for almost all t ∈ (0, T) and thus u_δ is uniformly bounded in L^∞(0, T;H¹(Ω)). From (23) and (24) it also follows that ^∂u_∂t^δ

δ>0 is uniformly bounded in L²(0, T;L²(Ω)). Hence (19) is proven.

Now we bound ¹_δφ(uˆ _δ) in L².

Setting χ≡ ¹_δP^Pφ(uˆ _δ) in (12) gives

ε

δ(^∂u_∂t^δ,P^Pφ(uˆ _δ)) + ^γε_δ(∇u_δ,∇P^Pφ(uˆ _δ)) + _εδ^γ(P^P(Dψ_δ(u_δ)),P^RP^Pφ(uˆ _δ)) = 0.

Using

N

P

i=1

∇(u_δ)_i = 0 a.e. which follows from (13) andP^P(Dψ_δ(u_δ))≡P^P(¹_δφ(uˆ _δ))−

P^P(Au_δ) we compute

ε

δ(^∂u_∂t^δ,P^Pφ(uˆ _δ)) + ^γε_δ(∇u_δ,∇φ(uˆ _δ)) + _δ^γ2ε(P^Pφ(uˆ _δ),P^RP^Pφ(uˆ _δ))

= _εδ^γ(P^PAu_δ,P^RP^Pφ(uˆ _δ)).

Noting that (P^Rv,R

Ω

−v) = 0 for any v ∈L²(Ω) and using (14) we obtain

γε

δ (∇uδ,∇φ(uˆ δ)) + _δ^γ2εkP^RP^Pφ(uˆ δ)k²

≤ ^ε_δ|(^∂u_∂t^δ,P^RP^Pφ(uˆ δ))|+ _εδ^γ|(P^PAuδ,P^RP^Pφ(uˆ δ))|.

Since ˆφ is non-decreasing we have that RT

0 (∇u_δ,∇φ(uˆ _δ))dt≥ 0 and hence Young’s inequality and the uniform estimates onu_δ ∈L^∞(0, T;H¹(Ω)) and∂_tu_δ ∈L²(0, T;L²(Ω)) yield

1

δ²kP^RP^Pφ(uˆ _δ)k²_L2(0,T;L²(Ω)) ≤C. (26) Choosing χ=u_δ in (12) and using (13) we obtain

0 = ε(^∂u_∂t^δ,u_δ) +γε(∇u_δ,∇u_δ) + ^γ_ε(P^PDψ_δ(u_δ),P^Ru_δ)

=ε(^∂u_∂t^δ,u_δ) +γε(∇u_δ,∇u_δ) + ^γ_ε(Dψ_δ(u_δ),P^Ru_δ). (27) From (27) and (16) it follows for any constantξ ∈R^N and for a.e. t∈(0, T) that

γ

ε(Dψδ(uδ),ξ−R

Ω

−uδ) = ^γ_ε(Dψδ(uδ),ξ−uδ)−γεk∇uδk²−ε(^∂u_∂t^δ,uδ)

≤ ^γ_ε(ψ_δ(ξ)−ψ_δ(u_δ),1) + ^γσmax_2ε^(A)kξ−u_δk²+εk^∂u_∂t^δkku_δk.

Setting now ξ ≡ (R

Ω

− u_δ)± βeⁿ = m ± βeⁿ, where eⁿ is the n-th unit vector, n = 1, ..., N, and β ∈ (0,1) such that β1 <R

Ω

−u₀ =m < 1 we have ˆψ(ξ) = 0 and obtain

γ

δε(φ(uˆ _δ),±βeⁿ) ≤ ^γ_ε(Au_δ,±βeⁿ)−_2ε^γ(A(m±βeⁿ),m±βeⁿ) + _2ε^γ(Au_δ,u_δ)

−_δε^γ( ˆψ(uδ),1) + ^γk|Ak|_2ε km±βeⁿ−uδk²+εk^∂u_∂t^δkkuδk

≤ ^βγk|Ak|

√

|Ω|

ε ku_δk+ ^3γk|Ak|_2ε km±βeⁿk²+^3γk|Ak|_2ε ku_δk² +εk^∂u_∂t^δkku_δk

where we used that ˆψ(u_δ)≥0. The above estimate gives the existence of a constant C which does not depend onδ such that for almost all (t∈(0, T)

1 δ|R

Ω

− {(φ(uˆ _δ)}_n| ≤C(1 +ku_δk²+k^∂u_∂t^δkku_δk)

for alln = 1, ..., N and hence taking the L²(Ω)-norm of the constanct vectorR

Ω

−φ(uˆ δ) we obtain

1 δkR

Ω

−φ(uˆ _δ)k ≤C(1 +ku_δk²+k^∂u_∂t^δkku_δk) (28) for almost all t, where C depends on ε, γ, N,A,|Ω| and β but not on δ. Squaring (28) gives after integration overt ∈(0, T) and noting (19) that

k¹_δR

Ω

−φ(uˆ _δ)k²_L2(0,T;L²(Ω)) ≤C. (29) Combining (26) and (29) gives that

1

δkP^Pφ(uˆ _δ)k_L²_(0,T;L²_(Ω))≤C. (30)

From (13) it follows that for almost all (x, t)∈Ω_T there exists ann(x, t)∈ {1, ..., N} such that {u_δ(x, t)}_n(x,t) ≥0 and hence {φ(uˆ _δ(x, t))}_n(x,t)= 0. This implies

(P^Pφ(uˆ _δ(x, t)))_n(x,t) =−P−φ(uˆ _δ(x, t)) and noting (30) we obtain that

1 δkP

−φ(uˆ _δ)k_L²_(0,T;L²_(Ω)) ≤ ¹_δkP^Pφ(uˆ _δ)k_L²_(0,T;L²_(Ω))≤C.

Together with (30) we obtain

1

δkφ(uˆ _δ)k_L²_(0,T;L²_(Ω)) ≤ ¹_δkP^Pφ(uˆ _δ)k_L²_(0,T;L²_(Ω))+ ¹_δkP

− φ(uˆ _δ)1k_L²_(0,T;L²_(Ω)) ≤C.

We refer to [3] where similar arguments have been used in the context of a Cahn- Hilliard system with a logarithmic free energy. Finally, the fact that theL²(0, T;H²(Ω))- norm ofu_δ is uniformly bounded in δ, see (22), follows from (12), (19), (21) and by applying elliptic regularity theory, see [26], on time slices.

We are now in a position to show an existence and uniqueness theorem for the original problem.

Theorem 2.2 LetΩ⊂R^d be a bounded domain and assume that eitherΩis convex or fulfills ∂Ω ∈ C^1,1. Let u₀ ∈ H¹(Ω) with u₀ ≥ 0,0 < R

Ω

− u₀ = m < 1 and

N

P

i=1

(u₀)_i = 1 a.e. in Ω. Then there exists a unique solution (u,µ,λ,Λ) to (5) - (8) with the following properties:

u(x,0) =u₀(x) for almost all x∈Ω, (31) u∈L^∞(0, T;H¹(Ω))∩H¹(0, T;L²(Ω))∩L²(0, T;H²(Ω)), (32) µ∈L²(0, T;L²(Ω)), (33) λ∈L²(0, T) and

N

X

i=1

λ_i = 0 for almost all t∈(0, T), (34) Λ∈L²(0, T;L²(Ω)). (35) Proof: As the bounds (19), (21) and (22) are independent ofδ, it follows that there exists au∈L^∞(0, T;H¹(Ω))∩H¹(0, T;L²(Ω))∩L²(0, T;H²(Ω)), and a subsequence {u_δ⁰} which converges to u as δ⁰ →0

a) in L^∞(0, T;H¹(Ω)) weak-star,

b) in H¹(0, T;L²(Ω))∩L²(0, T;H²(Ω)) weakly, c) in L²(0, T;L²(Ω)) strongly,

d) almost everywhere in Ω×(0, T)

(36)

where c) follows from a) and b), see [38].

Since u_δ ∈L²(0, T;H²(Ω)) we can use the strong formulation of (12) and obtain ε^∂u_∂t^δ −γε∆u_δ− ^γ_εAu_δ− ¹_εµ_δ− ¹_εΛ_δ1− ¹_ελ_δ =0 (37) where Λ_δ :=γDψ_δ(u_δ) = ^γ_δP−φ(uˆ _δ)−γP−Au_δ,λ_δ :=γR

ΩP^PDψ_δ(u_δ) = ^γ_δR

Ω

−φ(uˆ _δ)− γR

Ω

−Au_δ−R

Ω

−Λ_δ1 and µ_δ :=−^γ_δφ(uˆ _δ)≥0.

Note that Λ_δ ∈ L²(0, T;L²(Ω)), λ_δ ∈ L²(0, T) and µ_δ ∈ L²(0, T;L²(Ω)) are uni- formly bounded, see (19) and (21). Hence there exist Λ ∈ L²(0, T;L²(Ω)),λ ∈ L²(0, T) and µ∈L²(0, T;L²(Ω)) such that for a subsequence

Λ_δ*Λ inL²(0, T;L²(Ω)) as δ⁰ →0 λ_δ *λ inL²(0, T;L²(Ω)) as δ⁰ →0 µ_δ⁰ *µ inL²(0, T;L²(Ω)) as δ⁰ →0.

Consequently, equation (5) holds for the limit. Furthermore, from (36), (13) and (14) it follows that

N

P

i=1

u_i = 1 and _dt^d R

Ω

u=0 and hence R

Ω

−u =m. Taking the limit δ⁰ →0 in (20) gives that [u]−=0 and thus u(x, t)≥0 fora.e. (x, t)∈ΩT.

Since µ is the weak limit of functions which are componentwise non-negative we obtain µ≥0 a.e. in Ω_T. In order to show that (µ,u) = 0 we first note that

(µ_δ,u_δ) =−^γ_δ(φ(uˆ _δ),u_δ)≤0,

and using that u_δ → u and µ_δ * µ in L²(0, T,L²(Ω)) it follows that (µ,u) ≤ 0.

However, since u ≥ 0 and µ≥ 0 we have that (u,µ) = 0 a.e. in (0, T). Note that

N

P

i=1

(λ_δ)_i = 0 and hence

N

P

i=1

λ_i = 0, i.e. λ∈S.

It remains to show uniqueness. Assume that there are two solutions (u¹,µ¹,λ¹,Λ¹) and (u²,µ²,λ²,Λ²). Then we define ¯u = u¹ −u²,µ¯ = µ¹−µ². Multiplying the difference of the equation (5) for u¹ and u² with ¯u gives, after integration and using 1·u¯ = 0 and R

−u=0 , that ε_dt^d

Z

Ω

|¯u|² +γε Z

Ω

|∇u|¯² −¹_ε Z

Ω

(µ¹−µ²)·(u¹−u²)≤ ^γσmax_ε^(A) Z

Ω

|u|¯ ². (38)

The complementarity conditions (7)-(8) imply that (µ¹ −µ²)·(u¹ −u²) ≤ 0 and hence we deduce that

ε_dt^d Z

Ω

|u|¯ ²+γε Z

Ω

|∇¯u|² ≤ ^γσmax_ε^(A) Z

Ω

|u|¯ ².

Using a Gr¨onwall argument now gives uniqueness ofu. Henceµ+λ+Λ1is uniquely given through equation (5).

Now we show the uniqueness of the Lagrange multipliers λ,Λ and µ. For what follows we fix t ∈ (0, T) such that u(t) ∈ H²(Ω) and define the inactive sets I_i :=

{x∈Ω|u_i(x, t)>0}, the interface between phasesi and j as I_ij :=I_i∩I_j and the measure|I_ij| of I_ij.

We claim: λ_i−λ_j is uniquely defined for all pairs (i, j) with |I_ij|>0.

Using (5) and recalling that e^k is the k-th unit vector we obtain for (i, j) with

|I_ij|>0 that

(ε^∂u_∂t −γε∆u−^γ_εAu−¹_ελ)·(eⁱ−e^j) = 0 on Iij (39) where we have used that µ_i =µ_j = 0 onI_ij. We hence conclude

λ_i−λ_j = _|I¹

ij|

Z

Iij

h

ε^2∂u_∂tⁱ −γε²∆u_i−γ(Au)_i−ε^2∂u_∂t^j +γε²∆u_j +γ(Au)_ji .

This implies that the differenceλ_i−λ_j is uniquely defined if there exists an interface between phases iand j, i.e. if |I_ij|>0.

We now define a graph over {1, ..., N} with the edges E ={{i, j}:|Iij|>0}. If the graph is connected, which we show in the following, the differencesλ_i−λ_j are for all i, j ∈ {1, ..., N} uniquely defined. Together with the condition

N

P

i=1

λ_i = 0 we obtain the uniqueness of λ.

In order to show that the graph is connected, we define the following sets of indices L ={i∈ {1, ..., N}: there is a path from 1 to i} and M={1, ..., N} \ L.

We need to show that M=∅ and therefore we assume M 6=∅. We set

v =X

i∈L

u_i and w= X

j∈M

u_j

and note that v ≥0, w ≥0 and v+w= 1. Now one observes that the set A :={x∈Ω :v(x)>0 and w(x)>0}

has measure zero. This is true because otherwise there exist i∈ L and j ∈ M such that |I_ij|>0 which contradicts the definition ofL and the assumption M 6=∅. We hence obtain that v only attains the values 0 and 1. Since M 6= ∅ we obtain that v is not constant. Since an H¹-function that attains finitely many values has to be constant we obtain a contradiction. Hence M=∅ and the graph is connected.

Now we show the uniqueness of Λ. Since PN

i=1u_i = 1 and u ≥ 0, we can find for any x₀ ∈ Ω an i ∈ {1, ..., N} such that x₀ ∈ I_i and |I_i| > 0. On I_i we know that µ_i = 0 and hence we can define

Λ = (ε^2∂u_∂t −γε²∆u−γAu−λ)_i on I_i. (40)

The Lagrange multiplier Λ is well defined and unique since

(ε^∂u_∂t −γε∆u−^γ_εAu− ¹_ελ)_i(x, t) = (ε^∂u_∂t −γε∆u− ^γ_εAu−¹_ελ)_j(x, t)

for almost every x ∈ I_ij and since for almost all x ∈ Ω, there is an i ∈ {1, . . . , N} such thatui(x, t)>0, i.e. x∈Ii. Having shown uniqueness ofu,λand Λ uniqueness of µ follows from equation (5).

Remark 2.1 i) It can be shown that a solution to(P_m)is unique. Hence we can

conclude that for all solutions to(Pm)there exist Lagrange multipliersµ,λ,Λ such that (5) - (8) hold. Furthermore, problem (5)-(8) is equivalent to (P_m).

ii) In [35] the existence of a solution for vector-valued Allen-Cahn variational in- equalities without volume constraints is shown by using a representation of the Lagrange multipliers which cannot be used directly for a numerical approach.

There u_N is substituted by 1−PN−1

i=1 u_i and a system of parabolic variational inequalities in R^N−1 is considered.

3 Primal-dual active set approach

For the numerical approximation of solutions u of (P_m) we introduce a primal- dual active set method or equivalently a semi-smooth Newton method [7, 27]. Both are well known in the context of optimization with partial differential equations as constraints. We present a time discretization of the Allen-Cahn system and reformulate the complementarity conditions using primal-dual active sets. Finally, even though the method is not applicable to the time discretized problem, we present for ease of understanding the idea of the resulting iterative solution procedure for the time discretized problem, which will be applied to the fully discretized problem in the next section.

We denote the time step byτ, which can be a variable time step,t₀ = 0,t_n:=tn−1+τ and uⁿ⁻¹ :=u(., tn−1). For simplicity we denote by u the time discrete solution at time t_n. Then possible time discretizations of (P_m) are given as (semi-)implicit or explicit Euler-discretizations. Explicit Euler discretizations for Allen-Cahn obstacle problems have been used for example in [12, 20, 23, 24]. Numerical analysis for (semi-) implicit discretizations of the Allen-Cahn model has been performed in the papers [15, 22, 31, 32, 33, 34] and in works cited in these papers. Fully implicit discretizations are the most accurate, see e.g. [9], but due to the non-monotonicity for large time steps they can be either very expansive or they are not uniquely solvable. It will turn out that this is not the case for the primal dual active set approach as for well developed interfaces also larger time steps can be used, see e.g.

Remark 4.3. In this paper we focus on the implicit discretization of the vector-valued Allen-Cahn obstacle problem leading to the following formulation:

(P^τ_m)Given uⁿ⁻¹ ∈ G^m find u=uⁿ ∈ G^m such that

ε

τ(u−uⁿ⁻¹,η−u) +γε(∇u,∇η− ∇u)≥ ^γ_ε(Au,η−u) ∀η∈ G^m. (41) This discretization can also be seen as the Euler-Lagrange equation of an implicit time discretization of the L₂ gradient flow of the energyE, which is given as

u∈Gmin^mE(u) :=

Z

Ω

_γε

2 |∇u|²+^γ_εψ(u) dx+_2τ^εku−uⁿ⁻¹k²_L2. (42) As in Lemma 2.1 one can reformulate (P^τ_m) by using scaled Lagrange-multipliersµ∈ L²(Ω) for the inequality constraintu≥0,λ∈S for the constraintP_S(R

Ω

−u−m) =0 and Λ∈L²(Ω) for the sum constraint PN

i=1u_i = 1 to obtain:

ε²

τ (u−uⁿ⁻¹)−γε²∆u−γAu−µ−λ−Λ1 =0 a.e. in Ω, (43)

∂u

∂ν =0 a.e. on∂Ω, (44) P_S(R

Ω

−u−m) =0 (45)

N

P

i=1

u_i = 1 (46)

together with the complementarity conditions

u ≥0 a.e. in Ω, µ≥0 a.e. in Ω, (µ,u) = 0. (47) Now the idea is to reformulate the complementarity conditions using active sets based on the primal variable u and the dual variable µ. Then, for any c > 0, (47) is equivalent to the following: For all i∈ {1, ..., N}:

ui = 0 a.e. inAi; µi = 0 a.e. in Ii := Ω\ Ai; (48) where the primal-dual active sets are given by

A_i ={x∈Ω|cu_i(x)−µ_i(x)<0}. (49) If the sets Ai are known, we can determine the solution as follows. First we set u_i = 0 onA_i and µ_i = 0 on I_i, see (48). It now turns out that u_i only needs to be determined in points inI_i in which another component is inactive, i.e. on

D_i =I_i∩([

j6=i

I_j).

OnI_i\ D_i we observe thatiis the only inactive component and hence the constraint PN

j=1u_j = 1 leads to

u_i = 1 on I_i\ D_i. (50)

Figure 1: The computational effort in the primal-dual active set method is restricted to the diffuse interface.

Defining now the total diffuse interface region, see Figure 1, as D:=

N

[

i=1

D_i ⊇ {x∈Ω| there exists i∈ {1, . . . , N} s.t. 0 < u_i(x)<1}

we need to solve the following system: Find, for i = 1, . . . , N, u_i on D_i, Λ on D, λ∈R^N such that fori= 1, . . . , N

−γε²∆u_i+^ε_τ²(u_i−uⁿ⁻¹_i )−λ_i−Λ−γ(Au)_i = 0, a.e. on D_i, (51)

∂u_i

∂ν = 0 a.e. on ∂D_i∩∂Ω, u_i = 0 on∂D_i∩∂A_i, (52)

N

X

i=1

u_i = 1 on D and P_S(R

Ω

−u−m) = 0,

N

X

i=1

λ_i = 0. (53)

Then we have to determine Λ on Ω\ D. Since in each point x ∈ Ω at least one component is inactive there exists for a given x ∈ Ω\ D an i such that x ∈ I_i\ D and given µ_i = 0 onI_i we set

Λ =−γε²∆u_i+^ε_τ²(u_i−uⁿ⁻¹_i )−λ_i−γ(Au)_i. (54) Then Λ is completely determined. Finally we set

µ_i =−γε²∆u_i+^ε_τ²(u_i−uⁿ⁻¹_i )−λ_i−γ(Au)_i−Λ onA_i. (55) This leads to the idea of the Primal-Dual Active Set (PDAS) algorithm:

Given initial active sets A⁰_i for alli∈ {1, ..., N} iterate the following steps fork≥0 (where we define I_i^k,D_i^k,D^k analogous to the discussion above)

i) Set u^k_i = 0 on A^k_i and µ^k_i = 0 on I_i^k for all i∈ {1, ..., N}, u^k_i = 1 on I_i^k\ D^k_i. ii) Solve the coupled system (51)-(53) for λ^k, Λ^k on D and u^k_i on D_i^k for all

i∈ {1, ..., N}.

iii) Determine Λ^k on Ω\ D using (54).

iv) Determineµ^k_i on A^k_i using (55) for all i∈ {1, ..., N}.

v) Determine the new active sets A^k+1_i = {x ∈ Ω|cu^k_i(x)−µ^k_i(x) < 0} for all i∈ {1, ..., N}.

vi) Stop the iteration if A^k+1_i =A^k_i for all i∈ {1, ..., N}, otherwise set k =k+ 1 and goto 1.

Except for the sign conditions for u and µ all conditions (44)-(48) hold in each it- eration. As mentioned in the beginning of this section we cannot apply the method to the time discretized Allen-Cahn variational inequality. The reason is that al- though one can show the existence of the Langrange-multipliers and the regularity µ ∈ L²(Ω) this regularity does in general not hold in each iteration of the PDAS- algorithm. Then the multipliers may still exist but are only measures. This effect is also known for obstacle problems, see [29], and is discussed in more detail for Cahn-Hilliard problems in [8]. Therefore, the pointwise definition of the active sets A^k_i is not possible. However we show in the next section that the application of the PDAS-method to the fully discretized problem is possible.

The feature that the method is not applicable in function space may lead to mesh dependence of the PDAS iteration numbers for a fixed time step. Analysis for mesh independence is still lacking and needs further research. However, our numerical in- vestigations clearly indicate mesh independence if the time and space discretizations are reduced simultaneously, see Section 5.

4 Finite element approximation

For space discretization we employ a finite element approximation which we present in this section. Furthermore, we present the PDAS-algorithm for the fully discretized system.

4.1 Notation

For simplicity we assume that Ω is a polyhedral domain. LetTh be a regular triangu- lation of Ω into disjoint open simplices, i.e. in particular Ω =∪_T∈T_hT. Furthermore, we define h:= max_T∈T_h{diam T} to be 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 S_h :=n

ϕ ∈C⁰(Ω) ϕ

T

∈P₁(T) ∀ T ∈ T_ho

⊂H¹(Ω),

where we denote by P₁(T) the set of all affine linear functions on T. Furthermore we denote the standard nodal basis functions of Sh by χj for all j ∈ J and we set S_h = (S_h)^N. Then u_j ∈ R^N for j ∈ J denotes the coefficients of the basis representation of u_h inS_h which is given by u_h =P

j∈J u_jχ_j. In order to derive a discretization of the Allen-Cahn model we set

G_h :={χ∈S_h| χ≥0 and

N

X

i=1

(χ_i)_j = 1 ∀j ∈ J } and

G_h^m :={η ∈ G_h| R

Ω

−η =m}.

Here (χi)j denotes the i-th component χi of χ at thej-th node. We introduce also the lumped mass semi-inner product (f, g)_h = R

ΩI_h(f g) instead of (f, g), where I_h :C⁰(Ω)→S_h is the standard interpolation operator such that (I_hf)(p_j) =f(p_j) for all nodes j ∈ J.

Defining m_j := (1, χ_j) we have R

Ω

− u_i = P

j∈J m_j(u_i)_j/P

j∈J m_j where u_i ∈ S^h and i ∈ {1, ..., N}. Moreover, we define the stiffness matrix as S := (s_ij) with s_ij = (∇χ_j,∇χ_i), the mass matrix M := ((χ_j, χ_i)_h) = diag(m_j) and the mass vector m := (m_j)j∈J. Also we denote the entries of A by a_ij, i, j = {1, . . . , N}.

Recall that the total spatial amount to be conserved is denoted bym= (mⁱ)_i=1,...,N which should not be confused with the mass vector m.

4.2 Finite element approximation and the PDAS-algorithm

We now introduce the following finite element approximations of (P^τ_m) given by (41).

In the following we consider a fixed time stepτ =t_n−tn−1 and omit in some places the superscript n:

(P^τ_m,h) Given uⁿ⁻¹_h ∈ G_h^m find u_h =uⁿ_h ∈ G_h^m such that

(_τ^ε(u_h−uⁿ⁻¹_h )− ^γ_εAu_h,χ−u_h)_h+γε(∇u_h,∇(χ−u_h))≥0 ∀χ∈ G_h^m. (56)

Due to the use of piecewise linear finite elements and nodal basis functions the reformulation of (P^τ_m,h) with Lagrange multipliersµ_h ∈S_h can be stated as follows: