43. Splitting Methods for the time-dependent convection diffusion equation

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

We now want to solve a time-dependent convection diffusion equation:

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

• we use a DG method for the transport term, and an HDG method for the diffusion term

• we treat the transport term explicit. This allows to treat non-linear terms as in the Navier-Stokes equation easily. Also, solving non-symmetric systems is more difficult

• we treat the diffusion term implicit. An explizit treatment of the second-order term would lead to a severe time-step restriction. It is a symmetric term, for which fast solvers are available.

The equation for the time-step becomes:

\[\begin{split} \frac{1}{\tau} \left( \begin{array}{cc} M & 0 \\ 0 & 0 \end{array} \right) \left( \begin{array}{c} u^n-u^{n-1} \\ \widehat u^n - \widehat u^{n-1} \end{array} \right) + A^{diff} \left( \begin{array}{c} u^n \\ \widehat u^n \end{array} \right) + \left( \begin{array}{cc} A^{conv} & 0 \\ 0 & 0 \end{array} \right) \left( \begin{array}{c} u^{n-1} \\ \widehat u^{n-1} \end{array} \right) = f \end{split}\]

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

order = 3
fes1 = L2(mesh, order=order)
fes2 = FacetFESpace(mesh, order=order, dirichlet=".*")
fes = fes1*fes2

b = CoefficientFunction( (y-0.5, 0.5-x) )

tau = 1e-3
eps = 1e-4
tend = 10

u,uhat = fes.TrialFunction()
v,vhat = fes.TestFunction()

h = specialcf.mesh_size
n = specialcf.normal(2)
alpha = 2
dS = dx(element_boundary=True)

diffusion = eps*grad(u)*grad(v)*dx \
    + eps*(-n*grad(u)*(v-vhat)-n*grad(v)*(u-uhat)) * dS \
    + eps*alpha*(order+1)**2/h*(u-uhat)*(v-vhat) * dS

uup = IfPos(b*n, u, u.Other(bnd=0))
convection = -b*u*grad(v)*dx + b*n*uup*v * dS

adiff = BilinearForm(diffusion)

aconv = BilinearForm(fes, nonassemble=True)
aconv += convection

mstar = BilinearForm(u*v*dx + tau*diffusion)

f = LinearForm(fes)

mstar.Assemble()
adiff.Assemble()
f.Assemble();

gfu = GridFunction(fes)
gfu.components[0].Set(exp(-10**2*((x-0.5)*(x-0.5) +(y-0.75)*(y-0.75))))
scene = Draw(gfu.components[0], min=0, max=1, order=3, autoscale=False)

convu = gfu.vec.CreateVector()
w = gfu.vec.CreateVector()
r = gfu.vec.CreateVector()

inv = mstar.mat.Inverse(fes.FreeDofs(), inverse="sparsecholesky")
t = 0
cnt = 0
SetNumThreads(4)
while t < tend:
    t += tau
    aconv.Apply(gfu.vec, convu)
    r.data = f.vec - convu - adiff.mat * gfu.vec
    w.data = inv * r
    gfu.vec.data += tau*w
    cnt += 1
    if cnt % 10 == 0:
        scene.Redraw()

For higher order IMplicit-EXplicit (IMEX) methods see:

Implicit-explicit Runge-Kutta methods for time-dependent pdes by Uri M. Ascher, Steven J. Ruuth, and Raymond J. Spiteri, Applied Numerical Mathematics, 25, 1997