Remarks on the history of Arbitrary LagrangianEulerian finitedifference
Remarks on the history of Arbitrary Lagrangian-Eulerian finite-difference methods Rainer Bleck NASA Goddard Institute for Space Studies, NYC NOAA Earth System Research Lab, Boulder, CO October 2016
What is ALE? • ALE is a technique for blurring the distinction between fixed (“Eulerian”) and flow-following (“Lagrangian”) computational meshes used in solving partial differential equations. • “L” is generally understood to include adaptive (i. e. not necessarily Lagrangian) grids. • “A” stands for arbitrary.
Bibliography - general • Hirt, C. W. , A. A. Amsden, and J. L. Cook, 1974: An arbitrary Lagrangian-Eulerian computing method for all flow speeds. J. Comput. Phys. , 14, 227 -253. • Donea, J. , A. Huerta, J. -P. Ponthot, and A. Rodriguez. Ferran, 2004: Arbitrary Lagrangian-Eulerian methods. Encycl. Comput. Mech. , E. Stein et al. , Ed. , Wiley & Sons.
Bibliography -specific • Bleck, R. , and D. Boudra, 1981: Initial testing of a numerical ocean circulation model using a hybrid (quasi-isopycnic) vertical coordinate. J. Phys. Ocean. , 11, 755 -770. (Apologies for tooting my own horn. ) • Lin, S. -J. , 2004: A ‘‘vertically Lagrangian’’ finitevolume dynamical core for global models. Mon. Wea. Rev. , 132, 2293 -2307.
Hirt et al. (1974): J. G. TRULIO, Air Force Weapons Laboratory, Kirtland Air Force Base Report No. AFWL- TR-66 -19, 1966.
From Hirt et al. (1974)
Another 2 -D detonation example, this one from Donea et al. (2004): ALE Lagrangian
References to ALE literature in Donea et al. (2004):
Done with generalities. Focus now on ocean modeling… • Most ALE schemes were developed to solve 2 D and 3 -D initial value problems. • In geophysical fluid models, the horizontal mesh is typically fixed, so ALE is applied in the vertical only. • We use ALE in the vertical for a number of reasons…
First, some historical remarks • The idea to use a material vertical coordinate in geophysical fluid analysis and modeling goes back to -at least- the 1930 s (Shaw, Montgomery, Rossby). • In numerical weather prediction, its main attraction was the elimination of vertical transport terms (Starr, Shuman). • Arbitrarily chosen material surfaces were soon found to fail due to folding (Shuman). • But isentropic coordinates worked, presumably because their slope is tied to available potential energy (Eliassen). • Later, aligning coordinate surfaces with neutrally buoyant surfaces was proposed as an antidote to “false” (i. e. numerically induced) diabatic transport & diffusion. Today this is a major motivation for using generalized vertical coordinates.
The 3 -D grid in geophysical fluids • The distinguishing feature of geophysical fluid modeling is the pre-eminence of gravity. • Mixing of vertical and horizontal forces in the dynamic equations must therefore be avoided at all cost. • This is accomplished by making the z axis coincide with the direction of gravity. • The other 2 coordinate directions must be perpendicular to the z axis. This way, gravity will not contaminate the force balance there. • (ALE could, in principle, be applied in x, y direction, but grid columns would have to remain vertical at all times. )
The 3 -D grid in geophysical fluids (cont’d) • If ALE is applied in z direction, differential vertical motion of (Lagrangian) coordinate surfaces will inevitably destroy the “horizontal” alignment of neighboring grid points. • Once that happens, the horizontal force balance must be interpolated from information carried on sloping coordinate surfaces. • This interpolation, necessarily taking place in the direction of gravity, is the Achilles Heel of generalized-vertical coordinate modeling of geophysical fluids.
Donea et al. (2004) on the subject of grid maintenance: Both concepts are employed in the ALE ocean model HYCOM and its forerunners: regularization in the mixed layer, adaptation in the interior.
“adaptation” “regularization” ALE in early models of Bleck & Boudra (1981):
Continuity equation in generalized (“s”) coordinates (zero in fixed grids) (zero in material coord. ) (known)
Recapping… • Original concept (Hirt et al. , 1974): maintain Lagrangian character of coordinate but “re-grid” intermittently to keep grid points from fusing. • In HYCOM, we apply ALE in the vertical only and regrid for 2 reasons: – (1) to maintain minimum layer thickness (“regularization”); – (2) to nudge an entropy-related thermodynamic variable toward a prescribed layer-specific “target” value by importing water from above or below (“adaptation”). • (2) renders the grid quasi-isentropic. • (1) overrides (2).
Back to the “Achilles heel”: the PGF issue a) Aligning coordinate surfaces with isentropic surfaces avoids irreversible entropy production by numerically induced dispersion. b) Principal challenge: finding the horizontal pressure gradient force (PGF) on tilted coordinate surfaces: c) 2 -term PGF reduces to numerically stable 1 -term expression on surfaces of constant in-situ density:
The PGF issue (cont’d) d) In HYCOM we try to achieve a) by making potential density (rpot) the vertical coordinate. e) We try to achieve c) by introducing 2 mutually compatible approximations [Spiegel & Veronis, 1960]: (1) treat the fluid as incompressible (2) replace in-situ density r by rpot. f) This creates problems with thermobaricity [Sun et al. , 1999]. g) 2 nd problem: globally referenced rpot is not necessarily monotonic in the real ocean.
What Spiegel and Veronis (1960) say: The equations governing convection in a perfect gas are formally equivalent to those for an incompressible fluid if the static temperature gradient is replaced by its excess over the adiabatic … and the following approximations are valid: (a) the vertical dimension of the fluid is much less than any scale height, and (b) the motion-induced fluctuations in density and pressure do not exceed. . . the total static variations of these quantities.
The Adcroft et al. (1) approach • Goal: use realistic eqn. of state to overcome problems f) and g) in HYCOM. • Replace 2 -term finite-difference PGF expression by a line integral of “contact pressure” along the perimeter of each grid cell in the vertical plane – a finite-volume approach. • Doesn’t work well with traditional HYCOM grid generator (problem: gradually appearing thickness anomalies in middepth layers). • Works better with alternate grid generator originally developed for atmospheric applications (FIM(2)). Adcroft, A. , Hallberg, R. and M. Harrison, 2008: A finite volume discretization of the pressure gradient force using analytic integration. Ocean Modelling, 22, (3 -4), 106 -113. (2) http: //fim. noaa. gov (1)
Hycom 2 (=hycom with compressible eqn. of state) Things look good at the surface – for a while. However, looking below the surface …
Hycom 2 (=hycom with compressible eqn. of state) zonal velocity, yr 3. Ideal-age tracer, Yr 3.
Hycom 2 (=hycom with compressible eqn. of state) zonal velocity, yr 10. Ideal-age tracer, yr 10.
More successful attempts to use the full eqn of state are expected be presented in this meeting. Improving the grid “adaptation” algorithm likely is the key. Finally a news item from the numerical weather prediction community:
Lin’s (2004) atmospheric circulation model, which is on track to become the U. S. Weather Service’s “Next Generation Global Prediction System” (NGGPS), happens to be an ALE model:
unused slides
HYCOM limitations relevant for climate modeling • Aligning coordinate surfaces with isentropic surfaces avoids irreversible entropy production related to numerically induced dispersion. • Principal challenge: finding the horizontal pressure gradient force (PGF) on tilted coordinate surfaces. • 2 -term PGF reduces to a 1 -term expression (highly desirable!) on surfaces of constant in-situ density. • PGF is “almost” a 1 -term expression on surfaces of constant potential density. • Potential density comes much closer than in situ-density to being isentropic. Hence…. • Potential density is chosen as vertical coordinate. • Potential density can be multivalued in the vertical. • The pot. -density reference pressure should be chosen to avoid multivaluedness in critical regions (present choice: 10 KPa ~ 1 km).
Lagrangian versus Eulerian grids • Lagrangian methods can be capricious, and thus are generally avoided. They are good for solving certain initial value problems (shock formation, frontogenesis). • Extension to longer integration times (typically associated with boundary value problems) is problematic. • A generally accepted solution: periodic regridding => ALE (Arbitrary Lagrangian-Eulerian). • Applying ALE in a single dimension is relatively easy, but there are pitfalls. A bit of capriciousness remains.
The PGF issue a) Aligning coordinate surfaces with isentropic surfaces avoids irreversible entropy production by numerically induced dispersion. b) Principal challenge: finding the horizontal pressure gradient force (PGF) on tilted coordinate surfaces. c) 2 -term PGF reduces to numerically stable 1 -term expression on surfaces of constant in-situ density. d) In HYCOM we try to achieve a) by making potential density (rpot) the vertical coordinate. e) Giving the fluid in each coord. layer the density rpot and treating it as incompressible achieves goal c). f) This creates problems with thermobaricity. g) Potential density (globally referenced) can be multivalued in the vertical.
- Slides: 30