Advection-reaction (Hyperbolic)¶

This benchmark shows how to solve a steady-state 1^st order hyperbolic PDE with discontinuous boundary conditions.

Model problem¶

Equation solved: Steady-state advection-reaction equation

(1) $\vec\beta\cdot\nabla u + cu = f, \vec\beta = ( 10y^2 - 12x + 1, 1 + y ),\quad c = f = 0.$

Domain of interest: Square $[0, 1]\times [0, 1]$ .

Boundary conditions: Dirichlet, prescribed on the inflow parts of the domain boundary (i.e. where $\vec\beta\cdot\vec{n} < 0$ )

(2) $\begin{cases} 0 & \mbox{ for } (x = 0 \land 0.5 < y \leq 1)\,\lor\,(0.5<x\leq 1 \land y = 0),\\ 1 & \mbox{ for } (x = 0 \land 0 < y \leq 0.5)\,\lor\,(0 \leq x\leq 0.5 \land y = 0),\\ \sin^2(\pi y) & \mbox{ for } x = 1 \land 0\leq y \leq 1. \end{cases}$

Exact solution¶

Exact solution has been obtained by the method of characteristics, but it cannot be stated in a fully analytic form. Using Mathematica 7.0, the expression for the solution has been derived symbolically without any numerical approximation and its numerical evaluation at any point $(x,y)\in\Omega$ only amounts to calulating roots of 14^th order polynomials with coefficients defined by $x,y$ , evaluating standard analytical functions (exponential, sine) and performing standard algebraic operations. All these operations can be realized by Mathematica in machine precision and hence the solution may be considered exact for the purposes of comparison with the approximate results of Hermes. The source notebook is included [1].

For illustration purposes, the plots below were produced by the Matlab script plotu.m from file sol_101x101.map, which contains the solution values on a uniform grid of $101\times101$ points in $\Omega$ .

Methods overview¶

According to the classical theory, the advection term in eq. (1) propagates boundary data into the domain along the characteristics, i.e. curves $\vec x(s)$ satisfying $\mathrm d\vec x/\mathrm ds = \vec\beta$ . Hence, (1) with jump-discontinuous boundary conditions (2) has a unique solution with jump-discontinuities normal to the characteristics passing through the point of discontinuity at the boundary. It is well known that the classical conforming FEM solution realized by globally continuous basis functions is not likely to produce satisfactory results – unlike the problems with steep but continuous solution gradients, which could be resolved by sufficiently high h-refinement, spurious Gibbs oscillations will be introduced in the vicinity of the jumps and spread to the more distant parts of the domain – see the figure below.

Solution by the continuous Galerkin method without stabilization.

The two mostly used Galerkin methods which circumvent these stability issues will be briefly introduced now.

Stabilized continuous Galerkin methods. These methods effectively add artificial diffusion to the equation, changing its behavior to that of an advection-reaction-diffusion equation with a globally continuous solution. Typical methods from this category include the Streamline upwind Petrov-Galerkin (SUPG), Galerkin least squares (GLS) or Subgrid scale (SGS) methods (see e.g. [C98] and the references therein). The discretization procedure (carried out in the standard $H^1(\Omega)$ space) leads to much smaller linear algebraic systems comparing to the methods from the other category. However, it is usually very difficult to tune the amount of added diffusion so that the solution does not diverge too much from that of the original 1^st order problem, particularly for the higher-order multidimensional schemes.
Discontinuous Galerkin (DG) methods. In the DG methods, discretization is carried out in the $L^2(\Omega)$ space, using basis functions which are smooth inside each element but discontinuous across element interfaces. Relaxation of the interelement continuity allows for capturing the solution jumps and results in a stable scheme without the need for additional tuning, but at the expense of a bigger algebraic system to solve and more difficult assembling.

This benchmark implements the SUPG and DG methods (and the classical continuous FEM for comparison). Note that both methods may be used to ensure stability, i.e. that numerical oscillations will be fully contained in a close vicinity of the point where they appear, but they do not per se prevent the oscillations from actually occuring. In order to do so, monotonicity of the scheme has to be ensured as well by proper discontinuity capturing techniques, which are however not implemented in this benchmark.

Weak forms¶

Streamline upwind Petrov-Galerkin¶

The bilinear form for the SUPG discretization of problem (1), (2) is composed of three parts. The first one is formally obtained by multiplying eq. (1) by an $L^2$ -integrable test function and integrating over the whole domain:

$\int_\Omega (\vec\beta\cdot\nabla u + c)v \,\mathrm{d}x.$

Note that we do not apply the Green’s theorem and seek the strong solution, which lies in $L^2(\Omega)$ together with its streamline derivative $\vec\beta\cdot\nabla u$ [2]. Space of such functions contains $H^1(\Omega)$ and in particular its finite-dimensional subspace of piecewise continuous polynomials up to a specified order, which we use for the practical implementation.

The second part reads

(3) $\int_\Omega (\vec\beta\cdot\nabla u + cu - f)\, \tau \,\vec\beta\cdot\nabla v \,\mathrm{d}x,$

where $\tau$ must be judiciously chosen to ensure the stability. Note that since (3) contains the whole residual of (1), consistency is not broken by this additional contribution to the whole SUPG bilinear form.

Appropriate choice of parameter $\tau$ , so that the scheme is neither over-stabilized nor under-stabilized, is the major concern when implementing the SUPG method. Theoretically justified rules are known mostly for first order accurate, one dimensional, advection-reaction-diffusion problems. Therefore, for our isotropic finite element grid, the classical expression (see e.g. [C98]) is extended by letting the diffusion coefficient vanish and taking a square of the result, with the 1D advection coefficient replaced in each element $K$ by $||\vec\beta||_{L^2(K)\times L^2(K)}$ . The result is

$\tau = \frac{\mathrm{diam}(K)^2}{4 ||\vec\beta||^2_{L^2(K)\times L^2(K)}}$

and is working reasonably well for the current problem. Nevertheless, the reader is encouraged to derive and experiment with his own expressions.

The final part of the SUPG bilinear form, together with the corresponding linear form, enforces the Dirichlet boundary conditions on the inflow boundaries in an $L^2$ -integral sense. Although this form has been traditionally used in literature rather for mathematical analysis than for practical computation, we have found it advantageous for the latter purpose as well since with open quadrature rules currently used in Hermes, the problematic evaluation of the boundary condition and approximate solution at the vertices of discontinuity is hence avoided.

Discontinuous Galerkin¶

There are several possibilities how to formulate the DGM. We choose that presented and analysed in [BMS04].

The weak solution is well defined in the broken Sobolev space of functions $u\in L^2(\Omega)$ such that $u\in H^1(K)$ for every element $K$ , whose finite-dimensional subspace suitable for the FE discretization is represented in Hermes by class L2Space. Since $u$ is not expected to be continuous across element interfaces, the Green’s theorem has been applied element-wise. The consequence is the presence of surface integrals of basis and test functions, or more precisely of their arithmetic averages and jumps across element interfaces. If the surface form representing these integrals is added with the special marker H2D_DG_INNER_EDGE.

The values of the traces of the shape functions from both sides of an interface may then be obtained at the quadrature points along the interface.

The jump and average operators defined as in the macros AVG and JUMP above. Note that in order to apply the Green’s theorem, a transition to the conservative form of (1) has been performed using the product rule for derivatives (utilizing differentiability of $\vec\beta$ ), eventually leading to the term $-\nabla\cdot\vec\beta$ added to the reaction term $c$ in the final weak form.

Sample results¶

In this section, we present graphs of the solutions we obtained using the various adaptivity schemes. All approaches started from the unrefined, 4-element mesh, used the following heuristic setting

STRATEGY = 0;
THRESHOLD = 0.20;
CONV_EXP = 1.0;

and were ended when the refined finite element mesh contained more than 90000 dof. Manual setting of the error weighting applied during refinement also proved to be neccessary in order to prioritize h-refinement and hence better resolve the discontinuity regions.

Streamline upwind Petrov-Galerkin¶

h-adaptivity, P = 1 uniformly¶

h-adaptivity, P = 2 uniformly¶

hp-adaptivity¶

Discontinuous Galerkin¶

h-adaptivity, P = 0 uniformly¶

h-adaptivity, P = 1 uniformly¶

hp-adaptivity¶

Convergence comparisons¶

Below we compare the convergence of the various adaptive methods using two metrics.

Integral value of the weighted flux at the outflow boundary (the top edge of the square $\Omega$ ):

$\int_{\Gamma_\mathrm{out}} \vec\beta\cdot\vec n uw\,\mathrm{d}s,\quad \Gamma_\mathrm{out} = \{(x,y):\,\vec\beta(x,y)\cdot\vec n(x,y) > 0\},$

where the weighting function has been chosen as in [HRS00]: $w(x,y) = \sin(\pi x/2)\;$ for $\;(x,y)\in [0,1]\times {1}$ .
Relative $L^2(\Omega)$ error w.r.t. the exact (semi-analytic) solution:

$\frac{||u_{\mathrm{ex}} - u_h||_{L^2(\Omega)}}{||u_\mathrm{ex}||_{L^2(\Omega)}}$

In order to calculate this quantity, the exact solution has been evaluated at the $(50+51)\times (50+51)$ nodal points of the two-dimensional 50^th-order Gauss quadrature rule with Kronrod extension and saved together with the corresponding quadrature weights to file sol_GaussKronrod50.map. There is a class SemiAnalyticSolution responsible for loading the file and repeatedly calculating the norm, but be warned that since the latter operation involves a call to Solution::get_pt_value, computation of this metric considerably prolongates each adaptation step (particularly when there are many small low-order elements).

References¶

[BMS04]

F. Brezzi, L. D. Marini, and E. Suli: Discontinuous Galerkin methods for first-order hyperbolic problems. http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.4.333

[HRS00]

P. Houston, R. Rannacher, E. Süli: A posteriori error analysis for stabilised finite element approximations of transport problems. Comput. Meth. Appl. Mech. Engrg. 190 (2000), pp. 1483-1508.

[C98]

(1, 2) R. Codina: Comparison of some finite element methods for solving the diffusion-convection-reaction equation. Comput. Meth. Appl. Mech. Engrg. 156 (1998), pp. 185-210.

Footnotes

[1]	If you do not have Mathematica installed, a limited view is possible by the free Mathematica player.

[2]	This expresses the fact that the solution may be possibly discontinuous across certain characteristic curves.

Advection-reaction (Hyperbolic)¶

Model problem¶

Exact solution¶

Methods overview¶

Weak forms¶

Streamline upwind Petrov-Galerkin¶

Discontinuous Galerkin¶

Sample results¶

Streamline upwind Petrov-Galerkin¶

h-adaptivity, P = 1 uniformly¶

h-adaptivity, P = 2 uniformly¶

hp-adaptivity¶

Discontinuous Galerkin¶

h-adaptivity, P = 0 uniformly¶

h-adaptivity, P = 1 uniformly¶

hp-adaptivity¶

Convergence comparisons¶

References¶

Table Of Contents

Previous topic

Next topic

Navigation

Advection-reaction (Hyperbolic)¶

Model problem¶

Exact solution¶

Methods overview¶

Weak forms¶

Streamline upwind Petrov-Galerkin¶

Discontinuous Galerkin¶

Sample results¶

Streamline upwind Petrov-Galerkin¶

h-adaptivity, P = 1 uniformly¶

h-adaptivity, P = 2 uniformly¶

hp-adaptivity¶

Discontinuous Galerkin¶

h-adaptivity, P = 0 uniformly¶

h-adaptivity, P = 1 uniformly¶

hp-adaptivity¶

Convergence comparisons¶

References¶

Table Of Contents

Previous topic

Next topic

Quick search

Navigation