A Self-Adaptive Numerical Method to Solve Convection-Dominated Diffusion Problems

Convection-dominated diffusion problems usually develop multiscaled solutions and adaptive mesh is popular to approach high resolution numerical solutions. Most adaptive mesh methods involve complex adaptive operations that not only increase algorithmic complexity but also may introduce numerical dissipation. Hence, it is motivated in this paper to develop an adaptive mesh method which is free from complex adaptive operations. The method is developed based on a range-discrete mesh, which is uniformly distributed in the value domain and has a desirable property of self-adaptivity in the spatial domain. To solve the time-dependent problem, movement of mesh points is tracked according to the governing equation, while their values are fixed. Adaptivity of the mesh points is automatically achieved during the course of solving the discretized equation. Moreover, a singular point resulting from a nonlinear diffusive term can be maintained by treating it as a special boundary condition. Serval numerical tests are performed. Residual errors are found to be independent of the magnitude of diffusive term. The proposed method can serve as a fast and accuracy tool for assessment of propagation of steep fronts in various flow problems.


Introduction
It is well known that convection-diffusion equations arise in a variety of important science and engineering fields, for example, thermodynamics, fluid mechanics, and chemical reactions [1,2].In many applications recognized as convectiondominated problems, diffusion may be quite small in contrast to convection, so that the solution will develop steep moving fronts that are nearly shocks.Efficient and accurate assessment of propagations of moving fronts is important.However, standard Finite Difference Methods (FDM) tend to introduce spurious oscillations and stabilized methods often resort to upwind schemes, which introduce excessive numerical dissipations that will oversmooth the fronts [3,4].Consequently, an extremely fine mesh is needed to increase the resolution, making the convection-diffusion equation a most challenging one to be solved numerically.
There are mainly two lines of research in numerical methods for achieving higher resolution with limited computation cost.One is the development of higher-order methods trying to reduce the numerical diffusion, such as high-order FDM with flux limiters [5][6][7][8][9][10] and the method of characteristics [3,[11][12][13][14].The other important approach is using adaptive meshes.That is, mesh points are concentrated in the region where the solution varies steeply and is less concentrated in slowly varying regions.
Different adaptive strategies, for example, Adaptive Mesh Refinement method and moving mesh methods, have shown success in solving convection-diffusion problems [15][16][17][18][19][20][21][22].As an adaptive scheme is designed to yield highly nonuniform mesh, discretization of the governing equation on general meshes should be settled first [23].As for how to locate mesh points adaptively, such methods often appeal to equidistributing a monitor function or solving mesh equations [19].As a whole, adaptive methods generally consist of two parts: one is to locate mesh points adaptively and the other is to solve the governing equation on the nonuniform mesh.As a result, computation complexity is certainly increased.Moreover, a local interpolation is usually needed to project the information to new meshes.Such a remeshing operation often brings numerical diffusion.So, it is desirable to develop a numerical method that solves the problem adaptively but without complex adaptive operations.
In this paper, we propose a self-adaptive numerical method to capture the steep fronts accurately.We first consider the following convection-diffusion equation to illustrate the main ideas of our method [1,2]: where the diffusive coefficient () is positive.Instead of a traditional spatial-discretized mesh, the proposed method is based on a "range-discrete mesh" which is uniformly distributed in the value domain [24].It is found that the rangediscrete mesh has an adaptive nature because the density of mesh points is proportionate to the spatial gradient of the solution.During the computation process, positions of the mesh points are determined by the governing equation, while their values are fixed.In this way, movement of the mesh points not only solves the equation but also locates mesh points adaptively.In other words, the proposed method combines the two parts of an adaptive method, the adaptive strategy and solving the equation, in one procedure.Therefore, neither complex adaptive operations nor any interpolations are involved in the proposed method, leading to a simple and efficient method for convection-diffusion problems.
It must be emphasized that, with a linear diffusive coefficient, the solution of ( 1) is smooth, though there exist steep gradients, which is the most concerned situation in the literature.However, with a nonlinear diffusive coefficient, though small, spatial gradient of the solution can be discontinuous or even infinite [25][26][27].This singularity phenomenon is common in many applications but often smeared out by traditional methods.By simply treating it as a boundary condition, such a singular point can be maintained in our method.
It should be noted that, for monotone solutions, the proposed method is similar to the method of Fayers and Sheldon [25].Unfortunately, as the unknown is no longer single valued in their Lagrangian form, their method does not work for nonmonotone solutions.The proposed method overcomes this disadvantage by discretizing the governing equation over dynamic control volumes.Furthermore, how to deal with different kinds of boundary conditions with the proposed method is discussed.This paper is organized as follows.Section 2 introduces the range-discrete strategy to obtain an adaptive mesh.Next, the numerical scheme based on the range-discrete mesh, including boundary conditions, will be illustrated in detail in Sections 3, 4, and 5. Serval numerical tests are conducted to validate the numerical method in Section 6 and conclusions are drawn in Section 7.
For a constant piece in the solution curve, the governing equation is satisfied trivially.The solutions on both sides of this constant piece do not interact with each other until the constant piece disappears.Consequently, each endpoint of the constant piece can be treated as a boundary condition for each side, respectively.With two adjacent mesh points assigned to the same value and marked as boundary points, for example,  1 and  2 shown in Figure 1, a constant piece is sufficiently described.
If  = () is nonmonotone, serval intersections of the line  =   and the curve  = () exist, for example, points  1 ,  2 ,  3 shown in Figure 1.It is intuitive to employ all these intersections, assigned with the same discrete value   , as different mesh points to describe the nonmonotone feature.Moreover, a mesh point should also be set at the position of an extremum of  = (), for example, point   in Figure 1, because, without it, the solution between  2 and  3 will be identified as a constant piece.
In contrast to spatial-discrete mesh used by FDM and many other methods, the discrete mesh obtained using the rules above is named as a "range-discrete" mesh.We should state that, though Δ  can be different from each other, we focus our attention on the equidistance situation for the sake of simplicity.
If there exit  mesh points in the spatial area [  ,   ], the density of mesh points can be defined as  = /|  −   |.For a range-discrete mesh, variation of the solution over the area [  ,   ] can be defined as where () [  ,  ] /|  −   | can be interpreted as the spatial variation speed of the solution.Thus, we can conclude from (2) that  is proportionate to the spatial variation speed of the solution, so that the density of mesh points is large where the solution varies quickly and vice versa, indicating an adaptive nature, which can also be seen from Figure 1.As a special case, the density of mesh points for a constant piece is actually zero, being quite computational effort saving.
In fact, as lim , we can observe an equidistribution of the solution gradient |  ()| in the range-discrete grid.This is an obvious adaptive criterion widely used in many adaptive mesh methods [19,28,29].While complex adaptive operations have to be involved in these methods, it is done by the range-discrete mesh simply via a equidistance discretization in value domain.
With the technique illustrated above, a range-discrete mesh can be obtained for general solutions.Next, the numerical scheme based on the range-discrete mesh for solving (1) will be present.

Numerical Scheme Based on the Range-Discrete Mesh
For the initial-boundary value problem of (1), one can get an initial range-discrete mesh using the technique illustrated in Section 2. Once a range-discrete mesh is obtained, a dynamic control volume can be formed around each mesh point.Take [ 1 (),  2 ()] in Figure 2 as an example, and let ( 1 (), ) ≡  1 = (  +  −1 )/2 and ( 2 (), ) ≡  2 = ( +1 +   )/2, making the control volume moving with the flow.Integrating (1) on this control volume and a small time interval where the integration of the convective and diffusive term can be done directly as As for the transient term, according to Appendix A, if the solution curve (, ) is monotone in [ 1 (),  2 ()], we have Consequently, the integration form of ( 1) is The item on the left-hand side of ( 6) is the variation of the total quantity of the control volume over the time interval [ 0 ,  1 ]; ( 1 )−( 1 ) 1 and ( 2 )−( 2 ) 2 can be interpreted as the left and right boundary flux, respectively.Thus, ( 6) is actually a natural result of mass conservation principle.It should be emphasized that ( 6) is exact as it is derived from (1) without any approximation.Using linear approximation in the FDM manner, we can obtain first-order discrete form of ( 6) as where If the items on the right-hand side of (7) are the values at , the discrete equation is explicit; and if they are values at  + Δ, it is implicit.One may notice that ( 7) is quite similar to the discrete equations of FDM and can be solved in the same way Figure 3: The dynamic control volume around an extremum point. of FDM, for example, Newton's method.Generally speaking, other kinds of discrete equations and higher-order forms can also be obtained in the same manner as FDM.
However, we should emphasize that the unknown in ( 7) is   , not   , which is the essential difference between the proposed method and FDM.It means that the movements of the mesh points are tracked while their values are fixed throughout the computation process.When new positions of the mesh points are located, the approximation of the solution is obtained.Moreover, adaptivity is obtained simultaneously as a benefit of the range-discrete strategy.
We should note that though derived in different manners, (7) is similar to that in the method of Fayers and Sheldon [25].And also, it only fits for monotone solutions in [ 1 (),  2 ()].
To make the range-discrete grid work for nonmonotone solutions, adjustments are needed.

Adjustments for the Extremum Point
As has been mentioned in Section 2, a mesh point, for example,   in Figure 1 or Figure 3, is set at the extremum if the solution is nonmonotone.For the control volume of such an extremum point, (6) does not fit and adjustments are needed.
The dynamic control volume [ 1 (),  2 ()] for an extremum mesh point   is shown in Figure 3, where  1 =  2 =  −1 + Δ/2.According to Appendix A, The items inside the brackets in (3) for the control volume of an extremum point are So, we have With an parabolic approximation of the solution curve between [ 1 (),  2 ()], we obtain Consequently, we have Thus, the discrete equation for an extremum point is obtained.
One can conclude from (10) that, with a positive (), lefthand side of (10) always decreases as time goes on.It means that the value of the extremum   is approaching  1 , making the extremum point a special case: its value and position are always changing and the control volume is shrinking.
To avoid too small a control volume, the mesh points  −1 and  +1 will be removed from the range-discrete mesh if |  − 1 | is less than a threshold value, for example, |  − 1 | < Δ/4.Then, let  1 =  −2 + Δ/2 and a new control volume is formed immediately to continue the computation.In this way, the extremum point in the initial mesh is always an extremum point and represents the movement of the extremum in the solution.

Dirichlet Boundary Condition.
A Dirichlet boundary condition, that is, (  , ) =   , can be used directly in the discrete equation in the same way as FDM.That is, one can consider a discrete point   , set at the boundary to provide the boundary information, with its position and value always fixed as   =   and   =   during the computing process.

Moving Boundary Condition.
As has been mentioned in Section 2, either endpoint of a constant piece in the solution is treated as a boundary condition.The value of this boundary is fixed, which is quite similar to the Dirichlet boundary condition, but the position of the boundary varies.Thus, we name it a moving boundary condition.
At a moving boundary, we have ((), ) =   and (()  )| () = 0 (see Appendix B).Thus, it can be involved in the range-discrete mesh by a boundary point   with   =   and (()  )|   = 0.The boundary condition is not involved in a discrete manner as what was done in [25], but in a more precise manner, making it more promising for high precision.
One can see that the position of the moving boundary point is not really needed in the range-discrete equation.
However, a constant piece in the solution might disappear and the solution on two sides will then interact with each other.To illustrate this, we need to merge the two boundary points into a normal mesh point if their distance is less than a threshold  min .The position of the boundary point is only needed here and can be interpolated using () =  0 ( − )  by two nearby mesh points, where  can be settled using the method illustrated in Appendix B.
As a benefit of such a special boundary condition, the singular point, where the gradient of the solution might be discontinuous or infinite, as indicated in Appendix B, can be maintained by the proposed method.Numerical tests in Section 6 also confirm this.

Neumann Boundary Condition.
Given a Neumann boundary condition   (  ) =   , a boundary mesh point   will also be set.The discrete value of   is determined dynamically both by the value of   , which is the point nearest the boundary, and by the sign of   , that is,   =   − sgn(  )Δ for the left boundary and   =   + sgn(  )Δ for the right boundary.The position of   can be calculated by an linear interpolation: or other higher-order interpolations.As we can see, the position of a Neumann boundary point determined by ( 13) can be either in or out of the calculation area.So, we design the following rules: (1) If   is out the calculation area, computing continuous on.
(2) If   is in the calculation area, this boundary point becomes a normal mesh point and a new boundary point will be inserted into the range-discrete mesh.
(3) If   is out of the calculation area,   will be removed immediately and the information of   is determined by the new nearest boundary mesh point.

A Travelling Wave Solution.
With () =  2 /2 and () being a constant, we have Burgers' equation from (1): It is easy to verify that this equation has the travelling wave solution: where  is the speed of the wave.With  = 1, (15) is the solution of ( 14) with an initial condition (, 0) = 2 −/ /(1 +  −/ ). Figure 4 and Table 1 show the solutions and error tables at different time using the proposed method when  = 10 −1 , 10 −2 , 10 −3 .It can be seen from Figure 4 that mesh points are located adaptively according to the spatial gradient.Moreover, the adaptivity follows the movement of the travelling front.As a benefit, with only a few mesh points, the front is solved with good accuracy.One should notice that, with the diffusive coefficient  decreasing, the front becomes steeper, tending to form a moving shock.While it is hard to maintain the accuracy near the front with the same number of mesh points using FDM, the accuracy of the proposed method for different diffusive coefficients is uniform.

A Nonmonotone Case.
In this numerical example, the solution of Burgers' equation is designed to be nonmonotone to test the numerical method.The initial condition and boundary condition are set to be (, 0) = sin(), (0, ) = 0 and (1, ) = 0, respectively.According to Caldwell and Smith [30], the solution is where The numerical solutions at different times with diffusive coefficients  = 10 −1 , 10 −2 , 10 −3 are given by Figure 5 and the errors are Table 2.The results also show good performance of the proposed method for nonmonotone solutions.Again, the adaptivity functions well and tracks the front in an accurate manner.

Nonlinear Diffusive Term.
How a nonlinear diffusive term effects the distribution of the front is tested here.Consider Burgers' equation with a nonlinear diffusive term: with two Dirichlet boundary conditions (0, ) = 1, (1, ) = 0 and a lamp initial condition, which is (0 ≤  ≤ 0.2,  = 0) = 1 − 5 and (0.2 <  ≤ 1,  = 0) = 0. Diffusive coefficients set as () = 0.1  with  = 0.5, 1, 2 are tested.Solutions at  = 0.5 are shown in Figure 6.It is shown by Figure 6 that the gradient of the solution with  = 0.5 tends to be zero at the front and thus a smooth solution.However, the gradient with  = 2 tends to be infinite at the front, resulting in a continuous but not     phase, for example, water, and the region  > 0 is occupied by a nonwetting phase, for example, oil, initially.The governing equation can be obtained by Buckley-Leverett theory [31]: where  is the saturation of the wetting phase, () is the total infiltration rate,   () represents the capillary pressure, and () =     /(    +     ) is the flux coefficient of the wetting phase.In this example, () = 1 and porosity and absolute permeability of the layer are nondimensionalized to be  = 1 and  = 1, respectively.Relative permeability of the wetting phase (w) and the nonwetting phase (n) are   =  4  and   = (1 + )(1 − ) 3 , respectively, and their viscosity is   =   = 1.The capillary pressure is the -function   = 0.1 −0.5 .The convective term is typically nonconvex.Another problem is that the diffusive term is nonlinear, which will cause infinite saturation gradients at the front of each phase.Most methods will smear this singularity out because of their numerical diffusion, no matter how small.With the rangediscrete mesh method, this singular feature can be maintained by dealing the front as a moving boundary condition.
The numerical solutions at different time with Δ = 0.05 are shown in Figure 8.It is obvious that there exists not only a steep front where saturation varies quickly but also two singular points at  → 0 and  → 1, where the saturation gradient tends to be infinite.As a benefit of the range-discrete mesh and the moving boundary condition, the steep front is described preciously and two singular points are maintained intactly with only 20 mesh points.

Conclusion
Achievements of this paper are summarized as follows: (i) Traditional adaptive mesh methods involve complex adaptive operations and introduce numerical dissipation by interpolation.
(ii) We have developed a self-adaptive numerical method based on uniformly distributed range-discrete mesh for convection-diffusion equation.The governing equation is discretized by integration on the dynamic control volume of each mesh point.The goal of adaptivity and solving the equation are achieved simultaneously by solving the discretized equation, which moves the mesh points but does not change their values.
(iii) The method was verified against the analytical solution of serval numerical examples.Steep fronts were tracked accurately with only a few mesh points.The residual error is found to be independent of the magnitude of diffusion.
(iv) With a moving boundary condition, singular points at the fronts were maintained by the proposed method.It is found that, despite its magnitude, a nonlinear diffusion coefficient () =   ( ≥ 1) can produce singularities where the gradient of the solution can be discontinuous or infinite.
In all, the proposed method has shown advantages in solving convection-diffusion equations with high efficiency.The method can serve as a fast numerical tool for a wide range of front propagation problems, such as underground water flow, chemical floods of oil fields, and air pollution problems.On the other hand, it is obvious that, for higher dimensions, the elements in the range-discrete mesh become contour lines or contour planes and the solution can be achieved by tracking their movements.This work is in our next research scope.

A. Integration of the Transient Term in (3)
In (3), we have the integration of the transient term over the dynamic control volume [ 1 (),   So, integration of the transient term in (3) is The items inside the brackets in (A.2) are In short, the integration of the transient item is ( (,  0 ) −  2 ) d, for non-monotone  (, ) . (A.7)

B. The Moving Boundary Condition
Assume that there is a piece of constant solution at time , a moving boundary condition is formed at either end of the constant piece.In this section, without loss of generality, the constant piece is  > (), which means that ( ≥ (), ) =   , and  = (t) is the so-called moving boundary.
A dynamic control volume [ 1 ,  2 ], where  1 () = () and  2 is an arbitrary point in the constant piece, is chosen as shown in Figure 10 However, () is not a given condition but rather a result of the solution.It can be calculated by an interpolation as one of the two following cases: It can be seen from (B.5) that when  > 1, the gradient of the solution tends to be infinite when  → ().The situation when  = 2 graphed in Figure 7 shows us an example.
With the information of the two neighboring mesh points known, () can be obtained approximately using (B.5) or (B.8).

Figure 1 :
Figure 1: Description of the range-discrete mesh.

Figure 2 :
Figure 2: A dynamic control volume around a mesh point.

Figure 4 :
Figure 4: Numerical solutions at  = 0.5 and  = 1 of the travelling wave with different diffusive coefficients.Solid line represents exact solutions and squares, circles, and triangles in - figures are numerical solutions.The - figures below each - figure show the spatial distribution of the mesh points at different time.

Figure 5 :Figure 8 :
Figure 5: Numerical solutions obtained by the proposed method at  = 0.5 and  = 1 of Burgers' equation with different diffusive coefficients.Solid line represents exact solutions and squares, circles, and triangles in - figures are numerical solutions.The - figures below each - figure show the spatial distribution of mesh points at different time.

Figure 9 :
Figure 9: Geometric description of the integration.

Figure 10 :
Figure 10: The description of a moving boundary condition.

Table 1 :
Error table of the numerical solution of Burgers' equation with different diffusive coefficients at different time. is the number of mesh points.

Table 2 :
Error table of the numerical solution of the burgers' equation with different diffusive coefficients at different time. is the number of mesh points.