Cahn-Hilliard

In this script we show how to use the \(C^1\) non-conforming virtual element space to solve the Cahn-Hilliard equation. We use a fully implicit scheme here. You can find more details in [DH21]

try:
    import dune.vem
except:
    print("This example needs 'dune.vem'")
    import sys
    sys.exit(0)
import matplotlib
matplotlib.rc( 'image', cmap='jet' )
from matplotlib import pyplot
import random
from dune.grid import cartesianDomain, gridFunction
import dune.fem
from dune.fem.plotting import plotPointData as plot
from dune.fem.function import integrate, discreteFunction

from ufl import *
import dune.ufl

dune.fem.threading.use = 4

Grid and space construction - we use a cube grid here:

order        = 3
polyGrid = dune.vem.polyGrid( cartesianDomain([0,0],[1,1],[30,30]), cubes=True)
ncC1testSpaces = [ [0], [order-3,order-2], [order-4] ]
space = dune.vem.vemSpace(polyGrid, order=order, testSpaces=ncC1testSpaces)

order=3 using 3: [10,6,3,0,0]   edge: [3,2,2] 2 1 max size of edge set: 3
max stabilization factor: 0.68232

To define the mathematical model, let \(\psi\colon{\mathbb R} \rightarrow \mathbb{R}\) be defined as \(\psi(x) = \frac{(1-x^2)^2}{4}\) and let \(\phi(x) = \psi(x)^{\prime}\). The strong form for the solution \(u\colon \Omega \times [0,T] \rightarrow {\mathbb R}\) is given by \begin{align*} \partial_t u - \Delta (\phi(u)-\epsilon^2 \Delta u) = 0 \quad &\text{in} \ \Omega \times [0,T] ,\\ u(\cdot,0) = u_0(\cdot) \quad &\text{in} \ \Omega,\\ \partial_n u = \partial_n \big( \phi(u) - \epsilon^2\Delta u \big) = 0 \quad &\text{on} \ \partial \Omega \times [0,T]. \end{align*}

We use a backward Euler discretization in time and will fix the constant further down:

t     = dune.ufl.Constant(0,"time")
tau   = dune.ufl.Constant(0,"dt")
eps   = dune.ufl.Constant(0,"eps")
df_n  = discreteFunction(space, name="oldSolution") # previous solution
x     = SpatialCoordinate(space)
u     = TrialFunction(space)
v     = TestFunction(space)

H     = lambda v: grad(grad(v))
laplace = lambda v: H(v)[0,0]+H(v)[1,1]
a     = lambda u,v: inner(H(u),H(v))
b     = lambda u,v: inner( grad(u),grad(v) )
W     = lambda v: 1/4*(v**2-1)**2
dW    = lambda v: (v**2-1)*v

equation = ( u*v + tau*eps*eps*a(u,v) + tau*b(dW(u),v) ) * dx == df_n*v * dx

dbc = [dune.ufl.DirichletBC(space, 0, i+1) for i in range(4)]

# energy
Eint  = lambda v: eps*eps/2*inner(grad(v),grad(v))+W(v)

Next we construct the scheme providing some suitable expressions to stabilize the method

biLaplaceCoeff = eps*eps*tau
diffCoeff      = 2*tau
massCoeff      = 1
scheme = dune.vem.vemScheme(
                       [equation, *dbc],
                       solver=("suitesparse","umfpack"),
                       hessStabilization=biLaplaceCoeff,
                       gradStabilization=diffCoeff,
                       massStabilization=massCoeff,
                       boundary="derivative") # only fix the normal derivative = 0

To avoid problems with over- and undershoots we project the initial conditions into a linear lagrange space before interpolating into the VEM space:

def initial(x):
    h = 0.01
    g0 = lambda x,x0,T: conditional(x-x0<-T/2,0,conditional(x-x0>T/2,0,sin(2*pi/T*(x-x0))**3))
    G  = lambda x,y,x0,y0,T: g0(x,x0,T)*g0(y,y0,T)
    return 0.5*G(x[0],x[1],0.5,0.5,50*h)
initial = dune.fem.space.lagrange(polyGrid,order=1).interpolate(initial(x),name="initial")
df = space.interpolate(initial, name="solution")

Finally the time loop:

t.value = 0
eps.value = 0.05
tau.value = 1e-02

fig = pyplot.figure(figsize=(40,20))
count = 0
pos = 1
while t.value < 0.8:
    df_n.assign(df)
    info = scheme.solve(target=df)
    t.value += tau
    count += 1
    if count % 10 == 0:
        df.plot(figure=(fig,240+pos),colorbar=None,clim=[-1,1])
        energy = integrate(polyGrid,Eint(df),order=5)
        print("[",pos,"]",t.value,tau.value,energy,info,flush=True)
        pos += 1

[ 1 ] 0.09999999999999999 0.01 0.21103596308802303 {'converged': True, 'iterations': 2, 'linear_iterations': 2}
[ 2 ] 0.20000000000000004 0.01 0.1994196973749465 {'converged': True, 'iterations': 3, 'linear_iterations': 3}
[ 3 ] 0.3000000000000001 0.01 0.15198154606026518 {'converged': True, 'iterations': 3, 'linear_iterations': 3}
[ 4 ] 0.4000000000000002 0.01 0.13251053958143946 {'converged': True, 'iterations': 3, 'linear_iterations': 3}
[ 5 ] 0.5000000000000002 0.01 0.11181937957616332 {'converged': True, 'iterations': 3, 'linear_iterations': 3}
[ 6 ] 0.6000000000000003 0.01 0.10911562145620904 {'converged': True, 'iterations': 2, 'linear_iterations': 2}
[ 7 ] 0.7000000000000004 0.01 0.10660773284453659 {'converged': True, 'iterations': 2, 'linear_iterations': 2}
[ 8 ] 0.8000000000000005 0.01 0.10399486868805612 {'converged': True, 'iterations': 2, 'linear_iterations': 2}

For the final coarsening phase we increase the time step a bit:

tau.value = 4e-02
fig = pyplot.figure(figsize=(30,10))
count = 0
pos = 1
while t.value < 2.0:
    df_n.assign(df)
    info = scheme.solve(target=df)
    t.value += tau
    count += 1
    if count % 10 == 0:
        df.plot(figure=(fig,140+pos),colorbar=None,clim=[-1,1])
        energy = integrate(polyGrid,Eint(df),order=5)
        print("[",pos,"]",t.value,tau.value,energy,info,flush=True)
        pos += 1

[ 1 ] 1.2000000000000008 0.04 0.09160398897458233 {'converged': True, 'iterations': 3, 'linear_iterations': 3}
[ 2 ] 1.6000000000000012 0.04 0.07006415046565458 {'converged': True, 'iterations': 3, 'linear_iterations': 3}
[ 3 ] 2.0000000000000013 0.04 0.05882233784407939 {'converged': True, 'iterations': 2, 'linear_iterations': 2}

This page was generated from the notebook chimpl_nb.ipynb and is part of the tutorial for the dune-fem python bindings