1D upwind in DG advection-diffusion

Question

Hello,

I am trying to write a 1D code for the advection-diffusion equation using DG method.
I was looking at the demo

https://github.com/FEniCS/dolfin/blob/master/demo/undocumented/dg-advection-diffusion/python/demo_dg-advection-diffusion.py

trying to modify it.
I changed the dimensions and defined my velocity field as a constant:

from dolfin import *

# Dirichlet boundary
class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0]*(1 - x[0]), 0)

# Mesh
mesh = UnitIntervalMesh(50)

# Function spaces
V_dg = FunctionSpace(mesh, "DG", 1)
V_cg = FunctionSpace(mesh, "CG", 1)
V_u  = FunctionSpace(mesh, "CG", 2) 

# Velocity 
u_const = Constant(0.5)
u = interpolate(u_const, V_u)

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

# Diffusivity
kappa = Constant(0.0)

# Source term
f = Constant(0.0)

# Penalty term
alpha = Constant(5.0)

# Mesh-related functions
n = FacetNormal(mesh)
h = CellSize(mesh)
h_avg = (h('+') + h('-'))/2

# Upwind 
un = (u*n - abs(u*n))/2

# Mark facets of the mesh
boundaries = FacetFunction('size_t', mesh, 0)
DirichletBoundary().mark(boundaries, 1)

#Bilinear form
a_int = inner(v.dx(0), kappa*phi.dx(0) - u*phi)*dx \


a_fac = kappa*(alpha/h('+'))*inner(jump(v, n), jump(phi, n))*dS \
         - kappa*inner(avg(grad(v)), jump(phi, n))*dS \
         - kappa*inner(jump(v, n), avg(grad(phi)))*dS 


a_vel = jump(v)* (un('+')*phi('+') - un('-')*phi('-'))*dS  + (v*un*phi)*ds

When i try to run it, I have the error message

Can only integrate scalar expressions. The integrand is a tensor expression with value shape (1,) and free indices with labels ().
Traceback (most recent call last):
File "dg-advection-diffusion.py", line 62, in
a_vel = jump(v)* (un('+')phi('+') - un('-')phi('-'))dS + (vunphi)ds
File "/usr/lib/python2.7/dist-packages/ufl/measure.py", line 429, in rmul
(integrand.ufl_shape, integrand.ufl_free_indices))
File "/usr/lib/python2.7/dist-packages/ufl/log.py", line 171, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Can only integrate scalar expressions. The integrand is a tensor expression with value shape (1,) and free indices with labels ().

I think it's a problem with the upwind, but I'm not able to find it out.

What's wrong?

Thanks in advance

Answer 1

Hi,

the problem is that 'un' has rank 1 (i.e. it is a vector, since you are using the FacetNormal vector when defining 'un') whereas 'v' and 'phi' are scalar quantities.
I guess your code won't return an error when using

uv = as_vector((u,))
un = (dot(uv,n) - abs(dot(uv,n)) )/2

Now 'un' will be a scalar quantity and your form will compile (without any other guarantees of the correctness of your formulation).

Hope this helps?

Comments

Yes, thanks! Now it works