The gradient discretisation method for nonlinear porous media equations

# The gradient discretisation method for nonlinear porous media equations

Jérôme Droniou School of Mathematical Sciences, Monash University, Clayton, Victoria 3800, Australia. jerome.droniou@monash.edu  and  Kim-Ngan Le School of Mathematical Sciences, Monash University, Clayton, Victoria 3800, Australia. ngan.le@monash.edu
July 2, 2019
###### Abstract.

The gradient discretisation method (GDM) is a generic framework for designing and analysing numerical schemes for diffusion models. In this paper, we study the GDM for the porous medium equation, including fast diffusion and slow diffusion models. Using discrete functional analysis techniques, we establish a strong -convergence of the approximate gradients and a uniform-in-time convergence for the approximate solution, without assuming non-physical regularity assumptions on the data or continuous solution. Being established in the generic GDM framework, these results apply to a variety of numerical methods, such as finite volume, (mass-lumped) finite elements, etc. The theoretical results are illustrated, in both fast and slow diffusion regimes, by numerical tests based on mass-lumped conforming finite elements.

###### 2010 Mathematics Subject Classification:
65M08, 65M12, 65M60, 76S05
\definecolor

labelkeyrgb0.6,0,1 \definecolorlabelkeyrgb0.6,0,1 \definecolorvioletrgb0.580,0.,0.827 \definecolorshadecolorgray0.92 \definecolorTFFrameColorgray0.92 \definecolorTFTitleColorrgb0,0,0 \DeclareDocumentCommand\RPiD OD O,0 Π_#1(X_#1#2)

## 1. Introduction

In this paper, we study nonlinear porous media equations of the form

 ∂tu−Δβ(u)=fin ΩT:=(0,T)×Ω, u(0,⋅)=u0in Ω, (1) β(u)=0on (0,T)×∂Ω,

where with , is an open bounded domain of () with boundary , and .

In the case , the equation (1) is the standard model of diffusion of a gas in porous media, which is also called the slow diffusion model. The case describes the flow of an ideal gas in porous media while that of a diffusion of a compressible fluid through porous media. Other choices of exponent appear in different physical situations, such as thermal propagation in plasma () or plasma radiation (). The slow diffusion equation is not uniformly elliptic as it degenerates at unknown points where . An interesting feature is that, if is compactly supported, then at any the support of has a free boundary with a finite speed of propagation, see e.g. [5, 18].

In the fast diffusion case , the equation (1) is relevant in the description of plasma physics, the kinetic theory of gas or fluid transportation in porous media [3, 6, 19]. Since the modulus of ellipticity blows up whenever vanishes, the fast diffusion equation is a singular equation. An essential difference between the fast and slow diffusion cases is that for fast diffusion the solution decays to zero in some finite time depending on the initial data while for slow diffusion the solution decays to zero in infinite time like an inverse power of , see e.g. [4, 7].

If we formally represent (1) as

 ∂tu−div(β′(u)∇u)=f, (2)

where , then the above equation is a classical nonlinear diffusion equation. The term in (2) induces the nonlinearity that raises many challenges in the analysis of the porous media equations. These problems have been studied extensively both in theory, see e.g. [4, 17, 19], numerical analysis, see e.g. [11, 12, 13, 1], as well as in numerical approximations (without convergence proof) [14].

In particular, authors in [11] study equation (2) with variable exponent of nonlinearity, i.e. is replaced by where (i.e. in our case). In order to deal with the degeneracy in the problem, an approximate regularized problem is investigated. A space-time discretization scheme using the finite element method in space and the discontinuous Galerkin method in time is proposed for the regularized model (not the original problem). Furthermore, error estimates are obtained with strong regularity assumptions on the solution of the regularized model, which are not expected to be satisfied by classical solutions to (1) such as Barrenblatt solutions [15]. A space and time dependent exponent is studied in [1] using the same method as in [11]. The slow diffusion case () is studied in [12] where a fully discrete scheme is proposed and error estimates are proved with strong assumptions on the solution and the pressure . In [13], the author studies the fast and slow diffusion cases in which a fully discrete Galerkin approximation is considered for the transformed Richards model

 ∂tψ(w)−Δw=0,

where and . Error estimates in non-standard quasi norms and rates of convergence are proved with strong regularity assumptions on the solution of (1). In a recent preprint [16], a general nonlinear diffusion equation on the whole space is considered. The authors propose monotone schemes of finite difference type on Cartesian meshes and obtain -convergence for these schemes.

In this paper, we use the gradient discretisation method to approximate (1) and discrete functional analysis techniques to obtain an -convergence result without assuming non-physical regularity assumptions on the data. The gradient discretisation method (GDM) is a generic framework for the design and analysis of numerical schemes for diffusion models. Using only a few discrete elements (a space, and a function and gradient reconstruction), it describes a variety of numerical schemes – such as finite volumes, finite elements, discontinuous Galerkin, etc. – and identifies three key properties of the discrete elements that ensure the convergence for linear and nonlinear models. Analyses carried out within the GDM therefore yield convergence theorems that directly apply to all the methods that fit into this framework. We refer to the monograph [10] for a detailed presentation of the GDM, and of the methods it contains.

The paper is organised as follows. Section 2.1 recalls notations of the GDM and introduce an implicit-in-time gradient scheme for (1). The definition of weak solutions and the main convergence result are stated in Section 2.2. We provide a priori estimates for the approximate solution in Section 3. Section 4 contains the initial weak convergence of the gradient scheme. The main convergence result, including the uniform-in-time convergence, is proved in Section 5. Numerical results are provided in Section 6 to illustrate the generic convergence results; to this purpose, we choose a particular gradient scheme, the mass-lumped conforming method, and we evaluate its accuracy when approximating the Barrenblatt solution. These tests are presented for a variety of exponents , including both fast diffusion and slow diffusion ranges. A brief conclusion is presented in Section 7, and technical results are gathered in an appendix (Section 8).

## 2. The gradient discretisation method and the main convergence result

We recall here the notions of the gradient discretisation method. The idea of this general analysis framework is to replace, in the weak formulation of the problem, the continuous space and operators by discrete ones; the set of discrete space and operators is called a gradient discretisation (GD), and the scheme obtained after substituting these elements into the weak formulation is called a gradient scheme (GS). The convergence of the obtained GS can be established based on only a few general concepts on the underlying GD. Moreover, different GDs correspond to different classical schemes (finite elements, finite volumes, etc.). Hence, the analysis carried out in the GDM directly applies to all these schemes, and does not rely on the specificity of each particular method.

###### Definition 2.1.

is a space-time gradient discretisation for homogeneous Dirichlet boundary conditions, with piecewise constant reconstruction, if

1. the set of discrete unknowns is a finite dimensional real vector space,

2. the linear map is a piecewise constant reconstruction operator in the following sense: there exists a basis of and a family of disjoint subsets of such that, for all , it holds , where is the characteristic function of ,

3. the linear mapping gives a reconstructed discrete gradient. It must be chosen such that is a norm on ,

4. is an interpolation operator,

5. .

We then let and .

For any , we define the piecewise-constant-in-time functions , and by: For , for any , for a.e.

 ΠDv(0,x):=ΠDv(0)(x), ΠDv(t,x):=ΠDv(n+1)(x), ∇Dv(t,x):=∇Dv(n+1)(x), δDv(t)=δ(n+12)Dv:=ΠDv(n+1)−ΠDv(n)δt(n+12)∈L2(Ω).

Note that is defined everywhere, including at . This will be required to state a uniform-in-time convergence result on this reconstructed function.

If and satisfies , we define . The piecewise constant feature of then shows that

 ∀v∈XD,0,ΠDg(v)=g(ΠDv). (3)

Once a GD has been chosen, an implicit-in-time gradient scheme for (1) is defined the following way.

###### Algorithm 2.2 (GS for (1)).

Set and let satisfy:

 (4)

for all ‘test function’ .

In order to establish the stability and convergence of the GS (4), sequences of space-time gradient discretisations are required to satisfy consistency, limit-conformity and compactness properties [10]. The consistency is slightly adapted here to account for the nonlinearity we consider. In the following, we let .

###### Definition 2.3 (Consistency).

A sequence of space-time gradient discretisations in the sense of Definition 2.1 is said to be consistent if

• for all , letting

 ^SDl(ϕ):=minw∈XDl(∥ΠDlw−ϕ∥L1+ˆm(Ω)+∥∇Dlw−∇ϕ∥L2(Ω)),

we have as ,

• for all , in as

• as .

###### Definition 2.4 (Limit-conformity).

A sequence of space-time gradient discretisations in the sense of Definition 2.1 is said to be limit-conformity if, for all , letting

 WDl(ϕ):=maxv∈XDl∖{0}∣∣∣∫Ω(∇Dlv(x)⋅ϕ(x)+ΠDlv(x)divϕ(x))dx∣∣∣∥∇Dlv∥L2(Ω),

we have as .

###### Definition 2.5 (Compactness).

A sequence of space-time gradient discretisations in the sense of Definition 2.1 is said to be compact if

 limξ→0supl≥1TDl(ξ)=0,

where

 TDl(ξ):=maxv∈XDl∖{0}∥ΠDlv(⋅+ξ)−ΠDlv∥L2(Rd)∥∇Dlv∥L2(Ω),∀ξ∈Rd,

and has been extended by outside .

A sequence of GDs that is compact or limit-conforming also satisfies another important property: the coercivity [10, Lemmas 2.6 and 2.10].

###### Lemma 2.6 (Coercivity of sequences of GDs).

If a sequence of space-time gradient discretisations in the sense of Definition 2.1 is compact or limit-conforming, then it is coercive: there exists a constant such that

 CDl:=maxv∈XDl∖{0}∥ΠDlv∥L2(Ω)∥∇Dlv∥L2(Ω)≤Cp,∀l≥1.

### 2.2. The main convergence result

To state the main result of this paper, we introduce a space of continuous functions for the weak topology of , and we define a weak solution to (1).

 C([0,T];Lm+1(Ω)w):= space of functions v:[0,T]→Lm+1(Ω) that are continuous for the weak topology of Lm+1(Ω).

For a given , we set . We also set

 ζ(z):=∫z0β(s)ds=1m+1|z|m+1 for all z∈R. (5)
###### Definition 2.7 (Weak solution to (1)).

Assume that , and . A weak solution to (1) is a function such that

1. and in ,

2. , ,

3. , and for any

 ∫T0H−1⟨∂t¯u(t),ϕ(t)⟩H10dt+⟨∇β(¯u),∇ϕ⟩L2(ΩT)=⟨f,ϕ⟩L2(ΩT). (6)

The existence of a weak solution to (1) will be obtained as a by-product of our convergence analysis.

###### Theorem 2.8 (Convergence of the gradient scheme).

Let be a sequence of gradient discretisations that is consistent, limit-conforming and compact. Then, for each there exists solution to the gradient scheme (4) with . Moreover, there exists a weak solution to (1) in the sense of Definition 2.7 such that, up to a subsequence as ,

• strongly in ,

• strongly in .

## 3. A priori estimates

We first provide a priori estimates for the solution to (4), and then deduce its existence in Corollary 3.2. For legibility, we drop the index in sequences of gradient discretisations, and we simply write instead of .

###### Lemma 3.1.

Let be defined by (5), and be a solution of (4). Then, for any integer number , we have

 ∫Ωζ(ΠDu(k))(x)dx+ ∥∇Dβ(u)∥2L2(Ωtk) (7) ≤∫Ωζ(ΠDu(0))(x)dx+⟨f,ΠDβ(u)⟩L2(Ωtk).

Consequently, there exists a constant depending only on , and such that

 ∥ΠDu∥L∞(0,T;Lm+1(Ω))+∥ΠDζ(u)∥L∞(0,T;L1(Ω))+∥∇Dβ(u)∥L2(ΩT)≤C. (8)
###### Proof.

We choose the test function in (4) to write

 (9)

For any and , we estimate the first term in the left hand side of (9), starting from

 (10)

Since is increasing, is convex and thus above its tangent line, which implies for all . Applying this inequality with and , it follows from (10) that

 1δt(n+12)(ζ(ΠDu(n+1))−ζ(ΠDu(n)))≤δDu(t)ΠDβ(u(n+1)). (11)

Plugging that in (9) and using a telescopic sum yields (7). The a priori estimates (8) follow from this relation, by writing

 ∫Ωζ(ΠDu(k))(x)dx +∥∇Dβ(u)∥2L2(Ωtk)≤∫Ωζ(ΠDu(0))(x)dx+⟨f,ΠDβ(u)⟩L2(Ωtk) ≤∫Ωζ(ΠDu(0))(x)dx+C2p2∥f∥2L2(Ωtk)+12C2p∥ΠDβ(u)∥2L2(Ωtk).

Noting that for all and recalling that

 ∥ΠDβ(u)∥L2(Ωtk)≤CD∥∇Dβ(u)∥L2(Ωtk)d≤Cp∥∇Dβ(u)∥L2(Ωtk)d, (12)

we obtain the estimates in (8). ∎

The following corollary guarantees the existence of a solution to (4).

###### Corollary 3.2.

If is a gradient discretisation in the sense of Definition 2.1 then there exists at least one solution to the gradient scheme (4).

###### Proof.

The result is obtained using (8), the properties and , and the same arguments as in [9, Corrollary 4.2]

In the following lemma, we estimate the dual norm of the discrete time derivative . The dual norm on is defined by

 ∀v∈ΠD(XD,0), |v|∗,D:=sup{∫Ωv(x)ΠDϕ(x)dx:ϕ∈XD,0,∥∇Dϕ∥L2(Ω)=1}.
###### Lemma 3.3.

There exists a constant depending on , and such that

 ∫T0|δDu(t)|2∗,Ddt≤C.
###### Proof.

We first fix and then, for any , take in (4) such that and for all . Using the Cauchy–Schwarz inequality and the definition of we deduce

 δt(k+12)⟨δ(k+12)Du,ΠDξ⟩L2(Ω) ≤δt(k+12)∥∇Dβ(u(k+1))∥L2(Ω)∥∇Dξ∥L2(Ω) +δt(k+12)∥f(k+1)∥L2(Ω)CD∥∇Dξ∥L2(Ω),

where . Taking the supremum over all satisfying , we deduce

 δt(k+12)|δ(k+12)Du|∗,D≤δt(k+12)∥∇Dβ(u(k+1))∥L2(Ω)+δt(k+12)∥f(k+1)∥L2(Ω)CD.

Divide by , square, and sum over (and use Jensen or Cauchy–Schwarz to estimate the term involving ) to get

 ∫T0|δDu(t)|2∗,Ddt≤2∥∇Dβ(u)∥2L2(ΩT)d+2C2D∥f∥2L2(ΩT).

This together with (8) complete the proof of the lemma. ∎

To deal with the lack of global Lipschitz estimates on , we introduce cutoff functions. The definition of these functions is different in the case (fast diffusion) and (slow diffusion), because each of these cases correspond to a certain breakdown in the Lipschitz properties of : for fast diffusion, is not Lipschitz-continuous at 0, whereas for slow diffusion, the Lipschitz constant of explodes at . We therefore define and for any as

 for 01/k, (13)

and

 for m>1,βsk(r):=⎧⎨⎩kmfor r≥kβ(r)for −k

The cutoff functions and are globally Lipschitz continuous with respective Lipschitz constants and . Our goal here is to estimate the time-translates of . To achieve this, we will be using the cutoff functions, as well as the following relation.

###### Lemma 3.4.

Recalling that , for any we have

 (βik(b)−βik(a))2≤ˆmLik(b−a)(β(b)−β(a)),% for i=s,f.
###### Proof.

By noting that for and

 (βfk)′(r):={k1−m% for −1/k1/k,
 (βsk)′(r):={0for r<−k or r>k,β′(r)for −k

we obtain

 0≤(βfk)′(r) ≤1mβ′(r)∀r∈R∖{−1/k,0,1/k} 0≤(βsk)′(r) ≤β′(r)∀r∈R∖{−k,0,k}.

The above inequalities imply that for any ,

 ∣∣∣∫ba(βik)′(s)ds∣∣∣≤ˆm∣∣∣∫baβ′(s)ds∣∣∣,for i=s,f,

or, equivalently,

 ∣∣βik(b)−βik(a)∣∣≤ˆm∣∣β(b)−β(a)∣∣,for i=s,f. (14)

Together with the global Lipschitz property of and the monotonicity of , this yields

 (βik(b)−βik(a))2 ≤Lik|b−a|∣∣βik(b)−βik(a)∣∣ ≤Lik|b−a|ˆm∣∣β(b)−β(a)∣∣=ˆmLik(b−a)(β(b)−β(a)),

which completes the proof. ∎

We can now estimate the time translates of .

###### Lemma 3.5.

Let be a gradient discretisation in the sense of Definition 2.1 and be a solution to scheme (4). Then, there exists a constant depending only on , , and such that:

• for ,

 ∥∥ΠDβ(u)(⋅+τ,⋅)−ΠDβ(u)∥∥2L2(ΩT−τ)≤C(δtD+τ)2m/(m+1), (15)
• for ,

 ∥∥ΠDβ(u)(⋅+τ,⋅)−ΠDβ(u)∥∥2L1(ΩT−τ)≤C(δtD+τ)2/(m+1). (16)
###### Proof.

In this proof, denotes a generic constant that has the same dependencies as in the lemma. For and any we have

 ΠDβ(u) (⋅+τ,⋅)−ΠDβ(u)=Γi1,k+Γi2,k+Γi3,k (17)

with

 Γi1,k:=ΠD(βik(u)−β(u)),Γi2,k:=ΠD(β(u)−βik(u))(⋅+τ,⋅),

and

 Γi3,k:=ΠDβik(u)(⋅+τ,⋅)−ΠDβik(u).

To estimate , define

 Ωfk:= {(t,x)∈ΩT−τ:|ΠDu(t,x)|<1/k}, Ωsk:= {(t,x)∈ΩT−τ:|ΠDu(t,x)|>k}.

The commutativity property (3) shows that and thus, when ,

 ∥∥Γf1,k∥∥2L2(ΩT−τ) =∫Ωfk∣∣ΠD(|u|m−1u−k1−mu)∣∣2(t,x)dxdt ≤2∫Ωfk|ΠDu|2m(t,x)+k2−2m|ΠDu|2(t,x)dxdt ≤4k−2m∣∣Ωfk∣∣≤4∣∣ΩT∣∣k−2m. (18)

By the Hölder and Chebyschev inequalities, (3) and (8), we estimate :

 ∥∥Γs1,k∥∥L1(ΩT−τ) ≤∫ΩT−τ\mathbbm1ΩskΠD(|u|m+km)(t,x)dxdt ≤2∫ΩT−τ\mathbbm1Ωsk|ΠDu|m(t,x)dxdt ≤2k−1∥ΠDu∥m+1Lm+1(ΩT−τ)≤Ck−1. (19)

Similarly, we also have estimates for :

 ∥∥Γf2,k∥∥2L2(ΩT−τ)≤Ck−2m,∥∥Γs2,k∥∥L1(ΩT−τ)≤Ck−1. (20)

We estimate using the same arguments in [9, Lemma 4.4]. Thanks to Lemma 3.4 we have

 ∥∥Γi3,k∥∥2L2(ΩT−τ)≤ˆmLik∫T−τ0A(t,τ)dt, (21)

with

 A(t,τ):=∫Ω(ΠDu(t+τ,x)−ΠDu(t,x))(ΠDβ(u)(t+τ,x)−ΠDβ(u)(t,x))dx.

For any there exists such that . We rewrite by expressing as the sum of its jumps in time, and use the definition of dual semi-norm to infer

 A(t,τ) =nt+τ∑j=nt+1δt(j+12)∫Ωδ(j+12)Du(x)ΠD(β(u(nt+τ+1))−β(u(nt+1)))(x)dx ≤∥∇D(β(u(nt+τ+1))−β(u(nt+1)))∥L2(Ω)nt+τ∑j=nt+1δt(j+12)|δ(j+12)Du|∗,D =∥∇Dβ(u)(t+τ)−∇Dβ(u)(t)∥L2(Ω)nt+τ∑j=nt+1δt(j+12)|δ(j+12)Du|∗,D.

Together with the Cauchy–Schwarz inequality (on the sums and the integrals) and (8), this implies

 ∫T−τ0 ×[∫T−τ0(nt+τ∑j=nt+1δt(j+12)|δ(j+12)Du|∗,D)2dt]1/2 ≤C[∫T−τ0(nt+τ∑j=nt+1δt(j+12))(nt+τ∑j=nt+1δt(j+12)|δ(j+12)Du|2∗,D)dt]1/2. (22)

For any