Advection in periodic domain using DG

Question

This is related to this and this.

Has the problem of periodic boundary conditions and DG been solved?

With version 1.3 of Fenics, the code below compiles and works. However as soon as the function rho reaches the end of the domain, it leaves the boundary and does not "wrap around" as expected for periodic boundary conditions (see images below)

Am I missing something simple here? Or is the problem still around?

from dolfin import *

cpp.set_log_level(cpp.WARNING)

class PBC1D(SubDomain):
    def __init__(self, xmin=0, xmax=1):
        SubDomain.__init__(self)
        assert(xmax>xmin)
        self.xmin=xmin
        self.xmax=xmax
    # Left boundary is "target domain" G
    def inside(self, r, on_boundary):
        # return True if on left boundary
        return bool(near(r[0], self.xmin) and on_boundary)
    def map(self, r, s):
        s[0] = r[0] - self.xmax + self.xmin

mesh = IntervalMesh(100, 0.0, 1.0)
pbc = PBC1D(0.0, 1.0)

# Define function spaces and basis functions
V_dg = FunctionSpace(mesh, "DG", 1, constrained_domain=pbc)
Mv = VectorFunctionSpace(mesh, "CG", 1, constrained_domain=pbc)

# advecting velocity
v_hat = Expression(("-1.0",))
u = interpolate(v_hat, Mv)

# Test and trial functions
phi = TestFunction(V_dg)
rho = TrialFunction(V_dg)

# Mesh-related functions
n = FacetNormal(mesh)

# ( dot(v, n) + |dot(v, n)| )/2.0
un = (dot(u, n) + abs(dot(u, n)))/2.0

# parameters
dt = 0.001

# Bilinear form
a_mass = phi*rho*dx
a_int = dot(grad(phi), -u*rho)*dx
a_flux = (dot(jump(phi), un('+')*rho('+') - un('-')*rho('-') )*dS  
          + dot(phi, un*rho)*ds)

M = assemble(a_mass)
arhs = -dt*(a_int + a_flux)

rho0 = Expression("exp(-pow((x[0]-0.5),2)/0.01)")
rho  = interpolate(rho0, V_dg)

drho1 = Function(V_dg)
rho1  = Function(V_dg)
drho2 = Function(V_dg)
rho2  = Function(V_dg)
drho3 = Function(V_dg)


t = 0.0
T = 1.0

plot(rho, rescale=False)
while(t < (T-dt/2)):
    L = assemble(action(arhs,rho))
    solve(M,drho1.vector(),L)

    rho1.vector()[:] = rho.vector().copy()
    rho1.vector().axpy(1.0, drho1.vector())
    L = assemble(action(arhs,rho1))
    solve(M,drho2.vector(),L)

    rho2.vector()[:] = rho.vector().copy()
    rho2.vector().axpy(0.25, drho1.vector())
    rho2.vector().axpy(0.25, drho2.vector())
    L = assemble(action(arhs,rho2))
    solve(M,drho3.vector(),L)

    rho.vector().axpy((1.0/6.0), drho1.vector())
    rho.vector().axpy((1.0/6.0), drho2.vector())
    rho.vector().axpy((2.0/3.0), drho3.vector())
    t +=dt
    plot(rho)

interactive()

Comments

Any update on this....?

I want to know whether it is my mistake or this just does not work in FEniCS yet. Would really appreciate any help.

---

This was 1.5 years ago, but any news on the topic?

Answer 1

This is of interest to me also. Have you tried to implement the periodic bc through the fluxes. Mark the left and right boundaries as "0" and "1". Your advection speed = -1. So at right boundary you need a flux of the form

phi*u*rho(0)*ds(1)

I may be a bit wrong with the sign of the terms but the idea should work.

rho(0) gives the value of rho at x=0. I think if you assemble the rhs like this the periodic bc should work.

Comments

Thanks for the suggestion. I tried to implement what you said. But it does not seem to work. Changing the sign in front of the additional flux term also does not seem to change anything.

from dolfin import *

cpp.set_log_level(cpp.WARNING)

class PBC1D(SubDomain):
    def __init__(self, xmin=0, xmax=1):
        SubDomain.__init__(self)
        assert(xmax>xmin)
        self.xmin=xmin
        self.xmax=xmax
    # Left boundary is "target domain" G
    def inside(self, r, on_boundary):
        # return True if on left boundary
        return bool(near(r[0], self.xmin) and on_boundary)
    def map(self, r, s):
        s[0] = r[0] - self.xmax + self.xmin

mesh = IntervalMesh(100, 0.0, 1.0)
pbc = PBC1D(0.0, 1.0)

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)

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)

# Define function spaces and basis functions
V_dg = FunctionSpace(mesh, "DG", 1, constrained_domain=pbc)
Mv = VectorFunctionSpace(mesh, "CG", 1, constrained_domain=pbc)

# advecting velocity
v_hat = Expression(("-1.0",))
u = interpolate(v_hat, Mv)

# Test and trial functions
phi = TestFunction(V_dg)
rho = TrialFunction(V_dg)

# Mesh-related functions
n = FacetNormal(mesh)

# ( dot(v, n) + |dot(v, n)| )/2.0
un = (dot(u, n) + abs(dot(u, n)))/2.0

# parameters
dt = 0.001

# Bilinear form
a_mass = phi*rho*dx
a_int = dot(grad(phi), -u*rho)*dx
a_flux = (dot(jump(phi), un('+')*rho('+') - un('-')*rho('-') )*dS  
          + dot(phi, un*rho)*ds - phi*un*rho(0)*ds(1))

M = assemble(a_mass)
arhs = -dt*(a_int + a_flux)

rho0 = Expression("exp(-pow((x[0]-0.5),2)/0.01)")
rho  = interpolate(rho0, V_dg)

drho1 = Function(V_dg)
rho1  = Function(V_dg)
drho2 = Function(V_dg)
rho2  = Function(V_dg)
drho3 = Function(V_dg)


t = 0.0
T = 1.0

plot(rho, rescale=False)
while(t < (T-dt/2)):
    L = assemble(action(arhs,rho))
    solve(M,drho1.vector(),L)

    rho1.vector()[:] = rho.vector().copy()
    rho1.vector().axpy(1.0, drho1.vector())
    L = assemble(action(arhs,rho1))
    solve(M,drho2.vector(),L)

    rho2.vector()[:] = rho.vector().copy()
    rho2.vector().axpy(0.25, drho1.vector())
    rho2.vector().axpy(0.25, drho2.vector())
    L = assemble(action(arhs,rho2))
    solve(M,drho3.vector(),L)

    rho.vector().axpy((1.0/6.0), drho1.vector())
    rho.vector().axpy((1.0/6.0), drho2.vector())
    rho.vector().axpy((2.0/3.0), drho3.vector())
    t +=dt
    plot(rho)

interactive()