# Correct implementation of AllenCahn

**URL:** <https://fenicsproject.discourse.group/t/correct-implementation-of-allencahn/13487>\
**Category:** variational formulation\
**Tags:** dolfinx\
**Created:** [February 2, 2024, 10:04pm UTC](https://fenicsproject.discourse.group/t/correct-implementation-of-allencahn/13487 "2024-02-02T22:04:30Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![AymaneGr](https://avatars.discourse-cdn.com/v4/letter/a/ba9def/32.png) [@AymaneGr](https://fenicsproject.discourse.group/u/AymaneGr)\
**Post date:** [February 2, 2024, 10:04pm UTC](https://fenicsproject.discourse.group/t/correct-implementation-of-allencahn/13487/1 "2024-02-02T22:04:31Z")

</div>

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 :

```auto
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):
    return phi_prev*v*dx-dt* ufl.dot(ufl.grad(phi_prev),ufl.grad(v))*dx - dt*phi_prev*v*dx + dt*phi_prev**3*v*dx

##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 for your help

---

<div class="post-metadata">

**Author:** ![francesco-ballarin](https://avatars.discourse-cdn.com/v4/letter/f/f05b48/32.png) [@francesco-ballarin](https://fenicsproject.discourse.group/u/francesco-ballarin)\
**Post date:** [February 3, 2024, 5:33am UTC](https://fenicsproject.discourse.group/t/correct-implementation-of-allencahn/13487/2 "2024-02-03T05:33:12Z")

</div>

Please read [Read before posting: How do I get my question answered?](https://fenicsproject.discourse.group/t/read-before-posting-how-do-i-get-my-question-answered/21). We cannot comment unless you post the full code.

---

<div class="post-metadata">

**Author:** ![AymaneGr](https://avatars.discourse-cdn.com/v4/letter/a/ba9def/32.png) [@AymaneGr](https://fenicsproject.discourse.group/u/AymaneGr)\
**Post date:** [February 3, 2024, 8:26am UTC](https://fenicsproject.discourse.group/t/correct-implementation-of-allencahn/13487/3 "2024-02-03T08:26:15Z")

</div>

Thank you very much! i have update the post

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [February 6, 2024, 8:15pm UTC](https://fenicsproject.discourse.group/t/correct-implementation-of-allencahn/13487/4 "2024-02-06T20:15:55Z")

</div>

Please note that you should first:

> [@AymaneGr](#):
>
> ` problem = LinearProblem(a(u,v), L(v),mpc)`

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

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [March 20, 2024, 8:23pm UTC](https://fenicsproject.discourse.group/t/correct-implementation-of-allencahn/13487/5 "2024-03-20T20:23:02Z")

</div>

This is what I would do:

```python
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 * ufl.dot(ufl.grad(phi_prev), ufl.grad(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

 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/5/53cd1a39c232850ee77f1fbe581d8624562d61e2.jpeg)

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

---

<div class="post-metadata">

**Author:** ![RemDelaporteMathurin](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/remdelaportemathurin/32/4543_2.png) [@RemDelaporteMathurin](https://fenicsproject.discourse.group/u/RemDelaporteMathurin)\
**Post date:** [April 10, 2024, 2:44pm UTC](https://fenicsproject.discourse.group/t/correct-implementation-of-allencahn/13487/6 "2024-04-10T14:44:41Z")

</div>

Hi @AymaneGr @dokken

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

 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/1/11d4d0f5e6b4edb37f3238939f04595b8967d26c.jpeg)  
Taken from [this paper](https://doi.org/10.1016/j.ijhydene.2023.11.270)