# 38. Stationary Transport Equation#

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

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

are given on the inflow boundary

## 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:

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

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.

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

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 elementsIn 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

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\):

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:

However, this leads only to

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

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.