# DG for advection equation with periodic boundary condition

**URL:** <https://fenicsproject.discourse.group/t/dg-for-advection-equation-with-periodic-boundary-condition/272>\
**Category:** Uncategorized\
**Created:** [February 23, 2019, 5:24pm UTC](https://fenicsproject.discourse.group/t/dg-for-advection-equation-with-periodic-boundary-condition/272 "2019-02-23T17:24:44Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![JAEKWANG\_KIM](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/jaekwang_kim/32/5263_2.png) [@JAEKWANG\_KIM](https://fenicsproject.discourse.group/u/JAEKWANG_KIM)\
**Post date:** [February 23, 2019, 5:24pm UTC](https://fenicsproject.discourse.group/t/dg-for-advection-equation-with-periodic-boundary-condition/272/1 "2019-02-23T17:24:44Z")

</div>

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(2_pi_x[0])”,degree=2, domain=mesh)

Dbc\_values = Expression((“sin(2_pi_x[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 - dt_dot(grad(v),b\*phi)_dx + dt_a1  
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

---

<div class="post-metadata">

**Author:** ![tormod](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/tormod/32/27_2.png) [@tormod](https://fenicsproject.discourse.group/u/tormod)\
**Post date:** [February 24, 2019, 6:03pm UTC](https://fenicsproject.discourse.group/t/dg-for-advection-equation-with-periodic-boundary-condition/272/2 "2019-02-24T18:03:17Z")

</div>

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