GPU-accelerated numerical solution to the Richards Equation:performance and prospects

Catchment-scale hydrological simulations typically require numerical solutions to the Richards equation, which describes variably-saturated water flow in porous media. In recent years, with the widespread use of high-performance computing (HPC), the simulation speed of hydrological models have been significantly enhanced. However, existing numerical schemes for the Richards equation show different performance under different HPC configurations. It remains unclear if the serial Richards solvers scale well in parallel. 

In this work, four popular numerical schemes (the fully explicit scheme, the predictor-corrector scheme, the iterative Picard scheme, and the fully implicit Newton's scheme) are implemented to solve the three-dimensional Richards equation, aiming at investigating the performance of these schemes on both multi-core CPUs and GPUs. The codes are built under the Kokkos framework to achieve performance portability between CPU and GPU. Two infiltration problems with analytical solutions available are chosen to evaluate the model performance in terms of accuracy and efficiency. 

As expected, the simulation results indicate that the optimal solution schemes on CPU and GPU could be different. Moreover, the numerical scheme, the linear system solver, and the soil properties all have influence on scaling. A hybrid scheme is promising for minimizing the computational cost under various simulation conditions. The findings of this work will guide the development of the subsurface flow module of the performance-portable, multi-physics model, SERGHEI.