New approach to software design for moving mesh method

We intend to design our code in a ``plug and play'' fashion which is easy to use by the non-adaptive mesh user. We must emphasize here that the grid generated by the moving mesh is totally non-uniform; thus legacy codes written for the uniform grid, if they want to employ the advantages of moving mesh, should be rewritten for a nonuniform grid. This is in contrast to the adaptive mesh refinement (AMR) method, in which the uniformity can be preserved during the refinement and hence the uniform legacy codes can be used without any difficulties.

Our designs use a modular tree structure,

• Adaptive Mesh Method (AMM)
  • Static Rezone (SR)
    • Global Static Rezone (GSR)
      • find new grid
      • interpolation
      • warm restart
    • Adaptive Mesh Refinement (AMR), we will discuss this method in later chapters.
  • Moving Mesh (MM)
    • Time variation approach (TV)
    • Equidistribution principle (EP)
      • Dorfi and Drury method (DD)
      • Local smoothing method (LSM), MMPDE6 with two times Hamming smoothing
      • Other kinds of equidistributing mesh methods,
    • Regularization on spacing
    • Other kinds of moving mesh strategies
  • A simple h-r method
  • Mesh (or monitor) Function Evaluation(MF_EVA)
    • Arclength (ARC)
    • Curvature (CURVTR)
    • Arc-curvature, optimal combination of two
    • Smooth (SMTH)
    • Regularization (REGUL)
    • Other kinds of mesh functions
  • Initial Mesh Generator
    • Static rezone
    • Moving Mesh
  • Utility routines
    • Scaling
    • Smoothing
    • Interface to DASSL and DASPK

We now illustrate each part in more detail.

Mesh function evaluation

This is a component used by all of our adaptive mesh methods. For different adaptive mesh methods, the mesh function (or monitor function) might be the same or different. Generally, for equidistributing moving mesh methods, the monitor function can be arclength (most recommended) or curvature. For global static rezone methods based on the equidistribution principle and inverse interpolation, the curvature or arc-curvature monitors are recommended. For local mesh refinement or AMR, the mesh function can be much simpler (see Chapter 7 for detail). A smooth mesh function includes local smoothing or global smoothing. The regularization mesh function is only for a moving mesh to control the minimum spacing, which is most suitable for explicit moving mesh methods. The user can also supply other kinds of mesh function by supplying the ``other'' routine.

Static rezone

In this part, we provide two kinds of static rezone methods: global static rezone based on the equidistribution principle, and local mesh refinement method (adaptive mesh refinement (AMR)). The AMR method will be discussed in Chapter 7. For the global static rezone method, the initial number of nodes must be given. The user can choose a method with a fixed number of nodes, or a method of variable number of nodes, by giving an error-tolerance. There are two parts of the global static rezone. The first part is to find a new grid. Given the number of points in the new grid and the value of the mesh function for each cell of the old grid, the new grid can be generated via the dual-connecting inverse interpolation method discussed in [9]. The second part is to interpolate the solution from the old grid to the new one. The ``warm restart'' is also included in this part. The user can choose the ``warm restart'' for multistep time integrations or the ``full restart'' for single step time integrations. We supply the interfaces of the ``warm restart'' only to the DASSL and DASPK solver in our package. We have four interpolation routines: linear, cubic Hermite, monotone cubic Hermite, and quintic Hermite. For most smooth problems, cubic Hermite interpolation is sufficient. For problems with shocks or very steep wave fronts, the monotone cubic Hermite interpolation is recommended.

Moving mesh method

In this part, we supply several moving mesh strategies: the time variation approach, the Dorfi and Drury method, and the local smoothing method based on MMPDE6 and two Hamming smoothings. Regularizations for minimum and maximum spacing are also given (see Chapter 9 for details).

Initial mesh generator

Two kinds of initial mesh generator are provided. The first one computes the initial equidistributing mesh by the static rezone method. The second computes the initial equidistributing mesh by the moving mesh method. If the initial solution is given as a smooth function, the equidistributing mesh is computed only by the moving mesh equations. If the initial solution is given as discrete values, the equidistributing mesh is computed by the whole moving mesh system with u_t=0. The computation terminates when the mesh velocity is small enough (less than the error tolerance given by the user).

Structure of our adaptive mesh packages

Each module is a stand-alone routine and independent of the applications, and can be called without invoking the others. Expandability is achieved by the ``Other'' routine, which is supplied by the user. The user can also use a combination of two or three methods supplied by our packages, for example, ``0.5*ARC + 0.5*CURVTR + 2*SMTH'' means a mesh function of combination of arclength and curvature with two times smoothing. For moving mesh (MM) methods, the equation formulation is different for implicit and explicit integrations.

We have a friendly interface between the adaptive mesh modules with the other parts of the PDE solver. If a user has a code that works for a non-adaptive grid, and wants to incorporate the adaptive mesh modules into the code, the user does not need to recode but to set some parameters in the interface to indicate which adaptive strategies or which mesh function to use. The user can also turn off the adaptation at any time. If the user wants to test a different adaptive strategy or a different mesh function, he can supply the ``Other'' routine and does not need to code the adaptive part.

Our codes are originally designed for PDELIB in Los Alamos National Laboratory. So we use the PDELIB data structure and conventions to define a PDE solver. It is not difficult to adjust other PDE solvers to use our modules. Indeed, our adaptive modules are independent of any PDE solvers.

As in many other moving mesh codes, there are several parameters that should be specified by the user in the setting-up routine. Most of them are insensitive to different applications and can be chosen by default. The most critical parameter in equidistributing mesh moving strategies is the time scale $\tau$. A wrong choice of $\tau$ can result in starting problems or a grid lagging behind. If $\tau$ is very large and tends to infinity, the grid will not move and cannot be adapted. If $\tau$ is extremely small and tends to zero, the index of the DAE system may be high and there may be numerical problems. We have experienced that if the $\tau$ is changed, the number of time steps and the number of Jacobian evaluations are also changed. The optimal $\tau$ is also related to the size and smoothness of the monitor function. Many heuristic guidelines for how to choose $\tau$ have been proposed in the past decade, but none of them seems to work for all of the applications. This is especially true for explicit integration, where a small $\tau$ can result in mesh crossing and failures in integrations. For implicit integration methods, changing $\tau$ can result in a new evaluation of the Jacobian for the Newton iteration; hence the $\tau$should not change frequently and might be kept as a constant. In next section, we will give some heuristic suggestions on how to choose $\tau$ for general problems.