DG for advection equation with periodic boundary condition

Hi, all

I am a new to Fenics and I am learning this amazing tools these days.

I have a question on implementation of DG with periodic boundary condition.

The most simplified problem (minimal example) that describes my interest is

At domain x in [0,1] Solve for

(dP)/(dt) + u \codt (dP)/(dx) = 0

with given and fixed u=(1,0) and initial condition of P

P(x,t=0)=sin(2 pi x)

If I implement periodic boundary condition P(x=0)=P(x=1) correctly,

my solution will be the traveling sin wave.

The scratch of code is attached…

I marked one line with red on my understanding is not clear.

The main problem is that my solution at each time dissipates to zero

both P.min() and P.max()…which was initially -1 and 1 respectively.

Can anyone suggest any idea on what point I am missing here?..

I will be very appreciated that.


Best regards,

Jaekwang Kim

Code lines

import time
import os
import math
from dolfin import *
import matplotlib.pyplot as plt
from matplotlib import animation
import subprocess

t_end = 3.00
dt = 0.05

Boundary conditions

class DirichletBoundary(SubDomain):
def inside(self, x, on_boundary):
return bool((x[1] < DOLFIN_EPS or x[1] > (1.0 - DOLFIN_EPS))
and on_boundary)

##To create periodic boundary
class PeriodicBoundary(SubDomain):

def inside(self, x, on_boundary):
return bool(x[0] < DOLFIN_EPS and x[0] > -DOLFIN_EPS and on_boundary)

def map(self, x, y):
y[0] = x[0] - 1.0
y[1] = x[1]

pbc = PeriodicBoundary()

implement subdomains for periodic BC

class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0)

class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 1.0)

Create mesh and define function space

mesh = RectangleMesh(Point(0,0),Point(1.0,1.0),40,40,“right/left”)
V = FunctionSpace(mesh, “DG”, 1,constrained_domain=pbc)
ic= Expression(“sin(2pix[0])”,degree=2, domain=mesh)

Dbc_values = Expression((“sin(2pix[0])”),degree=2, domain=mesh)
dbc = DirichletBoundary()

sub_domains = MeshFunction(“size_t”, mesh, mesh.topology().dim() - 1)
left = Left()
left.mark(sub_domains, 0)
right = Right()
right.mark(sub_domains, 1)

Given Flow field b=(1,0)

b = Expression((“1”,“0”),degree=2, domain=mesh)

Define unknown and test function(s)

v = TestFunction(V)
phi = TrialFunction(V)

phi_c = Function(V) #current_solution
phi_0 = Function(V) #old_time_solution
phi_0 = interpolate(ic,V)

n = FacetNormal(mesh)
bn = (dot(b,n)+abs(dot(b,n)))/2.0

Define (Internal and boundary) Face Integral

def a(phi,v) :
a_internal_f = dot(jump(v), bn(’+’)*phi(’+’) - bn(’-’)*phi(’-’) )dS
a_external_f = dot(v,bn
return a_internal_f+a_external_f

Define variational forms

a1=a(phi,v) ##implicit

Implicit method

A_form = inner(v,phi)dx - dtdot(grad(v),b*phi)dx + dta1
A = assemble(A_form)

For periodic boundary condition…I am not sure on this

   b_form = inner(v,phi_0)*dx +  dt * inner(v,bn*phi_0)*ds(0) 

bvector = assemble(b_form)



while t < t_end:

print(‘t=’ “%.2f” % t)
solve(A,phi_c.vector(), bvector, ‘lu’,‘none’)

bvector = assemble(b_form)

fname = ‘tmp/tmp%03d.png’ % n

print(‘phi max:’,phi_c.vector().get_local().max())
print(‘phi min:’,phi_c.vector().get_local().min())

vtkfile = File(‘final_solution.pvd’)
vtkfile << phi_c

FEniCS does not currently support periodic BCs for DG function spaces as far as I know. There was talk about redoing the periodic BC support to be an integral part of the mesh definition, which would make it work, but I do not think that is done yet. If it is implemented it is only in DOLFIN-X, the experimental version where many other things are missing, so it is probably not quite ready for normal users yet.

Your options are:

  • Implement the periodic facet integrals “manually”. Beware that you probably must extend the non-zero locations in the sparse matrix to include the coupling terms you need
  • Fix DOLFIN so that periodic BCs work for DG function spaces. Potentially easier than the first option, I have not really looked into what needs to be done
  • Do not use DG + periodic BCs