Performance improvements to modern hydrological models via lookup table optimizations

Highlights

• The Canadian Hydrological Model (CHM) is profiled and expensive mathematical functions identified.

• FunC was used to replace the expensive mathematical functions in CHM with lookup tables.

• The run-time performance of CHM was improved by approximately 20% on two realistic simulations.

• A general methodology for using FunC to replace expensive mathematical functions with lookup tables is given.

Abstract

Distributed hydrological models predict the spatial variability in processes that govern observed mass and energy fluxes. A challenge associated with the use of these models is the computational burden associated with representing the Earth's (sub)surface via millions of computational elements. This burden is exacerbated as more complex process representations are included because their parameterizations involve computationally intensive mathematical functions. Lookup tables (LUTs) approximate a mathematical function by interpolating precomputed values of the function. Highly accurate approximations are possible for substantially reduced computational costs. In this work, a general methodology using the C++ LUT library FunC is applied to identify and replace computationally intensive mathematical function evaluations in the Canadian Hydrological Model (CHM). The use of LUTs introduces a pointwise relative error below ${10}^{-8}$ and provides a reduction in run time of almost 20%. This work shows how LUTs can be implemented with relatively little pain and yield significant computational savings for distributed hydrological models.

Introduction

Computer simulation is an increasingly important tool to address both research and operational water resource issues related to the hydrological cycle. These issues include improving our understanding of complex interactions affecting streamflow, helping to inform water management decisions, quantifying the effects of changes in climate and landcover, determining the changing frequency of extreme events, predicting flood impacts, and designing future infrastructure (Freeze and Harlan, 1969; Mote et al., 2005; Milly et al., 2008; Nazemi et al., 2013; Debeer et al., 2015; Wheater, 2015; Musselman et al., 2017). Distributed hydrological models represent the surface of the Earth using computational elements defined on regular or irregular grids (e.g., rasters or triangles) (Marsh et al., 2020b). The use of distributed models can provide key insights into the spatial variability of hydrological processes that govern the spatial heterogeneity in mass and energy fluxes (Dornes et al., 2008; Fatichi et al., 2016; Kumar et al., 2013). However, as the spatial extents investigated increase or the spatial scales of the computational elements decrease, there are limitations to the practicality of distributed models, partly due to computational costs associated with the simulations. Certainly, there are open questions regarding equifinality (Beven, 1993), parameter estimation (Beven, 1993; Savenije, 2009), initial and boundary conditions, and uncertainty in distributed hydrological models; however, model complexity and structure decisions are often influenced by computational cost (Fatichi et al., 2016).

The use of distributed hydrological models in practice is particularly challenging in the complex topography of mountainous regions. In these ‘‘water towers of the world’’ (Viviroli et al., 2007), seasonal snowpacks store a significant portion of the annual precipitation. The melting of these spatially variable snowpacks generates the spring freshet that can result in the largest discharge of the year (Gray and Male, 1981; Davies et al., 1987). This discharge provides a critically important source of fresh water to many downstream regions (Groisman and Davies, 2001; Mote et al., 2005; Stewart et al., 2008; Woo et al., 2008), including ecosystems (Davies et al., 1987; Groisman and Davies, 2001; Stewart et al., 2008; Woo et al., 2008; Boyd and Banzhaf, 2007) and agricultural, industrial, and municipal users (Harder et al., 2019; Nazemi et al., 2013). In order to accurately predict spring snowcovers, there is substantial motivation to simulate the mid-winter transport, both horizontal (blowing snow) (Pomeroy and Bernhardt, 2017; Marsh et al., 2020a; Mott et al., 2018) and vertical (avalanche) (Bernhardt et al., 2012) that, in conjunction with substantial horizontal variability in energy fluxes, leads to variability in the spring snowcover ablation. However, sub-hyper-resolution (sub-1 km) (Wood et al., 2011) and snow-drift-resolving scales (1 m–100 m) (Pomeroy and Bernhardt, 2017) are required to resolve this process heterogeneity, and these requirements impose a computational burden that limits the practical use of such models. Furthermore, ‘‘one-and-done’’ model runs are no longer generally acceptable, and uncertainty estimations require running a model tens, hundreds, or even thousands of times (Razavi et al., 2018). Thus, there is a critical need to adopt code optimization techniques to reduce the computational burden of distributed hydrological models; many such techniques are common in other computational sciences and can be directly applied to hydrological modelling.

Lookup Tables (LUTs) are an optimization approach that replaces the direct evaluation of mathematical functions with approximations based on sampling the input space of the function, precomputing and storing function values at the sample points, and then using interpolation to provide function evaluations. LUTs need not generally be applied to the entire model. Rather, they can be strategically applied to computationally costly mathematical functions within the model, e.g., an empirical parametrization that utilizes many mathematical functions. The benefit of LUTs is that fewer floating-point operations (FLOPs) are used, generally leading to decreased model execution times, albeit at the cost of decreased accuracy and increased data transfer and storage.

There are generally two categories of LUTs: software-based LUTs and hardware-based LUTs. Software-based LUTs make use of the general architecture of a CPU for their implementation. They find use in areas of numerical computing such as heart simulation (Cooper et al., 2006; Green et al., 2019) and atmospheric radiative transfer (Buehler et al., 2005). In contrast, hardware-based LUTs make use of specific hardware for storage and evaluation and find use in image processing (Battiato and Lukac, 2008; Gonzales and Woods, 2002; Pharr and Fernando, 2005), field programmable gate arrays (FPGAs) (Kuon et al., 2008), and hardware neural networks (Kumar Meher, 2010; Reis et al., 2014; Dias et al., 2014). Although it is true that hardware-based LUTs typically outperform their software-based counterparts, in practice, the use of software-based LUTs is an efficient and effective optimization for computationally intensive and repeatedly called mathematical functions (Green et al., 2019) and can lead to significant decreases in model execution time when run on a standard CPU. In this study, the focus is exclusively on software-based LUTs.

Historically, the implementation of LUTs required significant manual development work and resulted in software that was difficult to maintain (Loh et al., 2005). Wilcox et al. (2011) developed the MESA software package for the autonomous generation of LUTs. The MESA package is capable of generating LUTs according to an evaluation error tolerance, but it does not allow users to change the type of interpolation that is used — it uses only a basic piecewise constant interpolation over all subintervals. It can be shown that when using piecewise constant interpolation methods, the table sizes required to represent continuous functions to high precision can easily be on the order of terabytes or even greater (Wilcox et al., 2011; Green et al., 2019). Another downside of MESA is that it must be run external to the application where the LUT is to be used. A more modern approach to the generation of LUTs, such as the Function Comparator library (FunC) (Green et al., 2019) for C++, addresses these limitations by providing optimized implementations of higher-order interpolation methods as well as dynamic table generation facilities that can be incorporated directly into the desired C++ application.

FunC is a C++ library that facilitates the generation and assessment of LUT implementations for univariate function evaluation over bounded input spaces. FunC provides i) LUT implementations that use uniformly spaced sample points with interpolation schemes up to third order, ii) a LUT generation tool that can determine the required table parameters for each implementation for a given function, domain, and evaluation tolerance, and iii) benchmarking tools for evaluating and comparing the performance of the various LUT implementations against each other and against the direct evaluation of the original function.

Investigation and implementation of LUTs is thus no longer a significant obstacle for univariate functions. With guarantees over total error introduced as well as automatic generation and estimates for performance gains, these modern approaches allow LUTs to be readily included in existing software with unprecedented convenience.

In this work, the Canadian Hydrological Model (CHM) (Marsh et al., 2020b) is used in conjunction with FunC to determine whether CHM would benefit from LUT optimizations. Specifically: a) are there univariate mathematical functions in CHM that are computationally expensive and would benefit from LUT optimization?, b) do these benefits offer sufficient performance gains so as to warrant inclusion?, and c) are there general criteria for LUT use in hydrological models?

The remainder of the paper proceeds as follows. Section 2 describes CHM and its configuration for this study as well as FunC and its LUT capabilities. A general methodology is also given for the incorporation of LUTs in application software through the use of code profiling and FunC. Candidate functions are identified for replacement by LUTs, and details on LUT construction and use are given. Section 3 gives the results of applying this methodology within CHM on two realistic simulations. Section 4 discusses these results and offers general takeaways from this study to other applications. Finally, section 5 summarizes the conclusions from the study.