170 likes | 303 Vues
This lecture by Nils Wedi, hosted in room 007, delves into the semi-Lagrangian technique for solving the advection equation, building on prior materials by Mariano Hortal and Agathe Untch. The method focuses on avoiding quadratic terms by calculating time evolution along trajectories. Key topics include stability analysis, interpolation methods like cubic splines and quasi-monotone interpolation, and the nuances of applying semi-Lagrangian schemes in different dimensions. This presentation also addresses challenges in stability, trajectory calculations, and the treatment of the Coriolis term in atmospheric modeling.
E N D
Numerical methods III(Advection: the semi-Lagrangian technique) by Nils Wedi (room 007; ext. 2657) Based on previous material by Mariano Hortal and Agathe Untch
x x x x x x x x x Advection: The semi-Lagrangian technique material time derivative or time evolution along a trajectory thus avoiding quadratic terms; disadvantage: not flux-form! x x x x x x From a regular array of points we end up after Δt with a non-regular distribution x x x Semi-Lagrangian: (usually)tracking back Solution of the one-dimensional advection equation: computing the origin point via trajectory calculation origin point interpolation
e.g. McDonald (1987) Stability in one dimension Linear advection equation without r.h.s. p Origin of parcel at j: X*=Xj-U0Δt x j α “multiply-upstream” p: integer Linear interpolation α is not the CFL number except when p=0, then=> upwind Von Neumann: |λ|≤1 if 0 ≤α≤1 (interpolation from two nearest points) Damping!
Cubic spline interpolation S(x) is a cubic polynomial - S(xj)=φj at the neighbouring grid points - ∂S(x)/ ∂x is continuous - ∫d2S/dx2 dx is minimal Then: S(x)=Dj-1(xj-x)2(x-xj-1)/(Δx)2-Dj(x-xj-1)2(xj-x) /(Δx)2 + +φj-1(xj-x)2[2(x-xj-1)+ Δx] /(Δx)3+ φj(x-xj-1)2[2(xj-x)+ Δx] /(Δx)3 where (Dj-1+4Dj+Dj+1)/6=(φj+1- φj-1)/2 Δx
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 and quasi-monotone interpolation - Spline or Hermite interpolation derivatives modified derivatives interpolation - Quasi-monotone interpolation
Interpolation in the IFS semi-Lagrangian scheme y x x x x x with the weights x x x x x x x x x x x x x x x x x 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 on computations: cubic interpolation only for nearest neighbour rows, linear interpolation for rest => “quasi-cubic interpolation” => 7 cubic + 10 linear in 3 dimensions.
o o x x 3-t-l Semi-Lagrangian schemes in 2-D L: linear operatorN: non-linear function • Interpolating G x x x x x x x x x x x x I x Two interpolations needed • Ritchie scheme U=U*+U’ V=V*+V’ G x x x x x x x x x x x x I’ 2V*Δt o’ 2V’Δt The three of them are second-order accurate in space-time • Non-interpolating Average the non-linear terms between points G and o’
Stability of 2-D schemes • In the linear advection equation the interpolating scheme is stable provided the interpolations use the nearest points • In the linear shallow-water equations (with rotation), treating the linear terms implicitly, the stability limit is • In the two non-interpolating schemes the stability is given by Δt f ≤1 Coriolis term (kU’+lV’) Δt ≤1 Which is always true due to the definition of U’ and V’
Treatment of the Coriolis term In three-time-level semi-Lagrangian: • treated explicitly with the rest of the rhs In two-time-level semi-Lagrangian: Extrapolation in time to the middle of the trajectory leads to instability (Temperton (1997)) Two stable options: • Advective treatment: (Vector R here is the position vector.) • Implicit treatment : • Helmholtz eqs partially coupled for individual spectral components => tri-diagonal system to be solved.
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 (because 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. Tangent plane projection Trajectory calculation
Iterative trajectory calculation V1Δt V0Δt x x x x x x x x x x x x x x x x x x x x r0 r1 can be taken as trajectory straight line or great circle assumed constant during 2Δt or implicit
Iterative trajectory computation (1 dimension) r(n+1)=g-V0(n)Δt Where, for simplicity, we have taken a 2-time-level scheme and taken the velocity at the departure point of the trajectory Assume that V varies linearly between grid-points V=a+b.r b = r (n+1) = g - aΔt - Δt b r(n) For this procedure to converge, it must have a solution of the form r = λn + K; (| λ| < 1) |b|Δt < 1 Substituting, we get K=(g - a Δt)/(1 + b Δt) and λ = -b Δt more generally for three dimensions this translates to the determinant of a matrix: Lipschitz condition (Smolarkiewicz and Pudykiewicz, 1992) Physical meaning: To prevent trajectory intersections !!! It is in general less restrictive than the CFL condition and it does not depend on the mesh size.
Mass and tracer conservation • Semi-Lagrangian formulations are based on the advective form of the equations but can be made mass conserving (e.g. Zerroukat 2007; Kaas2008, Tolstykh 2012) • 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. • Cheaper solutions are global or local fixers applied.
• Three-time-level schemes - centered (second-order accurate) scheme - split in time (first-order accurate) - R at the departure point (first-order accurate) Semi-Lagrangian advection with rhs the r.h.s. R can be evaluated by interpolation to the middle of the trajectory or averaged along the trajectory: RM(t)={RD(t)+RA(t)}/2
with Extrapolation in time to middle of time interval Semi-Lagrangian advection with rhs Two-time-level (second-order accurate) schemes : Unstable! => noisy forecasts Forecast of temperature at 200 hPa (from 1997/01/04)
With and Forecast 200 hPa T from 1997/01/04 using SETTLS Stable extrapolating two-time-level semi-Lagrangian (SETTLS): Taylor expansion to second order
Trajectory computation with SETTLS Mean velocity during Δt The Taylor expansion from which we started is: which represents a uniformly accelerated movement Note: The middle of the trajectory is not the average between the departure and the arrival points.