On Lagrangian schemes for porous medium type generalized diffusion equations: a discrete energetic variational approach

TL;DR

This paper introduces a systematic Lagrangian scheme for porous medium equations using a discrete energetic variational approach, ensuring properties like energy dissipation and mass conservation, and demonstrates its effectiveness through numerical experiments.

Contribution

It develops a novel framework for deriving energy-dissipative Lagrangian schemes for generalized diffusion equations, applicable to multidimensional porous medium equations.

Findings

The scheme accurately captures free boundaries in PME.

It effectively estimates waiting times in simulations.

Numerical results confirm the scheme's stability and energy dissipation properties.

Abstract

In this paper, we present a systematic framework to derive a Lagrangian scheme for porous medium type generalized diffusion equations by employing a discrete energetic variational approach. Such discrete energetic variational approaches are analogous to energetic variational approaches in a semidiscrete level, which provide a basis of deriving the "semi-discrete equations" and can be applied to a large class of partial differential equations with energy-dissipation laws and kinematic relations. The numerical schemes derived by this framework can inherit various properties from the continuous energy-dissipation law, such as conservation of mass and the dissipation of the discrete energy. As an illustration, we develop two numerical schemes for the multidimensional porous medium equations (PME), based on two different energy-dissipation laws. We focus on the numerical scheme based on the energy-dissipation law with $\frac{1}{2} \int_{\Omega} |\mathbf{u}|^2 \mathrm{d} \mathbf{x}$ as the dissipation. Several numerical experiments demonstrate the accuracy of this scheme as well as its ability in capturing the free boundary and estimating the waiting time for the PME in both 1D and 2D.

Tables

Table 1. Table 4.1 : The convergence rate of numerical solutions for ρ 0 ​ ( X ) = B α ​ ( X , 1 ) subscript 𝜌 0 𝑋 subscript 𝐵 𝛼 𝑋 1 \rho_{0}(X)=B_{\alpha}(X,1) at the finite time T = 1 for α = 2 𝛼 2 \alpha=2 and α = 4 𝛼 4 \alpha=4 .

$\alpha =3$		scheme 1				scheme 2
N	$\tau$	Error at $X\left(0\right)$	Order	$L_m^2$-error	Order	Error at $X\left(0\right)$	Order	$L_m^2$-error	Order
$51$	1/100	2.6421e-04		0.0036		9.0421e-05		3.5109e-04
$101$	1/400	7.0619e-05	1.9036	0.0016	1.1699	2.2692e-05	1.9945	9.5494e-05	1.8784
$201$	1/1600	1.9303e-05	1.8712	7.0440e-04	1.1836	5.7058e-06	1.9917	2.5714e-05	1.8929
$\alpha =4$		scheme 1				scheme 2
N	$\tau$	Error at $X\left(0\right)$	Order	$L_m^2$-error	Order	Error at $X\left(0\right)$	Order	$L_m^2$-error	Order
$51$	1/100	2.6925e-04		0.0051		2.3969e-04		5.0786e-04
$101$	1/400	7.6531e-05	1.8148	0.0027	0.9175	6.0120e-05	1.9953	1.4894e-04	1.7697
$201$	1/1600	2.2855e-05	1.7435	0.0014	0.9475	1.5073e-05	1.9959	4.3149e-05	1.7873

Table 2. Table 4.2 : Numerical Error of right interface at T = 1 𝑇 1 T=1 for initial data ρ 0 ​ ( X ) = B α ​ ( X , 1 ) subscript 𝜌 0 𝑋 subscript 𝐵 𝛼 𝑋 1 \rho_{0}(X)=B_{\alpha}(X,1) with different α 𝛼 \alpha ( α = 4 , 5 , 6 𝛼 4 5 6 \alpha=4,5,6 and 8 8 8 ) at T = 1 𝑇 1 T=1 . The numerical solutions are computed by scheme 2.

$\alpha$	$N=51,\tau =1/100$	$N=101,\tau =1/400$	Order
4	3.7241e-04	8.9717e-05	2.0534
5	4.4021e-04	1.0629e-04	2.0502
6	4.6867e-04	1.1329e-04	2.0486
8	4.8008e-04	1.1619e-04	2.0468

Table 3. Table 4.3 : The convergence rate of numerical solutions with ρ 0 ​ ( 𝐗 ) = B α ​ ( 𝐗 , 1 ) subscript 𝜌 0 𝐗 subscript 𝐵 𝛼 𝐗 1 \rho_{0}(\mathbf{X})=B_{\alpha}(\mathbf{X},1) at T = 0.1 𝑇 0.1 T=0.1 for the PME with α = 2 𝛼 2 \alpha=2 and α = 4 𝛼 4 \alpha=4 .

$\alpha =2$				$\alpha =4$
N	$\tau$	${ℒ}^2$-error	Order	N	$\tau$	${ℒ}^2$-error	Order
132	1/100	6.5388e-04		135	1/100	0.0066
524	1/400	1.6053e-04	2.0262	516	1/400	5.5299e-04	1.7807
2103	1/1600	3.9867e-05	2.0096	2124	1/1600	1.6518e-04	1.9982

Equations

$${\rho }_t=\nabla \cdot (D(\rho ,\mathbf{x})\nabla \rho ),\mathbf{x}\in \Omega \subset {\mathbb{R}}^d,\alpha >1,t>0,$$

$$\rho (\mathbf{x},0)={\rho }_0(\mathbf{x}),\mathbf{x}\in \Omega ,$$

$$\frac{\partial \rho }{\partial \mathbf{n}}=0,\mathbf{x}\in \partial \Omega ,t>0,$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{E}^{\text{total}}(t)=-2\mathcal{D}(t),$$

$$\delta \mathcal{A}={\int }_0^T{\int }_{{\Omega }_t}(f_{\text{inertial force}}-f_{\text{conservative force}})\cdot \delta \mathbf{x}\mathrm{d}\mathbf{x}\mathrm{d}t,$$

$$\delta \mathcal{D}={\int }_{{\Omega }_t}{f}_{\text{dissipative force}}\cdot \delta {\mathbf{x}}_t\mathrm{d}\mathbf{x}.$$

$$\frac{\delta A}{\delta \mathbf{x}}=\frac{\delta \mathcal{D}}{\delta {\mathbf{x}}_t}.$$

$$\mathbf{x}(\mathbf{X},0)=\mathbf{X},$$

$$\mathbf{u}(\mathbf{x}(\mathbf{X},t),t)={\mathbf{x}}_t(\mathbf{X},t),$$

$$F(\mathbf{X},t)={\nabla }_{\mathbf{X}}\mathbf{x}(\mathbf{X},t),$$

$${\int }_{{\Omega }_t}\rho (\mathbf{x}(t),t)\mathrm{d}\mathbf{x}={\int }_{{\Omega }_0}{\rho }_0(\mathbf{X})\mathrm{d}\mathbf{X},$$

$$\rho (\mathbf{x}(\mathbf{X},t),t)={\rho }_0(\mathbf{X})/\det F(\mathbf{X},t).$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{{\Omega }_t}\rho (\mathbf{x}(\mathbf{X},t),t)\mathrm{d}\mathbf{x}$$

$$={\int }_{{\Omega }_t}\left({\rho }_t+\nabla \rho \cdot \mathbf{u}+\rho (\nabla \cdot \mathbf{u})\right)\mathrm{d}\mathbf{x}=0.$$

$${\rho }_t+\nabla \cdot (\rho \mathbf{u})=0.$$

$$\rho (\mathbf{x},t)={\rho }_0\tilde{\circ }{\mathbf{X}}^{-1}(\mathbf{x},t)=\frac{{\rho }_0\circ {\mathbf{X}}^{-1}(\mathbf{x},t)}{\det F(\mathbf{X},t)},$$

$$\varphi (\mathbf{x},t)={\varphi }_0\circ {\mathbf{X}}^{-1}(\mathbf{x},t).$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{\Omega }\omega (\rho )\mathrm{d}\mathbf{x}=-{\int }_{\Omega }\eta (\rho )|\mathbf{u}{|}^2\mathrm{d}\mathbf{x},$$

$$\delta \mathcal{A}=-\delta {\int }_0^T{\int }_{{\Omega }_0}\omega ({\rho }_0(X)/\det F)\det F\mathrm{d}\mathbf{X}$$

$$=-{\int }_0^T{\int }_{{\Omega }_0}\left(-{\omega }_{\rho }\left(\frac{{\rho }_0(X)}{\det F}\right)\frac{{\rho }_0(X)}{\det F}+\omega \left(\frac{{\rho }_0(X)}{\det F}\right)\right)(F^{-\mathrm{T}}:{\nabla }_{\mathbf{X}}\delta \mathbf{x})\det F\mathrm{d}\mathbf{X},$$

$$\delta \mathcal{A}$$

$$f_{\text{dissipative force}}=\frac{\delta \mathcal{D}}{\delta {\mathbf{x}}_t}=\eta (\rho ){\mathbf{x}}_t.$$

$$\eta (\rho )\mathbf{u}=-\nabla p,$$

$${\rho }_t=\nabla \cdot \left(\frac{\rho }{\eta (\rho )}\nabla p(\rho )\right),$$

$${\rho }_t=\Delta {\rho }^{\alpha },\alpha >1,$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{\Omega }\frac{1}{\alpha -1}{\rho }^{\alpha }=-{\int }_{\Omega }\rho |\mathbf{u}{|}^2\mathrm{d}x,$$

$$\rho \mathbf{u}=-\nabla {\rho }^{\alpha },$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{\Omega }\omega (\rho )\mathrm{d}\mathbf{x}=-{\int }_{\Omega }|\mathbf{u}{|}^2\mathrm{d}\mathbf{x},\omega (\rho )=\{\begin{array}{ll}& 2\rho \ln \rho ,\alpha =2,\\ & \frac{\alpha }{(\alpha -1)(\alpha -2)}{\rho }^{\alpha -1},\alpha >2,\end{array}$$

$$\mathbf{u}=-\nabla \left(\frac{\alpha }{\alpha -1}{\rho }^{\alpha -1}\right),$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{\Omega }\frac{1}{2}|\nabla \phi {|}^2+\frac{1}{4{ϵ}^2}({\phi }^2-1{)}^2=-{\int }_{\Omega }{\varphi }^2|\mathbf{u}{|}^2\mathrm{d}\mathbf{x},$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{E}_h\left(\mathbf{\Xi }(t)\right)=-2{\mathcal{D}}_h(\mathbf{\Xi }(t),{\mathbf{\Xi }}^{\prime }(t)),$$

$${\mathcal{A}}_h(\mathbf{\Xi }(t))=-{\int }_0^T{E}_h(\mathbf{\Xi }(t))\mathrm{d}t$$

Full text

1 Introduction

Many mathematical models in physics, chemistry and biology can be viewed as generalized diffusions with different choices of free energy and dissipation functional from an energetic variational viewpoint [22]. Examples include the porous medium equation (PME) [57], Keller-Segel model [30], Poisson-Nernst-Planck (PNP) [19, 63] and Cahn-Hilliard equations [8, 36].

One simple example of a generalized diffusion equation is given by

[TABLE]

where $\rho$ is a non-negative function, $D(\rho,\bm{x})$ is the diffusion coefficient, $\Omega$ is a bounded domain and $\mathbf{n}$ is the external normal direction in $\partial\Omega$. When $D(\rho,\bm{x})=c\rho^{\alpha-1}$ ($\alpha>1$), (1.1) is the porous medium equation (PME)[57].

There has been an increasing interesting in designing Lagrangian methods for diffusion equations [9, 60, 10, 28, 12, 39, 11]. The Lagrangian methods is particularly suitable for problems involving sharp interface and free boundary, especially for those with singularity. However, it is a common challenge to solve Lagrangian schemes numerically, especially in high spatial dimensions. Moreover, most of previous Lagrangian methods start with the nonlinear PDEs for the Lagrangian maps [10, 12, 18], it might be difficult to choose the proper weak form to preserve the original variational structure.

In this paper, we present a systematic framework to construct Lagrangian schemes to generalized diffusion equations (1.1) by a discrete energetic variational approach, which can be easily applied to multiple spatial dimensions. A discrete energy variational approach is an analogue to energetic variational approaches [35, 22] in a semidiscrete level, which provides a general framework to derive the “semi-discrete equations”, a system of ordinary differential equations in time, from a discrete energy-dissipation law. Our approach is quite close to the methods based on gradient flow structure in the $L_{2}$-Wasserstein metric [9, 60, 10, 28, 12, 39, 11], but various energetic variational forms can be utilized in such a framework. As an illustration, we develop two numerical schemes for the multidimensional porous medium equation (PME), based on two different energy-dissipation laws.

The PME has a wide application in physical and biological problems, such as the flow of an ideal gas though a porous medium [32], radiative transfer theory [33], biological aggregation [56], population dynamics [62], and tumor growth [50]. Theoretical studies [44, 29, 31, 61, 3, 57] have shown that the existence of the solutions to the PME, which will have a compact support at any time $t>0$ if the initial data has a compact support. The interface between the compact support and zero-region, called the free boundary, will moves outward in a finite speed, known as the finite speed propagation. Unlike the heat equation, the solution of the PME could become non-smooth even if the initial data is smooth. Moreover, for a certain initial data, the interface will not move until a finite positive time, called the waiting time [31, 3].

From a computational perspective, its lack of regularity and free boundary [1, 3, 57] poses challenges for numerical simulations for the PME. For instance, numerical solutions by standard numerical approaches, such as PCSFE (Predictor-Correction Algorithm and Standard Finite Element) method, may contain oscillations near the free boundary, which cannot be removed by raising the degree of finite element space and/or by refining spatial meshes [65]. On the other hand, it is difficult to track the movement of free boundary and estimate the waiting time of the PME in high accuracy. During the past, quite a number of numerical methods have been proposed for the PME [55, 23, 17, 24, 51, 52, 6, 25, 40, 26, 13, 65, 37, 7, 5, 4, 42, 43], most of them are Eulerian methods and focus on the one-dimensional case.

A commonly used numerical approach is the interface tracking algorithm [55, 23, 17, 24, 6, 25, 40], in which the interface equation are solved in Lagrangian coordinate, while the solution inside the numerical support is updated in Eulerian coordinate. However, it is not easy to apply such method to the PME in high spatial dimensions and the case with complex support. In order to eliminate the oscillations around the free boundary, some numerical methods from hyperbolic conservation law, such as relaxation scheme [26, 13], locally discontinuous Galerkin (LDG) [65] and WENO [37], have been applied to the PME successfully. Since the solution of the PME always steeper near the interface, some adaptive moving mesh methods have proved to be useful in solving the PME [7, 5, 4, 42, 43], especially in multiple space dimensions. A typical moving mesh method updates the computational mesh according to the monitor function defined by the numerical solution at current time, which might effect the dynamic of evolution. Hence, it might be difficult to track the free boundary and estimate the waiting time in high accuracy within the adaptive moving mesh methods.

The rest of the paper is organized as follows. In section 2, we briefly review the energetic variational approaches for porous medium type generalized diffusions in a continuous level. In section 3, we introduce a discrete energetic variational approach and apply it to construct two numerical schemes for the PME based on two different energy-dissipation laws. The numerical experiments for the PME in one and two spatial dimensions are presented in sections 4. We focus on the numerical scheme based on the energy-dissipation law with $\frac{1}{2}\int_{\Omega}|\mathbf{u}|^{2}\mathrm{d}\bm{x}$ as the dissipation. These numerical results demonstrate the accuracy of our Lagrangian scheme as well as its ability in capturing the free boundary and estimating the waiting time for the PME in both 1D and 2D.

2 Energetic Variational Approaches

For a given energy-dissipation law and the kinematic (transport) relations, an energetic variational approach provides a general framework to determine the dynamics of system in a unique and well-defined way, through two distinct variational processes: Least Action Principle (LAP) and Maximum Dissipation Principle (MDP) [35, 22]. This approach is originated from pioneering work of Onsager [45, 46] and Rayleigh [53], and has been successfully applied to build up many mathematical models for complex fluids [35, 54, 19, 22].

The starting point of an energetic variational approach is the first and second laws of thermodynamics [22], which yields the following energy-dissipation law

[TABLE]

for an isothermal closed system. $E^{\text{total}}$ is the total energy, which is the sum of the Helmholtz free energy $\mathcal{F}$ and the kinetic energy $\mathcal{K}$. $2\mathcal{D}$ is the rate of energy dissipation, which is related to the entropy production in thermodynamics. The Least Action Principle states that the dynamics of a conservative system is determined as a critical point of the action functional $\mathcal{A}(\bm{x})=\int_{0}^{T}\mathcal{K}-\mathcal{F}\mathrm{d}t$ with respect to $\bm{x}$ (the trajectory, if applicable) [2, 22], i.e.,

[TABLE]

where $\Omega_{t}$ is the physical domain at time $t$. On the other hand, for a dissipative system ($\mathcal{D}\geq 0$), according to Onsager [45, 46], the dissipative force can be obtained by minimization of the dissipation functional $\mathcal{D}$ with respect to the “rate” $\bm{x}_{t}$, known as Maximum Dissipation Principle (MDP), i.e.,

[TABLE]

Then, according to the force balance (Newton’s second law, in which the inertial force plays role of $ma$), we have

[TABLE]

We refer the reader to [22] for more detailed description of energetic variational approaches. Here we only focus on the modeling of generalized diffusions by an energetic variational approach.

Generalized diffusions are concerned with a conserved quantity $\rho(\bm{x},t)$. An energetic variational approach to a diffusion based on a Lagrangian description to the system, in which the system is characterized by a flow map $\bm{x}(\mathbf{X},t):\Omega_{0}\rightarrow\Omega_{t}$ satisfies

[TABLE]

where $\mathbf{X}$ are the Lagrangian coordinates and $\bm{x}$ are Eulerian coordinates. For a given flow map $\bm{x}(\mathbf{X},t)$, the velocity $\mathbf{u}(\bm{x}(\mathbf{X},t),t)$ of the system is defined by

[TABLE]

and the deformation matrix (the deformation gradient) $F(\mathbf{X},t)$ is the Jacobian matrix of the flow map $\bm{x}(\mathbf{X},t)$:

[TABLE]

which determines the kinematic relations of physical quantities in Lagrangian coordinates.

For a conserved quantity $\rho(\bm{x},t)$, let $\rho_{0}(\mathbf{X})$ be the initial mass density, then the conservation of mass

[TABLE]

indicates

[TABLE]

In Eulerian coordinates, the conservation of mass (2.7) is equivalent to

[TABLE]

Hence, $\rho(\bm{x}(\mathbf{X},t),t)$ satisfies a transport equation

[TABLE]

One can view (2.8) as a kind of composition of the flow map $\bm{x}(\mathbf{X},t)$ and the initial density $\rho_{0}(\mathbf{X})$, i.e.,

[TABLE]

where $\mathbf{X}^{-1}:\Omega_{t}\rightarrow\Omega_{0}$ is the inverse of the flow map $\bm{x}(\mathbf{X},t)$.

Remark 2.1**.**

The scalar transport equation $\phi_{t}+(\mathbf{u}\cdot\nabla)\phi=0$ is equivalent to $\phi(\bm{x}(\mathbf{X},t),t)=\phi_{0}(\mathbf{X})$ in the Lagrangian coordinates, which also can be viewed as a composition

[TABLE]

The scalar transport equation is a suitable kinematic for designing a Lagrangian methods for the gradient flow system, such as Allen-Cahn equation.

Remark 2.2**.**

The “initial data”, or “reference data”, $\rho_{0}(\mathbf{X})$ ($\phi_{0}(\mathbf{X})$) carries many information of the solutions, which may not be available for general problems. In practice, one may obtain $\rho_{0}(\mathbf{X})$ ($\phi_{0}(\mathbf{X})$) from other methods, such as those of Eulerian approaches.

In a porous medium type generalized diffusion considered in this paper, the overall energy-dissipation law of the system is given by

[TABLE]

where the kinetic energy $\mathcal{K}$ is neglected, $\omega(\rho)$ is the free energy density, and the energy dissipation $\mathcal{D}$ is taken to be $\frac{1}{2}\int_{\Omega}\eta(\rho)|\mathbf{u}|^{2}$ like the Darcy’s law (the friction to the resting media) [22].

We first perform the LAP, i.e., compute the variation of $\mathcal{A}=\int_{0}^{T}-\mathcal{F}\mathrm{d}t$ with respect to $\bm{x}(\mathbf{X},t)$. Direct computation results in

[TABLE]

where $\delta\bm{x}$ is the test function. It can be noticed that even for this simple case, the variational result in Lagrangian coordinates is quite complicated, which involves $F^{-1}$ and $\det F$.

Pull (2.13) back to Eulerian coordinates and apply the integration by parts, one can get

[TABLE]

where the boundary term vanishes due to the choice of $\delta\bm{x}$.

The MDP can be done by taking variation of $\mathcal{D}$ with respect to $\bm{x}_{t}$, one can easily obtain that

[TABLE]

By the force balance ($\frac{\delta\mathcal{A}}{\delta\bm{x}}=\frac{\delta\mathcal{D}}{\delta\bm{x}_{t}}$), an energetic variational approach leads to

[TABLE]

where $p=\omega_{\rho}\rho-\omega$. The force balance equation (2.16) defines the velocity in the transport equation (2.9). Combining (2.9) and (2.16), we get a generalized diffusion equation

[TABLE]

which satisfies the energy-dissipation law (2.12).

For the PME

[TABLE]

a commonly used energy-dissipation law is

[TABLE]

which is consistence with the Wasserstein gradient flows [27, 47]. Since $\omega(\rho)=\frac{1}{\alpha-1}\rho^{\alpha}$ and $\eta(\rho)=\rho$, the force balance gives

[TABLE]

which in turns yields the original PME (2.18) by (2.17).

It should be remarked that same governing equations can be obtained by using different $\omega(\rho)$ and $\eta(\rho)$ [see [18] for examples in the PME]. Correspondingly, different numerical schemes can be derived based on different energy-dissipation laws. Besides the classical energy-dissipation law (2.19) and the two used in ([18]), we can employ another energy-dissipation law

[TABLE]

for the PME with $\alpha\geq 2$ to construct a numerical scheme. Following the above computation, one can verify that, under the energy-dissipation law (2.21), the force balance results in

[TABLE]

which is equivalent to (2.20) on the compact support of $\rho_{0}(\mathbf{X})$. We will show that in the following sections the numerical scheme derived from of (2.21) has an advantage in tracking the free boundary in the PME, especially in multiple spatial dimensional situation.

Remark 2.3**.**

Although the energy-dissipation law (2.21) is only defined for the PME with $\alpha\geq 2$, the force balance (2.22), i.e., the trajectory equation, is well defined for $\alpha>1$. In the later section, we will show that the numerical scheme derived from of (2.21) also works well for the case that $1<\alpha<2$, although it can not be interpreted through (2.21) in such cases.

Remark 2.4**.**

In this paper, we only consider the diffusions satisfy the energy-dissipation law (2.12). The free energy in (2.12) can be generalized to more complicated form. For instance, the famous Cahn-Hilliard equation can be viewed as a $H^{1}$-diffusion satisfies the energy-dissipation law [22, 36]

[TABLE]

where $\varphi$ is the conserved quantity satisfies $\varphi_{t}+\nabla\cdot(\varphi\mathbf{u})=0$.

3 A discrete energetic variational approach and numerical schemes

In this section, we first introduce the abstract framework of a discrete energetic variational approach, and apply it to construct some Lagrangian schemes for the porous-medium type generalized diffusion equations. As an example, we construct two numerical schemes for the PME based on energy-dissipation laws (2.19) and (2.21).

As mentioned in the beginning, a discrete energetic variational approach is an analogue to an energetic variational approach in a semidiscrete level. For a given continuous energy-dissipation law, we can write down a discrete energy-dissipation, that is

[TABLE]

by introducing a spatial discretization, where $\bm{\Xi}(t)=\left(\Xi_{1}(t),\Xi_{2}(t),\ldots,\Xi_{K}(t)\right)\in\mathbb{R}^{K}$ is the state variable for the discrete energy-dissipation law (3.1). In the following, we only discuss the case that the kinetic energy $\mathcal{K}=0$ in the continuous energy-dissipation law.

Similar to a continuous energetic variational approach stated in sect. 2, we can get the governing equation of $\bm{\Xi}(t)$ by performing LAP and MDP. In a semidiscrete level, LAP corresponds to taking variation of the discrete action functional

[TABLE]

with respect to $\bm{\Xi}(t)$, and MDP corresponds to taking variation of the discrete dissipation functional $\mathcal{D}_{h}(\bm{\Xi},\bm{\Xi}(t))$ with respect to $\bm{\Xi}^{\prime}(t)$. Hence, the force balance results in:

[TABLE]

which is a system of $K$ nonlinear ordinary differential equations, i.e., “semi-discrete equations”, of $\bm{\Xi}(t)$. Due to the assumption that the dissipation is quadratic in “rate” $\bm{\Xi}^{\prime}(t)$, equations (3.2) can be written as

[TABLE]

where ${\sf D}\left(\bm{\Xi}(t)\right)$ is a $K\times K$ matrix. A numerical scheme can be obtained by introducing a proper temporal discretization to (3.3). If ${\sf D}\left(\bm{\Xi}(t)\right)$ is an identical matrix (which is not true in general), (3.3) can be viewed as a fast descent on each component $\Xi_{k}$.

A discrete energetic variational approach follows the “discretize-then-variation” strategy [20, 15]. The idea of “discretize-then-variation”, or “discretize-then-minimize”, has been successfully applied to study numerous problems [20, 15, 16, 58, 59, 9, 10, 64]. In a recent published book [20], Furihata and Matsuo show that “discretize-then-variation” is a systematic way to derive a structure-preserving numerical methods for a large class of PDE. By a discrete energetic variational approach, we are able to apply this idea to more general problems with energy-dissipation laws and kinematic relations. In general, “discrete-then-variation” and “variation-then-discrete” may give us different numerical schemes. It is important to notice that the force balance (2.4) uses the strong form of the variational results [22], since the test functions may be in different space. Hence, the numerical scheme derived by a discrete energetic variational approach normally gives a better approximation to the original energy-dissipation law.

3.1 Lagrangian schemes for generalized diffusions

Now we are ready to construct some Lagrangian schemes for porous medium type generalized diffusions with energy-dissipation law (2.12) by a discrete energetic variational approach. The procedure can be extended to more complicated type of diffusions, such as Cahn-Hilliard and PNP equations.

For generalized diffusions considered in this paper, we can discretize the energy-dissipation law by approximating the flow map $\bm{x}(\mathbf{X},t)$ directly, and the value of $\rho(\bm{x},t)$ is determined by the kinematic relation (transport equation) (2.8), i.e., a composition of the flow map $\bm{x}(\mathbf{X},t)$ and the initial data $\rho_{0}(\mathbf{X})$:

[TABLE]

This is the main difference between our numerical approach and most of traditional Eulerian approaches to diffusion equations.

Remark 3.1**.**

Although the original free energy possesses certain convexity property with respect to $\rho$, it is not clear whether the convexity is preserved in Lagrangian coordinates, i.e., with respect to the flow map $\bm{x}(\mathbf{X},t)$ in high dimensional situation. Hence, for the general case, the evolution of the flow map $\bm{x}(\mathbf{X},t)$ may approach to a “local” minimum, which is not necessarily to be the global minimum in theory.

In current study, we will discretize $\bm{x}(\mathbf{X},t)$ by a piecewise linear map. An advantage of a piecewise linear approximation to the flow map $\bm{x}(\mathbf{X},t)$ is that the deformation matrix $F(\mathbf{X},t)$ is piecewise constant (matrix) for give $t$, so are $\det F(\mathbf{X},t)$ and $F(\mathbf{X},t)^{-1}$.

One way to construct a piecewise linear approximation to the flow map $\bm{x}(\mathbf{X},t)$ is to use a finite element method [11]. To be more precisely, let $\mathcal{T}_{h}$ be triangulation of $\Omega_{0}\in\mathbb{R}^{d}$. $\mathcal{T}_{h}$ consists of a set of simplexes $\{\tau_{e}~{}|~{}e=1,\ldots M\}$ and a set of nodal points $\mathcal{N}_{h}=\{\mathbf{X}_{1},\mathbf{X}_{2},\ldots,\mathbf{X}_{N}\}$. Define the finite element space by

[TABLE]

which is a linear finite element space. Then the flow map $\bm{x}(\mathbf{X},t)$ can be approximated by

[TABLE]

where $\phi_{i}(\mathbf{X}):\mathbb{R}^{d}\rightarrow\mathbb{R}$ is the hat function satisfies $\phi_{i}(\mathbf{X}_{j})=\delta_{ij}$, and $\bm{\xi}_{i}(t)\in\mathbb{R}^{d}$ are coefficients to be determined. Hence,

[TABLE]

where $K=N\times d$, which is the state variable in (3.1).

Since $\bm{x}_{h}(\mathbf{X}_{i},t)=\bm{\xi}_{i}$, one can view $\bm{\xi}_{i}(t)$ as coordinates in $\Omega_{t}$ and $\bm{x}_{h}(\mathbf{X}_{i},t)$ as a trajectory of a “particle” $\mathbf{X}_{i}$. The finite element method enables us to compute the deformation matrix $F_{h}(\mathbf{X})=\nabla_{X}\bm{x}_{h}$ explicitly for given $\bm{\Xi}$. We can write down the deformation matrix $F_{h}(\mathbf{X})$ as a $d\times d$-matrix-valued function of $\bm{\Xi}$ on each element $\tau_{e}$, denoted by

[TABLE]

The admissible set of $\bm{\Xi}$ is

[TABLE]

Correspondingly, the admissible set for $\bm{x}_{h}$ is defined by

[TABLE]

The non-negativity of $\rho(\bm{x}_{h}(\mathbf{X},t),t)$ is naturally preserved as if $\bm{x}_{h}(\mathbf{X},t)$ is in the admissible set $F_{ad}^{\bm{x}_{h}}$. It can be noticed that $F_{ad}^{\bm{\Xi}}$ is not a convex set when $d\geq 2$, which imposes difficulty in numerical analysis.

In the following, we restrict ourselves in a two dimensional case, the numerical schemes in other spatial dimensions follow easily from this. In order to simplify the notation, we let $\bm{\xi}_{i}(t)=(a_{i}(t),b_{i}(t))$, and denote

[TABLE]

Hence, $\bm{\Xi(t)}=(\bm{a}(t),\bm{b}(t))$. Let $N(i)$ be all the indices $e$ such that $\mathbf{X}_{i}$ is contained in $\tau_{e}$ for given $\mathbf{X}_{i}\in\mathcal{N}_{h}$. The support of $\phi_{i}$ is denote by $G(i)=\cup_{e\in N(i)}\tau_{e}$. For each element $\tau_{e}$, the nodes of $\tau_{e}$ is denoted by $\mathbf{X}^{e}_{1}$, $\mathbf{X}^{e}_{2}$, $\mathbf{X}^{e}_{3}$, the global index of the nodes of $\tau_{e}$ are denoted by $en(e,l)$ ($l=1,2,3$).

Substituting (3.6) into (2.12), we get a discrete action functional

[TABLE]

and the discrete dissipation

[TABLE]

where

[TABLE]

is the deformation matrix $F(\mathbf{X},t)$ on each element $\tau_{e}$ at $t$, which can be written down as a function of $a_{en(e,l)}(t)$ and $b_{en(e,l)}(t)$ ($l=1,2,3$) explicitly [see Appendix A].

By taking variation of (3.8) with respect to $a_{i}(t)$ and $b_{i}(t)$, we have

[TABLE]

where

[TABLE]

on each element $\tau_{e}$, $\chi=a_{en(e,l)}$ or $b_{en(e,l)}$.

On the meanwhile, taking variations of (3.14) with respect to $a_{i}^{\prime}(t)$ and $b_{i}^{\prime}(t)$ results in

[TABLE]

where Einstein summation convention is used, $\eta(\rho)=\eta\left(\rho_{0}(\mathbf{X})/\det F_{e}(t)\right)$, and $M_{ij}(\bm{a}(t),\bm{b}(t))$ is defined by

[TABLE]

Hence, the force balance equation (3.2) results in the “semi-discrete equations”

[TABLE]

for $i=1,2,\ldots,N$, which can be discretized in time by using some numerical method for systems of ordinary differential equations.

Remark 3.2**.**

In this special case, ${\sf D}(\bm{\Xi}(t),\bm{\Xi}^{\prime}(t))$ in (3.3) is given by

[TABLE]

which is a symmetric matrix. From a computational point of view, ${\sf D}\left(\bm{\Xi}(t),\bm{\Xi}^{\prime}(t)\right)$ plays a important role in maintaining the integrity of the flow map in the evolution process.

Although both (3.10) and (3.12) involve the numerical integration over each element, they can be computed out by centroid method (known as midpoint method in 1D). As $F(\mathbf{X},t)$, $\det F(\mathbf{X},t)$ and $F(\mathbf{X},t)^{-1}$ are approximated in a piecewise constant manner, using high-accuracy numerical quadrature over each element cannot improve the numerical accuracy. This is another advantage of the piecewise linear approximation to the flow map, which is actually quadrature-free and can be applied to a high spatial dimensional case. The computational cost is roughly proportional to the number of nodes (“particle”). One can view our numerical approach as a type of cell-centered Lagrangian scheme, where the momentum is defined at the nodes and the other variables (density, pressure, and specific internal energy) are cell-centered [38]. In the following, we denote

[TABLE]

where $\mathbf{X}_{c}^{e}$ is the centroid of $\tau_{e}$.

In order to get a numerical scheme, we need to introduce a proper temporal discretization to the “semi-discrete equation” (3.14). One can use explicit Euler scheme, and the numerical scheme can be written as

[TABLE]

where

[TABLE]

Although the explicit Euler scheme is simple in the numerical implementation, one have to choose $\tau$ to be significantly small to ensure $\bm{\Xi}^{n+1}_{h}\in F_{ad}^{\bm{\Xi}}$ and the dissipation of the discrete energy.

A better approach is to adopt a backward Euler scheme for the temporal discretization, i.e.,

[TABLE]

where $M_{ij}^{*}(\bm{a}^{n},\bm{b}^{n})$ is defined by

[TABLE]

and $\gamma_{i}=n$ or $n+1$. Here we have the choice of taking $\eta(\rho_{e}(t))$ and $\det F_{e}(t)$ in (3.13) explicitly or implicitly, such that ${\sf M}^{*}$ depends $\bm{a}^{n},\bm{b}^{n}$, but is independent with $\bm{a}^{n+1},\bm{b}^{n+1}$. The theoretical analysis of the choice are in the progress.

Remark 3.3**.**

In general, we cannot have an explicit form for the “semi-discrete equation” (3.14) and the numerical scheme (3.18) in high dimensional situations, as (3.10) and (3.13) depend on the triangulation $\mathcal{T}_{h}$. In practice, $M_{ij}(\bm{a}(t),\bm{b}(t))$ and $\frac{\delta\mathcal{A}_{h}}{\delta a_{i}}(\bm{a}(t),\bm{b}(t))$ can be assembled using the standard technique in the finite element methods, that is, summing the results on each element over the mesh [34]. [See Appendix B for the procedure of numerical implementation.]

The numerical scheme (3.18) can be written as

[TABLE]

where

[TABLE]

is the discrete energy, ${\sf D}^{*}_{n}={\sf D}^{*}(\bm{\Xi}^{n})$ is defined by (3.15) with ${\sf M}={\sf M}^{*}$. Although (3.20) is highly nonlinear equation in

[TABLE]

as an advantage of our discrete energetic variational approach, we can get a solution of (3.20) by solving minimization problem in $F_{ad}^{\bm{\Xi}}$:

[TABLE]

where

[TABLE]

Moreover, we can prove the following theorem for the numerical scheme (3.20):

Theorem 3.1**.**

If ${\sf D}^{*}_{n}$ is a symmetric positive-definite matrix, and $\omega$ satisfies

[TABLE]

then for given $\bm{\Xi}^{n}\in F_{ad}^{\bm{\Xi}}$, there exists a solution $\bm{\Xi}^{n+1}$ to numerical scheme (3.20) such that the following discrete energy dissipation law holds, i.e.,

[TABLE]

Proof.

Since a minimizer of the minimization problem (3.22) is the solution of (3.20), we only need to show that there exists a minimizer of $J(\bm{\Xi})$ in

[TABLE]

Since for $\forall\bm{\Xi}\in\partial F_{ad}^{\bm{\Xi}}$, $\bm{\Xi}=\infty$, following the proof in the Lemma 3.1 in [11], the existence of a minimizer can be obtained by showing the set

[TABLE]

is a non-empty compact subset of $\mathbb{R}^{K}$. Obviously, $\bm{\Xi}^{n}\in\mathcal{S}$, so $\mathcal{S}$ is non-empty.

In order to show $\mathcal{S}$ is compact, we first show that $S$ is bounded. Since ${\sf D}^{*}_{n}$ is positive-definite, there exists $\lambda_{1}>0$ such that $\forall\bm{\Xi}\in\mathcal{S}$

[TABLE]

Hence, $\mathcal{S}$ is bounded. Next we show $\mathcal{S}$ is closed in $\mathbb{R}^{K}$, it suffices to show that the limit $\widetilde{\bm{\Xi}}$ for any sequence $(\bm{\Xi}^{(k)})_{k=1}^{\infty}\subset\mathcal{S}$ is in $F_{ad}^{\bm{\Xi}}$. For $\forall e\in\{1,2,\ldots M\}$ and all $k$

[TABLE]

where $|\tau_{e}|$ is the area of element $\tau_{e}$. Since $\omega\left(\frac{\rho_{e}^{0}}{\det F_{e}(\bm{\Xi}^{(k)})}\right)\det F_{e}(\bm{\Xi}^{(k)})\rightarrow\infty$ if $\det F_{e}\rightarrow 0$, we can conclude that $\det F_{e}(\bm{\Xi}^{(k)})>0$ is uniformly bounded away from zero. So $\det F_{e}(\bar{\bm{\Xi}})>0$ for all $e$, which means $\widetilde{\bm{\Xi}}\in F_{ad}^{\bm{\Xi}}$.

If $\bm{\Xi}^{n+1}$ is a minimizer in $\mathcal{S}$, we have

[TABLE]

which completes the proof. ∎

Remark 3.4**.**

In one-dimensional case, due to $\det F=F$, it is easy to show that $E_{h}(\bm{\Xi})$ is convex. Hence, we can have a stronger result similar to that in [18]:

[TABLE]

where the first inequality follows the convexity of $E_{h}(\bm{\Xi})$. Moreover, since $\bar{F}_{ad}^{\bm{\Xi}}$ is the closed convex set in 1D [18], we can have a uniquely solvable result of numerical scheme (3.20).

*In high dimensional situation ($d\geq 2$), we can not get a uniquely solvable result of numerical scheme (3.20) due to the lack of convexity of $J(\bm{\Xi})$ and $E(\bm{\Xi})$ [11]. *

Remark 3.5**.**

Although there exists a (local) minimizer in $\mathcal{S}\subset F_{ad}^{\bm{\Xi}}$ for $J(\bm{\Xi})$, which is a solution of numerical scheme (3.18) decreases the discrete energy $E_{h}$, we still need to choose a proper optimization methods to find a minimizer $\bm{\Xi}^{n+1}\in\mathcal{S}$. For the PME in 1D or 2D, numerical tests show that a standard damped Newton’s method with fixed step-size is adequate to this purpose. For more general problem, we can adopt optimizaition methods wih line search to find such a minimizer. Since we may only find a local minimizer of $J(\bm{\Xi})$ in $\mathcal{S}$ due to lack of convexity for the general problem, the dynamical evolution of a flow-map based Lagrangian methods might be different from Eulerian methods.

Remark 3.6**.**

The above approach start with the spatial discretization to the flow map (3.6), we can also begin with introducing a temporal discretization to the continuous energy-dissipation law (2.12) by

[TABLE]

where $\bm{x}^{n}=\bm{x}(\mathbf{X},t_{n})$, and $\gamma_{i}=n$ or $n+1$, as in (3.19). The RHS of (3.31) can be viewed as an approximation to $\frac{\mathrm{d}}{\mathrm{d}t}E(\bm{x}(\mathbf{X},t^{*}))$ for some $t^{*}\in[t^{n},t^{n+1}]$, we obtain (3.31) by further approximating

[TABLE]

Using the spatial discretization (3.6), we will have the following discrete energy-dissipation law:

[TABLE]

where ${\sf D}^{*}_{n}={\sf D}^{*}(\bm{\Xi}_{n})$ is same to that in (3.20). Note that

[TABLE]

for some $\bm{\Xi}^{c}=(1-c)\bm{\Xi}^{n+1}+c\bm{\Xi}^{n}$. It is straightforward to show that if $\bm{\Xi}^{n+1}\in F_{ad}^{\bm{\Xi}}$ satisfies

[TABLE]

then $\bm{\Xi}^{n+1}$ satisfies (3.33).

In our scheme (3.20), we approximate $\nabla_{\bm{\Xi}}E_{h}(\bm{\Xi}_{c})$ by $\nabla_{\bm{\Xi}}E_{h}(\bm{\Xi}^{n+1})$, which causes the difference between (3.25) and (3.33).

Remark 3.7**.**

We can design a numerical scheme that satisfies (3.33) exactly. Let

[TABLE]

where $\bm{e}_{i}$ is the standard orthonormal basis in $\mathbb{R}^{K}$. Note the RHS of (3.33) can be written as

[TABLE]

On the other hand, the LHS of (3.33) can be written as

[TABLE]

where

[TABLE]

Hence, direct computation shows that $\bm{\Xi}^{n+1}=\bm{\Xi}^{n}_{(K)}$ satisfies (3.33) if

[TABLE]

which give us a numerical scheme

[TABLE]

*The scheme (3.41) preserves the discrete energy-dissipation law (3.33), and might be useful when the variation of $E(\bm{\Xi})$ ($\mathcal{A}(\bm{\Xi})$) cannot be computed efficiently. However, to our knowledge, the numerical analysis and experiments for such type of scheme is lacking. We will explore this type of scheme in the future work. *

Remark 3.8**.**

One can also adopt some high-order temporal discretization to the “semi-discrete equations” (3.14). Since the resulting numerical scheme might be complicated to deal with, we will study the high-order temporal discretization in the future work.

Remark 3.9**.**

In general, discrete-then-variation and variation-then-discrete may give us different numerical scheme. In order to get the numerical schemes (3.16) and (3.18), one should substitute (3.6) into a particular weak form of the force balance (2.16) (strong form of the variational results), and introduce a proper approximation and temporal discretization. A weak form of (2.16) can be written as

[TABLE]

where $\bm{y}$ is a test function, $F=\nabla_{\mathbf{X}}\bm{x}$. One can get (3.18) by taking the test function $\bm{y}=\phi_{i}$ ($i=1,\ldots N$) and approximating (3.42) by

[TABLE]

It should be remarked that we might need to approximate $\det F$ in both side of (3.42) in different manners (explicitly or implicitly) to get back to (3.18). From an energetic variational approach point of view, the test functions may be in different space for the LAP and the MDP. Hence, starting with the force balance (2.16) may not give us a structure-preserving Lagrangian scheme without using a proper weak form.

At the end of this subsection, we briefly talk about the post-process after we obtain the discrete flow map $\bm{x}_{h}(\mathbf{X},\bm{x})$. According to the kinematic relation, $\rho(\bm{x})$ can be computed by

[TABLE]

Hence, we can compute the density on each element, i.e, the density in the centroid of each element is by

[TABLE]

where $\bm{x}_{c}^{e}$ is the centroid of $\bm{x}_{h}(\tau_{e})$, while $\mathbf{X}_{c}^{e}$ is the centroid of $\tau_{e}$. For each node $X_{i}$, since the determinant of the deformation matrix $F$ may be different in different elements contain $\mathbf{X}_{i}$, we can compute the $\rho_{h}(\bm{x}_{h}(\mathbf{X}_{i}),t)$ in each nodes by

[TABLE]

3.2 Application in the PME

The above framework works for any porous medium type generalized diffusions has the energy-dissipation law (2.12) and the kinematic relation (2.9). Next, we apply such framework to develop two numerical schemes for the PME based on energy-dissipation law (2.19) and (2.21). The RHS of (3.18) can be computed from (3.10) and (3.11) for a given $\omega(\rho)$. For the LHS in (3.18), we take

[TABLE]

for the PME with energy-dissipation law (2.19), while for the energy-dissipation law (2.21), we take

[TABLE]

We call the numerical scheme based on energy-dissipation law (2.19) as the scheme 1 [(3.18) with (3.47)] and the numerical scheme based on energy-dissipation law (2.21) as scheme 2 [(3.18) with (3.48)]. One can develop other numerical schemes for the PME by other energy-dissipation laws, such as two used in Ref. [18].

Due to the degeneracy of the PME at $\rho=0$, the semi-discrete equation (3.2), corresponding to the energy-dissipation law (2.19), is also degenerate at the region that $\rho_{0}(\mathbf{X})=0$, that is to say, if $\mathbf{X}_{k}\in\mathbb{R}^{2}\backslash\overline{\Omega}_{c}$, where $\Omega_{c}$ is the support of $\rho_{0}(\mathbf{X})$, then

[TABLE]

which means the nodes (“particle”) in such region have no well-defined velocity. In the meantime, we can only derive the PME from the energy-dissipation law (2.21) on the compact support of $\rho_{0}(\mathbf{X})$. Hence, in the following, we always assume that $\Omega_{0}$ is the compact support of $\rho_{0}(\mathbf{X})$ or $\rho_{0}(\mathbf{X})>0$ on $\Omega_{0}$, which also guarantees the conditions of Theorem 3.1 hold.

It is a commonly used strategy that solving the PME only within the solution support [55, 17, 24, 6, 7, 5, 43, 18], known as the non-embedding approach. The main challenge of this approach is that the evolution of free boundary has to be tracked explicitly [43]. Our cell-centered Lagrangian schemes enable us to treat the movement of free boundary in a uniform way, as the velocity of the free boundary is also well-defined, which is a a major advantage in high dimensional situations. In later section, we will show that our schemes, especially the scheme based on energy-dissipation law (2.21) can capture the free boundary of the PME, in both 1D and 2D, without explicitly track the movement of the free boundary.

Remark 3.10**.**

Since the corresponding trajectory equation (2.22) of (2.21) is exactly the equation for the movement of free boundary, we expect the numerical scheme based on the energy-dissipation law (2.21) has an advantage in tracking the movement of the free boundary. Although the energy-dissipation law (2.21) is only valid for the PME with $\alpha\geq 2$, our numerical tests show that the scheme 2 also works well for $1<\alpha<2$.

Numerical schemes in 1D: We can write down our two numerical schemes, based on energy-dissipation law (2.19) and (2.21), explicitly in one-dimensional case. Let $\Omega_{0}=[\xi_{l},\xi_{r}]$ be the compact support of $\rho_{0}(X)$, and $\xi_{l}=X_{1}<X_{2}<\ldots<X_{N}=\xi_{r}$ be nodes in $\Omega_{0}$. We can approximate the flow map $x(X,t)$ by

[TABLE]

where $\phi_{i}(X)$ is a hat function satisfies $\phi_{i}(X_{j})=\delta_{ij}$. Let $\bm{a}(t)=(a_{1}(t),a_{2}(t),\ldots,a_{N}(t))^{\prime}$ and $h_{i}=X_{i+1}-X_{i}$, in one-dimensional case, the admissible set of $\bm{a}$ is simply as

[TABLE]

The discrete action functional and the discrete dissipation can be written as

[TABLE]

where

[TABLE]

For the energy-dissipation law (2.19), direct computation results in

[TABLE]

and

[TABLE]

Hence, in the one-dimensional case, our scheme 1 based on energy-dissipation law (2.19) can be written as

[TABLE]

where ${\sf M}$ is a triangular matrix, given by

[TABLE]

where we define $X_{1/2}=0$, $X_{N+1/2}=0$ and $h_{0}=h_{N}=h$ to simplify the notation.

Remark 3.11**.**

In Ref. [18], the authors develop several numerical schemes for the one-dimensional PME based on different energy-dissipation laws. For the energy-dissipation law (2.19), their numerical scheme can be written as (scheme 0 in [18])

[TABLE]

Formally, with the equidistant node, our scheme (3.55) only differs from theirs in the temporal discretization for all inner points. A drawback of their temporal discretization is that, in order to prevent the solution to escape from the admissible set, they should use a specific and inefficient damped Newton methods, in which the step-size is roughly proportional to $\min(\rho_{0}(X_{i}))$. Our temporal discretization follow the maximum dissipation principle, in which $M(\bm{a})$ control the movement of each nodes. A standard damped Newton with a fixed step-size can prevent the solution to escape from the admissible set in our scheme. Their scheme can only be applied to the region in which $\rho_{0}>0$ (a different scheme is used for the movement of free boundary). Moreover, the starting point of the numerical approach in Ref. [18] is the force balance (2.20) in strong form, it is difficult to extend their approach to a high dimensional situation. Strictly speaking, their numerical schemes may not preserve the original energy-dissipation law in general.

Similarly, in the one-dimensional case, our scheme 2, which is based on energy-dissipation law (2.21), can be written as

[TABLE]

where

[TABLE]

where we define $X_{1/2}=0$, $X_{N+1/2}=0$ and $h_{0}=h_{N}=h$ to simplify the notation.

4 Numerical Experiments

In this section, we present some numerical results to the PME in both 1D and 2D to demonstrate the accuracy of our numerical methods. We’ll focus on our scheme 2, which is based on the energy-dissipation law (2.21). The error in space of a numerical solution is measured in $L^{2}$-norms defined by

[TABLE]

where $e_{h}(\bm{x})$ is the difference between the numerical solution $\rho_{h}(\bm{x})$ and the exact solution $\rho(\bm{x})$. We compute the numerical integration in (4.1) by the centroid method (known as midpoint method in 1D), which defines the discrete $\mathcal{L}^{2}$-norm, i.e.,

[TABLE]

where $\bm{x}_{c}^{e}$ is the centroid of $\bm{x}_{h}(\tau_{e})$, $\rho_{h}(\bm{x}_{c}^{e})$ is computed by (3.45).

4.1 One-dimensional problems

4.1.1 Barenblatt-Pattle solution

To verify the accuracy of our numerical methods, we first consider a benchmark solution, the Barenblatt-Pattle solution, which is a exact weak solution for the PME established by Barenblatt [21] and Pattle [48]. The one-dimensional Barenblatt-Pattle solution is given by,

[TABLE]

where $k=(\alpha+1)^{-1}$ and $u_{+}=\mathrm{max}(u,0)$. For any time $t>0$, this solution has a compact support $[-\xi_{\alpha}(t),\xi_{\alpha}(t)]$ with the interface $|x|=\xi_{\alpha}$ moving outward at a finite speed, where

[TABLE]

We take the Barenblatt solution at $t=1$ , i.e., $B_{\alpha}(x,1)$, as the initial data, and compare our numerical solution at time $T$ with $B_{\alpha}(x,T+1)$. Fig. 4.1 shows the numerical and exact solutions for $\alpha=4$ at $T=1$ and $T=10$, where the numerical solutions are computed by scheme 1 (3.55) [shown by blue-square] and scheme 2 (3.58) [shown by red-circle] with $\Omega_{0}$ to be the compact support of the initial data. The results demonstrate that both our numerical schemes can approximate the exact solutions well without oscillation. The numerical solutions by scheme 2 (3.58) approximate the exact solutions better near the interface.

The converge rate of Barenblatt solutions with $\alpha=3$ and $4$ for both numerical schemes is shown in Table 4.1. The error at $X=0$ and the error in $\mathcal{L}^{2}$-norm is presented. For both schemes, the converge rate of the error at $X=0$ is second order in space since the solution is smooth far away from the interface. The scheme 1 has a first-order converge rate, while the scheme 2 can achieve second-order in $\mathcal{L}^{2}$-norm at $T=1$. When $\alpha$ becomes larger, the converge rate of scheme 2 in $\mathcal{L}^{2}$-norm can keep in second order. As expected, the numerical error of scheme 2 is smaller than that of scheme 1, as it track the movement of free boundary better. We can reduce the numerical error of our scheme 1 by tracking the movement of the free boundary explicitly, i.e., replacing the equations of $X_{b}\in\partial\Omega_{0}$ with the equation of free boundary. In the following, we’ll focus on our scheme 2, which based on energy-dissipation law (2.21). All the following numerical solutions are computed by scheme 2.

Fig. 4.2(a) shows the trajectory of each node for Barenblatt solution with $\alpha=4$ [$N=51$, $\tau=0.01$, scheme 2]. It can be noticed that the final grid is almost uniform. This is because the solution doesn’t become steeper during the the time evolution. We plot in Fig. (4.2)(b) the evolution of the numerical interface for the Barenblatt solution [$N=51$, $\tau=0.01$, scheme 2], with four different parameters $\alpha=4,5,6$ and $8$, from $T=0$ to $1$, in which the solid line is the position of the exact interface, and the circle indicates the position of the numerical interface. The results show that the exact interface can be approximated in high accuracy.

Quantitatively, we show the error of the right interface of our scheme 2 with different $\alpha$ ($\alpha=4,5,6$ and $8$) at $T=1$ in Table. 4.2, which shows that our scheme 2 can track the movement of the free boundary in second-order, even for large $\alpha$.

Remark 4.1**.**

Although the scheme 2 is derived by the energy-dissipation law (2.21), which is valid for $\alpha\geq 2$, numerical tests show that this scheme also works well for the case with $1<\alpha<2$. Fig. 4.3 shows the numerical and exact solutions for $\alpha=5/3$ and $\alpha=1.1$ at $T=1$ with $\rho_{0}(X)=B_{\alpha}(x,1)$.

4.1.2 Waiting Time

Now, we study the waiting time phenomenon of the PME by our numerical methods. The waiting time phenomenon occurs for a certain type of initial data. For instance, if $\rho_{0}(X)$ satisfies

[TABLE]

then there exist a positive waiting time for $0\leq\theta\leq 1$ [3]. According the theoretical result [3], for $0\leq\theta\leq 1/4$, the waiting time for the initial date (4.5) is $\frac{1}{2(\alpha+1)(1-\theta)}$.

Most of previous studies estimate the waiting time from the trajectory of numerical interface [55, 6, 43], but the numerical waiting time may not be clearly estimated in some cases. As an alternative approach, Nakaki and Tomoeda estimate the waiting time for the one-dimensional PME by transforming the original equation into another equation whose solution will blow up at the waiting time of the original PME [41]. Recently, Duan et al. proposed an elegant criterion to determine the waiting time in one-dimensional case, and they manually set the velocity of interface to be zero before the numerical waiting time. We can easily adopt the criterion proposed in [18] into our numerical scheme in one-dimensional case. Fig. 4.4 shows the numerical solutions for the initial data (4.5) with $\theta=0$ and $\alpha=4$ at various time, in which we manually refine the first and last element to improve the numerical accuracy ($N=107,\tau=10^{-4}$). The numerical waiting time we obtained is $t=0.1054$, which is consistent with theoretical results ($t=0.1$).

A natural question is whether we can estimate the waiting time from our numerical interface without manually setting the velocity of interface to be zero, as it is difficult to apply similar criterion to the PME in high spatial dimensions. Fig. 4.5 (a) shows the right numerical interface computed by our numerical method with and without manually setting the velocity of interface to be zero before the numerical waiting time for the initial data (4.5) with $\theta=0$ and $\alpha=4$. Surprisingly, these two line almost coincide with each other after the numerical waiting time. Moreover, we noticed that the compact support will shrink at beginning due to the numerical approximation. These motivate us to define the numerical waiting time by the time when the numerical support begins to expand, i.e.,

[TABLE]

where $\xi_{b}(t)$ is the position of left and right interface at $t$ and $\xi_{b}^{0}$ is the initial position of left and right interface. Although no theoretical justification is available at this point, numerical experiments show that this criterion give a clear estimation to the waiting time for all cases that we tested. Fig. 4.5 (b) shows another example for the initial data (4.5) with $\theta=\frac{1}{2}$ and $\alpha=7$, in which the numerical interface for both approaches are shown. Both approaches indicate that the waiting time in this case is about $t^{*}\approx 0.124$.

4.2 2D Simulations

We now apply our numerical scheme 2, which is based on the energy-dissipation law (2.21), to the PME in two dimensions. Although the explicit form of the numerical scheme cannot be given in a general mesh, using the explicit form of $F_{e}$ as a function of $a_{en(e,l)}$ and $b_{en(e,l)}$ on each element $\tau_{e}$ (shown in the Appendix), we can compute $\frac{\delta\mathcal{A}_{h}}{\delta a_{i}}(\bm{a}^{n+1},\bm{b}^{n+1})$ ($\frac{\delta\mathcal{A}_{h}}{\delta b_{i}}(\bm{a}^{n+1},\bm{b}^{n+1})$) and ${\sf M}^{*}(\bm{a}^{n},\bm{b}^{n})$ by Algorithm 1 and Algorithm 2 in remark 3.6, during the computer implementation.

4.2.1 Barenblatt-Pattle solution

We first validate our numerical scheme by studying the 2D Barenblatt-Pattle solution. The Barenblatt-Pattle solution in $d$-dimensions is given by

[TABLE]

where $k=(\alpha-1+2/d)^{-1}$, $C_{0}$ is a constant, related to the initial mass. This solution is radially symmetric, self-similar, and has compact support $|\bm{x}|\leq\xi_{\alpha}(t)$ for any finite time, where

[TABLE]

We take Barenblatt solution (4.7) for $C_{0}=0.1$ at $t=1$ as the initial data, and compare the numerical solution at $T$ with the exact solution $B_{\alpha}(\bm{x},T+1)$. Since the Barenblatt solution is steeper near the interface, we use the non-uniform mesh in $\Omega_{0}$ in our 2D simulation, which can largely reduce computational cost. The initial non-uniform mesh in $\Omega_{0}$ is generated by DistMesh [49].

Fig. 4.6 shows the numerical solutions for $\alpha=4$ at $T=1$ in a non-uniform mesh with $N=516$. The numerical error at $T=1$ are plot in Fig. 4.6(b), which indicates that the $L^{\infty}$-norm of the numerical solution is $O(10^{-4})$ for $N=516$. The numerical and exact interfaces at $T=1$ are shown in Fig. 4.6(c), which demonstrated that our numerical scheme can track the movement of free boundary in 2D in high accuracy with a small $N$.

We now test the converge rate of 2D Barenblatt solution. The error for numerical solutions with $\rho_{0}(\mathbf{X})=B_{\alpha}(\mathbf{X},1)$ at $T=0.1$ in $\mathcal{L}^{2}$-norm for both $\alpha=2$ and $\alpha=4$ is shown in Table 4.3. It shows that we can achieve a second-order convergence rate in space for both $\alpha=2$ and $\alpha=4$ in $\mathcal{L}^{2}$-norm.

4.2.2 Waiting time

Now, we apply our numerical scheme (scheme 2) to the PME with the initial data has a waiting-time phenomenon in two-dimensional situations. We take the initial value as:

[TABLE]

which has a positive waiting time.

The numerical solutions for this initial data at various time are shown in Fig. 4.9 [$N=2105$ and $\tau=10^{-3}$]. We use the non-uniform mesh in $\Omega_{0}$, which is dense around the free boundary. In order to validate our numerical results, we also compute the same initial date within the axisymmetric assumption, which can reduce the problem into a one-dimensional problem and enable us to apply the criterion in [18] to estimate the waiting time. The numerical solutions obtained within the axisymmetric assumption for $N=201$ and $\tau=10^{-4}$ are shown in Fig. 4.9 (d)-(f).

The location of free boundary by the 2D simulation [blue solid line] and 1D simulation with the axisymmetric assumption [red dashed line] is plotted in Fig. 4.8(a), where we define the location of free boundary in 2D simulation results by

[TABLE]

where $\xi_{0}$ is the radius of the initial support. Both results indicate the wait time is around $0.128$. Fig. 4.8(b) plot the the numerical interface at $t=0.5$ for both approaches, in which the free boundary in the 2D simulation result is shown by blue-square and the result obtained by 1D simulation is shown in red solid-line . It can be concluded that the numerical solutions obtained by both approaches are consistent, which validates our numerical scheme 2 in 2D, and we can have a clear estimation to the waiting time in 2D from the numerical interface obtained by the scheme 2.

4.2.3 Complex Support

Next we consider examples with complex compact supports in 2D. We first take the initial data as

[TABLE]

which has a partial donut-shaped support. Similar examples are studied in [5, 43].

Fig. 4.9 shows the numerical solutions of the PME with $\alpha=3$ for the initial condition (4.10) at various time, which are computed by scheme 2 with $N=910$ and $\tau=10^{-2}$. One can see that our method works well for the concave domain. However, as pointed out in [5], a Lagrangian method can not handle the topological change automatically, which is also a limitation of our numerical approach. For this example, since the domain will evolve to reach a point where two ends of the “horseshoe” intersecting each other, the tangling of mesh cannot be avoided after this point. In order to get solutions beyond this point, one can manually interpolate the solution on to a new mesh, which can be viewed as an update to $\rho_{0}(\mathbf{X})$ (along with the mesh).

In our last numerical example, we study a peaks merge problem for the PME with $\alpha=4$, in which the initial data has two peaks, connected by a thin layer of mass. Similar numerical experiments are studied in [42, 11].

Let $\Omega=[-1,1]^{2}$, consider the initial data

[TABLE]

The boundary condition on $\partial\Omega$ is the Neumann boundary condition. In our numerical simulation, we take $\Omega_{0}=\Omega$ and manually set the velocity of the nodes on $\partial\Omega_{0}$ to be zero for the Neumann boundary condition. Fig. 4.10 shows the numerical solutions by scheme 2 at various time ($t=0.1,1$ and $5$) with a uniform initial mesh in $\Omega_{0}$ [$N=841$, $\tau=10^{-2}$]. It can be seen that the mesh is concentrated around the “interface”, in which the solution is steep, during the evolution. One can view our approach as a kind of moving mesh method. Unlike most of traditional moving mesh approach, in which both updating the mesh and solving the equation in the new mesh are required, we only need to solve the equation of the flow map, i.e., the equation of the mesh, the numerical solutions are determined by the kinematic relation. Similar to [11], it is a surprise that our numerical method can handle the situation of “peaks merge”, although local coarsening and remeshing is still needed in order to get better numerical solutions.

5 Summary

Structure-preserving and adaptive are two important aspects in designing efficient numerical methods for PDEs arising in numerous physical and biological modeling. It is usually a difficult task to construct a numerical scheme to a conservative or dissipative PDE that retain the conservation/dissipation properties in a discrete sense [20]. In this paper, we proposed a general framework to derive an efficient structure-preserving numerical scheme by employing a discrete energetic variational approach. A discrete energetic variational approach provides basis of deriving the “semi-discrete equations” after introducing a proper spatial discretization to the given energy-dissipation law, and can be applied to a large class of partial differential equations with energy-dissipation laws and kinematic relations, such as generalized diffusion equations, phase-field equations, and equations of liquid crystal. Within a piecewise linear approximation to the flow map, our approach is capable of handling high spatial dimensional situations. As an application, we develop two numerical schemes for the PME based on different energy-dissipation laws. By performing numerical experiments in both 1D and 2D, we show that the numerical scheme based on the energy-dissipation law (2.21) can better capture the free boundary and estimate the waiting time for the PME, without explicitly tracking the movement of the free boundary. Although our Lagrangian scheme performs well in numerical tests, a detailed numerical analysis of such schemes are required. We are working hard in this direction. As mentioned in the sect. 3, there is essential difference between Lagrangian schemes in 1D and high dimensions since both $\det F$ and the admissible set are not convex, which impose challenges in numerical analysis.

On the notion of our numerical approach, large deformations and topological changes can be difficult to handle by the dynamics of the flow map, usually with higher nonlinearity, degeneracy, or possible singularities. Furthermore, for general problems, a proper “initial data” may not be available. In order to solve these problems, we will employ a hybrid method, which solves the original equation in Eulerian and Lagrangian coordinates alternatively, in the ongoing work. In other word, during the evolution of the flow map, we can update the “initial data” (“reference data”) and the mesh, where the “initial data” is obtained by some Eulerian methods.

Finally, in the current approach, our piecewise linear approximation to the flow map is based on a finite element method, which provides us a simple framework to compute the deformation matrix $F$ on each element explicitly, but it is not obvious to incorporate the local coarsening and remeshing into it. One possible idea is to incorporate the particle methods [14] into our framework, in which remeshing for particle distortion can be easily dealt with, although how to compute the deformation matrix $F$ might still be a challenge.

Acknowledgment

The authors acknowledge the partial support of NSF (Grant DMS-1759536). We thank Prof. Cheng Wang, Prof. Xingye Yue, Dr. Chenghua Duan, Dr. Pei Liu and Dr. Qing Cheng for helpful discussions. Y. Wang would also like to thank Department of Applied Mathematics at Illinois Institute of Technology for their generous support and for a stimulating environment.

Appendix A Explicit form of $F_{e}$ in 2D

In the appendix, we give the explicit form of $F_{e}$ as a function of $a_{en(e,l)}$ and $b_{en(e,l)}$ ($l=1,2,3$) on each element $\tau_{e}$.

For a given element $\tau_{e}$, we denote the nodes of it by $\mathbf{X}^{e}_{l}$, and we let $\bm{\xi}^{e}_{l}=(a_{l}^{e},b_{l}^{e})\triangleq\bm{\xi}_{en(l,e)}$ ($l=1,2,3$). The discrete flow map on $\tau_{e}$ can be written as

[TABLE]

where $\lambda^{e}_{l}\in P_{1}(\tau_{e})$ is a nodal basis on $\tau_{e}$, which satisfies $\lambda^{e}_{l}(\mathbf{X}^{e}_{m})=\delta_{lm}$ ($l,m=1,2,3$).

We can compute $\lambda_{l}^{e}(\mathbf{X})$ by mapping $\tau_{e}$ into the reference triangle $\tau_{s}=\{\hat{\mathbf{X}}=(\hat{X},\hat{Y})\in\mathbb{R}^{2}~{}:~{}\hat{X}\geq 0,\hat{X}+\hat{Y}\leq 1\}$ with nodes $\hat{\mathbf{X}}_{1}^{s}=(0,0),\hat{\mathbf{X}}_{2}^{s}=(1,0)$, and $\hat{\mathbf{X}}_{3}^{s}=(0,1)$. The map between $\tau_{s}$ to $\tau_{e}$ is given by

[TABLE]

where

[TABLE]

Hence,

[TABLE]

Since the nodal basis on $\tau_{s}$ is given by

[TABLE]

we have

[TABLE]

where

[TABLE]

Then, we can compute $F_{e}$ as a function of $a_{l}^{e}$ and $b_{l}^{e}$ on $\tau_{e}$ directly by (1.1), which is

[TABLE]

where $X_{l-m}^{e}=X_{\l}^{e}-X_{m}^{e}$, $Y_{l-m}^{e}=Y_{l}^{e}-Y_{m}^{e}$ ($l,m=1,2,3$).

Appendix B Numerical Implementation in 2D

In this appendix, we present the detailed numerical implementation for the numerical scheme (3.18), which cannot be written down in an explicit form in a general mesh as $\frac{\delta\mathcal{A}}{\delta a_{i}}({\bm{a}}^{n+1},{\bm{b}}^{n+1})$ ($\frac{\delta\mathcal{A}}{\delta b_{i}}({\bm{a}}^{n+1},{\bm{b}}^{n+1})$) and ${\sf M}^{*}({\bm{a}}^{n},{\bm{b}}^{n})$ depend on the triangulation $\mathcal{T}_{h}$. Since we have the explicit form of $F_{e}$ as a function of $a_{en(e,l)}$ and $b_{en(e,l)}$ on each element $\tau_{e}$ (shown in the Appendix A), we can compute $\frac{\delta\mathcal{A}}{\delta a_{i}}({\bm{a}}^{n+1},{\bm{b}}^{n+1})$ ($\frac{\delta\mathcal{A}}{\delta b_{i}}({\bm{a}}^{n+1},{\bm{b}}^{n+1})$) and ${\sf M}^{*}({\bm{a}}^{n},{\bm{b}}^{n})$ in the numerical implementation by using the standard technique in finite element methods, that is, summing the results on each element over the mesh [34]. For instance, $\frac{\delta\mathcal{A}_{h}}{\delta a_{i}}(\bm{a}^{n+1},\bm{b}^{n+1})$ ($\frac{\delta\mathcal{A}_{h}}{\delta b_{i}}(\bm{a}^{n+1},\bm{b}^{n+1})$) and ${\sf M}^{*}(\bm{a}^{n},\bm{b}^{n})$ in our scheme 2 [scheme (3.18) with (3.48)] can be computed by Algorithm 1 and Algorithm 2. We can have an explicit form of $\frac{\delta\mathcal{A}_{h}}{\delta a_{i}}(\bm{a}^{n+1},\bm{b}^{n+1})$ ($\frac{\delta\mathcal{A}_{h}}{\delta b_{i}}(\bm{a}^{n+1},\bm{b}^{n+1})$) and ${\sf M}^{*}(\bm{a}^{n},\bm{b}^{n})$ in a uniform triangulation in a rectangular domain by applying Algorithm 1 and Algorithm 2.