Steady-state continuity equation on an extruded mesh

This demo showcases the use of extruded meshes, including the new regions of integration and the construction of sophisticated finite element spaces.

We now consider the equation

\[\begin{split}\nabla\cdot(\vec{u}q) &= 0 \\ q &= q_\mathrm{in} \quad \text{on} \ \Gamma_\mathrm{inflow},\end{split}\]

in a domain \(\Omega\), where \(\vec{u}\) is a prescribed vector field, and \(q\) is an unknown scalar field. The value of \(q\) is known on the ‘inflow’ part of the boundary \(\Gamma\), where \(\vec{u}\) is directed towards the interior of the domain. \(q\) can be interpreted as the steady-state distribution of a passive tracer carried by a fluid with velocity field \(\vec{u}\).

We apply an upwind DG method, as we saw in the previous example. Denoting the upwind value of \(q\) on interior facets by \(\widetilde{q}\), the full set of equations are then

\[-\int_\Omega \! q \vec{u_0} \cdot \nabla \phi \, \mathrm{d} x + \int_{\Gamma_{\mathrlap{\mathrm{ext, outflow}}}} \! \phi q \vec{u} \cdot \vec{n} \, \mathrm{d} s + \int_{\Gamma_\mathrm{int}} \! (\phi_+ \vec{u} \cdot \vec{n}_+ + \phi_- \vec{u} \cdot \vec{n}_-) \widetilde{q} \, \mathrm{d} S \quad = \quad -\int_{\Gamma_{\mathrlap{\mathrm{ext, inflow}}}} \phi q_\mathrm{in} \vec{u} \cdot \vec{n} \, \mathrm{d} s \quad \forall \ \phi \in V,\]

We will take the domain \(\Omega\) to be the cuboid \(\Omega = [0,1] \times [0,1] \times [0,0.2]\). We will use the uniform velocity field \(\vec{u} = (0, 0, 1)\). \(\Gamma_\mathrm{inflow}\) is therefore the base of the cuboid, while \(\Gamma_\mathrm{outflow}\) is the top. The four vertical sides can be ignored, since \(\vec{u} \cdot \vec{n} = 0\) on these faces.

We use an extruded mesh, where the base mesh is a 20 by 20 unit square, divided into triangles, with 10 evenly-spaced vertical layers. This gives prism-shaped cells.

from firedrake import *
m = UnitSquareMesh(20, 20)
mesh = ExtrudedMesh(m, layers=10, layer_height=0.02)

We will use a simple piecewise-constant function space for the unknown scalar \(q\):

V = FunctionSpace(mesh, "DG", 0)

Our velocity will live in a low-order Raviart-Thomas space. The construction of this is more complicated than element spaces that have appeared previously. The horizontal and vertical components of the field are specified separately. They are combined into a single element which is used to build a FunctionSpace.

# RT1 element on a prism
W0_h = FiniteElement("RT", "triangle", 1)
W0_v = FiniteElement("DG", "interval", 0)
W0 = HDivElement(TensorProductElement(W0_h, W0_v))
W1_h = FiniteElement("DG", "triangle", 0)
W1_v = FiniteElement("CG", "interval", 1)
W1 = HDivElement(TensorProductElement(W1_h, W1_v))
W_elt = W0 + W1
W = FunctionSpace(mesh, W_elt)

As an aside, since our prescibed velocity is purely in the vertical direction, a simpler space would have sufficed:

# Vertical part of RT1 element
# W_h = FiniteElement("DG", "triangle", 0)
# W_v = FiniteElement("CG", "interval", 1)
# W_elt = HDivElement(TensorProductElement(W_h, W_v))
# W = FunctionSpace(mesh, W_elt)

Or even:

# Why can't everything in life be this easy?
# W = VectorFunctionSpace(mesh, "CG", 1)

Next, we set the prescribed velocity field:

velocity = as_vector((0.0, 0.0, 1.0))
u = project(velocity, W)

# if we had used W = VectorFunctionSpace(mesh, "CG", 1), we could have done
# u = Function(W)
# u.interpolate(velocity)

Next, we will set the boundary value on our scalar to be a simple indicator function over part of the bottom of the domain:

x, y, z = SpatialCoordinate(mesh)
inflow = conditional(And(z < 0.02, x > 0.5), 1.0, -1.0)
q_in = Function(V)
q_in.interpolate(inflow)

Now we will define our forms. We use the same trick as in the previous example of defining un to aid with the upwind terms:

n = FacetNormal(mesh)
un = 0.5*(dot(u, n) + abs(dot(u, n)))

We define our trial and test functions in the usual way:

q = TrialFunction(V)
phi = TestFunction(V)

Since we are on an extruded mesh, we have several new integral types at our disposal. An integral over the cells of the domain is still denoted by dx. Boundary integrals now come in several varieties: ds_b denotes an integral over the base of the mesh, while ds_t denotes an integral over the top of the mesh. ds_v denotes an integral over the sides of a mesh, though we will not use that here.

Similiarly, interior facet integrals are split into dS_h and dS_v, over horizontal interior facets and vertical interior facets respectively. Since our velocity field is purely in the vertical direction, we will omit the integral over vertical interior facets, since we know \(\vec{u} \cdot \vec{n}\) is zero for these.

a1 = -q*dot(u, grad(phi))*dx
a2 = dot(jump(phi), un('+')*q('+') - un('-')*q('-'))*dS_h
a3 = dot(phi, un*q)*ds_t  # outflow at top wall
a = a1 + a2 + a3

L = -q_in*phi*dot(u, n)*ds_b  # inflow at bottom wall

Finally, we will compute the solution:

out = Function(V)
solve(a == L, out)

By construction, the exact solution is quite simple:

exact = Function(V)
exact.interpolate(conditional(x > 0.5, 1.0, -1.0))

We finally compare our solution to the expected solution:

assert max(abs(out.dat.data - exact.dat.data)) < 1e-10

This demo can be found as a script in extruded_continuity.py.