DG formulation of two-dimensional convection-diffusion equation

Question

Hi,

I have written the following code for the formulation of a DG problem. The discretisation of the diffusion part of the equation is the SIP method, and the convection part is discretized by upwinding. The result for the solution that I get is -nan, which is very unphysical.

Thanks for your help!

b = Expression(('gamma_in*(vbar_in-x[0])', '(beta_in-alpha_in*x[1])'), gamma_in = gamma, vbar_in =vbar, beta_in = beta, alpha_in = alpha_1)

f = Constant(-1.0)

n = FacetNormal(mesh)
h = CellSize(mesh)

of = Function(V_dg)

kappa = Constant(-1.0)

alpha = 32.0

#bilinear form B(u,v) with both SIP DG diffusion term and convection term

#diffusion part
a_convection = +dot(grad(v), dot(kappa, grad(u)))*dx + v*inner(b, grad(u))*dx \

-kappa*dot(jump(v, n), avg(grad(u)) )*dS \

-kappa*dot(avg(grad(v)), jump(u, n))*dS \

-kappa*dot(inner(v, n), grad(u))*ds \

-kappa*dot(grad(v), inner(u, n))*ds \

+kappa*alpha/h('+')*dot(jump(v, n), jump(u, n))*dS \

+kappa*alpha/h*v*u*ds

#convection part
#exterior
- dot(b, n)*(u('+'))*v('+')*ds\

#interior
- dot(b, n)*(u('+') - u('-'))*v('+')*dS\

L = f*v*dx


# Solution function
solution = Function(V_dg)

Comments

In the code above you never assemble/solve your problem. Can you provide a complete example?

---

Yes, sure. Sorry for the omission.

from dolfin import *
from numpy import *

gamma = 1.0
vbar = 9.0
beta = 10.0
alpha_1 = 1.0

epsilon = 1.0

set_log_level(10)

mesh = Mesh("D2_small_29072014.xml.gz")


# Defining the function spaces

V_dg = FunctionSpace(mesh, "DG", 1)


# Define Dirichlet boundary
def boundary(x):
    return abs(x[1] -(x[0]))< DOLFIN_EPS

u0 = Constant(0.0)
bc = DirichletBC(V_dg, u0, boundary, "geometric")


v = TestFunction(V_dg)
u = TrialFunction(V_dg)


b = Expression(('gamma_in*(vbar_in-x[0])', '(beta_in-alpha_in*x[1])'), gamma_in = gamma, vbar_in =vbar, beta_in = beta, alpha_in = alpha_1)

f = Constant(-1.0)

n = FacetNormal(mesh)
h = CellSize(mesh)

of = Function(V_dg)

kappa = Constant(-1.0)


alpha = 32.0

a_convection = +dot(grad(v), dot(kappa, grad(u)))*dx + v*inner(b, grad(u))*dx \

-kappa*dot(jump(v, n), avg(grad(u)) )*dS \

-kappa*dot(avg(grad(v)), jump(u, n))*dS \

-kappa*dot(inner(v, n), grad(u))*ds \

-kappa*dot(grad(v), inner(u, n))*ds \

+kappa*alpha/h('+')*dot(jump(v, n), jump(u, n))*dS \

+kappa*alpha/h*v*u*ds

#exterior
- dot(b, n)*(u('+'))*v('+')*ds\

#interior
- dot(b, n)*(u('+') - u('-'))*v('+')*dS\

L = f*v*dx


# Solution function
solution = Function(V_dg)

# Assemble and apply boundary conditions
A_convection = assemble(a_convection)
b_assembly = assemble(L)
bc.apply(A_convection, b_assembly)
# Solve system
solve(A_convection, solution.vector(), b_assembly)


# Plot solution
plot(solution, interactive = True)

---

For completeness, the PDE I want to solve is

$
\kappa \Delta u(x,y) + \mathbf{b}(\mathbf{x})\cdot \nabla u = -1,
$

with the choice of $\mathbf{b}$ given above.

Many thanks!

---

I am confused by the tense product.

Answer 1

Hi, first of all doing the import

from dolfin import *
from numpy import *

guarantees trouble as you override many DOLFIN objects with those from numpy, e.g. dot, inner just to name a few. Once you fix this I suggest you take a look at dg-advection-diffusion demo.

Comments

Hi, many thanks!

As far as I can see, the equation in the demo you mentioned is for the equation

$$
-\kappa \Delta u + \nabla \cdot (\mathbf{b} u) = f.
$$

Is that right? Please note that in my case, $\nabla \cdot \mathbf{b} \neq 0$, so I think I cannot use this formulation in its present form.

Thanks, Wilhelm