Consider a different scalar PDE: Convection-Diffusion equations [Mika Meyer, Kolja Straatman, Pascal Irmer]¶

Consider the (reaction-)convection-diffusion equation $$ -\nabla \cdot (\alpha \nabla u) + \beta \cdot \nabla u + \gamma u = f $$ in $\Omega$, equipped with suitable boundary conditions.

Implement a FEM solver for this problem and carry out numerical experiments to investigate the influence of the convection term on the solution.

Consider the following standard test problem:

  • Set $ \gamma = 0$, and $\beta = (1,1)$, $\Omega = (0,1)^2$
  • Use homogeneous Dirichlet boundary conditions: $u|_{\partial \Omega} = 0$, see Dirichlet demo notebook.
  • Consider the so-called mesh Péclét number $Pe = \frac{2 \| \beta \|_\infty}{\gamma h}$, where $h$ is the (initially uniform) mesh size.
  • Take the details as in this paper, page 19 ff., i.e. take the exact solution and compute (verify) the corresponding r.h.s. and compute the numerical solutions.
  • Implement the (missing) convection integral for $\int_\Omega \beta \cdot \nabla u v \, dx $ in the bilinear form (including tests and notebook(s))
  • Do numerical studies over $h$ and/or $Pe$ with different FE spaces
  • Document the results in proper notebook(s)
  • Use strechted grids (by introducing a mapping in the construction of the mesh) with increasingly small elements towards $x=1$ and $y=1$ and investigate the influence on the solution. (including documentation in notebook(s))

Weak Formulation: $$\int_\Omega \alpha \nabla u \nabla v dx + \int_\Omega \beta \nabla u v dx +\int_\Omega \gamma u v dx=\int_\Omega fv dx + \int_{\partial\Omega} \alpha \nabla u v ds$$ convection term: $$\int_\Omega \beta \nabla u v dx$$ implementation

In [1]:
from import_hack import *
from methodsnm.mesh_2d import *
from methodsnm.visualize import *
from methodsnm.fes import *
from methodsnm.forms import *
from methodsnm.formint import *
from numpy import pi, cos, sin, arctan, tanh, tan, e

beta=np.array([0,0])
beta1=np.array([1,1])
beta2=np.array([-1,-1])
gamma = 0
alpha = 0.1
N = 5
uh = Solver(N,beta,alpha,gamma,lambda x: 1,function=lambda x,y: (x,y))
uh1 = Solver(N,beta1,alpha,gamma,lambda x: 1,function=lambda x,y: (x,y))
uh2 = Solver(N,beta2,alpha,gamma,lambda x: 1,function=lambda x,y: (x,y))
source module for methodsNM imported.
Warning: Using numerical differentiation for deriv evaluation in <class 'methodsnm.fe_2d.P3_Triangle_FE'> object.
In [2]:
DrawFunction2D(uh, contour=False)
DrawFunction2D(uh, contour=True)
No description has been provided for this image
No description has been provided for this image
In [3]:
DrawFunction2D(uh1, contour=False)
DrawFunction2D(uh1, contour=True)
DrawFunction2D(uh2, contour=False)
DrawFunction2D(uh2, contour=True)
No description has been provided for this image
No description has been provided for this image