Presentation Transcript

1. Numerical methods IV(time stepping)

   by MichailDiamantakis(room 012; ext. 2402) In part based on previous material by Nils Wedi,AgatheUntch, Mariano Hortal
2. What makes a good numerical scheme? The usual essential requirements for any application: Good stability Large ∆t can be used without the model blowing up Accuracy and not stability should be what restricts ∆t Ideally “unconditional stability” – no ∆t restrictions! Good accuracy: At least 1st order accurate i.e. global error → 0 as ∆t → 0 (convergence = at least zero stability + order≥1) Higher order is desirable: fast convergence to the “analytical” solution as ∆t is decreased In NWP and CFD very high order is unnecessary
3. Enhancing stability Simplifying assumptions in the equation model: Hydrostatic approximation Filters fast acoustic modes ⇒CFL restriction less severe. Sufficient for relatively coarse global resolutions (≥ 10km) Anelastic approximation also filters sound waves but allows buoyancy (prognostic equation for vertical motion) Full equations can be used WITH An unconditionally stable semi-implicitscheme retaining the slow advectivecomponent undistorted but reducing the propagation speed of sound waves and fastest GW modes Split-explicitintegration of the full equations, since explicit NOT practical (much slower)
4. Accuracy Slower solution modes such as Rossby & slow gravity waves must be resolved – important for good weather predictions Fast acoustic & gravity waves can be supressed (damped) – these carry little energy and are not important: Ideally a numerical scheme is needed which damps these fast modes while is not damping the other significant waves Very high orders of accuracy in time-stepping are not necessary: In any model, errors from other model components such as parametrizations are usually large enough to make use of very high integration order not a practical option. Usually order 2-3 is good enough.
5. Choices for numerical implementation Do we want to avoid solving an elliptic equation (pressure solver)? Go for splitting into fractional steps eg. split-explicitSkamarock and Klemp (MWR, 1992); Durran (1999) If we want to solve an elliptic equation Projection method; Durran (1999) Semi-implicit Durran (1999); Cullen et. al.(QJRMS, 1990); Benard et al. (QJRMS, 2005) Preconditioned conjugate-residual solvers (eg. GMRES) or multigrid methods for solving the resulting Poisson or Helmholtz equations; Skamarock et. al. (MWR,1997) Direct Methods; an example is IFS using the spectral transform method
6. Spectralsemi-Lagrangiansemi-implicit (compressible) viable ? Action: Establish if okay to slow-down acoustic waves Action: Establish if spectral viable at ultra-high resolution yes no Action: Establish if semi-Lagrangian viable for CFL larger than 1 with permitted/resolved convection Action: Firmly establish efficiency of 3D elliptic solvers on MPP computers (linear with number of gridpoints) Action: Invest in grid/finite-volume methods no yes Action: Can advantage for SL schemes be maintained on MPP computers and issues with lower boundary overcome ? yes no Improve on existing NH system: Investigate methods to reduce cost of nonlinear equation solver andspectral transform explicit Options can be explored: -Compressible semi-implicit Anelastic/pseudo-incompressible Coupled model equations no Fast/slow split Action: Invest in flux-form Eulerian formulation Single step, horizontally explicit/vertically implicit
7. Time splitting Split slow fast in two sub-problems: to obtain a solution expressed as the product between a slow and a fast term:
8. Split-explicit integration Skamarock and Klemp(MWR, 1992); Durran (1999); Domsand Schäettler (DWD Tech. Rep. 1999) Based on previous analysis, a discrete scheme can be applied as follows: Solve slow part of problem Solve fast part of problem doing M-successive integrations with ∆t/M and φs initial conditions e.g. implemented in popular limited-area models: DWD Lokal Modell, WRF model Disadvantage: splitting error is an additional source of error
9. Example: 1d gravity wave equations Fluid mean depth Linearized: Perturbation from mean depth Seeking linear analytic solution of the form: implies a phase speed:
10. 2D shallow water equations : advection : gravity-wave (or sometimes called ‘adjustment’) term CFL≤1⇒ ∆x= ∆y=∆s advection adjustment combined In the atmosphere sec in synoptic-scale models
11. Explicit time-stepping • Leap-frog explicit scheme un-staggered grid arrangement xxxxxxx staggered grid arrangement xoxoxoxoxoxox Neutral (no damping) + 2nd order phase + dispersion errors Von Neuman stability: (for un-staggered mesh) Solution is a combination of a physical and a “parasitical” computational mode which can be damped by use of a time filter, e.g. Asselin filter
12. Increasing the allowed timestep • Forward-backward scheme forward backward Doubles the leapfrog timestep Neutral (no damping) • Runge-Kutta schemes such as RK3 used in WRF A high order FD scheme used to estimate derivatives Almost doubles (1.62) Δt when 3rdorder spatial discretization used [Also worth mentioning HEVI and promising IMEX Runge-Kutta (Pareschi & Russo 2005)]
13. Split-explicit time-stepping (index m for small t-steps starting from slow solution) Fast Slow Slow part canbe integrated with M-times bigger ∆t Restriction not severe as it deals with slow advective part Severe ∆t restriction deals with fast part Potential drawbacks: splitting errors, conservation However recent advances for NH NWP suggested in (Klempet al 2007) Note: The fast solution may be computed implicitly.
14. x x x x x x x x x Semi-implicit on 2D shallow water without rotation N Substitute first two equations to 3rd to obtain: E W j S a Helmholtz equation needs to be solved! Stability: now onlylimited by the advection terms coupling with Semi-Lagrangian lifts any stability ∆t restriction
15. Design of semi-implicit methods Treat all terms involving the fastest propagation speeds implicitly (acoustic waves, gravity waves) Assume that the energy in those components is negligible Consider the solvability of the resulting implicit system, which is typically an elliptic equation with coefficients: constant e.g. IFS(ECMWF)/Arpege/Aladin NH constant in time e.g. MC2 model non-constant e.g. UK Met Office NH model, EULAG model Large time-steps typically used by semi-implicit methods result in to phase errors in fast waves of the solution
16. Example: compressible Euler equations Simplified & modified version of Davies et al (QJRMS, 2005): inviscid& adiabatic w/o rotation, SL continuity, fully-interpolating SL θ advection on xyz plane
17. Semi-Lagrangian semi-implicit discretization 1/2≤ α≤ 1: semi-implicit (off-centring) weight - controls damping α=1/2 Crank-Nicholson ⇒ neutral (no damping) However, often a small off-centring amount may be useful to control numerical noise (damping unwanted oscillations artifact of the numerical discretization)
18. Deriving Helmholtz equation Substitute in momentum equations, calculate using and substitute in continuity to derive a Helmholtz equation for : Divide continuity by and substitute above expression Variable coefficient Helmholtz linear system of equations for Option to iterate this non-linear solver in some models (for convergence)
19. Semi-implicit time integration in IFS Choice of which terms in RHS to treat implicitly is guided by the knowledge of which waves cause instability because they are too fast (violate the CFL condition). In a hydrostatic model, fastest waves are horizontally propagating external gravity waves (long surface gravity waves), Lamb waves (acoustic wave not filtered out by the hydrostatic approximation) and long internal gravity waves => implicit treatment of the adjustment terms L= linearization of part of RHS (i.e. terms supporting the fast modes) => obtain a linear system of equations that can be solved analytically in a spectral model with respect to the unknown variables representing field values at future times.
20. Define: “semi-implicit averaging term” 2-TL SLSI integration of hydrostatic PE ε=0=> semi-implicit ε=1=> fully implicit ε=-1=> explicit X : advected variable F: right-hand side L: linearized part of RHS containing fast terms which will be treated implicitly “0”: indicates value at dep. point “1/2”:indicates value at mid-point t+Δt/2 “+”: indicates value at arrival point t+Δt Here, for simplicity, we have assumed advective treatment of Coriolis i.e.
21. Deriving constant coefficient Helmholtz equation Extract linearized fast terms L and discretize: Rewrite using definition for : integration matrices [Ritchie et al (1995) MWR 123] Take divergence of momentum Eliminate to obtain Helmholtz equation for
22. One equation for each because in spectral space Solving Helmholtz equation Vertically coupled set of Helmholtz equations.Coupling through Decouple equations by diagonalizing Once D+ has been computed, it is easy to compute the other variables at “+” (back-substitution). Forecast at t+Δt complete!
23. Time stepping summary Currently two dominant approaches in atmospheric modelling: Split-explicit (also HEVI: Horizontally explicit, vertically implicit) Semi-implicit Due to future H/W requirements there is currently a lot of research on what will be the best (scalable) time-stepping technique on exascale machines.