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,bnphi)*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