Stationary Transport Equation

30. Stationary Transport Equation#

\(\DeclareMathOperator{\opdiv}{div}\) We consider the stationary transport equation (convection 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 \}. \]

30.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} \operatorname{div} (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 the 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))

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

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

u = fes.TrialFunction()
v = fes.TestFunction()

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)*(y-0.5))*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);
Contents