Home Anisotropic Adaptive Finite Elements for a p-Laplacian Problem
Article Open Access

Anisotropic Adaptive Finite Elements for a p-Laplacian Problem

  • Paride Passelli EMAIL logo and Marco Picasso ORCID logo
Published/Copyright: June 26, 2024

Abstract

The p-Laplacian problem - ( ( μ + | u | p - 2 ) u ) = f is considered, where μ is a given positive number. An anisotropic a posteriori residual-based error estimator is presented. The error estimator is shown to be equivalent, up to higher order terms, to the error in a quasi-norm. The involved constants being independent of μ, the solution, the mesh size and aspect ratio. An adaptive algorithm is proposed and numerical results are presented when p = 3 . From this model problem, we propose a simplified error estimator and use it in the framework of an industrial application, namely a nonlinear Navier–Stokes problem arising from aluminium electrolysis.

MSC 2020: 65N30; 65N50

1 Introduction

Adaptive meshes demonstrated their efficiency to solve PDEs at a reduced computational cost, for a given level of accuracy. When boundary or internal layers are involved, anisotropic meshes, that are meshes with large aspect ratio, should be considered [24, 3]. The adaptive criteria is, when possible, based on an a posteriori error estimates, see for instance [1, 4, 6, 18, 40, 47, 49] in the isotropic framework and [9, 10, 16, 21, 22, 23, 27, 35, 36, 39, 42, 43, 50] when anisotropic meshes are considered.

A posteriori error estimates for the p-Laplacian problem - ( ( μ + | u | p - 2 ) u ) = f were proposed in [14, 7, 32, 13, 8, 17] for isotropic finite elements. A quasi-norm was introduced to obtain efficient and reliable error estimates. In [8, 17] convergence of an adaptive finite elements method for the p-Laplace equation is proven, when isotropic meshes and recursive bisection refinement are considered.

In this paper the model problem is studied with anisotropic finite elements. Since the anisotropic meshes are not nested, the framework of [8, 17] cannot be used. A residual-based anisotropic error estimator is derived for the p-Laplacian problem and an anisotropic adaptive algorithm is presented. We show that the error estimator is equivalent to the error in a quasi-norm, up to higher order terms. The involved constants are independent of μ, the solution, the mesh size and aspect ratio. We derive an upper bound for the classical W 0 1 , p norm and a quasi-norm. A lower bound for the error in a quasi-norm is proved. Numerical experiments confirm that when μ is small, the error estimator is sharp with respect to the error in the quasi-norm only. When μ is large, the error estimator is sharp with respect to both the W 0 1 , 3 and the quasi-norm.

The outline is the following: in Section 2 the problem and the numerical method are presented. In Section 3, we introduce the anisotropic finite element setting and Clément interpolation estimates in L p norms. In Section 4, we introduce the new error estimator, we prove an upper bound for the quasi-norm and the W 0 1 , p norm and a lower bound for the quasi-norm only. In Section 5, we present numerical experiments with non-adapted meshes. In Section 6, the adaptive algorithm and numerical experiments with adapted anisotropic meshes are presented and confirm the theoretical findings. In Section 7, an application to the Navier–Stokes equation with the so-called Smagorinsky turbulence model [45] in the context of aluminium electrolysis is presented. This industrial problem has some multi-scale features (from meters to millimeters), thus anisotropic finite elements are particularly efficient.

2 Problem Statement and Numerical Method

Consider a polygonal domain Ω d for d = 2 , 3 and a function u : Ω solution of the following problem

(2.1) { - ( ( μ + | u | p - 2 ) u ) = f in  Ω , u = 0 on  Ω ,

where μ 0 , p d and | | denotes the Euclidian norm in d . We assume throughout this paper that f L 2 ( Ω ) . The weak formulation of problem (2.1) consists in finding u W 0 1 , p ( Ω ) such that

(2.2) Ω ( μ + | u | p - 2 ) u v = Ω f v for all  v W 0 1 , p ( Ω ) .

When μ = 0 , existence and uniqueness of a solution u W 0 1 , p ( Ω ) have been proved in [31, 25]. The case μ > 0 can be proved in a similar way.

To solve problem (2.2), we use classical continuous piecewise linear finite elements. Consider, for any 0 < h < 1 , a conforming triangulation 𝒯 h of Ω ¯ having triangles K of diameter h K h . Let V h W 0 1 , p ( Ω ) be the classical finite element space of continuous piecewise linear functions on triangles K 𝒯 h with zero value on the boundary Ω . To approximate the solution of (2.2), we are looking for u h V h such that

(2.3) Ω ( μ + | u h | p - 2 ) u h v h = Ω f v h for all  v h V h .

Existence and uniqueness of u h V h can be again proved similarly as in [31, 25]. To solve this nonlinear problem, Newton’s method is advocated.

3 Anisotropic Finite Elements and Clément’s Interpolation Estimates

Anisotropic finite elements, that is meshes with possibly large aspect ratio, are considered. The framework of [22, 23] is used. For any K 𝒯 h , we note T K : K ^ K one of the affine transformation mapping the reference triangle K ^ into K defined by

(3.1) 𝐱 = T K ( 𝐱 ^ ) = M K 𝐱 ^ + 𝐭 K ,

where M K d × d and 𝐭 K d . Since M K is invertible, it admits a singular value decomposition M K = R K T Λ K P K , where R K and P K are orthogonal matrices and Λ K is a diagonal matrix

Λ K = ( λ 1 , K λ d , K ) , λ 1 , K λ d , K > 0 ,

and

R K = ( 𝐫 1 , K T 𝐫 d , K T ) .

In the above notations, for a given K 𝒯 h , 𝐫 1 , K , , 𝐫 d , K are the unit vectors corresponding to the stretching directions and λ 1 , K , , λ d , K their respective amplitude. In the context of anisotropic meshes, the classical minimal angle condition is not necessary anymore. However some additional assumptions are necessary. From now on, we will assume the two following geometrical assumptions:

  1. For each K, the cardinality of Δ K , that is the union of triangles sharing a vertex with K, is uniformly bounded from above, independently of the mesh geometry.

  2. For each K, the diameter of Δ K ^ = T K - 1 ( Δ K ) is uniformly bounded from above, independently of the mesh geometry.

Both assumptions guarantee that the constants involved in the interpolation estimates are mesh aspect ratio independent, which is necessary when working with anisotropic meshes. In particular, the second assumption excludes too distorted reference patches, an example of acceptable and non-acceptable patch is presented in Figure 1. The two assumptions seem to be fulfilled when using the anisotropic mesh generators in our numerical experiment [22, 23, 38, 43]. The Clément’s interpolant [15] is an important tool for deriving residual based a posteriori error estimates in the H 0 1 norm. Note that when d = 2 and p > 2 , since W 1 , p ( Ω ) C 0 ( Ω ¯ ) , the Lagrangian interpolant could also be used. However, in Section 7, we use our results when d = 3 and p = 3 , thus we only consider Clément’s interpolant in this paper.

Figure 1 
               Top: example of acceptable patch: the size of 
                     
                        
                           
                              Δ
                              ⁢
                              
                                 K
                                 ^
                              
                           
                        
                        
                        {\Delta\hat{K}}
                     
                   is independent of the aspect ratio 
                     
                        
                           
                              H
                              /
                              h
                           
                        
                        
                        {H/h}
                     
                  . Bottom: example of non-acceptable patch where the size of 
                     
                        
                           
                              Δ
                              ⁢
                              
                                 K
                                 ^
                              
                           
                        
                        
                        {\Delta\hat{K}}
                     
                   depends of the aspect ratio 
                     
                        
                           
                              H
                              /
                              h
                           
                        
                        
                        {H/h}
                     
                  .
Figure 1

Top: example of acceptable patch: the size of Δ K ^ is independent of the aspect ratio H / h . Bottom: example of non-acceptable patch where the size of Δ K ^ depends of the aspect ratio H / h .

Let R h : W 1 , p ( Ω ) V h be Clément’s interpolant with p 2 . Proceeding as in [22, 23] for the case d = 2 and p = 2 , we obtain the following.

Proposition 1 (Anisotropic Clément Interpolation Error Estimate in the L p Norm).

For p 2 , there exists a constant C ^ > 0 depending only on the reference triangle K ^ and on p such that for all v W 1 , p ( Ω ) , for all K T h , we have for i = 1 , , d + 1 ,

(3.2) ( j = 1 d λ j , K | K i | ) 1 p v - R h ( v ) L p ( K i ) + v - R h ( v ) L p ( K ) C ^ ω p , K ( v ) ,

where

(3.3) ( ω p , K ( v ) ) p = i = 1 d λ i , K p v 𝐫 i , K L p ( Δ K ) p ,

Δ K is the union of triangles sharing a vertex with those of K and K i is an edge (face if d = 3 ) of K.

Proof.

Let v ^ : K ^ be defined by

v ^ ( 𝐱 ^ ) = v ( M K 𝐱 ^ + 𝐭 𝐤 ) .

Using the equivalence between the p-norm and the 2-norm in d and the orthogonality of P K there exists C ^ depending only on p such that

^ v ^ L p ( K ^ ) p = M K T v L p ( K ^ ) p = K ^ P K T Λ K R K v p p C ^ K ^ Λ K R K v p p = C ^ i = 1 d λ i , K p K ^ v 𝐫 i , K p p .

After the change of variable (3.1), since det M K = j = 1 d λ i , K , we obtain

(3.4) ^ v ^ L p ( K ^ ) p C ^ i = 1 d ( λ i , K p j = 1 d λ j , K v 𝐫 i , K L p ( K ) p ) .

Moreover,

v - R h ( v ) L p ( K ) p = ( j = 1 d λ j , K ) v ^ - R ^ h ( v ^ ) L p ( K ^ ) p
C ^ h K ^ p ( j = 1 d λ j , K ) ^ v ^ L p ( Δ K ^ ) p
= C ^ ( j = 1 d λ j , K ) T Δ K ^ ^ v ^ L p ( T ) p .

Using (3.4), we therefore obtain

v - R h ( v ) L p ( K ) p C ^ ( ω p , K ( v ) ) p .

Proceeding in a similar way we have for i = 1 , , d + 1 ,

v - R h ( v ) L p ( K i ) p C ^ | K i | v ^ - R ^ h ( v ^ ) L p ( K ^ i ) p C ^ | K i | ^ v ^ L p ( Δ K ^ ) p .

Using again (3.4), we obtain

v - R h ( v ) L p ( K i ) p | K i | j = 1 d λ j , K ( ω p , K ( v ) ) p ,

which concludes the proof. ∎

Remark 1.

When p = 2 , we have as in [22, 23]

( ω 2 , K ( v ) ) 2 = i = 1 d λ i , K 2 ( 𝐫 i , K T G K ( v ) 𝐫 i , K ) ,

with G K ( v ) d × d defined as ( G K ( v ) ) i j = Δ K v x i v x j for i , j = 1 , , d .

4 The Anisotropic Error Estimator

Let u W 0 1 , p ( Ω ) be the weak solution of (2.2) and u h V h the solution of the finite element approximation (2.3). We derive an estimator for the error in the norm

μ ( u - u h ) L 2 ( Ω ) 2 + ( u - u h ) L p ( Ω ) p

and also in the quasi-norm

Ω | ( u - u h ) | 2 ( μ + ( | ( u - u h ) | + | u | ) p - 2 ) .

We first focus on the quasi-norm.

Proposition 2.

Let μ 0 , for w W 0 1 , p ( Ω ) and v W 0 1 , p ( Ω ) let

v ( w ) p = Ω | v | 2 ( μ + ( | v | + | w | ) p - 2 ) .

Then ( w ) is a quasi-norm on W 0 1 , p ( Ω ) .

The proof is as in [7]. We now state and prove a slight modification of [34, Lemma 2.1], where we set α 1 = α 2 = 0 , δ = 0 and μ has been added.

Lemma 1.

Let μ 0 , p 2 , d = 2 , 3 and let k C 1 ( 0 , ) C [ 0 , [ be defined by k ( t ) = μ + t p - 2 . Then there exist two constants C , M > 0 such that for all μ 0 and for all x , y R d we have

(4.1) | k ( | 𝐱 | ) 𝐱 - k ( | 𝐲 | ) 𝐲 | C | 𝐱 - 𝐲 | ( μ + ( | 𝐱 | + | 𝐲 | ) p - 2 )

and

(4.2) ( k ( | 𝐱 | ) 𝐱 - k ( | 𝐲 | ) 𝐲 ) ( 𝐱 - 𝐲 ) M | 𝐱 - 𝐲 | 2 ( μ + ( | 𝐱 | + | 𝐲 | ) p - 2 ) .

Proof.

We start by proving (4.1). First, observe that for all s , t > 0 ,

(4.3) | k ( t ) t - k ( s ) s | C | t - s | ( μ + ( t + s ) p - 2 ) .

Indeed, assume for instance 0 s t ; by the mean value theorem there exists s x t such that

k ( t ) t - k ( s ) s = ( t - s ) ( ( p - 2 ) x p - 3 x + ( μ + x p - 2 ) )
( t - s ) ( p - 1 ) ( μ + x p - 2 )

and (4.3) follows directly. Proceeding as in [34, Lemma 2.1], we can obtain (4.1). To prove (4.2), we observe that there exists M > 0 for all t s 0 such that

k ( t ) t - k ( s ) s M ( t - s ) ( μ + ( t + s ) p - 2 ) .

Indeed, we have

( t - s ) ( μ + ( t + s ) p - 2 ) = ( t - s ) μ + t exp ( ( p - 2 ) ln ( t + s ) ) - s exp ( ( p - 2 ) ln ( t + s ) )
( t - s ) μ + 2 p - 2 ( t exp ( ( p - 2 ) ln ( t ) ) - s exp ( ( p - 2 ) ln ( s ) ) )
2 p - 2 ( t ( μ + t p - 2 ) - s ( μ + s p - 2 ) ) .

Again following the idea of [34, Lemma 2.1], we prove (4.2). ∎

We now recall [32, Lemma 5.2] when θ = 1 .

Lemma 2.

For all σ 1 , σ 2 , σ 3 0 we have

( σ 3 + σ 1 ) p - 2 σ 1 σ 2 ( σ 3 + σ 1 ) p - 2 σ 1 2 + ( σ 3 + σ 2 ) p - 2 σ 2 2 .

As in [33] we set

a ( u , v ) = Ω ( μ + | u | p - 2 ) u v ,

where it must be noticed that a ( , ) is not linear with respect to the first variable.

Following [32] and [7] we can prove the following result.

Proposition 3.

Let C and M be the constants of Lemma 1. Let u W 0 1 , p ( Ω ) . Then for all v , w W 0 1 , p ( Ω ) we have

(4.4) a ( u , u - v ) - a ( v , u - v ) M 2 u - v ( u ) p

and

(4.5) | a ( u , w ) - a ( v , w ) | C ( u - v ( u ) p + w ( u ) p ) .

Proof.

First we prove equation (4.4) using relation (4.2) of Lemma 1. We have

a ( u , u - v ) - a ( v , u - v ) = Ω ( ( μ + | u | p - 2 ) u - ( μ + | v | p - 2 ) v ) ( u - v )
M Ω | ( u - v ) | 2 ( μ + ( | u | + | v | ) p - 2 ) .

Using the fact that for all 𝐱 , 𝐲 d we have

2 ( | 𝐱 | + | 𝐲 | ) | 𝐱 | + | 𝐱 - 𝐲 | ,

we obtain

a ( u , u - v ) - a ( v , u - v ) M 2 Ω | ( u - v ) | 2 ( μ + ( | u | + | ( u - v ) | ) p - 2 )
= M 2 u - v ( u ) p .

We now prove (4.5). We have

| a ( u , w ) - a ( v , w ) | Ω | ( ( μ + | u | p - 2 ) u - ( μ + | v | p - 2 ) v ) w | .

We apply (4.1) of Lemma 1 and the fact that for all 𝐱 , 𝐲 d we have

| 𝐱 | + | 𝐲 | 2 | 𝐱 | + | 𝐱 - 𝐲 |

to obtain

| a ( u , w ) - a ( v , w ) | C Ω | ( u - v ) | ( μ + ( | u | + | ( u - v ) | ) p - 2 ) | w | .

Using Lemma 2, we obtain

| a ( u , w ) - a ( v , w ) | C ( Ω | ( u - v ) | 2 ( μ + ( | u | + | ( u - v ) | ) p - 2 ) + Ω ( μ + ( | u | + | w | ) p - 2 ) | w | 2 ) ,

which concludes the proof. ∎

The next result follows directly taking w = u - v in (4.4)–(4.5).

Proposition 4.

Let u W 0 1 , p ( Ω ) and v W 0 1 , p ( Ω ) . Then there exist two constants M , C > 0 independent of μ, u, v such that

(4.6) 1 2 C ( a ( u , u - v ) - a ( v , u - v ) ) u - v ( u ) p 2 M ( a ( u , u - v ) - a ( v , u - v ) ) .

We can also prove the following proposition.

Proposition 5.

There exists α > 0 such that for all μ 0 and u , v W 0 1 , p ( Ω ) we have

(4.7) α ( μ ( u - v ) L 2 ( Ω ) 2 + ( u - v ) L p ( Ω ) p ) a ( u , u - v ) - a ( v , u - v ) .

Proof.

It suffices to prove that there exists α such that for all 𝐱 , 𝐲 d we have

(4.8) ( ( μ + | 𝐱 | p - 2 ) 𝐱 - ( μ + | 𝐲 | p - 2 ) 𝐲 ) ( 𝐱 - 𝐲 ) α ( μ | 𝐱 - 𝐲 | 2 + | 𝐱 - 𝐲 | p ) .

Indeed, we have

( ( μ + | 𝐱 | p - 2 ) 𝐱 - ( μ + | 𝐲 | p - 2 ) 𝐲 ) ( 𝐱 - 𝐲 ) = μ | 𝐱 - 𝐲 | 2 + ( | 𝐱 | p - 2 𝐱 - | 𝐲 | p - 2 𝐲 ) ( 𝐱 - 𝐲 )

so that (4.8) follows using [25, Lemma 5.1]. ∎

We can now state the main result of the paper.

Theorem 1.

Let u W 0 1 , p ( Ω ) be a solution of (2.2), u h V h a solution of (2.3) and e = u - u h . Let p 2 . Then there exists C ^ 1 > 0 independent of μ, u, the mesh size and the aspect ratio such that

(4.9) μ ( u - u h ) L 2 ( Ω ) 2 + ( u - u h ) L p ( Ω ) p + u - u h ( u ) p C ^ 1 ( K 𝒯 h η 2 , K + ϵ 1 ) .

Moreover, if there exists a constant C ^ , dependent only on the reference triangle K ^ , such that for all K T h and for i = 1 , , d - 1 ,

(4.10) λ i , K 2 ( 𝐫 i , K T G K ( u - u h ) 𝐫 i , K ) C ^ λ d , K 2 ( 𝐫 d , K G K ( u - u h ) 𝐫 d , K )

and assuming λ i , K vary smoothly around K, then there exists C ^ 2 > 0 independent of μ , u , the mesh size and the aspect ratio such that

(4.11) K 𝒯 h η 2 , K C ^ 2 ( u - u h ( u ) p + j = 1 4 ϵ j ) .

Here

(4.12)

η 2 , K = ( ( ( μ + | u h | p - 2 ) u h ) + Π K f L 2 ( K )
+ 1 2 i = 1 d + 1 ( | K i | j = 1 d λ j , K ) 1 2 [ ( μ + | u h | p - 2 ) u h 𝐧 ] L 2 ( K i ) ) ω 2 , K ( u - u h )

and

ϵ 1 = K 𝒯 h f - Π K f L 2 ( K ) ω 2 , K ( e ) , ϵ 2 = K 𝒯 h λ d , K f - Π K f L 2 ( K ) 2 ,
ϵ 3 = K 𝒯 h K 𝒫 K K | ( μ + | u h | ) | K p - 2 - ( μ + | u h | ) | K p - 2 | | e | 2 , ϵ 4 = K 𝒯 h λ d , K e L 2 ( Δ K ) 2 .

Here ω 2 , K ( ) is defined in (3.3), n stands for the unit outer normal to triangle K, [ ] denotes the jump across the edges (or faces if d = 3 ) K i of K for i = 1 , , d + 1 ( [ ] = 0 if K Ω ). For all K T h , Π K f is the L 2 ( K ) projection of f onto the set of piecewise constant functions defined by

Π K f = 1 | K | K f .

Moreover, we denote

𝒫 K = Δ K ( i = 1 d + 1 Δ K i )

with K i the ith element sharing a facet K i with K.

Remark 2.

Let 2 p p and let q be the Hölder conjugate of p . Then the upper bound (4.9) can be generalized to

C ^ 1 ( K 𝒯 h η p , K + K 𝒯 h f - Π K f L q ( K ) ω p , K ( u - u h ) ) ,

where

η p , K = ( ( ( μ + | u h | p - 2 ) u h ) + Π K f L q ( K )
+ 1 2 i = 1 d + 1 ( | K i | j = 1 d λ j , K ) 1 p [ ( μ + | u h | p - 2 ) u h 𝐧 ] L q ( K i ) ) ω p , K ( u - u h ) .

Remark 3.

Assume f H 1 ( Ω ) . Then we have

f - Π K f L 2 ( K ) 2 C ^ i = 1 d λ i , K 2 f 𝐫 i , K L 2 ( K ) 2 .

In the isotropic case this yields ϵ 1 , ϵ 2 , ϵ 4 = O ( h 3 ) , that are high order terms compared to u - u h ( u ) p which should be, according to [7], O ( h 2 ) . In the anisotropic context if, for instance, f depends only on x 2 , d = 2 and for all K 𝒯 h , 𝐫 1 , K = ( 0 , 1 ) , then ϵ 1 , ϵ 2 , ϵ 4 = O ( ( max K 𝒯 h λ 2 , K ) 3 ) . Numerical experiments performed in Sections 5 and 6 confirm previous predictions and show that ϵ 3 is also of higher order with respect to u - u h ( u ) p , see Figures 2 and 3. Thus when max K 𝒯 h λ d , K 2 0 , then u - u h ( u ) p and K 𝒯 h η 2 , K are equivalent up to higher order terms.

Remark 4.

Observe that local error estimator (4.12) is not standard since the exact solution is present in the term ω 2 , K ( u - u h ) . In practice post-processing techniques can be applied in order to approximate the gradient of the exact solution. For instance Zienkiewicz–Zhu (ZZ) post-processing can be applied [2, 51, 52]. We will replace

( u - u h ) x i by Π h Z Z u h x i - u h x i , i = 1 , , d ,

where for any v h V h , Π h Z Z v h x i is an approximate L 2 ( Ω ) projection of v h x i onto V h . The use of ZZ post-processing can be justified theoretically whenever superconvergence occurs. For elliptic equations and structured meshes, superconvergence of the ZZ recovery is presented in [12, 52]. Then, the post-processing is asymptotically exact, that is,

lim h 0 Π h Z Z u h - u h L 2 ( Ω ) u - u h L 2 ( Ω ) = 1 .

In [28] equivalence between Π h Z Z u h - u h L 2 ( Ω ) and u - u h L 2 ( Ω ) on general meshes has been proven, while the superconvergence of the ZZ gradient recovery has been shown for unstructured anisotropic meshes in [11]. Moreover, the efficiency of the ZZ post-processing has already been numerically demonstrated in [42, 43, 41, 26, 37, 44, 10] in the case of anisotropic meshes for elliptic, parabolic and hyperbolic equations. Notice that post-processing techniques for the p-Laplacian problem have already been studied in [14, 33].

Proof of the Upper Bound (4.9) of Theorem 1.

For all v W 0 1 , p ( Ω ) and v h V h we have using (2.2) and (2.3),

a ( u , v ) - a ( u h , v ) = Ω f ( v - v h ) - Ω ( μ + | u h | p - 2 ) u h ( v - v h )
= K 𝒯 h ( K ( f + ( μ + | u h | p - 2 ) u h ) ( v - v h ) + 1 2 K [ ( μ + | u h | p - 2 ) u h 𝐧 ] ( v - v h ) ) ,

where we integrated by parts. We set v = u - u h and v h = R h ( u - u h ) , adding and subtracting Π K f , using the Cauchy–Schwarz inequality and Proposition 1, we obtain

| a ( u , v ) - a ( u h , v ) |
K 𝒯 h ( Π K f + ( ( μ + | u h | p - 2 ) u h ) L 2 ( K ) + f - Π K f L 2 ( K ) ) u - u h - R h ( u - u h ) L 2 ( K )
    + 1 2 K 𝒯 h [ ( μ + | u h | p - 2 ) u h 𝐧 ] L 2 ( K ) u - u h - R h ( u - u h ) L 2 ( K )
C K 𝒯 h ( Π K f + ( ( μ + | u h | p - 2 ) u h ) L 2 ( K ) + 1 2 i = 1 d + 1 ( | K i | j = 1 d λ j , K ) 1 2 [ ( μ + | u h | p - 2 ) u h 𝐧 ] L 2 ( K i )
       + f - Π K f L 2 ( K ) ) ω 2 , K ( u - u h ) ,

which, together with Propositions 4 and 5, yields (4.9). ∎

We now aim to prove the lower bound (4.11). The classical standard bubble functions [5, 48], adapted to the anisotropic case [43, 20] and modified to take into account the nonlinearity, are involved.

Proposition 6.

Let e = u - u h . Under the same assumptions of Theorem 1 there exist a function φ W 1 , p ( Ω ) and a constant C ^ > 0 (that depends only on the reference triangle K ^ and on p) such that for any K T h and i = 1 , , d + 1 ,

K i [ ( μ + | u h | p - 2 ) u h 𝐧 ] φ = 1 2 ( | K i | 1 2 ( j = 1 d λ j , K ) 1 2 ω 2 , K ( e ) + | K i | 1 2 ( j = 1 d λ j , K i ) 1 2 ω 2 , K i ( e ) )
(4.13) × [ ( μ + | u h | p - 2 ) u h 𝐧 ] L 2 ( K i ) ,
(4.14) K ( Π K f + ( ( μ + | u h | p - 2 ) u h ) ) φ = Π k f + ( ( μ + | u h | p - 2 ) u h ) L 2 ( K ) ω 2 , K ( e ) ,
(4.15) ( j = 1 d λ j , K 2 φ 𝐫 j , K L 2 ( K ) 2 ) 1 2 C ^ ( ω 2 , K ( e ) + i = 1 d + 1 ω 2 , K i ( e ) ) ,
(4.16) K | φ | 2 C ^ ( ω 2 , K 2 ( e ) λ d , K 2 + i = 1 d + 1 ω 2 , K i 2 ( e ) λ d , K i 2 ) ,
(4.17) K | φ | p Z C ^ e L p ( 𝒫 K p ) .

We denote for i = 1 , , d + 1 , K i the ith element sharing facet K i with K.

Proof.

Following the proof of [43], we claim that

φ = K 𝒯 h C K Ψ K + 1 2 K 𝒯 h i = 1 d + 1 C K i Ψ K i ,

where C K , C K i are constants Ψ K , Ψ K i are the usual bubble functions over K and its edges (faces if d = 3 ) K i , i = 1 , , d + 1 . We set C K i = 0 if K i Ω . First we compute C K i for i = 1 , , d + 1 . We require that the constants associated to the same facet shared by two elements are equal. Using the fact that Ψ K , Ψ K j are zero over K i for all i j , the fact that ( Π K f + ( ( μ + | u h | p - 2 ) u h ) and [ ( μ + | u h | p - 2 ) u h 𝐧 ] are constants over K and K i respectively, | K | = C ^ j = 1 d λ j , K and (4.14) and (6), we obtain

C K i = ± 1 2 ( ( | K i | 2 j = 1 d λ j , K ) 1 2 ω 2 , K ( e ) + ( | K i | 2 j = 1 d λ j , K i ) 1 2 ω 2 , K i ( e ) ) 1 K i Ψ K i

and

C K = ± 1 K Ψ K ( C ^ ( j = 1 d λ j , K ) 1 2 ω 2 , K ( e ) i = 1 d + 1 1 2 ( ( | K i | 2 j = 1 d λ j , K ) 1 2 ω 2 , K ( e ) + ( | K i | 2 j = 1 d λ j , K i ) 1 2 ω 2 , K i ( e ) ) K Ψ K i K i Ψ K i ) ,

where the signs are chosen again accordingly to the signs of [ Π K i μ u h 𝐧 ] and Π K f + ( ( μ + | u h | p - 2 ) u h ) ) . By construction φ satisfies thus both equations (6)–(4.14). We now prove (4.17). We have

K | φ | p C ^ K ( | C K | p | Ψ K | p + i = 1 d + 1 | C K i | p | Ψ K i | p ) .

Using the change of variables (3.1), we have

K Ψ K = | K | K ^ Ψ ^ K ^ = C ^ j = 1 d λ j , K ,
K Ψ K i = | K | K ^ Ψ ^ K ^ i = C ^ j = 1 d λ j , K , i = 1 , , d + 1 ,
K i Ψ K i = | K i | K ^ i Ψ ^ K ^ i = C ^ | K i | , i = 1 , , d + 1 .

Thus we obtain the following bounds:

| C K i | p C ^ ( ω 2 , K p ( e ) ( j = 1 d λ j , K ) p 2 + ω 2 , K i p ( e ) ( j = 1 d λ j , K i ) p 2 ) ,
| C K | p C ^ ( ω 2 , K p ( e ) ( j = 1 d λ j , K ) p 2 + i = 1 d + 1 ω 2 , K i p ( e ) ( j = 1 d λ j , K i ) p 2 ) ,

hence

K | φ | p C ^ ( ω 2 , K p ( e ) ( j = 1 d λ j , K ) p 2 + i = 1 d + 1 ω 2 , K i p ( e ) ( j = 1 d λ j , K i ) p 2 ) ( K | Ψ K | p + j = 1 d + 1 K | Ψ K j | p ) .

Using assumption (4.10), we have

ω 2 , K p ( e ) C ^ ( λ d , K 2 e L 2 ( Δ K ) 2 ) p 2 .

Proceeding in the same way as we did for (3.4), we can show that there exists a constant C ^ > 0 such that for all v W 0 1 , p ( Ω ) ,

(4.18) v L p ( K ) p C ^ i = 1 d ( j = 1 d λ j , K λ i , K p ) ^ v ^ 𝐩 i , K L p ( K ^ ) p .

Since λ 1 , K λ d , K > 0 , we obtain from (4.18) that

Ψ K L p ( K ) p C ^ j = 1 d λ j , K λ d , K p ^ Ψ ^ K ^ L p ( K ^ ) p

and the same occurs for Ψ K i . Altogether we obtain using Hölder’s inequality

K | φ | p C ^ ( j = 1 d λ j , K ) 2 - p 2 ( ( Δ K | e | 2 ) p 2 + i = 1 d + 1 ( Δ K i | e | 2 ) p 2 )
C ^ ( j = 1 d λ j , K ) 2 - p 2 ( | Δ K | p - 2 2 Δ K | e | p + i = 1 d + 1 | Δ K i | p - 2 2 Δ K i | e | p ) ,

and using assumption (1) of Section 3 and the fact that λ i , K vary smoothly around K, we obtain (4.17). To prove (4.15), we proceed in the same way to get

K | φ 𝐫 i , K | 2 C ^ ( ω 2 , K 2 ( e ) j = 1 d λ j , K + l = 1 d + 1 ω 2 , K l 2 ( e ) j = 1 d λ j , K l ) ( K | Ψ K 𝐫 i , K | 2 + j = 1 d + 1 K | Ψ K j 𝐫 i , K | 2 ) .

Using again (3.4) (which is an equality when p = 2 ), we obtain the desired result. Finally, we proceed in a similar way to prove (4.16) ∎

We are now ready to prove (4.11).

Proof of the Lower Bound (4.11) of Theorem 1.

Using the definition of η 2 , K and identities (6)–(4.14), one can write

K 𝒯 h η 2 , K = K 𝒯 h K ( Π K f + ( ( μ + | u h | p - 2 ) u h ) φ ) + 1 2 K [ ( μ + | u h | p - 2 ) u h 𝐧 ] φ ,

where φ is the bubble function introduced in proposition 6. Using (2.2) and integration by parts we obtain

K 𝒯 h η 2 , K = K 𝒯 h K ( ( μ + | u | p - 2 ) u - ( μ + | u h | p - 2 ) u h ) φ - K 𝒯 h K ( f - Π K f ) ( φ - Π K φ )
| a ( u , φ ) - a ( u h , φ ) | + K 𝒯 h K | f - Π K f | | φ - Π K φ | .

Using Proposition 3 together with the Cauchy–Schwarz inequality, we get

(4.19) K 𝒯 h η 2 , K C ^ ( e ( u h ) p + φ ( u h ) p ) + K 𝒯 h f - Π K f L 2 ( K ) φ - Π K φ L 2 ( K ) .

Using the change of variable (3.1) and the Poincaré–Wirtinger inequality, we can show that

φ - Π K φ L 2 ( K ) 2 = j = 1 d λ j , K K ^ ( φ ^ - Π ^ K ^ φ ^ ) 2 C ^ j = 1 d λ j , K K ^ | ^ φ ^ | 2 ,

which together with (3.4) and (4.15) gives

φ - Π K φ L 2 ( K ) C ^ ( ω 2 , K ( e ) + i = 1 d + 1 ω 2 , K i ( e ) ) .

We obtain thus, using Young’s inequality,

K 𝒯 h f - Π K f L 2 ( K ) φ - Π K φ L 2 ( K ) C ^ ( ϵ 1 + ϵ 2 + ϵ 4 ) .

We have

φ ( u h ) p = K 𝒯 h K ( μ + | u h | + | φ | ) p - 2 | φ | 2 K 𝒯 h C ^ K ( μ + | u h | ) p - 2 | φ | 2 + C ^ K 𝒯 h K | φ | p
C ^ K 𝒯 h 𝒫 K ( μ + | u h | p - 2 ) | K | e | 2 + | e | p ,

where we used (4.16)–(4.17) and (4.10). We have

C ^ K 𝒯 h 𝒫 K ( μ + | u h | ) | K p - 2 | e | 2 = C ^ K 𝒯 h 𝒫 K ( μ + | u h | ) p - 2 | e | 2
+ C ^ K 𝒯 h K 𝒫 K K ( ( μ + | u h | ) | K p - 2 - ( μ + | u h | ) | K p - 2 ) | e | 2 .

Altogether, we have

K 𝒯 h η 2 , K C ^ ( e ( u h ) p + ϵ 1 + ϵ 2 + ϵ 4 ) + C ^ K 𝒯 h K 𝒫 K K ( ( μ + | u h | ) | K p - 2 - ( μ + | u h | ) | K p - 2 ) | e | 2 .

Therefore, applying the triangle inequality e ( u h ) p C ^ e ( u ) p , we get the desired result. ∎

5 Numerical Experiments with Non-Adapted Meshes

The goal of this section is to check numerically the sharpness of the error estimator (4.12) derived in the previous section. We are interested by the case p = 3 , since it corresponds to the so called Smagorinsky model [46, 45] discussed in Section 7. We define

e Q N = u - u h ( u ) 3 , e 3 = ( u - u h ) L 3 ( Ω ) 3 and e 2 = μ ( u - u h ) L 2 ( Ω ) 2 .

We also define the following effectivity indices:

(5.1) ei N = K 𝒯 h η 2 , K e 2 + e 3 ,
(5.2) ei Q N = K 𝒯 h η 2 , K e Q N ,
ei Z Z = u h - Π h Z Z u h L 2 ( Ω ) ( u - u h ) L 2 ( Ω ) .

Let Ω = ( 0 , 1 ) 2 and consider problem (2.1). Results with various mesh sizes ( h 1 and h 2 being the mesh sizes in directions x 1 and x 2 ) and various values of μ are reported in Table 1 for the exact solution

(5.3) u ( x , y ) = 4 ( 1 - e - α x - ( 1 - e - α ) x ) y ( 1 - y )

as in [23], with α = 50 . Similarly, in Tables 2 and 3 results are reported for the exact solution

(5.4) u ( x , y ) = tanh ( x - 0.5 ϵ )

with ϵ = 0.1 and ϵ = 0.05 respectively. Note that u defined by (5.4) is not zero on the boundary, thus the theory does not apply as is anymore. Nevertheless, numerical results indicate that the error estimator is still sharp.

We observe that the error estimator is accurate for both the W 0 1 , 3 norm and the quasi-norm when μ is large. When μ is large and h is small enough the effectivity indices are around 9 which corresponds to the linear Laplace problem [20, 43]. When μ is small, only the quasi-norm should be considered. Indeed, when h is small enough, the corresponding effectivity indices remain between 11 and 18 uniformly with respect to μ and h. We can therefore conclude that the upper bound (4.9) seems to be sharp with p = 3 . We now discuss the higher order terms of the upper and lower bounds, ϵ 1 , ϵ 2 , ϵ 3 and ϵ 4 (see Theorem 1). Consider again the first numerical experiment when μ = 0 and u is given by (5.3) with α = 50 , let h 1 = 1 N x and h 2 = 10 h 1 . In Figure 2 left we can observe the order of convergence of the error estimator K 𝒯 h η 2 , K = O ( h 2 ) , the quasi-norm error e Q N = O ( h 2 ) and the W 0 1 , 3 error e 3 = O ( h 3 ) . Moreover, in Figure 2 right, we observe ϵ 1 , ϵ 2 , ϵ 4 = O ( h 3 ) and ϵ 3 = O ( h 5 ) , which are higher order terms with respect to e Q N as discussed in Remark 3.

Figure 2 
               Left: error estimator, quasi-norm error, 
                     
                        
                           
                              W
                              0
                              
                                 1
                                 ,
                                 3
                              
                           
                        
                        
                        {W^{1,3}_{0}}
                     
                   norm error with corresponding slopes, when u is given by (5.3), 
                     
                        
                           
                              μ
                              =
                              0
                           
                        
                        
                        {\mu=0}
                     
                   and 
                     
                        
                           
                              α
                              =
                              50
                           
                        
                        
                        {\alpha=50}
                     
                  . Right: Higher order terms 
                     
                        
                           
                              ϵ
                              1
                           
                        
                        
                        {\epsilon_{1}}
                     
                  , 
                     
                        
                           
                              ϵ
                              2
                           
                        
                        
                        {\epsilon_{2}}
                     
                  , 
                     
                        
                           
                              ϵ
                              3
                           
                        
                        
                        {\epsilon_{3}}
                     
                  , 
                     
                        
                           
                              ϵ
                              4
                           
                        
                        
                        {\epsilon_{4}}
                     
                   with corresponding slopes, when u is given by (5.3), 
                     
                        
                           
                              μ
                              =
                              0
                           
                        
                        
                        {\mu=0}
                     
                   and 
                     
                        
                           
                              α
                              =
                              50
                           
                        
                        
                        {\alpha=50}
                     
                  .
Figure 2 
               Left: error estimator, quasi-norm error, 
                     
                        
                           
                              W
                              0
                              
                                 1
                                 ,
                                 3
                              
                           
                        
                        
                        {W^{1,3}_{0}}
                     
                   norm error with corresponding slopes, when u is given by (5.3), 
                     
                        
                           
                              μ
                              =
                              0
                           
                        
                        
                        {\mu=0}
                     
                   and 
                     
                        
                           
                              α
                              =
                              50
                           
                        
                        
                        {\alpha=50}
                     
                  . Right: Higher order terms 
                     
                        
                           
                              ϵ
                              1
                           
                        
                        
                        {\epsilon_{1}}
                     
                  , 
                     
                        
                           
                              ϵ
                              2
                           
                        
                        
                        {\epsilon_{2}}
                     
                  , 
                     
                        
                           
                              ϵ
                              3
                           
                        
                        
                        {\epsilon_{3}}
                     
                  , 
                     
                        
                           
                              ϵ
                              4
                           
                        
                        
                        {\epsilon_{4}}
                     
                   with corresponding slopes, when u is given by (5.3), 
                     
                        
                           
                              μ
                              =
                              0
                           
                        
                        
                        {\mu=0}
                     
                   and 
                     
                        
                           
                              α
                              =
                              50
                           
                        
                        
                        {\alpha=50}
                     
                  .
Figure 2

Left: error estimator, quasi-norm error, W 0 1 , 3 norm error with corresponding slopes, when u is given by (5.3), μ = 0 and α = 50 . Right: Higher order terms ϵ 1 , ϵ 2 , ϵ 3 , ϵ 4 with corresponding slopes, when u is given by (5.3), μ = 0 and α = 50 .

Table 1

True errors and effectivity indices for various non-adapted meshes and various choices of μ when solving problem (2.1) with f such that the exact solution is given by (5.3), with α = 50 .

h 1 h 2 ei N ei Q N e Q N e 3 e 2 e i Z Z
μ = 0 0.025 0.025 8.43 2.71 81.50 26.21 0 0.59
0.0125 0.0125 20.08 6.14 15.32 3.11 0 0.82
0.00625 0.00625 85.17 9.66 3.12 3.53e - 01 0 0.93
0.003125 0.003125 208.47 11.99 7.11e - 01 4.07e - 02 0 0.98
0.01 0.1 55.12 10.07 7.75 1.40 0 0.96
0.005 0.05 111.62 12.10 1.86 2.02e - 01 0 0.98
0.0025 0.025 235.82 13.89 4.60e - 01 2.70e - 02 0 0.99
0.00125 0.0125 510.447 15.33 1.02e - 01 3.05e - 03 0 0.99
μ = 1 0.025 0.025 8.13 2.72 82.71 26.06 1.64 0.59
0.0125 0.0125 27.28 6.13 15.70 3.12 4.02e - 01 0.82
0.00625 0.00625 68.39 9.57 3.21 3.54e - 01 9.60e - 02 0.93
0.003125 0.003125 134.81 11.84 7.34e - 01 4.09e - 02 2.36e - 02 0.98
0.01 0.1 44.00 9.94 8.17 1.40 4.42e - 01 0.97
0.005 0.05 73.06 11.83 1.98 2.01e - 01 1.20e - 01 0.98
0.0025 0.025 111.91 13.48 4.92e - 01 2.70e - 02 3.22e - 02 0.99
0.00125 0.0125 148.42 14.78 1.10e - 01 3.07e - 03 7.91e - 03 0.99
μ = 100 0.025 0.025 4.16 3.23 222.81 23.33 149.50 0.63
0.0125 0.0125 7.22 5.65 53.55 3.00 38.89 0.84
0.00625 0.00625 9.40 7.36 12.54 3.51e - 01 9.47 0.94
0.003125 0.003125 10.69 8.36 3.05 4.09e - 02 2.35 0.98
0.01 0.1 8.66 7.60 51.04 1.37 43.42 0.97
0.005 0.05 9.08 7.99 13.72 1.99e - 01 11.87 0.98
0.0025 0.025 9.45 8.34 3.67 2.69e - 02 3.21 0.99
0.00125 0.0125 9.66 8.58 8.92e - 01 3.06e - 03 7.90e - 01 0.99
Table 2

True errors and effectivity indices for various non-adapted meshes and various choices of μ when solving problem (2.1) with f such that the exact solution is given by (5.4), with ϵ = 0.1 .

h 1 h 2 ei N ei Q N e Q N e 3 e 2 e i Z Z
μ = 0 0.1 0.1 36.63 8.08 8.13 1.79 0 1.15
0.05 0.05 85.12 11.59 1.30 1.77e - 01 0 1.04
0.025 0.025 159.87 13.17 2.73e - 01 2.25e - 02 0 1.02
0.0125 0.0125 309.63 13.82 6.11e - 02 2.73e - 03 0 1.00
0.00625 0.00625 643.61 14.22 1.51e - 02 3.33e - 04 0 1.00
0.003125 0.003125 1280.79 14.34 3.74e - 03 4.19e - 05 0 1.00
0.01 0.1 378.17 15.30 3.61e - 02 1.46e - 03 0 1.00
0.005 0.05 772.78 15.81 8.26e - 03 1.69e - 04 0 0.99
0.0025 0.025 1547.62 15.57 1.92e - 03 1.94e - 05 0 1.00
0.00125 0.0125 3098.95 15.62 4.79e.04 2.41e - 06 0 1.00
0.002 0.2 2574.62 18.56 1.25e - 03 9.02e - 06 0 1.00
0.001 0.1 4704.70 17.58 3.16e.04 1.18e - 06 0 0.99
0.0005 0.05 8910.96 17.36 7.57e - 05 1.47e - 07 0 0.99
μ = 1 0.1 0.1 27.07 8.15 8.76 1.70 9.37e - 01 1.16
0.05 0.05 45.73 11.27 1.47 1.73e - 01 1.88e - 01 1.06
0.025 0.025 59.82 12.46 3.15e - 01 2.24e - 02 4.33e - 02 1.02
0.0125 0.0125 71.31 12.95 7.12e - 02 2.73e - 03 1.02e - 02 1.01
0.00625 0.00625 80.49 13.29 1.76e - 02 3.33e - 04 2.57e - 03 1.00
0.003125 0.003125 85.29 13.39 4.39e - 03 4.10e - 05 6.47e - 04 1.00
0.01 0.1 79.92 14.25 4.20e - 02 1.46e - 03 6.04e - 03 1.00
0.005 0.05 88.81 14.67 9.68e - 02 1.69e - 04 1.43e - 03 1.00
0.0025 0.025 91.69 14.42 2.26e - 03 1.94e - 05 3.36e - 04 1.00
0.00125 0.0125 95.00 14.46 5.61e - 04 2.41e - 06 8.30e - 05 1.00
0.002 0.2 111.81 17.23 1.47e - 03 9.03e - 06 2.17e - 04 1.00
0.001 0.1 108.85 16.32 3.70e - 04 1.18e - 06 5.44e - 05 0.99
0.0005 0.05 107.82 16.09 8.88e - 05 1.47e - 07 1.31e - 05 0.99
μ = 100 0.1 0.1 8.28 7.80 88.35 1.36 81.91 1.17
0.05 0.05 8.82 8.33 18.90 1.63e - 01 17.70 1.08
0.025 0.025 8.66 8.19 4.56 2.22e - 02 4.29 1.03
0.0125 0.0125 8.55 8.09 1.08 2.73e - 03 1.02 1.01
0.00625 0.00625 8.65 8.19 2.73e - 01 3.34e - 04 2.58e - 01 1.00
0.003125 0.003125 8.64 8.17 6.83e - 02 4.21e - 05 6.46e - 02 1.00
0.01 0.1 8.88 8.40 6.37e - 01 1.46e - 03 6.00e - 01 1.01
0.005 0.05 8.92 8.44 1.51e - 01 1.79e - 04 1.42e - 01 1.00
0.0025 0.025 8.69 8.22 3.55e - 02 1.94e - 05 3.36e - 02 1.00
0.00125 0.0125 8.68 8.21 8.78e - 03 2.42e - 06 8,30e - 03 1.00
0.001 0.1 10.07 9.52 5.74e - 03 1.17e.06 5.42e - 03 0.99
0.0005 0.05 9.75 9.22 1.38e - 03 1.46e - 07 1.31e - 03 0.99
Table 3

True errors and effectivity indices for various non-adapted meshes and various choices of μ when solving problem (2.1) with f such that the exact solution is given by (5.4), with ϵ = 0.05 .

h 1 h 2 ei N ei Q N e Q N e L 3 e L 2 e i Z Z
μ = 0 0.05 0.05 14.42 4.66 53.93 17.41 0 0.70
0.025 0.025 74.61 11.67 4.93 7.70e - 01 0 1.06
0.0125 0.0125 153.99 13.08 1.05 8.89e - 02 0 1.00
0.00625 0.00625 324.70 13.90 2.50e - 01 1.07e - 02 0 1.01
0.003125 0.003125 650.11 14.17 6.12e - 02 1.33e - 03 0 1.00
0.01 0.1 198.78 15.46 6.62e - 01 5.15e - 02 0 0.90
0.005 0.05 384.54 14.91 1.38e - 01 5.36e - 03 0 0.89
0.0025 0.025 766.47 15.20 3.13e - 02 6.22e - 04 0 0.99
0.00125 0.0125 1544.29 15.53 7.71e - 03 7.76e - 05 0 1.00
0.002 0.2 1184.92 17.28 2.01e - 02 2.93e - 04 0 1.00
0.001 0.1 2189.84 17.27 5.26e - 03 4.14e - 05 0 1.00
0.0005 0.05 4364.82 17.12 2.21e - 03 4.77e - 06 0 1.00
μ = 1 0.05 0.05 14.14 5.03 51.76 14.25 4.19 0.76
0.025 0.025 53.56 11.46 5.25 7.64e - 01 3.58e - 01 1.08
0.0125 0.0125 83.35 12.73 1.13 8.85e - 02 8.33e - 02 1.02
0.00625 0.00625 115.70 13.44 2.70e - 01 1.07e - 02 2.07e - 02 1.01
0.003125 0.003125 138.78 13.68 6.64e - 02 1.33e - 03 5.21e - 03 1.00
0.01 0.1 97.92 14.94 7.01e - 01 4.96e - 02 5.74e - 02 0.95
0.005 0.05 121.83 14.53 1.48e - 01 5.23e - 03 1.24e - 02 0.95
0.0025 0.025 150.69 14.62 3.40e - 02 6,22e - 04 2.68e - 03 1.00
0.00125 0.0125 167.63 14.92 8.39e - 03 7.76e - 05 6.68e - 04 1.00
0.002 0.2 179.26 16.61 2.18e - 02 2.94e - 04 1.73e - 03 1.00
0.001 0.1 191.54 16.58 5.71e - 03 4.14e - 05 4.53e - 05 1.00
0.0005 0.05 196.31 16.44 1.32e - 03 4.77e - 06 1.06e - 05 1.00
μ = 100 0.05 0.05 8.22 7.36 185.30 6.02 159.75 1.16
0.025 0.025 9.29 8.34 38.88 7.36e - 01 34.19 1.09
0.0125 0.0125 9.35 8.39 9.16 8.78e - 02 8.14 1.03
0.00625 0.00625 9.47 8.50 2.31 1.07e - 02 2.06 1.01
0.003125 0.003125 9.47 8.50 5.80e - 01 1.34e - 03 5.19e - 01 1.00
0.01 0.1 9.83 8.83 5.49 4.72e - 02 4.88 1.02
0.005 0.05 9.76 8.77 1.24 5.16e - 03 1.11 1.00
0.0025 0.025 9.51 8.53 2.98e - 01 6.22e - 04 2.67e - 01 1.00
0.00125 0.0125 9.58 8.60 7.45e - 02 7.77e - 05 6.68e - 02 1.00
0.002 0.2 10.91 9.79 1.93e - 01 2.93e - 04 1.73e - 01 1.00
0.001 0.1 10.72 9.61 5.05e - 02 4.12e - 06 4.52e - 02 1.00
0.0005 0.05 10.58 9.49 1.18e - 02 4.75e - 05 1.09e - 01 0.99

In the next section, we introduce an adaptive algorithm aiming to control the relative error in the quasi norm ( u ) considering error estimator presented in Theorem 1 when p = 3 .

6 Adaptive Algorithm and Numerical Experiments with Adapted Meshes

The adaptive algorithm is similar to the one presented in [19, 43, 41, 20]. The goal of the adaptive algorithm is to build a sequence of meshes with possibley large aspect ratio, such that the relative estimated error is closed to a prescribed tolerance TOL, i.e.

(6.1) 0.75 TOL ( K 𝒯 h η 2 , K u h ( u ) 3 ) 1 2 1.25 TOL .

In the term

u h ( u ) 3 = Ω | u h | 2 ( μ + | u | + | u h | ) ,

u is in practice replaced by its approximate L 2 ( Ω ) projection onto V h , as discussed in Remark 4. Condition (6.1) holds whenever, for all K 𝒯 h ,

(6.2) 1 N K 0.75 2 TOL 2 u h ( u ) 3 η 2 , K 1 N K 1.25 2 TOL 2 u h ( u ) 3 ,

where N K is the total number of elements of 𝒯 h .

Assumption (4.10) suggests that the error estimator should be equidistributed in all stretching directions. We recall that λ i , K is the length of the stretching direction 𝐫 i , K of the triangle K and ω 2 , K ( v ) and G K ( v ) are defined in Remark 1. We define for all v W 0 1 , 2 ( Ω ) and for i = 1 , , d ,

(6.3) ( ω 2 , K , i ( v ) ) 2 = λ i , K 2 𝐫 i , K T G K ( v ) 𝐫 i , K

and the local error estimator in direction x i

(6.4) η 2 , K , i = ( ( ( μ + | u h | ) u h ) + Π K f L 2 ( K ) + 1 2 l = 1 d + 1 ( | K l | j = 1 d λ j , K ) 1 2 [ ( μ + | u h | ) u h 𝐧 ] L 2 ( K l ) ) ω 2 , K , i ( u - u h ) ,

so that

ω 2 , K ( v ) = ( i = 1 d ( ω 2 , K , i ( v ) ) 2 ) 1 2

and

η 2 , K = ( i = 1 d η 2 , K , i 2 ) 1 2 .

Thus in order to satisfies (6.2), we require that for all K 𝒯 h and for i = 1 , , d ,

(6.5) 1 d N K 2 0.75 4 TOL 4 u h ( u ) 6 η 2 , K , i 2 1 d N K 2 1.25 4 TOL 4 u h ( u ) 6 .

In this way we will equidistribute the local error estimator in all directions 𝐫 1 , K , , 𝐫 d , K . The adaptive strategy consists in aligning each triangle K (i.e. its directions of anisotropy 𝐫 1 , K , , 𝐫 d , K ) with the eigenvectors of the matrix G K and to define a new mesh size based on (6.5), details can be found in [19, 43, 41, 20] (we use the mesh generator BL2D [30] in two dimensions and the 3D Precise Mesh software [53] in three dimensions).

Let Ω = ( 0 , 1 ) 2 and consider again problem (2.1) with exact solution (5.3) with α = 100 . We apply the adaptive algorithm with a starting mesh of size h = 0.1 . Starting with a tolerance TOL = 0.5 , 40 iterations of the adaptive algorithm are performed. The tolerance is then halved and the process is repeated. We report in Table 4 the results obtained for different values of μ. We observe that, for any choice of μ, the value of ei Q N does not increase together with the average aspect ratio av ar and the maximum aspect ratio max ar . For μ = 100 , we report in Figure 3 the graph of convergence with respect to TOL = 0.5 × 2 - N of the errors (left) and the higher order terms ϵ 1 , ϵ 2 , ϵ 3 and ϵ 4 (right). Theoretical predictions of Remark 3 are again numerically verified.

Table 4

True errors, effectivity indices, number of vertices and aspect ratio for different values of tolerance TOL and μ, when solving problem (2.1) with f such that the exact solution is given by (5.3) with α = 100 .

TOL ei N ei Q N e Q N ei Z Z N v av ar max ar
μ = 0 0.5 33.52 10.25 20.50 0.97 89 25 116
0.25 70.84 12.97 3.48 0.98 249 32 272
0.125 98.70 13.67 8.17e - 01 0.97 838 33 239
0.0625 338.71 17.13 1.61e - 01 1.01 2975 36 421
0.03125 756.82 18.82 3.63e - 02 1.01 11582 37 825
0.015625 1539.65 18.33 9.43e - 03 1.00 44454 39 962
μ = 1 0.5 38.02 11.01 14.23 0.99 95 26 131
0.25 61.19 13.65 3.38 0.99 254 29 305
0.125 83.56 15.21 7.28e - 01 0.99 863 33 257
0.0625 112.37 16.89 1.58e - 01 1.01 3182 36 607
0.03125 135.67 17.64 3.93e - 02 1.00 11936 37 695
0.015625 150.25 17.90 9.84e - 03 1.00 46032 41 670
μ = 100 0.5 9.94 8.37 40.89 0.93 171 18 202
0.25 10.39 8.91 10.16 0.95 502 23 225
0.125 10.65 9.18 2.42 0.95 1713 26 253
0.0625 11.85 10.15 5.23e - 01 0.99 6649 26 436
0.03125 12.44 10.63 1.26e - 02 1.00 25536 27 635
0.015625 13.07 11.13 3.02e - 02 1.00 98424 29 1064
Figure 3 
               Left: error estimator, quasi-norm error, 
                     
                        
                           
                              W
                              0
                              
                                 1
                                 ,
                                 3
                              
                           
                        
                        
                        {W^{1,3}_{0}}
                     
                   and 
                     
                        
                           
                              W
                              0
                              
                                 1
                                 ,
                                 2
                              
                           
                        
                        
                        {W^{1,2}_{0}}
                     
                   norm error with corresponding slope, when u is given by (5.3), 
                     
                        
                           
                              μ
                              =
                              100
                           
                        
                        
                        {\mu=100}
                     
                   and 
                     
                        
                           
                              α
                              =
                              100
                           
                        
                        
                        {\alpha=100}
                     
                  . Right: Higher order terms 
                     
                        
                           
                              ϵ
                              1
                           
                        
                        
                        {\epsilon_{1}}
                     
                  , 
                     
                        
                           
                              ϵ
                              2
                           
                        
                        
                        {\epsilon_{2}}
                     
                  , 
                     
                        
                           
                              ϵ
                              3
                           
                        
                        
                        {\epsilon_{3}}
                     
                  , 
                     
                        
                           
                              ϵ
                              4
                           
                        
                        
                        {\epsilon_{4}}
                     
                   with corresponding slopes, when u is given by (5.3), 
                     
                        
                           
                              μ
                              =
                              100
                           
                        
                        
                        {\mu=100}
                     
                   and 
                     
                        
                           
                              α
                              =
                              100
                           
                        
                        
                        {\alpha=100}
                     
                  .
Figure 3 
               Left: error estimator, quasi-norm error, 
                     
                        
                           
                              W
                              0
                              
                                 1
                                 ,
                                 3
                              
                           
                        
                        
                        {W^{1,3}_{0}}
                     
                   and 
                     
                        
                           
                              W
                              0
                              
                                 1
                                 ,
                                 2
                              
                           
                        
                        
                        {W^{1,2}_{0}}
                     
                   norm error with corresponding slope, when u is given by (5.3), 
                     
                        
                           
                              μ
                              =
                              100
                           
                        
                        
                        {\mu=100}
                     
                   and 
                     
                        
                           
                              α
                              =
                              100
                           
                        
                        
                        {\alpha=100}
                     
                  . Right: Higher order terms 
                     
                        
                           
                              ϵ
                              1
                           
                        
                        
                        {\epsilon_{1}}
                     
                  , 
                     
                        
                           
                              ϵ
                              2
                           
                        
                        
                        {\epsilon_{2}}
                     
                  , 
                     
                        
                           
                              ϵ
                              3
                           
                        
                        
                        {\epsilon_{3}}
                     
                  , 
                     
                        
                           
                              ϵ
                              4
                           
                        
                        
                        {\epsilon_{4}}
                     
                   with corresponding slopes, when u is given by (5.3), 
                     
                        
                           
                              μ
                              =
                              100
                           
                        
                        
                        {\mu=100}
                     
                   and 
                     
                        
                           
                              α
                              =
                              100
                           
                        
                        
                        {\alpha=100}
                     
                  .
Figure 3

Left: error estimator, quasi-norm error, W 0 1 , 3 and W 0 1 , 2 norm error with corresponding slope, when u is given by (5.3), μ = 100 and α = 100 . Right: Higher order terms ϵ 1 , ϵ 2 , ϵ 3 , ϵ 4 with corresponding slopes, when u is given by (5.3), μ = 100 and α = 100 .

Consider again problem (2.1) with exact solution (5.4) with ϵ = 0.05 and Ω = ( 0 , 1 ) 2 . In Table 5 we can observe results obtained applying the adaptive algorithm with a starting mesh of size h = 0.1 . Starting with a tolerance TOL = 0.5 , 40 iterations of the adaptive algorithm are performed. The tolerance is then halved and the process is repeated. Results are again reported for different choices of μ.

The obtained results confirm that the effectivity index ei Q N is independent from the aspect ratio, but slightly depends on μ. In fact, when μ varies between 0 and 100 and the average aspect ratio varies between 1 and 15000, then the effectivity index ei Q N is between 11 and 24 for TOL small enough. In Figure 4 left, an example of obtained mesh is reported.

Table 5

True errors, effectivity indices, number of vertices and aspect ratio for different values of tolerance TOL and μ, when solving problem (2.1) with f given such that the exact solution is given by (5.4) with ϵ = 0.05 .

TOL ei N ei Q N e Q N ei Z Z N v av ar max ar
μ = 0 0.5 58.15 11.68 3.34 0.91 52 20 61
0.25 139.26 18.36 6.68e - 01 0.97 70 47 142
0.125 225.39 17.53 1.64e - 01 0.96 150 88 325
0.0625 412.32 20.95 4.00e - 02 0.99 246 189 865
0.03125 558.11 20.95 9.93e - 03 0.97 492 351 1365
0.015625 1282.96 22.44 2.34e - 03 0.98 863 845 4025
0.0078125 2426.23 23.76 5.73e - 04 0.99 1703 1661 6348
0.00390625 3384.70 22.91 1.45e - 04 0.98 3434 3286 16889
0.00195312 7852.33 23.20 3.61e - 05 0.99 6821 6640 38352
0.000976562 14905.10 23.23 9.12e - 06 0.99 13510 13250 61063
μ = 1 0.5 51.02 12.56 2.88 0.93 80 14 46
0.25 72.59 15.53 6.46e - 01 0.96 76 51 171
0.125 104.87 19.87 1.50-01 0.99 137 104 359
0.0625 95.51 19.04 4.37e - 02 0.98 272 166 575
0.03125 108.03 20.36 1.11e - 02 0.99 488 367 1444
0.015625 107.27 20.01 2.63e - 03 0.99 967 788 3325
0.0078125 104.77 20.58 6.76e - 04 0.97 1828 1644 7717
0.00390625 120.28 21.25 1.66e - 04 1.00 3520 3315 15048
0.00195312 116.67 20.94 4.10e - 05 0.99 7169 6587 40360
0.000976562 119.16 21.27 1.03e - 05 1.00 14026 13530 66011
μ = 100 0.5 13.03 12.03 17.13 0.97 59 23 90
0.25 10.08 9.42 6.18 0.91 105 45 131
0.125 11.39 10.53 1.18 0.97 210 92 364
0.0625 13.48 12.51 2.91e - 01 0.97 401 180 971
0.03125 12.51 11.61 6.89e - 02 0.99 732 424 2046
0.015625 13.04 12.07 1.79e - 02 1.00 1354 872 4121
0.0078125 13.14 12.15 4.50e - 03 0.99 2622 1735 7951
0.00390625 13.48 12.45 1.10e - 03 1.00 5080 3770 17808
0.00195312 13.54 12.48 6.82e - 05 1.00 9976 15286 77916
0.000976562 13.60 12.54 7.01e - 05 1.00 19640 14987 71481

In [29] it is proven, up to higher order terms, that edge residuals dominate for the Laplace problem and linear finite elements method on anisotropic meshes. We numerically compare now, error estimator (4.12) with the simplified error indicator, where only edge residuals are considered

(6.6) η ~ 2 , K = ( 1 2 i = 1 d + 1 ( | K i | j = 1 d λ j , K ) 1 2 [ ( μ + | u h | p - 2 ) u h 𝐧 ] L 2 ( K i ) ) ω 2 , K ( u - u h ) .

Let again Ω = ( 0 , 1 ) 2 and consider again problem (2.1) with exact solution (5.4) with ϵ = 0.05 and μ = 0 . As before, we apply the adaptive algorithm to a starting mesh of size h = 0.1 . The starting tolerance is TOL = 0.5 , 40 iterations of the adaptive algorithm are performed. The tolerance is then halved and the process is repeated. We compare results obtained in Table 5 when adapting with the local error estimator (4.12) and results obtained when considering (6.6). We define as in (5.1)–(5.2) the effectivity index of the modified error indicator

ei ~ Q N = K 𝒯 h η ~ 2 , K e Q N .

In Table 6 results obtained for both error indicators η 2 , K and η ~ 2 , K are reported. Both error indicators give similar results, in particular the effectivity index ei Q N is around 23 while ei ~ Q N is around 9. In Figure 4 two adapted meshes with their respective solution can be observed. Error indicators η 2 , K and η ~ 2 , K seem to be equivalent with different magnitudes.

Table 6

True errors, effectivity indices for different values of tolerance TOL, when solving problem (2.1) with f such that the exact solution is given by (5.4) with ϵ = 0.05 and μ = 0 . Left: Adapting with respect to the local error estimator (4.12). Right: Adapting with respect to the local error indicator (6.6).

TOL ei Q N e Q N N v av ar ei ~ Q N e Q N N v av ar
0.03125 20.95 9.93e - 03 492 350 8.24 2.47e - 02 306 241
0.015625 22.45 2.34e - 03 863 844 8.80 5.90e - 03 587 491
0.0078125 23.76 5.73e - 04 1703 1661 9.03 1.49e - 03 1133 1008
0.00390625 22.91 1.46e - 04 3425 3286 9.37 3.68e - 04 2220 2040
Figure 4 
               Left: Adapted mesh obtained when 
                     
                        
                           
                              TOL
                              =
                              0.25
                           
                        
                        
                        {\mathrm{TOL}=0.25}
                     
                   using the local error estimator (4.12), when u is given by (5.4) with 
                     
                        
                           
                              ϵ
                              =
                              0.05
                           
                        
                        
                        {\epsilon=0.05}
                     
                   and 
                     
                        
                           
                              μ
                              =
                              0
                           
                        
                        
                        {\mu=0}
                     
                  . Right: Adapted mesh obtained when 
                     
                        
                           
                              TOL
                              =
                              0.125
                           
                        
                        
                        {\mathrm{TOL}=0.125}
                     
                   using the local error indicator (6.6), when u is given by (5.4) with 
                     
                        
                           
                              ϵ
                              =
                              0.05
                           
                        
                        
                        {\epsilon=0.05}
                     
                   and 
                     
                        
                           
                              μ
                              =
                              0
                           
                        
                        
                        {\mu=0}
                     
                  .
Figure 4 
               Left: Adapted mesh obtained when 
                     
                        
                           
                              TOL
                              =
                              0.25
                           
                        
                        
                        {\mathrm{TOL}=0.25}
                     
                   using the local error estimator (4.12), when u is given by (5.4) with 
                     
                        
                           
                              ϵ
                              =
                              0.05
                           
                        
                        
                        {\epsilon=0.05}
                     
                   and 
                     
                        
                           
                              μ
                              =
                              0
                           
                        
                        
                        {\mu=0}
                     
                  . Right: Adapted mesh obtained when 
                     
                        
                           
                              TOL
                              =
                              0.125
                           
                        
                        
                        {\mathrm{TOL}=0.125}
                     
                   using the local error indicator (6.6), when u is given by (5.4) with 
                     
                        
                           
                              ϵ
                              =
                              0.05
                           
                        
                        
                        {\epsilon=0.05}
                     
                   and 
                     
                        
                           
                              μ
                              =
                              0
                           
                        
                        
                        {\mu=0}
                     
                  .
Figure 4

Left: Adapted mesh obtained when TOL = 0.25 using the local error estimator (4.12), when u is given by (5.4) with ϵ = 0.05 and μ = 0 . Right: Adapted mesh obtained when TOL = 0.125 using the local error indicator (6.6), when u is given by (5.4) with ϵ = 0.05 and μ = 0 .

7 An Industrial Application: Aluminium Electrolysis

Our objective is to apply our adaptive strategy to a fluid flow simulation arising from aluminium electrolysis [46]. Let Ω 3 be the domain occupied by two fluids, namely the liquid electrolyte and the liquid aluminium, which are immiscible and have different densities. Let 𝐮 : Ω 3 and p : Ω be the velocity and the pressure of the fluids in the cell. Let 𝐅 denote the gravity and electromagnetic forces. For a given interface between the two fluids, we are looking for 𝐮 , p such that

(7.1) { ρ ( 𝐮 ) 𝐮 - ( 2 ( μ ( | ϵ ( 𝐮 ) | ) ) ϵ ( 𝐮 ) ) + p = 𝐅 in  Ω , 𝐮 = 0 in  Ω .

On Ω no slip boundary conditions are enforced. Here

ϵ i j ( 𝐮 ) = 1 2 ( u i x j + u j x i )

and ρ is the density, which is piecewise constant in the liquid electrolyte and the liquid aluminium. Smagorinsky’s model [46, 45] is considered

μ ( | ϵ ( 𝐮 ) | ) = μ L + ξ turb ρ l 2 | ϵ ( 𝐮 ) | ,

where 𝝁 L , ξ turb and l are positive numbers and

| ϵ ( 𝐮 ) | = ( i , j 3 ( ϵ i j ( 𝐮 ) ) 2 ) 1 2 .

At the interface between the two fluids a non-penetration condition is imposed: 𝐮 𝐧 = 0 so as zero tangential forces [ ( 2 μ ϵ ( 𝐮 ) - p I d ) 𝐧 𝐭 i ] = 0 for i = 1 , 2 . Our goal is to reduce the CPU time using an anisotropic adapted mesh. The numerical results obtained using this adapted mesh will be compared to those obtained with a reference industrial mesh. The classical mini element is used, leading to an approximation 𝐮 h and p h of (7.1) as in [46]. Starting from a given mesh, we apply the adaptive strategy of Section 6. We denote 𝒯 h a mesh of Ω in tetrahedron K. From discussion of the previous section, the following simplified error indicator with d = 3 is considered:

K 𝒯 h 1 λ 3 , K 1 2 [ ( μ L + ξ turb ρ l 2 | ϵ ( 𝐮 h ) | ) ϵ ( 𝐮 h ) 𝐧 ] L 2 ( K ) ω 2 , K ( 𝐮 - 𝐮 h ) .

Here u h is the linear part of the finite element velocity, for 𝐯 ( H 0 1 ( Ω ) ) 3 we have defined

( ω 2 , K ( 𝐯 ) ) 2 = i = 1 3 ( ω 2 , K ( v i ) ) 2 ,

where v i are the components of 𝐯 and ω 2 , K ( v i ) is defined in (3.3). We choose a starting mesh of 10693 vertices, we run ten iterations of the adaptive algorithm with TOL = 1 , ten additional iterations with TOL = 0.5 and 20 more iterations with TOL = 0.375 . The final mesh has 24255 vertices and the CPU time of the whole iterative process is 4.08 hours.

Once the anisotropic adapted mesh is obtained, it is used to solve a free interface problem [46] to equilibrate the normal forces at the electrolyte-aluminium surface [ ( 2 μ ϵ ( 𝐮 ) - p I d ) 𝐧 𝐧 ] = 0 . This requires to compute ten more Navier–Stokes approximated solutions, the CPU time needed being 0.96 hours. On the other hand, the reference (non-adapted) industrial mesh has 326099 vertices and the computational time to solve the free-surface problem is 5.78 hours. Therefore, once the adapted mesh is available, solving the interface Navier–Stokes problem requires five times less CPU than when using the standard mesh. It should be mentioned that the same adapted mesh can be used for several computations. In Figures 5, 6 and 7, results obtained with the reference industrial mesh and the adapted one are shown. In Figures 8 and 9, we report a plot over two different lines of the velocity magnitude. Similar accuracy between the adapted mesh and the reference one can be observed.

Figure 5

View of the computational domain with corresponding mesh.

(a) 
                  Reference industrial mesh (326099 vertices).
(a)

Reference industrial mesh (326099 vertices).

(b) 
                  Adapted mesh (24255 vertices).
(b)

Adapted mesh (24255 vertices).

Figure 6

Cut at x = - 6 of the fluid domain.

(a) 
                  Reference industrial mesh (326099 vertices) and solution after solving the free-surface problem.
(a)

Reference industrial mesh (326099 vertices) and solution after solving the free-surface problem.

(b) 
                  Adapted mesh (24255 vertices) and solution after solving the free-surface problem.
(b)

Adapted mesh (24255 vertices) and solution after solving the free-surface problem.

Figure 7

Cut at z = 0.3 of the fluid domain, view from below.

(a) 
                  Reference industrial mesh (326099 vertices) and solution after solving the free-surface problem.
(a)

Reference industrial mesh (326099 vertices) and solution after solving the free-surface problem.

(b) 
                  Adapted mesh (24255 vertices) and solution after solving the free-surface problem.
(b)

Adapted mesh (24255 vertices) and solution after solving the free-surface problem.

Figure 8 
               Plot along y at 
                     
                        
                           
                              x
                              =
                              
                                 -
                                 6.0
                              
                           
                        
                        
                        {x=-6.0}
                     
                   and 
                     
                        
                           
                              z
                              =
                              0.3
                           
                        
                        
                        {z=0.3}
                     
                  . The velocity magnitude, computed solving the free interface problem, for the adapted mesh and the reference one is reported.
Figure 8

Plot along y at x = - 6.0 and z = 0.3 . The velocity magnitude, computed solving the free interface problem, for the adapted mesh and the reference one is reported.

Figure 9 
               Plot along z at 
                     
                        
                           
                              x
                              =
                              
                                 -
                                 6.0
                              
                           
                        
                        
                        {x=-6.0}
                     
                   and 
                     
                        
                           
                              y
                              =
                              
                                 -
                                 1.0
                              
                           
                        
                        
                        {y=-1.0}
                     
                  . The velocity magnitude, computed solving the free interface problem, for the adapted mesh and the reference one is reported.
Figure 9

Plot along z at x = - 6.0 and y = - 1.0 . The velocity magnitude, computed solving the free interface problem, for the adapted mesh and the reference one is reported.

8 Conclusion and Perspectives

An anisotropic error estimator for the approximation of a p-Laplacian problem is introduced. The constants involved in the upper and lower bounds are shown to be independent from the aspect ratio, the solution and μ, up to higher order terms. An adaptive algorithm based on the derived error estimator (4.12) is presented. Numerical experiments for p = 3 show that the proposed error estimator is accurate for both the W 0 1 , 3 ( Ω ) norm and the quasi-norm, when μ is large and only for the quasi-norm when μ is small, as already shown for isotropic finite elements. Effectivity indices for the quasi-norm remain in [ 11 , 24 ] interval, when μ ranges from 0 to 100 and when adapted anisotropic meshes are used when TOL is small enough. Finally, a simplified error estimator, which is used to obtain an adapted mesh in the framework of aluminium electrolysis is presented. The results show a reduction of the computational time, while keeping a comparable accuracy with respect to a reference industrial mesh. An extension to higher order finite elements should be considered first when p = 2 , then when p 2 .

Funding statement: Paride Passelli is financed by Rio Tinto Aluminium LRF Research Center at Saint Jean de Maurienne (EPFL industrial grant).

Acknowledgements

Alexei Lozinski is acknowledged for advice. Maude Girardin, Samuel Dubuis and Jacques Rappaz are acknowledged for discussions. The referees are acknowledged for fruitful suggestions.

References

[1] M. Ainsworth and J. T. Oden, A posteriori error estimation in finite element analysis, Comput. Methods Appl. Mech. Engrg. 142 (1997), no. 1–2, 1–88. 10.1016/S0045-7825(96)01107-3Search in Google Scholar

[2] M. Ainsworth, J. Z. Zhu, A. W. Craig and O. C. Zienkiewicz, Analysis of the Zienkiewicz–Zhu a posteriori error estimator in the finite element method, Internat. J. Numer. Methods Engrg. 28 (1989), no. 9, 2161–2174. 10.1002/nme.1620280912Search in Google Scholar

[3] F. Alauzet and A. Loseille, A decade of progress on anisotropic mesh adaptation for computational fluid dynamics, Comput.-Aided Des. 72 (2016), 13–39. 10.1016/j.cad.2015.09.005Search in Google Scholar

[4] I. Babuška, J. Chandra and J. E. Flaherty, Adaptive Computational Methods for Partial Differential Equations, Society for Industrial and Applied Mathematics, Philadelphia, 1983. Search in Google Scholar

[5] I. Babuška, R. Durán and R. Rodríguez, Analysis of the efficiency of an a posteriori error estimator for linear triangular finite elements, SIAM J. Numer. Anal. 29 (1992), no. 4, 947–964. 10.1137/0729058Search in Google Scholar

[6] E. Bänsch, Local mesh refinement in 2 and 3 dimensions, Impact Comput. Sci. Engrg. 3 (1991), no. 3, 181–191. 10.1016/0899-8248(91)90006-GSearch in Google Scholar

[7] J. W. Barrett and W. B. Liu, Quasi-norm error bounds for the finite element approximation of a non-Newtonian flow, Numer. Math. 68 (1994), no. 4, 437–456. 10.1007/s002110050071Search in Google Scholar

[8] L. Belenki, L. Diening and C. Kreuzer, Optimality of an adaptive finite element method for the p-Laplacian equation, IMA J. Numer. Anal. 32 (2012), no. 2, 484–510. 10.1093/imanum/drr016Search in Google Scholar

[9] E. Boey, Y. Bourgault and T. Giordano, Anisotropic space-time adaptation for reaction-diffusion problems, preprint (2017), https://arxiv.org/abs/1707.04787. Search in Google Scholar

[10] Y. Bourgault and M. Picasso, Anisotropic error estimates and space adaptivity for a semidiscrete finite element approximation of the transient transport equation, SIAM J. Sci. Comput. 35 (2013), no. 2, A1192–A1211. 10.1137/120891320Search in Google Scholar

[11] W. Cao, Superconvergence analysis of the linear finite element method and a gradient recovery postprocessing on anisotropic meshes, Math. Comp. 84 (2015), no. 291, 89–117. 10.1090/S0025-5718-2014-02846-9Search in Google Scholar

[12] C. Carstensen, All first-order averaging techniques for a posteriori finite element error control on unstructured grids are efficient and reliable, Math. Comp. 73 (2004), no. 247, 1153–1165. 10.1090/S0025-5718-03-01580-1Search in Google Scholar

[13] C. Carstensen and R. Klose, A posteriori finite element error control for the p-Laplace problem, SIAM J. Sci. Comput. 25 (2003), no. 3, 792–814. 10.1137/S1064827502416617Search in Google Scholar

[14] C. Carstensen, W. Liu and N. Yan, A posteriori FE error control for p-Laplacian by gradient recovery in quasi-norm, Math. Comp. 75 (2006), no. 256, 1599–1616. 10.1090/S0025-5718-06-01819-9Search in Google Scholar

[15] P. Clément, Approximation by finite element functions using local regularization, Rev. Franç. Automat. Inform. Rech. Opér. Sér. Rouge Anal. Numér. 9 (1975), no. R-2, 77–84. 10.1051/m2an/197509R200771Search in Google Scholar

[16] T. Coupez, Metric construction by length distribution tensor and edge based error for anisotropic adaptive meshing, J. Comput. Phys. 230 (2011), no. 7, 2391–2405. 10.1016/j.jcp.2010.11.041Search in Google Scholar

[17] L. Diening and C. Kreuzer, Linear convergence of an adaptive finite element method for the p-Laplacian equation, SIAM J. Numer. Anal. 46 (2008), no. 2, 614–638. 10.1137/070681508Search in Google Scholar

[18] W. Dörfler, A convergent adaptive algorithm for Poisson’s equation, SIAM J. Numer. Anal. 33 (1996), no. 3, 1106–1124. 10.1137/0733054Search in Google Scholar

[19] S. Dubuis, Adaptive algorithms for two fluids flows with anisotropic finite elements and order two time discretizations, Ph.D. Thesis, Ecole polytechnique fédérale de Lausanne, 2019. Search in Google Scholar

[20] S. Dubuis, P. Passelli and M. Picasso, Anisotropic adaptive finite elements for an elliptic problem with strongly varying diffusion coefficient, Comput. Methods Appl. Math. 22 (2022), no. 3, 529–543. 10.1515/cmam-2022-0036Search in Google Scholar

[21] L. Formaggia, S. Micheletti and S. Perotto, Anisotropic mesh adaption in computational fluid dynamics: Application to the advection-diffusion-reaction and the Stokes problems, Appl. Numer. Math. 51 (2004), no. 4, 511–533. 10.1016/j.apnum.2004.06.007Search in Google Scholar

[22] L. Formaggia and S. Perotto, New anisotropic a priori error estimates, Numer. Math. 89 (2001), no. 4, 641–667. 10.1007/s002110100273Search in Google Scholar

[23] L. Formaggia and S. Perotto, Anisotropic error estimates for elliptic problems, Numer. Math. 94 (2003), no. 1, 67–92. 10.1007/s00211-002-0415-zSearch in Google Scholar

[24] P. J. Frey and F. Alauzet, Anisotropic mesh adaptation for CFD computations, Comput. Methods Appl. Mech. Engrg. 194 (2005), no. 48–49, 5068–5082. 10.1016/j.cma.2004.11.025Search in Google Scholar

[25] R. Glowinski and A. Marrocco, Sur l’approximation, par éléments finis d’ordre un, et la résolution, par pénalisation-dualité, d’une classe de problèmes de Dirichlet non linéaires, Rev. Franç. Automat. Inform. Rech. Opér. Sér. Rouge Anal. Numér. 9 (1975), no. R-2, 41–76. 10.1051/m2an/197509R200411Search in Google Scholar

[26] O. Gorynina, A. Lozinski and M. Picasso, Time and space adaptivity of the wave equation discretized in time by a second-order scheme, IMA J. Numer. Anal. 39 (2019), no. 4, 1672–1705. 10.1093/imanum/dry048Search in Google Scholar

[27] G. Kunert, An a posteriori residual error estimator for the finite element method on anisotropic tetrahedral meshes, Numer. Math. 86 (2000), no. 3, 471–490. 10.1007/s002110000170Search in Google Scholar

[28] G. Kunert and S. Nicaise, Zienkiewicz–Zhu error estimators on anisotropic tetrahedral and triangular finite element meshes, M2AN Math. Model. Numer. Anal. 37 (2003), no. 6, 1013–1043. 10.1051/m2an:2003065Search in Google Scholar

[29] G. Kunert and R. Verfürth, Edge residuals dominate a posteriori error estimates for linear finite element methods on anisotropic triangular and tetrahedral meshes, Numer. Math. 86 (2000), no. 2, 283–303. 10.1007/PL00005407Search in Google Scholar

[30] P. Laug and H. Borouchaki, The BL2D mesh generator – Beginner’s guide, user’s and programmer’s manual, Technical Report RT-0194, Institut National de Recherche en Informatique et Automatique (INRIA), Rocquencourt-Le Chesnay, 1996. Search in Google Scholar

[31] J.-L. Lions, Optimal Control of Systems Governed by Partial Differential Equations, Grundlehren Math. Wiss. 170, Springer, New York, 1971. 10.1007/978-3-642-65024-6Search in Google Scholar

[32] W. Liu and N. Yan, Quasi-norm local error estimators for p-Laplacian, SIAM J. Numer. Anal. 39 (2001), no. 1, 100–127. 10.1137/S0036142999351613Search in Google Scholar

[33] W. Liu and N. Yan, Some a posteriori error estimators for p-Laplacian based on residual estimation or gradient recovery, J. Sci. Comput. 16 (2001), no. 4, 435–477. Search in Google Scholar

[34] W. B. Liu and J. W. Barrett, Quasi-norm error bounds for the finite element approximation of some degenerate quasilinear elliptic equations and variational inequalities, RAIRO Modél. Math. Anal. Numér. 28 (1994), no. 6, 725–744. 10.1051/m2an/1994280607251Search in Google Scholar

[35] A. Loseille and F. Alauzet, Continuous mesh framework part I: Well-posed continuous interpolation error, SIAM J. Numer. Anal. 49 (2011), no. 1, 38–60. 10.1137/090754078Search in Google Scholar

[36] A. Lozinski, M. Picasso and V. Prachittham, An anisotropic error estimator for the Crank–Nicolson method: Application to a parabolic problem, SIAM J. Sci. Comput. 31 (2009), no. 4, 2757–2783. 10.1137/080715135Search in Google Scholar

[37] S. Micheletti and S. Perotto, Reliability and efficiency of an anisotropic Zienkiewicz–Zhu error estimator, Comput. Methods Appl. Mech. Engrg. 195 (2006), no. 9–12, 799–835. 10.1016/j.cma.2005.02.009Search in Google Scholar

[38] S. Micheletti, S. Perotto and M. Picasso, Stabilized finite elements on anisotropic meshes: A priori error estimates for the advection-diffusion and the Stokes problems, SIAM J. Numer. Anal. 41 (2003), no. 3, 1131–1162. 10.1137/S0036142902403759Search in Google Scholar

[39] J.-M. Mirebeau, Optimally adapted meshes for finite elements of arbitrary order and W 1 , p norms, Numer. Math. 120 (2012), no. 2, 271–305. 10.1007/s00211-011-0412-1Search in Google Scholar

[40] M. Picasso, Adaptive finite elements for a linear parabolic problem, Comput. Methods Appl. Mech. Engrg. 167 (1998), no. 3–4, 223–237. 10.1016/S0045-7825(98)00121-2Search in Google Scholar

[41] M. Picasso, An anisotropic error indicator based on Zienkiewicz–Zhu error estimator: Application to elliptic and parabolic problems, SIAM J. Sci. Comput. 24 (2003), no. 4, 1328–1355. 10.1137/S1064827501398578Search in Google Scholar

[42] M. Picasso, Numerical study of the effectivity index for an anisotropic error indicator based on Zienkiewicz–Zhu error estimator, Comm. Numer. Methods Engrg. 19 (2003), no. 1, 13–23. 10.1002/cnm.546Search in Google Scholar

[43] M. Picasso, Adaptive finite elements with large aspect ratio based on an anisotropic error estimator involving first order derivatives, Comput. Methods Appl. Mech. Engrg. 196 (2006), no. 1–3, 14–23. 10.1016/j.cma.2005.11.018Search in Google Scholar

[44] M. Picasso, Numerical study of an anisotropic error estimator in the L 2 ( H 1 ) norm for the finite element discretization of the wave equation, SIAM J. Sci. Comput. 32 (2010), no. 4, 2213–2234. 10.1137/090778249Search in Google Scholar

[45] S. B. Pope, Turbulent flows, Meas. Sci. Technol 12 (2001), Paper No. 11. 10.1088/0957-0233/12/11/705Search in Google Scholar

[46] J. Rochat, Approximation numérique des coulements turbulents dans des cuves d’électrolyse de l’aluminium, Ph.D. Thesis, Ecole polytechnique fédérale de Lausanne, 2016. Search in Google Scholar

[47] R. Verfürth, A posteriori error estimation and adaptive mesh-refinement techniques, J. Comput. Appl. Math. 50 (1994), no. 1–3, 67–83. 10.1016/0377-0427(94)90290-9Search in Google Scholar

[48] R. Verfürth, A Review of A Posteriori Error Estimation and Adaptive Mesh-Refinement Techniques, Wiley-Teubner, New York, 1996. Search in Google Scholar

[49] R. Verfürth, A Posteriori Error Estimation Techniques for Finite Element Methods, Oxford University, Oxford, 2013. 10.1093/acprof:oso/9780199679423.001.0001Search in Google Scholar

[50] J. Xu and Z. Zhang, Analysis of recovery type a posteriori error estimators for mildly structured grids, Math. Comp. 73 (2004), no. 247, 1139–1152. 10.1090/S0025-5718-03-01600-4Search in Google Scholar

[51] O. C. Zienkiewicz and J. Z. Zhu, A simple error estimator and adaptive procedure for practical engineering analysis, Internat. J. Numer. Methods Engrg. 24 (1987), no. 2, 337–357. 10.1002/nme.1620240206Search in Google Scholar

[52] O. C. Zienkiewicz and J. Z. Zhu, The superconvergent patch recovery and a posteriori error estimates. I. The recovery technique, Internat. J. Numer. Methods Engrg. 33 (1992), no. 7, 1331–1364. 10.1002/nme.1620330702Search in Google Scholar

[53] Spatial Corp Headquarters, Broomfield, 3D Precise Mesh, https://www.3ds.com/. Search in Google Scholar

Received: 2022-10-12
Revised: 2024-06-04
Accepted: 2024-06-06
Published Online: 2024-06-26
Published in Print: 2025-04-01

© 2024 Walter de Gruyter GmbH, Berlin/Boston

This work is licensed under the Creative Commons Attribution 4.0 International License.

Downloaded on 16.9.2025 from https://www.degruyterbrill.com/document/doi/10.1515/cmam-2022-0205/html
Scroll to top button