# Correct implementation of AllenCahn

Hello everyone,

I am new to FeniCS and trying to learn it for a school project by implementing the Allen-Cahn equation:

\frac{\partial \varphi}{\partial t } = \nabla^2 \varphi +\varphi -\varphi^3

Its variational formulation is as follows:

\int_{\Omega}\varphi_{n+1} v \,d\Omega = \int_{\Omega}\varphi_{n} v \,d\Omega + dt \left[-\int_{\Omega} \nabla\varphi_n \nabla v \, d\Omega \, + \int_{\Omega} (\varphi_n - \varphi_n^3)v \, d\Omega \right]

Which i implemented as a linear problem in DolfinX :

import dolfinx
from dolfinx.io.gmshio import model_to_mesh
from mpi4py import MPI
import numpy as np
import gmsh
import pyvista
import dolfinx.fem as fem
import dolfinx.mesh as mesh
import dolfinx.io as io
import dolfinx.plot as plot
import ufl
from dolfinx.fem import Expression
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
import dolfinx_mpc.utils
from dolfinx_mpc import LinearProblem, MultiPointConstraint

##Generate square mesh##
L,W = 1,1
nx,ny = 50,50
mesh = dolfinx.mesh.create_rectangle(
MPI.COMM_WORLD,
[np.array([0, 0]), np.array([L, W])],
[nx, ny], dolfinx.mesh.CellType.triangle
)

##Define scalar field
element_a= ufl.FiniteElement('CG',mesh.ufl_cell(),degree=1)
V_a = fem.FunctionSpace(mesh, element_a)
phi_next = fem.Function(V_a, name="phase")
phi_prev = fem.Function(V_a, name="phaseold")

##initial conditions : Random distribution between -1 and 1
def randm(x):
return [np.random.random()*2-1 for i in range(len(x[0]))]
phi_prev.interpolate(randm)

## Define variational problem
v = ufl.TestFunction(V_a)
u = ufl.TrialFunction(V_a)
from ufl import Measure
dx = ufl.Measure("dx", domain=mesh)
dt = 1e-6
def a(u,v):
return u*v*dx

def L(v):

##Define periodic boundary conditions along x, it need sto be along y too but corner points are wrong
def periodic_boundary_x(x):
return np.isclose(x[0], 1, atol=1e-6)

def periodic_relation_x(x):
out_x = np.zeros_like(x)
out_x[0] = x[0]-1
out_x[1] = x[1]
out_x[2] = x[2]
return out_x
mpc = MultiPointConstraint(V_a)
mpc.create_periodic_constraint_geometrical(V_a, periodic_boundary_x, periodic_relation_x,[])
mpc.finalize()

#Main solver loop
ffile = dolfinx.io.XDMFFile(MPI.COMM_WORLD, "AC_Random.xdmf", "w")
ffile.write_mesh(mesh)
ffile.write_function(phi_prev, 0)
Nstep=1000
for n in range(1,Nstep):
print("Timestep ",n)
problem = LinearProblem(a(u,v), L(v),mpc)
ph = problem.solve()
ph.name = "phase"
phi_prev.x.array[:] = ph.x.array
ffile.write_function(ph, n)

ffile.close()

The code is running and i am getting a solution but it behaves in the wrong way, as if the \varphi-\varphi³ term is not being considered. In fact, I initialize the field with a a random noise between -1 and 1 and it evolves to a smooth variation around 0, which is not the right physical behavior because the last term is supposed to penalize this and separate phases into regions of -1 others of +1

Is my implementation correct?

Thank you very much! i have update the post

Please note that you should first:

Move this outside of your temporal loop, as it does not change within the loop.

This is what I would do:

import dolfinx
from mpi4py import MPI
import numpy as np
import gmsh
import dolfinx.fem as fem
import dolfinx.mesh as mesh
import dolfinx.io as io
import dolfinx.plot as plot
import ufl
from mpi4py import MPI
from dolfinx_mpc import LinearProblem, MultiPointConstraint

## Generate square mesh##
L, W = 1, 1
nx, ny = 50, 50
mesh = dolfinx.mesh.create_rectangle(
MPI.COMM_WORLD,
[np.array([0, 0]), np.array([L, W])],
[nx, ny], dolfinx.mesh.CellType.triangle
)

# Define scalar field
element_a = ufl.FiniteElement('CG', mesh.ufl_cell(), degree=1)
V_a = fem.FunctionSpace(mesh, element_a)

# initial conditions : Random distribution between -1 and 1

def randm(x):
return [np.random.random()*2-1 for i in range(len(x[0]))]

# Define periodic boundary conditions along x, it need sto be along y too but corner points are wrong

def periodic_boundary_x(x):
return np.isclose(x[0], 1, atol=1e-6)

def periodic_relation_x(x):
out_x = np.zeros_like(x)
out_x[0] = x[0]-1
out_x[1] = x[1]
out_x[2] = x[2]
return out_x

mpc = MultiPointConstraint(V_a)
mpc.create_periodic_constraint_geometrical(
V_a, periodic_boundary_x, periodic_relation_x, [])
mpc.finalize()

# Get new MPC
phi_next = fem.Function(V_a, name="phase")
phi_prev = fem.Function(V_a, name="phaseold")
phi_prev.interpolate(randm)

# Define variational problem
v = ufl.TestFunction(V_a)
u = ufl.TrialFunction(V_a)
dx = ufl.Measure("dx", domain=mesh)
dt = 1e-6

# Main solver loop
ffile = dolfinx.io.XDMFFile(MPI.COMM_WORLD, "AC_Random.xdmf", "w")
ffile.write_mesh(mesh)
Nstep = 1000

a = u*v*dx
L = phi_prev*v*dx
L += dt * 10000 * (phi_prev-phi_prev**3)*v*dx

problem = LinearProblem(a, L, mpc)
ph = problem.u
ph.name = "phase"

ph.interpolate(phi_prev)
ffile.write_function(ph, 0)

for n in range(1, Nstep):
print("Timestep ", n)
problem.solve()
phi_prev.x.array[:] = ph.x.array
ffile.write_function(ph, n)

ffile.close()

Note that I turned up the effect of your term (as dt is very small as you are using an explicit scheme). Then I get

Reducing the magnitude of the penalization term to 1000 dfiffusion takes over.

I’m kinda new to this (phase- ) field.
Can we use this model to simulate grain growth in a microstructure?
Like this:

Taken from this paper