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)
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)

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.211035963088022 {'converged': True, 'iterations': 2, 'linear_iterations': 2, 'timing': [0.9910155230000001, 0.608492848, 0.38252267500000003]}
[ 2 ] 0.20000000000000004 0.01 0.19941969737503615 {'converged': True, 'iterations': 3, 'linear_iterations': 3, 'timing': [1.919479224, 1.339368585, 0.5801106390000002]}
[ 3 ] 0.3000000000000001 0.01 0.1519815460612822 {'converged': True, 'iterations': 3, 'linear_iterations': 3, 'timing': [1.629681892, 1.133977905, 0.49570398699999996]}
[ 4 ] 0.4000000000000002 0.01 0.1325105395823258 {'converged': True, 'iterations': 3, 'linear_iterations': 3, 'timing': [1.303063359, 0.897980981, 0.405082378]}
[ 5 ] 0.5000000000000002 0.01 0.11181937957615219 {'converged': True, 'iterations': 3, 'linear_iterations': 3, 'timing': [0.81255796, 0.5387415050000001, 0.27381645499999996]}
[ 6 ] 0.6000000000000003 0.01 0.10911562145598705 {'converged': True, 'iterations': 2, 'linear_iterations': 2, 'timing': [0.49923511600000003, 0.30994939099999996, 0.18928572500000007]}
[ 7 ] 0.7000000000000004 0.01 0.10660773284430594 {'converged': True, 'iterations': 2, 'linear_iterations': 2, 'timing': [0.49076930900000004, 0.30500042800000005, 0.185768881]}
[ 8 ] 0.8000000000000005 0.01 0.1039948686878176 {'converged': True, 'iterations': 2, 'linear_iterations': 2, 'timing': [0.532539046, 0.31969493800000004, 0.21284410799999998]}

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.09160398897424067 {'converged': True, 'iterations': 3, 'linear_iterations': 3, 'timing': [0.75154695, 0.457609778, 0.293937172]}
[ 2 ] 1.6000000000000012 0.04 0.0700641504662594 {'converged': True, 'iterations': 3, 'linear_iterations': 3, 'timing': [0.722812607, 0.45562290000000005, 0.267189707]}
[ 3 ] 2.0000000000000013 0.04 0.058822337844102145 {'converged': True, 'iterations': 2, 'linear_iterations': 2, 'timing': [0.541193616, 0.30652525900000005, 0.234668357]}

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