Stabilized advection-diffusion equation

Stabilized advection-diffusion equation#

In the sequel we consider two different stabilization methods, leading to formulations of the form

\[ \begin{align} \text{find $u \in V$ such that} \quad a(u,v) + a_s(u,v) = l(v) + l_s(v) \quad \forall v \in V, \end{align} \]

where \(a_s\) and \(l_s\) are additional terms depending on the stabilization method.

The streamline upwind Petrov-Galerkin (SUPG) method is given by

\[ \begin{align} a_s(u,v) = \int_{\Omega} \tau ~ L(u) ~ \left( \mathbf{b} \cdot \nabla v \right) ~d\Omega \quad \mbox{and} \quad l_s(v) = \int_{\Omega} \tau ~ f ~ \left( \mathbf{b} \cdot \nabla v \right) ~d\Omega && \text{[SUPG]} \end{align} \]

while the Galerkin least squares (GLS) method is as follows,

\[ \begin{align} a_s(u,v) = \int_{\Omega} \tau ~ L(u) ~ L(v) ~d\Omega \quad \mbox{and} \quad l_s(v) = \int_{\Omega} \tau ~ f ~ L(v) ~d\Omega && \text{[GLS]} \end{align} \]

where we introduced the differential operator \(L(u) := - \kappa \nabla^2 u + \mathbf{b} \cdot \nabla u\).

Notice that in these formulations \(\tau\) stands for a piecewise function that is constant on each element of the tessellation associated to the computational domain \(\Omega\).

Formal Model#

Build a manifactured solution#

from sympde.expr import BilinearForm, LinearForm, integral
from sympde.expr     import find, EssentialBC, Norm, SemiNorm
from sympde.topology import ScalarFunctionSpace, Square, element_of
from sympde.calculus import grad, dot, laplace
from sympde.core     import Constant

from sympy import pi, sin, Tuple

domain = Square()

x,y = domain.coordinates

kappa = Constant('kappa', is_real=True)
b1 = 1.
b2 = 0.
b = Tuple(b1, b2)

L = lambda w: -kappa * laplace(w) + dot(b,grad(w))

ue = sin(pi*x)*sin(pi*y)
f = L(ue)

You can have the full expression of f by calling the TerminalExpr as follows

from sympde.expr import TerminalExpr

print(TerminalExpr(f, domain))

2*kappa*pi**2*sin(pi*x1)*sin(pi*x2) + 1.0*pi*sin(pi*x2)*cos(pi*x1)

Let’s go back now to our formal model.

V = ScalarFunctionSpace('V', domain)

u,v = [element_of(V, name=i) for i in ['u', 'v']]

# bilinear form
expr = kappa * dot(grad(v), grad(u)) + dot(b, grad(u)) * v
a = BilinearForm((u,v), integral(domain, expr))

# linear form
l = LinearForm(v, integral(domain, f*v))

tau   = Constant('tau', is_real=True)

# ...
s_supg = BilinearForm((v,u), integral(domain, tau * L(u) * dot(b, grad(v))))

l_supg = LinearForm(v, integral(domain, tau * f * dot(b, grad(v))))
# ...

# ...
s_gls = BilinearForm((v,u), integral(domain, tau * L(u) * L(v)))

l_gls = LinearForm(v, integral(domain, tau * f * L(v)))
# ...

# Dirichlet boundary conditions
bc = [EssentialBC(u, 0, domain.boundary)]

# Variational problem
equation_supg = find(u, forall=v, lhs=a(u, v) + s_supg(u,v), rhs=l(v) + l_supg(v), bc=bc)
equation_gls  = find(u, forall=v, lhs=a(u, v) + s_gls(u,v), rhs=l(v) + l_gls(v), bc=bc)

Discretization#

We shall need the discretize function from PsyDAC.

from psydac.api.discretization import discretize

degree = [2,2]
ncells = [8,8]

# Create computational domain from topological domain
domain_h = discretize(domain, ncells=ncells, comm=None)

# Create discrete spline space
Vh = discretize(V, domain_h, degree=degree)

# Discretize equation SUPG
equation_supg_h = discretize(equation_supg, domain_h, [Vh, Vh])

# Discretize equation GLS
equation_gls_h = discretize(equation_gls, domain_h, [Vh, Vh])

Solving the PDE#

Since the problem is not symmetric, we shall use gmres for the linear solver.

equation_supg_h.set_solver('gmres', info=False, tol=1e-8)
equation_gls_h.set_solver('gmres', info=False, tol=1e-8)

uh_supg = equation_supg_h.solve(kappa=1.e-1, tau=1.e-3)
uh_gls  = equation_gls_h.solve(kappa=1.e-1, tau=1.e-3)

Computing the error norm#

Computing the \(L^2\) norm#

u = element_of(V, name='u')

# create the formal Norm object
l2norm = Norm(u - ue, domain, kind='l2')

# discretize the norm
l2norm_h = discretize(l2norm, domain_h, Vh)

# assemble the norm
l2_error = l2norm_h.assemble(u=uh_supg)

# print the result
print('>>> SUPG = ', l2_error)

# assemble the norm
l2_error = l2norm_h.assemble(u=uh_gls)

# print the result
print('>>> GLS = ', l2_error)

>>> SUPG =  0.002139370600719921
>>> GLS =  0.00021873265193415864
Contents