Free outflow in DG with dolfin_dg

1

I am using the dolfin_dg library (https://bitbucket.org/nate-sime/dolfin_dg/src/master/) to solve the convection diffusion equation

\mathbf{b} \cdot \nabla t - \kappa \nabla^2 t = 0

The domain is a simple unit square. The advection field is a horizontal vector \mathbf{b} = (1.0). There are two Dirichlet BC in the lower left side of the domain. Zero outflow in the top and bottom and I’d like to have free outflow in the right hand side. What I imagine is having very little diffusion and the advection field transporting the Dirichlet BC from the left to the right uniformly. What happens now though when I impose zero flow everywhere is that there is some sort of accumulation of the field t in the upper side instead of leaving freely. This is because I am marking all the boundary edges in

CompiledSubdomain("on_boundary").mark(inlets_mshfunc, OUTFLOW)

If I change it to

CompiledSubdomain("on_boundary && !(x[0] > (1.0 - DOLFIN_EPS)").mark(inlets_mshfunc, OUTFLOW)

The right side is not marked, but the problem does not converge. I think I might be missing some fundamental concept on how to pose boundary conditions on hyperbolic problems

Here is the code

from dolfin import *
from dolfin_dg import *
import numpy as np

__author__ = 'njcs4'

parameters['form_compiler']["cpp_optimize"] = True
parameters['form_compiler']["optimize"] = True
parameters['form_compiler']['representation'] = 'uflacs'
parameters["ghost_mode"] = "shared_facet"

ele_n = 100
p = 1

mesh = UnitSquareMesh(ele_n, ele_n, 'left/right')

V = FunctionSpace(mesh, 'DG', p)
v = TestFunction(V)

class Inlet1(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] < DOLFIN_EPS and between(x[1], (0.2, 0.4))

class Inlet2(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] < DOLFIN_EPS and between(x[1], (0.0, 0.2))

inlets_mshfunc = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
inlets_mshfunc.set_all(0)

INLET1, INLET2, OUTFLOW = 1, 3, 30

CompiledSubDomain("on_boundary && !(x[0] > (1.0 - DOLFIN_EPS))").mark(inlets_mshfunc, OUTFLOW)

inlets = Inlet1()
inlets.mark(inlets_mshfunc, INLET1)
inlets = Inlet2()
inlets.mark(inlets_mshfunc, INLET2)

ds = Measure('ds', domain=mesh, subdomain_data=inlets_mshfunc)

gD1 = Constant(100)
gD2 = Constant(1)
gN = Constant(0.0)

u = Function(V)
b = Constant((10, 0))
f = Constant(0.0)
n = FacetNormal(mesh)

def F_c(u):
    return b*u

H = LocalLaxFriedrichs(lambda u, n: dot(b, n))

kappa = 1e-8
def F_v(u, grad_u):
    return kappa*grad_u

bcs = [DGNeumannBC(ds(OUTFLOW), gN),
        DGDirichletBC(ds(INLET1), gD1),
        DGDirichletBC(ds(INLET2), gD2)]
ho = HyperbolicOperator(mesh, V, bcs, F_c, H)
eo = EllipticOperator(mesh, V, bcs, F_v)

residual = ho.generate_fem_formulation(u, v) \
+ eo.generate_fem_formulation(u, v) \
- f*v*dx

du = TrialFunction(V)
J = derivative(residual, u, du)
solve(residual == 0, u, [], J=J, solver_parameters={"newton_solver" : {"absolute_tolerance": 1e-6}})
File("dolfin_dg_adv_diff.pvd") << u
2

So I think the issue you’re having is that your problem is ilposed when kappa=0.0.

You have a region of the mesh, x[1] > 0.4 with no numerical flux because no boundary condition is prescribed for information to “flow” through the mesh. So the system is singular.

It will work when kappa > 0 because numerical flux is allowed to “diffuse” into the x[1] > 0.4 region.

You should think carefully about your boundary conditions relative to the physical problem you want to solve.