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))
In [2]:
DrawFunction2D(uh, contour=False)
DrawFunction2D(uh, contour=True)
In [3]:
DrawFunction2D(uh1, contour=False)
DrawFunction2D(uh1, contour=True)
DrawFunction2D(uh2, contour=False)
DrawFunction2D(uh2, contour=True)