Below you can find a link to the MATLAB code calculating the advection-diffusion problem by the finite element method with stabilization by minimizing the residual for the model Eriksson-Johnson problem.
First, let us recall the advection-diffusion problem, the Eriksson-Johnson model problem (described, for example, in [1] defined on a square domain \( \Omega = [0,1]^2 \) in the following way: We seek \( u \) such as
\( a(u,v)=l(v) \forall v \) where
\( a(u,v) =\int_{\Omega} \beta_x(x,y) \frac{\partial u(x,y) }{\partial x } dxdy + \int_{\Omega} \beta_y(x,y) \frac{\partial u(x,y)}{\partial y } dxdy \\ +\epsilon \int_{\Omega} \frac{\partial u(x,y) }{\partial x} \frac{\partial v(x,y)}{\partial x } dxdy +\epsilon \int_{\Omega} \frac{\partial u(x,y)}{\partial y } \frac{\partial v(x,y)}{\partial y } dxdy \)
\( l(v) = \int_{\partial \Omega } f(x,y) v dxdy \)
\( f(x,y)=sin(\pi y)(1-x) \) is an extension of the Dirichlet boundary condition to the entire area, while
\( \beta = (1,0) \)
represents the wind blowing from left to right, while \( \epsilon = 10^{-2} \) denotes diffusion coefficient. The quantity \( Pe=1/ \epsilon = 100 \) is called the Peclet number, and it defines the sensitivity of the advection-diffusion problem.
In line 830 we specify the diffusion coefficient ϵ=10−2. Its inverse is the Peclet number \( Pe=\frac{1}{\epsilon} \).
We construct the knot vector \( knot\_trial\_x = [0 \quad 0 \quad 0 \quad 1 \quad 2 \quad 2 \quad 2] \) for the approximation space, and vector of point along x axis, between od 0 and 1, \( points\_x = [0 \quad 0.5 \quad 1] \). At these points, we construct the basis of B-spline functions for approximation of the Eriksson-Johnson problem. In a similar way, we define the knot vector for the test space \( knot\_test\_x = [0 \quad 0 \quad 0 \quad 0 \quad 1 \quad 2 \quad 2 \quad 2 \quad 2] \).
According to the idea of the residual minimization method, there must be more testing functions than approximating functions. In our example, we use the quadratic B-spline function with C1 continuity for the approximation space (trial) and the third degree B-spline function with C1 continuity for the test space (test).
The node vectors for the approximating and testing space, and the points on which we span the node vectors, are defined separately for the x axis and for the y axis.
The code is activated by opening it in Octave and typing a command
\( advection\_igrm\_adapt \)
After a moment of calculation, the code opens an additional window and draws a numerical and exact solution in it. In the second window, the code draws a residuum, computed additionally by the residuum minimization method. It is a measure of the local solution error.
The code calculates the residuum error \( Residual norm: L2 0.015808, H1 0.100417 \)
The code computes also the L2 and H1 norms of the solution
\( Error: L2 45.90 percent, H1 61.85 percent \)
and compares to L2 and H1 norm of the exact solution.
\( Best possible: L2 2.02 percent, H1 14.20 percent \)
You can see that in order to improve the accuracy of the solution, it is necessary to increase the mesh.
Exercise 1: Uniform increase of trial grid and test
Content of exercise:
Please increase the number of points to five in x and y directions, keeping the degrees of B-spline polynomials in the trial and test space. Please check how this operation improves the accuracy of the solution.
similarly in the x direction.
Residual norm: L2 0.009895, H1 0.080162
Error: L2 37.83percent H1 50.92percent
Best possible: L2 0.41percent H1 36.77percent
Exercise 2: Adaptive enlargement of the trial grid and test
Content of exercise:
Please increase adaptively the number of points towards the right bank [0 0.5 0.9 0.95 1] where the coastal layer is, keeping the degrees of B-spline polynomials in the trial and test space. Please check how this operation improves the accuracy of the solution.
Exercise 3: Adaptive increasing the trial grid and test by breaking in half the element closest to the singularity for a large Peclet number
Content of exercise:
Please modify the problem so that the Peclet number is one million. Please start from the grid [0 0.5 1] in the x direction and [0 0.25 0.5 0.75 1] in the y direction, and increase adaptively the number of points towards the right edge where the boundary layer is, keeping the degrees of B-spline polynomials in the trial space and test. Please check how this operation improves the accuracy of the solution.
Solution:
We enter \( \epsilon=10^{-6} \) and we modify the grid according to the sequence