Construction of a p-Adaptive Continuous Residual Distribution Scheme

Abstract

A p-adaptive continuous residual distribution scheme is proposed in this paper. Under certain conditions, primarily the expression of the total residual on a given element K into residuals on the sub-elements of K and the use of a suitable combination of quadrature formulas, it is possible to change locally the degree of the polynomial approximation of the solution. The discrete solution can then be considered non continuous across the interface of elements of different orders, while the numerical scheme still verifies the hypothesis of the discrete Lax–Wendroff theorem which ensures its convergence to a correct weak solution. We detail the theoretical material and the construction of our p-adaptive method in the frame of a continuous residual distribution scheme. Different test cases for non-linear equations at different flow velocities demonstrate numerically the validity of the theoretical results.

Appendices

Appendix 1: Implicit Numerical Solver

We need to solve the system of equations (20), written in a compact way as:

$$\begin{aligned} R({\mathbf {u}}_{h})=0. \end{aligned}$$

(60)

This problem is first relaxed as:

$$\begin{aligned} \frac{d{\mathbf {u}}_{h}}{dt}=-R({\mathbf {u}}_{h}). \end{aligned}$$

(61)

To approximate this time derivative, we use the backward Euler formula, and the problem becomes:

$$\begin{aligned} \frac{{\mathbf {u}}^{n+1}_{h}-{\mathbf {u}}^{n}_{h}\mathbf {{}}_{}}{\Delta t^{n}}=-R({\mathbf {u}}^{n+1}_{h}), \ n=0,1,2,\ldots \end{aligned}$$

(62)

Thus, we have to solve a non linear problem at each time step $n$. We use for that a Newton method, that when applied to (60) reads:

$$\begin{aligned} {\mathbf {u}}^{k+1}_{h}={\mathbf {u}}^{k}_{h} -J^{-1}R({\mathbf {u}}^{k}_{h}),\;\;k=0,1,2, \end{aligned}$$

(63)

where

$$\begin{aligned} J=\frac{\partial R({\mathbf {u}}^{k}_{h})}{\partial {\mathbf {u}}} \end{aligned}$$

(64)

is the Jacobian of R. In practice, the Jacobian J is computed from a first order scheme with the stabilization term. And so, with our Newton method applied to (62) , the linear system we have to solve at each time step $n$ reads:

$$\begin{aligned}&[\frac{1}{\Delta t^{n}}+J({\mathbf {u}}^{n}_{h_{k}})]\Delta {\mathbf {u}}^{n}_{h_{k}}=-R({\mathbf {u}}^{n}_{h_{k}}) \end{aligned}$$

(65)

$$\begin{aligned}&{\mathbf {u}}^{n}_{h_{k+1}}={\mathbf {u}}^{n}_{h_{k}}+\Delta {\mathbf {u}}^{n}_{h_{k}},\;\;k=0,1,2\ldots \end{aligned}$$

(66)

In practice, for each time step $n$, we use only one iteration on $k$, which seems enough to reach convergence of the nodal residuals to the zero machine.

Appendix 2: Comparison with Constrained Approximation for the Finite Element Method

In [8] the authors present a mesh adaptation method called “constrained approximation“ in the frame of a Petrov–Galerkin finite-element method. The idea behind the method is to use a non regular mesh (non conforming and of various order) and to constrain the value of the discrete solution at the interface shared by elements of different size (h-adaptation) and different order (p-adaptation) to keep the solution continuous across the interface.

Fig. 28

From left to right: initial mesh, h-adaptation and p-adaptation

In Fig. 28 inspired from [8], the authors give an example of h-refinement and p-refinement. In the h-refinement case, the two triangles share a common edge with two different sizes, which generates a hanging node in the middle of the edge. In the p-refinement case, the triangles are of two different polynomial orders (${\mathbb {P}}1$ and ${\mathbb {P}}2$) and so a hanging node is generated in the middle of the edge. In both cases, the solution is constrained at the hanging node so that it remains continuous.

The authors use a hierarchical basis of shape functions [24]. For example, in dimension two this basis is constituted of functions associated with the vertices, the edges and the element and are called accordingly vertex, edge and bubble functions. They are built incrementally with respect to the desired order of the approximation, by using kernel functions in the definition of the edge and bubble functions. The main advantages of this basis over the classical Lagrange basis functions, is that the polynomial order of the approximation can be changed easily (by adding or removing kernel functions) and in the case of p-adaptation, like in Fig. 28, the hanging node is naturally suppressed, by changing to one the polynomial order of the shape function associated with the edge at the interface in the ${\mathbb {P}}_2$ element.

Once the constrains are established, the next step is to apply these constrains. The authors present a complete set of rules to impose the continuity of the solution in the case of h-adaptation and p-adaptation. The authors suggest to follow a rule called ”1-irregularity“ rule (only one hanging node between two elements) in order to avoid too complex situations, at least for h-adaptation.

We believe that this method is very different from the method we have proposed in this paper. The first difference is that the generation of hanging nodes makes the mesh non conforming. In a RD scheme, the question of how this hanging node should be treated does not seem evident. The distribution schemes (like for example (21) and (35)) would probably have to be redefined to take into account the hanging nodes, so that the scheme remains conservative (which is an essential ingredient for the convergence) and still exhibit the same behavior (for example (non) diffusive or (non) oscillatory). Even if it is still possible to use a classical Lagrange basis, it is advised in [8] to use a hierarchical basis of functions, as such a basis is more suited for p-adaptation. So, the compatibility of a hierarchical basis function with RD schemes should be analyzed. Finally, there is the question as to whether the hanging nodes should be taken into account in the system of equations (20), and for the implementation of an implicit RD scheme, what the consequences would be for the computation of the Jacobian matrix. We are not saying that the constrained approximation method is not compatible with RD schemes (indeed, we have obtained some preliminary numerical results with a prototype of a RD scheme based upon the constrained approximation method, in the case of scalar equations with Lagrange elements and an explicit scheme). But the questions we have raised should first be carefully assessed before trying to make an RD scheme compatible with the adaptive method of [8]. To the contrary, the p-adaptive method for RD schemes we have proposed here is designed for RD schemes by construction and as such avoids all the problems evoked above.