31. Instationary Transport Equation

\(\DeclareMathOperator{\opdiv}{div}\)

The time dependent transport equation is to find an \(u = u(x,t)\) such that

\[ \frac{\partial u}{\partial t} + \opdiv bu = f \]

with an given initial condition \(u(x,t=0) = u_0(x)\), and a boundary condition on \(\Gamma_{in}\).

The explicit Euler time-discretization method is to find \(u^n \approx u(t_n)\):

\[ \frac{1}{t_n - t_{n-1}} \int_\Omega (u^n - u^{n-1}) v + A^{DG} (u^{n-1},v) = f(v) \]

Here, \(A^{DG}(.,.)\) is the bilinear-form of the stationary transport equation.

To compute the new coefficient vector for the new \(u^n\), we have to solve an equation with the mass matrix

\[ M (u^n-u^{n-1}) = \tau \, (f - A^{DG} u^{n-1}) \]

Here, we realize the second advantage of DG methods: Since basis functions are defined element by element, the mass matrix is block diagonal. Thus, it is cheap to invert it.

from ngsolve import *
from netgen.geom2d import unit_square
from ngsolve.webgui import Draw
mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))

b = CoefficientFunction( (y-0.5, 0.5-x) )
Draw (b, mesh, "wind", vectors={"grid_size" : 20 });

fes = L2(mesh, order=3)

u = fes.TrialFunction()
v = fes.TestFunction()
a = BilinearForm(fes, nonassemble=True)

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

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

f = LinearForm(fes)
f.Assemble()

gfu = GridFunction(fes)
gfu.Set(exp(-10**2*((x-0.5)*(x-0.5) +(y-0.75)*(y-0.75))))

scene = Draw(gfu, min=0, max=1, order=3, autoscale=False)

tau = 0.001
tend = 50
t = 0
cnt = 0

w = gfu.vec.CreateVector()
# the mass matrix (implemented matrix-free)
invm = fes.Mass(rho=1).Inverse()
# SetNumThreads(3)
with TaskManager():
    while t < tend:
        # print (t)
        # apply the transport operator 
        a.Apply (gfu.vec, w)
        gfu.vec.data -= tau * invm * w
        if cnt%20 == 0:
            scene.Redraw()    
        t += tau
        cnt +=1