Automatic Mesh Refinement

The ability to dynamically alter the mesh during the solution of a numerical model has been a significant feature of all our finite element-based software since their release in the year 2000. We remain the only geotechnical software company worldwide to offer this advanced feature.

Since there is often a significant amount of confusion regarding AMR, this webpage is designed to answer each question and alleviate the confusion.

In particular, the following questions will be answered:

  1. Why is AMR important?
  2. What are the different mesh refinement methods?
  3. Are smearing errors ever encountered?
  4. Does the mesh actually move?
  5. In what situations is AMR beneficial?
    1. Penetration of sharp wetting fronts during infiltration in initially dry soil
    2. Sharp concentration fronts during convection-dominant solute transport problems
  6. Does each timestep need a different mesh?
  7. Should meshes be the same between applications?
  8. Why not just make the mesh dense everywhere?



1. Why is AMR important?

“Insufficient spatial resolution is commonly recognized to be the primary source of errors in numerical solutions of partial differential equations (PDEs) for water flow and solute transport in the vadose zone.” (Yeh, 2000)

Proper mesh assignment is critical to valid numerical solutions. But what is our traditional methodology? Typically the user i) creates outlining geometry and then ii) picks a preliminary grid that seems to suite model geometry. Most fixed-grid software companies will argue that when a pre-generated, fixed grid is well-chosen for a given problem, conventional numerical solution methods will provide valid results.

This approach is inherently flawed as noted by Tannehill, et al. 1997.

“One difficulty in solving a PDE with this approach (i.e., using a pre-generated, fixed grid) is that the grid is constructed and points are distributed in the physical domain before details of the solution are actually known. As a consequence, the grid may not be the best one for the particular problem.” (Tannehill, et al. 1997)

Other comments from literature note the following:

“Scientific and engineering computation has become so complex that traditional numerical computation on uniform meshes is generally not possible or too expensive” (Bern et al., 1999)

“Although faster computers have emerged in recent years, faster computers tend to be used to solve even more difficult problems.” (Finlayson, 1992)

“Therefore, developing efficient numerical methods that accurately solve a problem with a minimum of computer time continues to provide a worthy challenge. This is particularly important for three-dimensional simulations of large geographical areas of the vadose zone and associated groundwater zone.” (Mansell, 2002)

2. What are the different mesh refinement methods?

There are three main categories of LAGR methods

  • Mesh refinement (h-methods)
  • Moving mesh (r-methods)
  • Subspace enrichment (p-methods)

h-methods refine the mesh to increase nodal density in regions of the flow domain that have large errors. They can be very effective in producing near-optimal meshes for given error tolerances (Mansell, 2002). r-methods relocate nodes within a region to increase nodal densities in regions with large errors. r-methods have the downside in that they can generate smearing errors. p-methods increase the degree of interpolation in high-error regions.

Our software makes use of h-methods which do not have any associated smearing errors and are highly effective in producing near-optimal meshes. Error estimates are formed in each cell, and cells with high error are split into two new cells, resulting in a new computation mesh and a repeat of the solution process.  The exact algorithm is proprietary, but similar techniques are discussed in the open literature.

Reference: The importance of mesh resolution

3. Are smearing errors ever encountered?

No – for explanation see point #2.

4. Does the mesh actually move?

No – it is merely a subdivision and node-releasing process. Nodes can be added and taken away but never actually move in the mesh.

5. In what situations is AMR beneficial?

The benefit of AMR as pertaining to geotechnical numerical modeling problems is well-documented in open literature with respect to the following types of models:

a) Penetration of sharp wetting fronts during infiltration in initially dry soil,
b) Sharp concentration fronts during convection-dominant solute transport problems.

5.a) Penetration of sharp wetting fronts during infiltration in initially dry soil

Critical cases of water flow include evaporation near the soil surface and infiltration into initially dry soil profiles. These types of scenarios are common in numerical models of cover design and other climate-influenced models. Highly nonlinear relationships between hydraulic conductivity and pressure head contribute to very steep wetting fronts during infiltration into initially dry soil. In the vicinity of the wetting front for initially dry soil, small values of hydraulic conductivity require very large gradients to move even a small amount of water (Pan and Wierenga, 1995). Insufficient local resolution for such cases of water flow can result in numerical oscillation and numerical smearing (Mansell, 2002).

5.b) Sharp concentration fronts during convection-dominant solute transport problems.

Steep moving fronts of solute concentration may occur for convective-dispersive transport during water flow in porous media for conditions when the convection term dominates or the concentration moves along in a wave (Finlayson, 1992). Accurate numerical solutions for solute transport involving steep concentration fronts create special requirements for optimum resolution (Mansell, 2002). For convection-dominant transport of solutes, sharp peaks and valleys of concentration, as well as details of the solution controlling transport between nodes after each time step are truncated by interpolation between nodes. Such truncation is the source of all numerical problems in simulations of convective-dispersive solute transport (Yeh, 2000).

Truncation of peaks and valleys is an important cause of peak clipping, numerical spreading, and spurious oscillation Incorporation of local adaptive grid refinement algorithms in numerical models provides opportunity to enhance the accuracy of numerical approximations by automated adjustment of local spatial resolution for such cases (Mansell, 2002).

The benefit of AMR is of particular note in the solution of the classic “advection only” example model. In this scenario, a contaminant is released from the left side of the “pipe” and allowed to flow to the right only under the process of advection. A perfect solution is therefore a step function. Most software packages will not yield a step function, however and what appears to be diffusion is actually due to a truncation error. Standard finite difference method presents an non-ideal solution. The Method of Characteristics (MOC) and the Modified Method of Characteristics (MMOC) methods have subsequently been presented in literature along with the total variation diminishing (TVD) method for reducing the errors. The TVD produces reasonable results at the expense of not maintaining accurate mass balance in the formulation.

SVFlux was able to match the best results of the TVD method using the standard Galerkin finite element method along with AMR. In this case, the AMR was set such that it followed the moving contaminant front.

6. Does each time-step need a different mesh?

ChemFlux Automatic Mesh RefinementIt is entirely possible and reasonable to expect that each individual time-step in a transient analysis may need a separate mesh. This is particularly true in a infiltration or contaminant transport model where the wetting or contaminant front may be moving with time. If the front is moving with time then the zone of high gradients is constantly spatially changing and likely to require significant modifications to the mesh over time in order to accurately model the changing conditions.

The requirement for a separate mesh with each time-step is particularly apparent in the problem of a contaminant released from a box (or underground storage tank). In this model a contaminated “zone” is created at the top of the numerical model. The contaminant is then “released” at time=0 and the contaminant is allowed to flow down under both the processes of diffusion and a downward advective gradient.

It can be seen from the figure that the result is highly dependant on the mesh density used for the problem.

7. Should meshes be the same between applications?

Definitely not. The Galerkin finite element process solves based on an error estimate of the correctness of the solution relative to a primary variable. For example, head (h) is typically the formulated dominant variable in solution of the Richard’s based flow equation.

The reason it is not good to assume the same mesh for each type of analysis (i.e. flow, heat flow, air flow, stress/deformation, contaminant transport) is that the primary variable for which the finite element solver is solving for is completely different in each analysis. Using the same mesh for different analysis (as is done in other geotechnical software suites) may achieve a reasonable answer but is not a good practice and can lead to an unstable and non-ideal solution.

Let’s consider the following example of the modeling of a cross-section of an asphalt roadway. In this example there are three processes at work, namely; i) rain is falling on the shoulders, ii) the sun is coming out and warming the surface, and iii) traffic is moving down the roadway. Therefore these three processes may be modeled with seepage (SVFlux), thermal (SVHeat), and stress/deformation (SVSolid) software.

Let's take a look at the following example. In this case a road is modeled with a small stream to the left. In the stress analysis the user is interested in stress created by vehicle loading. Vehicle loading is applied to the asphalt road surface.



In the stress portion of the analysis we can see that mesh density is focused around the areas right under the tire loads. This is logical as these are the zones of high stress change.



In the seepage portion of the analysis, there is rainfall applied to all upper boundaries except for the road surface. It can be seen from the seepage analysis that mesh refinement is required on the upper boundary.




In the thermal aspect of the analysis there is heat generated from radiation on the asphalt surface. There is also a cooling effect from the stream next to the road. This thermal scenario again requires a seperate mesh for proper solution.

8. Why not just make the mesh dense everywhere?

A fixed grid with a fine mesh may provide accurate numerical solutions for certain problems.  High computational costs, however, may be prohibitive for complex multi-dimensional simulations of large-scale problems.

Speed is largely related to the number of nodes in the model. If the user proceeds with adding 5x more nodal density to the model then the solution time would be expected to grow by at least 5x.

“Scientific and engineering computation has become so complex that traditional numerical computation on uniform meshes is generally not possible or too expensive” (Bern et al., 1999)

“Although faster computers have emerged in recent years, faster computers tend to be used to solve even more difficult problems.” (Finlayson, 1992)

“Therefore, developing efficient numerical methods that accurately solve a problem with a minimum of computer time continues to provide a worthy challenge. This is particularly important for three-dimensional simulations of large geographical areas of the vadose zone and associated groundwater zone.” (Mansell, 2002).

In the example shown to the right, we take the standard problem of water infiltration into a dry soil, the model creates a wetting front in the soil which requires increased nodal accuracy in the zone of high gradient changes as may be seen in the adjacent figure.

As an experiment, the mesh for the example was then increased to make the mesh dense everywhere at the same level of density selected by the mesh refinement algorithm at the wetting front. If the mesh density was globally increased to this density, the model would then take more than 30 times longer to solve in order to achieve the same level of accuracy as an AMR-based solution.


REFERENCES

R. S. Mansell,* Liwang Ma, L. R. Ahuja, and S. A. Bloom - 2002
"Adaptive Grid Refinement in Numerical Models for Water Flow and Chemical
Transport in Soil: A Review"
Vadose Zone Journal, 1:222–238

Bern, M.W., J.E. Flaherty, and M. Luskin (ed.) 1999. Grid generation and adaptive algorithms. The Institute for Mathematics and its Applications. Vol. 113. Springer-Verlag, New York, NY.

Finlayson, B.A. 1992. Numerical methods for problems with moving fronts.
Ravenna Park Publ., Seattle, WA.

Pan, L., and P.J. Wierenga. 1995. A transformed pressure head-based approach to solve Richards’ equation for variably saturated soils. Water Resour. Res. 31:925–931.

Tannehill, J.C., D.A. Anderson, and R.H. Pletcher. 1997. Computational fluid mechanics and heat transfer. 2nd ed. Taylor and Francis Publishers, Washington, DC.

Yeh, G.-T. 2000. Computational subsurface hydrology: Reactions, transport, and fate. Kluwer Academic Publ., Boston, MA.