基于PETSc的海岸带地下水数值模型MARUN的并行优化

海岸带地区地下水系统复杂，通常利用数值模拟仿真技术对地下水的流动和盐分运移过程进行定量刻画与研究。然而，随着所模拟的海岸带含水层尺度和维数的增加，传统的数值模拟方法面临着计算效率低、内存占用大等问题。开发基于并行技术的程序被认为是解决这类问题的有效技术手段。目前，已有大量的研究工作将并行技术在地下水模拟领域进行应用，但却极少探讨不同数值算法对求解海岸带地下水流-溶质运移耦合模型并行性能的影响。

本研究通过分析海岸带地下水数值模拟程序MARUN的结构框架以及并行框架PETSc的优越性，确定了以PETSc为基础，从矩阵、向量构造组装和线性系统求解两个方面对MARUN程序实现并行优化。通过采用两种Krylov子空间方法和五种预处理方法，分别对三个不同的算例进行强可扩展性测试或弱可扩展性测试，并分析对比了不同算法的并行性能。

研究结果表明，在使用较多处理器求解固定规模问题时，不同算法导致的计算性能差异较小。然而在处理核心数较少时，不同迭代算法和预处理方法对计算时间的影响显著，需要谨慎选取求解算法。此外，使用HYPRE库默认的BoomerAMG预处理器或BiCGSTAB子空间方法均能在用时最多的水头求解部分表现出较高的计算性能和可扩展性。对于流速和盐度的求解，SOR预处理器或Block-Jacobi和ILU结合的PETSc默认预处理器则是较为稳妥的选择。

Owing to the complexity of coastal groundwater systems, numerical simulation techniques are usually utilized to quantitatively study the flow system and salt transport processes. However, as the scale and dimension of the simulated coastal aquifer increase, traditional numerical simulation methods face challenges (e.g., low computational efficiency and large memory consumption). Optimizing programs based on parallel techniques is considered an effective way to solve such problems. At present, parallel techniques have been widely applied in the field of groundwater simulation, but there is little exploration on how different numerical algorithms impact the parallel performance of solving the coastal groundwater flow and solute transport coupling models.

This study elaborates on the structure of the coastal groundwater numerical program MARUN, and the superiority of the parallel framework PETSc. It is determined that the parallel optimization of MARUN program is based on PETSc from two aspects (i.e., matrix-vector construction and linear system solving). By using two Krylov subspace methods and five preconditioning techniques respectively, strong or weak scaling tests are conducted on three different examples. The parallel performance of different algorithms is therefore compared.

The results indicate that when using multiple processors to solve fixed-scale problems, the difference in computational performance caused by different algorithms is relatively small. However, when dealing with fewer processors, different Krylov subspace methods and preconditioning techniques have a significant impact on the runtime. The solving algorithms should be selected with caution. In addition, both the BiCGSTAB subspace method and the default BoomerAMG preconditioner from HYPRE can achieve high computational performance and scalability in the most time-consuming head-solving part. As for the solving process of velocity and salinity, SOR preconditioner or the default preconditioner of PETSc which combines Block-Jacobi with ILU is a more reliable choice.

[1] 周树荃. 有限元结构分析并行计算[M]. 科学出版社, 1994.
[2] Essink G H P. Mathematical models and their application to salt water intrusion problems[C]//TIAC’03: Coastal aquifers intrusion technology: Mediterranean countries. 2003: 57-77.
[3] Neumann B, Vafeidis A T, Zimmermann J, et al. Future coastal population growth and exposure to sea-level rise and coastal flooding-a global assessment[J]. PloS one, 2015, 10(3): e0118571.
[4] Verkaik J, Van Engelen J, Huizer S, et al. Distributed memory parallel computing of three-dimensional variable-density groundwater flow and salt transport[J]. Advances in Water Resources, 2021, 154: 103976.
[5] 苗新强, 金先龙, 丁峻宏. 结构静力有限元分层并行计算方法[J]. 力学学报, 2014, 46(4): 611-618.
[6] 余德浩. 大规模科学计算研究[J]. 中国基础科学, 2001年1期: 19-25.
[7] Eijkhout V. Introduction to high performance scientific computing[M]. North Carolina: Lulu Press, 2010.
[8] 张林波, 迟学斌, 莫则尧, 等. 并行计算导论[M]. 北京: 清华大学出版社, 2006.
[9] Delsman J R, Mulder T, Verastegui B R, et al. Reproducible construction of a high-resolution national variable-density groundwater salinity model for the Netherlands[J]. Environmental Modelling & Software, 2023, 164: 105683.
[10] Hwang H T, Park Y J, Sudicky E A, et al. A parallel computational framework to solve flow and transport in integrated surface–subsurface hydrologic systems[J]. Environmental Modelling & Software, 2014, 61: 39-58.
[11] Abdelaziz R, Le H H. MT3DMSP–A parallelized version of the MT3DMS code[J]. Journal of African Earth Sciences, 2014, 100: 1-6.
[12] Huang J, Christ J A, Goltz M N. An assembly model for simulation of large-scale ground water flow and transport[J]. Groundwater, 2008, 46(6): 882-892.
[13] Quarteroni A, Sacco R, Saleri F. Numerical mathematics[M]. Berlin: Springer Science & Business Media, 2006.
[14] Barlow R H, Evans D J. Parallel algorithms for the iterative solution to linear systems[J]. The Computer Journal, 1982, 25(1): 56-60.
[15] Gander W, Gander M J, Kwok F. Scientific computing-An introduction using Maple and MATLAB[M]. Berlin: Springer Science & Business Media, 2014.
[16] Saad Y, Van Der Vorst H A. Iterative solution of linear systems in the 20th century[J]. Journal of Computational and Applied Mathematics, 2000, 123(1-2): 1-33.
[17] Alyahya H, Mehmood R, Katib I. Parallel iterative solution of large sparse linear equation systems on the intel MIC architecture[J]. Smart Infrastructure and Applications: Foundations for Smarter Cities and Societies, 2020: 377-407.
[18] Xiao G, Yin C, Zhou T, et al. A survey of accelerating parallel sparse linear algebra[J]. ACM Computing Surveys, 2023, 56(1): 1-38.
[19] Sun Q, Zhang C, Wu C, et al. Bandwidth reduced parallel SpMV on the SW26010 many-core platform[C]//Proceedings of the 47th International Conference on Parallel Processing. 2018: 1-10.
[20] Fei X, Zhang Y. Regu2D: Accelerating vectorization of SpMV on intel processors through 2D-partitioning and regular arrangement[C]//Proceedings of the 50th International Conference on Parallel Processing. 2021: 1-11.
[21] Sadi F, Sweeney J, Low T M, et al. Efficient spmv operation for large and highly sparse matrices using scalable multi-way merge parallelization[C]//Proceedings of the 52nd Annual IEEE/ACM International Symposium on Microarchitecture. 2019: 347-358.
[22] Parravicini A, Sgherzi F, Santambrogio M D. A reduced-precision streaming SpMV architecture for Personalized PageRank on FPGA[C]//Proceedings of the 26th Asia and South Pacific Design Automation Conference. 2021: 378-383.
[23] Chen Y, Xiao G, Wu F, et al. tpSpMV: A two-phase large-scale sparse matrix-vector multiplication kernel for manycore architectures[J]. Information Sciences, 2020, 523: 279-295.
[24] Xiao G, Li K, Chen Y, et al. Caspmv: A customized and accelerative spmv framework for the sunway taihulight[J]. IEEE Transactions on Parallel and Distributed Systems, 2019, 32(1): 131-146.
[25] Bienz A, Gropp W D, Olson L N. Node aware sparse matrix–vector multiplication[J]. Journal of Parallel and Distributed Computing, 2019, 130: 166-178.
[26] Niu Y, Lu Z, Dong M, et al. TileSpMV: A tiled algorithm for sparse matrix-vector multiplication on GPUs[C]//2021 IEEE International Parallel and Distributed Processing Symposium (IPDPS). IEEE, 2021: 68-78.
[27] Li K, Yang W, Li K. Performance analysis and optimization for SpMV on GPU using probabilistic modeling[J]. IEEE Transactions on Parallel and Distributed Systems, 2014, 26(1): 196-205.
[28] Miller C T, Dawson C N, Farthing M W, et al. Numerical simulation of water resources problems: Models, methods, and trends[J]. Advances in Water Resources, 2013, 51: 405-437.
[29] Mo Z, Zhang A, Cao X, et al. JASMIN: a parallel software infrastructure for scientific computing[J]. Frontiers of Computer Science in China, 2010, 4: 480-488.
[30] Heroux M A, Bartlett R A, Howle V E, et al. An overview of the Trilinos project[J]. ACM Transactions on Mathematical Software (TOMS), 2005, 31(3): 397-423.
[31] Balay S, Brown J, Buschelman K, et al. PETSc users manual[R]. Tech. Rep. ANL-95/11—Revision 3.3, Argonne Natl. Lab., Argonne, Ill., 2012.
[32] Zhang K, Wu Y S, Bodvarsson G S. Parallel computing simulation of fluid flow in the unsaturated zone of Yucca Mountain, Nevada[J]. Journal of Contaminant Hydrology, 2003, 62-63: 381-399.
[33] White S P, Creighton A L, Bixley P E, et al. Modeling the dewatering and depressurization of the Lihir open-pit gold mine, Papua New Guinea[J]. Geothermics, 2004, 33(4): 443-456.
[34] Katz R F, Knepley M G, Smith B, et al. Numerical simulation of geodynamic processes with the Portable Extensible Toolkit for Scientific Computation[J]. Physics of the Earth and Planetary Interiors, 2007, 163: 52-68.
[35] White M D, Oostrom M. STOMP Subsurface Transport Over Multiple Phases, version 4.0, user's guide[R]. PNNL-15782, Pacific Northwest National Laboratory, Richland, 2006.
[36] White M D, Oostrom M, Rockhold M L, et al. Scalable modeling of carbon tetrachloride migration at the Hanford site using the STOMP simulator[J]. Vadose Zone Journal, 2008, 7(2): 654-666.
[37] Lu B, Wheeler M F. Iterative coupling reservoir simulation on high performance computers[J]. Petroleum Science, 2009, 6(1):43-50.
[38] Yamamoto H, Zhang K, Karasaki K, et al. Numerical investigation concerning the impact of CO2 geologic storage on regional groundwater flow[J]. International Journal of Greenhouse Gas Control, 2009, 3(5): 586-599.
[39] Kollet S J, Maxwell R M, Woodward C S, et al. Proof of concept of regional scale hydrologic simulations at hydrologic resolution utilizing massively parallel computer resources[J]. Water Resources Research, 2010, 46(4): W04201.
[40] Le P, Kumar P, Valocchi A J, et al. GPU-based high-performance computing for integrated surface–sub-surface flow modeling[J]. Environmental Modelling & Software, 2015, 73: 1-13.
[41] Sarris T S, Scott D M, Close M E. Parallel simulation in subsurface hydrology: evaluating the performance of modelling computers[J]. Groundwater, 2021, 59:109-116.
[42] Verkaik J, Hughes J D, van Walsum P E V, et al. Distributed memory parallel groundwater modeling for the Netherlands Hydrological Instrument[J]. Environmental Modelling & Software, 2021, 143: 105092.
[43] Horgue P, Renard F, Gerlero G S, et al. porousMultiphaseFoam v2107: An open-source tool for modeling saturated/unsaturated water flows and solute transfers at watershed scale[J]. Computer Physics Communications, 2022, 273:108278.
[44] Lu C, Lichtner P C. PFLOTRAN: massively parallel 3-D simulator for CO2 sequestration in geologic media[C]//Fourth Annual Conference on Carbon Capture and Sequestration DOE/NETL. 2005.
[45] Hammond G E, Lichtner P C, Mills R T. Evaluating the performance of parallel subsurface simulators: An illustrative example with PFLOTRAN[J]. Water Resources Research, 2014, 50(1):208-228.
[46] Wang W, Fischer T, Zehner B, et al. A parallel finite element method for two-phase flow processes in porous media: OpenGeoSys with PETSc[J]. Environmental Earth Sciences, 2015, 73: 2269-2285.
[47] Wang W, Kolditz O, Nagel T. Parallel finite element modelling of multi-physical processes in thermochemical energy storage devices[J]. Applied Energy, 2016, 185(pt.2):1954-1964.
[48] Su D, Mayer K U, Macquarrie K T B. Parallelization of MIN3P-THCm: A high performance computational framework for subsurface flow and reactive transport simulation[J]. Environmental Modelling and Software, 2017, 95: 271-289.
[49] Liu X, Zhang Q, Cheng T. Accelerating contaminant transport simulation in MT3DMS using JASMIN-based parallel computing[J]. Water, 2020, 12(5):1480.
[50] Boufadel M C, Suidan M T, Venosa A D. A numerical model for density-and-viscosity-dependent flows in two-dimensional variably-saturated porous media[J]. Journal of Contaminant Hydrology, 1999, 37(1-2): 1-20.
[51] Li H, Boufadel M. Long-term persistence of oil from the Exxon Valdez spill in two-layer beaches[J]. Nature Geoscience, 2010, 3: 96–99.
[52] Huyakorn P S, Pinder G F. Computational methods in subsurface flow[M]. New York: Academic Press, 1983.
[53] CRC Handbook of Chemistry and Physics[M]. 63rd ed. Florida: CRC Press, 1982: 261.
[54] Van Genuchten M T. A closed‐form equation for predicting the hydraulic conductivity of unsaturated soils[J]. Soil Science Society of America Journal, 1980, 44(5): 892-898.
[55] Frind E O. Simulation of long-term transient density-dependent transport in groundwater[J]. Advances in Water Resources, 1982, 5: 73-88.
[56] Celia M A, Bouloutas E T, Zarba R L. A general mass-conservative numerical solution for the unsaturated flow equation[J]. Water Resources Research, 1990, 26: 1483-1496.
[57] Balay S, Buschelman K, Eijkhout V, et al. PETSc: home page[EB/OL]. http://www.mcs.anl.gov/petsc.
[58] Knepley M G, Katz R F, Smith B. Developing a geodynamics simulator with PETSc[C]//Lecture Notes in Computational Science and Engineering. Springer Berlin Heidelberg, 2006.
[59] Sleijpen G L G, Fokkema D R. BiCGstab (l) for linear equations involving unsymmetric matrices with complex spectrum[J]. Electronic Transactions on Numerical Analysis, 1993, 1(11): 2000.
[60] Saad Y. Iterative methods for sparse linear systems[M]. Society for Industrial and Applied Mathematics, 2003.
[61] 李晓梅, 吴建平. 数值并行算法与软件[M]. 北京: 科学出版社, 2007.
[62] Stoeckl L, Houben G. Flow dynamics and age stratification of freshwater lenses: Experiments and modeling[J]. Journal of Hydrology, 2012, 458: 9-15.
[63] Stoeckl L, Walther M, Graf T. A new numerical benchmark of a freshwater lens[J]. Water Resources Research, 2016, 52(4): 2474-2489.
[64] Fulton S R, Ciesielski P E, Schubert W H. Multigrid methods for elliptic problems: A review[J]. Monthly Weather Review, 1986, 114(5): 943-959.
[65] Goswami R R, Clement T P. Laboratory‐scale investigation of saltwater intrusion dynamics[J]. Water Resources Research, 2007, 43(4): W04418.
[66] Povich T J, Dawson C N, Farthing M W, et al. Finite element methods for variable density flow and solute transport[J]. Computational Geosciences, 2013, 17: 529-549.
[67] Henry H R. Effects of dispersion on salt encroachment in coastal aquifers[J]//Seawater in Coastal Aquifers. US Geological Survey Water Supply Paper, 1964: 70-84.
[68] Segol G, Pinder G F, Gray W G. A Galerkin-finite element technique for calculating the transient position of the saltwater front[J]. Water Resources Research, 1975, 11(2): 345-346.