Anisotropic Mesh Adaptation for Finite Element Solution of Anisotropic Porous Medium Equation
AAnisotropic Mesh Adaptation for Finite Element
Solution of Anisotropic Porous Medium Equation
Xianping Li ∗ Anisotropic Porous Medium Equation (APME) is developed as an extension of thePorous Medium Equation (PME) for anisotropic porous media. A special analytical so-lution is derived for APME for time-independent diffusion. Anisotropic mesh adaptationfor linear finite element solution of APME is discussed and numerical results for two di-mensional examples are presented. The solution errors using anisotropic adaptive meshesshow second order convergence.
AMS 2010 Mathematics Subject Classification.
Key words. porous medium equation, anisotropic mesh adaptation, adaptive mesh, anisotropicdiffusion, finite element, moving mesh
In this paper, we extend the porous medium equation (PME) to the Anisotropic Porous MediumEquation (APME) that takes into consideration the anisotropic physical properties of the porousmedia such as permeability. Then we study the linear finite element solution of APME. We considerthe following problem u t = ∇ · ( u m D ∇ u ) , in Ω T = Ω × ( t , T ] u ( x , t ) = 0 , on ∂ Ω × ( t , T ] ,u ( x ,
0) = u ( x ) , in Ω × { t = t } (1)where u = u ( x , t ) is a nonnegative scalar function, m ≥ ⊂ R d is aconnected polygonal or polyhedral domain of d -dimensional space, t ≥ T > u ( x ) ≥ D = D ( x ) is a general symmetricand strictly positive definite matrix-valued function on Ω T that takes both isotropic and anisotropicdiffusion as special cases. For simplicity, we consider only time-independent diffusion matrix D in thiswork. The principles can also be applied to the time-dependent situation with minor modifications.Porous medium equation (PME) arises in many fields of science and engineering such as fluid flowin porous media, heat transfer or diffusion, image processing, and population dynamics [36]. Thegeneral form is given as u t = ∇ · ( ∇ u m +1 ) , (2) ∗ Department of Mathematics and Statistics, the University of Missouri-Kansas City, Kansas City, MO 64110, U.S.A.( [email protected] ) a r X i v : . [ m a t h . NA ] A ug r in the modified form u t = ∇ · ( u m ∇ u ) . (3)For gas flow in porous media, m is the heat capacity ratio, u represents the density, u m is the pressure,and −∇ u m is the velocity.Mathematically, the parameter m in (3) can take any real value. Specifically, when m = 0, the PME(2) or (3) reduces to the heat equation. When m >
0, the PME becomes a nonlinear evolution equationof parabolic type that has attracted interests of both theoretical and computational mathematicians.Particularly, when m = 1, the PME is called the Boussinesq’s equation that models groundwater flowin a porous stratum.The nonlinear term u m in (3) induces the so-called nonlinear diffusion that brings up many challengesin the analysis of the PME. However, this nonlinear diffusion is not the physical property of the porousmedia such as permeability. In this paper, we generalize PME to Anisotropic Porous Medium Equation(APME) that takes into account the anisotropy of the physical property of the porous media. In themean time, APME can also be viewed as a nonlinear anisotropic diffusion problem.General PME (2) or (3) has been studied extensively both in theory [28, 1, 19, 12, 13, 2, 34, 36]and numerical approximations [31, 6, 5, 37, 10, 38, 30, 11, 33, 7, 8, 22, 26]. In particular, thenonlinear diffusion term and the sharp gradient near the free boundary make it difficult to achievehigh order convergence of the numerical solutions. For example, finite difference moving mesh methodhas been developed in [33, 22] for PME in one-dimensional space and second order convergence hasbeen observed; however, no result is provided for PME in 2D. Error estimates have been developedin [30] for finite volume discretization of PME in 2D, which shows only first order convergence. Forfinite element discretization using quasi-uniform meshes [31, 27, 32, 9, 37], the convergence is at mostfirst order for m = 1 and decreases for larger m . For one dimensional PME, high order convergencerate was achieved on a uniform mesh by using a high order local discontinuous Galerkin finite elementmethod [38].On the other hand, adaptive meshes and moving meshes are great choices to improve computationalefficiency and accuracy by concentrating more mesh elements in the regions where solution changessignificantly. In particular, Ngo and Huang [26] have studied moving mesh finite element solutionof PME and demonstrated the advantages of moving mesh over quasi-uniform meshes. Their resultsshow second-order convergence of solution errors using moving meshes.However, no result is currently available for APME. Interesting features of PME also appear inAPME such as finite propagation, free boundaries and waiting time phenomenon. Moreover, withthe anisotropy of the porous media, satisfaction of maximum principle becomes more challenging andspecial mesh adaptation is needed, see [21, 23, 24] and the references therein. This paper serves asa starting effort about APME and its numerical solutions. Anisotropic mesh adaptation techniqueis applied in the numerical computations to improve efficiency and accuracy. Different than movingmesh method that keeps the connectivity of the elements, our anisotropic mesh adaptation techniquecan change the connectivity as well as the number of elements as desired. Therefore, a coarse uniformmesh can be used as the initial mesh in our adaptation, while a fine initial mesh is usually needed formoving mesh method in order to capture the sharp change of solution on the initial free boundary.The outline of this paper is as follows. In Section 2, the anisotropic porous medium equation(APME) is developed and the exact solution for a special case is discussed. Section 3 gives a briefsummary of linear finite element solution of the APME, and Section 4 introduces the anisotropic meshadaptation methods. Numerical examples are presented in Section 5 to show the different behaviorbetween APME and PME and also demonstrate the advantages of adaptive anisotropic meshes overother meshes. Finally, some conclusions are drawn in Section 6.2 Anisotropic Porous Medium Equation (APME)
Firstly, we derive the model for fluid flow through anisotropic porous media as follows. Let u be thedensity of the fluid. The flow is governed by the following three equations.(I) Conservation of mass (continuity equation) εu t = −∇ · ( u v ) , (4)where ε ∈ (0 ,
1) is the porosity of the media and v is the velocity.(II) Darcy’s law in anisotropic porous media v = − µ D ∇ p, (5)where µ is the viscosity of the fluid, D is the permeability matrix of the porous media, and p is thepressure. In most cases, the porous media is anisotropic, thus D has different eigenvalues. If D varieswith location, then it also represents heterogeneous media.(III) The equation of state p = p u m , (6)where p is the reference pressure and m ≥ u t = p εµ ∇ · ( u D ∇ u m )= p εµ ∇ · ( D u · m · u m − ∇ u )= mp εµ ∇ · ( D u m ∇ u ) (7)= mp ( m + 1) εµ ∇ · ( D ∇ u m +1 ) . (8)Scale time t with the constant mp ( m +1) εµ in (8) or mp εµ in (7), we obtain the anisotropic porous mediumequation (APME) as follows u t = ∇ · ( D ∇ u m +1 ) (9)or u t = ∇ · ( u m D ∇ u ) . (10)Combining u m D in (10) together as the diffusion term, APME can be viewed as an anisotropic diffusionequation with induced nonlinearity from the solution. Thus, available results for anisotropic diffusionproblems [4, 35, 21, 23, 24] are useful for the study of APME.Note that, V´azquez has mentioned the extension of PME in non-homogeneous media (NHPME) in[36] as ε ( x , t ) u t = ∇ · ( c ( x , t ) ∇ u m +1 ) , (11)where ε ( x , t ) and c ( x , t ) are nonnegative functions. However, the anisotropy of the porous media hasnot been considered in NHPME. In this sense, both PME (2) and NHPME (11) are special casesof APME (10), where both the anisotropy and heterogeneity are taken into account in the diffusionmatrix D . 3ext, we derive a special solution for (1). We start with the Barenblatt-Pattle solution for PME(2) developed by Barenblatt [3] and Pattle [29] independently. Suppose the initial solution u ( x ) (attime t ) is compact-supported in a region given by r > u ( x ) = max (cid:40)(cid:18) − x T x r (cid:19) m , (cid:41) . (12)The solution at time t ≥ t for PME is given by u ( x , t ) = max (cid:40) κ d (cid:18) − x T x r κ (cid:19) m , (cid:41) , (13)where κ = (cid:18) tt (cid:19) β , β = 1 dm + 2 , and t = 12 βmr . (14)For APME, we consider the mesh T h in the physical domain Ω as a uniform mesh T c in the compu-tational domain Ω c specified by the metric M . By choosing M = D − , the diffusion can be consideredas isotropic in the computational domain Ω c with mesh T c [17, 23]. We apply the Barenblatt-Pattlesolution for PME in Ω c and then map the solution to Ω for APME. Denote the vertices in T c by ˜ x ,we have ˜ x = M x = D − x . (15)We also consider the initial solution in a region given by r > c as u ( x ) = max (cid:18) − ˜ x T ˜ x r (cid:19) m , = max (cid:40)(cid:18) − x T D − x r (cid:19) m , (cid:41) . (16)Then the solution at time t ≥ t for APME is given by u ( x , t ) = max κ d (cid:18) − ˜ x T ˜ x r κ (cid:19) m , = max (cid:40) κ d (cid:18) − x T D − x r κ (cid:19) m , (cid:41) (17)where κ , β , and t are the same as those defined in (14). R emark . If the initial solution of APME is taken as (12) instead of (16), the boundary ofthe compact support is not a circle in the computational domain Ω c . Therefore, the evolution of thesolution cannot be described by the solution (17). The difference will be demonstrated in Examples5.1 and 5.2. In this section, we briefly describe the linear finite element formulation of IBVP (1). Assume that anaffine family of simplicial triangulations {T h } is given for the physical domain Ω, and define U = { v ∈ H (Ω) | v | ∂ Ω = 0 } . T h by U h . Then a linear finite elementsolution u h ( t ) ∈ U h for t ∈ (0 , T ] to IBVP (1) is defined by (cid:90) Ω ∂u h ∂t v h d x + (cid:90) Ω ( ∇ v h ) T (cid:16) ( u h ) m D (cid:17) ∇ u h d x = 0 , ∀ v h ∈ U h , t ∈ ( t , T ] . (18)Denote the numbers of the elements, vertices, and interior vertices of T h by N , N v , and N vi , respec-tively. Assume that the first N vi vertices are the interior vertices. Then U h and u h can be expressedas U h ( t ) = span { φ , · · · , φ N vi } ,u h = N vi (cid:88) j =1 u j ( t ) φ j + N v (cid:88) j = N vi +1 u j ( t ) φ j , (19)where φ j is the linear basis function associated with the j th vertex, x j , at time t . The boundary andinitial conditions in (1) are approximated as u j ( t ) = 0 , j = N vi + 1 , ..., N v , (20)and u j ( t ) = u ( x j ) , j = 1 , ..., N v . (21)Substituting (19) into (18), taking v h = φ i ( i = 1 , ..., N vi ), and combining the resulting equationswith (20), we obtain the linear algebraic system M d u dt + A ( u h ) u = , (22)where u = ( u , ..., u N vi , u N vi +1 , ..., u N v ) T is the unknown vector and M and A are the mass andstiffness matrices, respectively. The entries of the matrices are given as follows. For j = 1 , ..., N v , m ij = (cid:82) Ω φ j φ i d x = (cid:80) K ∈T h (cid:82) K φ j φ i d x , i = 1 , ..., N vi , i = N vi + 1 , ..., N v (23) a ij = (cid:82) Ω ( ∇ φ i ) T ( u h ) m D ∇ φ j d x = (cid:80) K ∈T h (cid:82) K ( ∇ φ i ) T ( u h ) m D ∇ φ j d x , i = 1 , ..., N vi δ ij , i = N vi + 1 , ..., N v . (24)The system (22) is solved using the fifth-order Radau IIA method with a two-step error estimator[14]. The relative and absolute tolerances are chosen as 10 − and 10 − , respectively.Next, we apply the anisotropic mesh adaptation method to generate the mesh T nh for time t = t n ( n = 0 , , · · · ). In this section we briefly introduce the anisotropic mesh adaptation method used in the computationsfor this paper. More details about the method can be found in [16, 17, 23, 18]. We take the so-called M -uniform mesh approach where an adaptive mesh is viewed as a uniform one in the metric specified5y a tensor M = M ( x ) that is assumed to be symmetric and uniformly positive definite on Ω. It isshown in [17] that an M -uniform mesh T h satisfies | K | det( M K ) = σ h N , ∀ K ∈ T h (25)1 d tr (cid:0) ( F (cid:48) K ) T M K F (cid:48) K (cid:1) = det (cid:0) ( F (cid:48) K ) T M K F (cid:48) K (cid:1) d , ∀ K ∈ T h (26)where M K = 1 | K | (cid:90) K M ( x ) d x , σ h = (cid:88) K ∈T h det( M K ) | K | . (27)Condition (25) is called as the equidistribution condition which determines the size of K from det( M K ) ,while (26) is called the alignment condition which controls the shape and orientation of K .In this paper, we consider three choices of metric tensors for adaptive meshes. The first choice isbased on minimization of the H semi-norm of linear interpolation error and is given in [16] by M adap ( K ) = (cid:13)(cid:13)(cid:13) I + 1 α h | H K ( u h ) | (cid:13)(cid:13)(cid:13) det (cid:18) I + 1 α h | H K ( u h ) | (cid:19) − (cid:20) I + 1 α h | H K ( u h ) | (cid:21) , (28)where u h is the finite element solution, H K ( u h ) is a recovered Hessian of u h over K , | H K ( u h ) | is theeigen-decomposition of H K ( u h ) with the eigenvalues being replaced by their absolute values, and α h is a positive regularization parameter.The second and third choices are related to the diffusion matrix and are defined in [23] as M DMP ( K ) = D − K , ∀ K ∈ T h (29)and M DMP + adap ( K ) = (cid:18) α h B K (cid:19) d +2 det ( D K ) d D − K , (30)where D K = 1 | K | (cid:90) K D ( x ) d x ,B K = det ( D K ) − d (cid:107) D − K (cid:107) · (cid:107) D K | H K ( u h ) |(cid:107) ,α h = | Ω | (cid:88) K ∈T h | K | B dd +2 K d +2 d . The M -uniform meshes associated with M DMP and M DMP + adap satisfy the alignment condition (26),and the mesh elements are aligned along the principal diffusion direction that is the direction of theeigenvector corresponding to the largest eigenvalue of the diffusion matrix D .If only fixed mesh is used in the computations, a very fine initial mesh is needed in order tocapture the sharp change of solution at the initial free boundary. Even for moving mesh methods, aninitial fine mesh is usually needed since the connectivity of the mesh elements is fixed. For adaptivemeshes, however, they only need to concentrate elements around the initial free boundary or movingboundary. Therefore, an anisotropic initial mesh will significantly increase the efficiency and accuracyof the computation. The anisotropic initial mesh can be generated by adapting a uniform Delaunay6 dapted meshat t = t Update timeSolve PDE Check if t = T Yes Stop computationOutput resultsNoCompute metric tensorGenerate new meshInterpolate solutionfrom old mesh
Figure 1: Procedures for anisotropic mesh adaptation based on M -uniform mesh approachmesh according to a metric tensor M . The adaptation can iterate a few times to generate an initialmesh with good quality.Starting with the adapted initial mesh, the procedures for anisotropic mesh adaptation based on M -uniform mesh is shown in Fig. 1. Some numerical results are presented in the next section. R emark . In the mesh adaptation procedures shown in Fig. 1, we use linear finite elementinterpolation to map the solution from old mesh to the new adaptive mesh for the current time t n .Then the new mesh is used to solve PDE for the next time t n +1 . Higher order interpolations can alsobe used if needed. R emark . The adaptation for the mesh at a given time t n can be iterated a few times. However,the solution also needs to be interpolated correspondingly, which may reduce the accuracy of thesolution. Therefore, the mesh is only adapted once for each time step in our computations. R emark . Comparing to uniform meshes, the condition numbers of mass matrix M in (23)and stiffness matrix A in (24) for anisotropic meshes are affected by mesh non-uniformity. However,the condition numbers are still bounded. In anisotropic meshes, elements are usually concentrated ina small portion of the physical domain. With Jacobi preconditioning, the condition numbers of M and A for M -uniform meshes are comparable with Delaunay meshes [20] and the impact on iterativeconvergence is not significant. In this section we present numerical results obtained in two dimensions for three examples to demon-strate the different behavior between PME and APME as well as the significance of mesh adaptationin the computations. For comparison purpose, we consider four types of meshes. One is a fixed meshthat does not change through out the computations. The other three meshes are M adap mesh, M DMP mesh, and M DMP + adap mesh that have been introduced in Section 4.The fixed mesh is generated by splitting the rectangle domain into sub-rectangles, then each sub-rectangles are divided into four triangles by the diagonal lines. The adaptive meshes are generatedusing the Bidimensional Anisotropic Mesh Generator (BAMG) developed by Hecht [15] based on theDelaunay triangulation and local node movement. E xample . The first example is in the form of IBVP (1) for m = 1 on domain Ω = [ − , withdiffusion D = (cid:18) . . . . (cid:19) and u given in (16) with r = 0 .
5. The initial time is t = 0 . M adap mesh, if not specified otherwise, the regularization parameter in (28) is chosenas α h = 0 .
01. As can be seen, the elements in both M DMP mesh and M DMP + adap mesh are alignedalong the north-east direction. -3 -2 -1 0 1 2 3-3-2-10123 (a): Fixed mesh, N = 6 , -3 -2 -1 0 1 2 3-3-2-10123 (b): M adap mesh, N = 3 , -3 -2 -1 0 1 2 3-3-2-10123 (c): M DMP mesh, N = 3 , -3 -2 -1 0 1 2 3-3-2-10123 (d): M DMP + adap mesh, N = 4 , t = t and the final solution at t = 0 . M adap and M DMP + adap meshes at different times are shown in Fig. 4. Comparing with the initial meshes inFig. 2(b) and (d), it is clear that the free boundary moves outward in the shape of an ellipse in thephysical domain and the meshes are adapted accordingly. R emark . From the definition of metric tensors (28) and (30), the smaller α h (hence larger1 /α h ), the more elements will be concentrated around the region with sharp change of solution. Fig.5 shows M adap meshes obtained using α h = 1 at t = t and t = 0 .
06. Comparing with Fig. 2(b) and8 S o l u t i on y x
01 -1 -1-2 -2-3 -3 (a): t = t S o l u t i on y x
01 -1 -1-2 -2-3 -3 (b): t = 0 . α h = 0 .
01 provides better adaptation than α h = 1 for M adap .For convergence of solution errors, we choose the final time as T = ( t + 0 . / . L norm of the solution errors obtained from different meshes are plotted in Fig. 6, where M adap, and M adap, . represent the results from M adap meshes obtained using α h = 1 and α h = 0 .
01, respectively.For comparison purpose, the first order and second order reference rates are also plotted in Fig. 6and computed as 0 . / √ N and 1 /N , respectively.As can be seen from Fig. 6, the solution errors obtained from fixed mesh and M DMP mesh havefirst order convergence, with M DMP mesh providing smaller errors than fixed mesh. On the otherhand, the errors obtained from M adap (both M adap, and M adap, . ) and M DMP + adap meshes havesecond order convergence, with M adap, . mesh providing smallest errors.It is interesting to observe that with slight adaptation using M adap meshes with α h = 1 (that is, M adap, ), the convergence rate of the solution errors has already been improved to second order. Withbetter adaptation using α h = 0 .
01, the solution errors are not only of second order convergence butalso much smaller than others.The result demonstrates that adaptive meshes help the numerical solutions to achieve higher orderconvergence. In the mean time, M DMP + adap mesh also performs well as a combination of M DMP and M adap , and the errors are between those from M adap, and M adap, . meshes. R emark . Similar to the effect of α h for M adap meshes, the regularization parameter α h in (30)for M DMP + adap meshes can also be adjusted to improve the numerical solutions. E xample . The second example is the same as Example 5.1 except that the initial solution u is defined as in (12) where the initial free boundary is a circle with radius r in the physical domainΩ. The purpose of this example is to show the different behavior between our APME and generalPME. For PME, as we know from the exact solution (13), the free boundary will move outward in ashape of a circle. However, it is not the case for APME.Fig. 7 shows the initial and final meshes adapted from metric tensors M adap with α h = 0 .
01 and M DMP + adap . The solutions at different times are shown in Fig. 8. As can be seen from the meshes inFig. 7, the free boundary varies from a circle gradually to an ellipse. At a given time, the ellipticalboundary is different than the ones in Example 5.1 as shown in Fig. 4 (b) and (d). In particular, theeccentricity of the elliptical boundary is smaller than that in Example 5.1. The solution at t = 0 . (a): t = 0 . M adap mesh -3 -2 -1 0 1 2 3-3-2-10123 (b): t = 0 . M adap mesh -3 -2 -1 0 1 2 3-3-2-10123 (c): t = 0 . M DMP + adap mesh -3 -2 -1 0 1 2 3-3-2-10123 (d): t = 0 . M DMP + adap meshFigure 4: Example 5.1. M adap and M DMP + adap meshes at different times.For investigation purpose, the results for m = 2 using the same initial solutions are presented inFig. 9. The observations are similar to the results for m = 1 except that the free boundary has notmoved as much as that for m = 1 during the same time period. E xample . In the third example, we choose m = 6 and consider an initial solution with two-isolated support in domain Ω = [ − , . We also consider heterogeneous anisotropic diffusion. Theinitial solution is defined at t = 0 as u = , x ∈ Ω = (1 , × (0 , , , x ∈ Ω = ( − , − × (0 . , . , x ∈ Ω \ (Ω ∪ Ω ) . (31)10 (a): t = t -3 -2 -1 0 1 2 3-3-2-10123 (b): t = 0 . M adap meshes obtained with α h = 1 at different times. L n o r m o f e rr o r N: number of elements fixed M adap, M adap, . M DMP M DMP + adap first ordersecond order L n o r m o f e rr o r N: number of elements fixed M adap, M adap, . M DMP M DMP + adap first ordersecond order Figure 6: Example 5.1. L norm of solution errors obtained from different meshes.The end of time evolution is chosen as T = 120. The diffusion matrix D = D ( x ) is taken as D ( x ) = (cid:18) cos θ − sin θ sin θ cos θ (cid:19) (cid:18)
50 00 1 (cid:19) (cid:18) cos θ sin θ − sin θ cos θ (cid:19) , (32)where θ = θ ( x, y ) = π sin(0 . x ) cos(0 . y ) is the angle between the direction of the principle eigenvector11 (a): M adap mesh at t = t -3 -2 -1 0 1 2 3-3-2-10123 (b): M adap mesh at t = 0 . -3 -2 -1 0 1 2 3-3-2-10123 (c): M DMP + adap mesh at t = t -3 -2 -1 0 1 2 3-3-2-10123 (d): M DMP + adap at t = 0 . M adap and M DMP + adap meshes at different times, m = 1.of D and the positive x -axis. With this choice of θ , the primary diffusion direction changes at differentlocations.The initial and final solutions are shown in Fig. 10. The initial M adap and M DMP + adap meshes aredisplayed in Fig. 11. As can be seen, the elements of the M DMP + adap mesh are aligned along theprinciple diffusion directions at different places.Fig. 12 and Fig. 13 show the M adap mesh and the corresponding numerical solutions at differenttimes, respectively. It is observed that the two free boundaries start merging at around t = 60. R emark . Due to the challenges of the anisotropic and nonlinear diffusion in the APME, somenegative values may occur in the numerical solutions during the computations, which violates thediscrete maximum principle (DMP). The cut-off method [25] is applied to force all numerical solutionsto be nonnegative. Satisfaction of DMP for APME and the effect of cut-off are under investigation.12 S o l u t i on y x
01 -1 -1-2 -2-3 -3 (a): t = t S o l u t i on y x
01 -1 -1-2 -2-3 -3 (b): t = 0 . m = 1. -3 -2 -1 0 1 2 3-3-2-10123 (a): M adap mesh -3 -2 -1 0 1 2 3-3-2-10123 (b): M DMP + adap meshFigure 9: Example 5.2. M adap and M DMP + adap meshes at t = 0 . m = 2. S o l u t i on y x (a): t = 0 S o l u t i on y x (b): t = 120Figure 10: Example 5.3. Initial and final solutions, m = 6.13 (a): M adap mesh -3 -2 -1 0 1 2 3-3-2-10123 (b): M DMP + adap meshFigure 11: Example 5.3. Initial M adap and M DMP + adap meshes, m = 6. In the previous sections, we have generalized the Porous Medium Equation (PME) (2) or (3) toAnisotropic Porous Medium Equation (APME) (10) that takes into account the anisotropy and het-erogeneity of the physical properties of the porous media such as permeability. A special exact solutionfor APME is developed in (17) based on the Barenblatt-Pattle solution (13) for PME.Meanwhile, anisotropic mesh adaptation technique has been applied to obtain the finite elementsolutions of APME, which have helped improving both accuracy and efficiency of the computation.For general quasi-uniform meshes, the convergence rate for the solution error can be at mostly firstorder for m = 1. However, with adaptive M adap mesh or M DMP + adap mesh, we have attained secondorder convergence rate, as demonstrated by Example 5.1. Furthermore, for a same metric tensor,better adaptation and smaller errors can be achieved by adjusting the regularization parameter α h inthe computation.Our result is comparable with those obtained using moving mesh methods in [26] for PME butcan achieve the same accuracy with less number of mesh elements. Different than the moving meshmethod that keeps the connectivity of the mesh elements, the adaptive meshes, including M adap and M DMP + adap , not only change the connectivity but also can change the number of elements. Therefore,we can start with a coarse initial uniform mesh and adapt it to concentrate more elements around thefree boundaries. The adapted mesh is used as a better initial mesh for later computations.Numerical results also show different behavior between APME and PME, as demonstrated byExamples 5.1 and 5.2. For PME, the diffusion is the same in all directions, however, for APME,the diffusion is more significant along the principle diffusion direction. The merger between isolatedfree boundaries also occurs along the principle diffusion directions as demonstrated by Example 5.3.The challenges for general anisotropic diffusion problems also apply to APME such as satisfaction ofmaximum principle [21, 23, 24]. More investigation on properties of APME is needed. Acknowledgment.
This work was partially supported by the UMRB grant KDM56 (Li) from theUniversity of Missouri Research Board. 14 (a): t = 30 -3 -2 -1 0 1 2 3-3-2-10123 (b): t = 60 -3 -2 -1 0 1 2 3-3-2-10123 (c): t = 90 -3 -2 -1 0 1 2 3-3-2-10123 (d): t = 120Figure 12: Example 5.3. M adap meshes at different times, m = 6. References [1] D.G. Aronson. Regularity properties of flows through porous media.
SIAM J. Appl. Math. ,17:461–467, 1969.[2] D.G. Aronson, L.A. Caffarelli, and S. S. Kamin. How an initially stationary interface begins tomove in porous medium flow.
SIAM J. Appl. Math. , 14:639–658, 1983.[3] G.I. Barenblatt. On some unsteady motions of a liquid or a gas in a porous medium.
Prikl. Mat.Mekh. , 16 (1):67–78, 1952.[4] J. Brandts, S. Korotov, and M. Kˇr´ıˇzek. The discrete maximum principle for linear simplicial finiteelement approximations of a reaction-diffusion problem.
Linear Algebra and its Applications ,429:2344–2357, 2008. 15 (a): t = 30 -3 -2 -1 0 1 2 3-3-2-10123 (b): t = 60 -3 -2 -1 0 1 2 3-3-2-10123 (c): t = 90 -3 -2 -1 0 1 2 3-3-2-10123 (d): t = 120Figure 13: Example 5.3. Numerical solutions using M adap meshes at different times, m = 6.[5] C. Budd, G. Collins, W. Huang, and R.D. Russell. Self-similar numerical solutions of the porousmedium equation using moving mesh methods. Phil. Trans. R. Soc. Lond. A , 357:1047–1078,1999.[6] E. DiBendedetto and D. Hoff. An interface tracking algorithm for the porous medium equation.
Trans. Amer. Math. Soc. , 284 (2):463–500, 1984.[7] J.C.M. Duque, R.M.P. Almeida, and S.N. Antontsev. Convergence of the finite element methodfor the porous media equation with variable exponent.
SIAM J. Numer. Anal. , 51 (6):3483–3504,2013.[8] J.C.M. Duque, R.M.P. Almeida, and S.N. Antontsev. Application of the moving mesh method tothe porous medium equation with variable exponent.
Math. Comput. Simulation , 118:177–185,2015. 169] C. Ebmeyer. Error estimates for a class of degenerate parabolic equations.
SIAM J. Numer.Anal. , 35(3):10951112, 1998.[10] C. Ebmeyer and W.B. Liu. Finte element approximation of the fast diffusion and the porousmedium equations.
SIAM J. Numer. Anal. , 46(5):2393–2410, 2008.[11] E. Emmrich and D. Siska. Full discretization of the porous medium/fast diffusion equation basedon its very weak formulation.
Commun. Math. Sci. , 10:1055–1080, 2012.[12] B.H. Gilding and L.A. Peletier. On a class of similarity solutions of the porous media equation.
J. Math. Anal. Appl. , 55:351–364, 1976.[13] B.H. Gilding and L.A. Peletier. On a class of similarity solutions of the porous media equationii.
J. Math. Anal. Appl. , 57:522–538, 1977.[14] S. Gonz´alez-Pinto, J.I. Montijano, and S. P´erez-Rodr´ıguez. Two-step error estimators for implicitRunge-Kutta methods applied to stiff systems.
ACM Trans. Math. Softw. , 30 (1):1–18, 2004.[15] F. Hecht. Bidimensional anisotropic mesh generator software (bamg). , 2010.[16] W. Huang. Metric tensors for anisotropic mesh generation.
J. Comput. Phys. , 204:633–665, 2005.[17] W. Huang. Mathematical principles of anisotropic mesh adaptation.
Comm. Comput. Phys. ,1:276–310, 2006.[18] W. Huang and R.D. Russell.
Adaptive Moving Mesh Methods . Springer, New York, 2011.[19] A.S. Kalashnikov. The propagation of disturbances in problems of nonlinear heat conductionwith absorption.
USSR Comp. Math. and Math. Phys. , 14:70–85, 1974.[20] L. Kamenski, W. Huang, and H. Xu. Conditioning of finite element equations with arbitraryanisotropic meshes.
Math. Comput , 83:2187–2211, 2014.[21] D. Kuzmin, M.J. Shashkov, and D. Svyatskiy. A constrained finite element method satisfying thediscrete maximum principle for anisotropic diffusion problems.
J. Comput. Phys. , 228:3448–3463,2009.[22] T.E. Lee, M.J. Baines, and S. Langdon. A finite difference moving mesh method based onconservation for moving boundary problems.
Journal of Computational and Applied Mathematics ,288:1–17, 2015.[23] X. Li and W. Huang. An anisotropic mesh adaptation method for the finite element solution ofheterogeneous anisotropic diffusion problems.
J. Comput. Phys. , 229:8072–8094, 2010.[24] X. Li and W. Huang. Maximum principle for the finite element solution of time-dependentanisotropic diffusion problems.
Numer. Meth. PDEs , 29:1963–1985, 2013.[25] C. Lu, W. Huang, and E. VanVleck. The cutoff method for the numerical computation ofnonnegative solutions of parabolic PDEs with application to anisotropic diffusion and lubrication-type equations.
J. Comput. Phys. , 242:24–36, 2013.1726] C. Ngo and W. Huang. A study on moving mesh finite element solution of the porous mediumequation.
J. Comput. Phys. , 331:357–380, 2017.[27] R.H. Nochetto and C. Verdi. Approximation of degenerate parabolic problems using numericalintegration.
SIAM J. Numer. Anal. , 25(2):784–814, 1988.[28] O.A. Oleinik, A.S. Kalashnikov, and Y.-I. Chzou. The cauchy problem and boundary problemsfor equations of the type of unsteady filtration.
Izv. Akad. Nauk SSR Ser. Math. , 22:667–704,1958.[29] R.E. Pattle. Diffusion from an instantaneous point source with concentration dependent coeffi-cient.
Quart. Jour. Mech. Appl. Math. , 12:407–409, 1959.[30] I.S. Pop, M. Sep´ulveda, F.A. Radu, and O.P. Vera Villagr´an. Error estimates for the finitevolume discretization for the porous medium equation.
Journal of Computational and AppliedMathematics , 234:2135–2142, 2010.[31] M.E. Rose. Numerical methods for flows through porous media. i.
Mathematics of Computation ,40 (162):435–467, 1983.[32] J. Rulla and N.J. Walkington. Optimal rates of convergence for degenerate parabolic problemsin two dimensions.
SIAM J. Numer. Anal. , 33(1):56–67, 1996.[33] G. Scherer and M.J. Baines. Moving mesh finite difference schemes for the porous mediumequation.
Mathematics Report Series 1/2012 , pages Department of Mathematics and Statistics,University of Reading, UK, 2012.[34] S.I. Shmarev. Interfaces in multidimensional diffusion equations with absorption terms.
NonlinearAnal. , 53:791–828, 2003.[35] V. Thom´ee and L.B. Wahlbin. On the existence of maximum principles in parabolic finite elementequations.
Math. Comput. , 77:11–19, 2008.[36] J.L. V´azquez.
The porous medium equation - Mathematical theory . Oxford Mathematical Mono-graphs, The Clarendon Press Oxford University Press, Oxford, 2007.[37] D. Wei and L. Lefton. A priori l p error estimates for Galerkin approximations to porous mediumand fast diffusion equations. Mathematics of Computation , 68:971–989, 1999.[38] Q. Zhang and Z.-L. Wu. Numerical simulation for porous medium equation by local discontinuousGalerkin finite element method.