THIRD ORDER MAXIMUM-PRINCIPLE-SATISFYING DG SCHEMES Third Order Maximum-Principle-Satisfying DG schemes for Convection-Diffusion problems with Anisotropic Diffusivity DIFFUSIVITY
Hui Yu Hailiang Liu
math.NA cs.NA 65M06, 65M60, 65M12

Maximum-principle-satisfying discontinuous Galerkin schemes]Third order maximum-principle-satisfying DG schemes for convection-diffusion problems with anisotropic diffusivity




Tsinghua University, Yau Mathematical Sciences Center, Beijing, China 100084 huiyu@tsinghua.edu.cn
Iowa State University, Mathematics Department, Ames, IA 50011 hliu@iastate.edu
65M06, 65M60, 65M12
 


 
For a class of convection-diffusion equations with variable diffusivity, we construct third order accurate discontinuous Galerkin (DG) schemes on both one and two dimensional rectangular meshes.
The DG method with an explicit time stepping can well be applied to nonlinear convection–diffusion equations.  
It is shown that under suitable time step restrictions, the scaling limiter proposed in [Liu and Yu, SIAM J. Sci. Comput. 36(5): A2296–A2325, 2014] when coupled with the present DG schemes preserves the solution bounds indicated by the initial data, i.e., the maximum principle, while maintaining uniform third order accuracy.
These schemes can be extended to rectangular meshes in three dimension.
The crucial for all model scenarios is that  an effective test set can be identified to verify the desired bounds of numerical solutions. 
This is achieved mainly by taking advantage of the flexible form of the diffusive flux  and the adaptable decomposition of weighted cell averages. 
Numerical results are presented to validate the numerical methods and demonstrate their effectiveness.


[
    
Hui Yu and Hailiang Liu

    December 30, 2023
=============================






§ INTRODUCTION
 
In this paper, we construct and analyze third order maximum-principle-satisfying (MPS) discontinuous Galerkin (DG) schemes for the following problem,
 
    {[ M(x) ∂_t u + ∇· f(u) = ∇·(A∇ u),                    x∈ℝ^d,   t>0,;                  u(0,x)=u_0 (x),                           x∈ℝ^d. ].

This equation can be considered as a model for numerous physical problems. 
In this equation,  u=u(t, x) is a time-dependent unknown scalar function, f(u) is the smooth vector flux,  and d is the spatial dimension. We assume that the diffusion tensor A=A(x, u) is a symmetric and nonnegative definite matrix, and M(x) is strictly positive scalar function. In this paper, we present our results and analysis for d=1, 2, extension to three dimension can be done as well. 

Our concern for (<ref>) arises from several typical scenarios. 
The first case with f=0 is a model for heat conduction in a non-uniform body with M(x)=c(x)/ρ(x), where c(x) is the specific heat and ρ(x) the mass density. 
The heat flux is proportional to the temperature difference -A ∇ u, known as the Fourier law of heat conduction. 
Here A measures the ability of the material to conduct heat, called the thermal conductivity. 
In the second case with  M=1, we have the usual convection-diffusion equation, 

    ∂_t u + ∇· f(u)  = ∇·(A(x, u) ∇ u).

There are many interpretations and derivations from fluid dynamics and other application areas that motivate the convection-diffusion equation, 
such as Navier-Stokes equations and the porous medium equation.  The third case is the Fokker-Planck equation, 

    ∂_t ρ =∇·(∇ρ +∇ V(x) ρ),

which can be rewritten in terms of u=ρ e^V(x) as (<ref>) with A(x) = M(x)=e^-V(x), and f=0. This setting allows for many variants.  
	
The above three types of problems can be formulated under the model class (<ref>). 
One important solution to (<ref>) is the one that is bounded by two constants dictated by the initial data, leading to the so-called Maximum Principle (MP). 
In other words, if

    c_1 = min_x u_0,	   c_2 = max_x u_0,

then u(t,x) ∈ [c_1, c_2]  for any x∈ℝ^d and t>0.  

From the analytical viewpoint, the maximum principle is quite general but very significant due to its physical implications.  From the numerical viewpoint, it is widely recognized that maximum principle provides a valuable tool in proving solvability results (existence and uniqueness of discrete solutions), enforcing numerical stability, and deriving convergence results (a priori error estimates) for the sequence of approximate solutions.
Design of a high order scheme to preserve the maximum principle  is known a challenging task. Our goal is to better understand how a high order DG scheme can be constructed for (<ref>)  to respect the MPS property.
The main difficulty in the anisotropic case with a variable weight M(x) is the derivation of suitable sufficient conditions so that the weighted cell 
average stays in [c_1,c_2] during the time evolution. Such weighted cell average is essentially used for limiting the numerical solution into [c_1, c_2], 
without destroying accuracy.   



 §.§ Related work
 
Early discussion of the discrete maximum principle for the convection-diffusion equations includes the linear finite element solutions for parabolic equations <cit.>,  and recent developments <cit.>, as well as <cit.> by the Petrov–Galerkin finite element method to solve convection dominated problems. However, they are under a different framework. 

The present investigation involves the choice of  numerical fluxes and monotonicity of weighted cell averages in terms of point values. In the largest sense, the origins of  these ideas  go all the way back to monotone schemes for hyperbolic conservation laws.  

Indeed, for scalar conservation laws, i.e., (<ref>) with A=0,  many first order classical schemes can be shown to be MPS (other names of this sort include bound-preserving,  positivity preserving, or maximum-principle-preserving) since such low order accurate  schemes are usually monotone. On the other hand, the Godunov Theorem states that a linear monotone scheme is at most first order accurate for the convection equation <cit.>. 
To construct high order accurate MPS schemes for scalar convection, weak monotonicity in finite volume type schemes including DG methods was first used in <cit.>. 
Here by weak monotonicity it means that each cell average  is monotone with respect to point values in that cell, see e.g.  <cit.>.  
The main idea in their work  is to find sufficient conditions to preserve the desired bounds of cell averages by repeated convex combinations.  A simple and efficient local MPS limiter can then be used to control the solution values at test points without affecting accuracy and conservation. Together with strong stability preserving (SSP) Runge–Kutta or multistep methods <cit.>,  which are convex combinations of several formal forward Euler steps, a high order accurate finite volume or DG scheme can be rendered MPS with the limiter. 

For diffusion,  a  linear finite volume scheme can only be up to second order accurate in order to preserve the weak monotonicity, unless a non-conventional discretization is used in the scheme construction such as that in <cit.>. 
For DG methods, in general only second order accuracy can be obtained to feature the MPS property; see <cit.> for solving (<ref>) on triangular meshes. 

The only DG method known to  satisfy the weak monotonicity up to third order accuracy is the direct DG  (DDG) method introduced 
in <cit.>.   
Indeed, the special method  parameters in the DDG discretization allowed us to design in <cit.>  a third order MPS method for the linear Fokker-Planck equation (<ref>). 
A key idea in  <cit.> is the use of the non-logrithimic Landau formulation 

    M∂_t u= ∇·(M∇ u)    with  M(x) = e^-V(x) and  u=ρ/M,

so that the corresponding maximum principle on ρ(t,x): 

    c_1 e^-V≤ρ(0, x) ≤ c_2 e^-V⟹ c_1 e^-V≤ρ(t, x) ≤ c_2 e^-V  ∀ t > 0,

reduces to  

    c_1 ≤ u(0, x) ≤ c_2   ⟹ c_1 ≤ u(t, x) ≤ c_2   ∀ t >0.

With this reformulation, one can show that each weighted cell average is monotone in terms of point values under appropriate CFL conditions. The result in <cit.> is directly applicable to multi-dimensional diffusion on rectangular meshes.  However, it gets subtle to ensure the MPS property on unstructured meshes; we refer to <cit.> for a third order such DDG method to solve diffusion equations on unstructured triangular meshes.

Another approach towards a positivity-preserving scheme with high order accuracy is to use the local DG (LDG) method <cit.>, combined with some special positivity-preserving fluxes.
 Such an effort was first made in <cit.> for constructing high order accurate positivity-preserving DG schemes for compressible Navier–Stokes equations. Concerning further developments in this direction, we refer to <cit.> for solving convection–diffusion problems.   
 An MPS third-order LDG method  using overlapping meshes has been recently proposed in <cit.> for convection-diffusion equations.

One noteworthy alternative to enforce positivity in high order schemes is to take a convex combination of high order flux with a first order positivity-preserving one; the method has been applied to various high order schemes including finite difference,  finite volume,  and DG methods <cit.>, while rigorous justification of accuracy for such methods seems difficult. 



 §.§ Present investigation
 The main distinction of our present investigation from  the  above mentioned works is the use of weighted cell averages for both ensuring the weak monotonicity and applying the scaling limiter, with special attention on difficulty caused by the  anisotropic diffusivity.  
 
 The spatial discretization explored in this work is the DDG method introduced by Liu and Yan in <cit.>.  Besides the usual advantages of a DG method (see e.g. <cit.>), one main feature of the DDG method lies in numerical  flux choices for the solution gradient, which involve second order derivatives evaluated crossing cell interfaces (see (<ref>) below). With this choice, the obtained schemes are provably stable and optimally convergent as well as superconvergent for β_1≠0 <cit.>. Such method has also been successfully extended to various application problems, including Fokker-Planck type equations <cit.>, and the three dimensional compressible Navier–Stokes equation <cit.>. 

Built upon the work <cit.>, we present third order DG methods for solving the initial value problem (<ref>), with an application to nonlinear convection-diffusion equations of form (<ref>). 
Let us illustrate the main ideas via a simple one-dimensional equation subject to periodic boundary condition:

    M(x)∂_t u  = ∂_x(A(x)∂_xu).
  
The computational cell is denoted by I_j = [x_j-1/2,x_j+1/2],  where x_j+1/2's are the grid points.
Let τ and h be the time and space steps for a uniform mesh and n the index for the time stepping. 
The update on ⟨ u_h^n ⟩_j, the weighted cell average of the numerical approximation u_h^n, is given by 

    ⟨ u^n+1_h⟩_j =  ⟨ u^n_h ⟩_j - τ/h. A∂_x u_h^n|_∂ I_j   where ⟨ u ⟩_j:= 1/h∫_I_jM(x)u(x) dx.

Note that ∂_x u_h is the approximation to ∂_x u at ∂ I_j, the cell interface of I_j, given by 

    ∂_x u    =β_0/h[u] +{∂_x u} +β_1 h[∂_x^2 u].

As mentioned above, this form of the diffusive flux was originally introduced in <cit.> as part of the DDG method for solving the diffusion problem. 
If the flux parameters satisfy 

    β_0≥ 1    and   1/8≤β_1≤1/4,

then the procedure developed in <cit.> can be extended to the present setting to conclude that, under a suitable CFL condition, the simple Euler forward will keep the cell average  u̅^n_j =⟨ u_h^n⟩_j/⟨ 1 ⟩_j∈  [c_1, c_2] in each time step,  and the validity of the maximum principle when combined with a scaling limiter. 

For a two dimensional problem with A = ([ a c; c b ]), the DDG method on shape-regular Cartesian meshes with κ^-1≤Δ x/Δ y≤κ can be rendered MPS if

    β_0 ≥ 1+κ|c| L(L-1)/2 min{a, b}   and   1/8≤β_1 ≤1/4.

Here L is the number of Gauss-Lobatto points used in the numerical evaluation of involved integrals. 
With f(u) and A = A(x,y,u), the MPS DDG schemes are analyzed in Section <ref> where the parameter range and the CFL conditions are established.

The main conclusion is as follows: by applying the weighted MPS limiter introduced in <cit.>  to the DDG scheme designed here for (<ref>), with the time discretization by an SSP Runge-Kutta method (see <cit.>), we obtain a third  order accurate scheme solving (<ref>) satisfying the strict maximum principle in the sense that the numerical solution never goes out of the range [c_1, c_2] as indicated by the initial data. 



 §.§  Organization of the paper

The rest of the paper is organized as follows. In Section <ref>, we design the numerical method for one dimensional problems.  
We first formulate the  DDG scheme to solve the model heat equation, and prove the MPS property of the third order fully discretized DDG scheme, we then apply the result to show the MPS property for nonlinear convection-diffusion equations. 
Section <ref> is organized similarly for two dimensional problems on Cartesian meshes. 
In Section <ref>,  we present the MPS limiter, with which the algorithm is complete.  
In Section <ref>, numerical tests for the DDG method are reported, including examples from the heat equation, porous media equation 
the Buckley-Leverett equation, and two dimensional diffusion with the anisotropic diffusion.  
Concluding remarks are given in Section <ref>.




§ MPS SCHEMES IN ONE DIMENSION

We first investigate the MPS property for third order DDG schemes to weighted diffusion equations, and show how to apply the scheme obtained to nonlinear convection-diffusion equations. 
 


 §.§ The diffusion equation
 We begin with the heat equation of the form 

    M(x)∂_t u=∂_x (A(x)∂_x u),

with M(x)>0 and A(x)≥ 0 on the spatial domain Ω,  subject to initial data u_0(x) and the periodic boundary condition.  It is known that  the following maximum principle holds:

    if  c_1 ≤ u_0(x) ≤ c_2   ∀ x ∈Ω,  then  c_1 ≤ u(t, x) ≤ c_2   ∀ x ∈Ω, t>0.

In general, the problem can be either defined on a connected compact domain with proper boundary conditions, or it can involve the whole real line with solutions vanishing at the infinity.  
In our numerical scheme, we will always choose the spatial domain to be a connected interval. 
For simplicity, the periodic boundary conditions are applied. 

We partition the domain Ω by regular cells such that  Ω=⋃_j=1^NI_j with I_j = [x_j-1/2, x_j+1/2]. 
Denote h_j=x_j+1/2-x_j-1/2 and h=max_j h_j.  
We seek numerical solutions in the discontinuous piecewise polynomial space, 

    V_h={v ∈ L^2(Ω) |   v|_I_j∈ P^k(I_j),    j=1⋯, N}.

Here P^k(I_j) is the space of k-th order polynomials on I_j. Note that the functions in V_h can be double-valued at cell interfaces.  Hence notations v^- and v^+ are used for the left limit and right limit of v.  The jump of these two values, v^+-v^-, is denoted by [v], and the average by {v}. 

Throughout this paper we adopt the DDG numerical flux of the form 

    ∂_x v    =β_0/h_j+1/2[v] +{∂_x v} +β_1 h_j+1/2[∂_x^2 v]    with  h_j+1/2 = h_j + h_j+1/2,

when crossing the cell interface x_j+1/2, and (β_0, β_1) are in the range to be specified
so that the underlying third order scheme can weakly satisfy the maximum-principle.  
The parameter range was first identified in <cit.> for a third order DDG scheme to feature the MPS property for linear Fokker-Planck equations.  

For model equation (<ref>), we consider a (k+1)th-order DG scheme: to find u_h∈ V_h such that for any test function v∈ V_h,

    ∫_I_j M(x) ∂_t u_h v dx     =  -∫_I_j A(x) ∂_x u_h ∂_x v dx + . A [∂_x u_hv +(u_h -{u_h})∂_x v]|^x_j+1/2_x_j-1/2,

with the diffusive flux ∂_x u_h as defined in (<ref>). 
This DDG scheme with interface correction was proposed in <cit.> for the diffusion problem, as an improved version of that in <cit.>. Based on the DDG scheme in <cit.>,  we will have 

    ∫_I_j M(x) ∂_t u_h v dx     =  -∫_I_j A(x) ∂_x u_h ∂_x v dx + . A ∂_x u_hv|^x_j+1/2_x_j-1/2.


With the forward Euler time discretization, the weighted cell average for either of two DDG schemes evolves as

    ⟨ u^n+1_h ⟩_j=  ⟨ u^n_h ⟩_j + μ h. A ∂_x u_h^n|^x_j+1/2_x_j-1/2,

where μ = τ/h^2 is the mesh ratio.  For a concise presentation, a uniform mesh is assumed. Here and in what follows, we denote the time step length as τ. 
Note that for periodic boundary conditions considered in this paper, we take  

    u_h,1/2^-=u^-_h, N+1/2,    u_h, N+1/2^+=u^+_h, 1/2

in the DDG numerical flux formula (<ref>) when j=1/2, N+1/2.

We also use the notation 

    ⟨ q(ξ) ⟩_j= 1/2∫_-1^1M(x_j+h/2ξ)q(ξ)dξ,  for any  q(ξ)  on [-1, 1].

And the cell average u̅_j is
u̅_j = ⟨ u_h ⟩_j/⟨ 1 ⟩_j.
For j=1,⋯, N, we define 

    a_j= ⟨ξ-ξ^2 ⟩_j/⟨ 1-ξ⟩_j,    
    b_j=⟨ξ +ξ^2 ⟩_j/⟨ 1+ξ⟩_j,

and

    ω̃^1_j =    ⟨γ-ξ(1+γ)+ξ^2⟩_j /2(1+γ),
    ω̃^2_j =   ⟨ 1-ξ^2⟩_j /1-γ^2, 
    ω̃^3_j =   ⟨ -γ +ξ(1-γ)+ξ^2⟩_j /2(1-γ)

with the weight M(x)|_I_j=M(x_j+h/2ξ). We recall the following key result. 
<cit.>  ω̃^i_j >0 for i=1, 2,3 if and only if 

    γ∈ (a_j, b_j),

where a_j, b_j satisfy -1<a_j<b_j<1. 
 
 Here we only show a_j< b_j, which means that selection of γ is always ensured. A direct calculation gives 
 
    b_j - a_j = 2⟨ 1⟩_j⟨ξ^2⟩_j - ⟨ξ⟩_j^2/⟨ 1+ξ⟩_j⟨ 1-ξ⟩_j≥ 0,

where the numerator can be reformulated as 
 
    ∫_-1^1∫_-1^1(η^2-ξη)M̃(ξ) M̃(η)dξ dη +∫_-1^1∫_-1^1(ξ^2-ηξ)M̃(η)M̃(ξ)dη dξ
    
    =   ∫_-1^1∫_-1^1(ξ-η)^2 M̃(ξ)M̃(η)dξ dη >0,

where  M̃(ξ):= M(x_j+h/2ξ) has been used.
 

We thus have the following result.  
(k=2) The scheme (<ref>) 
with

    β_0≥ 1
      and  1/8≤β_1 ≤1/4

is maximum-principle-satisfying, namely, c_1≤u̅^n+1_j ≤ c_2 if u_h^n(x) is in [c_1, c_2] on  {S_j}_j=1^N, where

    S_j = x_j +h/2{-1, γ, 1}

with γ satisfying

    a_j <γ < b_j  and   |γ|≤  8β_1-1,

under the CFL condition
μ≤μ_0,
where μ_0 is given in (<ref>) below. 

  
Step 1.  Weighted integral decomposition. 
Define 

    p_j(ξ)=u_h(x_j+h/2ξ)  for ξ∈ [-1,1],

we see that in the case of  p_j(ξ)∈P^2[-1, 1], for any γ∈ (-1,1),  the unique interpolation of p_j at three points {-1, γ, 1} gives the following
 
    p_j(ξ)=(ξ-1)(ξ -γ)/2(1+γ)p_j(-1)+(ξ -1)(ξ +1)/(γ -1)(γ +1)p_j(γ)+(ξ +1)(ξ -γ)/2(1-γ)p_j(1).

This yields the following identity for the weighted average,

    ⟨ u_h⟩_j =ω̃^1_j p_j(-1) + ω̃^2_j p_j(γ) + ω̃^3_j p_j(1),

where ω̃^i_j given in (<ref>) are ensured positive by Lemma 2.1.

Step 2.  Flux representation.  A direct calculation gives

    h .∂_xu_h|_x_j+1/2
    =   α_3(-γ) p_j+1(-1)+α_2(-γ)p_j+1(γ)+α_1(-γ)p_j+1(1)
       - [α_1(γ) p_j(-1)+α_2(γ)p_j(γ)+α_3(γ)p_j(1)],

where

    α_1(γ)= 8β_1-1+γ/2(1+γ),   α_2(γ)=21-4β_1/1-γ^2,   α_3(γ)=β_0 +8β_1-3+γ/2(1-γ),

are all positive due to (<ref>) and (<ref>). 

Step 3. Monotonicity under some CFL condition. We now substitute (<ref>) and (<ref>) into (<ref>) to obtain

    ⟨ u_h^n+1⟩_j=R_j^n(M(·), μ, h, A)

with 

    R_j^n(M(·), μ, h, A)=    ⟨ u^n_h ⟩_j + μ(.Ah∂_xu_h^n|_x_j+1/2 - .Ah ∂_xu_h^n|_x_j-1/2) 
    
    =   [ω̃_j^1-μ(α_3(-γ) A_j-1/2 +α_1(γ) A_j+1/2) ]p_j(-1) 
       +[ω̃_j^2 - μ(α_2(-γ) A_j-1/2 +α_2(γ) A_j+1/2)]p_j(γ)
       +[ω̃_j^3 - μ(α_1(-γ) A_j-1/2 +α_3(γ) A_j+1/2) ]p_j(1) 
       +μ  A_j+1/2[α_3(-γ) p_j+1(-1)+α_2(-γ)p_j+1(γ)+α_1(-γ)p_j+1(1) ]
       +μ  A_j-1/2[α_1(γ) p_j-1(-1)+α_2(γ)p_j-1(γ)+α_3(γ)p_j-1(1) ].

Here it is understood that p_0(ξ):=p_N(ξ+|Ω|) and p_N+1(ξ+|Ω|)=p_1(ξ) for incorporating the periodic boundary conditions. Note also that ∑_i^3α_i(γ)=β_0= ∑_i^3α_i(-γ).   
Using the fact  ω̃_j^3(γ)=ω̃_j^1(-γ) and formulas for
α_2(γ),  ω̃_j^2, we see that if  μ is chosen to be smaller than μ_0 where 

    μ_0 = 
    	(max_1≤ j≤ N A(x_j+1/2)  )^-1min_1≤ j≤ N{ω̃_j^1(±γ)/α_3(∓γ) 
    +α_1(±γ) ,  ⟨ 1-ξ^2⟩_j/4(1-4β_1)},

(<ref>) is nondecreasing in the point values p_j(±1), p_j(γ), p_j± 1(±1), p_j± 1(γ), hence when these values are replaced with the lower and upper bounds c_1 and c_2 respectively, we have

    c_1∑_i=1^3ω̃_j^i ≤⟨ u_h^n+1⟩_j ≤ c_2 ∑_i=1^3ω̃_j^i,

since the terms with α_i's are cancelled out. Moreover the sum of ω̃_j^i is ⟨ 1⟩_j. Therefore

    c_1⟨ 1⟩_j ≤⟨ u_h^n+1⟩_j ≤ c_2⟨ 1⟩_j   ⇒   c_1 ≤u̅_j^n+1≤ c_2.



 
For other types of boundary conditions, the boundary flux needs to be modified.  A  similar result to Theorem <ref> may be established as long as the PDE problem satisfies a maximum principle.


The use of γ is essential for the success of our schemes, in particular when M(x) is a function.  
For the special case M(x)=1, we have γ∈ (-1/3, 1/3), hence γ=0, which corresponds to the usual Gauss quadrature point, is admissible.  
The CFL number given in (<ref>) may be optimized by carefully tuning γ∈ (a_j, b_j), but not in a linear fashion.  
Nevertheless, it was observed in <cit.> that  the larger |γ| is, the better the scheme's performance for the heat equation. 




 §.§ Application to nonlinear convection-diffusion equation

In this section we will demonstrate how to apply our MPS DG method in section 2.1 to the nonlinear convection-diffusion equation, 

    ∂_t u +∂_x f(u) =∂_x (A(x, u)∂_x u),

where f(u) is a smooth function and diffusion coefficient A(x, u)≥ 0, subject to initial data u(0, x)=u_0(x), and periodic boundary conditions.  By applying the DG approximation, we obtain the following scheme. We seek u_h∈ V_h such that for any test function v∈ V_h,

    ∫_I_j∂_t u_h v dx     = ∫_I_j f(u_h)∂_xv dx -f̂(u_h^-, u_h^+) v |^x_j+1/2_x_j-1/2
          -∫_I_j A_h ∂_x u_h ∂_x v dx + .{A_h}[∂_x u_hv +(u_h -{u_h})∂_x v]|^x_j+1/2_x_j-1/2.
 
For diffusion part,  we adopt the DDG diffusive flux (<ref>) 
For convection, any monotone numerical flux can be used, i.e., f̂(u, v) is Lipschitz continuous, nondecreasing in u and nonincreasing in v, consistent with f(u) in the sense that f̂(u, u)=f(u). For example, the global Lax-Friedrichs flux 

    f̂(u_h^-, u^+_h)=1/2 (f(u^-_h) +f(u^+_h) -σ (u^+_h-u^-_h)), 
      σ=max_u ∈ [c_1, c_2]|f'(u)|.


We consider the first order Euler forward temporal discretization of (<ref>) to obtain 

    ∫_I_ju_h^n+1-u_h^n/τ v dx     = ∫_I_j f(u^n_h)∂_xv dx 
    	- f̂((u_h^n)^-, (u_h^n)^+) v |^x_j+1/2_x_j-1/2
          -∫_I_j A_h^n ∂_x u^n_h ∂_x v dx + .{A_h^n}[∂_x u^n_hv +(u^n_h -{u^n_h})∂_x v]|^x_j+1/2_x_j-1/2.

By taking the test function v = 1 on I_j and 0 elsewhere, we obtain the evolutionary update for the cell average, 

    u̅^n+1_j=  u̅^n_j - λf̂^n |^x_j+1/2_x_j-1/2  
     + μ h {A_h^n}∂_xu_h^n|^x_j+1/2_x_j-1/2,

where λ = τ/h  and μ = τ/h^2 are the mesh ratios.  Assuming that u̅^n_j ∈ [c_1, c_2] for all j's, we would like to derive some sufficient conditions such that 
u̅^n+1_j ∈ [c_1, c_2] under certain restrictions on λ and μ. 

For piecewise quadratic polynomials, the main result can be stated as follows.  
(k=2) The scheme (<ref>) with

    β_0≥ 1
      and  1/8≤β_1 ≤1/4

is maximum-principle-satisfying; namely, u̅_j^n+1∈ [c_1, c_2]  if u_h^n(x) ∈ [c_1, c_2] on the set S_j's where

    S_j = x_j +h/2{-1, γ, 1}

with γ satisfying

    -1/3 <γ <  1/3  and   |γ|≤  8β_1-1,

under the CFL condition

    λ≤λ_0,   μ≤μ_0

for some λ_0 and μ_0 defined in (<ref>) and  (<ref>), respectively. 

 We present the proof in four steps:  

Step 1.  Split: we split the average u̅_j^n into two halves so that   

    u̅_j^n+1=1/2𝒞_j^n +1/2𝒟_j^n,

where the convection term is 

    𝒞_j=u̅_j -2λf̂ |^x_j+1/2_x_j-1/2

and the diffusion term is 

    𝒟_j =u̅_j + 2μ h {A_h^n}∂_xu_h |^x_j+1/2_x_j-1/2.

This split is for convenient presentation and may not lead to an optimal CFL condition. 

Step 2. The integral decomposition. From (<ref>)  it follows 

    u̅_j = ω^1 p_j(-1)+ω^2 p_j(γ)+ω^3 p_j(1),

where ω^i=ω̃^i, ⟨ u_h⟩_j = u̅_j for M(x)≡1, and  

    ω^1 = 1+3γ/6(1+γ),  ω^2 = 2/3(1-γ^2),  ω^3 = 1-3γ/6(1-γ).

These coefficients are positive for γ satisfying (<ref>). 


Step 3.   The convection term.

Using the cell average decomposition and the flux formula,  we rewrite (<ref>) as  

    𝒞_j =    u̅^n_j -2 λ(.f|_x_j+1/2 - .  f|_x_j-1/2) 
    
    =   ω^3 p_j(1) - 2λf̂(p_j(1), p_j+1(-1)) + ω^2 p_j(γ)  + ω^1 p_j(-1) 
    + 2λf̂(p_j-1(1), p_j(-1)) 
    
    =    : G(p_j-1(1), p_j(-1), p_j(γ), p_j(1), p_j+1(-1)).

For a monotone flux f̂(u, v) being Lipschitz continuous with Lipschitz constant ℒ, 
G is non-increasing in all the four arguments provided the following condition is met, 

    2 λℒ≤min{ω^1, ω^3}.
 
Note that for the Lax-Friedrichs flux, ℒ = max_u∈[c_1,c_2] |f'(u)|.  
Moreover, the consistency of the flux f̂(u, u)=f(u) implies that G(u, u, u, u)=u.  
Hence we have 

    𝒞_j ∈ [G(c_1, c_1, c_1, c_1, c_1), G(c_2, c_2, c_2, c_2, c_2)]=[c_1, c_2],

as long as the involved values are in [c_1, c_2]. It suffices to take 

    λ_0= 1/2ℒmin{ω^1, ω^3}=1-3γ/12ℒ(1-γ),

where we have used the fact that ω^1(γ)=ω^3(-γ).

Step 4. The  diffusion term. We apply the result in Theorem <ref> to the case with M(x)≡1 and μ replaced by 2μ to conclude that 
𝒟_j ∈ [c_1, c_2] if u_h^n(x) ∈ [c_1, c_2] on S_j and μ≤μ_0 with 

    μ_0     = 1/2(max_j {A_h^n}_j+1/2)^-1min_1≤ j≤ N{ω^1(±γ)/α_3(∓γ) 
    +α_1(±γ) ,  1/3(1-4β_1)}
        =1/12(max_1≤ j≤ N{A_h^n}_j+1/2)^-1min_1≤ j≤ N{1± 3γ/β_0(1±γ)+8β_1-2, 2/1-4β_1}.



The above analysis can be readily carried over to (<ref>) in one dimension, i.e.,   

    M(x)∂_t u +∂_x f(u) =∂_x (A(x, u)∂_x u).

We summarize the result in the following. 
(k=2) The scheme (<ref>) when applied to (<ref>) with

    β_0≥ 1
      and  1/8≤β_1 ≤1/4

is maximum-principle-satisfying; namely, u̅_j^n+1∈ [c_1, c_2]  if u_h^n(x) ∈ [c_1, c_2] on the set S_j's where

    S_j = x_j +h/2{-1, γ, 1}

with γ satisfying

    γ∈ (a_j, b_j)  as in (<ref>) and   |γ|≤  8β_1-1,

under the CFL condition

    λ≤λ_0,   μ≤μ_0

for some λ_0 and μ_0 defined in (<ref>) and  (<ref>), respectively. 

 
The proof is entirely analogous to the proof of Theorem <ref>. One only needs to replace ω^i by ω̃^i for i=1, 2, 3, and  u̅_j by ⟨ u_h⟩_j, respectively, so to obtain different CFL conditions.  More precisely,  (<ref>) in Step 3 needs to be  replaced by 

    λ_0= 1/2ℒmin{ω̃^1(γ), ω̃^3(γ)}= 1/2ℒmin{ω̃^1(±γ)},

and (<ref>) by a factor of 1/2 gives    

    μ_0 = 1/2(max_1≤ j≤ N A(x_j+1/2)  )^-1min_1≤ j≤ N{ω̃_j^1(±γ)/α_3(∓γ) 
    +α_1(±γ) ,  ⟨ 1-ξ^2⟩_j/4(1-4β_1)}.

 



§ MPS SCHEMES IN TWO DIMENSIONS

In this section, we design an MPS DG scheme to solve two dimensional problems on Cartesian meshes.
Consider the model equation

    M(x,y)∂_t u = ∇·(A∇ u)  with  
    A = ([ a c; c b ]),    for  (x, y)∈Ω⊂ℝ^2,    t>0,

subject to the initial data u_0(x, y) and periodic boundary conditions.  Here M(x, y) > 0 is a given function, 
a, b and c are constant parameters so that A is nonnegative definite.
The domain Ω =I × J is a rectangle given by two intervals I and J in x and y direction, respectively.
Let 
∪_i=1^N_x∪_j=1^N_y K_ij be a partition of the domain Ω, with K_ij=I_i× J_j, where 

    I_i=[x_i-1/2, x_i+1/2],    J_j = [y_j-1/2, y_j+1/2].

The finite element space is defined as 

    V_h={v∈ L^2(Ω),  v|_K_ij∈ Q^k(K_ij), i=1, ⋯, N_x, j=1,⋯, N_y}.

Here Q^k(K_ij) is the tensor product space of P^k(I_i) and P^k(J_j).   Hence the DDG scheme can be formulated 
as follows: to find u_h∈ V_h such that for all v∈ V_h, 

    ∫_K_ij Mu_h^n+1v dxdy 
    =    ∫_K_ij Mu_h^nv dxdy 
    - τ∫_K_ij A∇ u_h^n·∇ v dxdy 
       + τ∫_∂ K_ijA∇ u_h^n·ν v ds 
     + τ∫_∂ K_ijA∇ v·ν (u_h^n - {u_h^n}) ds,

where ν is the outward unit normal to the cell boundary ∂ K_ij, and the numerical flux 

    ∇ u_h= (∂_x u_h, ∂_y u_h)^⊤

is defined as follows, 

    .∂_x u_h|_(x_i+1/2,y)    =β_0/Δ x[u_h] +{∂_x u_h} +β_1Δ x [∂_x^2 u_h],  .∂_y u_h|_(x_i+1/2,y)  ={∂_y u_h}, 
    .∂_y u_h|_(x,y_j+1/2)    =β_0/Δ y[u_h] +{∂_y u_h} +β_1Δ y [∂_y^2 u_h],  .∂_x u_h|_(x,y_j+1/2) ={∂_x u_h},

where β_0, β_1 are flux parameters to be determined to ensure the desired MPS property.  
Note that in ∂_y u_h|_(x_i+1/2,y), jump terms do not show up since along interface 
 x = x_i+1/2 and y ∈ [y_j-1/2, y_j+1/2], there is no jump of polynomials  in y direction. This argument applies to ∂_x u_h|_(x,y_j+1/2) as well.
Here for a concise expression of the numerical flux, a uniform mesh has been used, with Δ x = x_i+1/2-x_i-1/2 and Δ y=y_j+1/2-y_j-1/2.  

To proceed, we  recall some conventions similar to the one-dimensional case.  The weighted cell average is defined as 

    ⟨ u^n+1⟩_ij = ∫_K_ijMu_h dxdy/Δ x Δ y
    	= _J_j_I_iMu_h dxdy,

where  denotes the average integral.  We also define the weighted interval average in x and y, respectively,

    ⟨ϕ(ξ)  ⟩_i(y)=_-1^1M(x_i+Δ x/2ξ, y)ϕ(ξ)dξ,   ⟨ϕ(η) ⟩_j(x)=_-1^1M(x, y_j+Δ y/2η)ϕ(η)dη.

The cell average 

    u_ij = ∫_K_ijMu_h dxdy/∫_K_ijM dxdy=⟨ u_h⟩_ij/⟨ 1⟩_ij

update can be obtained from  (<ref>) as

    ⟨ u^n+1_h⟩_ij =   ⟨ u^n⟩_ij 
    + μ_xΔ x ._J_j(a∂_x u_h^n + c∂_y u_h^n) dy|_∂ I_i
    + μ_yΔ y ._I_i(b∂_y u_h^n + c∂_x u_h^n) dx|_∂ J_j,

where μ_x = τ/(Δ x)^2 and μ_y = τ/(Δ y)^2. 
Let μ = μ_x + μ_y and decompose ⟨ u^n⟩_ij as

    ⟨ u^n⟩_ij = μ_x/μ⟨ u^n⟩_ij  + μ_y/μ⟨ u^n⟩_ij,
 
so that (<ref>) can be rewritten as 

    ⟨ u^n+1⟩_ij =   μ_x/μ_J_jH_1(y) dy + μ_y/μ_I_i H_2(x) dx  +B,

where 

    H_1(y)  =_I_iM(x,y)u_h^n dx
    	+ μΔ x. a∂_x u_h^n|_∂ I_i,
        H_2(x)  =_J_jM(x,y)u_h^n dy
    	+ μΔ y. b∂_y u_h^n|_∂ J_j,
        B =  cτ/|K_ij|[ ∫_J_j{∂_y u_h}dy |_∂ I_i +∫_I_i{∂_x u_h}dx |_∂ J_j].

Notice that B can be expressed as a combination of point values of u_h^n at four vertices of K_ij as 

    B    =cτ/2Δ xΔ y(2u_h^n(x_i+1/2^-, y_j+1/2^-) - 2u_h^n(x_i+1/2^-, y_j-1/2^+)
    -2u_h^n(x_i-1/2^+, y_j+1/2^-) + 2u_h^n(x_i-1/2^+, y_j-1/2^+)) 
        + cτ/2Δ xΔ y(u_h^n(x_i+1/2^+, y_j+1/2^-) +u_h^n(x_i+1/2^-, y_j+1/2^+)  - u_h^n(x_i+1/2^-, y_j-1/2^-)
    -u_h^n(x_i+1/2^+, y_j-1/2^+))
        + cτ/2Δ xΔ y(u_h^n(x_i-1/2^+, y_j-1/2^-)+u_h^n(x_i-1/2^-, y_j-1/2^+) -u_h^n(x_i-1/2^+, y_j+1/2^+) 
    -u_h^n(x_i-1/2^-, y_j+1/2^-) ).

The two integrals in (<ref>) can be approximated by the Gauss-Lobatto quadrature rule with sufficient accuracy. 
Let us assume that we use an L-point Gauss-Lobatto quadrature rule with L≥k+3/2 points, which has accuracy of at least O(h^k+2).
Let

    Ŝ_i^x = {x_i-1/2 = x̂^1_i < ⋯ < x̂_i^σ < ⋯ < x̂^L_i=x_i+1/2}, 
       Ŝ_j^y = {y_j-1/2 = ŷ^1_j < ⋯ < ŷ_j^σ < ⋯ < ŷ^L_j=y_j+1/2}

denote the quadrature points on I_i and J_j, respectively, 
and ω̂^σ's be the associated quadrature weights so that

    ∑_σ=1^L ω̂^σ=1.

Using the quadrature rule on the right-hand side of (<ref>), we obtain the following scheme

    ⟨ u^n+1_h ⟩_ij= μ_x/μ∑_σ=1^L ω̂^σH_1(ŷ_j^σ) + μ_y/μ∑_σ=1^L ω̂^σH_2(x̂_i^σ)+B.

Also set 

    S^x_i     = {x_i-1/2, x_i^γ, x_i+1/2} = x_i + Δ x/2{-1, γ^x, 1}, 
    
     S^y_j     = {y_j-1/2, y_j^γ, y_j+1/2} = y_j + Δ y/2{-1, γ^y, 1}

 to denote the test sets on [x_i-1/2, x_i+1/2] and [y_j-1/2, y_j+1/2], respectively,  with γ^x, γ^y satisfying

    ⟨ξ-ξ^2⟩_i/⟨ 1- ξ⟩_i (y_j^σ) < γ^x < ⟨ξ +ξ^2⟩_i/⟨ 1+ξ⟩_i (y_j^σ),    |γ^x|≤ 8β_1-1,
       ⟨η-η^2⟩_j/⟨ 1- η⟩_j (x_i^σ) < γ^y < ⟨η +η^2⟩_j/⟨ 1+η⟩_j (x_i^σ),     |γ^y|≤ 8β_1-1.

We use ⊗ to denote the tensor product and define

    S_ij=(S_i^x⊗Ŝ^y_j)∪ (Ŝ_i^x ⊗ S_j^y).

The main result can now be stated in the following.
(k=2)
Consider the two dimensional DDG scheme (<ref>) on rectangular meshes.
Assume the mesh is regularly shaped, i.e., κ^-1≤Δ x/Δ y≤κ,  for some constant κ > 0, and the flux parameters (β_0, β_1)  satisfy 

    β_0 ≥ 1+κ|c|/2 ω̂^1 min{a, b}.
       and   1/8≤β_1 ≤1/4

If u_h^n(x, y)∈ [c_1, c_2] for all (x, y)∈ S_ij, then there exists μ_0>0 such that if μ≤μ_0 the cell average
u̅^n+1_ij∈ [c_1, c_2]. 
More precisely, we have 

    μ_0 = min_i,jω_ijmin{ω̂^1/ω̂^1max{a,b}(β_0 + 8β_1-2/1+γ) + κ|c|,
    	1-γ^2/4max{a,b}(1-4β_1)},

where ω_ij>0 is defined in equation (<ref>), and γ = max{|γ^x|, |γ^y|}≤ 8β_1-1.


The proof is relegated to Appendix A.  We proceed to deal with nonlinear diffusion equations in the next  subsection.




 §.§ Application to nonlinear diffusion equations
 

This section is devoted to application to nonlinear diffusion equations of the form 

    ∂_t u = ∇·(A∇ u)  with  
    A(x,y,u) = ([ a c; c b ])    for  (x, y)∈Ω⊂ℝ^2,    t>0,

subject to initial data u_0(x, y) and periodic boundary conditions, and A(x,y,u) is nonnegative definite. 
This type of model arises in a wide range of applications. 


Hence the DDG scheme can be formulated as follows: to find u_h∈ V_h such that for all v∈ V_h, 

    ∫_K_ij u_h^n+1v dxdy 
    =    ∫_K_ij u_h^nv dxdy 
    - τ∫_K_ij A_h^n ∇ u_h^n·∇ v dxdy 
       + τ∫_∂ K_ij{A_h^n}∇ u_h^n·ν v ds 
     + τ∫_∂ K_ij{A_h^n}∇ v·ν (u_h^n - {u_h^n}) ds,

where ν is the outward unit normal to the cell boundary ∂ K_ij, the numerical flux 
∇ u_h
is defined in (<ref>), and A_h^n=A(x, y, u_h^n). 

The  cell average evolves according to 

    u̅^n+1_ij = u̅^n_ij 
    + μ_xΔ x ._J_j({a_h^n}∂_x u_h^n + {c_h^n}∂_y u_h^n) dy|_∂ I_i
    + μ_yΔ y ._I_i({b_h^n}∂_y u_h^n + {c_h^n}∂_x u_h^n) dx|_∂ J_j.

That is 

    u̅^n+1_ij  =  μ_x/μ_J_jH_1(y) dy + μ_y/μ_I_i H_2(x) dx  +B,

where  μ_x = τ/(Δ x)^2,  μ_y = τ/(Δ y)^2 and μ = μ_x + μ_y,  with 

    H_1(y)  =_I_iu_h^n dx
    	+ μΔ x. {a_h^n}∂_x u_h^n|_∂ I_i,
        H_2(x)  =_J_ju_h^n dy
    	+ μΔ y. {b_h^n}∂_y u_h^n|_∂ J_j,
        B =  τ/|K_ij|[ ∫_J_j{c_h^n}{∂_y u_h}dy |_∂ I_i +∫_I_i{c_h^n}{∂_x u_h}dx |_∂ J_j].


The main result  can be stated in the following.
(k=2) Consider the two dimensional DDG scheme (<ref>) on rectangular meshes.
Assume the mesh is regularly shaped, i.e., κ^-1≤Δ x/Δ y≤κ,  for some constant κ > 0, and the flux parameters (β_0, β_1)  satisfy 

    β_0 ≥ 1+2κc_∞ L(L-1) /(1-γ) min{a, b},
       and   1/8≤β_1 ≤1/4

and γ = max{|γ^x|, |γ^y|}≤ 8β_1-1,
where

    c_∞ = max_(x,y) ∈Ω, u∈[c_1, c_2]|c(x,y,u)|,   min{a, b} = min_(x,y) ∈Ω, u∈[c_1, c_2]{a(x,y,u), b(x,y,u)}.

If u_h^n(x, y)∈ [c_1, c_2] for all (x, y)∈ S_ij, then there exists μ_0>0 such that if μ≤μ_0 the cell average
u̅^n+1_ij∈ [c_1, c_2]. 
More precisely, we have 

    μ_0 = min{1-3γ/6max{a,b}(β_0 + 8β_1-2/1+γ)(1-γ) + 12κc_∞ L(L-1),
    	1/6max{a,b}(1-4β_1)},

where 

    max{a, b} = max_(x,y) ∈Ω, u∈[c_1, c_2]{a(x,y,u), b(x,y,u)}.


The proof is relegated to Appendix B. 



§ SCALING LIMITER AND THE MPS ALGORITHM



 §.§ Scaling limiter
  The one dimensional result in Theorem <ref> and the two dimensional result  in Theorem <ref>
tell us that for the DDG scheme with forward Euler time discretization, we need to modify u_h^n such that it is in [c_1, c_2] on the test set S=S_j or S_ij.  
In one dimensional case, we can use the following  scaling limiter

    ũ_h(x) = θ(u_h(x) - u̅_j) + u̅_j    with θ = min{1,  |u̅- c_1/u̅_j - m_1|,  
    |c_2-u̅_j/m_2- u̅_j| },

where

    m_1 = min_x∈ S_ju_h(x),      m_2 = max_x∈ S_ju_h(x).

In the two dimensional case, 

    ũ_h(x, y) = θ(u_h(x, y)-u̅_ij) + u̅_ij   where θ = min{1,  | u̅_ij-c_1/u̅_ij - m_1|, c_2-u̅_ij/m_2-u̅_ij| },

 where 
 
    m_1 = min_(x, y)∈ S_ij u_h(x, y)  and    m_2=max_(x, y)∈ S_ij u_h(x, y).

The modified polynomials  are indeed in [c_1, c_2] and preserve the cell average.  
Moreover, following <cit.> it can be shown that  the above scaling limiters do not destroy the accuracy. We summarize this for two dimensional case only.  
 If u̅_ij∈ (c_1, c_2),  then the modified polynomial (<ref>) is as accurate as u_h(x,y) to approximate u(t,x,y) for the same t 
in the following sense:

    |u_h(x, y)-ũ_h(x, y)| ≤ C_k u_h(·, ·) - u(t, ·, ·)_∞,

where C_k is a constant depending on the polynomial degree k and the weight function M(x,y).




 §.§ Algorithm

The fact that we only require u_h^n be in the desired range [c_1, c_2] at certain points in ∪_i,j S_ij can be used to reduce the computational cost in a great deal. 

Given the weighted L^2 projection u_h^0 computed from the initial data u_0(x, y), the algorithm is stated below:


  * Initialization.

Obtain u_h^0∈ V_h using the standard piecewise L^2 projection 

    ∫_I_ij u_0(x, y)ϕ dx dy =∫_I_ij u_h^0ϕ dx dy    for  ϕ∈ V_h.



  * Time evolution.

For n = 0, 1, 2, ...,


  (a) Check the point values of u_h^n on the test set S_ij. 
If one of them goes outside of [c_1, c_2], reconstruct û_h^n using the formula (<ref>) and (<ref>) and set u_h^n = û_h^n.


  (b) Use the scheme (<ref>) to compute u_h^n+1.

End


This algorithm is guaranteed to produce numerical solutions within the range with uniform third order accuracy for smooth exact solutions.

The algorithm with  forward Euler time discretization can be extended to high order ODE solves, such as the strong stability preserving Runge-Kutta methods, since they are a convex linear combination of the forward Euler; see <cit.>. 
The desired MPS property can be  ensured as long as the proper time step restrictions are respected.  



§ NUMERICAL TESTS
 
In this section, we present the results of numerical tests using our third-order maximum-principle-preserving DDG schemes. 
The numerical integration is computed using the Gaussian quadrature rule. Since M(x) and A(x) can be complex functions, sufficient number of quadrature points will guarantee the desired order of accuracy and induce small numerical errors.  Hence we take 16 quadrature points in each cell through all the examples.
As for the time stepping, we implement the SSP(3,3) scheme as in <cit.> for the strong-stability-preservation in time.
The one-dimensional error is measured by the discrete norms:

    e^h_p(t) =(∑_j=1^N_xu_h(t,·) - u(t,·)^p_L^p(I_j))^1/p,  for  p = 1  or  2,
 
and u(t,x) is taken as the exact solution or the reference solution given by greatly refined spacial discretization. 
If neither of them is available, we can compute the consecutive errors between u_h(t,x) and u_h/2(t,x), where the subindex indicates the mesh size h and h/2.
We introduce e^h_∞(t) to demonstrate the discrete MPS property at the test points in S_j: 

    e^h_∞(t) = max_x ∈ S_j, 0 ≤ j ≤ N{c_1 - u_h(t, x), u_h(t, x) - c_2},

If e^h_∞(t)>0, then the discrete MPS property is violated.



 §.§ One dimensional numerical tests
 
In our numerical tests we choose scheme parameters as 

    β_0 = 2,  β_1 = 0.16,  γ = 0.1.



 §.§ Accuracy test

We construct a linear problem of form (<ref>) to demonstrate the third order accuracy of our numerical schemes:

    {[ M(x)∂_t u = ∂_x(A(x) ∂_xu),          x ∈ [1, 3],   t>0,;        u(0, x)=sin (x^2-1),                  x ∈ [1, 3] ].

with M(x) = 4xe^-x^2+1, A(x) = e^-x^2+1/x. The exact slution is given by 

    u(t,x) = exp(-t)sin(x^2-1-t).

The boundary condition is imposed by using the exact solution.
We take the final time t=0.1 with different mesh sizes N_x = 16, 32, 64 and 128. 
Fig. <ref> shows the logarithm of the L^2 error for u_h, denoted by circles.
One can observe the third order of accuracy for u_h.





 §.§ Porous medium equation
 The porous medium equation

    ∂_t u =∂_x^2 (u^m),    m>1

is known to admit the Barenblatt solution of the form 

    B_m(t, x)=t^-α[ 1-α(m-1)/2m|x|^2/t^2α]_+^1/m-1,    with α = 1/m+1,

which is compactly supported.  

Fig. <ref> provides the accuracy test on (<ref>) with m=5 at t = 0.1.
Removing the two non-smooth corners, one can observe the third-order accuracy for the numerical solutions of u_h.

We compute the numerical solution with initial data B_m(1,x) subject to zero boundary conditions for m=2, up to final time t=3. 
From the numerical results in Fig. <ref>(a) with the MPS limiter we see a sharp resolution of discontinuities, and keeping the solution strictly within the initial bounds everywhere for all time.  
Fig. <ref>(b) shows a zoom-in at the nonsmooth corner for x ∈ [-6, -4] where the solution u is well simulated.
In contrast, without MPS limiter, it brings in significant overshoots near the upper bound of the exact solution, as evidenced by oscillations already appearing at t=1.0025 in Fig. <ref>(c).  In addition, the scheme without the MPS limiter will blow up in a short time.



 §.§ The Buckley-Leverett equation
 The convection-diffusion Buckley-Leverett equation of the form 

    ∂_t u +∂_x f(u) =ε∂_x (ν(u) ∂_x u),

is a model often used in reservoir simulations; see <cit.>. Here ε>0 is a small parameter,  f has an s-shape:

    f(u)=u^2/u^2+(1-u)^2,

and 

    ν(u)=4 u(1-u)1_0≤ u≤ 1.

We numerically solve (<ref>) with ε=0.01, subject to the following initial and boundary conditions 

    u(0, x)=(1-3x)1_0≤ x≤ 1/3,    u(t, 0)=1,    u(t, 1)=0.

The exact solution is not available.  With numerical convergence we demonstrate that our numerical scheme is capable of simulating the sharp corner of the solution that is moving in time. From the results in Fig. <ref> we observe the numerical convergence when the spacial mesh is refined. 
Moreover,  the lower bound  of u_h(t,x) is well preserved around the corner point x = 0.5.
We note that the numerical solution here is comparable to that obtained in <cit.> by the second order central scheme. 





 §.§ Two dimensional numerical tests

We consider the convection-diffusion equation:

    ∂_t u + ∇· f(u) = ∇·(A∇ u)  for  (x, y) ∈ [-1, 1]×[-1,1]  and  t>0,

where

    f(u) = uv⃗,   v⃗:=(0.01,0.01)^⊤

and A is a symmetric, positive definite matrix. 
For such an equation, exact solutions can be found of the form u(t,x,y)=a(t)exp(-ξ^⊤ B(t) ξ) with ξ=x-v⃗t, x:= (x,y)^⊤, provided a'=-a tr (AB) and B'=-2B^⊤ AB. 
Here a(0) can be chosen small enough to ensure that the periodic boundary condition adopted is reasonable and accurate at a finite time.  
In fact, if we set σ_0 = 0.01 and a 2×2 matrix σ(t) be such that

    σ(t) = σ_0^2Id + 2At    with Id being the identity matrix.

Then the function

    u(t,x,y) = σ_0^2/|(σ)|^1/2exp(-(x-v⃗t)^Tσ^-1(x-v⃗t)/2),   x:= (x,y)^⊤

 as given in <cit.>, is a solution to equation (<ref>).

In our numerical tests,  we take three choices of the tensor A as 

    A = 
    ([ 1 0; 0 1 ]),
    ([ 1 0; 0 2 ])
     or ([ 1 1; 1 2 ]),

which are usually denoted as isotropic, diagonally anisotropic and fully anisotropic diffusion tensors.  

We begin to demonstrate the necessity of the MPS limiter using the first isotropic problem.
Table <ref> shows the minumum and maximum of the numerical solutions u_h(t,x,y) using the MPS limiter over all the points in the test sets S_ij.
They are well bounded in the interval [0,1], satisfying the maximum principle.
Table <ref> shows the error e_∞^h(t) in (<ref>) without the MPS limiter.
One can observe that the numerical solutions u_h(t,x,y) violates the maximum principle and the simulation will break down after some time.
Larger β_0 helps to suppress the overshoot or undershoot, but cannot realize the bound preservation ideally.


In the first two cases, we take β_0 = 2, β_1 = 0.16, γ = 0.1 for both variables x and y. 
For the fully anisotropic one, we take β_0 = 4. For smaller β_0,  small oscillations develop in time and the scheme can become unstable.
We observe a nice agreement of the numerical solution to the exact solution for all three cases of A. 





For the equation (<ref>) with the third choice of the matrix A, we take care of the anisotropy by using a larger β_0 = 4, and adpative mesh size and time steps.
In the beginning, the solution is highly concentrated around the origin and requires a good resolution.
Therefore, for t ∈ [0, 10^-5], we employ the discretization with Δ x = Δ y = 0.005, Δ t = 10^-7.
Afterwards, a coarser mesh with Δ x = Δ y = 0.01, Δ t = 10^-6 is used.
For time t = 0.01, Fig. <ref>(a-b) show the contours of the exact solution and the numerical solution with the MPS limiter.
Fig. <ref>(c) shows a slice of the exact solution (circles) and the numerical solution (crosses) for y = 0.
One can observe a good agreement of the two solutions.





§ CONCLUDING REMARKS

In this paper, we present third order accurate DDG schemes which can be proven maximum-principle-satisfying for a class of diffusion equations with  variable diffusivity in terms of spatial variables and/or the unknown, in both one and two dimensional settings. 
Through careful theoretical analysis and numerical tests, we show that under suitable CFL conditions, with a simple scaling limiter involving little additional computational cost, the numerical schemes satisfy the strict maximum principle while maintaining uniform third order accuracy. 
The methodology extends to three dimensional rectangular meshes as well.
The effectiveness of the maximum-principle-satisfying DG schemes has been demonstrated through extensive numerical examples.

The presence of variable diffusivity has led to slightly more restrictive CFL conditions in order to preserve the maximum principle. 
In two dimensional case, the CFL condition seems a main factor for the slow numerical time evolution.  It would be interesting to extend the present result to implicit schemes so to improve the computational efficiency. 



§ PROOF OF THEOREM <REF>

Upon regrouping, we decompose the right-hand side of (<ref>) as

    ⟨ u^n+1⟩_ij    =  μ_x/μ∑_σ=2^L-1ω̂^σH_1(ŷ_j^σ) + μ_y/μ∑_σ=2^L-1ω̂^σH_2(x̂_i^σ) 
          +μ_x ω̂_1/μ(H_1(y_j-1/2)+H_1(y_j+1/2) ) + μ_yω̂_1/μ(H_2(x_i-1/2)+H_2(x_i+1/2) ) + B.

Here ω̂_1 = ω̂_L =1/L(L-1) is used. From (<ref>) we see that 

    H_1(ŷ_j^σ)=R_i(M(·,ŷ_j^σ), μ, Δ x, a),    H_2(x̂_i^σ)=R_j(M(x̂_i^σ, ·), μ, Δ y, b)

for 1≤σ≤ L.
Notice that the terms in B, induced from the nontrivial c, involve only polynomials values at four vertices of K_ij, 
we proceed to regroup and combine them with H_1(y_j±1/2) and H_2(x_i±1/2), respectively, 
in the following  way: 

    B_1     = cτ/2Δ xΔ y(u_h(x_i-1/2^-, y_j-1/2^+) + u_h(x_i-1/2^+, y_j-1/2^+) - u_h(x_i+1/2^-, y_j-1/2^+) - u_h(x_i+1/2^+, y_j-1/2^+))   with  H_1(y_j-1/2),  
    
    B_2     = cτ/2Δ xΔ y(u_h(x_i+1/2^-, y_j+1/2^-) - u_h(x_i-1/2^+, y_j+1/2^-) + u_h(x_i+1/2^+, y_j+1/2^-) - u_h(x_i-1/2^-, y_j+1/2^-))  with  H_1(y_j+1/2), 
    
    B_3     = cτ/2Δ xΔ y(u_h(x_i-1/2^+, y_j-1/2^+) - u_h(x_i-1/2^-, y_j+1/2^-) + u_h(x_i-1/2^+, y_j-1/2^-) - u_h(x_i-1/2^+, y_j+1/2^+))  with   H_2(x_i-1/2), 
    
    B_4     = cτ/2Δ xΔ y(u_h(x_i+1/2^-, y_j-1/2^+) - u_h(x_i+1/2^-, y_j+1/2^-) + u_h(x_i+1/2^-, y_j-1/2^-) - u_h(x_i+1/2^-, y_j+1/2^+))  with   H_2(x_i+1/2).

We shall also use the following notations.  

    [   ω̃_i^x,1(γ^x,y) = ⟨γ^x-ξ(1+γ^x) + ξ^2⟩_i(y)/2(1+γ^x),   ω̃_j^y,1(γ^y,x) = ⟨γ^y-η(1+γ^y) + η^2⟩_j(x)/2(1+γ^y),;            ω̃_i^x,2(γ^x,y) = ⟨ 1 - ξ^2⟩_i(y)/1-(γ^x)^2,            ω̃_j^y,2(γ^y,x) = ⟨ 1 - η^2⟩_j(x)/1-(γ^y)^2,; ω̃_i^x,3(γ^x,y) = ⟨ -γ^x+ξ(1-γ^x) + ξ^2⟩_i(y)/2(1-γ^x), ω̃_j^y,3(γ^y,x) = ⟨ -γ^y+η(1-γ^y) + η^2⟩_j(x)/2(1-γ^y). ]

For the first group, we have 

    H_1(y_j-1/2) + μ/μ_xω̂_1B_1 = R_i(M(·, y_j-1/2), μ, Δ x, a) + μ/μ_xω̂_1B_1 
    
    =   [ω̃_i^x,1(γ^x, y_j-1/2) - μ a (α_3(-γ^x) +α_1(γ^x)) + cτμ/2Δ xΔ yμ_xω̂_1]u_h(x_i-1/2^+, y_j-1/2^+) 
       +[ω̃_i^x,2(γ^x, y_j-1/2) - μ a (α_2(-γ^x) +α_2(γ^x))]u_h(x_i^γ,y_j-1/2^+)
       +[ω̃_i^x,3(γ^x, y_j-1/2) - μ a (α_1(-γ^x) +α_3(γ^x)) - cτμ/2Δ xΔ yμ_xω̂_1]u_h(x_i+1/2^-, y_j-1/2^+) 
       +(μ aα_3(-γ^x) - cτμ/2Δ xΔ yμ_xω̂_1)u_h(x_i+1/2^+, y_j-1/2^+)
          +μ a[α_2(-γ^x)u_h(x_i+1^γ,y_j-1/2^+)+α_1(-γ^x)u_h(x_i+3/2^-, y_j-1/2^+) ]
       +μ a[α_1(γ^x) u_h(x_i-3/2^+, y_j-1/2^+)+α_2(γ^x)u_h(x_i-1^γ,y_j-1/2^+) ]
        + (μ aα_3(γ^x) + cτμ/2Δ xΔ yμ_xω̂_1)u_h(x_i-1/2^-, y_j-1/2^+),

and the second group reduces to  

    H_1(y_j+1/2) + μ/μ_xω̂_1B_2 = R_i(M(·, y_j+1/2), μ, Δ x, a) + μ/μ_xω̂_1B_1 
    
    =   [ω̃_i^x,1(γ^x, y_j+1/2) - μ a (α_3(-γ^x) +α_1(γ^x)) - cτμ/2Δ xΔ yμ_xω̂_1]u_h(x_i-1/2^+, y_j+1/2^-) 
       +[ω̃_i^x,2(γ^x, y_j+1/2) - μ a (α_2(-γ^x) +α_2(γ^x))]u_h(x_i^γ,y_j+1/2^-)
       +[ω̃_i^x,3(γ^x, y_j+1/2) - μ a (α_1(-γ^x) +α_3(γ^x)) + cτμ/2Δ xΔ yμ_xω̂_1]u_h(x_i+1/2^-, y_j+1/2^-) 
       +(μ aα_3(-γ^x) + cτμ/2Δ xΔ yμ_xω̂_1)u_h(x_i+1/2^+, y_j+1/2^-) 
        +μ a[α_2(-γ^x)u_h(x_i+1^γ,y_j+1/2^-)+α_1(-γ^x)u_h(x_i+3/2^-, y_j+1/2^-) ]
       +μ a[α_1(γ^x) u_h(x_i-3/2^+, y_j+1/2^-)+α_2(γ^x)u_h(x_i-1^γ,y_j+1/2^-) ]
        +
       (μ aα_3(γ^x) - cτμ/2Δ xΔ yμ_xω̂_1)u_h(x_i-1/2^-, y_j-1/2^+).

From the above two groups we see that all coefficients of solution values involved  are nonnegative if 


    aα_3(±γ^x) - Δ x|c|/2Δ yω̂_1≥ 0, 
       μ≤min{ω̃_i^x,1(±γ^x, y_j±1/2)/a(α_1(±γ^x) + α_3(∓γ^x)) + Δ x |c|/Δ yω̂_1,
    ω̃_i^x,2(γ^x, y_j±1/2)/2α_2(γ^x)}.


Here we used ω̃_i^x,3(γ^x, y) = ω̃_i^x,1(-γ^x, y) and α_2(γ^x) = α_2(-γ^x).

Similarly,  all coefficients of solution values in 

    H_2(x_i-1/2) + μ/μ_yω̂_1B_3 = R_j(M(x_i-1/2, ·), μ, Δ y, b) + μ/μ_yω̂_1B_3

and 

    H_2(x_i+1/2) + μ/μ_yω̂_1B_4 = R_j(M(x_i+1/2, ·), μ, Δ y, b) + μ/μ_yω̂_1B_4

are also nonnegative if


    bα_3(±γ^y) - Δ y|c|/2Δ xω̂_1≥ 0,
       μ≤min{ω̃_j^y,1(±γ^y, x_i±1/2)/b(α_1(±γ^y) + α_3(∓γ^y)) + Δ y |c|/Δ xω̂_1,
    ω̃_j^y,2(γ^y, x_i±1/2)/2α_2(γ^y)}.


Here we used ω̃_j^y,3(γ^y, x) = ω̃_j^y,1(-γ^y,x) and α_2(γ^y) = α_2(-γ^y).

Since ω̃_i^x,σ, ω̃_j^y,σ for σ = 1, 2, 3 only depend on M(x,y)|_K_ij and γ, therefore bounded  from below. We set such bound as 

    ω_ij = min_γ^x, γ^y, x∈Ŝ_i^x, y∈Ŝ_j^y{ω̃_i^x,1(±γ^x, y), 
    ω̃_i^x,2(γ^x,y), 
    ω̃_j^y,1(±γ^y, x), 
    ω̃_j^y,2(γ^y, x)}.

Notice also that for γ = max{|γ^x|, |γ^y|},  using β_1≤ 1/4, we have 

    α_1(±γ^x) + α_3(∓γ^x) = β_0 + 8β_1-2/1±γ^x≤β_0 + 8β_1-2/1-γ,

which also hods when γ^x is replaced by γ^y.  Hence both (<ref>b) and  (<ref>b) are ensured  by  (<ref>).

Observe that (<ref>a) and  (<ref>a) are implied by 

    β_0 +min_s∈ [-γ, γ]8β_1-3+s/2(1-s)≥|c|κ/2ω̂_1 min{a, b}.

Using β_1 ≤ 1/4 we see that  the minimum on the left hand side is -1, obtained at s=γ when γ = 8β_1-1, hence this relation gives the lower bound in (<ref>). 



§ PROOF OF THEOREM <REF> 
 
For simplicity of presentation, in the following we consider  A=A(x, y) instead of A(x, y, u), for which we only need to use {A_h} on interfaces, so that c=c(x, y), a=a(x, y) and b=b(x, y). 
 
Since u_h(x,y) is quadratic in terms of x and y respectively,  we can use the formula (<ref>) to obtain 

    ṗ(η)=   ω̇^1(η)p(-1) + ω̇^2(η)(γ) + ω̇^3(η)p(1) 
    
    =   2η -1-γ/2(1+γ)p(-1)+2η/(γ^2 -1)p(γ)+2η +1-γ/2(1-γ)p(1).

Here η is the variable in the reference element [-1,1], 
and ω̇^σ(η) = d/dηω^σ(η).
Then 

    ∫_J_j c(x,y){∂_y  u_h(x, y)}dy    = {u_h(x, y_j-1/2^+)}∫_-1^1 c(x,y_j+Δ y/2η)ω̇^1(η) dη
          +{u_h(x, y^γ_j)}∫_-1^1 c(x, y_j+Δ y/2η)ω̇^2(η) dη
          +{u_h(x, y_j+1/2^-)}∫_-1^1 c(x, y_j+Δ y/2η)ω̇^3(η) dη.

Here the average {·} is taken with respect to x^- and x^+.
Using the quadrature rule on the right-hand side of (<ref>), we obtain the following scheme

    u̅^n+1_ij= μ_x/μ∑_σ=1^L ω̂^σH_1(ŷ_j^σ) + μ_y/μ∑_σ=1^L ω̂^σH_2(x̂_i^σ)+B.

The test set now reduces to 

    S^x_i = {x_i-1/2, x_i^γ, x_i+1/2} = x_i + Δ x/2{-1, γ, 1},  
        S^y_j = {y_j-1/2, y_j^γ, y_j+1/2} = y_j + Δ y/2{-1, γ, 1}

with γ satisfying

    |γ| ≤1/3 and  |γ|≤ 8β_1-1.

It can be shown that there exists L ≥2k+3/2 and γ satisfying  (<ref>)  such that x_i^γ∈Ŝ^x_i and y_j^γ∈Ŝ^y_j, 
the corresponding quadrature weight is denoted by ω̂^*.
Upon regrouping,  the right-hand side of (<ref>) may be decomposed  as 

    u̅^n+1_ij
     =      μ_x/μ∑_σ≠ 1, *, Lω̂^σH_1(ŷ_j^σ) + μ_y/μ∑_σ≠1,*,Lω̂^σH_2(x̂_i^σ) 
        + μ_x/μ(ω̂^1 H_1(y_j-1/2)+ ω̂^* H_1(y_j^γ)+ ω̂^1 H_1(y_j+1/2) ) 
        + μ_y/μ(ω̂^1 H_2(x_i-1/2)+ ω̂^* H_2(x_i^γ)+ ω̂^1 H_2(x_i+1/2) ) + B.

Here ω̂^1 = ω̂^L is used. 
The terms in B will be regrouped correspondingly.  More precisely, the first term involving integration on J_j is regrouped in terms involving solution values at y_j-1/2^+, y_j^γ and y_j+1/2^-,  and combined with H_1(y_j-1/2), H_1(y_j^γ) and H_1(y_j+1/2), respectively.
We check upon the first group only: 

    H_1(y_j-1/2) + μ/μ_xω̂^1B_1 = R_i(M(·, y_j-1/2), μ, Δ x, a) + μ/μ_xω̂^1B_1 
    
    =   [ω^1 - μ a (α_3(-γ) +α_1(γ)) - μΔ x/2Δ yω̂^1∫_-1^1c(x_i-1/2, y_j+Δ y/2η)ω̇^1(η) dη]u_h(x_i-1/2^+, y_j-1/2^+) 
       +[ω^2 - μ a (α_2(-γ) +α_2(γ))]u_h(x_i^γ,y_j-1/2^+)
       +[ω^3 - μ a (α_1(-γ) +α_3(γ)) +μΔ x/2Δ yω̂^1∫_-1^1c(x_i+1/2, y_j+Δ y/2η)ω̇^1(η) dη]u_h(x_i+1/2^-, y_j-1/2^+) 
       +(μ aα_3(-γ) +μΔ x/2Δ yω̂^1∫_-1^1c(x_i+1/2, y_j+Δ y/2η)ω̇^1(η) dη)u_h(x_i+1/2^+, y_j-1/2^+) 
       +μ a[α_2(-γ)u_h(x_i+1^γ, y_j-1/2^+)+α_1(-γ)u_h(x_i+3/2^-,y_j-1/2^+) ] 
       +μ a[α_1(γ) u_h(x_i-3/2^+, y_j-1/2^+)+α_2(γ)u_h(x_i-1^γ, y_j-1/2^+) ] 
       +(μ aα_3(γ) -μΔ x/2Δ yω̂^1∫_-1^1c(x_i-1/2, y_j+Δ y/2η)ω̇^1(η) dη)
       u_h(x_i-1/2^-, y_j-1/2^+),

From the three groups of H_1(y_j-1/2), H_1(y_j^γ) and H_1(y_j+1/2), we see that all the coefficients of solutions values are nonnegative if 

    aα_3(±γ) - Δ x/2Δ y min{ω̂_1, ω̂^*}max_x ∈ S_i^x; σ|g_j^σ(x) | ≥ 0,
       μ≤min{min{ω^1,ω^3}/a(α_1(±γ) + α_3(∓γ)) + Δ x/2Δ y min{ω̂_1, ω̂^*}max_x ∈ S_i^x; σ|g_j^σ(x) |,
    ω^2/2aα_2(γ)},

where g_j^σ(x)=∫_-1^1 c(x, y_j+Δ y/2η)ω̇^σ(η) dη.  The rest of terms in B when combined with H_2(y) for y∈ S_j^y led to the following conditions 

    b α_3(±γ) - Δ y/2Δ x min{ω̂_1, ω̂^*}max_y ∈ S_j^y; σ|∫_-1^1 c(x_i+Δ x/2ξ, y)ω̇^σ(ξ) dξ| ≥ 0,
       μ≤min{min{ω^1,ω^3}/b(α_1(±γ) + α_3(∓γ)) + Δ y/2Δ x min{ω̂_1, ω̂^*}max_y ∈ S_j^y; σ|g_i^σ(y)|,
    ω^2/2aα_2(γ)},

where g_i^σ(y)=∫_-1^1 c(x_i+Δ x/2ξ, y))ω̇^σ(ξ) dξ. 
Note that both |∫_-1^1c(x, y_j+Δ y/2η)ω̇^σ(η) dη| and |∫_-1^1c(x_i+Δ x/2ξ, y))ω̇^σ(ξ) dξ| are bounded from above by 

    c_∞|∫_-1^1ω̇^σ(η) dη| ≤4c_∞/1-γ.

Using the fact that κ^-1≤Δ x/Δ y≤κ, ω^3 ≤ω^1, 1/L(L-1)=ω̂_1 < ω̂^*, 
and α_3(±γ) ≥β_0 - 1, as well as 

    α_1(±γ) + α_3(∓γ) ≤β_0 + 8β_1 - 2/1+γ.

We see that all the above needed constraints are implied by the more severe constraints, including  (<ref>) for β_0, 
and the CFL condition (<ref>).





§ ACKNOWLEDGMENTS
  Liu's research was partially supported by the National Science Foundation under Grant DMS1812666 and by NSF Grant RNMS (Ki-Net) 1107291. 






10

BR97
F. Bassi and  S. Rebay,
A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier–Stokes equations.
 J. Comput. Phys.,  131(2):  267–279, 1997. 

CHY16
Z. Chen, H. Huang and J. Yan.
Third order Maximum-principle-satisfying direct discontinuous Galerkin 
methods for time dependent convection diffusion equations on unstructured triangular meshes.
Journal of Computational Physics, 308:198–217, 2016.

CLZ17  
W.-X. Cao, H. Liu and Z.-M. Zhang. 
Superconvergence of the direct discontinuous Galerkin method for
convection-diffusion equations. 
Numer Methods Partial Differential Eq., 33:290–317, 2017.

CSG16
Z. Chai, B. Shi and Z. Guo.
A multip-relaxation-time Lattice Boltzmann model for general  nonlinear anisotropic convection–diffusion equations.
J Sci Comput.   69:355–390, 2016.

CS98 
B. Cockburn and C.-W. Shu.
 The local discontinuous Galerkin method for time-dependent convection–diffusion systems. 
 SIAM J. Numer. Anal.,  35(6):2440–2463, 1998. 

CYLLL16
J. Cheng, X. Yang, X. Liu, T.-G. Liu and H. Luo.
  A direct discontinuous Galerkin method for the compressible Navier-Stokes equations on arbitrary grids. 
 J. Comp. Phys.,  327: 484–502,  2016.

DY19 
J. Du and Y. Yang. 
  Maximum-principle-preserving third-order local discontinuous Galerkin method for convection-diffusion equations on overlapping meshes. 
 J. Comp. Phys.,  377: 117–141, 2019.

Fu73
H. Fujii.
Some remarks on finite element analysis of time-dependent field problems.
Theory and Practice in Finite Element Structural Analysis, University of Tokyo Press, Tokyo,  91–106, 1973. 

FH06
 I. Faragó and R. Horváth.
  Discrete maximum principle and adequate discretizations of linear parabolic problems.
 SIAM J. Sci. Comput., 28: 2313–2336, 2006.

FHK05 
I. Faragó, R. Horváth, and S. Korotov.
 Discrete maximum principle for linear parabolic problems solved on hybrid meshes.
Appl. Numer. Math., 53:249–264, 2005. 

FK12
 I. Faragó, J. Karátson and S. Korotov.
 Discrete maximum principles for nonlinear parabolic PDE systems.
IMA Journal of Numerical Analysis, 32(4):1541–1573, 2012. 
 
 GKS09
S. Gottlieb, D.I. Ketcheson and C.-W. Shu.
High order strong stability preserving time discretizations.
Journal of Scientific Computing, 38:251–289, 2009.

GKS11
S. Gottlieb, D.I. Ketcheson, C.-W. Shu.
Strong Stability Preserving RungeKutta and Multistep Time Discretizations.
World Scientific, 2011.

HW07
J. S. Hesthaven and T. Warburton.
Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and
  Applications.
Springer Publishing Company, Incorporated, 1st edition, 2007.

JX13 
Y. Jiang and Z. Xu.
Parametrized maximum principle preserving limiter for finite difference
WENO schemes solving convection-dominated diffusion equations.
SIAMJ. Sci. Comput., 35:A2524A2553, 2013.

KT00
A. Kurganov and E. Tadmor.
New High-Resolution Central Schemes for Nonlinear Conservation Laws and Convection-Diffusion Equations.
Journal of Computational Physics, 160(1):241–282, 2000.

Le92
R.J. LeVeque.
Numerical Methods for Conservation Laws.
Birkhäuser, Basel, 1992.

Liu15  
H. Liu. 
Optimal error estimates of the direct discontinuous Galerkin method for convection–diffusion equations.
Math. Comp.,  84:2263–2295, 2015.

LW16
H. Liu and Z.-M. Wang. 
An entropy satisfying discontinuous Galerkin method for nonlinear Fokker–Planck equations. 
J Sci Comput., 68:1217–1240, 2016.

LW17
 H. Liu and Z.-M. Wang. 
 A free energy satisfying discontinuous Galerkin method for Poisson–Nernst–Planck
systems. J. Comput. Phys., 238: 413–437, 2017.

LY09
H. Liu and J. Yan.
The Direct Discontinuous Galerkin (DDG) methods for diffusion
  problems.
SIAM Journal on Numerical Analysis, 47(1):675–698, 2009.

LY10
H. Liu and J. Yan.
The direct discontinuous Galerkin (DDG) method for diffusion with
  interface corrections.
Commun. Comput. Phys., 8(3):541–564, 2010.

LY12
H. Liu and H. Yu.
An entropy satisfying conservative method for the Fokker-Planck
  equation of FENE dumbbell model for polymers.
SIAM Journal on Numerical Analysis, 50(3):1207–1239, 2012.

LY14
H. Liu and H. Yu.
Maximum-principle-satisfying third order discontinuous Galerkin schemes for Fokker–Planck equations.
SIAM Journal on Scientific Computing, 36 (5):A2296–A2325, 2014.

LY15
H. Liu and H. Yu.
The entropy satisfying discontinuous Galerkin method for Fokker–Planck equations
Journal of Scientific Computing, 62(3):803–830, 2015. 

MH85
A. Mizukami and T. J. Hughes.
A Petrov-Galerkin finite element method for convection-dominated
flows: an accurate upwinding technique for satisfying the maximum principle.
Comput. Meth. Appl. Mech. Engrg., 50:181–193, 1985.

Ri08
B. Riviére. 
Discontinuous Galerkin methods for solving elliptic and parabolic equations. 
SIAM series: Frontiers in Applied Mathematics ,   2008.

Sh09 
C.-W. Shu. 
Discontinuous Galerkin methods: general approach and stability. 
In Numerical solutions of partial differential equations, Adv. Courses Math. CRM Barcelona, pages 149-201. 
Birkhauser, Basel, 2009.

SCS18
Z. Sun, J.A. Carrillo, and C.-W. Shu.
A discontinuous Galerkin method for nonlinear parabolic equations and gradient flow problems with interaction potentials.
 Journal of Computational Physics, 382:76–104, 2018.

SPZ18 
S. Srinivasana, J. Poggiea, and X. Zhang. 
A positivity-preserving high order discontinuous Galerkin scheme for convection–diffusion equations.
 Journal of Computational Physics, 366: 120–143, 2018.

TW08
V. Thomee and L.B. Wahlbin.
On the existence of maximum principles in parabolic finite element equations.
Mathematics of Computation,  77: 11–19,  2008.

VKH10 
T. Vejchodský, S. Korotov, and A. Hannukainen.
Discrete maximum principle for parabolic problems solved by prismatic finite elements.
Math. Comput. Simulation, 80: 1758–1770, 2010.

XQX15
T. Xiong, J.-M. Qiu and Z. Xu. 
High order maximum-principle-preserving discontinuous Galerkin method for convection–diffusion equations.
SIAM J. Sci. Comput.,  37(2):  A583–A608, 2015.

YXQX16 
P. Yang, T. Xiong, J.-M. Qiu and Z. Xu.
High order maximum principle preserving finite volume method for convection dominated problems.
 J. Sci. Comput.,  67(2):  795–820, 2016.


Zh17
X. Zhang.
On positivity-preserving high order discontinuous Galerkin schemes for compressible Navier–Stokes equations.
 J. Comput. Phys.,  328:301–343, 2017.

ZS10
X. Zhang and C.-W. Shu.
On maximum-principle-satisfying high order schemes for scalar conservation laws.
J. Comput. Phys., 229(9):3091–3120, 2010.

ZS11 
X. Zhang, C.-W. Shu.
Maximum-principle-satisfying and positivity-preserving high-order schemes for conservation laws: survey and new developments.
 in: Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, vol. 467, 
 The Royal Society, 2752–2776, 2011.

ZLS12
X. Zhang, Y. Liu, and C. Shu.
Maximum-principle-satisfying high order finite volume weighted
  essentially nonoscillatory schemes for convection-diffusion equations.
SIAM Journal on Scientific Computing, 34(2):A627–A658, 2012.

ZZS13
Y. Zhang, X. Zhang, and C.-W. Shu.
Maximum-principle-satisfying second order discontinuous Galerkin
  schemes for convection-diffusion equations on triangular meshes.
J. Comput. Phys., 234:295–316, 2013.



