Transcript Recent developments in Numerical Methods with a report from the ECMWF annual seminar on ‚Recent developments in numerical methods for atmosphere and.

Recent developments in Numerical Methods with a report from the ECMWF annual seminar on
‚Recent developments in numerical methods for
atmosphere and ocean modelling 2013‘
29th WGNE meeting
10-14 March 2014, Melbourne
M. Baldauf (DWD)
M. Baldauf (DWD)
1
ECMWF annual seminar 2013 on ‚Recent developments in
numerical methods for atmosphere and ocean modelling‘
2-5 Sept. 2013, ECMWF, Reading, UK
organizer: Nils Wedi (ECMWF)
invited talks also from many globally modelling centers
http://www.ecmwf.int/newsevents/meetings/annual_seminar/2013/index.html
(last annual seminars on numerical methods: 2004, 1998, 1991, …)
horizontal grid:
N. Wood, G. Zängl, B. Skamarock, H. Tomita, J. Thuburn
vertical discret.:
G. Zängl, M. Hortal, M. Baldauf
Finite element based: C. Cotter, F. Giraldo, M. Baldauf
time integration:
P. Benard, R. Klein, S.-J. Lock
unstructured meshes: J. Szmelter
reduced equation systems: P. Smolarkiewicz, R. Klein
Physics-dynamics-coupling: S. Malardel
M. Baldauf (DWD)
2
Classical requirements for a dynamical core:
• Stability
• Accuracy (up to a certain order; at least 2)
• Efficiency
Today, additional requirements:
• conservation of certain variables
• „mimetic“ properties (discretization reproduces some exact analytical properties)
• well-balancing
• scalability
• efficient on (quite) different computer architectures (CPU, GPU,…)
M. Baldauf (DWD)
3
Development and projection of HPC over the years (source: www.top500.org)
1 EFlops
3*106 cores
1 PFlops
~104 cores
estimated
~108 cores
1 TFlops
factor 1.8 per year
= factor 2 per 14 months
1 GFlops
2013
2003
1993
 code scalability is crucial
M. Baldauf (DWD)
4
Since the clock rate of supercomputers does not increase
 increase of computer power only by parallelisation
 scalability of the code is very important
Many current global models (still) use a lat-lon grid
 clustering of grid points around the pole
 get rid of the strong time step restriction by combination of
Semi-Lagrangian-advection and semi-implicit time integration
 SL for high Courant numbers requires intensive communication
 detrimental for scalability
 look for horizontally more uniform grids
(lat-lon is not a problem for limited area models)
M. Baldauf (DWD)
5
Horizontal grids
review of horizontal grids on the sphere: Staniforth, Thuburn (2012) QJRMS
Why bothering with all this horizontal grid stuff?
• more uniform grids (see above)
• grid staggering
•  degrees of freedom (DoFs) per cell and/or variable
see also a general overview about ‚computational modes‘ by J. Thuburn
M. Baldauf (DWD)
6
Why grid staggering?
Dispersion relation of the linear 1D wave equation
Frequency 
Phase-, group-velocity
 x/c
=0 for
2x waves
vph
unstaggered
negative
group
velocity
vgr
k x
grid
k x
 x/c
tstagg =
½ tunstagg
vph
staggered
vgr
k x
k x
M. Baldauf (DWD)
8
M. Baldauf (DWD)
9
N. Wood
(UK MetOffice)
M. Baldauf (DWD)
10
N. Wood
(UK MetOffice)
unfortunately low
accuracy of
discretizations
M. Baldauf (DWD)
11
N. Wood
(UK MetOffice)
B. Skamarock
(NCAR)
M. Baldauf (DWD)
12
Currently used grids for operational models
• UK MetOffice: currently lat-lon grid (‚New Dynamics‘, ‚ENDGame‘)
GungHo: more uniform grids (example see above)
• NCAR: MPAS uses icosahedron  dual grid  hexagonal C-grid
• DWD: currently (GME): icosahedron (however, regular grid on 10 diamonds)
new: ICON uses icosahedron with triangle refinement, triangular C-grid
• Canada: currently lat-lon
plans: Yin-Yang grid in 2015
• China (CMA): FD on regular lat-lon grid  SISL scheme
new: (multi-moment) FD; Yin-Yang grid or cubed sphere.
• …
M. Baldauf (DWD)
13
B. Skamarock
(NCAR)
this is still a
conformal grid
(i.e. no hanging nodes)
M. Baldauf (DWD)
14
alternatively: hide most of the problems with the horizontal grid
behind a spectral representation 
• ECMWF: horizontally spectral (vertically FE)
• MeteoFrance: currently spectral models (share dyn. core with ECMWF)
however, vision for FE discretization horizontally
• JMA: GSM uses spectral approach
• Brazil: spectral / lat-lon  SI integration
• …
M. Baldauf (DWD)
15
The Integrated Forecasting System (IFS)
technology applied at ECMWF for the last 30 years …
A spectral transform, semi-Lagrangian, semi-implicit
(compressible) (non-)hydrostatic model
Big obstacle could be removed:
Legendre transform is O(N^3), fast Legendre Transform O(N^2 log N)
(Wedi et al. (2013), Tygert (2008, 2010) )
Similar statements hold for Aladin (biperiodic Fourier series)
N. Wedi (ECMWF)
16
N. Wedi (ECMWF)
M. Baldauf (DWD)
17
by G. Mozdzynski
M. Baldauf (DWD)
18
Vertical grids
issues:
• vertical finite element representation
• proper vertical averages in vertically staggered grid for FD are necessary
M. Baldauf (DWD)
19
Improvement of the vertical discretization
Vertical (Lorenz)-staggering in COSMO (and a few other models):
Half levels (w-positions) are defined by a stretching function zk = f(k);
Main levels (p‘, T‘-pos.) lie in the middle of two half levels
Arithmetic average from half levels to main level:
Weighted average from main levels to half level
Derivatives always by centered differences (appropriate average used before)
The above mentioned staggering requires proper averaging; example …
M. Baldauf (DWD)
20
Truncation error in stretched grids
Divergence – grid stretching variant B
Divergence with weighted average
of u (and v) to the half levels:
s·dz
dz
1/s ·dz
Divergence with only arithmetic average of u (and v) to the half levels:
not a consistent discretization if s1 !
M. Baldauf (DWD)
21
Time integration
• Meteo France (P. Benard): 3-time level Semi-implicit Leapfrog based integration
+ iteration step (ICI-scheme) (on spectral model)
• DWD (ICON): 2-time level Predictor-Corrector –scheme (HEVI in each step)
COSMO: 2-timelevel Runge-Kutta (split explicit)
• NCAR: 2-timelevel Runge-Kutta (split explicit) for WRF and MPAS
• UK MetOffice: currently SI-SL scheme
outlook: IMEX Runge-Kutta under consideration
• JMA: GSM uses SI-SL scheme (for the spectral model)
• …
M. Baldauf (DWD)
22
M. Baldauf (DWD)
23
M. Baldauf (DWD)
24
M. Baldauf (DWD)
25
Advection schemes
For tracer advection beyond sufficient accuracy, two ingredients are
seen as important:
• conservation  finite volume discretization
• consistency with continuity equation
(=the advection equation for q should reproduce those of  if q=1)
easy, if the same advection scheme is used for  as for all tracers;
however, NICAM, ICON, …use different time steps for continuity
equation and tracers  use time averaged mass fluxes! (talks by
H. Tomita, G. Zängl)
M. Baldauf (DWD)
27
Semi-Lagrangian advection (overview given by M. Diamantakis)
as used for the dynamical core in:
•
•
•
•
•
•
•
•
•
ARPEGE(MeteoFrance)/
IFS(ECMWF)/ALADIN
UM(UKMO)
HIRLAM
SL-AV(Russia)
GEM(Environment Canada)
GFS(NCEP)
GSM(JMA)
…
For tracer advection, local conservation properties are not easy to achieve:
e.g. SLICE-3D (Zerroukat, Allen (2012) QJRMS) implemented in
ENDGame version of UM (but computationally expensive)
on unstructured grids: conservative SL probably too expensive.
(conservative) SL schemes on structured grids may be beneficial for many tracers
M. Baldauf (DWD)
28
Finite-element based methods
several Continuous Galerkin (CG) / spectral elements /
Discontinuous Galerkin (DG) developments:
•
•
•
•
•
NUMA (Naval postgraduate school, Monterey) F. Giraldo, …
HOMME (NCAR) R. Nair, …
CAM-SE (Sandia Nat. lab.) M. Taylor
COSMO (DWD) first steps by D. Schuster, …
…
general remarks in talk of C. Cotter
M. Baldauf (DWD)
29
linear 1D wave equation 
C. Cotter (Imperial Coll.)
M. Baldauf (DWD)
30
C. Cotter (Imperial Coll.)
 Ladyzhenskaya/Babuska/Brezzi (LBB) condition guarantees compatibility
M. Baldauf (DWD)
31
Discontinuous Galerkin (DG) methods in a nutshell
 weak formulation
(increases solution space)
Finite-element ingredient
Finite-volume ingredient
 talk by F. Giraldo
From Nair et al. (2011) in ‚Numerical techniques …
e.g.
Cockburn, Shu (1989) Math. Comput.
Cockburn et al. (1989) JCP
e.g. Legendre-Polynomials
Lax-Friedrichs flux
 ODE-system for q(k)
M. Baldauf (DWD)
32
1D wave expansion with a Discontinuous Galerkin (DG) discretization
Literature:
Hu, Hussaini, Rasetarinera (1999) JCP: 1D advection-, 2D wave-equation
Hu, Atkins (2002) JCP: non-uniform grids  k=k()
Ainsworth (2004) JCP: multi-dim. advection equation
M. Baldauf (DWD)
33
Numerical dispersion relation for the 1D wave equation
DG with p=0,1,2,3
(=c used)
p
Re  x/c
max ||x/c
0
1
1
3.9
2
7.51
3
11.83
4
16.86
5
22.58
6
28.96
10
60.75
15
113.68
k x
Im  x/c
 max || x/c  1 + 2.6 p + 0.33 p²
increases slightly stronger
than linear with p.
Choose not too large p!
k x
M. Baldauf (DWD)
34
F. Giraldo (NPS)
M. Baldauf (DWD)
35
Other remarks
• „WGNE table about computer ressources 2014“
up to now, I have received contributions from
Chiashi (Japan), Hoon (Korea), Xueshun (China), Ayrton (Canada),
Jean-Noel (ECMWF) and DWD
• shall we continue the survey about dynamical cores (by Mikhail Tolstyk) ?
M. Baldauf (DWD)
36