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):
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
) field.