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