A high-order finite volume remapping scheme for nonuniform grids: The piecewise quartic method (PQM)

doi:10.1016/j.jcp.2008.04.026

Journal of Computational Physics

Volume 227, Issue 15, 20 July 2008, Pages 7394-7422

https://doi.org/10.1016/j.jcp.2008.04.026 Get rights and content

Abstract

A hierarchy of one-dimensional high-order remapping schemes is presented and their performance with respect to accuracy and convergence rate investigated. The schemes are also compared based on remapping experiments in closed domains. The piecewise quartic method (PQM) is presented, based on fifth-order accurate piecewise polynomials, and is motivated by the need to significantly improve hybrid coordinate systems of ocean climate models, which require the remapping to be conservative, monotonic and highly accurate. A limiter for this scheme is fully described that never decreases the polynomial degree, except at the location of extrema. We assess the use of high-order explicit and implicit (i.e., compact) estimates for the edge values and slopes needed to build the piecewise polynomials in both piecewise parabolic method (PPM) and PQM. It is shown that all limited PQM schemes perform significantly better than limited PPM schemes and that PQM schemes are much more cost-effective.

Introduction

Remapping is a crucial component of most arbitrary Lagrangian–Eulerian (ALE) algorithms used in computational fluid dynamics [13], [17], [14]. These algorithms involve a regridding step, whereby a new grid is generated based on some criteria, and a remapping step, whereby the variables are remapped from the old grid onto the new grid (Fig. 1). It is generally required that remapping be both conservative and monotonic in the sense that no new extrema should be created nor existing ones amplified. This is particularly important in applications where boundedness of some variables must be guaranteed.

The present study is motivated by the growing need to improve vertical coordinate systems in ocean general circulation models used for climate predictions. Over the last four decades or so, the vast majority of ocean models have used a single coordinate system in the vertical, usually aiming at a better representation of selected physical processes. These ocean models, however, have difficulties resolving physical processes for which they were not primarily designed. Hybrid coordinate ocean models have thus naturally emerged where the vertical grid is built by combining different coordinate systems in different regions [16], [2], [6], [1]. Due to the dynamical nature of the ocean, these hybrid coordinate systems are adapted in the course of the simulations, which means that vertical remapping is a key component. The accuracy of remapping is a major research issue in hybrid coordinate ocean models that prevents the hybrid framework from being truly convincing; presently, third-order reconstruction is used, at best, in these models. Hence, there is a need to explore higher-order representations.

An essential element of any remapping scheme lies in the reconstruction. Based on a set of cell averages on a given grid, the objective of a reconstruction scheme is to accurately and conservatively represent the underlying, real data with piecewise functions. Other constraints, such as monotonicity, may also apply. In this paper, we limit ourselves to piecewise polynomials but other functions can be employed, such as rational functions, e.g. [20]. Once the polynomial degree is chosen (for example, degree two leading to piecewise parabolas), a number of degrees of freedom must be determined that will produce a unique polynomial over the cell. There are basically two ways of improving a polynomial reconstruction method: First, the accuracy of the estimates for the pending degrees of freedom (e.g., edge values) may be improved. Second, when a monotonic reconstruction is required, improving the limiter may yield more accurate results. The latter domain of improvement, in particular, has drawn most of the attention in the case of the piecewise parabolic method (PPM) [4], [3], [12], [18], [11]. Besides these efforts, piecewise polynomial methods have witnessed almost no improvement regarding the edge estimates or the use of higher-degree polynomials. Only recently has the parabolic spline method been introduced but it also relies on polynomials of degree two [21]. Daru and Tenaud [5] clearly showed that successively high-order schemes lead to more accurate solutions for both limited and unlimited finite volume fluxes. The close relationship between their third-order method and PPM inspired us to seek improvement by using higher-order polynomials.

The objective of this paper is twofold: First, we introduce the piecewise quartic method (PQM) that uses piecewise polynomials of degree four and which, to our knowledge, has not been presented before. A limiter is devised that ensures monotonicity of any PQM-based remapping scheme. Second, a range of explicit and implicit schemes to estimate the edge values and slopes is investigated on the basis of accuracy and convergence analysis.

The main part (Section 2) of the paper focuses on reconstruction. We present PQM and the wide range of explicit and implicit schemes to estimate the edge values and slopes. The associated PPM and PQM schemes are evaluated in terms of accuracy. A limiter for PQM is also described in detail. The treatment of boundaries is covered in Section 3 and some comments on computational costs are made in Section 4. Conclusions are given in Section 5.

Section snippets

Piecewise polynomial reconstruction schemes

Given a nonuniform grid of cell widths $h_{j}$ and cell averages ${\bar{u}}_{j}$ for $j = 1, 2, \dots, N$ , where N is the number of cells, the objective is to determine a piecewise polynomial reconstruction that accurately approximates the underlying data (Fig. 2). The polynomial over cell j is noted $R_{j}$ and, for convenience, use is made of a local coordinate $ξ \in [0, 1]$ such that the global coordinate x is given by $x = x_{j - \frac{1}{2}} + (x_{j + \frac{1}{2}} - x_{j - \frac{1}{2}}) ξ = x_{j - \frac{1}{2}} + h_{j} ξ,$ where the global coordinates of the left and right cell interfaces are $x_{j - \frac{1}{2}}$ and $x_{j}$

Treatment of boundaries

We now turn our attention on estimating the edge values and slopes at the boundaries. We limit our discussion to the following schemes: PPM $h_{4}$ , PPM ${ih}_{4}$ , PQM ${ih}_{4} / {ih}_{3}$ and PQM ${ih}_{6} / {ih}_{5}$ , because they turn out to be the most effective. The concepts presented here are readily applicable to other schemes.

For PPM $h_{4}$ , edge values are estimated by using fourth-order polynomials spanning four cells. Thus, this scheme is not directly applicable to the first and last two edges for which another approach must

Computational cost

The relative computational costs of PPM and PQM schemes are now briefly investigated on the basis of remapping experiments in closed domains, as described in Fig. 11(b) (high-order boundary conditions are used). The errors and elapsed computational times are reported in Table 5 for unlimited and limited remapping experiments consisting of 20,000 cycles between a uniform 100-cell grid and a nonuniform 90-cell grid. Both the unlimited and limited versions of PQM schemes are more cost-effective

Conclusion

We have presented a hierarchy of one-dimensional high-order remapping schemes and investigated their performance with respect to accuracy and convergence rate. The schemes have also been compared based on remapping experiments in closed domains. We have introduced the new PQM scheme that is based on fifth-order accurate piecewise quartics. A limiter for this scheme has been fully described that never decreases the polynomial degree, except at the location of extrema where piecewise constants

Acknowledgments

We thank Bob Hallberg and Steve Griffies for comments regarding the topics covered in this paper. We also thank two anonymous reviewers for valuable suggestions that helped significantly improve the manuscript. Laurent White is supported by the ECCO2 project (http://ecco2.org/), “Estimating the Circulation and Climate of the Ocean, Phase II: High Resolution Global-Ocean and Sea-Ice Data Synthesis”, NASA award number NNG06GC28G. Laurent White is an honorary postdoctoral researcher with the

References (21)

C.N. Barron et al.
Formulation, implementation and examination of vertical coordinate choices in the Global Navy Coastal Ocean Model (NCOM)
Ocean Model
(2006)
R. Bleck
An oceanic general circulation model framed in hybrid isopycnic-Cartesian coordinates
Ocean Model
(2002)
P. Colella et al.
The piecewise parabolic method (ppm) for gas-dynamical simulations
J. Comput. Phys.
(1984)
V. Daru et al.
High order one-step monotonicity-preserving schemes for unsteady compressible flow calculations
J. Comput. Phys.
(2004)
G.R. Halliwell
Evaluation of vertical coordinate and vertical mixing algorithms in the HYbrid-Coordinate Ocean Model (HYCOM)
Ocean Model
(2004)
M.H. Kobayashi
Marcelo H. Kobayashion a class of Padé finite volume methods
J. Comput. Phys.
(1999)
C. Lacor et al.
A finite volume formulation of compact central schemes on arbitrary structured grids
J. Comput. Phys.
(2004)
L.G. Margolin et al.
Second-order sign-preserving conservative interpolation (remapping) on general grids
J. Comput. Phys.
(2003)
J. Morrell et al.
A cell by cell anisotropic adaptive mesh ALE scheme for the numerical solution of the Euler equations
J. Comput. Phys.
(2007)
M. Shashkov et al.
The repair paradigm and application to conservation laws
J. Comput. Phys.
(2004)

There are more references available in the full text version of this article.

Cited by (38)

A vertically-Lagrangian, non-hydrostatic, multilayer model for multiscale free-surface flows
2020, Journal of Computational Physics
Citation Excerpt :
In multiple dimensions this is not a simple operation, in particular when conservation of the remapped quantities is required [42]. In the one-dimensional case we are interested in, this is much simpler and can be done using for example Piecewise Polynomial Reconstruction [43]. In our implementation we have used the quadratic piecewise remapping of the PPR library of Engwirda and Kelley (2016, https://github.com/dengwirda/PPR).
This article presents a semi-discrete, multilayer set of equations describing the three-dimensional motion of an incompressible fluid bounded below by topography and above by a moving free-surface. This system is a consistent discretisation of the incompressible Euler equations, valid without assumptions on the slopes of the interfaces. Expressed as a set of conservation laws for each layer, the formulation has a clear physical interpretation and makes a seamless link between the hydrostatic Saint-Venant equations, dispersive Boussinesq-style models and the incompressible Euler equations. The associated numerical scheme, based on an approximate vertical projection and multigrid-accelerated column relaxations, provides accurate and efficient solutions for all regimes. The same model can thus be applied to study metre-scale waves, even beyond breaking, with results closely matching those obtained using small-scale Euler/Navier–Stokes models, and coastal or global scale dispersive waves, with an accuracy and efficiency comparable to extended Boussinesq wave models. The implementation is adaptive, parallel and open source as part of the Basilisk framework and the documented source codes sufficient to reproduce all results and figures are provided.
High-accurate and robust conservative remapping combining polynomial and hyperbolic tangent reconstructions
2020, Computers and Fluids
Citation Excerpt :
Such slope limiters may not always be trivial to design even if there exists a vast literature about slope/flux limiters, some of them are now known after their discoverers (van Leer [31–33], van Albada [34], Sweby [35], Barth-Jespersen [36], Venkatakrishnan [37], Koren [38,39], etc.) or their particularity (minmod [40], superbee [40], monotonized central [32], etc.). Relatively few attempts have been made to increase the remapping accuracy to more than second order via high order polynomial reconstructions, for several examples see [41–43]. Any attempt has to face the difficulty of designing a ‘limiter’ for non-smooth solution.
In this article we present a 1D single-material conservative remapping method that relies on high accurate reconstructions: polynomial ( $P_{4},$ $P_{1}$ with slope limiter) and non-linear hyperbolic tangent (THINC) representations. Such remapping procedure is intended to be used pairwise with a cell-centered Lagrangian scheme along with a rezone strategy to build a so-called indirect Arbitrary-Lagrangian-Eulerian scheme. Most of practically used Lagrangian schemes are second-order accurate. The goal of this work is to handle with accuracy contact using THINC reconstructions. At the same time, the smooth part of the solution is dealt with quartic polynomials, resulting locally in fifth order accurate remapping method. To ensure robustness, TVD-like reconstructions ( $P_{1}$ with slope limiter) are employed otherwise. A simple feature tracking algorithm is designed to assign a reconstruction type per cell ( $P_{4},$ $P_{1}^{lim}$ or THINC). This tracking algorithm is based on the nature of the contact waves which are traveling at the fluid velocity, while the shocks are compressive and detectable by following a change of cell volumes. Numerical results assess the behavior of such a remapping method on pure remapping problems of a scalar quantity and in the context of the full hydrodynamics equations. The associated indirect cell-centered ALE numerical scheme is run and produces numerical results that are presented to assess the extreme accuracy gained by such a remapping procedure employing a mix of reconstruction types.
Solution property preserving reconstruction BVD+MOOD scheme for compressible euler equations with source terms and detonations: Solution property preserving reconstruction BVD+MOOD scheme
2020, Computers and Fluids
In “Solution Property Reconstruction for Finite Volume scheme: a BVD+MOOD framework”, Int. J. Numer. Methods Fluids, 2020, we have designed a novel solution property preserving reconstruction, so-called multi-stage BVD-MOOD scheme. The scheme is able to maintain a high accuracy in smooth profile, eliminate the oscillations in the vicinity of discontinuity, capture sharply discontinuity and preserve some physical properties like the positivity of density and pressure for the Euler equations of compressible gas dynamics. In this paper, we present an extension of this approach for the compressible Euler equations supplemented with source terms (e.g., gravity, chemical reaction). One of the main challenges when simulating these models is the occurrence of negative density or pressure during the time evolution, which leads to a blow-up of the computation. General compressible Euler equations with different type of source terms are considered as models for physical situations such as detonation waves. Then, we illustrate the performance of the proposed scheme via a numerical test suite including genuinely demanding numerical tests. We observe that the present scheme is able to preserve the physical properties of the numerical solution still ensuring robustness and accuracy when and where appropriate.
The numerics of hydrostatic structured-grid coastal ocean models: State of the art and future perspectives
2018, Ocean Modelling
Citation Excerpt :
Originally Hirt et al. (1974) defined the Arbitrary Lagrangian-Eulerian (ALE) method by a Lagrangian mesh movement (LVD step) followed by an optional readjustment of the mesh associated with grid-related advective transports (EVD step). The advection during the EVD step can be replaced by instantaneous “remapping” (monotone and conservative interpolation, see e.g. White and Adcroft, 2008) which relaxes the strict CFL constraint to a less severe Lipschitz type criteria of the form Δt max |∇u| < 1. Gräwe et al. (2015) demonstrated the clear advantages of adaptive coordinates compared to σ-coordinates in a realistic state-of-the-art model application with different dynamical regimes (tidal/non-tidal parts, seasonal thermocline / permanent halocline).
The state of the art of the numerics of hydrostatic structured-grid coastal ocean models is reviewed here. First, some fundamental differences in the hydrodynamics of the coastal ocean, such as the large surface elevation variation compared to the mean water depth, are contrasted against large scale ocean dynamics. Then the hydrodynamic equations as they are used in coastal ocean models as well as in large scale ocean models are presented, including parameterisations for turbulent transports. As steps towards discretisation, coordinate transformations and spatial discretisations based on a finite-volume approach are discussed with focus on the specific requirements for coastal ocean models. As in large scale ocean models, splitting of internal and external modes is essential also for coastal ocean models, but specific care is needed when drying & flooding of intertidal flats is included. As one obvious characteristic of coastal ocean models, open boundaries occur and need to be treated in a way that correct model forcing from outside is transmitted to the model domain without reflecting waves from the inside. Here, also new developments in two-way nesting are presented. Single processes such as internal inertia-gravity waves, advection and turbulence closure models are discussed with focus on the coastal scales. Some overview on existing hydrostatic structured-grid coastal ocean models is given, including their extensions towards non-hydrostatic models. Finally, an outlook on future perspectives is made.
Attribution of horizontal and vertical contributions to spurious mixing in an Arbitrary Lagrangian–Eulerian ocean model
2017, Ocean Modelling
We examine the separate contributions to spurious mixing from horizontal and vertical processes in an ALE ocean model, MOM6, using reference potential energy (RPE). The RPE is a global diagnostic which changes only due to mixing between density classes. We extend this diagnostic to a sub-timestep timescale in order to individually separate contributions to spurious mixing through horizontal (tracer advection) and vertical (regridding/remapping) processes within the model. We both evaluate the overall spurious mixing in MOM6 against previously published output from other models (MOM5, MITGCM and MPAS-O), and investigate impacts on the components of spurious mixing in MOM6 across a suite of test cases: a lock exchange, internal wave propagation, and a baroclinically-unstable eddying channel.
The split RPE diagnostic demonstrates that the spurious mixing in a lock exchange test case is dominated by horizontal tracer advection, due to the spatial variability in the velocity field. In contrast, the vertical component of spurious mixing dominates in an internal waves test case. MOM6 performs well in this test case owing to its quasi-Lagrangian implementation of ALE. Finally, the effects of model resolution are examined in a baroclinic eddies test case. In particular, the vertical component of spurious mixing dominates as horizontal resolution increases, an important consideration as global models evolve towards higher horizontal resolutions.
High-order accurate finite-volume formulations for the pressure gradient force in layered ocean models
2017, Ocean Modelling
Discretisation of the horizontal pressure gradient force in layered ocean models is a challenging task, with non-trivial interactions between the thermodynamics of the fluid and the geometry of the layers often leading to numerical difficulties. We present two new finite-volume schemes for the pressure gradient operator designed to address these issues. In each case, the horizontal acceleration is computed as an integration of the contact pressure force that acts along the perimeter of an associated momentum control-volume. A pair of new schemes are developed by exploring different control-volume geometries. Non-linearities in the underlying equation-of-state definitions and thermodynamic profiles are treated using a high-order accurate numerical integration framework, designed to preserve hydrostatic balance in a non-linear manner. Numerical experiments show that the new methods achieve high levels of consistency, maintaining hydrostatic and thermobaric equilibrium in the presence of strongly-sloping layer geometries, non-linear equations-of-state and non-uniform vertical stratification profiles. These results suggest that the new pressure gradient formulations may be appropriate for general circulation models that employ hybrid vertical coordinates and/or terrain-following representations.

View all citing articles on Scopus

View full text

A high-order finite volume remapping scheme for nonuniform grids: The piecewise quartic method (PQM)

Abstract

Introduction

Section snippets

Piecewise polynomial reconstruction schemes

Treatment of boundaries

Computational cost

Conclusion

Acknowledgments

Ocean Model

Ocean Model

J. Comput. Phys.

J. Comput. Phys.

Ocean Model

J. Comput. Phys.

J. Comput. Phys.

J. Comput. Phys.

J. Comput. Phys.

J. Comput. Phys.