Geosci. Model Dev., 12, 3135–3148, 2019 https://doi.org/10.5194/gmd-12-3135-2019 © Author(s) 2019. This work is distributed under the Creative Commons Attribution 4.0 License. How to use mixed precision in ocean models: exploring a potential reduction of numerical precision in NEMO 4.0 and ROMS 3.6 Oriol Tintó Prims1,2 , Mario C. Acosta1 , Andrew M. Moore3 , Miguel Castrillo1 , Kim Serradell1 , Ana Cortés2 , and Francisco J. Doblas-Reyes1,4 1 EarthSciences Department, Barcelona Supercomputing Center – Centro Nacional de Supercomputación, Barcelona, Spain 2 HPCA4SE research group, Computer Architecture and Operating Systems Department, Universitat Autònoma de Barcelona, Bellaterra, Spain 3 Ocean Sciences Department, University of California, Santa Cruz, CA, USA 4 ICREA, Barcelona, Spain Correspondence: Oriol Tintó Prims (oriol.tinto@bsc.es) Received: 21 January 2019 – Discussion started: 15 February 2019 Revised: 8 June 2019 – Accepted: 24 June 2019 – Published: 24 July 2019 Abstract. Mixed-precision approaches can provide substan- sion. The number of possible combinations is prohibitively tial speed-ups for both computing- and memory-bound codes large and therefore impossible to explore. The alternative of with little effort. Most scientific codes have overengineered performing a screening of the variables individually can pro- the numerical precision, leading to a situation in which mod- vide certain insight about the required precision of variables, els are using more resources than required without know- but, on the other hand, other complex interactions that in- ing where they are required and where they are not. Conse- volve several variables may remain hidden. Instead, we use quently, it is possible to improve computational performance a divide-and-conquer algorithm that identifies the parts that by establishing a more appropriate choice of precision. The require high precision and establishes a set of variables that only input that is needed is a method to determine which real can handle reduced precision. This method has been tested variables can be represented with fewer bits without affect- using two state-of-the-art ocean models, the Nucleus for Eu- ing the accuracy of the results. This paper presents a novel ropean Modelling of the Ocean (NEMO) and the Regional method that enables modern and legacy codes to benefit from Ocean Modeling System (ROMS), with very promising re- a reduction of the precision of certain variables without sac- sults. Obtaining this information is crucial to build an actual rificing accuracy. It consists of a simple idea: we reduce the mixed-precision version of the code in the next phase that precision of a group of variables and measure how it af- will bring the promised performance benefits. fects the outputs. Then we can evaluate the level of precision that they truly need. Modifying and recompiling the code for each case that has to be evaluated would require a prohibitive amount of effort. Instead, the method presented in this paper 1 Introduction relies on the use of a tool called a reduced-precision emula- tor (RPE) that can significantly streamline the process. Using Global warming and climate change are a great challenge the RPE and a list of parameters containing the precisions for human kind, and given the social (e.g., Kniveton et al., that will be used for each real variable in the code, it is possi- 2012, regarding climate refugees), economic (e.g., White- ble within a single binary to emulate the effect on the outputs man and Hop, 2013, regarding the trillion dollar problem) of a specific choice of precision. When we are able to emu- and environmental threat (e.g., Bellard et al., 2012, regard- late the effects of reduced precision, we can proceed with the ing mass extinctions) that they pose, any effort to understand design of the tests that will give us knowledge of the sensi- and fight them falls short. Better knowledge and greater ca- tivity of the model variables regarding their numerical preci- pacity to forecast how the climate will evolve can be a game- changing achievement, since it could help to justify ambi- Published by Copernicus Publications on behalf of the European Geosciences Union. 3136 O. Tintó Prims et al.: How to use mixed precision in ocean models: NEMO 4.0 and ROMS 3.6 tious policies that adapt to future scenarios (Oreskes et al., halved. For this reason, the use of smaller representations can 2010). The Earth system can be seen as an amalgamation of now provide performance benefits that justify the effort of many parts: the atmosphere, the hydrosphere, the cryosphere, optimizing the numerical precision (Baboulin et al., 2009), the land surface and the biosphere. All these elements are and while this is true with the actual hardware, the expected extremely rich in phenomena and are open and inter-related. potential is even bigger in future architectures that will in- Fluxes of mass, heat and momentum interchange in ways that clude not only 64- and 32-bit but also 16-bit arithmetic. are virtually endless, some of which are poorly understood We are now in a situation in which ESMs, as computer or unknown. The magnitude and complexity of these sys- codes of other domains, need to use computational resources tems make it difficult for scientists to observe and understand efficiently and in which mixed-precision approaches are them fully. For this reason, the birth of computational science emerging as a potential solution to help improve efficiency. was a turning point, leading to the development of Earth sys- The main risk of reducing precision is falling short, since tem models (ESMs) that allowed for the execution of exper- using less precision than needed can lead to numerical errors iments that were impossible until then. ESMs, despite be- that make model results inaccurate or simply wrong. The pre- ing incomplete, inaccurate and uncertain, have been a frame- cision required depends on many factors, so what is needed is work in which it is possible to build upon knowledge and a method to identify which variables can effectively use less have become crucial tools (Hegerl and Zwiers, 2011). Since precision without compromising the quality of the simula- their inception, the capability to mimic the climate system tions and which ones cannot. If the precision can be reduced, has increased and with it the capacity to perform useful fore- in many situations it is because the precision has been ov- casts (Bauer et al., 2015). The main developments that have erengineered. One example is the precision used to represent led to this improvement in model skill are the enhancement the input data provided to ESMs. The precision of the physi- of physical parameterizations, the use of ensemble forecasts, cal observations is limited by the instruments used to collect the improvement of the model initialization and increases in them. In the case of sea surface temperature measurements resolution (Bauer et al., 2015). Most of these contribute to a from Earth-orbiting satellites, this precision is a few tenths higher computational cost (see Randall et al., 2007). For this of a degree; 64-bit representations are often used when 16 reason, these developments are only possible with an increas- bits could potentially be enough. ing availability of computing power, a situation that will con- Recent work has demonstrated the potential benefits that tinue in the future. The motivation to make models efficient is mixed-precision approaches can provide to many different twofold. Firstly, developments that are considered crucial to kinds of codes, since it is possible to achieve substantial increase the skill of the models require more computational speed-ups for both computing- and memory-bound codes re- power. This is not just a matter of having a larger machine quiring little effort with respect to the code (Baboulin et al., since some of the issues that emerge are not trivial and re- 2009). The spectrum of studies goes from explicit code ma- quire additional developments (Dennis and Loft, 2011). Sec- nipulation in very specific algorithms to automatic modifi- ondly, the huge investment in computational resources that cation of binaries of any kind of code. Some studies have is necessary to perform simulations with ESMs implies that focused on the use and development of mixed-precision al- investing time in optimizing them will be of value. gorithms to obtain performance benefits without compromis- One research field that has gained momentum in recent ing the accuracy of the results (Baboulin et al., 2009 and years and that can improve model performance is the use of Haidar et al., 2018). There are also several automatic mixed- mixed-precision algorithms. Until not so long ago, the speed precision exploration tools (Graillat et al., 2016; Lam et al., of most computations was constrained by how fast a CPU 2012) that have been mainly tested on small benchmarks, could perform operations, with the speed of the memory be- usually C++ codes. These kinds of studies inspired Earth sci- ing fast enough to provide more data than the processor could ence groups working with ESMs that are willing to improve process. In addition, CPUs were designed in a way that they their computational performance to make bigger and more could virtually perform operations at the same speed no mat- ambitious experiments possible (Váňa et al., 2017; Düben ter whether they were operating with 32- or 64-bit floating- et al., 2014, 2017; Düben et al., 2017; Thornes, 2016). In- point representations. Therefore, the only benefit of using spired by previous work, we propose a method that automat- less precision was to reduce the memory requirements, and ically explores the precision required for the real variables the computational performance was not so much a motiva- used in state-of-the-art ESMs. tion. Mainly two factors have changed that scenario. First, The method emerged while trying to explore how to CPU speed increased at a faster rate than memory speed, achieve simulations as similar as possible to standard double- meaning that at some point many codes that were CPU- precision simulations while reducing the precision of some of bound before would become memory-bound, with the mem- the variables used. Our work extends the aforementioned re- ory bandwidth being insufficient to feed all the data that the search to achieve mixed-precision implementations for full- processor can process. Second, vector operations doubled the scale models. To do so we rely on a reduced-precision emu- number of floating-point operations per cycle that could be lator (RPE) (Dawson and Düben, 2017), which mimics the performed when the number of bits of the representation is effects of using an arbitrary number of significand bits to Geosci. Model Dev., 12, 3135–3148, 2019 www.geosci-model-dev.net/12/3135/2019/ O. Tintó Prims et al.: How to use mixed precision in ocean models: NEMO 4.0 and ROMS 3.6 3137 represent the real variables in a code and measure the im- tained at high precision except the one that is being tested. pact of a specific reduced-precision configuration in the out- While it is true that this approach can provide insight about put produced by the model. Minimizing user intervention by the precision needed by a particular variable, it may not automating all the tedious intermediate processes allows for reveal issues regarding more complex interactions between an analysis of models that would have otherwise required too more than one variable. Additionally, building a complete set much manpower. Although the tool is very convenient for of variables from the tests performed on individual variables exploring the impact of the bits used to represent the man- is not trivial at all; i.e., a set that consists of variables that, as tissa of the floating-point numbers, the effect of changing the individuals, can use lower precision may not behave as such number of bits devoted to representing the exponent is not when combined and produce inaccurate results. explored. Another alternative is a divide-and-conquer algorithm that To test the methodology, this work includes two case stud- identifies the sections of the code that cannot handle reduced ies. These cases correspond to two different ocean models precision and builds a complete set of variables that can. that are widely used worldwide: the Nucleus for European Starting from a variable set that contains all the real variables Modelling of the Ocean (NEMO) and the Regional Ocean that we want to analyze, the approach consists of evaluating Modeling System (ROMS). With these models we demon- what happens when the set uses reduced precision. If the re- strate how the methodology can be used with different appli- sults become inaccurate, we proceed to split the set in two cations, thus demonstrating its potential. parts that are evaluated individually. The process is recursive, and ideally a binary search is performed until the sets that are evaluated contain only a single variable. If a set containing a 2 Methods single variable ends up being inaccurate, this variable is kept in high precision in the preceding sets when they are reeval- In this section we will demonstrate how we can establish uated. The advantage of this approach is that it is cheaper to which real variables written in a Fortran code can effectively find all the variables that can individually compromise the use less precision than the de facto 64 bits. The reader will results, and in addition it is easier to rebuild a complete set find an explanation as to why and how we developed this by assimilating the results of the simulations of subsets that method, with the specific steps of the methodology detailed gave an accurate result. below. The basic idea behind the method is to perform simu- Nevertheless, there are some elements that prevent this lations with a Fortran model using a custom set of precisions approach from working properly. The nonlinearity of most and directly assess the outputs to see if the results are accu- Earth science codes implies that the differences between two rate enough. To do so we use the RPE tool, a Fortran library simulations performed using different numerical precisions that allows us to simulate the results of a floating-point oper- will not be constant. In many cases two accurate subsets re- ation given that the variables are using a specific number of sulted in an inaccurate set when combined. To increase the significand bits. This can be integrated into an actual code to confidence in the results, we propose reevaluating the sets mimic the possible consequences of using reduced precision whose results are accurate with different initial conditions. in certain variables. The advantage of the tool is the flexibil- This method is similar to the ensemble simulation used in ity that it offers once implemented into the code, allowing us ESMs (Palmer et al., 2005), which tries to assess the uncer- to easily test any given combination of precisions. The main tainty of the simulation outcomes by taking into consider- drawback is the considerable overhead added to the simula- ation the uncertainty in the model inputs. The method re- tions, increasing its cost. mains the same but adding an extra simulation initialized The objective of the method is to find a set of precisions with different initial conditions for the sets that show accu- that minimizes the numerical resources used while keeping rate results. According to our results, an ensemble of only the accuracy of the results. A set of precisions is a specific two members was required to solve most of the issues related combination of the precisions assigned to each variable, with to combinations of accurate subsets resulting in inaccurate 52n being the number of possible sets, where n is the number sets. of real variables used in the code and 52 the number of bits The steps of the methodology, which are discussed in the used to describe the significand in a double-precision repre- next subsections, consist of the following: sentation. This number makes it prohibitively expensive to explore all combinations. For any real-world code, this holds – implementing the emulator into the code, completing all true even when not considering all possible values between 0 the necessary actions to obtain a code that uses the em- and 52 but just considering double precision (64 bits), single ulator whereby it is possible to select the precision of precision (32 bits) and half-precision (16 bits), wherein the each real variable through a list of parameters; number of possible sets is still 3n . A feasible alternative is to perform a screening of the vari- ables individually. The idea is to simply perform n simula- – establishing a test that will determine if the results of a tions to observe what happens when all the variables are re- simulation are accurate enough; and www.geosci-model-dev.net/12/3135/2019/ Geosci. Model Dev., 12, 3135–3148, 2019 3138 O. Tintó Prims et al.: How to use mixed precision in ocean models: NEMO 4.0 and ROMS 3.6 – performing a precision analysis by launching the neces- – Although it is possible to adapt the RPE library to in- sary tests to obtain a set of variables that can effectively clude intrinsic functions (i.e., max, min) that can use use reduced precision to pass the accuracy test. rpe_var variables as arguments, there is still a problem when mixing rpe variables and real variables. The prob- 2.1 Implementing the emulator lem can be overcome by converting all the variables to the same type. The RPE is a Fortran library that allows the user to emulate floating-point arithmetic using a specific number of signifi- – When there is a call to an external library (i.e., NetCDF, cand bits. It has the capability to emulate the use of arbitrary message-passing interface – MPI), the arguments can- reduced floating-point precision. Its development was mo- not be rpe_var and must be reals. tivated by the need to explore mixed-precision approaches in large numerical models that demand more computational – Read and write statements that expect a real variable power, with weather and climate fields in mind, making the cannot deal with rpe_var variables. tool highly suitable for our purpose. The emulator is open source, can be accessed through GitHub and has documen- tation available, including a reference paper by Dawson and – Several other minor issues (pointer assignations, type Düben (2017) with more detailed information. conversions, the modules used, etc.) are also possible. Although the use of the emulator facilitates the testing In codes with hundreds to thousands of variables this task of different precision configurations without recompiling, in can represent months or years of work, and for this reason it large codes like ESMs the implementation of the emulator is worthwhile to automate the whole process. When all the can carry more work than expected. The length of the code, issues are solved the model should be able to compile and the large development time, the quantity of different devel- run. opers and the lack of a rigid style guide can result in a large number of exceptions that make it harder to fully automate 2.1.2 Selecting the precision the emulator implementation, requiring the user to solve the emerging issues. To specify at runtime the precision of each individual vari- Our implementation of the emulator has two different able, the method that we use is to create a new Fortran mod- parts: ule that includes an array of integers containing the precision 1. replacing the declaration of the real variables with the value of each one of the real variables of the model. The val- custom type rpe_var and ues of this array will be read and assigned at the beginning of the simulation. After implementing the emulator into a code, 2. introducing a method of selecting the precision of the one should be able to launch simulations individually spec- variables without requiring recompilation of the code. ifying the number of significand bits used for each variable and obtain outputs, having completed the most arduous part 2.1.1 Replace variable declarations of the proposed methodology. To use the emulator, the user has to replace the declarations 2.2 Designing the accuracy tests of the real variables with the custom type defined in the em- ulator library (see Dawson and Düben, 2017). Even though Once we are able to launch simulations, we must define how the idea is quite simple, the practical process in a complex we will verify the results, the kind of experiment that we want and large state-of-the-art ESM can present several minor is- to perform and a true–false test to perform on the outputs. sues that can add up to a considerable amount of work. This To define a test to verify results, we must define a function is a list of some of the specific issues that can be found when that can be applied to simulation output and determine if the implementing the emulator. outputs are correct. – All the real variables that were initialized at declaration To give an example, consider a simulation whose only out- time need to be initialized, providing a derived type of put is a single scalar value; then, we can consider a given test variable instead of a real, which requires modification to be accurate if the output and reference match to a certain of all these declarations. number of significant figures. Using the value of π as a ref- erence, we can consider a given simulation to be accurate – When a hard-coded real is used as a routine argument enough if the difference between the reference and the value in which the routine is expecting an rpe_var variable, obtained is smaller than a given threshold; i.e., the results it is necessary to cast this hard-coded variable into an coincide up to a specific number of significant figures. For rpe_var type of variable example, a result accurate to 6 significant figures would pass (i.e., “call routine(var, 1.0, var2, 0.1)” has to be if the required condition was to coincide with 4 significant “call routine(var, rpe_var(1.0), var2, rpe_var(0.1))”). figures, but it would fail if it had to coincide with 10. Geosci. Model Dev., 12, 3135–3148, 2019 www.geosci-model-dev.net/12/3135/2019/ O. Tintó Prims et al.: How to use mixed precision in ocean models: NEMO 4.0 and ROMS 3.6 3139 2.3 Performing the analysis that any given set has to follow in order to learn which ele- ments of the set can use reduced precision. We propose a recursive divide-and-conquer algorithm with a few slight modifications. For a given set of variables, we 3 Study cases generate a list of parameters that set these variables to use re- In this section we will present a proof of concept of duced precision and we launch a simulation. When the sim- the method using two state-of-the-art models, NEMO and ulation is completed, we proceed to apply the accuracy test ROMS. For each one of the two models we carried out two described in Sect. 2.2. If the simulation with the specific set different experiments. With NEMO, we performed two anal- passes the test we consider it safe to reduce the precision of yses using different accuracy tests, and for both cases the tar- this set. If this is not the case, the purpose of the algorithm get is to find which variables can use single precision (23-bit is to identify which part of the set is responsible for the inac- significand). With ROMS, we performed two analyses with curacy: we proceed by subdividing the set and evaluating its the same accuracy tests but having two different reduced- parts separately until we have identified one of the variables precision targets: single precision and half-precision (23- and that needs to preserve higher precision. The sets that yield 10-bit significand). The section is divided into three parts, inaccurate results require information from the subsets after with one part for each model, and finally a discussion of the these have been evaluated to be modified and reevaluated. results. This initial approach has some drawbacks. On the one hand, the effect of reducing the precision of a set of variables 3.1 NEMO can not be directly deduced from what happens to its individ- ual parts. The simplest case in which this can happen is when NEMO (Nucleus for European Modelling of the Ocean) is we deal with two variables that appear in the same arithmetic a state-of-the-art modeling framework of ocean-related en- operation: using the same logic as the actual processors, the gines. The physical core engines for the different parts of the emulator performs intermediate operations using the largest model are the OPA engine for the ocean, LIM for the sea ice precision between the variables involved, so it may happen and TOP–PISCES for biogeochemistry (Madec, 2008; Rous- that having any of the two variables in higher precision will set et al., 2015; Aumont et al., 2015). The range of applica- give an accurate enough result but not when both of the vari- tions includes oceanographic research, operational oceanog- ables use reduced precision. raphy, seasonal forecast and (paleo)climate studies. To illustrate this, let us consider the following example: Previous performance analyses of NEMO have shown that having two variables x and y, with x = 26 and y = 2−6 . If the most time-consuming routine was not even responsible we compute the sum z = x + y in real-number arithmetic, for 20 % of the total computation time (Tintó Prims et al., the result is z = 64.015625. Both x and y can be perfectly 2018). For this reason, any effort to improve the computa- represented using a 10-bit significand, but it is not the case tional performance of the model cannot target a single region for z, which requires more numerical precision. Following of the code but something that has to be applied along all the the processor logic, if either x or y uses a 52-bit significand, sections. the computation of z will be done using 52 bits, leading to As explained in Sect. 1, previous publications have a correct result, but if both variables are using 10 bits, the demonstrated the positive impact that the use of mixed- computation will be done using 10 bits and will yield a wrong precision approaches can have on the performance of scien- result (z = 26 ). tific models. Previous experiences in reducing the working On the other hand, numerical error in most algorithms has precision of NEMO from 64 to 32 bits demonstrated a sig- a stochastic component that, combined with the nonlinearity nificant change in the results (see Fig. B1), indicating that of these kind of models, can mean that a specific set of vari- blindly reducing the precision in the entire code was not an ables gives accurate results under certain conditions but gives option. inaccurate results under different conditions. To make the outcome of this work relevant for the model- For these two reasons we added an extra test into our work- ing community, the analysis has been performed using ver- flow that consists of assessing the results with different initial sion 4.0b of the code, which at the time of writing is the lat- conditions whenever the first evaluation results in a positive est version available. The configuration used was an ocean- outcome. While we cannot be sure that performing only a only simulation using the ORCA2 grid (about 220 km reso- single reevaluation will be enough, it can be sufficient for the lution near the Equator), and the objective of the analysis is most sensitive cases. to identify which set of variables can effectively use 32-bit Another slight modification that we added to the initial ap- floating-point representations instead of the 64-bit standard proach is the definition of stricter thresholds for smaller sub- while keeping the difference with the reference below a cho- sets to prevent errors from accumulating in bigger sets. sen threshold. The entire algorithm is described in the pseudo-code pre- sented in Appendix A, which also contains the instructions www.geosci-model-dev.net/12/3135/2019/ Geosci. Model Dev., 12, 3135–3148, 2019 3140 O. Tintó Prims et al.: How to use mixed precision in ocean models: NEMO 4.0 and ROMS 3.6 3.1.1 Emulator implementation compare and the kind of test. For this analysis, we are using as a reference a 64-bit real-variable version of the code. The We have developed a tool (see “Code availability” section) runs consisted of 10 d simulations that produced daily out- that not only modifies the source code to implement the em- puts of the 3-D temperature field, the 3-D salinity field, and ulator solving all the issues mentioned in the Methods sec- the column-integrated heat and salt content. Temperature and tion, but also creates a database with information about the salinity were selected because these are the two active trac- sources, including its modules, routines, functions, variables ers that appear in the model equations. Reference simulations and their relations. This database will have several uses after- were launched for each different initial condition that was wards, since it can be used to generate the list of precisions used later in the analysis. assigned to each variable and to process the results of the In order to define some meaningful thresholds, we per- analysis. formed a simulation using a halved time step to compute the The code largely relies on MPI and NetCDF. Since there is accuracy score against the reference simulation; this value no special interest in analyzing these routines, a simple solu- was used as a first guess to define our thresholds (see Ta- tion is to keep them unmodified. The selection of the source ble 1). To show how defining different thresholds leads to files that should not be modified requires user expertise. After different results and to emphasize how with this method it is that the tool handles all the necessary workarounds to ensure possible to keep arbitrarily small output differences, two dif- that the proper variable type is passed as routine arguments. ferent thresholds were defined (see Table 1). The first of the two thresholds is defined to be 1000 times smaller than the 3.1.2 Designing the accuracy tests differences obtained using a halved time step. The second is defined to be only 10 times smaller. We will refer to these Our approach was to define a metric to evaluate how similar cases as the tight case and the loose case, respectively. two simulations are and define a threshold above which we will consider simulations inaccurate. 3.1.3 Executing the tests The metric used to evaluate the similarity between two simulation outputs is the root mean square deviation divided To execute the tests, we implemented the rules described in by the interquartile range. It is computed for each time step, the Methods section by developing a Python workflow man- and the maximum value is retained. ager capable of handling the dynamic workflow by creating RMSD(t)) = the simulation scripts, launching them on a remote platform, v and then checking the status and adequacy of the results. u max max max refi,j,k (t) − xi,j,k (t)2 u Pi Pj Pk t i=1 j =1 k=1 To perform tests, the Marenostrum 4 supercomputer has , (1) been used (see Table 2), with each individual simulation tak- imax · jmax · kmax ing about 8 min 30 s on a single node. where i, j, k are the spatial axis indices, t is the time index, During the implementation of the emulator, the declara- refi,j,k (t) is the value of the reference simulation at a given tions of more than 3500 real variables were replaced with point i, j, k and a given time t, and xi,j,k (t) the value of the emulator variables. These variables can be scalars or arrays simulation that is being evaluated at the same point i, j, k and of up to five dimensions. They can be global variables or just a given time t. temporary variables used only in a single line of code inside a subroutine. In a large code like NEMO there are usually IQR(t) = Q3(t) − Q1(t), (2) parts of the code that are either used or not depending on the specific parameters employed. For this reason, even when where Q3 and Q1 are the values of the third quartile and first more than 3500 variables were identified, using our specific quartile, respectively, and t is the time index. configuration only 942 are used. The variables that are used So the final metric will be the accuracy score: are identified during the runtime. Given that considering the   RMSD(t) unused variables for the analysis would be useless and dan- Accuracy score = max . (3) gerous (as incorrect conclusions could be drawn about what ∀t IQR(t) precision is needed) those are simply omitted from the anal- The final accuracy test can thus be defined as ysis. The initial set used to start the analysis will therefore be the set containing the 942 variables that are actually used Accuracy test (accuracy score, threshold) = ( during our specific case. True if accuracy score < threshold (4) 3.1.4 Discussion False if accuracy score ≥ threshold. To be able to use the test defined above we need a reference Starting with the evaluation of the initial set containing all simulation and to establish the thresholds. To obtain the ref- the variables, the results show that a global reduction of the erence simulation, first we have to define what we want to precision changes the results beyond what is acceptable. Af- Geosci. Model Dev., 12, 3135–3148, 2019 www.geosci-model-dev.net/12/3135/2019/ O. Tintó Prims et al.: How to use mixed precision in ocean models: NEMO 4.0 and ROMS 3.6 3141 Table 1. Accuracy score of the simulation performed using a time step shorter than the reference and the different thresholds used for the different variables in the tight and loose cases. The thresholds for the tight case are defined by multiplying the accuracy score of the half- time-step simulations by 10−3 , and the thresholds for the loose case are defined by multiplying the accuracy score of the half-time step by 10−1 . Variable Accuracy score Thresholds Thresholds half-time step tight case loose case (AS) (AS ×10−3 ) (AS ×10−1 ) Temperature (3-D field) 1.49 × 10−3 1.49 × 10−6 1.49 × 10−4 Salinity (3-D field) 6.41 × 10−3 6.41 × 10−6 6.41 × 10−4 Heat content (2-D field) 2.74 × 10−3 2.74 × 10−6 2.74 × 10−4 Salt content (2-D field) 5.47 × 10−6 5.47 × 10−9 5.47 × 10−7 Table 2. Marenostrum 4 node specifications. Processor Two-socket Intel Xeon Platinum 8160 CPU with 24 cores each @ 2.10 GHz for a total of 48 cores per node Memory L1d 32K; L1i cache 32K; L2 cache 1024K; L3 cache 33792K 96 GB of main memory 1.880 GB per core Network 100 Gbit s−1 Intel Omni-Path HFI Silicon 100 Series PCI-E adapter 10 Gbit Ethernet Local storage 200 GB local solid-state drive available as temporary storage during jobs ter running the full analysis for both cases we can see that Table 3. The table presents the total number of variables that were the results differ. On the one hand, the tight case shows that included in the analysis split by the dimension of the arrays (0 for 652 variables (69.2 %) could use single precision and, on the scalars). Also shown are the number of variables that can safely use other hand, in the loose case the obtained solution contained single precision for the two cases discussed in Sect. 3.1 (tight and 902 variables (95.8 %), thereby keeping only 40 variables in loose) and the percentage of the total variables that this represents. Using tight constraints, 69.2 % of the variables can use single preci- higher precision. As expected, defining looser thresholds al- sion, although only 52.5 % of the 3-D variables are responsible for lows for the use of lower precision in a larger portion of the the most memory usage. Using loose conditions, 95.8 % of the vari- variables. Table 3 presents the results split by array dimen- ables can use single precision, including 95 % of the 3-D variables. sion. Also, the tighter case required more tests, with 1442 simulations required to arrive at the solution (consuming Tight case Loose case 204.3 node hours), while the loose case required only 321 (consuming 45.5 node hours). Dimension Total n % n % variables The method ensures that the differences between the reduced-precision simulations performed with the final vari- 0 374 304 81.3 366 97.9 able sets and the reference are below the determined thresh- 1 32 23 71.9 32 100.0 olds. Although the accuracy test used for this example was 2 322 216 67.1 307 95.3 not designed to ensure the conservation of global quantities, 3 200 105 52.5 190 95.0 Fig. 1 shows that the global heat and salt content of the sim- 4 13 3 23.1 6 46.2 ulations performed using the loose configuration resembles 5 1 1 100.0 1 100.0 those of the reference simulation performed fully in double Total 942 652 69.2 902 95.8 precision, which was not the case for a simulation fully per- formed using single precision. The Fig. 2 shows the expected impact on memory usage when using the set obtained with the loose accuracy test. 3.2 ROMS Since the results produced by both cases are so close to the reference simulation, Fig. 2 only demonstrates the loose case The Regional Ocean Modeling System (ROMS) is a free- because we anticipate that a mixed-precision implementation surface, hydrostatic, primitive-equation ocean model that that is accepted by the scientific community would be more uses stretched, terrain-following coordinates on the verti- similar to the loose case. cal and orthogonal curvilinear coordinates on the horizon- tal. It contains a variety of features, including high-order ad- vection schemes, accurate pressure gradient algorithms, sev- eral subgrid-scale parameterizations, atmospheric, oceanic, www.geosci-model-dev.net/12/3135/2019/ Geosci. Model Dev., 12, 3135–3148, 2019 3142 O. Tintó Prims et al.: How to use mixed precision in ocean models: NEMO 4.0 and ROMS 3.6 Figure 1. Daily averages of the global heat content (a) and global salt content (b) for 1-month simulations using NEMO 4. The plot of the global salt content (b) has an offset of 4.731 × 1022 g. The reference simulation was performed using double precision, the single-precision simulation was performed using single-precision arithmetic for all the variables and the mixed-precision simulation was performed keeping 40 of the 942 model real variables in double precision and using single precision for the remaining 902 variables. resolution and 30 levels on the vertical (http://www.myroms. org, last access: 18 October 2018). This configuration was selected because 4DVAR ROMS has a large community of users, there is an easy-to-follow tutorial to set up the config- uration and it involves linear models that make this an inter- esting choice to expand the results obtained in Sect. 3.1. To perform this analysis, we set up the tutorial available online that performs an I4DVAR data assimilation cycle that spans the period 3–6 January 2004. The observations assim- ilated into the model are satellite sea surface temperature, satellite sea surface height in the form of a gridded prod- uct from Aviso, and hydrographic observations of tempera- ture and salinity collected from Argo floats and during the GLOBEC/LTOP and CalCOFI cruises off the coast of Ore- gon and southern California, respectively. With the I4DVAR version of ROMS, most of the com- putational time is spent inside the data assimilation “inner loops” within the adjoint and tangent linear sub-models (AD and TL, respectively). To give some numbers, in a simulation with the WC13 configuration, 35 % of the time is spent in the Figure 2. Estimation of the memory usage impact in an ORCA2 TL model and 50 % in the AD model, while the time spent in configuration using the set obtained with the loose accuracy test. the nonlinear model is below 14 %. It was felt that since the The estimation has been performed using the dimensions of the TL and AD models are used to compute an approximation variables and the size of the domain. In red (darker) we have the memory that is occupied by variables in double precision (100 % in gradient of the I4DVAR cost function, further approxima- the original case), and in light green the proportion of the memory tions in the gradient resulting from lower precision will prob- occupied by variables in single precision; the difference between the ably be not so detrimental, leading to our starting hypothesis original case and the loose case implies a potential 44.9 % decrease that the AD and TL models are better targets than the non- in memory usage. linear model. However, for this hypothesis to be true the AD needs to keep being the transpose of the TL to within a cho- sen precision; otherwise, it might prevent the convergence of and benthic boundary layers, biological modules, radia- the algorithm. For this reason, we are trying to minimize the tion boundary conditions, and data assimilation (http://www. precision used in these regions of the code. In this case, the myroms.org, last access: 18 October 2018). target reduced precision is not limited to single precision (23- The experiments performed with ROMS were done by ap- bit significand), but the half-precision is also explored (10-bit plying the primal form of the incremental strong constraint significand). 4D-Var (4DVAR) (Moore et al., 2011). The configuration It is important to remark that in our experiments the num- used is the US west coast at 1/3◦ horizontal, referred to as ber of bits used to represent the exponent of the variables is WC13, a standard model test case. It has ∼ 30 km horizontal 11 (the same number as used in double precision), no matter Geosci. Model Dev., 12, 3135–3148, 2019 www.geosci-model-dev.net/12/3135/2019/ O. Tintó Prims et al.: How to use mixed precision in ocean models: NEMO 4.0 and ROMS 3.6 3143 what number of bits is used for the significand. In practice, it Table 4. The table presents the total number of variables that were means that to ensure that it is possible to reduce the precision included in the ROMS analysis split by the dimension of the arrays to half-precision we would need to ensure that the values do (0 for scalars). Also shown are the number of variables that can not exceed the range that can be represented with a smaller safely use reduced precision for the two cases discussed in Sect. 3.2 exponent. However, the results are still interesting because (single precision and half-precision) and the percentage of the to- tal variables that they represent. The single-precision case was triv- we can learn the potential use of a reduced number of bits ial since all the variables considered can use single precision while for the significand. keeping the accuracy. For the half-precision case, 80.7 % of the vari- ables could use half-precision, although the proportion of variables 3.2.1 Emulator implementation that can use reduced precision becomes lower as the dimension in- creases. One aspect that made this exercise different from that of NEMO in terms of implementing the emulator was that the 23 bits 10 bits interest was focused on specific regions of the code, i.e., the parts related to the tangent linear model and the adjoint Dimension Total n % n % variables model. 0 326 326 100.0 277 85.0 3.2.2 Designing the tests 1 82 82 100.0 71 86.6 2 450 450 100.0 385 85.6 When running a data assimilation experiment, the objective 3 173 173 100.0 135 78.0 is to obtain a coherent state of the ocean that minimizes the 4 78 78 100.0 44 56.4 difference between the model and the observations. Through 5 31 31 100.0 9 29.0 different forward–backward iterations using the TL and AD 6 4 4 100.0 2 50.0 models, the model should converge to a state in which the 40 2 2 100.0 2 100.0 cost function (function that describes the difference between Total 1146 1146 100.0 925 80.7 the model state and the observations) is minimum. We can consider a simulation to be accurate enough if the model con- verges to the same result as the double-precision reference. Through the different iterations, the solution should converge numerical precision. To test that, we proceeded with a new to a minimum. To set a threshold for the accuracy of the sim- analysis using a 10-bit significand. The first simulation of the ulations, we can look at the difference between the last two analysis crashed, not providing any output and showing that inner-loop iterations. In the reference simulation the value it is not possible to reduce the precision to 10 bits overall. of this difference computed as a cost function defined in the Using the analysis to find where this precision can be used, model is 1.77 × 10−1 . Defining a threshold 10 times smaller the results show that 80.7 % of the variables can use 10 bits (1.77 × 10−2 ) ensures that the impact of reducing the preci- instead of 52. Looking to the effects produced by the vari- sion will be smaller than changing the number of inner itera- ables that were not able to use 10-bit precision we could see tions performed. that only a single variable was preventing the model from completing the simulations, even though many others were 3.2.3 Executing the tests introducing too many numerical errors. In Table 4 the results are presented split by dimension. For ROMS, we use the same workflow manager developed for NEMO, simply using different templates for launching 3.3 Common discussion simulations and evaluating the results. 3.2.4 Discussion From the two exercises shown in Sect. 3 we can draw sev- eral conclusions. The most important is that the method can As in the case with NEMO, we first performed the reference provide a set of variables that can use reduced numerical pre- executions for all the different initial conditions using 64-bit cision while preserving the accuracy of the results, as was precision. observed in the four cases explored (two different accuracy The parts selected for the analysis contain 1556 real- tests for NEMO and two target precisions for ROMS). The variable declarations, 1144 of which are actually used in our results in the four cases showed that the potential to reduce case. Starting with this initial set and trying to use a 23-bit the numerical precision is considerable even in the most con- significand, the results showed that in fact all the variables strained cases (i.e., NEMO tight). In the NEMO case, the can use reduced precision at once while still obtaining accu- analysis covers the full model, while the ROMS case was fo- rate results. cused only on a part of the model, demonstrating the versa- This suggests that this model and in particular this specific tility of the method. It has also been shown that it is possible configuration are suitable for a more drastic reduction of the to achieve results that are arbitrarily close to the reference, www.geosci-model-dev.net/12/3135/2019/ Geosci. Model Dev., 12, 3135–3148, 2019 3144 O. Tintó Prims et al.: How to use mixed precision in ocean models: NEMO 4.0 and ROMS 3.6 while requiring more exact results will have the cost of leav- Code and data availability. Releases of the Reduced Precision ing more parts of the code in higher precision. Emulator library are available at https://github.com/aopp-pred/rpe (Dawson and Dueben, 2016, last access: January 2019). The tool developed to automatically implement RPE can be found at https:// 4 Conclusions earth.bsc.es/gitlab/otinto/AutoRPE (Tintó Prims, 2019). The work- flow manager used to manage the analysis simulations can be found Previous works suggested that the generalized use of double at https://earth.bsc.es/gitlab/otinto/AutoRPE (Tintó Prims, 2019). precision for Earth science modeling is a case of overengi- The source code of the NEMO model and the input data used for neering and that for this reason adapting the computational the experiments referred to in this work can be found at https: models to use the smallest numerical precision that is essen- //www.nemo-ocean.eu/ (last access: January 2019). The source tial would compensate in terms of improvement in perfor- code of the ROMS model can be found at http://www.myroms.org/ (The ROMS/TOMS Group, 2019) and can be freely accessed un- mance. However, an improper reduction can lead to accuracy der registration. The details and input data used for the experiments losses that may make the results unreliable. In this paper we involving ROMS can be found at https://www.myroms.org/wiki/ presented a method that was designed to solve this problem I4DVAR_Tutorial (last access: January 2019). by finding which variables can use a lower level of precision without degrading the accuracy of the results. The method was designed to be applied to computational ocean models coded in Fortran. It relies on the Reduced Pre- cision Emulator tool, which was created to help us under- stand the effect of reduced precision within numerical sim- ulations. It allows us to simulate the results that would be produced in the case of performing arithmetic operations in reduced precision using a custom number of significant bits for each floating-point variable. The proposed analysis algo- rithm finds which variables need to be kept in double preci- sion and which ones can use reduced precision without im- pacting the accuracy of the results. The method has been tested with two widely used state- of-the-art ocean models, NEMO and ROMS. The experi- ences with the two models pursued different objectives. With NEMO, the analysis covered all the routines and variables used within the ocean-only simulation; the target precision was 32 bits, and we explored how the selection of specific accuracy tests changes which variables can safely use re- duced precision. However, with ROMS, the analysis covered only the variables belonging to the adjoint and tangent linear models using a single accuracy test and examining how the method can be used to discover the viability of using numer- ical precisions below 32 bits. The results presented in this work allow us to draw some conclusions. It is shown that both models can use reduced precision for large portions of the code, proving the feasibil- ity of mixed-precision approaches, and the method described is able to find these. The method can also provide a configu- ration that can be arbitrarily close to the double-precision ref- erence, whereby the amount of variables that can use reduced precision will depend on how strict the imposed conditions are. Users might want to follow the method presented here to build a version of the model that can benefit from a reduc- tion of precision and can be used by the whole community without concerns. Geosci. Model Dev., 12, 3135–3148, 2019 www.geosci-model-dev.net/12/3135/2019/ O. Tintó Prims et al.: How to use mixed precision in ocean models: NEMO 4.0 and ROMS 3.6 3145 Appendix A www.geosci-model-dev.net/12/3135/2019/ Geosci. Model Dev., 12, 3135–3148, 2019 3146 O. Tintó Prims et al.: How to use mixed precision in ocean models: NEMO 4.0 and ROMS 3.6 Appendix B Figure B1. Difference in sea surface temperature monthly mean (◦ C) for the first month of simulation between a NEMO 4.0 simulation performed emulating single-precision arithmetic and a NEMO 4.0 simulation performed using double-precision arithmetic. The monthly mean of the sea surface temperature for the first month of simulation shows localized biases exceeding 1.5 ◦ C. Geosci. Model Dev., 12, 3135–3148, 2019 www.geosci-model-dev.net/12/3135/2019/ O. Tintó Prims et al.: How to use mixed precision in ocean models: NEMO 4.0 and ROMS 3.6 3147 Author contributions. OTP led the design of the method, the exper- Düben, P. D., Subramanian, A., Dawson, A., and Palmer, T. N.: iments and the writing of the paper. MCA supervised all the pro- A study of reduced numerical precision to make superparam- cesses and had a main role in the writing of the paper. AMM par- eterization more competitive using a hardware emulator in ticipated in the design of the experiments involving ROMS and the the OpenIFS model, J. Adv. Model. Earth Sy., 9, 566–584, writing of the paper. MC was involved in the discussion and design https://doi.org/10.1002/2016MS000862, 2017. of the method and the experiments and participated in the writing of Düben, P. D., Subramanian, A., Dawson, A., and Palmer, T. N.: the paper. KS, AC and FJD participated in the discussions and the A study of reduced numerical precision to make superparam- revision of the paper. eterization more competitive using a hardware emulator in the OpenIFS model, J. Adv. Model. Earth Sy., 9, 566–584, https://doi.org/10.1002/2016MS000862, 2017. Competing interests. The authors declare that they have no conflict Graillat, S., Fabienne, J., Picot, R., and Lathuili, B.: PROMISE: of interest. floating-point precision tuning with stochastic arithmetic, 17th international symposium on Scientific Computing, Computer Arithmetic and Verified Numerics (SCAN 2016), UPPSALA, Financial support. The research leading to these results has re- Sweden, 98–99, September 2016. ceived funding from the EU ESiWACE H2020 Framework Pro- Haidar, A., Tomov, S., Dongarra, J., and Higham, N. J.: Harness- gramme under grant agreement no. 823988, from the Severo ing GPU Tensor Cores for Fast FP16 Arithmetic to Speed up Ochoa (SEV-2011-00067) program of the Spanish Government and Mixed-Precision Iterative Refinement Solvers, in: Proceedings from the Ministerio de Economia y Competitividad under contract of the International Conference for High Performance Comput- TIN2017-84553-C2-1-R. ing, Networking, Storage, and Analysis, SC’18 (Dallas, TX), IEEE Press, Piscataway, NJ, USA, 47:1–47:11, available at: http://dl.acm.org/citation.cfm?id=3291656.3291719 (last access: July 2019), 2018. Review statement. This paper was edited by Steven Phipps and re- Hegerl, G. and Zwiers, F.: Use of models in detection and at- viewed by Peter Düben and Matthew Chantry. tribution of climate change, Wires Clim. Change, 2, 570–591, https://doi.org/10.1002/wcc.121, 2011. Kniveton, D. R., Smith, C. D., and Black, R.: Emerging migration flows in a changing climate in dryland Africa, Nat. Clim. Change, References 2, 444–447, https://doi.org/10.1038/nclimate1447, 2012. Lam, M. O., De Supinksi, B. R., Legendre, M. P., and Aumont, O., Ethé, C., Tagliabue, A., Bopp, L., and Gehlen, Hollingsworth, J. K.: Automatically adapting programs for M.: PISCES-v2: an ocean biogeochemical model for carbon mixed-precision floating-point computation, Proceedings – and ecosystem studies, Geosci. Model Dev., 8, 2465–2513, 2012 SC Companion: High Performance Computing, Net- https://doi.org/10.5194/gmd-8-2465-2015, 2015. working Storage and Analysis, SCC 2012, 1423–1424, Baboulin, M., Buttari, A., Dongarra, J., Kurzak, J., Lan- https://doi.org/10.1109/SC.Companion.2012.231, 2012. gou, J., Langou, J., Luszczek, P., and Tomov, S.: Ac- Madec, G.: NEMO ocean engine, Note du Pole de mod- celerating scientific computations with mixed precision elisatio 1288–1619, https://www.nemo-ocean.eu/wp-content/ algorithms, Comput. Phys. Commun., 180, 2526–2533, uploads/Doc_OPA8.1.pdf (last access: July 2019), 2008. https://doi.org/10.1016/j.cpc.2008.11.005, 2009. Moore, A., Arango, H., Broquet, G., S. Powell, B., Weaver, A., Bauer, P., Thorpe, A., and Brunet, G.: The quiet revolu- and Zavala-Garay, J.: The Regional Ocean Modeling System tion of numerical weather prediction, Nature, 525, 47–55, (ROMS) 4-dimensional variational data assimilation systems https://doi.org/10.1038/nature14956, 2015. Part I – System overview and formulation, Prog. Oceanogr., 91, Bellard, C., Bertelsmeier, C., Leadley, P., Thuiller, W., and Cour- 34–49, https://doi.org/10.1016/j.pocean.2011.05.004, 2011. champ, F.: Impacts of climate change on the future of biodi- Oreskes, N., Stainforth, D. A., and Smith, L. A.: Adap- versity, Ecol. Lett., 15, 365–377, https://doi.org/10.1111/j.1461- tation to Global Warming: Do Climate Models Tell Us 0248.2011.01736.x, 2012. What We Need to Know?, Philos. Sci., 77, 1012–1028, Dawson, A. and Dueben, P.: aopp-pred/rpe: v5.0.0 (Version v5.0.0), https://doi.org/10.1086/657428, 2010. Zenodo, https://doi.org/10.5281/zenodo.154483, 2016. Palmer, T., Shutts, G., Hagedorn, R., Doblas-Reyes, F., Jung, T., Dawson, A. and Düben, P. D.: rpe v5: an emulator for reduced and Leutbecher, M.: Representing model uncertainty in wheather floating-point precision in large numerical simulations, Geosci. and climate prediction, Annu. Rev. Earth Pl. Sc., 33, 163–193, Model Dev., 10, 2221–2230, https://doi.org/10.5194/gmd-10- https://doi.org/10.1146/annurev.earth.33.092203.122552, 2005. 2221-2017, 2017. Randall, D., Wood, R., Bony, S., Colman, R., Fichefet, T., Fyfe, J., Dennis, J. M. and Loft, R. D.: Refactoring Scientific Applications Kattsov, V., Pitman, A., Shukla, J., Srinivasan, J., Stouffer, R., for Massive Parallelism, in: Numerical Techniques for Global At- Sumi, A., and Taylor, K.: Climate Models and Their Evaluation, mospheric Models, 80, 539–556, https://doi.org/10.1007/978-3- in: IPCC, 2007: Climate Change 2007: the physical science ba- 642-11640-7_16, 2011. sis. contribution of Working Group I to the Fourth Assessment Düben, P. D., McNamara, H., and Palmer, T. N.: The use Report of the Intergovernmental Panel on Climate Change, 323, of imprecise processing to improve accuracy in weather 589–662, https://doi.org/10.1016/j.cub.2007.06.045, 2007. & climate prediction, J. Comput. Phys., 271, 2–18, https://doi.org/10.1016/j.jcp.2013.10.042, 2014. www.geosci-model-dev.net/12/3135/2019/ Geosci. Model Dev., 12, 3135–3148, 2019 3148 O. Tintó Prims et al.: How to use mixed precision in ocean models: NEMO 4.0 and ROMS 3.6 Rousset, C., Vancoppenolle, M., Madec, G., Fichefet, T., Flavoni, Tintó Prims, O., Castrillo, M., Acosta, M., Mula-Valls, O., S., Barthélemy, A., Benshila, R., Chanut, J., Levy, C., Masson, Sanchez Lorente, A., Serradell, K., Cortés, A., and Doblas- S., and Vivier, F.: The Louvain-La-Neuve sea ice model LIM3.6: Reyes, F.: Finding, analysing and solving MPI communica- global and regional capabilities, Geosci. Model Dev., 8, 2991– tion bottlenecks in Earth System models, J. Comput. Sci., 3005, https://doi.org/10.5194/gmd-8-2991-2015, 2015. https://doi.org/10.1016/j.jocs.2018.04.015, in press, 2018. The ROMS/TOMS Group: ROMS, available at: http://www. Váňa, F., Düben, P., Lang, S., Palmer, T., Leutbecher, M., Salmond, myroms.org/, last access: July 2019. D., and Carver, G.: Single Precision in Weather Forecasting Thornes, T.: Can reducing precision improve accuracy in Models: An Evaluation with the IFS, Mon. Weather Rev., 145, weather and climate models?, Weather, 71, 147–150, 495–502, https://doi.org/10.1175/MWR-D-16-0228.1, 2017. https://doi.org/10.1002/wea.2732, 2016. Whiteman, G. and Hope, C. W. P.: Vast costs of Arctic change, Na- Tintó Prims, O.: AutoRPE, available at: https://earth.bsc.es/gitlab/ ture, 499, 401–403, https://doi.org/10.1038/499401a, 2013. otinto/AutoRPE, last access: July 2019. Geosci. Model Dev., 12, 3135–3148, 2019 www.geosci-model-dev.net/12/3135/2019/