Periodic Boundary Conditions for Allen-Cahn

variational formulation
1

Howdy forum,
To preface, I’m still pretty new to FEniCS and finite elements, so I appreciate the patience.

I’m working on an Allen-Cahn Equation that has periodic boundary conditions that’s presented as follows:

u_t -0.0001u_{xx} + 4(u^3-u) = 0, x\in[-1,1], t\in[0,1]
u(0,x)=x^2cos(2\pi x),
u(t,-1)=u(t,1),
u_x(t,-1)=u_x(t,1)

I’ve looked through this example in the documentation, but for whatever reason, my solution to this particular problem isn’t periodic at the boundaries.
tmp

Here’s a minimal working example of my code:

from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
%matplotlib notebook

# parameters
#####################################
nx       = 1000        # mesh points
theta    = 0.5         # time discretization
dt       = 5.0e-3      # time step
(x0, xf) = (-1.0, 1.0) # boundaries
order    = 2           # mesh polynomial order

# Class for interfacing with the Newton solver
class AllenCahnEquation(NonlinearProblem):
    def __init__(self, a, L):
        NonlinearProblem.__init__(self)
        self.L = L
        self.a = a
    def F(self, b, x):
        assemble(self.L, tensor=b)
    def J(self, A, x):
        assemble(self.a, tensor=A)
        
# Sub domain for Periodic boundary condition
class PeriodicBoundary(SubDomain):

    # Left boundary is "target domain" G
    def inside(self, x, on_boundary):
        return bool(x[0] < DOLFIN_EPS and x[0] > -DOLFIN_EPS and on_boundary)

    # Map right boundary (H) to left boundary (G)
    def map(self, x, y):
        y[0] = x[0] - (xf - x0)

# formulate problem
#####################################
# create periodic boundary condition
pbc = PeriodicBoundary()

# setup mesh
mesh = IntervalMesh(nx, x0, xf)
V = FunctionSpace(mesh, "CG", order, constrained_domain=pbc)

# define test & trial functions
du = TrialFunction(V)
v = TestFunction(V)

# define functions
u = Function(V)
u0 = Function(V)

# initial conditions
u_init = Expression("pow(x[0],2)*sin(2*pi*x[0])", degree=2)
u.interpolate(u_init)
u0.interpolate(u_init)

mu_mid = (1.0-theta)*u0 + theta*u
gamma1 = 0.0001
gamma2 = 4.0 # {1 (easiest):4 (hardest)}

F = u*v*dx - u0*v*dx +                    \
    dt*gamma1*dot(grad(u), grad(v))*dx +  \
    dt*gamma2*(u**3 - u)*v*dx

J = derivative(F, u, du)

problem = AllenCahnEquation(J, F)
solver = NewtonSolver()

# map mesh vertices to solution DOFs
#####################################
dof_coordinates = V.tabulate_dof_coordinates()

# get indicies of sorted result
dofs      = np.squeeze(dof_coordinates)
asc_order = np.argsort(dofs)

# time stepping
#####################################
(t, T) = (0.0, 1.0)

plt.figure
plt.title(f"Allen-Cahn: $\gamma_2={gamma2}$")
plt.ylabel("$u(t)$")
plt.xlabel("$x$")
labels = []

while t < T:
    
    # compute current solution
    solver.solve(problem, u.vector())
    
    # update previous solution
    u0.vector()[:] = u.vector()
    
    # plot 6 solution snapshots
    if round(t/dt) % round(T/dt/6) == 0:
        plt.plot(dofs[asc_order], u.vector()[asc_order])
        labels.append(f"t = {t/T:.2f}")
        
    # increment time
    t += dt
    n += 1
    
plt.legend(labels)
plt.show()

From my understanding, the map() function should be enforcing that any available solutions must have equal values on both boundaries, which doesn’t seem to be the case based on my results.

Also, how would I go about enforcing equal derivative values at the boundaries, as it states in the problem? Is this possible?

Thanks.

2

After looking into the example a bit more, I realized that I was defining my inside() function incorrectly. It should’ve been:

def inside(self, x, on_boundary):
    return near(x[0], x0) and on_boundary

I’m not extremely clear on why we should only be targetting the left boundary as both sides should be periodic, but the problem seems to have been fixed.