Advection-diffusion equation

Advection-diffusion equation#

We consider the advection-diffusion problem consisting of finding a scalar-valued function \(u\) such that

\[\begin{split} \begin{align} \begin{cases} - \kappa \nabla^2 u + \mathbf{b} \cdot \nabla u = f &\text{in $\Omega$}, \\ u = 0 &\text{on $\partial \Omega$}. \end{cases} \end{align} \end{split}\]

The Variational Formulation#

The variational formulation reads

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

where

  • \(V \subset H_0^1(\Omega)\),

  • \(a(u,v) := \int_{\Omega} \left( \kappa \nabla u \cdot \nabla v + \left( \mathbf{b} \cdot \nabla u \right) v \right) d\Omega\),

  • \(l(v) := \int_{\Omega} f v~d\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))

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

# Variational problem
equation   = find(u, forall=v, lhs=a(u, v), rhs=l(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
equation_h = discretize(equation, domain_h, [Vh, Vh])

Solving the PDE#

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

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

uh = equation_h.solve(kappa=1.e-1)

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)

# print the result
print(l2_error)

0.00022156219679932395

Computing the \(H^1\) semi-norm#

# create the formal Norm object
h1norm = SemiNorm(u - ue, domain, kind='h1')

# discretize the norm
h1norm_h = discretize(h1norm, domain_h, Vh)

# assemble the norm
h1_error = h1norm_h.assemble(u=uh)

# print the result
print(h1_error)

0.013033907531543579

Computing the \(H^1\) norm#

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

# discretize the norm
h1norm_h = discretize(h1norm, domain_h, Vh)

# assemble the norm
h1_error = h1norm_h.assemble(u=uh)

# print the result
print(h1_error)

0.013035790553238334
Contents