180 likes | 296 Views
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
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 conservation • Semi-Lagrangian formulations are based on the advective form of the equations but can be made mass conserving (e.g. Zerroukat 2007; Kaas 2008) • 2 fundamental issues: • The iterative trajectory calculation should (but normally does not) involve the continuity equation. • The conservation properties of the interpolation operator are important.
• 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
Example Let us apply each of the above schemes to the equation whose analytical solution (with appropriate initial and boundary conditions) is: Z = Re( Ae-ikx eωt) with ω=ikU0-k2K WARNING: the three-time-level scheme applied to the diffusion eq. has an absolutely unstable numerical solution! With the values A=1, k=2π/100, K=10-2, the r.m.s. error with respect to the analytical solution (before the unstable numerical solution grows too much) grows linearly with time. After 200 sec of integration, the error is: 5×10-4Δt for the split treatment 5×10-4Δt for r.h.s. at departure point 5×10-8 (Δt)2 for the centered scheme
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.