Application of Periodic Boundary Conditions with Unstructured Mesh

General
, ,
1

Hi everyone,

I’m using FEniCS 2019.1.0 to simulate the advection of a moving peak - ie, I’m solving the PDE

\frac{\partial u}{\partial t} + \mathbf{a} \cdot \nabla u = 0

on the domain \Omega = (0,1) \times (0,1) with periodic boundary conditions in the x-direction.

I find that if I use a structured mesh, for example

mesh = fe.UnitSquareMesh(nx, ny)

the simulation proceeds exactly as expected. However, if I instead use an unstructured mesh, say

domain = mshr.Rectangle(fe.Point(0.0, 0.0), fe.Point(1, 1))
mesh = mshr.generate_mesh(domain, numPts)

I find that if numPts is too large, spurious numerical… artifacts - appear on the left and right ends of the domain and eventually consume all of the computational domain.

I’m trying to simulate a system (much more complex than this) wherein I need periodic boundary conditions and an unstructured mesh - is there some subtlety to applying periodic boundary conditions on an unstructured mesh that I’m unaware of?

Find below a MWE:

import fenics as fe
import mshr

comm = fe.MPI.comm_world
rank = fe.MPI.rank(comm)

fe.parameters["form_compiler"]["optimize"] = True
fe.parameters["form_compiler"]["cpp_optimize"] = True

T = 100      
L_x = 1
L_y = 1
aVec = fe.Constant((1.0, 0.0))

domain = mshr.Rectangle(fe.Point(0.0, 0.0), fe.Point(L_x, L_y))
mesh = mshr.generate_mesh(domain, 80)
#mesh = fe.UnitSquareMesh(150,150)
h = mesh.hmin()

solnFile = fe.XDMFFile(comm, "unstructuredMesh.xdmf")
solnFile.parameters["flush_output"] = True
solnFile.parameters["functions_share_mesh"] = True
solnFile.parameters["rewrite_function_mesh"] = False

dt = 0.01*h
numSteps = int(T/dt)

class PeriodicBoundary(fe.SubDomain):

    def inside(self, x, on_boundary):
        return bool(
            fe.near(x[0], 0.0)
            and on_boundary
            and not fe.near(x[0], L_x)
        )

    def map(self, x, y):
        y[0] = x[0] - L_x
        y[1] = x[1]
            
pbc = PeriodicBoundary()


V = fe.FunctionSpace(mesh, "Lagrange", 1, constrained_domain=pbc)


u_trial = fe.TrialFunction(V)
v = fe.TestFunction(V)
u_n = fe.Function(V)
u_nP1 = fe.Function(V)

xc = 0.5
yc = 0.5
a = 0.25
b = 0.25

initExpr = fe.Expression(
    "max(1.0 - fabs(x[0]-xc)/a - fabs(x[1]-yc)/b, 0.0)",
    degree=1,
    xc=xc,
    yc=yc,
    a=a,
    b=b
)

u_n = fe.interpolate(initExpr, V)

bilinForm = u_trial*v*fe.dx 

linForm = u_n*v*fe.dx - dt*fe.dot(aVec,fe.grad(u_n))*v*fe.dx\
    - 0.5*dt**2 * fe.dot(aVec, fe.grad(v)) * fe.dot(aVec, fe.grad(u_n))*fe.dx


sysMat = fe.assemble(bilinForm)


t = 0.0
for n in range(numSteps):
    
    t += dt 
    
    rhsVec = fe.assemble(linForm)
    
    fe.solve(sysMat, u_nP1.vector(), rhsVec)
    
    u_n.assign(u_nP1)
    
    if n % 500 == 0:

        solnFile.write(u_n, t)