# 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.

Thanks.

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=0
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()
bc=DirichletBC(V,Dbc_values,dbc)

sub_domains = MeshFunction(“size_t”, mesh, mesh.topology().dim() - 1)
sub_domains.set_all(2)
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
phi)*ds
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)
bc.apply(A)

## 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)
bc.apply(bvector)

## Time-stepping

n=0

while t < t_end:

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

bvector = assemble(b_form)
bc.apply(bvector)

plot(phi_c)
fname = ‘tmp/tmp%03d.png’ % n
plt.savefig(fname)

print(‘phi max:’,phi_c.vector().get_local().max())
print(‘phi min:’,phi_c.vector().get_local().min())
n+=1

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.