540 likes | 750 Vues
Development of a new fast waves solver for the Runge-Kutta scheme. M. Baldauf Deutscher Wetterdienst, Offenbach. COSMO-General Meeting 05.09.2011, Rome. Motivation.
E N D
Development of a new fast waves solverfor the Runge-Kutta scheme M. Baldauf Deutscher Wetterdienst, Offenbach COSMO-General Meeting 05.09.2011, Rome
Motivation • Internal DWD project 'COSMO-DE L65' (increase # of vertical levels 50 65):(goal: better representation of convection and its initiation)But several numerical problems occur (lower BC for w, vertical advection scheme, tracer advection scheme, ...)Speculation: discretizations of fast processes not yet entirely consistent • DWD-'Extramurale Forschung' project 'Numerische Raumdiskretisierungsverfahren hoher Ordnung für das COSMO Modell' (A. Will/J. Ogaja)Prerequisite: use of only centered finite difference and centered averaging operators
‚Fast waves‘ processes (p'T'-dynamics): artificial divergence damping sound buoyancy D = divv fu, fv, ... denote advection, Coriolis force and all physical parameterizations integration procedure: horizontally forward-backward, vertically implicit
Integration procedure: horizontally: forward-backward (Mesinger (1977) Contr. Phys. Atm.) vertically: implicit (Klemp, Wilhelmson (1978) MWR, Wicker, Skamarock (1998, 2002) MWR) With the possibility of a 3D divergence damping, stability of the whole RK3-scheme (p'T'-dynamics) was shown in Baldauf (2010) MWR
1. Improvement of the vertical discretization Averages from half levels to main level: Averages from main levels to half level with appropriate weightings (!): centered differences (2nd order if used for half levels to main level) G. Zängl could show the advantages of weighted averages in the explicit parts of the fast waves solver. New: application to all vertical operations (also the implicit ones)
2. 'Strong conservation form' of the divergence operator Divergence operator used up to now: Strong conservation form: ~ d/dt proper BC Discretization of metric terms more compact expressions in the strong conservation form Doms, Schättler (2002) COSMO Sci. Doc. (I), Prusa, Smolarkiewicz (2003) JCP
3. Isotropic divergence damping Gassmann, Herzog (2007) MWR recommend the use of the complete (=isotropic) 3D divergence damping more realistic dispersion relation for sound and gravity waves The following tests compare the current FW solver (fast_waves_rk.f90, version COSMO 4.18) with the new FW solver (fast_waves_sc.f90)
3.a Linear flow over mountains current FW FW new U0=10 m/s, N=0.01 1/s, hmax=25 m test case definition in: Schär et al. (2002) MWR Linear analytic solution (black contours) : Baldauf (2008) COSMO-Newsl.
3.b Nonlinear flow over mountains current FW FW new U0=10 m/s, N=0.01 1/s, hmax=750 m
3.b Nonlinear flow over mountains current FW FW new model abort after 6 h U0=10 m/s, N=0.01 1/s, hmax=900 m
Results of idealised test cases see: COSMO-user Seminar, March 2011 SRNWP-workshop Bad Orb, May 2011 • Sound wave expansion • Linear Gravity wave in a channel (Skamarock, Klemp, 1994) • Linear flow over mountains (compare with analytic solution) • Non-linear flow over a mountain • mountain in a steady atmosphere • moist warm bubble test (Weisman, Klemp, 1982) • dry cold bubble (Straka et al.,1993) All these idealised tests are simulated with either similar accuracyor slightly better.
Real simulations ... ... at the beginning caused several problems, mostly in connection with the divergence damping near the bottom.
Quasi-3D - divergence damping in terrain following coordinates Stability criterium: • in particular near the bottom (x, y >> z) a strong reduction of div is necessary! • This violates the requirement of not too small div in the Runge-Kutta-time splitting scheme (xkd~0.1 (Wicker, Skamarock, 2002), • in Baldauf (2010) MWR evenxkd~0.3 is recommended). • Otherwise divergence damping is calculated as an additive tendency (no operator splitting) a certain 'weakening' of the above stability criterium is possible
… additionally necessary for stability: • don't apply divergence damping in the first small time step!(J. Dudhia, 1993 (?)) • boundary treatment in the term (d/dt) / of the divergence in ‚strong conservation‘ form:the obvious d/dt|ke+1=0 leads to a divergence which is (near the ground) about a factor 10 larger than a numerical approximation of this term.
Boundary condition for the Euler equations = free slip condition at the bottom (and at the top) 1.) Extrapolation of u, v to the ground (ke+1/2):this improved the pressure bias in COSMO-EU Options now: extrapolation 0th order (= 'true' free slip BC), 1st or 2nd order. 2.) Discretisation of dh/dx, dh/dy: up to now: analogous to the order of the advection operator (upwind 5th order) now: centred diff. 4th order necessary ke-1/2 w ke u ke+1/2 u w
'COSMO-EU, 05.01.2011, 0 UTC', PMSL after 78 h (Exp. 8230) Diff. Exp. 8230 Routine during the simulation a negative pressure bias of about -0.5 hPa/3d develops (remark: the operational COSMO-EU was nearly bias free in Jan. 2011)
(12.01.2011, 0 UTC run, stand alone) pressure bias in COSMO-EU FW_current FW_new
(12.01.2011, 0 UTC run, stand alone) pressure bias in COSMO-EU FW_current FW_new 0., 1., 2. orderExtrapolation BC for w: extrapolation of u, v to the lower half level (=surface) does not have any influence to the mean surface pressure
(12.01.2011, 0 UTC run, stand alone) pressure bias in COSMO-EU FW_current FW_current ldyn_bbc=F FW_new 'dynamical bottom BC for p' helps in producing the correct surface pressure (at least the mean value) in the current FW solver
Boundary treatment of p'/ (or of p'/z) at the lower boundary: 1.) one sided finite difference (G. Zängl): p'/ = 0p'(ke) + 1p'(ke-1) + 2p'(ke-2) 2.) ‚dynamical bottom BC‘ (A. Gassmann, 2004, COSMO-Newsl.) From and derive a condition for p'/ .
(12.01.2011, 0 UTC run, stand alone) pressure bias in COSMO-EU FW_new with ldyn_bbc=T FW_newwith ldyn_bbc=F 'dynamical bottom BC for p' also improves pressure bias in the new fast waves - solver
influence of the FW-solver to 'indirect' variables like precipitation is quite small example: COSMO-DE, 12.01.11, 0 UTC Diff. FW current FW new
influence of the FW-solver to 'indirect' variables like precipitation is quite small example: COSMO-DE, 28.04.11, 0 UTC Diff. FW current FW new
COSMO-DE, 12.01.2011, 0 UTC run, ldyn_bbc=.TRUE. current FW solver new FW solver very strong inversion in an Alpine valley
COSMO-DE, 12.01.2011, 0 UTC run, ldyn_bbc=.FALSE. current FW solver new FW solver
Which benefits are currently visible ? • The original intention to develop the new fast wave solver • was to produce more consistent dynamic fields. • Indeed in some situations the stability seems to be slightly higher, examples: • a model crash of COSMO-2 at 16.06.2011, 0 UTC runs stable with the new FW solver • a model crash of COSMO-DE at 12. July 2011, 6 UTC run couldbe repaired by the use of Bott2_Strang, but also alternatively by theuse of the new FW solver • a simulation with high resolution of 0.01° only runs stable with thenew FW solver • the new solver also runs without crash during 01.-31. Jan. 2011in both a COSMO-EU and a COSMO-DE setup But: of course the fundamental difficulty of split explicit schemes with steep orography remain
crash with the operational COSMO-DE at 12. July 2011, 6 UTC model crash with the current FW stable simulation with the new FW unrealistic high value of qr in one grid box; strongly deformed wind field
simulated radar reflectivity COSMO-run with a resolution of 0.01° (~ 1.1km) 1700 * 1700 grid points model crash after 10 time steps with the current fast waves solver stable simulation with the new FW from Axel Seifert (DWD)
Efficiency: on NEC - SX9: new fast_waves_sc: reaches ~20 GFlops and needs about 30% more computation time than the current fast waves solver (~18 GFlops) a COSMO-EU run needs about 5% more time but it is not yet optimized for Intel processors (cache based): it takes about 80% more computing time (reason?) higher computation time can be expected due to - more vertical weightings - exchange of p' and div v
Summary • New fast waves solver with • improved vertical discretizations • strong conservation form of divergence • (optional: 3D isotropic divergence damping) • Idealised test cases (stationary/unstationary, linear/nonlinear, with/without orography) are simulated with either similar accuracyor slightly better • runs stable in all inspected cases (COSMO-EU, COSMO-DE duringthe whole Jan 2011; and in selected cases) • for both FW solvers holds: dynamical bottom boundary condition necessary to reduce pressure bias in COSMO-EU, but can have detrimental influence on COSMO-DE • satisfying optimization for NEC SX9 achieved
Outlook • extensive verification • efficiency: • on NEC-SX9: further increase probably only with help of NEC specialists • on Intel/cache based: better inspection tools necessary (valgrind, ...) • add currently available features:lateral radiation BC; p'-dynamics solver; lower BC for w via 'RK advection of height', ... • application to the COSMO-DE L65 setup • Closer examination of the influence of (an-)isotropic divergence damping
Neu gegenüber dem bisherigen fast waves -Löser: • Verbesserung in der vertikalen Diskretisierung: abstandsgewichtete Mittelungsoperatoren • Divergenz in strong conservation form • optional volle 3D Divergenzdämpfung (d.h. in allen 3 Bewegungsgleichungen) • zusätzlich einige 'technische' Verbesserungen; • Erhöhung der Lesbarkeit neue Version fast_waves_sc.f90 (basierend auf COSMO 4.17)
Eine bessere analytische Lösung für die lineare Schwerewelle im Kanal • Analytische Lösung von • Skamarock, Klemp (1994) MWR • ist Boussinesq-approximiert • keine perfekte Übereinstimmung • zwischen analytischer und simulierter Lösung möglich. Eine bessere analytische Lösung als Referenz wäre auch bei Modellvergleichen wünschenswert und wird z.Z. im Rahmen von durchgeführt
Eine bessere analytische Lösung für die lineare Schwerewelle im Kanal Eine analytische Lösung für kompressible, nicht-hydrostatische Euler-Gleichungen kann für eine isotherme (T0=const) Hintergrund-Atmosphäre gefunden werden! Anfangszustand: ruhend (bzw. u=const.) , p‘=0, vorgegebene T‘- (oder ‘-) Verteilung z.B. der Form:
Die zeitliche Entwicklung der Fouriertransformierten Felder lautet (analog für u, T‘, …) und sind gegebene Funktionen von kx, kz, T0, g, cp, cv
Die zeitliche Entwicklung der Fouriertransformierten Felder lautet (analog für u, T‘, …) und sind gegebene Funktionen von kx, kz, T0, g, cp, cv
1. Sound wave expansion test current FW FW new Isothermal atmosphere (T=250 K)
1. Sound wave expansion test current FW FW new ?
1. Sound wave expansion test current FW FW new ?
1. Sound wave expansion test current FW FW new ? Artificial stationary solution (?)
2. Linear gravity wave current FW FW new test definition and (anelastic approx.) analytic solution: Skamarock, Klemp (1994) MWR
2. Linear gravity wave current FW FW new test definition and (anelastic approx.) analytic solution: Skamarock, Klemp (1994) MWR
2. Linear gravity wave current FW FW new test definition and (anelastic approx.) analytic solution: Skamarock, Klemp (1994) MWR
2. Linear gravity wave current FW FW new test definition and (anelastic approx.) analytic solution: Skamarock, Klemp (1994) MWR
4. Mountain in a steady atmosphere with vertical grid stretching current FW FW new
4. Mountain in a steady atmosphere (b) with vertical grid stretching current FW FW new
4. Mountain in a steady atmosphere (b) with vertical grid stretching current FW FW new
5. Testcase Weisman, Klemp (1982) current FW FW new Umax=30 m/s, qv,0,max=13 g/kg