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]
[1]:
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:
[2]:
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:
[3]:
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
[4]:
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:
[5]:
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:
[6]:
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:
[7]:
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]}
