Stationary Transport Equation

38. Stationary Transport Equation#

\(\DeclareMathOperator{\opdiv}{div}\) We consider the stationary transport equation

\[ \opdiv bu = f \qquad \text{ on } \Omega \]

Here, the vector field \(b\) is the given wind. It is a model for continuously injecting some concentration into a given flow, for example milk into coffee.

Boundary conditions

\[ u = u_{in} \]

are given on the inflow boundary

\[ \Gamma_{in} = \{ x \in \partial \Omega : b \cdot n < 0 \}. \]

38.1. Discontinuous Galerkin method#

Finite volume methods are very popular for transport equations. Their strength is the upwind technique for a stable discretization. A DG method combines the advantages of finite elements and finite volumes. It can be seen as the extension of finite volume methods to higher order.

A DG method uses discontinuous trial and test finite element spaces, for example piecewise polynomials. It is derived as follows:

Multiply the equation by element-wise smooth test functions, integrate:

\[ \int_{\Omega} \opdiv (bu) v = \int_\Omega f v \]

and integrate by parts on evey element \(T\):

\[ \sum_{T \subset \Omega} - \int_T b u \nabla v + \int_{\partial T} b_n u v = \int_\Omega f v \]

Here, \(n\) is the outward going normal vector of the element. The true solution \(u\) is continuous over element facets (edges in 2D and faces in 3D). But, for a discontinuous trial function we have to use the \(u\) on the boundary from the left, or from the right element. The stable decision is to use \(u\) from the upwind element, i.e. the element where \(b_n > 0\), i.e. the wind blows out of the element.

\[\begin{split} u^{up} := \left\{ \begin{array}{cl} u_T & \text{ if } b_n > 0 & \text{element outflow boundary } \\ u_{T,other} & \text{ if } b_n < 0 & \text{ element inflow boundary } \end{array} \right. \end{split}\]

Thus, the DG discretization method is as follows: \(V_h\) a piecewise polynomial, discontinuous finite element space. Find \(u_h \in V_h\) such that

\[ \sum_{T \subset \Omega} - \int_T b u_h \nabla v_h + \int_{\partial T} b_n u_h^{up} v_h = \int_\Omega f v_h \]

for all test functions \(v_h \in V_h\).

from ngsolve import *
from ngsolve.webgui import Draw
mesh = Mesh(unit_square.GenerateMesh(maxh=0.05))

Choose some wavy wind:

b = CF( (1,0.5*sin(2*6.28*x)))
Draw (b, mesh, "wind", min=0.8, max=1.2, vectors = { "grid_size" : 20 });

The upwind term is specified in NGSolve as

n = specialcf.normal(2)
uup = IfPos(b*n, u, u.Other())
a += b*n*uup*v * dx(element_boundary=True)

  • The IfPos function checks whether the first argument is positive. If so, it returns the secod, if not, it returns the third argument.

  • The dx(element_boundary=True) integrates over all element boundaries. The loop in the assembling is over all elements

  • In this context, the normal vector is the outward-pointing normal vector to the element

The source term is integrated over traces of the L2 finite element space. This would not work with the usual ds. The flag skeleton=True gives also access to the volume element next to the surface element.

fes = L2(mesh, order=2, dgjumps=True)

u,v = fes.TnT()

a = BilinearForm(fes)
a += -b*u*grad(v)*dx

# the upwind-term:
n = specialcf.normal(2)
uup = IfPos(b*n, u, u.Other())
a += b*n*uup*v * dx(element_boundary=True)
a.Assemble()

f = LinearForm(fes)
f += exp(-20**2*(y-0.5)**2)*v * ds(skeleton=True,definedon=mesh.Boundaries("left"))
f.Assemble();

gfu = GridFunction(fes)
gfu.vec.data = a.mat.Inverse() * f.vec

Draw (gfu, order=3);

38.2. Stability of the method#

The bilinear-form is \(B(.,.) : V \times W \rightarrow {\mathbb R}\), with

\[ B(u,v) = \int u b \cdot \nabla v \, dx \]

We choose \(\| u \|_V = L_2\), and \(\| v \|_W = \| b \cdot \nabla v \|_{L_2}\). The later is a norm, if the wind \(b\) is such that the every point is connected with the inflow boundary.

One inf-sup condition is very simple: Choose \(u = b \cdot \nabla v\):

\[ \sup_u \frac{| B(u,v) | }{ \| u \| } \geq \frac{ \| b \cdot \nabla v \|^2 }{ \| b \cdot \nabla v \| } = \| v \|_W \]

With DG-methods one can mimic this inf-sup proof for the discrete spaces: choose \(u_h\) as element-wise derivative. This would not be possible by continuous finite elements:

\[ u_h^1 = b \cdot \nabla_h v_h \]

However, this leads only to

\[ B(u_h^1, v_h) \succeq \| b \cdot \nabla v_h \|^2 - \sum_E \| [ v_h ] \|_{L_2}^2, \]

with a bad, negative term on the right. Now we choose a second candidate \(u_h^2 = v_h\), and find

\[ B(u_h^2, v_h) \succeq \sum_E \| [ v_h ] \|_{L_2}^2, \]

Here the upwind-choice comes into the analysis. Finally, we take a combination \(u_h = u_h^1 + \alpha u_h^2\) such that the second choice compensates the negative term from the first choice. Details are found in the lecture notes.

Contents