Simulation of Thermal Processes in Permafrost: Parallel Implementation on Multicore CPU Elena N. Akimova1,2 , Mikhail Yu. Filimonov1,2 , Vladimir E. Misilov1,2 , and Nataliia A. Vaganova1,2 1 Ural Federal University, Ekaterinburg, Russia 2 Krasovskii Institute of Mathematics and Mechanics, Ural Branch of RAS, Ekaterinburg, Russia aen15@yandex.ru, fmy@imm.uran.ru, v.e.misilov@urfu.ru, vna@imm.uran.ru Abstract. The processes of seasonal changes in upper layers of per- mafrost soil are described. A thermal conductivity equation is considered and algorithm is presented for solving it. An approach of computational optimization is suggested. The proposed parallel algorithm is implemen- ted for multicore processor using OpenMP technology. Keywords: thermal conductivity, numerical simulation, permafrost, par- allel computing, MPI, OpenMP 1 Introduction Permafrost occupies 35 million km2 (about 25% of the total land area of the world [1, 2]. In Russia, more than 60% of its territory is occupied by permafrost, and 93% natural gas and 75% oil are produced in these regions. Human ac- tivity, production and transportation of oil and gas have a significant effect on permafrost. The warm oil heated the pipes in wells and pipelines and other pro- cesses may lead to permafrost degradation. Therefore, the problem of reducing the intensity of thermal interactions in the “heat source — permafrost” zones has a particular importance for solving problems of energy saving, environmen- tal protection, safety, cost savings, and enhance operational reliability of various engineering structures [3, 5, 4, 6]. We consider new three-dimensional model which allows one to describe the heat distribution in upper layers in permafrost soils with taking into considera- tion not only most significant climatic factors (seasonal changes in temperature and intensity of solar radiation due to geographic location of the field) and physical factors (different thermal characteristics of non-uniform ground that change over the time), but also engineering construction features of the produc- tion wells, and other types of technical systems such as ripraps, tanks, pipelines, flare systems, and others [7–10]. Considered numerical algorithms are justified and approved on the base of data for 12 Russian oil and gas fields. When the source of heat in the frozen ground was a well, a comparison was made between numerical data on the distribution of the boundary of the melting of frozen 2 ground (zero isotherm) and experimental data. The detailed calculations of the long-term prediction of the permafrost boundary changes demand considerable computing power, the parallel approach to solving such problems is also re- quired [11]. The complete simulation of all technical systems located in a well pad makes it necessary to solve such problems in a significally larger area with a three-dimensional computational grid. This leads to an essential encreasing time of computations (up to 100 hours). In numerical simulations with using finite-difference methods, an original problem is often reduced to numerical solution of systems of linear algebraic equations (SLAE). In papers [12–15], for solving SLAE with the block-three- and block-fivediagonal matrices, direct and iterative parallel algorithms are designed. Those are: the conjugate gradient method with preconditioner, the conjugate gradient method with regularization, the square root method, and the parallel matrix sweep algorithm. In papers [16–18], the iterative gradient methods were used for solving SLAE with dense matrices arising in inverse gravity problems. These algorithms are implemented on multi-core processor with good efficiency. THERMAL EXCHANGE SOLAR RADIATION EMISSIVITY THERMAL CONDUCTIVITY PHASE TRANSITION Fig. 1. A basic scheme of simulated area and heat fluxes However, the considered model of heat distribution is assumed a nonlinear boundary condition of heat fluxes balance at the soil–air boundary that increases the computational comlexity and does not allow the direct application of the parallel methods of SLAE solving. In this paper, an appoach of effective code implementation of a splitting method by spatial variables is presented. 2 Statement of problem Simulation of processes of heat distribution in the permafrost soil is reduced to solution of a three-dimensional diffusivity equation with non-uniform coefficients. 3 The equation includes the localized heat of phase transition; i.e., an approach to solve the problem of the Stefan type without the explicit separation of the phase transition in 3D area (Fig. 1) [3]. The equation for the temperature T (t, x, y, z) at the instant t at the point (x, y, z) has the form  ∂T ρ cν (T ) + kδ(T − T ∗ ) = ∇ (λ(T )∆T ), (1) ∂t with initial condition T (0, x, y, z) = T0 (x, y, z). (2) Here, ρ is the density [kg/m3 ], T ∗ is the temperature of the phase transi- tion [K], c1 (x, y, z), T < T ∗ ,  cν (T ) = is the specific heat [J/(kg · K)], c2 (x, y, z), T > T ∗ λ1 (x, y, z), T < T ∗ ,  λ(T ) = is the thermal conductivity coefficient [W/(m · K)], λ2 (x, y, z), T > T ∗ k = k(x, y, z) is the specific heat of phase transition, δ is the Dirac delta function. Balance of the heat fluxes at the surface z = 0 defines the corresponding nonlinear boundary conditions ∂T (x, y, 0, t) γq + b(Tair − T (x, y, 0, t)) = εσ(T 4 (x, y, 0, t) − Tair 4 )+λ . (3) ∂z To determine the parameters in the boundary condition (3), an iterative algorithm is developed that takes into account the geographic coordinates of the area, lithology of soil, and other features of the considered location [10]. 3 Method for solving the problem An effective scheme of through computations with smoothing the discontinuous coefficients was developed for the equation of thermal conductivity by temper- ature in the neighborhood of the phase transformation. This scheme is applied to numerical simulations. On the basis of ideas in [3], a finite difference method is used to solve the problem with splitting by the spatial variables in three- dimensional domain. We construct an orthogonal grid, uniform, or condensing near the ground surface or to the surfaces of internal boundaries. The original equation for each spatial direction is approximated by a locally additive implicit pattern, and to solve SLAE, a combination of the sweep and Newton method is used. At the upper boundary z = 0, there is an algebraic equation of the fourth degree, which is solved by the Newton method. Solvability of the same difference problems approximating (1)–(3) is proved in [5]. 4 INPUT Preparation of parameters of the iterations by the time domain, constructions, and soil Computation of the temperature field at the instant t at the mesh nodes: comp_on_mesh (59,92%) OUTPUT Fig. 2. Principal computations algorithm 4 Parallel algorithm A basic scheme of the computational algorithm is presented in Fig. 2 Let us consider an approach to construction of the parallel algorithm for multicore CPU using OpenMP technology. To analyse the serial program, the Visual Studio 2017 Performance Profiler tool was used. The results presented in Fig. 3 show that the most time-costly pro- cedures are computation of the temperature field ‘comp on mesh()‘ and prepara- tion of the domain parameters ‘ zad teplofiz grunt ()‘ and ‘ zad teplofiz skv ()‘. Fig. 3. Profiler results 5 The temperature computation procedure consists of three steps of forming and solving SLAEs for each spatial direction. We can solve each SLAE inde- pendently within one direction. For example, when forming and solving SLAEs along the X axis, we can split the Y Z domain into a number of horizontal bars and distribute them between a number of threads. We do this by using ‘#pragma omp parallel for‘ for outer of two nested loops. After all SLAEs are solved for one direction, we should synchronize threads by ‘#pragma omp barrier‘ to avoid update conflicts. Figs. 4 and 5 show the serial algorithm and the parallel variant, respectively. The procedures for domain parameters consist of three nested loops; so, we can just use ‘ parallel for‘ for the outer one. Profiling of parallel program shows that the proportion of the parallel code is 75%. x-direction for k for j 1) forming SLAE 2) solving SLAE Computation of the temperature DATA 3) updating T X(i,j,k) field at the instant t Y(i,j,k) for i Z(i,j,k) y-direction for k T(t,x,y,z) 1) forming SLAE 2) solving SLAE parameters 3) updating T for i z-direction for j 1) forming SLAE 2) solving SLAE 3) updating T Fig. 4. Serial algorithm 5 Numerical experiments As a model problem, let consider a seasonal freezing–thawing of the upper lay- ers of a permafrost soil. The permafrost temperature lower than the area of influence of the seasonal changes (lower than 10 meters) is -0.7◦ C. The basic 6 parallel for k x-direction for j for j for j for j for j … 1)f 1)f 1)f 1)f 1)f Computation of the temperature DATA 2)s 2)s 2)s 2)s 2)s 3)u 3)u 3)u 3)u 3)u X(i,j,k) field at the instant t Y(i,j,k) parallel for k Z(i,j,k) y-direction for i for i for i for i for i T(t,x,y,z) … 1)f 1)f 1)f 1)f 1)f parameters 2)s 2)s 2)s 2)s 2)s 3)u 3)u 3)u 3)u 3)u parallel for i z -direction for j for j for j for j for j … 1)f 1)f 1)f 1)f 1)f 2)s 2)s 2)s 2)s 2)s 3)u 3)u 3)u 3)u 3)u Fig. 5. Parallel algorithm thermal parameters of the soil are in the following: the thermal conductivity is 1.82 and 1.58 W/(m·K), the volumetric heat is 2130 and 3140 kJ/(m3 ·K) for frozen and melted soil, respectively, the volumetric heat of phase transition is 1.384· 105 kJ/(m3 ·K). Table 1 shows the data for the considered area. The other parameters in equa- tion (1) and conditions (2) and (3) are determined as a result of the geophysical research. Table 1. Basic annual climatic parameters Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec ◦ Aver. temperature, C -28.1 -26.3 -20.0 -10.8 -2.0 8.5 15.1 11.5 5.1 -6.5 -20.8 -26.2 Solar radiation, W/m2 0 112 282 567 809 865 889 639 355 122 34 0 The experiments are carried out using the six-cores AMD Ryzen 5 1600X CPU. The grid size for problem is 91 × 91 × 51. The time interval was 10 days with the step of 1 day. Table 2 shows the computing times Tp , speedup Sp = Tp /T1 , and efficiency Ep = Sp /p for solving the fixed problem size by various numbers p of OpenMP 7 threads (strong scaling). It also contains the theoretical speedup calculated using the Amdahl law [19]: 1 Sp = , (4) α + 1−α p where α is proportion of the serial code (25% for our program). Table 2. Results of strong scaling experiments Number Calculation time Theoretical Speedup Sp Efficiency Ep of threads p Tp , sec speedup Sp 1 620 - - - 2 390 1.6 1.59 0.79 4 310 2.3 2 0.5 6 250 2.7 2.5 0.41 Table 3 shows computing times and efficiency Wp = T1 /Tp for fixed problem size per thread (weak scaling). The grid size is adjusted automatically, so the dimensions are not exactly scaled. Table 3. Results of weak scaling experiments Number Calculation time Grid size Efficiency Wp of threads p Tp , sec 1 91 × 91 × 51 620 - 2 181 × 91 × 51 820 0.76 4 379 × 85 × 51 1210 0.51 6 585 × 85 × 51 1500 0.41 By extrapolating these calculation times, we can say that the serial program would take 6.3 hours to solve the problem for 1 year interval and 188 hours (8 days) for 30 years interval. Using the parallel computing, we can cut this time to approximately 2 days. 6 Conclusion Thus, the developed mathematical model and software allow one to carry out detailed numerical calculations on long-term forecasting the temperature field 8 changes from different technical systems in the near-surface layer of soil in the permafrost zone. Effective and quite simple OpenMP approach for a multicore system allows one to reach approximately 90% of the theoretical speedup. The results of numerical experiments show that by using the six-core processor, the computation time is reduced up to 2.5 times. Acknowledgments The work was partially supported by the Russian Foundation for Basic Research 16–01–00401 and by the Center of Excellence “Geoinformation technologies and geophysical data complex interpretation” of the Ural Federal University Pro- gram. References 1. Zhang, T., Barry, R. G., Knowles, K., Heginbottom, J. A., Brown, J.: Statistics and characteristics of permafrost and ground ice distribution in the Northern Hemi- sphere. Polar Geography, Vol. 23(2), 132–154 (1999) https://doi.org/10.1080/ 10889379909377670 2. Jorgenson, M. T., Romanovsky, V., Harden J., Shur Y., O’Donnell,J., Schuur, E. A. G., Kanevskiy, M.: Resilience and Vulnerability of Permafrost to Climate Change. Canadian J. Forest Research, Vol. 40(7), 1219–1236 (2010) https://doi. org/10.1139/X10-060 3. Filimonov, M. Yu., Vaganova, N. A.: Simulation of thermal stabilization of soil around various technical systems operating in permafrost. Applied Mathemat- ical Sciences, 7(144), 7151–7160, (2013) https://doi.org/10.12988/ams.2013. 311669 4. Filimonov M. Yu., Vaganova N. A.: Simulation of Technogenic and Climatic In- fluences in Permafrost.: Lecture Notes in Computer Science. Vol. 9045. 178–185. (2015) https://doi.org/10.1007/978-3-319-20239-6_18 5. Vaganova, N. A., Filimonov., M. Yu.: Simulation of Engineering Systems in Per- mafrost. Siberian Journal of Pure and Applied Mathematics. (2013) 6. Filimonov, M., Vaganova, N.: Permafrost thawing from different technical systems in Arctic regions. IOP Conf. Ser.: Earth Environ. Sci. Vol. 72. No. art. 012006 (2017) https://doi.org/10.1088/1755-1315/72/1/012006 7. Vaganova, N., Filimonov, M.: Simulation of freezing and thawing of soil in Arctic regions. IOP Conf. Ser.: Earth Environ. Sci. Vol. 72. No. art. 012005 (2017) https: //doi.org/10.1088/1755-1315/72/1/012005 8. Vaganova, N. A., Filimonov, M. Yu.: Computer simulation of nonstationary ther- mal fields in design and operation of northern oil and gas fields. AIP Conf. Proc. Vol. 1690, No. 020016 (2015) https://doi.org/10.1063/1.4936694 9. Filimonov, M., Vaganova, N.: Short and Long Scale Regimes of Horizontal Flare System Exploitation in Permafrost. CEUR Workshop Proceedings. Vol. 1662, 253– 260 (2016) http://ceur-ws.org/Vol-1662/mod3.pdf 10. Filimonov, M. Yu., Vaganova, N. A.: On Boundary Conditions Setting for Nu- merical Simulation of Thermal Fields Propagation in Permafrost Soils, CEUR Workshop Proceedings, Vol. 2109, 18–24, (2018) http://ceur-ws.org/Vol-2109/ paper-04.pdf 9 11. Filimonov, M. Yu., Vaganova, N. A.: Parallel splitting and decomposition method for computations of heat distribution in permafrost, CEUR Workshop Proceedings, Vol. 1513, 42–49, (2015) http://ceur-ws.org/Vol-1513/paper-05.pdf 12. Akimova, E. N., Belousov, D. V.: Parallel algorithms for solving lin- ear systems with block-tridiagonal matrices on multi-core CPU with GPU, Journal of Compu- tational Science. Vol. 3, No. 6, 445449 (2012) https://doi.org/10.1016/j.jocs. 2012.08.004 13. Akimova, E. N., Belousov, D. V., Misilov, V. E.: Algorithms for solving inverse geo- physical problems on parallel computing systems. Numerical Analysis and Appli- cations, Vol. 6(2), 98–110 (2013) https://doi.org/10.1134/S199542391302002X 14. Akimova, E. N.: A parallel matrix sweep algorithm for solving linear system with block-fivediagonal matrices. AIP Conf. Proc., Vol. 1648, 850028. (2015) https: //doi.org/10.1063/1.4913083 15. Akimova, E., Belousov, D.: Parallel algorithms for solving linear systems with block-fivediagonal matrices on multi-core CPU. CEUR Workshop Proceedings, Vol. 1729, 38–48 (2016) http://ceur-ws.org/Vol-1729/paper-06.pdf 16. Martyshko, P. S., Pyankov, V. A., Akimova, E. N., Vasin, V. V., Misilov, V. E.: On solving a structural gravimetry problem on supercomputer “Uran” for the Bashkir Predural’s area, GeoInformatics 2013 - 12th International Conference on Geoinformatics: Theoretical and Applied Aspects (2013) http://www.earthdoc. org/publication/publicationdetails/?publication=68121 17. Akimova, E. N., Martyshko, P. S., Misilov, V. E., Kosivets, R.A.: An efficient numerical technique for solving the inverse gravity problem of finding a lateral density. Applied Mathematics and Information Sciences, Vol. 10(5), 1681–1688 (2016) http://doi.org/10.1063/1.4992206 18. Akimova, E. N., Martyshko, P. S., Misilov, V. E.: On finding a density in a curvi- linear layer by biconjugate gradient type methods, AIP Conference Proceedings, Vol. 1863, 050009, (2017) http://doi.org/10.1063/1.4992206 19. Voevodin, V. V., Voevodin, V. V.: Parallel Computing. BHV-Petersburg, St. Pe- tersburg (2002)