A C 1 virtual element method for the Cahn-Hilliard equation with - - PowerPoint PPT Presentation

A C 1 virtual element method for the Cahn-Hilliard equation with polygonal meshes Marco Verani

MOX, Department of Mathematics, Politecnico di Milano Joint work with:

• P. F. Antonietti (MOX Politecnico di Milano, Italy)
• L. Beir˜

ao da Veiga (Universit` a di Milano, Italy)

• S. Scacchi (Universit`

a di Milano, Italy)

Polytopal Element Methods, Georgia Tech, October 26th, 2015

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 1 / 28

Outline

1 Cahn-Hilliard equation 2 C 1-VEM for Cahn-Hilliard 3 Numerical results 4 Conclusions

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 2 / 28

Cahn-Hilliard (CH) problem

Let Ω ⊂ R2 be an open bounded domain. Let ψ(x) = (1 − x2)2/4 and φ(x) = ψ′(x).      ∂tu − ∆

• φ(u) − γ2∆u(t)
• = 0

in Ω × [0, T], u(·, 0) = u0(·) in Ω, ∂nu = ∂n

• φ(u) − γ2∆u(t)
• = 0
• n ∂Ω × [0, T],

where ∂n denotes the (outward) normal derivative and γ ∈ R+, 0 < γ ≪ 1, represents the interface parameter.

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 3 / 28

Cahn-Hilliard (CH) problem

Let Ω ⊂ R2 be an open bounded domain. Let ψ(x) = (1 − x2)2/4 and φ(x) = ψ′(x).      ∂tu − ∆

• φ(u) − γ2∆u(t)
• = 0

in Ω × [0, T], u(·, 0) = u0(·) in Ω, ∂nu = ∂n

• φ(u) − γ2∆u(t)
• = 0
• n ∂Ω × [0, T],

where ∂n denotes the (outward) normal derivative and γ ∈ R+, 0 < γ ≪ 1, represents the interface parameter. CH models the evolution of interfaces specified as the level set of a smooth continuos function u exhibiting large gradients across the interface [Korteweg 1901, Cahn and Hilliard 1958, Ginzburg and Landau

1965, van der Waals 1979]

Many applications: phase separation in binary alloys, tumor growth,

• rigin of Saturn’s rings, separation of di-block copolymers, population

dynamics, image processing, clustering of mussels ...

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 3 / 28

Huge literature on numerical methods for CH...

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 4 / 28

Huge literature on numerical methods for CH... Here we limit ourvselves to FEM as the main properties (and limitations) of these schemes are instrumental to motivate the introduction of our new approach.

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 4 / 28

Huge literature on numerical methods for CH... Here we limit ourvselves to FEM as the main properties (and limitations) of these schemes are instrumental to motivate the introduction of our new approach. CH is a fourth order nonlinear problem ❀ C 1 FEM (see [Elliott Zheng, 1986], [Elliott French, 1987]).

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 4 / 28

Huge literature on numerical methods for CH... Here we limit ourvselves to FEM as the main properties (and limitations) of these schemes are instrumental to motivate the introduction of our new approach. CH is a fourth order nonlinear problem ❀ C 1 FEM (see [Elliott Zheng, 1986], [Elliott French, 1987]). To avoid the well known difficulty met in the implementation of C 1 FEM, another possibility is the use of non-conforming (see [Elliott French, 1989] or discontinuous (see [Wells Kuhl Garikipati, 2006]) methods.

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 4 / 28

Huge literature on numerical methods for CH... Here we limit ourvselves to FEM as the main properties (and limitations) of these schemes are instrumental to motivate the introduction of our new approach. CH is a fourth order nonlinear problem ❀ C 1 FEM (see [Elliott Zheng, 1986], [Elliott French, 1987]). To avoid the well known difficulty met in the implementation of C 1 FEM, another possibility is the use of non-conforming (see [Elliott French, 1989] or discontinuous (see [Wells Kuhl Garikipati, 2006]) methods. Alternatively, use mixed methods (see e.g. [Elliot French Milner, 1989], [Elliott Larsson, 1992] and [Kay Styles Suli, 2009] for the continuous and discontinuous setting, respectively).

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 4 / 28

Huge literature on numerical methods for CH... Here we limit ourvselves to FEM as the main properties (and limitations) of these schemes are instrumental to motivate the introduction of our new approach. CH is a fourth order nonlinear problem ❀ C 1 FEM (see [Elliott Zheng, 1986], [Elliott French, 1987]). To avoid the well known difficulty met in the implementation of C 1 FEM, another possibility is the use of non-conforming (see [Elliott French, 1989] or discontinuous (see [Wells Kuhl Garikipati, 2006]) methods. Alternatively, use mixed methods (see e.g. [Elliot French Milner, 1989], [Elliott Larsson, 1992] and [Kay Styles Suli, 2009] for the continuous and discontinuous setting, respectively). Recently, isogeometric analysis [G´

• mez Calo Bazilevs Hughes, 2008].
• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 4 / 28

Huge literature on numerical methods for CH... Here we limit ourvselves to FEM as the main properties (and limitations) of these schemes are instrumental to motivate the introduction of our new approach. CH is a fourth order nonlinear problem ❀ C 1 FEM (see [Elliott Zheng, 1986], [Elliott French, 1987]). To avoid the well known difficulty met in the implementation of C 1 FEM, another possibility is the use of non-conforming (see [Elliott French, 1989] or discontinuous (see [Wells Kuhl Garikipati, 2006]) methods. Alternatively, use mixed methods (see e.g. [Elliot French Milner, 1989], [Elliott Larsson, 1992] and [Kay Styles Suli, 2009] for the continuous and discontinuous setting, respectively). Recently, isogeometric analysis [G´

• mez Calo Bazilevs Hughes, 2008].

Here we propose a C 1-VEM scheme (see, e.g., [Beir˜ ao da Veiga Brezzi Cangiani Manzini Marini Russo, 2013] for an intro to VEM for second order pbs and [Brezzi Marini, 2013] for plate bending pb.)

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 4 / 28

Weak formulation of CH

a∆(v, w) =

• Ω

(∇2v) : (∇2w) dx ∀v, w ∈ H2(Ω), a∇(v, w) =

• Ω

(∇v) · (∇w) dx ∀v, w ∈ H1(Ω), a0(v, w) =

• Ω

v w dx ∀v, w ∈ L2(Ω),

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 5 / 28

Weak formulation of CH

a∆(v, w) =

• Ω

(∇2v) : (∇2w) dx ∀v, w ∈ H2(Ω), a∇(v, w) =

• Ω

(∇v) · (∇w) dx ∀v, w ∈ H1(Ω), a0(v, w) =

• Ω

v w dx ∀v, w ∈ L2(Ω), and the semi-linear form r(z; v, w) =

• Ω

φ′(z)∇v · ∇w dx ∀z, v, w ∈ H2(Ω).

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 5 / 28

Weak formulation of CH

a∆(v, w) =

• Ω

(∇2v) : (∇2w) dx ∀v, w ∈ H2(Ω), a∇(v, w) =

• Ω

(∇v) · (∇w) dx ∀v, w ∈ H1(Ω), a0(v, w) =

• Ω

v w dx ∀v, w ∈ L2(Ω), and the semi-linear form r(z; v, w) =

• Ω

φ′(z)∇v · ∇w dx ∀z, v, w ∈ H2(Ω). Let V =

• v ∈ H2(Ω) : ∂nu = 0 on ∂Ω
• The weak formulation reads as: find u(·, t) ∈ V such that
• a0(∂tu, v) + γ2a∆(u, v) + r(u; u, v) = 0

∀v ∈ V , u(·, 0) = u0(·).

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 5 / 28

Towards C 1-VEM discretization of CH

Ωh decomposition of Ω into polygons E

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 6 / 28

Towards C 1-VEM discretization of CH

Ωh decomposition of Ω into polygons E Intermediate local space (see [Brezzi Marini, 2013]):

• Vh|E =
• v ∈ H2(E) : ∆2v ∈ P2(E), v|∂E ∈ C 0(∂E), v|e ∈ P3(e) ∀e ∈ ∂E,

∇v|∂E ∈ [C 0(∂E)]2, ∂nv|e ∈ P1(e) ∀e ∈ ∂E

• ,
• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 6 / 28

Towards C 1-VEM discretization of CH

Ωh decomposition of Ω into polygons E Intermediate local space (see [Brezzi Marini, 2013]):

• Vh|E =
• v ∈ H2(E) : ∆2v ∈ P2(E), v|∂E ∈ C 0(∂E), v|e ∈ P3(e) ∀e ∈ ∂E,

∇v|∂E ∈ [C 0(∂E)]2, ∂nv|e ∈ P1(e) ∀e ∈ ∂E

• ,

Linear Operators: D1 ❀ evaluation of v at the vertexes of E; D2 ❀ evaluation of ∇v at the vertexes of E.

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 6 / 28

The operator Π∆

E

The projection operator Π∆

E :

Vh|E → P2(E) is defined by

• a∆

E (Π∆ E vh, q) = a∆ E (vh, q)

∀q ∈ P2(E) ( (Π∆

E vh, q)

)E = ( (vh, q) )E ∀q ∈ P1(E), for all vh ∈ Vh|E where ( (vh, wh) )E=

• ν vertexes
• f ∂E

vh(ν) wh(ν) ∀vh, wh ∈ C 0(E). The operator Π∆

E is uniquely determined by D1 and D2, i.e.

Π∆

E vh depends only on the values of vh and ∇vh at the vertices of E.

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 7 / 28

Virtual local space

Wh|E =

• v ∈

Vh|E :

• E

Π∆

E (vh) q dx =

• E

vh q dx ∀q ∈ P2(E)

• .

[Ahmad Alsaedi Brezzi Marini Russo, 2013]

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 8 / 28

Virtual local space

Wh|E =

• v ∈

Vh|E :

• E

Π∆

E (vh) q dx =

• E

vh q dx ∀q ∈ P2(E)

• .

Properties: The set of operators D1 and D2 constitutes a set of degrees of freedom for the space Wh|E.

[Ahmad Alsaedi Brezzi Marini Russo, 2013]

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 8 / 28

Virtual local space

Wh|E =

• v ∈

Vh|E :

• E

Π∆

E (vh) q dx =

• E

vh q dx ∀q ∈ P2(E)

• .

Properties: The set of operators D1 and D2 constitutes a set of degrees of freedom for the space Wh|E. P2(E) ⊆ Wh|E.

[Ahmad Alsaedi Brezzi Marini Russo, 2013]

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 8 / 28

Virtual local space

Wh|E =

• v ∈

Vh|E :

• E

Π∆

E (vh) q dx =

• E

vh q dx ∀q ∈ P2(E)

• .

Properties: The set of operators D1 and D2 constitutes a set of degrees of freedom for the space Wh|E. P2(E) ⊆ Wh|E. The L2-projection Π0

E : Wh|E → P2(E) defined by

a0

E(Π0 Evh, q) = a0 E(vh, q)

∀q ∈ P2(E), is computable (only) on the basis of the values of the degrees of freedom D1 and D2.

[Ahmad Alsaedi Brezzi Marini Russo, 2013]

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 8 / 28

Global VEM space

Wh =

• v ∈ V : v|E ∈ Wh|E

∀E ∈ Ωh

• Properties:

Wh ⊂ H2(Ω) ❀ conforming VEM solution. The global degrees of freedom are:

values of vh at the vertexes of the mesh Ωh; values of ∇vh at the vertexes of the mesh Ωh.

Dimension of Wh is three times the number of vertexes in the mesh.

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 9 / 28

The operator Π∇

E

We define Π∇

E : Wh|E → P2(E) by

   a∇

E (Π∇ E vh, q) = a∇ E (vh, q)

∀q ∈ P2(E)

• E

Π∇

E vh dx =

• E

vh dx. Since the bilinear form a∇

E (·, ·) has a non trivial kernel (given by the

constant functions) we added a second condition in order to keep the

• perator Π∇

E well defined.

The operator Π∇

E is uniquely determined by D1 and D2

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 10 / 28

Local discrete bilinear forms

We propose the following discrete (and symmetric) local forms a∆

h,E(vh, wh) = a∆ E (Π∆ E vh, Π∆ E wh) + h−2 E sE(vh − Π∆ E vh, wh − Π∆ E wh),

a∇

h,E(vh, wh) = a∇ E (Π∇ E vh, Π∇ E wh) + sE(vh − Π∇ E vh, wh − Π∇ E wh),

a0

h,E(vh, wh) = a0 E(Π0 Evh, Π0 Ewh) + h2 E sE(vh − Π0 Evh, wh − Π0 Ewh),

for all vh, wh ∈ Wh|E, where sE(vh, wh) =

• ν vertexes
• f ∂E
• vh(ν) wh(ν) + (hν)2 ∇vh(ν) · ∇wh(ν)
• being hν some characteristic mesh size lenght associated to the node ν (for

instance the maximum diameter among the elements having ν as a vertex).

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 11 / 28

Properties of the discrete bilinear forms

The symbol † stands for the symbol ∆, ∇ or 0. (Consistency) a†

h,E(p, vh) = a† E(p, vh)

∀p ∈ P2(E), ∀vh ∈ Wh|E, (Stability)1 There exist two positive constants c⋆, c⋆ independent of the element E ∈ {Ωh}h such that c⋆ a†

E(vh, vh) ≤ a† h,E(vh, vh) ≤ c⋆a† E(vh, vh)

∀vh ∈ Wh|E, The bilinear forms a†

h,E(·, ·) are continuous with respect to the

relevant norm: H2 (for † = ∆), H1 (for † = ∇) and L2 (for † = 0).

1 It holds under regularity assumption on the polygonal partition

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 12 / 28

Global discrete bilinear forms

The global discrete bilinear form is assembled as usual a†

h(vh, wh) =

• E∈Ωh

a†

h,E(vh, wh)

∀vh, wh ∈ Wh, with the usual multiple meaning of the symbol †.

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 13 / 28

Discrete semilinear form

Continuous semilinear form r(z; v, w) =

• E∈Ωh

rE(z; v, w) ∀z, v, w ∈ H2(Ω), rE(z; v, w) =

• E

(3z(x)2 − 1)∇v(x) · ∇w(x) dx ∀E ∈ Ωh.

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 14 / 28

Discrete semilinear form

Continuous semilinear form r(z; v, w) =

• E∈Ωh

rE(z; v, w) ∀z, v, w ∈ H2(Ω), rE(z; v, w) =

• E

(3z(x)2 − 1)∇v(x) · ∇w(x) dx ∀E ∈ Ωh. Idea: on each element E, we approximate the term z(x)2 with its average

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 14 / 28

Discrete semilinear form

Continuous semilinear form r(z; v, w) =

• E∈Ωh

rE(z; v, w) ∀z, v, w ∈ H2(Ω), rE(z; v, w) =

• E

(3z(x)2 − 1)∇v(x) · ∇w(x) dx ∀E ∈ Ωh. Idea: on each element E, we approximate the term z(x)2 with its average Discrete local semilinear form rh,E(zh; vh, wh) = φ′(zh)|E a∇

h,E(vh, wh)

∀zh, vh, wh ∈ Wh|E where

• φ′(zh)|E=3|E|−1a0

h,E(zh, zh) − 1

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 14 / 28

Discrete semilinear form

Continuous semilinear form r(z; v, w) =

• E∈Ωh

rE(z; v, w) ∀z, v, w ∈ H2(Ω), rE(z; v, w) =

• E

(3z(x)2 − 1)∇v(x) · ∇w(x) dx ∀E ∈ Ωh. Idea: on each element E, we approximate the term z(x)2 with its average Discrete local semilinear form rh,E(zh; vh, wh) = φ′(zh)|E a∇

h,E(vh, wh)

∀zh, vh, wh ∈ Wh|E where

• φ′(zh)|E=3|E|−1a0

h,E(zh, zh) − 1

The global form is then assembled as usual rh(zh; vh, wh) =

• E∈Ωh

rh,E(zh; vh, wh) ∀wh, rh, vh ∈ Wh.

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 14 / 28

C 1-VEM semi-discrete scheme

Discrete global space with boundary conditions: W 0

h = Wh ∩ V =

• v ∈ Wh : ∂nu = 0 on ∂Ω
• .

The semi-discrete problem is:      Find uh(·, t) in W 0

h such that

a0

h(∂tuh, vh) + γ2a∆ h (uh, vh) + rh(uh, uh; vh) = 0,

uh(0, ·) = u0,h(·), for all vh ∈ W 0

h being u0,h ∈ W 0 h a suitable approximation of u0.

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 15 / 28

Convergence result for the semi-discrete scheme

Theorem (Antonietti, Beirao da Veiga, Sacchi, V., 2015)

Let u be the (sufficiently regular) exact solution of the Cahn-Hilliard problem and uh the solution of the semi-discrete VEM problem. Then for all t ∈ [0, T] it holds u − uhL2(Ω) h2.

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 16 / 28

Numerical results

Exact solution: u(x, y, t) = t cos(2πx) cos(2πy) on Ω = (0, 1)2

Table: Test 1: H2, H1 and L2 errors and convergence rates α computed on four quadrilateral meshes discretizing the unit square.

h |uh − u|H2(Ω) α |uh − u|H1(Ω) α ||uh − u||L2(Ω) α 1/16 1.35e-1 – 8.57e-2 – 8.65e-2 – 1/32 5.86e-2 1.20 2.20e-2 1.96 2.20e-2 1.97 1/64 2.79e-2 1.07 5.53e-3 1.99 5.52e-3 1.99 1/128 1.38e-2 1.02 1.37e-3 2.01 1.37e-3 2.01 Parameters: γ = 1/10 and ∆t = 1e − 7.

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 17 / 28

Evolution of an ellipse (I)

(a) Quadrilateral mesh of 16384 = 128 × 128 elements (49923 dofs).

Figure: Snapshots at three temporal frames (t = 0, 0.5, 1).

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 18 / 28

Evolution of an ellipse (II)

(a) Triangular mesh of 8576 elements (13167 dofs).

Figure: Snapshots at three temporal frames (t = 0, 0.5, 1).

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 19 / 28

Voronoi polygonal meshes

Figure: Examples of Voronoi polygonal meshes (quadrilaterals, pentagons, hexagons) with 10 (left) and 100 (right) elements.

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 20 / 28

Evolution of a cross (I)

(a) Quadrilateral mesh of 16384 = 128 × 128 elements (49923 dofs).

Figure: Snapshots at three temporal frames (t = 0, 0.05, 1).

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 21 / 28

Evolution of a cross (II)

(a) Triangular mesh of 8576 elements (13167 dofs).

Figure: Snapshots at three temporal frames (t = 0, 0.05, 1).

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 22 / 28

Evolution of a cross (III)

(a) Voronoi polygonal mesh of 10000 elements (59490 dofs).

Figure: Snapshots at three temporal frames (t = 0, 0.05, 1).

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 23 / 28

Spinoidal decomposition

(a) Voronoi polygonal mesh of 10000 elements (59490 dofs).

Figure: Snapshots at three temporal frames (t = 0.01, 0.05, 5).

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 24 / 28

Conclusions

C 1-VEM of minimal degree for the approximation of CH equation. Theoretical L2-convergence of the semi-discrete scheme Numerical tests Reference

• P. F. Antonietti, L. Beirao da Veiga, S. Scacchi and M. Verani, A C 1

virtual element method for the Cahn-Hilliard equation with polygonal meshes, accepted on SINUM, 2015.

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 25 / 28

Spinoidal decomposition (I)

(a) Quadrilateral mesh of 16384 = 128 × 128 elements (49923 dofs).

Figure: Snapshots at three temporal frames (t = 0.01, 0.05, 5).

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 26 / 28

Spinoidal decomposition (II)

(a) Triangular mesh of 8576 elements (13167 dofs).

Figure: Snapshots at three temporal frames (t = 0.075, 0.25, 1.25).

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 27 / 28

VEM vs mixed-DG

DG VEM h dof |uh − u|H1(Ω) ||uh − u||L2(Ω) dof |uh − u|H1(Ω) ||uh − u||L2(Ω) 1/16 1536 7.21e-3 7.23e-5 867 2.56e-2 9.06e-4 1/32 6144 3.63e-3 1.82e-5 3267 6.70e-3 2.31e-4 1/64 24576 1.82e-3 4.57e-6 12675 1.74e-3 5.90e-5 1/128 98304 9.12e-4 1.14e-6 49923 4.46e-4 1.49e-5

CH with exact solution, γ = 1/10 and ∆t = 1e − 7. The simulation is run for 10 time steps. The H1 and L2 errors are computed at the final time step on four triangular (for the DG method) and quadrilateral (for the VEM) meshes discretizing the unit square.

• M. Verani (MOX - PoliMi)

VEM Cahn-Hilliard POEMs 2015 28 / 28