2.1 Hardware complexity

Significantly, after many years of “computational similarity” based on “multi-core” Intel or Intel-like X86 systems, the future looks far more heterogenous than ever. There are a plethora of new architectures available, from evolutions of X86 with more complex vector processing units, such as the now terminated Intel Many Integrated Core (MIC, or “many-core”) line, to new ARM server chips and/or radically different architectures based on graphical processing units (GPUs). More exotic processor families are emerging too: from field programmable gate arrays (FPGAs) to new technologies such as Google's Tensor Processing unit (Jouppi et al., 2017). At the time of writing, one vendor alone, Cray, is offering very different nodes in their HPC line: multi-core CPU alone, many-core nodes, or nodes which combine CPUs with GPUs.

This proliferation of architectures is a direct consequence of the limitations of physics and power constraints in silicon: substantial year-on-year increases in clock speed are no longer possible, and increased compute capacity has to come from innovation in what is on the chip rather than how fast it goes. An additional driver is that the technology choices which arise are driven by economics, not the needs of the HPC community. Further, for big machines, power constraints dominate thinking, with existing CPU architectures consuming too much power. The US exascale programme has set a goal of 20 MW for an exaflop machine, but this would not be achievable with current technology (Kogge, 2015), leading to more pressure on hardware change – change which is likely to require lower power processors and more of them.

With more cores per chip comes a requirement for more memory, which cannot be met by current technology trends, particularly memory with low latency and high bandwidth. The hardware solutions being offered include a range of types of memory, each with their own characteristics requiring new techniques to manage and control memory interference between different cores, agents, and applications that share the memory system (Mutlu, 2015). Whether caching algorithms will provide the appropriate performance or programmers will need tools to directly manipulate data in these systems is not known. Further, this “near memory” is now being coupled with large fast storage “close” to the computing (so-called “burst buffers”), leading to nested hierarchies in a spectrum from “memory” to “archive” and new problems (and opportunities) for the workflow around applications.

In most modelling groups, the “modeller” currently has to have full knowledge of how the entire software stack works with the hardware even if implementation details are buried in libraries. The new architectures present the “full-stack” modeller with many challenges: can one code to exploit the parallelism and the memory? Is there sufficient parallelism in the code itself to exploit the architecture? How can one obtain performance data upon which to make informed choices? Can one code to exploit multiple different architectures (and maintain that code)?

Figure 1Architectural heterogeneity: a schematic comparison between an Intel Broadwell CPU socket (two sockets to a node, joined locally by the quick path interconnect, QPI) and an Intel Knights Landing MIC node. Notable differences include the two-dimensional mesh linking cores and memory (near and far) on the MIC node (with up to 16 GB of near memory L3 cache) and the bidirectional full ring on the CPU linking the cores (and supporting up to 50 MB shared L3 cache). The most important difference, however, is the number of cores: approximately 20 to a socket or 40 to a CPU node and up to 72 on the MIC node. (Not shown in this diagram is the 256 MB L2 cache in the Broadwell cores and their internal vector processing units.)

Download

Figure 1 shows just one comparison of two relatively similar architectures, a multi-core evolution of the traditional Intel CPU processor and a many-core node. This schematic is intended to show both the similarities and differences: both consist of multiple cores, but there are different arrangements of cache and main memory, different ways of accessing vector processing units, and different methods for linking between the cores.

Today, at the most fundamental level, modellers rely on the compiler to hide most of these details, but with the advent of coding for GPUs (of which more below), it has become clear that the better one can express parallelism in the high-level (e.g. Fortran) code, the better job the compiler can do even on traditional CPUs. Traditionally (say 10 years ago), in a weather or climate model one had two levels of concurrency in mind associated with shared and distributed memory.

• Shared memory concurrency involves organising the code so the compiler could be pre-fetching memory while doing calculations, i.e. worrying about loop order and organising loops or array declarations to enable the efficient use of vector hardware. This requires “good coding practice”, either gleaned by direct knowledge about how compilers and hardware work or simply by edict: e.g. “always code nested loops with the innermost loop using the left-hand index¹, and never mind why”.

• Distributed memory concurrency typically requires explicit handling in the science code. Processes (typically on different nodes) cooperate on a common task utilising memory, which is partitioned across nodes (“domain decomposition”), and exchange updates using the message passing interface (MPIs;  Gropp et al., 1996, ideally using libraries that vendors have optimised for their interconnects).

With the advent of the multi-core CPU (moving from dual core, to quad core, to now 20 + core) the modeller had to become aware of strategies to exploit the fact that the processor local cores shared memory and were connected together with lower latency and higher bandwidth than the interconnect between processor sockets and nodes.

In most codes now we see one of two strategies: either MPI is still used across both nodes and node-local cores, or MPI is used in conjunction with OpenMP (Dagum and Menon, 1998). In the latter case OpenMP is used to run some parallel threads on a core (“hyper threading”) or parallel threads on a socket, but the threads are sharing memory, and directives are used to identify where this can occur. However, this strategy of mixing explicit library calls (MPI) and compiler directives (OpenMP) does not always work better than MPI alone, not least because vendor implementations of MPI can implement on-node “messages” with highly optimised libraries exploiting shared memory with performance that can be hard to beat. Conversely, MPI decomposition cannot be changed without paying a synchronisation and data movement penalty, whereas OpenMP can exploit the dynamic scheduling of work and better handle load balancing around (for example) physics parameterisations such as precipitation. Hence, the modeller following a hybrid approach has many issues to consider: how many parallel threads per node (between zero and the number of hyper-threads supportable on the node) and what parts of the code can take advantage of this technique? There is a balance to be struck which depends on knowledge of the code, the hardware, and the quality of the vendor-supplied compilers and libraries, and this balance needs to be re-evaluated for each architecture and model version.

Other programming models, such as the partitioned global address space (PGAS), try to address the problem of the locality of memory with a single abstraction to address both shared and distributed memory. Implementations include the C extension, Unified Parallel C, Coarray Fortran, Chapel, and the (Open)SHMEM library. However, all of these suffer from lack of adoption and lack of universal availability, both of which contribute to the lack of wide use.

OpenMP also provides task concurrency, through which pools of threads are allocated amongst a set of tasks that can be done concurrently. However, in climate and weather models, this feature is not widely used.

Now our hypothetical modellers are dealing with at least two levels of concurrency, probably a new directive syntax for declaring what parts of the code can use OpenMP or a completely different programming model (PGAS). They may also be having to deal with multiple languages (e.g. some combination of C$++$ or Fortran and C) to best use the hardware and libraries available.

Most processors have vector units, and the trend is for these to grow both in size (width) and number in new processors. This provides another hardware feature that needs to understood and exploited and further complicates the task of model development. Then of course there are GPUs, which many think of as “big fat vector processing units”, but ones for which the needed levels of data parallelism are very high and out-of-order memory fetches are very expensive. However, a cursory look at Fig. 2 shows that GPU architectures are themselves complex, with sophisticated memory hierarchies and onboard methods of thread management.

The big difference between the CPU architectures (whether multi-cores or many-cores such as MIC) and the GPUs is the amount of intelligence (complexity providing extra instruction-level parallelism and caches to hide memory latency) in the cores and on the package around the cores and how accessible that intelligence is to the programmer (or how easily it can be exploited by the compiler without exposing hardware concepts to the high-level code). In CPUs, much of the silicon is devoted to supporting complexity in the algorithm (branching multiple instructions, for example) and in GPUs, much of the silicon is devoted to doing the same thing to lots of data. However, the compiler directives that support this parallelism are still rather different and conform to differing, rapidly evolving, and sometimes poorly defined standards, so the programmer potentially has to understand not only both hardware architectures, but also how best to exploit them. Requiring all modellers to have detailed hardware knowledge – and keep such knowledge up to date – is obviously not a tenable proposition.

Compounding all these difficulties is the fact that future processors are likely to have to operate at both low power and high processor counts; both factors are likely to lead to increased node failure rates. Existing software tools expect perfectly reliable hardware, with checkpointing as the only viable strategy to deal with hardware failure, so there will be a need for new mechanisms for handling hardware failure, which will undoubtedly add further complexity to the hardware programming environment.

Figure 2High-level schematic of an NVIDA Kepler GPU (adapted from  NVIDIA, 2012). Up to 15 streaming multiprocessor (SMX) units share a 1.5 MB common dedicated L2 cache, and up to six memory controllers allow access to main memory. Thread control exists across the entire chip (the GigaThread engine) and in each SMX. The SMX units include 192 single precision cores, 64 double precision cores and other special function units, and an internal interconnect so that blocks of 32 threads can be managed internally. They also include 64 KB of L1 memory that can be configured in a combination of shared memory or L1 cache, and 48 KB of read-only L1 memory.

Download

2.2 Software complexity

Weather and climate models traditionally differ primarily in three areas:

1. the duration of a simulation (days vs. years),

2. component complexity (climate models typically include more components associated with longer-timescale processes), and

3. how they are initialised (weather models by complex data assimilation systems and climate models by very long duration “spin-up” runs for components with poor observational coverage).

These differences have become blurred as weather prediction has begun to exploit ocean, snow, and ice models and with the advent of seasonal prediction. But for the purposes of the discussion here, what they have in common is more important anyway:

1. the importance of executing the model at a satisfactory speed (e.g. a year per day);

2. multiple interconnected but often independently developed sub-models (oceans, atmospheres, wave models, etc.), sometimes with extra versions to provide nested grids and or local refinement;

3. an interconnection framework or coupler (to exchange information amongst model components);

4. at least one component with an expensive dynamical core (atmosphere and/or ocean), with

5. a significant number of physics routines (radiation, microphysics, clouds, gravity waves etc.);

6. an internal diagnostics system (to accumulate variables or calculate extra diagnostics beyond the prognostic variables carried in the model components);

7. systems to manage ensemble methods (for the extra simulations needed to address uncertainty); and

8. an I/O subsystem (to manage efficient input and output).

These common attributes can be linked together in very complex ways.

At the high level, Alexander and Easterbrook (2015) showed some permutations of how large components can be coupled together in eight Earth system models and at the same time indicated the relative complexity of the key components. In general, the most complexity was found in the atmospheres of these models, although some Earth system models of intermediate complexity (EMICs) were dominated by oceans. The relative distribution of complexity also differed significantly between models aiming at the chemistry thought to be significant in coming decades and those aiming at longer-timescale processes associated with the carbon cycle.

The interconnection between all the large-scale components requires a specialised framework and/or coupler to provide methods for exchanging fields between components. Valcke et al. (2012) provide a comparison of coupling technologies used in CMIP5, but it is interesting to note that of the six technologies discussed in this paper, only two (ESMF; Theurich et al., 2015, and OASIS, Craig et al., 2017) have any widespread use, and many large institutes feel it necessary to build their own coupling systems.

One key piece of software complexity arises from the a priori necessity to find parallelism within the large components such as the atmosphere and oceans of any but the lowest-resolution models. With only a handful of large-scale components and machines with hundreds to thousands of nodes, exploiting the machine efficiently requires splitting components across nodes exploiting “domain decomposition” – splitting the multi-dimensional (potentially all of longitude, latitude, and height) space into “domains”, each of which is integrated forward each time step on just one node (with the necessity of passing information between nodes after each time step so that “cross-domain” differences can be calculated).

Figure 3A programmer view of data flow and the interaction of dynamics and physics in the atmosphere of three prominent models: ECMWF IFS, GFDL HIRAM, and the Met Office UM. Red boxes show the start and finish of time steps, yellow boxes show dynamics, blue boxes increment physics, beige boxes physics with full values, and external coupling or library routines are shown in the white boxes. Not all prognostic variables and associated routines are shown!

Download

The relationship between the dynamics and physics within an atmosphere is also complex, with even more scope for different arrangements of data flow between internal components and processes. Figure 3 provides a glimpse into the interplay of dynamics and physics in three different atmosphere components from the ECMWF Integrated Forecast System, IFS (IFS 41R1; https://www.ecmwf.int/en/forecasts/documentation-and-support, last access: 29 April 2018), the GFDL HiRAM (2012 public release; https://www.gfdl.noaa.gov/hiram/, last access: 29 April 2018; Zhao et al. , 2009), and the Met Office UM (v6.6; Davies et al., 2005). This figure describes the flow of data associated with model prognostics within a model time step.

IFS is a spectral transform model. The key property of such a model is that due to the properties of the Laplacian operator, equations for vorticity, divergence and diffusion can be easily calculated in spectral space, with advection and physics calculated on a grid distributed in physical space. The model exploits two time levels and a semi-implicit, semi-Lagrangian time-stepping routine. The time step begins with variables held in spectral space, which are then converted back to grid point space for the physics to be executed sequentially, after which they are converted back to spectral space via a fast Fourier transform (FFT) in longitude and a Legendre transform (LT) in latitude to spectral space, for which the Helmholtz equation resulting from the semi-implicit time stepping is solved and horizontal diffusion is computed.

HiRAM, the GFDL high-resolution atmospheric model, is a grid point model based on finite volumes and the cubed-sphere grid. The dynamics uses a two-level split explicit time-stepping scheme. The physics components use an implicit algorithm in the vertical (including fast coupling with land and ocean surfaces) integrated using a tridiagonal solver which splits into a down and up sequence. The down, dry processes are calculated first with the increments from the physics and dynamics solver calculated in parallel and then summed. The down physics processes are calculated sequentially. These are then summed and fed into the sea ice and land surface, and then the up, moist processes are calculated; again the increments from the physics and dynamics solver are calculated in parallel and summed at the end of time step.

The UM is also a grid point model. In this case the physics is split into “slow” and “fast” processes calculated before and after the main dynamics, which also exploits a semi-implicit semi-Lagrangian time-stepping scheme. The slow physics processes use the state of the previous time step, but the faster processes are integrated using the latest state calculated sequentially using the incremented versions.

It can be seen that these three examples show rather different methodologies. The first two use semi-implicit time stepping and the last explicit time stepping; as a consequence they have quite different communications patterns. The data layouts are very different, at least for the dynamics: IFS uses spectral transforms which necessarily means data manipulation between spectral and grid point space. Each has their own take on how to organise the sequence of physical processes and involve many independent pieces of code, each of which could need configuration and compiler- or even hardware-specific optimisation to take advantage of any given hardware environment.

Figure 3 does not describe actual code blocks or expose details of the internal time splitting, and it is important to note that any implicit parallelism shown may not be mirrored in the model code. Nonetheless, there are important differences in information flow which will feedback on opportunities for concurrent parallelism. For example, the UM does not use the updated state of the variables for the slow physical parameterisations so they could be calculated in parallel, and some subcomponents, e.g. radiation, are expensive enough that they may not be executed at every time step. Such subcomponents are large enough to treat, in terms of concurrency, in the same way as full model components. Balaji et al. (2016) show results from experiments with the same code family as HiRAM, doing exactly that. Their results show that this extra concurrency can be exploited in a range of ways (to shorten the radiation time step or to speed up model integrations, albeit with a potential increase in resources used). However, the actual performance of the necessary changes, both in terms of the effect on simulation quality and speed, will differ with the different model approaches to information flow.

Model formulations also take a range of approaches to providing diagnostic, ensemble, and I/O services, with no common practice as far as we are aware beyond some emerging community use of the XIOS I/O server (http://forge.ipsl.jussieu.fr/ioserver/wiki, last access: 29 April 2018). I/O servers provide a method of offloading (primarily) output to one set of processors to happen in parallel with computations happening on other processors. XIOS provides an XML interface which allows not only output offloading, but also some diagnostic processing en route, and work is underway to extend it to support ensemble management as well.

2.3 Where is the concurrency?

The most obvious paradigm for concurrency in weather and climate models is the domain decomposition introduced above, exploiting MPI, but clearly there are multiple levels of parallelism being provided by the hardware and/or exploitable in the software:

1. vectorisation within a GPU core or accelerator (using compiler directive languages such as OpenACC), with or without

2. vectorisation within a CPU core (handled by the compiler, but only if the code is prepared appropriately, possibly by embedding explicit compiler directives);

3. shared parallelism across CPU and accelerators and/or GPUs (e.g. see Mittal and Vetter, 2015);

4. threading providing shared memory concurrency within nodes or sockets (e.g. by inserting OpenMP directives in code);

5. distributed memory concurrency across nodes, either by

   • utilising MPI (traditional “domain decomposition”; directly or with a library such as YAXT or GCOM) or

   • exploiting an implementation of a partitioned global address space, such as Coarray Fortran (e.g. Mozdzynski et al., 2012, or language-neutral interfaces, such as GPI; Grünewald and Simmendinger, 2013);

6. component concurrency, whether using a framework such as ESMF or handled manually with, for example, OpenMP tasks (e.g. Balaji et al., 2016);

7. coupled model components either

   • as part of coupled model executed independently and concurrently using a coupler such as OASIS or executed together and concurrently using a framework (Valcke et al., 2012) or

   • concurrent models running as part of a single executable ensemble, handcrafted (e.g. Bessières et al., 2016) or via a framework such as the GFDL flexible modelling system used in Zhang et al. (2007); and

8. I/O parallelism (using an I/O server such as XML I/O Server, XIOS, or a parallel I/O library such as parallel NetCDF; Li et al., 2003).

Most of these levels of parallelism require decisions as to where and how the code should be refactored and/or restructured to take advantage of them, and it may not be possible to utilise them all in one code stack, not least because on any given platform implementations may not exist or they may be poorly implemented (in terms of performance or feature coverage).

In the discussion to follow we use the nomenclature introduced by Balaji et al. (2016) to term the first three of these modes as “fine-grained” concurrency and the last four as “coarse-grained” concurrency. The key distinction between these modes of concurrency is that the former involve organising codes to take advantage of single instruction multiple data (SIMD) hardware (such as vector units), and the latter involve higher-level assemblages interacting, mainly using the multiple programme multiple data (MPMD) pattern (except for the distributed memory options, which utilise the single programme multiple data, or SPMD, paradigm). We do not categorise threading (4) as either fine- or coarse-grained concurrency, since in practice from a programmer point of view it can be used following either mode (e.g. to parallelise loops in the same way as one parallelises them for a vector processor or to provide extra models of parallelisation for domain decomposition).

The current best practice for addressing these modes of concurrency is to

1. code for hierarchies of parallelism (loops, blocks, domains),

2. use standard directives (OpenMP or OpenACC),

3. optimise separately for many-core GPUs, and

4. try to minimise the code differences associated with architectural optimisations.

However, this is no longer seen as a long-term successful strategy, at least on its own, because the strategies needed for different architectures may be fundamentally different (e.g. for memory optimisations). With the advent of exascale systems, entirely new programming models are likely to be necessary (Gropp and Snir, 2013), potentially deploying new tools, or even the same tools (MPI, OpenMP), to deliver entirely new algorithmic constructs such as thread pools and task-based parallelism (e.g. Perez et al., 2008). Weather and climate models will undoubtedly target such systems, so their disruptive influence is likely to necessitate the wider use of these new programming models.

3.1 Domain-specific languages

As discussed above, much parallelisation in weather and climate codes is delivered by individuals constructing or modifying code-using libraries (such as MPI) and directives (such as OpenMP); this is despite decades of research into parallelising tools. The dominant reason for this is that humans can inspect an algorithm, and with knowledge of the hardware, further exploit domain-specific knowledge to reason about how to improve performance, but the domain-specific knowledge is not available to generic parallelisation tools. The result is that humans currently deliver science code that performs much better than that produced by generic tools (although this is not the case with compilers, which generally produce vectorised assembly code which performs better than that produced by humans)

So what are the domain-specific characteristics which are being exploited by the humans? Not surprisingly, they depend on the domain – the narrower the domain, the more a priori knowledge is available. Some of the earliest examples of this approach include the use of kinetic preprocessors for atmospheric chemistry. These take reaction specifications and generate code for specific processors (e.g. Linford et al., 2009).

Extending these ideas more widely means taking advantage of the knowledge that weather and climate models typically employ finite-difference, finite-volume, or finite-element discretisations of the underlying equations. As a consequence,

• operations are performed over a mesh,

• the same operations are typically performed at each mesh point, volume, or element,

• many operations are independent of each other, allowing for data parallelism,

• operations at a mesh point, volume, or element can either be computed locally or depend on neighbouring elements (leading to nearby neighbour – stencil – communication patterns), and

• reductions are required for convergence and/or conservation (leading to requirements for global communication).

In addition, the meshes themselves typically

• have fixed topologies,

• are structured or semi-structured (quasi-uniform),

• have resolution in the vertical which is much smaller than the horizontal (leading to 2-D + 1-D meshes rather than full 3-D meshes), and

• have data structured or unstructured in the horizontal and structured in the vertical.

We can further note that

• dynamical cores are mostly independent in the vertical, and

• physics (parameterisations) are mostly independent in the horizontal.

These characteristics fundamentally affect how to efficiently parallelise a code. For example, one would parallelise a 3-D finite-element model of an engine, which would typically be unstructured in all dimensions, very differently to a finite-element model of an atmosphere. In the latter case optimisations can make use of the knowledge that one direction – the vertical – is typically more structured and has lower dimensionality.

Domain-specific languages (DSLs) provide constructs that are suitable for the domain in question, rather than being general purpose. This approach can lead to a more intuitive (higher-level) interface for scientific programmers, since they will be dealing with a much simpler, more constrained language (in comparison with general-purpose languages). It can also lead to a better separation of concerns between the algorithmic description (using the DSL) and optimisation. For example, weather and climate codes in the form described above have two main communication patterns: nearest neighbour (stencil) accesses and global reductions. Therefore, the optimisation scheme associated with such a DSL need only understand and reason about these two forms of communication pattern.

DSLs are not necessarily designed for entire communities; they can be targeted at particular requirements within communities. For example, there are two exemplars of the DSL approach currently in use by the weather and climate community.

• GridTools (Fuhrer et al., 2014a), which was initially for the COSMO consortium. The GridTools project includes a wide range of activities, some of which are discussed further below, but the DSL component of the work has evolved from an earlier DSL (Fuhrer et al., 2014b and Leutwyler et al., 2016).

• PSyclone, which was for a new UK model – LFRic. PSyclone has evolved from the Gung Ho dynamical core project (Ford et al., 2013).

Each has been expressly developed to separate the science specification from the hardware optimisation so that the science writer only needs to write the science code, with separate implementation and optimisation for parallelism. This allows the science to be written once, with (potentially multiple) implementations tailored to particular architectures in such as way that they do not intrude on the maintainability (and evolution) of the science code. Both groups are already exploiting this potential to separately target different architectures (see e.g. the discussion in Sect. 3.2 below).

Both GridTools and PSyclone make use of domain-specific embedded languages (DSELs); i.e. their languages are embedded within an existing language, C$++$ and Fortran respectively. In both cases this choice was made as it was felt that the combination of available tools and community familiarity bought important benefits, but clearly different choices as to those benefits were made.

They also differ in their original target domains. The origin of GridTools was a DSL designed for a finite-difference regular latitude–longitude local area model, whereas PSyclone was designed for a finite-element model with an irregular mesh. However, the GridTools formalism is being extended to support irregular finite-difference models and PSyclone is being extended to support regular finite-difference models.

The approach to optimisation in the two systems is also quite different. In GridTools a particular back end (optimised by an HPC expert) is chosen and the science description is mapped to that particular back end via C$++$ templates. In PSyclone, an HPC expert specifies a recipe of optimisations which are then used by Psyclone to process the DSL input to produce code optimised for a particular architecture.

Despite these differences, both approaches rely on HPC expertise to deliver architecture-dependent performance. The DSL approach expedites the separation of concerns in such a way that computational experts can concentrate on the optimisation for given hardware, allowing different individuals to concentrate on different parts of the problem (e.g. Schulthess, 2015).

The use of DSLs has also expedited the involvement of vendors who usually have insight into how best to exploit their hardware. However, their ongoing deep engagement will be dependent on the size of the market; clearly for any given DSL, the more widely it is used, the more incentive there is for vendor engagement. Both projects recognise the importance of “growing their market” to maintain this incentive, but this is a complex issue (see Sect. 4.2).

3.2 Weather and climate dwarfs

As already noted, many components within weather and climate models have very flat computational profiles, but one approach to dealing with this is to identify “mini-apps” (analogous to the so-called Berkely dwarfs; Asanovic et al., 2006), which are exemplars of patterns in the code (e.g. common communication behaviours), so that performance issues can be investigated (or predicted) by using a small code base, typically a fraction of the size of the larger code (Heroux et al., 2009). When combined with a DSL approach, lessons learned from investigating mini-apps could be exploited more widely within the codes with minimal human intervention.

The ESCAPE (Energy-efficient Scalable Algorithms for Weather Prediction at Exascale; http://www.hpc-escape.eu, last access: 29 April 2018) project is investigating this approach for weather and climate with a specific focus on energy efficiency and emerging processor technology. ESCAPE initially defined dwarf categories (Table 1) which address key aspects of such models by including key components of the ECMWF spectral model IFS, its semi-implicit time-stepping formulation, and the semi-Lagrangian advection scheme. These allowed for an investigation of aspects of both global and limited-area prediction in the context of wave, ocean, and land surface components.

The IFS model takes one approach to the discretisation of prognostic partial differential equations and their evolution in time defined on a grid. However, alternative discretisations are also realised with a finite-volume approach, which added dwarfs related to the advection scheme (Kühnlein and Smolarkiewicz, 2017) and a 3-D grid-point-based elliptic solver (Smolarkiewicz et al., 2016), explored with both structured and unstructured grids. Two further dwarfs have also been chosen as representatives of the (many) physics parameterisations in models – a cloud microphysics and radiation scheme (Mengaldo, 2014; Geleyn et al., 2017). Different dwarfs would probably be necessary for other dynamics formulations such as those introduced in Sect. 2.2.

The computational performance challenges are fundamentally different between dwarfs just like their functional separation within an Earth system model. For example, spectral transforms are very memory and communication bandwidth intensive. Three-dimensional elliptic solvers are compute intensive and affected to a greater extent by communication latency. The semi-Lagrangian advection scheme is also communication intensive and depending on implementation may suffer from load-balancing issues. All three of these dwarfs have limited scalability because they involve data transfer and computations of global or regional fields, which is opposite to the nearest-neighbour principle. By contrast, parameterisations for cloud microphysics and radiation are very compute intensive but scalable, as they are applied to individual vertical columns that can be distributed across both threads and tasks in chunks aligned with the available memory. Hence, performance enhancement options need to be explored differently for each dwarf based on the range of available processor types and considering that all these processes need to be present in an Earth system model.

Table 1The ESCAPE dwarfs – a set of mini-apps providing a sample of characteristic behaviours of weather and climate model subcomponents. The “D” category dwarfs are related to the dynamics, the “I” dwarf is general purpose, and the “P” dwarfs represent physics-related issues.

Download Print Version | Download XLSX

For each dwarf, ESCAPE targets performance on existing and emerging processor technologies (specifically Intel Xeon, Xeon Phi, and NVIDIA GPGPU and a novel technique employing optical interferometry particularly suitable for Fourier transforms), but it is also targeting programming methodologies. Of necessity this involves comparing modified and existing MPI and OpenMP implementations with new OpenACC- and DSL-based implementations (the latter using GridTools).

3.3 Data models

The earliest weather and climate models discretised equations using a grid with one data point per cell. This led to easy data layouts, with data contiguous in memory, and relatively straightforward algorithms to find vectorisation. Modern grids and meshes are more complicated, not always structured, and often with multiple points per cell. Finite-element methods add additional complexity: instead of points representing values associated with each cell, basis functions are defined within each cell, and field “values” are established by parameters per cell known as “degrees of freedom (DOFs)”, which for each basis function establish how the field varies within the cell. These DOFs are then associated with mesh entities such as volumes, surfaces, edges, and vertices.

Figure 4Object orientation and the Grid classes supported in Atlas. Specific model grids can be constructed as instances of one of these classes. For example, the current ECMWF high-resolution (9 km) octahedral reduced Gaussian grid is an instance of the StructuredGrid, but can be referred to by the name Ø1280, which it inherits from the Grid class.

Download

There is no established convention or standard for layout in memory of these more complicated grids and meshes or how the various cell entities are associated with each cell. The current best practice is tightly coupled to knowledge of the grid being used; the ECMWF IFS model, for example, has a data structure designed specifically for the reduced Gaussian grids which arise from exploiting the spectral transforms used in the dynamical core (Hortal and Simmons, 1991). But the spectral transform method may no longer be attractive on future architectures due to its inherent all-to-all global communication pattern, and this is recognised in the exploration of alternative approaches via the extra dwarfs discussed above. Some of these other approaches rely on nearest-neighbour communication patterns and may even be applied to unstructured grids, leading to the requirement for different data structures to get maximum performance, potentially with different variants depending on the target hardware!

The Atlas library (Deconinck et al., 2016, 2017) is an attempt to address these issues. It provides a general framework for representing model data structures on different grids and meshes relevant to NWP and climate simulation. It is an open-source object-oriented C$++$ library with a modern Fortran interface designed to be maintainable and non-intrusively adapt to emerging hardware technologies and programming models. This is achieved by a design with a clear separation of concerns that, for example, provides a clean interface for exploitation by a DSL and allows multiple back ends for actual field data storage: a CPU back end and a GPU back end with the capability to synchronise with each other seamlessly.

The Atlas library provides a toolkit of functionalities ranging from low-level field containers to high-level mathematical operators and remapping operators. One example of the Atlas object-orientated design is the way that the grid structures are designed and exploited. We have already noted in Sect. 3.1 that weather and climate models have common properties that can be exploited, and grid type is one such common property, although the grids used vary, with regular grids, reduced Gaussian, icosahedral, and cubed-sphere grids all in common use. As a consequence, Atlas provides a classification of grid objects and supports their configuration for different model variants at different resolution via simple resolution parameters. Figure 4 shows how such grid objects can be constructed as specialisations of common characteristics with resolution parameters. More specific grids may be added non-intrusively by exploiting the object-orientated design.

Figure 5Conceptual workflow on how to use Atlas with most applications.

Download

Figure 5 sketches how these Grid objects can be used in a workflow to configure how memory structures can be set up using Atlas. A MeshGenerator can, using a given Grid, generate a parallel distributed Mesh. A Mesh partition provides an unstructured ordering of nodes that connected through lines (also called edges) form two-dimensional elements (triangles or quadrilaterals). Connectivity tables essentially link elements, edges, and nodes. Atlas further provides functionality to grow each Mesh partition with a halo so that parallelisation strategies requiring halo exchanges can be formulated. Using a Grid or a Mesh, a FunctionSpace can be created, which sets up how fields are discretised and parallelised.

A NodeColumns FunctionSpace is available to describe for every node of a horizontal mesh a contiguous-in-memory structured vertical column. A Spectral FunctionSpace describes a field with spectral coefficients. Additional FunctionSpace types can be created to describe discretisation for continuous or discontinuous spectral element methods (e.g. Marras et al., 2016).

Within ESCAPE, Atlas is being used as the foundation for the evolution of the NWP and climate dwarfs. As a consequence Atlas acquired improved mesh-coarsening strategies, a GPU-capable Field data storage, and support for regional grids with projections and checkerboard domain decomposition. Clean interfaces in both Atlas and GridTools are being used to exploit Atlas memory storage to allow GridTools to compute stencil operations on unstructured meshes.

3.4 MPI libraries, coupling, and I/O

All methods of coarse-grained parallelisation which do not use global address spaces exploit the MPI paradigm, whether it is managing domain decomposition, coupling components, or exploiting I/O servers.

Exploiting MPI directly is difficult, and models typically exploit some sort of middleware library such as Generalised COMmunications (GCOM) used by the Met Office UM or YAXT (Yet Another eXchange Tool; https://doc.redmine.dkrz.de/yaxt/html/da/d3c/concept.html, last access: 29 April 2018) under testing for use in ICON. Such tools can simplify the usage of MPI and therefore make the performance potential of data layout and data independence accessible (see e.g. Fig. 6).

Figure 6The influence of using a library to manage MPI communications for domain decomposition: experiences using YAXT to replace the existing ECHAM usage of MPI in a T63L47 model. The YAXT optimised library (green) can lead to factors of 10 improvement in the time taken to do transpositions of data over simple implementations (blue) and factors of 2 over manual optimisation (red, which was effectively an aggregation of multiple transpositions into bigger messages). In this case the two main optimisations in YAXT were avoiding additional data copies by the usage of MPI data types and reordering of messages to avoid congestion, with the latter becoming more important with higher core counts.

Download

The complexity of data transpositions and MPI communications increases in the context of coupling, and these also need to be understood in terms of how best to obtain performance. Many groups are investing in coupling tools, such as the OASIS coupler widely exploited in Europe. Recognising that there is much to be learnt from these various activities, the community has held a series of workshops on coupling, which led to the identification of a range of parameters which influence coupling and can be used in a benchmarking strategy (Table 2).

With these characteristics in mind, the IS-ENES2 project established a programme of investigation into coupling performance by developing some stand-alone components suitable for coupling tests and carried out some initial investigations testing the impact of the number of cores per component, the grid sizes, the grid aspect ratios and orientations, and having different numbers of cores for different components.

As a proof of concept, these coupled test cases were run using five different coupling technologies in different configurations on three different platforms: Bullx at CINES in France, Cray XC40 at the UK Met Office, and the Broadwell partition of Marconi at CINECA in Italy. Thus far the work has demonstrated the utility and versatility of the benchmarking environment; future work will focus on learning lessons from comparisons of performance (Valcke et al., 2017).

Table 2A set of characteristics that influence coupler performance.

Download Print Version | Download XLSX

3.5 Tools to support porting, optimisation, and performance portability

Automatic code-parsing tools provide a route to code migration (e.g. Fu et al., 2016, as discussed above), but their efficiency depends on the internal models of how codes are and can be organised, in particular the presence or absence of pre-existing manual optimisations, which if not understood by the tool, can require significant manual effort to remove (before and/or after the use of the tool).

It is clear that achieving performance requires a comprehensive understanding of what existing model components or dwarfs are doing, and performance portability requires designing code that will still have performance on alternative architectures. With the complexity of modern codes, this can no longer be efficiently done with access to simple timing information, and more complicated information needs to be gathered and presented for analysis by developers. Such information is normally provided by performance tools.

Performance tools can provide the support and the freedom to the developer to view all the different levels of concurrency (from the hardware counters underpinning fine-grained concurrency to the load-balancing involved in coarse-grained concurrency). While it is possible to collect a vast amount of information by relinking entire codes to performance monitoring libraries, results are often enhanced by manually identifying hot spots and extracting additional trace information, and automatically generated call trees can be enhanced by additional instrumentation.

In particular, like automatic code migration tools, the tools which generate trace files, such as the EXTRAE library underpinning the generation of trace files for PARAVER or DIMEMAS, depend on an internal understanding of the programming models in use in the target codes. It is this understanding that allows them to automatically trap library calls and read hardware counters without code invasion, but it is precisely because of this dependence on generalised programming models that extra manual intervention can be necessary to capture domain-dependent information in trace files.

Figure 7Understanding model performance on existing machines generally requires instrumenting the model code to output “trace” information that can be analysed and visualised by a dedicated application. This information can also be used, when combined with machine information, to predict performance (issues) on alternative architectures. This figure depicts the use of PARAVER (Labarta et al., 2005) for collecting trace information and DIMEMAS (Gonzalez et al., 2011) for simulating performance on alternative architectures.

Download

Figure 8Example trace information generated with EXTRAE used to analyse high-resolution coupling in EC-Earth. This view shows MPI calls during the execution (one colour per call). Top: 512 MPI tasks for the IFS atmospheric model. Bottom: 128 MPI tasks for the NEMO ocean model. This fine-grained performance information yields information as to how to get better performance out of coarse-grained concurrency.

Download

The resulting trace files can be analysed directly to help improve performance on any given architecture and also be used to predict performance on alternative architectures (see Fig. 7). Any extra levels of invasive alteration of code add further levels of programming complexity, and even with intuitive analysis tools (see e.g. Fig. 8), considerable time and expertise is required to interpret the outputs. However, the regular use of such tools by DSL developers and those evaluating dwarfs is one route to maximising the benefit of these tools without the community needing to gather the relevant expertise.

3.6 Related work

Weather and climate computing is not alone in facing these problems. For example, the notion of DSLs as a solution has a tried and tested heritage – examples include the Kokkos array library (Edwards et al., 2012), which like GridTools uses C$++$ templates to provide an interface to distributed data which can support multiple hardware back ends, and from computational chemistry, sophisticated codes (Valiev et al., 2010) built on top of a toolkit (Nieplocha et al., 2006), which facilitates shared memory programming. Arguably DSLs are starting to be more prevalent because of the advent of better tooling for their development and because the code they generate can be better optimised by autotuning (Gropp and Snir, 2013). However, in our case, we still believe that human expertise has a primary role to play, but this expertise can be better deployed by supporting DSLs than entire code bases.

The mini-app or dwarf approach we describe above is also starting to become more widely prevalent, mainly because it provides a tractable route to both predict performance and develop and test performance optimisations for new hardware. Its applicability to weather and climate modelling follows directly from the fact that the source of the methodology was in disciplines for which similar flat multiphysics profiles were involved. The concept of attempting to represent a cross section of applications or application characteristics is becoming well recognised (e.g. Martineau et al., 2016), and others are investigating applicability to environmental models (e.g. Stone et al., 2012) and noting the importance of representativeness in the selection of the mini-apps.

Problems of establishing appropriate data structure, generating meshes, and simplifying the related interfaces and workflows are prevalent across all communities addressing the solutions to partial differential equations, whether on a sphere or not, with decades of relevant research for deployment on parallel computers. More recently, work on data structures has begun to address cache awareness and the use of space filling curves (e.g. Günther et al., 2006), and a wide range of mesh generators are available. The issue for weather and climate modelling is not so much what relevant work there is, but how best to establish domain-specific tooling which yields scientific and software development productivity across the community.

Further up the software stack, model couplers, frameworks, and I/O servers are all tools for which the domain specificity dominates, although many communities have analogous problems. Tools for evaluating code performance are common across communities, with important related work in the establishment of standards for trace file formats (e.g. the Open Trace Format; Knüpfer et al., 2006) and common notations for adding directives for extra performance-related information.

4.1 Requirements

Existing models are known to have high levels of software quality (Pipitone and Easterbrook, 2012), and maintaining quality models in the face of increasing scientific demands on what is simulated, ensemble size, and resolution will be crucial. It is these that define scientific quality – a model could be very high performance and portable, but not be suitable for the scientific objectives, and there is a spectrum between “not suitable” and “ideal”. In fact, some level of model quality is often compromised for performance; for example, resolution and complexity are routinely sacrificed to fit weather forecasts in a particular time window or to reach a requisite number of simulated years per day for climate.

It is not so well known that quality is also often compromised for productivity and portability. One of the reasons for the popularity of the Weather Research and Forecasting model (WRF; Skamarock et al., 2005) is that it is easy to modify and run, and this productivity–portability benefit outweighs the use of alternative models even when they are perceived to deliver better-quality weather simulations.

Until recently, given the dominance of X86 systems, performance has been relatively independent of platform; once compiler differences have been negotiated, getting the best performance has been mostly about understanding how to exploit processor decomposition given the interconnect and ignored by much of the modelling community, except by those with the most demanding applications (typically in weather predictions in tight wall-clock time windows or in long, complex, or high-resolution climate runs). Nonetheless, performance also involves the need to optimally fill machines (or execute large ensembles), so optimising performance can also mean addressing throughput rather than speed (Balaji et al., 2017), meaning that performance optimisation should not only be about good fast code, but also about evaluating speed for a variety of situations within a large potential parameter space (e.g. processor decomposition and placement).

Portability, which is the ability to migrate codes between architectures and achieve acceptable performance (with confidence that the code is delivering the same scientific results; e.g. Baker et al., 2015) is very important, whether to ensure that codes can be executed on many platforms by a wide community or to ensure that the codes will work on future architectures. In many groups it is seen as one of the most important considerations (e.g. Norman et al., 2017). However, as we have seen, obtaining ideal performance on multiple architectures with the same code is problematic, so performance portability is always a compromise. In fact, it has been said that “the real goal is probably better labelled maintainability than strict portability” (Norman et al., 2017, again).

The difficulty we now face is the fact that because of the oncoming heterogeneity, even vendors are suggesting that the community is facing a world in which one can only have two of the three factors of performance, portability, and productivity – a happy compromise is no longer possible. While this is again to some extent hyperbole, it is certainly true that much more effort needs to be devoted to achieving these simultaneously and that quite major changes in approach (e.g. the use of DSLs) are needed: so it is the amount of effort needed, the change in approaches required, and the timescales in play that lead to the “cliff edge chasm” motif for this paper.

4.2 Making progress

To make progress in the science, existing weather models are being extended to seasonal and decadal scales by running with additional components (ocean, sea, and land ice), and existing climate models are confronting the need to start from prescribed initial conditions and run at increased resolution. New models are being developed, but the existing timescales from development to production science serving large communities are currently of the order of a decade. All these models need to exploit all the avenues of concurrency outlined in Sect. 2.3 and at the same time support multiple architectures. How then can we make progress?

There are essentially four strategies to sustaining a modelling programme:

1. develop a model and all the supporting components from scratch,

2. build a particular piece or pieces of a model and utilise components and/or infrastructure from elsewhere,

3. utilise someone else's model, but tinker with some aspects of it, or

4. simply configure and run models developed elsewhere.

Which strategy is most appropriate depends primarily on the resources and skill sets available and the productivity possible. Only a handful of groups can afford the first strategy, and the further down the list one goes, the less control one has over the future trajectory that is possible for that modelling programme because of the external dependencies.

Risks range from not being able to move fast enough to meet science goals, to having key capability removed because the external group has chosen either to do something completely different or change key aspects of the interfaces upon which the model depends. In the worst case, having outsourced some capability, the internal group may no longer have the capability to replace that capacity. Mainly for these reasons, large (primarily) national modelling endeavours prefer to keep as much development as possible in house. However, looking forward, it is not obvious that even large national groups have the internal resources to both keep up a production line of incremental model improvements associated with meeting near-term scientific (and/or operational) requirements and identify and take the requisite steps necessary to develop codes which can hit quality + performance + portability + productivity requirements using next generation computing, particularly when the latter is unlikely to be achievable with small steps. This is recognised in some large institutions; for example, the Met Office is investing in OASIS, NEMO, and XIOS, all of which have replaced or will replace internal developments.

If the existing incremental development method will not provide more than modest improvements in performance at very low levels of computational efficiency, how then can progress be made? We would assert that the way forward has to involve the community making better use of the collective resources and more groups utilising components, particularly infrastructure components, developed elsewhere. It is such shared libraries and tools that we see as the “spans” that will need to be utilised to cross the “chasm” of model development ahead. That said, we fully recognise the importance of competitive evolution: the optimal number of libraries and tools is not one of each; two or three might be the right answer, but it should certainly not approach N, where N is the number of large modelling groups.

We have introduced a range of such tools, libraries, and methodologies for making progress such as the identification of dwarfs, the publication of data structures, and the elements necessary to establish a comparison of coupling technologies, and to evaluate performance. Within some of that work there is a clear aim to develop and specify interfaces that can be reused in different ways (for example, the interfaces which allow GridTools to interface directly with the Atlas data structures), but many of these (and similar) projects have still mainly been developed by small groups who have not had the time, mandate, or wherewithal to publicly promulgate requests for requirements or specifications and then build to documented interfaces.

Progressing at the community level will require improved methods to allow the community to discuss, specify, design, develop, maintain, and document the necessary libraries and tools. The weather and climate community does not have a great track record at sharing such tools, although in recent years necessity has begun to influence practice, with tools such as OASIS becoming more prominent in more places. One of the reasons for this lack of sharing is in part the lack of a commonly deployed structured approach to sharing, specifically one that maximises the delivery of requirements while minimising risk of future technical burden – this is the sort of approach that has delivered the MPI libraries upon which nearly all of high-performance computing depends. While a fully fledged standards track is probably beyond the will of the community at this point, it is certainly possible for the community to investigate and then engage in more collaborative working.

Such steps would begin with institutions and modelling groups recognising the scale of the problem ahead and recognising the following:

• business as usual, consisting of modest incremental steps, is unlikely to deliver the requisite next generation models;

• they do not have enough internal resources to take the leap to the next generation alone; and, most importantly,

• there are library or tool projects which they can exploit and through which they can collaborate rather than compete, some of which may be from outside the traditional communities of collaborators.

Institutions which are most likely to be able to seize these collaborative opportunities are most likely to share certain characteristics. They will do the following:

• understand the issue fully at the management level, the science level, and in the infrastructure teams;

• be able to reward individuals for innovation in and/or contributions to external projects;

• recognise the benefit of external scrutiny and avoid situations in which they are forever arguing that their code is “not ready to share”;

• have the courage to stop existing activities and pick up and use or integrate third party libraries and tools;

• have the ability to recognise the cost–benefit trade-off between “doing it themselves” and contributing intellectually and financially to third party solutions; and

• be ready to apply more sophisticated and complex software engineering techniques and encourage more computational science research.

Projects which would be suitable for collaborative working will share some common characteristics. They will do the following:

• be open source and have an open development process,

• have clear goals, scope, and where appropriate, deliver stable software interfaces,

• have a mechanism to understand and respond to the timescales of collaborators (that is, some sort of governance mechanism which assimilates and responds to requirements),

• potentially be able to accumulate and spend funds to provide user support, training, and documentation,

• be not initially disruptive of existing solutions, and ideally

• engage both the scientific community and vendors (compare with MPI for which vendor implementations are often key to enhanced MPI performance and reliability).

The advantages of good responsiveness to requirements and/or vendor engagement cannot be underestimated – where potential users see a project which is responsive to problems or one with good vendor support they are more likely to commit to dependencies on that software. Where they see dependencies on software with inadequate or obscure support, commitment is more difficult. There are examples of such governance for libraries; for example, both ESMF and OASIS (Sect. 2.2) have established credible community governance. However, the issue is exacerbated for toolchain software dependencies – if the software dependency is a library, one can at least contemplate a substitution; if it were, for example, a DSL, then the room to manoeuvre is far more limited. This probably means that the use of DSLs (and other critical toolchain infrastructure) is most likely to spread slowly based on closer collaborations rather than loosely coupled relationships until there is a visible cross-institutional developer community and/or they become commercially supported. Even then, the scope of such dependency rather depends on the scope of the DSL itself, and it may be more appropriate for the community to aspire to interoperability and exploitation of libraries within more narrowly focused DSLs which do not themselves support large communities.