JIT error and dot() dimensions

Question

I am writing a DG discretization of advection-diffusion-reaction, and have two outstanding problems (reproducing (3.4) in http://arturo.imati.cnr.it/marini/reports/final_sinum.pdf). One is that dolfin.dot() says the dimensions are wrong in the following:

dolfin.dot(dolfin.jump(phi), dolfin.avg(dolfin.grad(psi)))*dolfin.dS

and instead compiles if I replace avg(grad(psi)) with avg(psi). But the dimensions of that seem wrong to me?

I also get a JIT error trying to integrate jump terms over the boundary,

dolfin.dot(dolfin.jump(phi), dolfin.jump(psi))*dolfin.dS

works, but

dolfin.dot(dolfin.jump(phi), dolfin.jump(psi))*dolfin.ds(0)

doesn't, returning error `TypeError: 52978288 format: a number is required, not NoneType
Minimum example:

import dolfin
import numpy as np

theta = 3.0*dolfin.pi/16.0
Omega = dolfin.as_vector((dolfin.cos(theta),dolfin.sin(theta)))
mesh = dolfin.UnitSquareMesh(20, 20)

# Define inflow boundary
class DirichletBoundary(dolfin.SubDomain):
    def inside(self, x, on_boundary):
        return (dolfin.near(x[0],0.0) or dolfin.near(x[1],0.0)) and on_boundary

# Define outflow boundary (disjoint covering of boundary with inflow boundary)
#   + We assume that Omega*n is positive on Neumann boundary, i.e. it is strictly outflow
class NeumannBoundary(dolfin.SubDomain):
    def inside(self, x, on_boundary):
        return (dolfin.near(x[0],1.0) or dolfin.near(x[1],1.0)) and on_boundary

# Mark different parts of boundary, 0 for Dirichlet and 1 for Neumann
boundary_parts = dolfin.MeshFunction("size_t", mesh, mesh.topology().dim()-1)
Gamma_D = DirichletBoundary()
Gamma_D.mark(boundary_parts, 0)
Gamma_N = NeumannBoundary()
Gamma_N.mark(boundary_parts, 1)
ds = dolfin.Measure("ds", domain=mesh, subdomain_data=boundary_parts)

# Function spaces
DG = dolfin.FunctionSpace(mesh,"DG", 1)

# Trial and test functions, normal vector, mesh size
psi = dolfin.TrialFunction(DG)
phi = dolfin.TestFunction(DG)
nor = dolfin.FacetNormal(mesh)

# Positive part of Omega*n
pos_nor = (dolfin.dot(Omega, nor) + np.abs(dolfin.dot(Omega, nor)))/2.0
# Negative part of Omega*n
neg_nor = (np.abs(dolfin.dot(Omega, nor)) - dolfin.dot(Omega, nor))/2.0

# DG, Eq. (3.4) 
#   - 5th line, jit error
#   - 7th and 8th line, dimensions of dot() error. If I make avg(grad(psi)) --> avg(psi),
#     it compiles, but this seems wrong...
bilinear_form = (dolfin.inner(psi,phi) + \
                dolfin.inner(dolfin.grad(psi), dolfin.grad(phi)) - \
                dolfin.inner(psi, dolfin.dot(Omega,dolfin.grad(phi))))*dolfin.dx + \
                dolfin.dot(dolfin.jump(phi), dolfin.jump(psi))*dolfin.dS + \
                dolfin.dot(dolfin.jump(phi), dolfin.jump(psi))*dolfin.ds(0) + \
                dolfin.dot(dolfin.jump(phi), pos_nor('+')*psi('+') - neg_nor('+')*psi('-') )*dolfin.dS + \
                dolfin.dot(dolfin.jump(phi), dolfin.avg(dolfin.grad(psi)))*dolfin.dS + \
                dolfin.dot(dolfin.jump(phi), dolfin.avg(dolfin.grad(psi)))*dolfin.ds(0) + \
                pos_nor*phi*psi*dolfin.ds

A = dolfin.assemble(bilinear_form)

Any help is greatly appreciated!

jmmal · Answer 1 · Jun 23, 2017

Hi,

to answer your two questions:
1) The reason why dolfin complains about the wrong dimensions, is that

dolfin.jump(phi)

is a scalar, whereas

dolfin.avg(dolfin.grad(psi))

is a vector.
I guess you want to get the jump over the facets, using the facet normals (and expect to get a vector quantity in return)? In your code, this is as simple as explicitly stating:

n = dolfin.FacetNormal(mesh)
dolfin.dot(dolfin.jump(phi,n), dolfin.avg(dolfin.grad(psi)))*dolfin.dS

So note the difference in the definition of the jump.

2) The definition of the jump and the average is slightly different for exterior boundary parts compared to interior boundaries. In the ufl language, you have to stay close to the mathematical definition to get the exterior terms correct. So the exterior term you are looking for can be defined as:

dolfin.dot(u*n, dolfin.grad(v))*ds(0)

Hopefully, this helps!

JIT error and dot() dimensions

1 Answer