1 / 23

Numerical methods III (Advection: the semi- Lagrangian technique)

Numerical methods III (Advection: the semi- Lagrangian technique). by Michail Diamantakis (room 012; ext. 2402). In part based on previous material by Nils Wedi and Mariano Hortal and Agathe Untch. Solving the advection problem. “Advection” is a dominant atmospheric process:

terry
Download Presentation

Numerical methods III (Advection: the semi- Lagrangian technique)

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Numerical methods III(Advection: the semi-Lagrangian technique) by MichailDiamantakis(room 012; ext. 2402) In part based on previous material by Nils Wedi and Mariano Hortal and AgatheUntch

  2. Solving the advection problem • “Advection” is a dominant atmospheric process: • Transport of momentum, heat, mass around the globe • Eulerianapproach: the “observer” sits in a fixed point in space watching the flow evolving around – fluid parcels are always evaluated at fixed locations (in practice grid-points) • Time step restrictions for explicit schemes • nonlinear advective terms more complicated equations • Lagrangian approach: an “observer” travels with a fluid element i.e. parcels move freely in space and time (not on fixed points) • Distribution of parcels eventually becomes highly non-uniform and numerical approximations inaccurate

  3. What is a semi-Lagrangian method? • A numerical technique for solving transport problems which combines the Lagrangian “philosophy” with Eulerian “technology”: • It is like a Lagrangian method because it calculates parcel movements along a trajectory • But it is also “slightly” Eulerian as there is a fixed regular grid where data are always mapped on – important flow features can always be resolved at the same accuracy as there is a fixed grid • It gradually evolved to current form from schemes introduced in the ’50s, ’60s and 70s (Wiin-Nielsen, Krishnamurti, Sawyer, Leith, Purnel)

  4. Semi-Lagrangian: an NWP favourite • The semi-Lagrangian (SL) technique is a popular one for solving the NWP advection problem: • Efficient as it is unconditionally stable: time step restriction of explicit Eulerian schemes is lifted • advection time-step is only restricted by accuracy, not by CFL number • Nice property for global models, especially for those using lat-lon grids • A “simplifier”: reduces the non-linearity of the advection problem • Good phase speeds with little numerical dispersion

  5. The departure point Consider the advection equation without forcing: Use Lagrangian derivative approximation: • Finding the “departure point” is an essential part of the technique: • Solution at t+Δt is the field value at the d.p.which is obtained by interpolation of the existing field defined in the regular grid! • Nonlinear term is absorbed by the Lagrangian derivative: simpler linear problem is solved!

  6. Stability in one dimension p • j-p • j • j-p-1 (constant wind) • x α “multiply-upstream” p: integer α is not the CFL number except when p=0, then=> upwind Assuming linear interpolation conduct Von Neuman stability analysis: |λ|≤1 if 0 ≤α≤1 (interpolation from two nearest points) Damping!

  7. A simple Semi-Lagrangian algorithm Solve Given field at time , a grid (arrival pts) and wind field V: • First compute departure point (d.p.) location, e.g. for simple case with constant wind • Using field values at points surrounding interpolate given field to obtain solution at future time , i.e. interpolation operator Accurate calculation of d.p.and an accurate interpolation scheme are essential!

  8. Finding d.p. in “real” atmospheric flows In real flows wind field changes in space and time • To find departure points, solve trajectory equation: where the position and velocity vector (wind) • Second order mid-point rule is often used: • Departure point is computed iteratively two-time level scheme: Trajectory midpoint t-extrapolationneeded: 3-time level scheme: No t-extrapolation but expensive

  9. The trajectory iterative scheme • Consider two time-level (TL) scheme and assume that during a time-step parcels follow straight lines (great circle) trajectories • Define displacement vector: • Iterate discretized trajectory equation: is the known part i.e. the position vector of the “arrival” gridpoint Interpolate V at r-d(k-1)/2 (use linear interpolation) Extrapolate V at t+∆t/2 normally K=2 Smolarkiewicz & Pudikiewicz (J. Atmos. Sci.1992): Convergence requires parcels do not overtake each other (no trajectory crossing): Doesn’t depend on mesh size and less restrictive than CFL for atmospheric flows

  10. Time extrapolated winds and stability • Extrapolating V att+∆t/2 (when computing d.p.) can be a source of weak instability Cordero et al (QJRMS, 2005) • Iterative (expensive) approach: • Time-step once to obtain good estimate of V(t+∆t) • Time-step again, BUT, do trajectory iterations interpolating, i.e. V(t+∆t/2)=[V(t)+V(t+∆t)]/2 rather than extrapolating! • Use Stable Extrapolating Two Time-Level Semi-Lagranian (SETTLS) scheme by Hortal (QJRMS, 2002) T forecast 200 hPa (from 1997/01/04) SETTLS Standard extrapolation

  11. SETTLS for computing trajectories in 2-TL scheme Taylor expansion to second order: and AV: average value along SL trajectory Hence, Compute d.p. by iterating: Interpolate at

  12. SL advection with forcing • Three TL 2nd order midpoint SL scheme for : • Compute by interpolating to mid-point • Two TL trapezoidal (Crank-Nicholson) SL scheme • Three TL schemes are centred in time but are twice as expensive to achieve same accuracy! Outdated … F at midpoint: implicit term

  13. SETTLS and nonlinear terms in ECMWF Discretize equations using a 2nd order, fully centred time scheme separating linear and non-linear terms: Trajectory midpoint L: linear operatorN: non-linear terms Crank-Nicholson can be obtained by one of the two extrapolation schemes discussed: SL averaging SETTLS 2nd order estimate (can be verified by Taylor expansion)

  14. Z V j j j A i x A Y D D M i i X Semi-Lagrangian advection on the sphere Momentum eq. is discretized in vector form (a vector is continuous across the poles, components are not!) Interpolations at departure point are done for components u & v of the velocity vec- tor relative to the system of reference local at D. Interpolated values are to be used at A, so the change of reference system from D to A needs to be taken into account. Trajectories are arcs of great circles if constant (angular) velocity is assumed for the duration of a time step. Rotation matrix Tangent plane projection Trajectory calculation

  15. Treatment of the Coriolis term Extrapolating in time Coriolis terms is unstable! (Temperton in A. Robert memorial volume 1997) For the NWP equation set which includesCoriolisterms there are two stable options: •Incorporate Coriolis terms in the advected field: terms with velocity in momentum eqn • Implicit treatment: Helmholtz eqs partially coupled for individual spectral components => tri-diagonal system to be solved.

  16. Interpolation in the IFS semi-Lagrangian scheme y x x x x x x x x x x x x x x x x x x x x x x • Two important steps in SL algorithm: • Compute departure point (trajectory calculation) • Interpolate advected field to d.p. To obtain: • ECMWF model uses quasi-monotonequasi-cubic Lagrange interpolation , Cubic Lagrange interpolation: Number of 1D cubic interpolations in two dimensions is 5, in three dimensions 21! To save computations: cubic interpolation only for nearest neighbour rows, linear interpolation for rest => “quasi-cubic interpolation”=> 7 cubic + 10 linear in 3 dimensions.

  17. x x x x x x x φmax x φmin Shape-preserving interpolation • Creation of artificial maxima /minima x: grid points x x x x: interpolation point x x • Shape-preserving (quasi-monotone)interpolation - Spline or Hermite interpolation derivatives modified derivatives interpolation - Quasi-monotone cubic interpolation

  18. Non-interpolating SL • Ritchie (1986) developed SL scheme without interpolation: • Why? To eliminate small scale damping associated with interpolations which could affect long-range simulations • How? Split velocity vector in two: component from arrival point to the nearest to the d.p. grid-point and the residual: 1d: x xxx x xxx x xxx No need to interpolate: - Use available field value at nearest point However, - A residual advection needs to be done. This is done using an neutral (undamped) 3-TL scheme x Parcel crosses Consider 3-TL p grid-points scheme • Scheme has dispersive properties of its • Eulerian counterpart • Based on a neutral 3-TL (more expensive)

  19. General remarks • Semi-Lagrangian is particularly efficient for the multi-tracer advection problem: all you need is multiple interpolations on a single departure point! • Cubic Lagrange is a popular option for NWP & climate models • Quasi-monotone filters often used, particularly in tracer fields (humidity, water species, ozone, gases) • Interpolation order and accuracy: for a p-th order interpolation scheme global truncation error in linear advection (constant wind) case is: i.e. smaller timestep doesn’t necessarily improve accuracy! Ho(improves accuracy in the calculation of d.p.)

  20. Factors affecting accuracy of trajectories • Interpolation of velocities when computing trajectories: • Mac Donald (MWR, 1987) demonstrates that interpolation of 1 degree lower than the one used to interpolate to the d.p. (e.g. quadratic when cubic) is sufficiently accurate. • But in practice, when cubic is used for the advected field linear works well for trajectory interpolations and is cheap Temperton& Staniforth (QJRMS, 1987) • In most models a second order integration scheme is used for computing trajectories • 2 iterations are usually enough for converging the trajectory equations

  21. A summary on stability • Fully interpolating SL advection is unconditionally stable provided • Interpolations use the (nearest) points surrounding the departure point • When forcing terms are present, i.e. NWP equation set, implicit treatment is used for fast terms and Coriolis terms (or alternative advective treatment). • In non-interpolating SL schemes using explicit Eulerian advection for the residual velocities, the restriction applies. However, by construction, it is automatically satisfied!(Eulerian advection takes place in distance less than one grid-length)

  22. A note on mass conservation • Semi-Lagrangian formulations are based on the advective form of the equations but can be made mass conserving e.g. Zerroukat (SLICE),Kaas (LMCSL), Lauritzen (CSLAM) • 3 fundamental issues: • The iterative trajectory calculation should (but normally does not) involve the continuity equation. • The conservation properties of the interpolation operator are important. • Typical semi-Lagrangian schemes are point-value based but conservative schemes preserve volume/area and are volume/area-integrated average based. • Cheapest solutions are global mass fixers.

  23. Useful references • Staniforth and Cote (MWR 1990) review paper: “Semi Lagrangian schemes for Atmospheric models” • Hortal (QJRMS 2002): “The development and testing of a new two-time level semi-Lagrangian scheme (SETTLS) in the ECMWF forecast model” • Dale Durran’s book: “Numerical methods for Wave Equations in Geophysical Fluid Dynamics” (1999) • Other references scattered in the slides … • The training course notes

More Related