# Two Finite Element Approaches For The Porous Medium Equation That Are
Positivity Preserving And Energy Stable

Arjun Vijaywargiya  and  Guosheng Fu  Department of Applied and Computational
Mathematics and Statistics, University of Notre Dame, USA<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

In this work, we present the construction of two distinct finite element
approaches to solve the Porous Medium Equation (PME). In the first approach,
we transform the PME to a log-density variable formulation and construct a
continuous Galerkin method. In the second approach, we introduce additional
potential and velocity variables to rewrite the PME into a system of
equations, for which we construct a mixed finite element method. Both
approaches are first-order accurate, mass conserving, and proved to be
unconditionally energy stable for their respective energies. The mixed
approach is shown to preserve positivity under a CFL condition, while a much
stronger property of unconditional bound preservation is proved for the log-
density approach. A novel feature of our schemes is that they can handle
compactly supported initial data without the need for any perturbation
techniques. Furthermore, the log-density method can handle unstructured grids
in any number of dimensions, while the mixed method can handle unstructured
grids in two dimensions. We present results from several numerical experiments
to demonstrate these properties.

###### Key words and phrases:

Porous medium equation; Mixed finite element method; log-density formulation;
Entropy stability; Positivity perservation

###### 1991 Mathematics Subject Classification:

65N30, 65N12, 76S05, 76D07

This research is partially supported by NSF grant DMS-2012031.

##  1\. Introduction

The porous medium equation (PME) is given by

$\frac{\partial\rho}{\partial t}-\nabla\cdot\nabla\left(\rho^{m}\right)=0,$
(1)

where $x\in\Omega\subset\mathbb{R}^{d}$, $\rho(t,x)\geq 0$ is the unknown
density function, and $m\geq 1$. The PME can be found in many physical and
biological phenomena, including the flow of an ideal gas through porous media,
groundwater infiltration, the spread of viscous fluids, and boundary layer
theory [15, 35, 9]. The PME is an example of a non-linear parabolic PDE with
an underlying Wasserstein gradient flow structure [29] that guarantees that
any solution for initially positive data will remain so for all times, and
will be energy dissipative. Moreover, the PME has the peculiar finite speed of
propagation property which states that any compactly supported initial data
will remain so for all times, and thus, the solution at the domain interface,
referred to as the free boundary, will only propagate at a finite speed. The
PME admits a classical self-similar weak solution [4, 30] that is compactly
supported and shows this property. It is also a commonly known fact that for
certain initial data, the solution exhibits the waiting-time phenomenon in
which the free-boundary does not move, but the interior profile continues to
evolve until a certain positive finite time [1].

The existence of a free boundary due to the degeneracy of the PME makes
traditional parabolic numerical techniques ineffective. Various numerical
methods have been developed for the PME [16, 21, 38, 36, 11, 27, 12, 22, 17,
24, 13, 31, 26, 18, 8, 7]. We summarize some of the broad themes of these
methods:

  * •

Standard numerical routines that invoke standard finite element procedures for
spatial discretization and a predictor-corrector formulation for temporal
discretization (PCSFE method) can suffer from oscillations at the free
boundary which cannot be suppressed by increasing the polynomial degree or
refining the spatial mesh [38].

  * •

Discontinuous Galerkin [11], Local Discontinuous Galerkin [38], and WENO [24]
methods have been adapted to the PME to suppress these oscillations arising
near the free boundary. Particularly, [38] makes use of a non-negativity
preserving limiter, while [11] implements a maximum-principle-satisfying (MPS)
limiter. However, these schemes introduce additional numerical viscosity, due
to which they fail to accurately track the free boundary, and thus, to
accurately estimate the waiting time.

  * •

Several interface tracking schemes [26, 18, 8, 7] have been developed that
track the free boundary by solving the equation of the interface in the
Lagrangian coordinate. However, these schemes have limited applicability in
higher dimensions and to initial data with complex support due to the
complexity of implementation. Recent works in [10, 23] develop fully
Lagrangian schemes that can be applied to 2D problems. Nonetheless, being a
Lagrangian scheme, it may not handle well initial data with complex support.
If any “tangling” of the mesh occurs, the solution has to be manually
interpolated onto a new mesh.

  * •

Perturbation techniques have also been used to remove the degeneracy of the
PME [13, 31, 17]. In particular, the non-negative initial data is perturbed by
a small parameter $\varepsilon>0$ to make it positive everywhere. Using this
approach, [17] constructs a first-order scheme that is unconditionally energy
dissipative as well as provably bound preserving on structured tensor-product
grids. However, in addition to the perturbation by $\varepsilon$, this scheme
has the limitation that it cannot be applied to complex geometries.

In this work, we present two positive and energy-preserving nonlinear schemes
for the PME: one is based on a log-density formulation, and the other uses a
mixed formulation. Our log-density formulation scheme is closely related to
the work in [17] but can handle fully unstructured grids in any number of
dimensions. Additionally, both schemes do not require any perturbation of a
non-positive initial data.

  * •

The log-density based scheme uses a classical conforming piecewise linear
finite element space for the spatial discretization in combination with a
first-order semi-implicit time discretization. Global mass conservation, the
positivity of density, and energy dissipation for the energy
$\int_{\Omega}\rho(\log(\rho)-1)\,\mathrm{dx}$ are proven for this fully
discrete scheme on general unstructured meshes. Unique solvability of the
nonlinear equations in each time step is established when the mass matrix is
lumped to be a diagonal matrix. Moreover, under the condition of a Delaunay
triangulation along with a special edge-based discretization of the nonlinear
diffusion term [37], the scheme can be further proven to satisfy a discrete
maximum principle. While the theoretical results for this scheme are proven
under the condition that the initial density $\rho^{0}(x)>0$ is positive
everywhere, we show a practical implementation to handle compactly supported
initial data which simply deactivates the degrees of freedom if the associated
diagonal entry of the stiffness matrix is below a threshold value (e.g.,
$10^{-14}$).

  * •

On the other hand, the mixed method is based on a first-order reformulation of
the PME, where the unknowns are density, potential, and velocity. The
classical RT0-P0 finite element pair is used for spatial discretization, where
the velocity is approximated via the lowest order Raviart-Thomas finite
element space [32], and the density and potential are approximated by a
discontinuous piecewise constant space. A first-order semi-implicit time
discretization is then applied to this spatial discretization, which results
in a nonlinear fully discrete scheme. The fully discrete scheme is locally
mass conservative and positivity preserving under a classical CFL condition on
general unstructured meshes. Furthermore, we prove energy dissipation for the
physical energy $\int_{\Omega}\frac{\rho^{m}}{m-1}\,\mathrm{dx}$ of this
scheme when mass lumping [3] is applied to evaluate the velocity mass matrix.
Positivity of this lumped velocity mass matrix requires the triangulated mesh
to be Delaunay. We note that this mass lumping is simply a trapezoidal rule on
tensor-product meshes, on which one can further prove the solvability of the
nonlinear system and unconditional positivity preservation following similar
arguments as in the log-density formulation; see also [17]. The nonlinear
system in each time step can be efficiently solved by expressing velocity and
potential in terms of the density unknowns and then solving the parabolic
system of nonlinear equations for density alone using Newton’s method. This
scheme naturally handles a compactly supported initial density profile as zero
density/potential is allowed in the scheme.

The rest of the paper is organized as follows. In section 2, we describe the
construction of the log-density scheme and prove several properties like mass
conservation, energy stability, unique solvability, and bound preservation. In
section 3, we describe the mixed method and prove properties like local mass
dissipation, energy stability, and positivity. In section 4, we discuss
several numerical experiments and compare the results of the two schemes. In
section 5, we provide some concluding remarks.

##  2\. Log-Density Formulation

We consider the PME (1) expressed using the log-density variable
$u:=\log(\rho)$ on a bounded polyhedral domain
$\Omega\subset\mathbb{R}^{d},d=1,2,3$ with a homogeneous Neumann boundary
condition:

$\displaystyle\frac{\partial\exp(u)}{\partial
t}-\nabla\cdot\left(m\exp(m\,u)\nabla u\right)$
$\displaystyle=0,\quad{\text{in}}\quad\Omega,$ (2a)
$\displaystyle\frac{\partial u}{\partial n}$
$\displaystyle=0,\quad{\text{o}n}\quad\partial\Omega,$ (2b) and initial data
$\displaystyle u(0,x)=$ $\displaystyle\;\log(\rho^{0}(x))\text{ in }\Omega,$
(2c)

where the initial density $\rho^{0}(x)>0$ is assumed to be positive everywhere
in the domain. This log-density based formulation was first developed in [25]
for the Poisson-Nernst-Planck equations, see also [14], which is closely
related to the the entropy-stable schemes based on the entropy variables for
hyperbolic conservation laws and compressible flow in the CFD literature [19,
34, 20, 6]. In particular, the density $\rho=\exp(u)$ is guaranteed to stay
positive as long as the initial data $\rho^{0}(x)$ is positive. We note that
the recently introduced bound preserving and energy dissipative finite
difference scheme by Gu and Shen [17] is also closely related to this log-
density formulation.

Note that the PME has a special finite speed of propagation property [35]
which states that if the initial data $\rho^{0}$ has compact support, then the
solution to the Cauchy problem of the PME will also have compact support at
any other time, $t>0$. For this reason, special care has to be taken in the
case of compactly supported initial data as $u$ is negative infinity at
locations where $\rho^{0}$ is zero. In [17], the authors add a small
perturbation $\varepsilon$ to the initial data to make it positive everywhere.
We employ a different approach in which we deactivate degrees of freedom
wherever $\rho^{0}$ is close to zero; see more discussion at the end of this
section.

The PME (2) satisfy the following three important properties:

  * (i)

Mass conservation:

$\displaystyle\int_{\Omega}\rho(t,x)dx=\int_{\Omega}\rho^{0}(x)dx.$ (3a)

  * (ii)

Positivity:

If $\rho^{0}(x)>0$, then $\rho(t,x)>0$ for any $t>0$. (3b)

  * (iii)

Energy dissipation:

$\displaystyle\frac{d}{dt}E=-\int_{\Omega}m\rho^{m}|\nabla u|^{2}dx,$ (3c)

where the energy $E$ is given by

$\displaystyle
E:=\int_{\Omega}\rho(\log(\rho)-1)dx=\int_{\Omega}\exp(u)(u-1)dx.$ (3d)

Our goal is to design a numerical scheme that preserves these three
properties.

###  2.1. Spatial and Temporal Discretizations

We describe the method on a general unstructured simplicial mesh, although
quadrilateral/hexahedral meshes can also be used. Let
$\Omega_{h}:=\\{K_{i}\\}_{i=1}^{N_{K}}$ be a conforming simplicial
triangulation of the domain $\Omega$ with $N_{K}$ elements. Denote
$\mathcal{E}_{h}=\\{E_{i}\\}_{i=1}^{N_{E}}$ as the collection of $N_{E}$ edges
of $\Omega_{h}$, and $\mathcal{V}_{h}=\\{v_{i}\\}_{i=1}^{N_{V}}$ as its
collection of $N_{V}$ vertices. For any element $K\in\Omega_{h}$, denote
$\mathcal{E}_{K}:=\\{E\in\mathcal{E}_{h}:\;E\subset\bar{K}\\}$ as its edges,
and $\mathcal{V}_{K}:=\\{v\in\mathcal{V}_{h}:\;v\in\bar{K}\\}$ as its
vertices.

We shall use the $H^{1}$-conforming finite element space

$\displaystyle V_{h}:=\\{v_{h}\in H^{1}(\Omega):\;v_{h}|_{K}\in
P_{1}(K),\quad\forall K\in\Omega_{h}\\},$ (4)

where $P_{1}(K)$ is the space of linear polynomials on a simplex $K$. The
space $V_{h}$ is equipped with the standard nodal (hat) basis
$\\{\phi_{i}(x)\\}_{i=1}^{N_{V}}$ in which $\phi_{i}(v_{j})=\delta_{ij}$ where
$\delta_{ij}$ is the Kronecker delta function. Hence any function $w_{h}\in
V_{h}$ can be expressed as

$w_{h}=\sum_{i=1}^{N_{V}}w_{i}\phi_{i},$

where $\underline{w}:=[w_{1},\cdots,w_{N_{V}}]^{\prime}$ is the coefficient
vector satisfying $w_{i}=w_{h}(v_{i})$.

The spatial discretization for (2) then reads: find $u_{h}\in V_{h}$ such
that, for $t>0$,

$\displaystyle\int_{\Omega}\frac{\partial\exp(u_{h})}{\partial
t}v_{h}\,\mathrm{dx}+\int_{\Omega}m\exp(m\,u_{h})\nabla u_{h}\cdot\nabla
v_{h}\,\mathrm{dx}=$ $\displaystyle 0,\quad\forall v_{h}\in V_{h},$ (5)

with initial conditions

$u_{h}(0,v_{i})=\log\left(\rho^{0}_{h}(v_{i})\right),\quad\forall
v_{i}\in\mathcal{V}_{h}.$ (6)

We apply a first-order semi-implicit discretization for the ODE system (5) to
arrive at the following fully discrete scheme: Given data $u_{h}^{n-1}\in
V_{h}$ at time $t^{n-1}$ and time step size $\Delta t$, find $u_{h}^{n}\in
V_{h}$ at time $t^{n}=t^{n-1}+\Delta t$ such that

$\mathcal{M}\left(\frac{\exp(u_{h}^{n})-\exp(u_{h}^{n$-1$})}{\Delta
t},v_{h}\right)+\mathcal{A}\left(m\exp(mu_{h}^{n$-1$});\nabla u_{h}^{n},\nabla
v_{h}\right)=0,\quad\forall v_{h}\in V_{h},$ (7)

where the mass operator $\mathcal{M}$ and stiffness operator $\mathcal{A}$
read as follows:

$\displaystyle\mathcal{M}(\alpha,\beta):=$
$\displaystyle\;\int_{\Omega}\alpha\cdot\beta\mathrm{dx},\quad\quad\mathcal{A}(\gamma;\alpha,\beta):=\;\int_{\Omega}\gamma\nabla\alpha\cdot\nabla\beta\mathrm{dx}.$

Unique solvability of the scheme requires the mass matrix to be mass lumped,
which is achieved by applying the following vertex-based quadrature rule

$\mathcal{M}_{h}(\alpha,\beta):=\sum_{K\in\Omega_{h}}\sum_{v\in\mathcal{V}_{K}}\frac{|K|}{d+1}\alpha(v)\beta(v)=\sum_{i=1}^{N_{v}}\frac{|S_{i}|}{d+1}\alpha(v_{i})\beta(v_{i}),$
(8)

where $S_{i}:=\cup_{v_{i}\in\bar{K}}\\{\bar{K}\\}$ is the vertex patch of
$v_{i}$ and $|S_{i}|$ is its volume.

Furthermore, we make use of the following edge-based integration formula [37]
for the stiffness matrix, which will be used to prove uniform the boundedness
of the scheme:

$\begin{split}\mathcal{A}_{h}(\gamma;\alpha,\beta):=&\;\sum_{K\in\Omega_{h}}\sum_{E\in\mathcal{E}_{K}}\omega^{K}_{E}\tilde{\gamma}_{E}\delta_{E}(\alpha)\delta_{E}(\beta)\\\
=&\;\sum_{E\in\mathcal{E}_{h}}\omega_{E}{\tilde{\gamma}}_{E}\delta_{E}(u_{h})\delta_{E}(v_{h}),\end{split}$
(9)

where $\omega_{E}:=\sum_{K\supset E}\omega^{K}_{E}$. Here for an edge $E$ with
vertices $v_{i}$ and $v_{j}$, we have

$\delta_{E}(u_{h})=u_{h}(v_{i})-u_{h}(v_{j}).$ (10)

The quantity $\tilde{\gamma}_{E}$ is the following harmonic average on $E$,

$\tilde{\gamma}_{E}=\left[\frac{1}{|E|}\int_{E}\frac{1}{\gamma}\,ds\right]^{-1}.$
(11)

Also, the weights $\omega^{K}_{E}$ are given by the identity [37, 5]:

$\omega_{E}^{K}=\frac{1}{d(d-1)}\left|\kappa_{E}^{K}\right|\cot\theta_{E}^{K},$
(12)

with $d$ being the number of dimensions, $\theta_{E}^{K}$ being the angle
between faces not containing edge $E$, and $\kappa_{E}^{K}$ the ($d-2$)
dimensional simplex formed by their intersection.

###  2.2. Properties

In this section, we prove several important results for the scheme (7). We
begin by showing that the scheme satisfies discrete versions of the properties
in (3).

######  Theorem 2.1.

The fully discrete, semi-implicit scheme (7) conserves mass, preserves
positivity, and dissipates energy in the following form

$E_{h}^{n}-E_{h}^{n$-1$}\leq-\int_{\Omega_{h}}m\Delta
t\exp\left(mu_{h}^{n$-1$}\right)|\nabla u_{h}^{n}|^{2}dx,$ (13)

where the energy $E_{h}^{n}:=\mathcal{M}(\exp(u_{h}^{n}),u_{h}^{n}-1).$

###### Proof.

The positivity of the scheme is guaranteed by the log-density variable
formulation since $\exp(u_{h}^{n})>0$ always. The mass conservation can be
proved by picking $v_{h}=1$ in (7):

$\mathcal{M}\left(\frac{\exp(u_{h}^{n})-\exp(u_{h}^{n$-1$})}{\Delta
t},1\right)+0=0\implies\mathcal{M}\left(\exp(u_{h}^{n}),1\right)=\mathcal{M}\left(\exp(u_{h}^{n$-1$}),1\right).$
(14)

The only non-trivial property to prove is energy dissipation. Observe that by
Taylor expansion, we get

$(\exp(a)-\exp(b))a=\exp(a)(a-1)-\exp(b)(b-1)+\frac{1}{2}\exp(\xi)(a-b)^{2},$
(15)

where $\xi$ is a function between $a$ and $b$. Now, picking $v=u_{h}^{n}$ and
using the above Taylor expansion, we get

$\begin{split}\mathcal{M}\left(\exp(u_{h}^{n}),u_{h}^{n}-1\right)-\mathcal{M}\left(\exp(u_{h}^{n$-1$}),u_{h}^{n$-1$}-1\right)=&-\Delta
t\mathcal{A}\left(m\exp(mu_{h}^{n$-1$}),\nabla u_{h}^{n},\nabla
u_{h}^{n}\right);\\\
&\quad-\mathcal{M}\left(\exp(\xi),\frac{(u_{h}^{n}-u_{h}^{n$-1$})^{2}}{2}\right)\end{split}$
(16)

where $\xi$ is a function between $u_{h}^{n}$ and $u_{h}^{n$-1$}$. This
completes the proof. ∎

Next we prove unique solvability of the scheme (7) where mass lumping (8) is
used for the mass operator.

######  Theorem 2.2.

The fully discrete scheme (7) is uniquely solvable provided that mass lumping
(8) is used to evaluate the mass operator.

###### Proof.

We prove this result using matrix-vector notation. Denoting
$\underline{u}^{j}$ as the coefficient vector of solution $u_{h}^{j}\in
V_{h}$, the scheme (7) with mass lumping can then be written in the following
matrix-vector form:

$\mathbf{M}(\exp(\underline{u}^{n})-\exp(\underline{u}^{n$-1$}))+\Delta
t\mathbf{A}^{n$-1$}\underline{u}^{n}=0,$ (17)

where $\mathbf{M}$ is the diagonal mass matrix, by virtue of mass-lumping,
with entries $\mathbf{M}_{ii}=\frac{|S_{i}|}{(d+1)},$ and $\mathbf{A}^{n$-1$}$
is a symmetric positive semidefinite stiffness matrix with entries:

$\mathbf{A}^{n$-1$}_{ij}=\mathcal{A}\left(m\exp(mu_{h}^{n$-1$});\nabla\phi_{i},\nabla\phi_{j}\right).$
(18)

It is clear that the nonlinear system (17) is the Euler-Lagrange equation of
the following minimization problem:

$\underline{u}^{n}:=\mathrm{argmin}_{\underline{u}\in\mathbb{R}^{N_{V}}}F(\underline{u}),$

where the energy functional is

$F(\underline{u})=\underline{1}\cdot\mathbf{M}\exp(\underline{u})-\underline{u}\cdot\mathbf{M}\exp(\underline{u}^{n$-1$})+\frac{\Delta
t}{2}\underline{u}\cdot\mathbf{A}^{n$-1$}\underline{u},$ (19)

where $\underline{1}$ is the vector of ones of size $N_{V}$. Hence, unique
solvability of (17) is equivalent to the coercivity (existence) and strictly
convexity (uniqueness) of the functional (19). Both properties can be easily
verified using the positivity of the mass and stiffness matrices and
elementary calculation. We leave out the details.

∎

Finally, we prove the uniform boundedness of the scheme (7) when mass lumping
(8) is used for the mass operator, and the edge-based integration (9) is used
for the stiffness operator.

######  Theorem 2.3.

If the triangulation $\Omega_{h}$ is Delaunay, then the solution to the fully
discrete scheme (7) with mass lumping (8) being used for the mass operator and
edge-based integration (9) being used for the stiffness operator is uniformly
bounded. That is, given $0<\varepsilon_{1}<\varepsilon_{2}$ such that
$\varepsilon_{1}\leq\exp(u_{h}^{n$-1$})\leq\varepsilon_{2}$, we have
$\varepsilon_{1}\leq\exp(u_{h}^{n})\leq\varepsilon_{2}$.

###### Proof.

We only prove the lower bound, i.e., if
$\exp(u_{h}^{n$-1$})\geq\varepsilon_{1}$, then
$\exp(u_{h}^{n})\geq\varepsilon_{1}$, as the upper bound use the same
argument.

Taking a non-negative test function $v_{h}\in V_{h}$ in (7) such that its
coefficient $v_{i}=\max\\{\epsilon_{1}-\exp(u_{i}^{n}),0\\}$, and using the
matrix-vector form (17), we get

$\displaystyle\underline{v}\cdot\mathbf{M}(\exp(\underline{u}^{n})-\exp(\underline{u}^{n$-1$}))+\Delta
t\underline{v}\cdot\mathbf{A}^{n$-1$}\underline{u}^{n}=0.$ (20)

Using the fact the $\mathbf{M}$ is a diagonal positive matrix,
$\exp(\underline{u}^{n$-1$})-\varepsilon_{1}\geq 0$, and definition of
$\underline{v}$, we have

$\underline{v}\cdot\mathbf{M}(\exp(\underline{u}^{n})-\exp(\underline{u}^{n$-1$}))=\underbrace{\underline{v}\cdot\mathbf{M}(\exp(\underline{u}^{n})-\varepsilon_{1})}_{\leq
0}-\underbrace{\underline{v}\cdot\mathbf{M}(\exp(\underline{u}^{n$-1$}))-\varepsilon_{1})}_{\geq
0}\leq 0.$

Next, using the edge integration formula (9), we have

$\displaystyle\underline{v}\cdot\mathbf{A}^{n$-1$}\underline{u}^{n}=\sum_{E_{ij}\in\mathcal{E}_{h}}\omega_{E_{ij}}\tilde{\gamma}_{E_{ij}}(v_{i}-v_{j})(u_{i}-u_{j}),$
(21)

where $\gamma=m\exp(u_{h}^{n$-1$})$. By definition of $v_{h}$, we have

$(v_{i}-v_{j})(u_{i}-u_{j})=(v_{i}-v_{j})((u_{i}-\varepsilon_{1})-(u_{j}-\varepsilon_{1}))\leq
0.$

Hence, the term (21) is non-positive as long as $\omega_{E_{ij}}>0$ for all
$E_{ij}\in\mathcal{E}_{h}$, which is equivalent to the requirement that the
triangulation $\Omega_{h}$ is Delaunay; see [37]. In this case, we have

$\underline{v}\cdot\mathbf{M}(\exp(\underline{u}^{n})-\varepsilon_{1})=\sum_{i=1}^{N_{V}}\mathbf{M}_{ii}(\exp({u}_{i}^{n})-\varepsilon_{1})\max\\{\epsilon_{1}-\exp(u_{i}^{n}),0\\}=0$

thanks to (20). Hence, $\exp(u_{i}^{n})\geq\varepsilon_{1}$, which completes
the proof.

∎

######  Remark 2.1.

We remark that the edge-based integration in Theorem 2.3 above is only used to
prove the non-positivity of the term (21) on Delaunay triangulations. When the
mesh is a structured tensor-product grid, such bound preservation can be
easily proven with other standard numerical integration rules; see, e.g.,
[17].

Moreover, although we use the edge integration formula (9) to prove the
uniform boundedness result in Theorem 2.3, in practical implementation, we
simply use the mass-lumping quadrature (8) to compute both $\mathbf{M}$ and
$\mathbf{A}^{n$-1$}$. The resulting scheme is still very robust, and is
provable unique solvable and satisfy the properties in Theorem 2.1.

######  Remark 2.2.

We conclude this section by remarking on the practical implementation of the
scheme (7) with compactly supported initial data $\rho^{0}$.

In this case, we still interpolate the initial data using (6). Note that we
have $u_{i}^{0}=-\infty$ whenever $\rho^{0}(v_{i})=0$. Now, the $k$-th Newton
iteration for (17) takes the form

$\left(\mathbf{M}\mathbf{D}^{(k-1)}+\Delta
t\mathbf{A}^{n$-1$}\right)\underline{u}^{(k)}=\mathbf{M}\left(\mathbf{D}^{(k-1)}\underline{u}^{(k-1)}-\exp{(\underline{u}^{(k-1)})}+\exp{(\underline{u}^{n-1})}\right),$
(22)

where $\mathbf{D}^{(k-1)}$ is a diagonal matrix with diagonal entry
$\mathbf{D}^{(k-1)}_{ii}=\exp{u_{i}^{(k-1)}}$, and
$\underline{u}^{(0)}=\underline{u}^{n$-1$}$. To solve the above linear system
for $\underline{u}^{(k)}$, we only activate the $i$-th degree of freedom if
the diagonal entry $\left(\mathbf{M}\mathbf{D}^{(k-1)}+\Delta
t\mathbf{A}^{n$-1$}\right)_{ii}$ is greater than a small cutoff value (e.g.,
$10^{-14}$), and set the inactive degrees of freedom to $-\infty$.

##  3\. Mixed Method

In this section, we develop our mixed formulation to solve the PME (1).

The physical energy for the PME is given by $U(\rho)=\frac{\rho^{m}}{m-1}$. We
define a new potential variable equal to the derivative of the physical
energy:

$\mu=U^{\prime}(\rho)=\frac{m}{m-1}\rho^{m-1}.$ (23)

We also define a velocity variable, $\mathbf{u}$, and set it equal to the
negative gradient of the potential, i.e.,

$\mathbf{u}=-\nabla\mu.$ (24)

Notice that using this definition, we easily observe that

$\nabla\rho^{m}=\rho\nabla\mu=-\rho\mathbf{u}.$ (25)

We have now successfully reformulated the PME (1) into the following first-
order system

$\begin{cases}&\rho_{t}+\nabla\cdot(\rho\mathbf{u})=0,\\\
&\mathbf{u}+\nabla\mu=0,\\\ &\mu=\frac{m}{m-1}\rho^{m-1},\end{cases}$ (26)

where the unknowns are density, potential, and velocity. We again equip the
PME system (26) with the homogeneous Neumann boundary condition

$\mathbf{u}\cdot\mathbf{n}=0\text{ on }\partial\Omega.$

###  3.1. Spatial and Temporal Discretizations

Our energy-stable mixed method can be constructed on structured tensor product
meshes in any space dimension, or on general unstructured Delaunay triangular
meshes in two dimensions. Due to the use of velocity mass lumping as a key
tool to establish the energy stability result, the method does not work on
general simplicial meshes in three dimensions. Below we describe the method in
detail on 2D triangular meshes, following the meshing notation as in Section
2.1. For any triangular element $K\in\Omega_{h}$, we denote $\partial K$ as
its boundary whose associated outward unit normal is $\mathbf{n}_{K}$.

We shall use the following two finite element spaces:

$\displaystyle Q_{h}=\left\\{q_{h}\in
L^{2}(\Omega):\enspace\left.q_{h}\right|_{K}\in P_{0}(K),\enspace\forall
K\in\Omega_{h}\right\\},$ (27)
$\displaystyle\mathbf{V}_{h}=\left\\{\mathbf{v}_{h}\in\operatorname{H}(\operatorname{div};\Omega):\enspace\left.\mathbf{v}_{h}\right|_{K}\in
RT_{0}(K),\enspace\forall
K\in\Omega_{h},\;\;\mathbf{v}_{h}\cdot\mathbf{n}=0,\quad\text{ on
}\partial\Omega\right\\},$

where $P_{0}(K)$ is the constant space, and
$RT_{0}(K)=[P_{0}(K)]^{2}\oplus\mathbf{x}P_{0}(K)$ is the lowest-order
Raviart-Thomas space on the triangle $K$. Note that the Neumann boundary
condition is encoded in the velocity space $\mathbf{V}_{h}$.

The spatial discretization now reads as follows: find $\rho_{h},\mu_{h}\in
Q_{h}$ and $\mathbf{u}_{h}\in\mathbf{V}_{h}$ such that, for $t>0$,

$\displaystyle\sum_{K\in\Omega_{h}}\left[\int_{K}(\rho_{h})_{t}q_{h}\,\mathrm{dx}+\int_{\partial
K}\hat{\rho}_{h}\mathbf{u}_{h}\cdot\mathbf{n}_{K}q_{h}\,\mathrm{ds}\right]=$
$\displaystyle\;0,\qquad\forall\enspace q_{h}\in Q_{h},$ (28a)
$\displaystyle\sum_{K\in\Omega_{h}}\left[\int_{K}\mathbf{u}_{h}\cdot\mathbf{v}_{h}\,\mathrm{dx}-\int_{\partial
K}\mu_{h}\mathbf{v}_{h}\cdot\mathbf{n}_{K}\,\mathrm{ds}\right]=$
$\displaystyle\;0,\qquad\forall\enspace\mathbf{v}_{h}\in\mathbf{V}_{h},$ (28b)
$\displaystyle\sum_{K\in\Omega_{h}}\int_{K}(\mu_{h}-\tfrac{m}{m-1}\rho_{h}^{m-1})r_{h}dx=$
$\displaystyle\;0,\qquad\forall\enspace r_{h}\in Q_{h},$ (28c)

where $\hat{\rho}_{h}$ is the upwinding numerical flux, i.e.,

given an edge $E$ shared by two elements $K^{+}$ and $K^{-}$ with
$\mathbf{u}_{h}\cdot\mathbf{n}_{K^{-}}|_{E}\geq 0$, $\hat{\rho}_{h}$ takes
value from $K^{-}$:

$\left.\hat{\rho}_{h}\right|_{E}=(\rho_{h}|_{K^{-}})|_{E},\quad\text{ where
}\mathbf{u}_{h}\cdot\mathbf{n}_{K^{-}}|_{E}\geq 0.$ (29)

Analogous to the log-density approach, we apply a first-order semi-implicit
discretization of the ODE system (28) to obtain the following fully discrete
scheme: Given time step size $\Delta t$ and data $\rho_{h}^{n$-1$}\in Q_{h}$
at time $t^{n$-1$}$, find $\mu^{n}_{h},\rho^{n}_{h}\in Q_{h}$, and
$\mathbf{u}_{h}^{n}\in\mathbf{V}_{h}$, at time $t^{n}$, such that

$\displaystyle\mathcal{M}_{h}\left(\frac{\rho_{h}^{n}-\rho_{h}^{n-1}}{\Delta
t},q_{h}\right)+\mathcal{A}_{h}\left(\hat{\rho}_{h}^{n$-1$};\mathbf{u}_{h}^{n},q_{h}\right)=0,\quad$
$\displaystyle\forall q_{h}\in Q_{h},$ (30a)
$\displaystyle\bar{\mathcal{M}_{h}}(\mathbf{u}_{h}^{n},\mathbf{v}_{h})-\bar{\mathcal{A}_{h}}(\mu^{n}_{h},\mathbf{v}_{h})=0,\quad$
$\displaystyle\forall\mathbf{v}_{h}\in\mathbf{V}_{h},$ (30b)
$\displaystyle\mathcal{M}_{h}\left(\mu^{n}_{h}-\frac{m}{m-1}(\rho_{h}^{n})^{m-1},r_{h}\right)=0,\quad$
$\displaystyle\forall r_{h}\in Q_{h},$ (30c)

where the associated operators are defined as:

$\displaystyle\mathcal{M}_{h}\left(\alpha,\beta\right)=\sum_{K\in\Omega_{h}}\int_{K}\alpha\beta\mathrm{dx},$
$\displaystyle\bar{\mathcal{M}}_{h}(\bm{\alpha},\bm{\beta})=\sum_{K\in\Omega_{h}}\int_{K}\bm{\alpha}\cdot\bm{\beta}\mathrm{dx},$
(31a)
$\displaystyle\mathcal{A}_{h}\left(\gamma;\bm{\alpha},\beta\right)=\sum_{K\in\Omega_{h}}\int_{\partial
K}\gamma\bm{\alpha}\cdot\mathbf{n}_{K}\beta\mathrm{dx},$
$\displaystyle\bar{\mathcal{A}}_{h}\left(\alpha,\bm{\beta}\right)=\sum_{K\in\Omega_{h}}\int_{\partial
K}\alpha\bm{\beta}\cdot\mathbf{n}_{K}\mathrm{dx}.$ (31b)

Again, the numerical flux $\hat{\rho}_{h}^{n$-1$}$ in (30a) is taken to be the
upwind flux: given an edge $E$ shared by two elements $K^{+}$ and $K^{-}$ with
$\mathbf{u}_{h}^{n}\cdot\mathbf{n}_{K^{-}}|_{E}\geq 0$, we take

$\left.\hat{\rho}_{h}^{n-1}\right|_{E}=(\rho_{h}^{n-1}|_{K^{-}})|_{E},\quad\text{
where }\mathbf{u}_{h}^{n}\cdot\mathbf{n}_{K^{-}}|_{E}\geq 0,$ (32)

Standard one-point numerical integration rules are used in the evaluation of
the operators $\mathcal{M}_{h},\mathcal{A}_{h},$ and $\bar{\mathcal{A}}_{h}$.
However, the energy stability of the scheme requires that the velocity mass
matrix be mass lumped. As such, the operator $\bar{\mathcal{M}}_{h}$ must be
evaluated using an appropriate mass-lumping quadrature. On tensor-product
grids, this is readily achieved by the use of the trapezoidal rule. On
triangular meshes, mass-lumping is achieved by using the following formula
given in [3]:

$\bar{\mathcal{M}_{h}}\left(\bm{\alpha},\bm{\beta}\right):=\sum_{K\in\Omega_{h}}\sum_{E\in\mathcal{E}_{K}}\omega_{E}^{K}\varphi_{E}(\bm{\alpha})\varphi_{E}(\bm{\beta})=\sum_{E\in\mathcal{E}_{h}}\omega_{E}\varphi_{E}(\bm{\alpha})\varphi_{E}(\bm{\beta}),$
(33)

where $\omega_{E}:=\sum_{K\supset E}\omega_{E}^{K}$, and
$\varphi_{E}(\bm{\alpha}):=\bm{\alpha}\cdot\mathbf{n}_{E}$ denotes the normal
flux of $\bm{\alpha}$ through edge $E$. The weights $\omega_{E}^{K}$ are given
by

$\omega_{E}^{K}=\frac{1}{2}\cot\theta_{E}^{K},$ (34)

where $\theta_{E}^{K}$ is the angle opposite to edge $E$ in $K$. It is clear
that the mass matrix associated with the integration rule in (33) is a
diagonal matrix, whose diagonal entries are positive provided that the mesh is
Delaunay.

To efficiently solve the scheme (30), we first apply static condensation to
locally solve the potential and velocity variables $\mu_{h}^{n}$ and
$\mathbf{u}_{h}^{n}$ as functions of the density $\rho_{h}^{n}$ using
equations (30b)–(30c), and then use Newton’s method to solve the resulting
parabolic nonlinear system (30a) for density alone.

###  3.2. Properties

######  Theorem 3.1.

Provided a Delaunay triangulation $\Omega_{h}$, the fully discrete scheme (30)
is mass conservative and energy dissipative in the following forms

$\displaystyle\mathcal{M}_{h}\left(\rho_{h}^{n},1\right)$
$\displaystyle=\mathcal{M}_{h}\left(\rho_{h}^{n$-1$},1\right)$ (35a)
$\displaystyle\mathcal{M}_{h}\left(U(\rho^{n}),1\right)-\mathcal{M}_{h}\left(U(\rho^{n$-1$}),1\right)$
$\displaystyle\leq-\Delta
t\sum_{E\in\mathcal{E}_{h}}\omega_{E}\varphi_{E}(\mathbf{\alpha})\varphi_{E}(\mathbf{u}_{h}^{n})^{2},$
(35b)

where $U^{n}=\frac{m}{m-1}(\rho^{n})^{m-1}$ is the physical energy, and the
right-hand side in (35b) comes from the quadrature formula (33).

###### Proof.

By picking $q_{h}=1$ in (30a) and using the homegenous boundary condition
$\mathbf{u}_{h}^{n}\cdot\mathbf{n}|_{\partial\Omega}=0$, we can easily recover
the mass conservation (35a). Now, to prove energy stability, we take the test
function $\mathbf{v}_{h}$ in (30b) to be a function in $\mathbf{V}_{h}$ such
that its normal flux through edge $E$ is given by
$\varphi_{E}(\mathbf{v}_{h})=\varphi_{E}(\hat{\rho}_{h}^{n$-1$}\mathbf{u}_{h}^{n})$.
For this choice of $\mathbf{v}_{h}$, we have

$\mathcal{A}_{h}\left(\hat{\rho}_{h}^{n$-1$};\mathbf{u}_{h}^{n},\mu^{n}_{h}\right)=\bar{\mathcal{A}_{h}}(\mu^{n}_{h},\mathbf{v}_{h}).$
(36)

Additionally, through an application of (33), we obtain

$\bar{\mathcal{M}_{h}}(\mathbf{u}_{h},\mathbf{v}_{h})=\sum_{E\in\mathcal{E}_{h}}\omega_{E}\hat{\rho}_{h}^{n$-1$}\varphi_{E}(\mathbf{u}_{h}^{n})^{2}$
(37)

Next taking test function $r_{h}=\frac{\rho_{h}^{n}-\rho_{h}^{n$-1$}}{\Delta
t}$ in (30c), we get

$\begin{split}&\mathcal{M}_{h}\left(\frac{\rho_{h}^{n}-\rho_{h}^{n-1}}{\Delta
t},\mu_{h}^{n}\right)=\mathcal{M}_{h}\left(\frac{\rho_{h}^{n}-\rho_{h}^{n-1}}{\Delta
t},U^{\prime}(\rho_{h}^{n})\right)\\\ &=\frac{1}{\Delta
t}\mathcal{M}_{h}\left(U(\rho_{h}^{n}),1\right)-\frac{1}{\Delta
t}\mathcal{M}_{h}\left(U(\rho_{h}^{n$-1$}),1\right)+\frac{1}{\Delta
t}\mathcal{M}_{h}\left(\frac{U^{\prime\prime}(\xi)}{2}(\rho_{h}^{n}-\rho_{h}^{n$-1$})^{2}\right),\end{split}$
(38)

where we have used Taylor expansion with $\xi$ being a function between
$\rho^{n}_{h}$ and $\rho_{h}^{n$-1$}$. Finally, taking the test function
$q_{h}=\mu_{h}^{n}$ in (30a), and using the above relations, we get

$\begin{split}\mathcal{M}_{h}\left(U(\rho^{n}),1\right)-\mathcal{M}_{h}\left(U(\rho^{n$-1$}),1\right)&=-\Delta
t\sum_{E\in\mathcal{E}_{h}}\omega_{E}\hat{\rho}_{h}^{n$-1$}\varphi_{E}(\mathbf{u}_{h}^{n})^{2}-\mathcal{M}_{h}\left(\frac{U^{\prime\prime}(\xi)}{2}(\rho_{h}^{n}-\rho_{h}^{n$-1$})^{2}\right)\\\
&\leq-\Delta
t\sum_{E\in\mathcal{E}_{h}}\omega_{E}\hat{\rho}_{h}^{n$-1$}\varphi_{E}(\mathbf{u}_{h}^{n})^{2}\end{split}$
(39)

The right-hand side is non-positive if $\omega_{E}\geq 0$ which is guaranteed
if the triangulation $\Omega_{h}$ is Delaunay. This completes the proof. ∎

Next, we prove that the scheme (30) is positivity preserving under a usual CFL
time stepping constraint.

######  Theorem 3.2.

Given $\rho^{n$-1$}_{h}\geq 0$, the fully discrete scheme (30) is positivity
preserving under the CFL condition,

${\Delta
t}\sum_{E\in\mathcal{E}^{-}_{K}}|\mathbf{u}_{h}^{n}\cdot{\mathbf{n}_{K}}|_{E}\frac{|E|}{|K|}\leq
1,\quad\forall K\in\Omega_{h},$ (40)

where

$\mathcal{E}_{K}^{-}:=\\{E\in\mathcal{E}_{K}:\quad\mathbf{u}_{h}^{n}\cdot\mathbf{n}_{K}|_{E}\geq
0.\\}$

###### Proof.

Restricting the mass conservation equation (30a) to a single element
$K\in\Omega_{h}$, we have

$\frac{|K|}{\Delta
t}(\rho_{K}^{n}-\rho_{K}^{n$-1$})+\sum_{E\in\mathcal{E}_{K}}\hat{\rho}^{n$-1$}_{E}\left.(\mathbf{u}^{n}_{h}\cdot\mathbf{n}_{K})\right|_{E}|E|=0,$
(41)

where $\rho_{K}^{n}:=\rho_{h}^{n}|_{K}$ is the restriction of the function
$\rho_{h}^{n}\in Q_{h}$ to the element $K$. This implies that

$\displaystyle\rho_{K}^{n}=\rho_{K}^{n$-1$}-\frac{\Delta
t}{|K|}\sum_{E\in\mathcal{E}_{K}}\hat{\rho}^{n$-1$}_{E}\left.(\mathbf{u}^{n}_{h}\cdot\mathbf{n}_{K})\right|_{E}|E|,$

By the definition of $\mathcal{E}^{-}_{K}$, we have
$\hat{\rho}_{E}^{n-1}=\rho_{K}^{n-1}$ for all $E\in\mathcal{E}^{-}_{K}$.
Hence,

$\displaystyle\rho_{K}^{n}=\rho_{K}^{n$-1$}(1-\frac{\Delta
t}{|K|}\sum_{E\in\mathcal{E}^{-}_{K}}\left|\mathbf{u}^{n}_{h}\cdot\mathbf{n}_{K}\right|_{E}|E|)+\frac{\Delta
t}{|K|}\sum_{E\in\mathcal{E}_{K}\backslash\mathcal{E}^{-}_{K}}\hat{\rho}_{E}^{n-1}\left|\mathbf{u}^{n}_{h}\cdot\mathbf{n}_{K}\right|_{E}|E|,$

Both terms on the right hand side are nonnegative under the assumption that
$\rho_{h}^{n-1}\geq 0$ and (40). This completes the proof. ∎

######  Remark 3.1.

We note that for tensor product meshes, the stronger property of unconditional
positivity preservation can be proven following similar arguments as in the
proof of Theorem 2.3; see also [17]. The key is to locally eliminate velocity
and potential degrees of freedom to express the scheme as a finite volume
scheme for the piecewise constant density unknown only. In this case,
unconditionally positivity preservation holds for any consistent numerical
flux, that is, the upwinding flux is not needed for the positivity proof on
tensor-product meshes. We leave out the detailed derivation.

##  4\. Numerical Results

In this section, we present out numerical findings. All computations are
performed using the Python interface of the open-source library NGSolve [33].
The source code for these computations is available at the following git
repository: https://github.com/avj-jpg/pme.

|  | Log-density method | Mixed method
---|---|---|---
m | N | $\Delta t$ |  | Error in
---
$[-5,5]$
Order |  | Error in
---
$[-10,10]$
Order | $\Delta t$ |  | Error in
---
$[-5,5]$
Order |  | Error in
---
$[-10,10]$
Order
2 | 100 | 1/5 | 1.19e-01 | - | 3.72e-01 | - | 1/10 | 4.53e-02 | - | 8.48e-02 | -
| 200 | 1/20 | 3.04e-02 | 1.971 | 1.03e-01 | 1.860 | 1/20 | 2.27e-02 | 0.999 | 4.26e-02 | 0.992
| 400 | 1/80 | 7.57e-03 | 2.009 | 2.49e-02 | 2.043 | 1/40 | 1.13e-02 | 0.999 | 2.14e-02 | 0.994
| 800 | 1/320 | 1.88e-03 | 2.007 | 6.18e-03 | 2.010 | 1/80 | 5.67e-03 | 1.000 | 1.08e-02 | 0.992
3 | 100 | 1/5 | 6.63e-02 | - | 3.15e-01 | - | 1/10 | 1.94e-02 | - | 8.61e-02 | -
| 200 | 1/20 | 1.64e-02 | 2.012 | 9.36e-02 | 1.752 | 1/20 | 9.76e-03 | 0.995 | 4.69e-02 | 0.876
| 400 | 1/80 | 3.97e-03 | 2.050 | 2.55e-02 | 1.878 | 1/40 | 4.89e-03 | 0.997 | 2.65e-02 | 0.825
| 800 | 1/320 | 9.31e-04 | 2.091 | 8.93e-03 | 1.511 | 1/80 | 2.45e-03 | 0.998 | 1.54e-02 | 0.782
4 | 100 | 1/5 | 4.07e-02 | - | 2.33e-01 | - | 1/10 | 1.19e-02 | - | 1.05e-01 | -
| 200 | 1/20 | 1.02e-02 | 1.994 | 7.51e-02 | 1.634 | 1/20 | 5.95e-03 | 1.006 | 6.03e-02 | 0.793
| 400 | 1/80 | 2.47e-03 | 2.050 | 2.70e-02 | 1.478 | 1/40 | 3.01e-03 | 0.980 | 3.25e-02 | 0.892
| 800 | 1/320 | 4.43e-04 | 2.479 | 1.24e-02 | 1.126 | 1/80 | 1.49e-03 | 1.016 | 2.14e-02 | 0.606

Table 1. Convergence results for 1D Barenblatt initial data.

(a) Log density method.

(b) Mixed method.

Figure 1. Evolution of 1D Barenblatt initial data by the two schemes for
$m=3$, $\Delta t=0.05$, and $N=200$ elements.

###  4.1. 1D Barenblatt solution

The Porous Medium Equation admits an exact weak solution formulated by
Barenblatt [4] and Pattle [30]. In the one-dimensional case, the Barenblatt
solution is given by the equation:

$\rho_{B}(x,t)=(t+1)^{-k}\left(s_{0}-\frac{k(m-1)}{2m}\frac{x^{2}}{(t+1)^{2k}}\right)_{+}^{\frac{1}{m-1}},\quad
t>0,$ (42)

where $k=(m+1)^{-1}$, and $s_{0}$ denotes a scaling factor. Note that this
data is compactly supported in the interval $[-\eta_{m}(t),\eta_{m}(t)]$,
where the right boundary $\eta_{m}(t)$ moves as:

$\eta_{m}(t)=\sqrt{\frac{2m}{k(m-1)}}(t+1)^{k}.$ (43)

To verify the accuracy of the two schemes, we use (42) as the initial data,
with $x\in[-10,10]$, $s_{0}=3$, and $m\in\\{2,3,4\\}$, and conduct a spacetime
mesh-refinement study for the $L^{2}$-error of density at final time $t=1$. We
record the $L^{2}$-error in the entire domain $[-10,10]$ and in the interval
$[-5,5]$ away from the interface. We utilize a sequence of meshes consisting
of $\\{100\cdot 2^{i}\\}_{i=0}^{3}$ spatial elements and set the timestep size
to $\frac{1}{5\cdot 4^{i}}$ for the log-density method and to
$\frac{1}{10\cdot 2^{i}}$ for the mixed method, correspondingly. The results
of this convergence study are recorded in table 1. For all three values of
$m$, we observe that the log-density method is second-order accurate in space
and first-order accurate in time in the region $[-5,5]$ where the solution
remains smooth at final time. For larger values of $m$, the order of
convergence deteriorates when the $L^{2}$-error calculated in the entire
domain. This is anticipated since the error is expected to be greater at the
interface for larger values of $m$, owing to the decreased regularity of the
solution.

On the other hand, we observe that the mixed scheme has the expected first
order accuracy in both space and time in the interval $[-5,5]$. Similar to the
case of the log-density method, the order of convergence is observed to decay
with $m$ if the $L^{2}$-error is measured in the entire domain $[-10,10]$.
However, the decay in order is slower than in the case of the log-density
scheme.

Figure 1 displays plots of the numerically computed Barenblatt solution using
the two schemes at time $t=1$, for $m=3$ and $\Delta t=0.05$. These plots
illustrate that the initial profile is accurately evolved by both schemes,
with no oscillations emerging near the interface.

(a) Log-density method

(b) Mixed method

(c) Log-density method

(d) Mixed method

Figure 2. Observation of waiting time phenomenon for $m=3$ and $\theta=0$.
Panels (A) and (C) show the solutions at various times, and panels (B) and (D)
show the density over time at the node ($x_{R}^{0}$) corresponding to the
initial location of the right interface.

###  4.2. Waiting time phenomenon

The waiting time phenomenon is a well-known feature of the PME for initial
data of form

$\rho^{0}(x)=\begin{cases}\left(\frac{m-1}{m}\left((1-\theta)\cos^{2}x+\theta\cos^{4}x\right)\right)^{\frac{1}{m-1}},&\text{
if }-\frac{\pi}{2}\leq x\leq\frac{\pi}{2}\\\ 0,\quad&\text{ otherwise,
}\end{cases}$ (44)

with $\theta$ in the interval $[0,1]$. For this data, the interface does not
move until $t$ is greater than a certain waiting time $t^{*}$, even though the
internal profile continues to evolve. The theoretical waiting time for
$0\leq\theta\leq\frac{1}{4}$ is given by [1]

$t^{*}=\frac{1}{2(m+1)(1-\theta)}.$ (45)

In our numerical experiment, we set the parameters $m=3$ and $\theta=0$, and
solve the equation (44) using the two schemes on a sequence of meshes with
$\\{200\cdot 2^{i}\\}_{i=0}^{6}$ elements and a fixed time step $\Delta
t=0.001$. The simulation is performed until the final time $t=0.15$, past the
theoretical waiting time $t^{*}=0.125$. We denote by $x_{R}^{0}$ the node that
corresponds to the location of the right interface at time $t=0$ and track the
density value at this node for each time step.

The results of our experiment are presented in Figure 2. Panels (A) and (B)
show the density profiles at times $t=0$, $t=0.125$, and $t=0.15$, while (C)
and (D) plot the density at $x_{R}^{0}$ as a function of time for different
values of $N$, as obtained by the log-density and the mixed method,
respectively. In both (A) and (B), we observe that the interface remains
stationary at $t=0.125$ while the internal profile continues to evolve. At
time $t=0.15$, a clear movement of the interface is noticed. The panels (C)
and (D) reveal for both methods, when $t$ is equal to the theoretical waiting
time (indicated by a black vertical line), the density value at $x_{R}^{0}$
gradually decreases towards 0 as the mesh is refined by increasing $N$.
Furthermore, for each value of $N$, the mixed method gives a smaller density,
and thus a smaller error, at the node $x^{R}_{0}$ when $t=t^{*}$. The smaller
error can be ascribed to the piecewise constant spatial discretization of the
mixed method, which allows for a sharper capture of the interface, reducing
the error near it.

###  4.3. Higher dimensional Barenblatt solution

|  | Log-density method | Mixed method
---|---|---|---
m | N | dt |  | Error in
---
$[-3,3]$
Order |  | Error in
---
$[-6,6]$
Order | dt |  | Error
---
in $[-3,3]$
Order |  | Error in
---
$[-6,6]$
Order
2 | 32x32 | 1/5 | 2.35e-02 | - | 6.37e-02 | - | 1/10 | 1.29e-01 | - | 2.45e-01 | -
| 64x64 | 1/20 | 6.61e-03 | 1.828 | 2.97e-02 | 1.101 | 1/20 | 6.46e-02 | 1.004 | 1.25e-01 | 0.971
| 128x128 | 1/80 | 1.70e-03 | 1.956 | 1.33e-02 | 1.161 | 1/40 | 3.21e-02 | 1.007 | 6.38e-02 | 0.971
| 256x256 | 1/320 | 4.29e-04 | 1.989 | 5.77e-03 | 1.202 | 1/80 | 1.61e-02 | 0.999 | 3.25e-02 | 0.973
3 | 32x32 | 1/5 | 9.55e-03 | - | 1.98e-01 | - | 1/10 | 7.47e-02 | - | 3.77e-01 | -
| 64x64 | 1/20 | 2.99e-03 | 1.677 | 1.21e-01 | 0.706 | 1/20 | 3.72e-02 | 1.006 | 2.01e-01 | 0.905
| 128x128 | 1/80 | 7.81e-04 | 1.936 | 6.98e-02 | 0.795 | 1/40 | 1.85e-02 | 1.011 | 1.13e-01 | 0.835
| 256x256 | 1/320 | 1.99e-04 | 1.970 | 4.01e-02 | 0.801 | 1/80 | 9.24e-03 | 0.999 | 6.22e-02 | 0.859
4 | 32x32 | 1/5 | 5.97e-03 | - | 3.40e-01 | - | 1/10 | 4.39e-02 | - | 5.28e-01 | -
| 64x64 | 1/20 | 1.55e-03 | 1.949 | 2.33e-01 | 0.547 | 1/20 | 2.19e-02 | 1.002 | 2.86e-01 | 0.884
| 128x128 | 1/80 | 4.04e-04 | 1.936 | 1.46e-01 | 0.672 | 1/40 | 1.09e-02 | 1.009 | 1.72e-01 | 0.729
| 256x256 | 1/320 | 1.03e-04 | 1.977 | 8.92e-02 | 0.712 | 1/80 | 5.45e-03 | 0.998 | 1.02e-01 | 0.764

Table 2. Convergence results for 2D Barenblatt initial data.

(a) Density; triangular elements

(b) Error; triangular elements

(c) Density; quadrilateral elements

(d) Error; quadrilateral elements

Log-density method

(e) Density; triangular elements

(f) Error; triangular elements

(g) Density; quadrilateral elements

(h) Error; quadrilateral elements

Mixed Method

Figure 3. Numerical results at final time $t=1$ for 2D Barenblatt initial data
with $m=3$, $\Delta t=0.025$, and $N\approx 1024$ elements.

(a) Density; tetrahedral elements.

(b) Error; tetrahedral elements

(c) Density; hexahedral elements

(d) Density; hexahedral elements

Figure 4. Numerical results at final time $t=0.2$ using the log-density method
for 3D Barenblatt initial data with $m=3$, $\Delta t=0.025$, and $N\approx
4100$ elements.

The d-dimensional version of (42) is given by

$\rho_{B}(\bm{x},t)=(t+1)^{-k}\left(s_{0}-\frac{k(m-1)}{2dm}\frac{|\bm{x}|^{2}}{(t+1)^{2k/d}}\right)_{+}^{\frac{1}{m-1}},\quad
t>0.$ (46)

Analogous to the one dimensional case, we verify the accuracy of the two
schemes in two dimensions by using (46) as initial data with $x\in[-6,6]^{2}$,
$s_{0}=1$, and $m\in\\{2,3,4\\}$ and conduct a spacetime mesh-refinement study
for the $L^{2}$ error at final time $t=0.2$. The $L^{2}$-error is recorded in
the entire domain and in the box $[-3,3]^{2}$ away from the interface. We use
a sequence of spatial meshes with $\\{(32\times 32)\cdot 2^{i}\\}_{i=0}^{3}$
tensor-product elements and set the timestep size to $\frac{1}{5\cdot 4^{i}}$
for the log-density method and to $\frac{1}{10\cdot 2^{i}}$ for the mixed
method, correspondingly. The results of the study are summarized in Table 2.
Akin to the one dimensional case, we observe second-order spatial accuracy and
first-order temporal accuracy of the log-density method in the box
$[-3,3]^{2}$ where the solution remains smooth at final time. In the same
region, the mixed method is observed to be first-order accurate in space and
time as expected. Additionally, the order of convergence is observed to decay
with $m$ when the error is measured in the entire domain. This deterioration
is slower for the case of the mixed method.

In figure 3, we plot the density and error profiles obtained using the two
schemes at final time $t=0.2$ for the data (46) with $m=3$, $\Delta t=0.025$,
and $N\approx 1024$ elements. Results on both triangular and quadrilateral
meshes are shown. Panels (A)-(D) display the profiles computed using the log-
density method, while (E)-(F) display the profiles computed using the mixed
method.

We additionally evolve the three-dimensional version of 46 on tetrahedral and
hexahedral meshes of the domain $\Omega_{h}=[-6,6]^{3}$ using the log-density
method until final time $t=0.2$ with $\Delta t=0.025$ and $N\approx 4100$
elements. We set $m=3$ and $s_{0}=1$. The plots of this simulation are shown
in figure 4. In the region $[-3,3]^{3}$, the $L^{2}$-error is 0.0196 for
hexahedral elements and 0.0674 for tetrahedral elements. We remark that the
mixed method can also be used for this case, but only hexahedral or regular
tetrahedral meshes may be used due to the lack of a proper mass-lumping
quadrature formula for general tetrahedral elements.

###  4.4. Merging Gaussians

We further investigate the robustness of the two schemes on a popular test
case consisting of two initial Gaussian peaks that merge into a single peak
under the action of the PME [22, 27, 10]. The initial condition for this test
is given by

$\rho^{0}(x,y)=e^{-20\left((x-0.3)^{2}+(y-0.3)^{2}\right)}+e^{-20\left((x+0.3)^{2}+(y+0.3)^{2}\right)},$
(47)

where the domain is taken to be $\Omega=[-1,1]^{2}$. We set $m=3$, $\Delta
t=0.001$, and evolve the initial data (47) using the two schemes on a
triangular grid with 3750 elements. We record the solution at $t=0$, $t=0.15$,
and $t=0.3$, and plot the results of the simulation in figure 5, where one can
observe that in both cases the two peaks move towards each other and
eventually start merging.

(a) $t=0$

(b) $t=0.15$

(c) $t=0.3$

Log-density method

(d) $t=0$

(e) $t=0.15$

(f) $t=0.3$

Mixed method

Figure 5. Merging Gaussians test with $m=3,\Delta t=0.001$, and $N=3750$
triangular elements.

###  4.5. Complex support

(a) $t=0$

(b) $t=0.5$

(c) $t=1.0$

Log-density method

(d) $t=0$

(e) $t=0.5$

(f) $t=1.0$

Mixed method

Figure 6. Complex support test with $m=3,\Delta t=0.001$, and $N=3750$
triangular elements.

In the final numerical example, we consider the following initial data [22,
28, 2]

$\rho^{0}(x,y)=\left\\{\begin{aligned}
25\left(0.25^{2}-\left(\sqrt{x^{2}+y^{2}}-0.75\right)^{2}\right)^{\frac{3}{2(m-1)}},&\quad\sqrt{x^{2}+y^{2}}\in[0.5,1]\text{
and }(x<0\text{ or }y<0),\\\
25\left(0.25^{2}-x^{2}-(y-0.75)^{2}\right)^{\frac{3}{2(m-1)}},&\quad
x^{2}+(y-0.75)^{2}\leq 0.25^{2}\text{ and }x\geq 0,\\\
25\left(0.25^{2}-(x-0.75)^{2}-y^{2}\right)^{\frac{3}{2(m-1)}},&\quad(x-0.75)^{2}+y^{2}\leq
0.25^{2}\text{ and }y\geq 0,\\\ 0,&\quad\text{
otherwise},\end{aligned}\right.$ (48)

on the domain $[-2,2]$. The support of intial data (48) has the shape of a
horseshoe or a partial donut. We set $m=3$, $\Delta t=0.001$, and evolve the
initial data using the two schemes on a triangular grid consisting of 3750
elements. The solution is recorded at times $t=0$, $t=0.5$, and $t=1.0$, and
the computed profiles are presented in figure 6. The plots demonstrate that
both the schemes deliver reliable results, as the horseshoe ends are observed
to evolve towards each other before ultimately intersecting. Additionally, in
line with previous observations, the boundary of the support is captured more
sharply by the mixed method.

We note that unlike Lagarangian schemes, the two methods exhibit a robust
handling of the topology change without requiring any interpolation of the
solution on to a new mesh when this event occurs.

##  5\. Conclusion

We presented two distinct, first-order spacetime accurate, finite element
approaches to the PME. The log-density approach is constructed for a problem
with Neumann boundary conditions, and the properties of mass conservation,
energy stability, unique solvability, and bound preservation on unstructured
Delaunay meshes are proved. The scheme is shown to be second order in space
and first order in time. The mixed approach is constructed for a problem with
Dirichlet boundary conditions and is shown to be mass conservative, energy
stable, and positivity preserving (under a CFL condition). The mixed scheme is
shown to be first order in both space and time. Both schemes can handle
compactly supported initial data without the need of any perturbation.

Funding This research is partially supported by NSF grant DMS-2012031.

Data Availability Enquiries about data availability shall be directed to the
authors.

## Declaration

Conflict of interest The authors declare that there are no known conflicts of
interest assocaited with this work.

## References

  * [1] D. G. Aronson, L. A. Caffarelli, and S. Kamin, How an initially stationary interface begins to move in porous medium flow, SIAM Journal on Mathematical Analysis, 14 (1983), pp. 639–658.
  * [2] M. Baines, M. Hubbard, and P. Jimack, A moving mesh finite element algorithm for the adaptive solution of time-dependent partial differential equations with moving boundaries, Applied Numerical Mathematics, 54 (2005), pp. 450–469.  Selected papers from the 16th Chemnitz Finite Element Symposium 2003.
  * [3] Baranger, Jacques, Maitre, Jean-François, Oudin, and Fabienne, Connection between finite volume and mixed finite element methods, ESAIM: Mathematical Modelling and Numerical Analysis - Modélisation Mathématique et Analyse Numérique, 30 (1996), pp. 445–465.
  * [4] I. S. Barenblatt, On some unsteady motions of a liquid and gas in a porous medium, Akad. Nauk SSSR. Prikl. Mat. Mekh., 16 (1952), pp. 67–78.
  * [5] T. J. Barth, Aspects of unstructured grids and finite-volume solvers for the euler and navier-stokes equations, AGARD, special course on unstructured grid methods for advection dominated flows, (1992).
  * [6] T. J. Barth, Numerical methods for gas-dynamics systems on unstructured meshes, Lecture Notes in Computational Science and Engineering VOL 5, Springer, 1999.
  * [7] E. D. Benedetto and D. Hoff, An interface tracking algorithm for the porous medium equation, Transactions of the American Mathematical Society, 284 (1984), pp. 463–500.
  * [8] M. BERTSCH and R. D. PASSO, A numerical treatment of a superdegenerate equation with applications to the porous media equation, Quarterly of Applied Mathematics, 48 (1990), pp. 133–152.
  * [9] J. Boussinesq, Recherches théoriques sur l’écoulement des nappes d’eau infiltrés dans le sol et sur le débit de sources. comptes rendus acad, Sci. / J. Math Pures Appl., 10 (1903/04), pp. 5–78.
  * [10] J. A. Carrillo, B. Düring, D. Matthes, and D. S. McCormick, A Lagrangian scheme for the solution of nonlinear diffusion equations using moving simplex meshes, J. Sci. Comput., 75 (2018), pp. 1463–1499.
  * [11] Z. Chen, H. Huang, and J. Yan, Third order maximum-principle-satisfying direct discontinuous Galerkin methods for time dependent convection diffusion equations on unstructured triangular meshes, J. Comput. Phys., 308 (2016), pp. 198–217.
  * [12] C. Duan, C. Liu, C. Wang, and X. Yue, Numerical methods for porous medium equation by an energetic variational approach, J. Comput. Phys., 385 (2019), pp. 13–32.
  * [13] J. C. M. Duque, R. M. P. Almeida, and S. N. Antontsev, Convergence of the finite element method for the porous media equation with variable exponent, SIAM Journal on Numerical Analysis, 51 (2013), pp. 3483–3504.
  * [14] G. Fu and Z. Xu, High-order space-time finite element methods for the Poisson-Nernst-Planck equations: positivity and unconditional energy stability, arXiv: [math.NA].
  * [15] J. Gratton and C. Vigo, Evolution of self-similarity, and other properties of waiting-time solutions of the porous medium equation: the case of viscous gravity currents, European J. Appl. Math., 9 (1998), pp. 327–350.
  * [16] J. L. Graveleau and P. Jamet, A finite difference approach to some degenerate nonlinear parabolic equations, SIAM J. Appl. Math., 20 (1971), pp. 199–223.
  * [17] Y. Gu and J. Shen, Bound preserving and energy dissipative schemes for porous medium equation, Journal of Computational Physics, 410 (2020), p. 109378.
  * [18] M. E. GURTIN, R. C. MacCAMY, E. A. SOCOLOVSKY, and E. A. Socolousky, A coordinate transformation for the porous media equation that renders the free-boundary stationary, Quarterly of Applied Mathematics, 42 (1984), pp. 345–357.
  * [19] A. Harten, On the symmetric form of systems of conservation laws with entropy, J. Comput. Phys., 49 (1983), pp. 151–164.
  * [20] T. J. R. Hughes, L. P. Franca, and M. Balestra, A new finite element formulation for computational fluid dynamics. V. Circumventing the Babuška-Brezzi condition: a stable Petrov-Galerkin formulation of the Stokes problem accommodating equal-order interpolations, Comput. Methods Appl. Mech. Engrg., 59 (1986), pp. 85–99.
  * [21] S. Jin, L. Pareschi, and G. Toscani, Diffusive relaxation schemes for multiscale discrete-velocity kinetic equations, SIAM J. Numer. Anal., 35 (1998), pp. 2405–2439.
  * [22] C. Liu and Y. Wang, On Lagrangian schemes for porous medium type generalized diffusion equations: a discrete energetic variational approach, J. Comput. Phys., 417 (2020), pp. 109566, 27.
  * [23] H. Liu and W. Maimaitiyiming, Efficient, positive, and energy stable schemes for multi-D Poisson-Nernst-Planck systems, arXiv:2001.08350 [math.NA].
  * [24] Y. Liu, C.-W. Shu, and M. Zhang, High order finite difference weno schemes for nonlinear degenerate parabolic equations, SIAM Journal on Scientific Computing, 33 (2011), pp. 939–965.
  * [25] M. S. Metti, J. Xu, and C. Liu, Energetically stable discretizations for charge transport and electrokinetic models, J. Comput. Phys., 306 (2016), pp. 1–18.
  * [26] Monsaingeon, Léonard, An explicit finite-difference scheme for one-dimensional generalized porous medium equations: Interface tracking and the hole filling problem, ESAIM: M2AN, 50 (2016), pp. 1011–1033.
  * [27] C. Ngo and W. Huang, A study on moving mesh finite element solution of the porous medium equation, J. Comput. Phys., 331 (2017), pp. 357–380.
  * [28] C. Ngo and W. Huang, Adaptive finite element solution of the porous medium equation in pressure formulation, Numerical Methods for Partial Differential Equations, 35 (2019), pp. 1224–1242.
  * [29] F. Otto, The geometry of dissipative evolution equations: the porous medium equation, Comm. Partial Differential Equations, 26 (2001), pp. 101–174.
  * [30] R. E. Pattle, Diffusion from an instantaneous point source with a concentration-dependent coefficient, The Quarterly Journal of Mechanics and Applied Mathematics, 12 (1959), pp. 407–409.
  * [31] I. Pop and W. Yong, A numerical approach to degenerate parabolic equations, Numer. Math., 92, pp. 357–381.
  * [32] P.-A. Raviart and J. M. Thomas, A mixed finite element method for 2nd order elliptic problems, in Mathematical aspects of finite element methods (Proc. Conf., Consiglio Naz. delle Ricerche (C.N.R.), Rome, 1975), Lecture Notes in Math., Vol. 606, Springer, Berlin, 1977, pp. 292–315.
  * [33] J. Schöberl, C++11 Implementation of Finite Elements in NGSolve, 2014.  ASC Report 30/2014, Institute for Analysis and Scientific Computing, Vienna University of Technology.
  * [34] E. Tadmor, Skew-selfadjoint form for systems of conservation laws, J. Math. Anal. Appl., 103 (1984), pp. 428–442.
  * [35] J. L. Vázquez, The Porous Medium Equation, Oxford University Press, Oxford, 2007.
  * [36] M. Westdickenberg and J. Wilkening, Variational particle schemes for the porous medium equation and for the system of isentropic Euler equations, M2AN Math. Model. Numer. Anal., 44 (2010), pp. 133–166.
  * [37] J. Xu and L. Zikatanov, A monotone finite element scheme for convection-diffusion equations, Mathematics of Computation, 68 (1999), pp. 1429–1446.
  * [38] Q. Zhang and Z.-L. Wu, Numerical simulation for porous medium equation by local discontinuous Galerkin finite element method, J. Sci. Comput., 38 (2009), pp. 127–148.