DOLFIN

Periodic BC for 1-d unsteady inviscid burgers equation

Asked by Praveen C

Hello all

I wrote a DG code for 1-d burgers equation and I want to use periodic BC. The code below is doing forward euler time stepping. I cant figure out where to apply periodic BC.

Thanks
praveen

# Sub domain for Dirichlet boundary condition
class Boundary(SubDomain):
   def inside(self, x, on_boundary):
      return (x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS) and on_boundary

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

# Initial condition
u0 = Expression("sin(2*pi*x[0])")

# Load mesh
mesh = UnitInterval(100)

# Defining the function spaces
V = FunctionSpace(mesh, "DG", 1)

# Periodic boundary condition
bc = PeriodicBC(V, Boundary())

# Test and trial functions
phi = TestFunction(V)
uh = TrialFunction(V)

uold = Function(V) # previous solution
u = Function(V) # current solution

# Set initial condition
u = project(u0, V)

# Mesh-related functions
n = FacetNormal(mesh)
un = dot(u, n)

# Numerical flux function - Roe flux
Hn = 0.25*(u('+')**2 + u('-')**2) - 0.5*abs(un('+')+un('-'))*(u('-') - u('+'))

# Mass matrix
M = uh*phi*dx
MMat = assemble(M)

dt = 0.001 # time step
T = 0.05 # final time
t = 0.0 # time counter

# RHS
R = uold*phi*dx + dt*0.5*u**2*phi.dx(0)*dx - dt*Hn*jump(phi)*dS

ufile = File("u.pvd")
ufile << u

while t < T + DOLFIN_EPS:
   uold.assign(u)
   rhs = assemble(R)
   solve(MMat, u.vector(), rhs)
   t += dt
   print t
   ufile << u

Question information

Language:
English Edit question
Status:
Expired
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Praveen C (cpraveen) said : 		#1

I have to add the boundary terms also

R = uold*phi*dx + dt*0.5*u**2*phi.dx(0)*dx - dt*Hn*jump(phi)*dS - dt*Hn*phi*ds

For the flux Hn at left boundary, u('-') is available but u('+') has to come from periodicity, i.e., it comes from the last cell which is adjacent to right boundary. I dont know how to put this in the form.

Revision history for this message
Garth Wells (garth-wells) said : 		#2

It's not clear what your question is. If you send code, it needs to be shorter.

Revision history for this message
Praveen C (cpraveen) said : 		#3

I have tried to simplify the problem. The function u is periodic on [0,1]. I need to evaluate a form R which consists of interior face and boundary face terms. This contains a numerical flux function Hn evaluated at every face, which takes u('+') and u('-') as arguments. But at x=0 and x=1, I dont know how to define Hn, since I have to use the fact that u is periodic. Is it possible to specify that u is periodic, so that the definition of Hn below will use that information to evaluate the ds terms ?

# Load mesh
mesh = UnitInterval(100)

# Defining the function spaces
V = FunctionSpace(mesh, "DG", 1)

# Test and trial functions
phi = TestFunction(V)
u = Function(V) # current solution

u0 = Expression("sin(2*pi*x[0])")

# Set initial condition
u = project(u0, V)

# Mesh-related functions
n = FacetNormal(mesh)
un = dot(u, n)

# Numerical flux function - Roe flux
Hn = 0.25*(un('+')**2 + un('-')**2) - 0.5*abs(un('+')+un('-'))*(u('-') - u('+'))

dt = 0.005 # time step

# RHS
R = - dt*Hn*jump(phi)*dS - dt*Hn*phi*ds

rhs = assemble(R)

Revision history for this message
Launchpad Janitor (janitor) said : 		#4

This question was expired because it remained in the 'Open' state without activity for the last 15 days.