1 / 18

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

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

trapper
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 Nils Wedi (room 007; ext. 2657) Based on previous material by Mariano Hortal and Agathe Untch

  2. 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

  3. 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!

  4. 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

  5. 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

  6. 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.

  7. 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’

  8. 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’

  9. 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.

  10. 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

  11. 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

  12. 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.

  13. 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.

  14. • 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

  15. 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

  16. 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)

  17. 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

  18. 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.

More Related