The gradient discretisation method for nonlinear porous media equations
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 uniformintime convergence for the approximate solution, without assuming nonphysical 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, (masslumped) finite elements, etc. The theoretical results are illustrated, in both fast and slow diffusion regimes, by numerical tests based on masslumped conforming finite elements.
2010 Mathematics Subject Classification:
65M08, 65M12, 65M60, 76S05labelkeyrgb0.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
(1)  
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
(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 spacetime 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
where and . Error estimates in nonstandard 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 nonphysical 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 implicitintime 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 uniformintime 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 masslumped 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
2.1. Gradient scheme
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 spacetime gradient discretisation for homogeneous Dirichlet boundary conditions, with piecewise constant reconstruction, if

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

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 ,

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

is an interpolation operator,

.
We then let and .
For any , we define the piecewiseconstantintime functions , and by: For , for any , for a.e.
Note that is defined everywhere, including at . This will be required to state a uniformintime convergence result on this reconstructed function.
If and satisfies , we define . The piecewise constant feature of then shows that
(3) 
Once a GD has been chosen, an implicitintime 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 spacetime gradient discretisations are required to satisfy consistency, limitconformity 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 spacetime gradient discretisations in the sense of Definition 2.1 is said to be consistent if

for all , letting
we have as ,

for all , in as

as .
Definition 2.4 (Limitconformity).
A sequence of spacetime gradient discretisations in the sense of Definition 2.1 is said to be limitconformity if, for all , letting
we have as .
Definition 2.5 (Compactness).
A sequence of spacetime gradient discretisations in the sense of Definition 2.1 is said to be compact if
where
and has been extended by outside .
A sequence of GDs that is compact or limitconforming 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 spacetime gradient discretisations in the sense of Definition 2.1 is compact or limitconforming, then it is coercive: there exists a constant such that
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).
For a given , we set . We also set
(5) 
Definition 2.7 (Weak solution to (1)).
Assume that , and . A weak solution to (1) is a function such that

and in ,

, ,

, and for any
(6)
The existence of a weak solution to (1) will be obtained as a byproduct of our convergence analysis.
Theorem 2.8 (Convergence of the gradient scheme).
Let be a sequence of gradient discretisations that is consistent, limitconforming 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.
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
(11) 
Plugging that in (9) and using a telescopic sum yields (7). The a priori estimates (8) follow from this relation, by writing
Noting that for all and recalling that
(12) 
we obtain the estimates in (8). ∎
The following corollary guarantees the existence of a solution to (4).
Corollary 3.2.
Proof.
In the following lemma, we estimate the dual norm of the discrete time derivative . The dual norm on is defined by
Lemma 3.3.
There exists a constant depending on , and such that
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
where . Taking the supremum over all satisfying , we deduce
Divide by , square, and sum over (and use Jensen or Cauchy–Schwarz to estimate the term involving ) to get
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 Lipschitzcontinuous at 0, whereas for slow diffusion, the Lipschitz constant of explodes at . We therefore define and for any as
(13) 
and
The cutoff functions and are globally Lipschitz continuous with respective Lipschitz constants and . Our goal here is to estimate the timetranslates 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
Proof.
By noting that for and
we obtain
The above inequalities imply that for any ,
or, equivalently,
(14) 
Together with the global Lipschitz property of and the monotonicity of , this yields
which completes the proof. ∎
We can now estimate the time translates of .
Lemma 3.5.
Proof.
In this proof, denotes a generic constant that has the same dependencies as in the lemma. For and any we have
(17) 
with
and
To estimate , define
The commutativity property (3) shows that and thus, when ,
(18) 
By the Hölder and Chebyschev inequalities, (3) and (8), we estimate :
(19) 
Similarly, we also have estimates for :
(20) 
We estimate using the same arguments in [9, Lemma 4.4]. Thanks to Lemma 3.4 we have
(21) 
with
For any there exists such that . We rewrite by expressing as the sum of its jumps in time, and use the definition of dual seminorm to infer
Together with the Cauchy–Schwarz inequality (on the sums and the integrals) and (8), this implies
(22) 
For any