The processes of seasonal changes in upper layers of permafrost 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 implemented for multicore processor using OpenMP technology.
Permafrost occupies 35 million km2 (about 25% of the total land area of the
world [
We consider new three-dimensional model which allows one to describe the
heat distribution in upper layers in permafrost soils with taking into
consideration not only most signi cant climatic factors (seasonal changes in temperature
and intensity of solar radiation due to geographic location of the eld) and
physical factors (di erent thermal characteristics of non-uniform ground that
change over the time), but also engineering construction features of the
production wells, and other types of technical systems such as ripraps, tanks, pipelines,
are systems, and others [7{10]. Considered numerical algorithms are justi ed
and approved on the base of data for 12 Russian oil and gas elds. 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
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
required [
In numerical simulations with using nite-di erence 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- vediagonal 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 e ciency.
THERMAL EXCHANGE
SOLAR RADIATION
EMISSIVITY THERMAL CONDUCTIVITY
PHASE TRANSITION
However, the considered model of heat distribution is assumed a nonlinear boundary condition of heat uxes 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 e ective code implementation of a splitting method by spatial variables is presented. 2
Simulation of processes of heat distribution in the permafrost soil is reduced to
solution of a three-dimensional di usivity equation with non-uniform coe cients.
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) [
T ) with initial condition
T (0; x; y; z) = T0(x; y; z):
Here, is the density [kg/m3], T is the temperature of the phase transition [K], (1) (2) c (T ) = (T ) = c1(x; y; z); T < T ; c2(x; y; z); T > T 1(x; y; z); T < T ; 2(x; y; z); T > T is the speci c heat [J/(kg K)]; is the thermal conductivity coe cient [W/(m K)]; k = k(x; y; z) is the speci c heat of phase transition, is the Dirac delta function.
Balance of the heat uxes at the surface z = 0 de nes the corresponding nonlinear boundary conditions q + b(Tair
T (x; y; 0; t)) = " (T 4(x; y; 0; t)
Ta4ir) + : (3)
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 [
An e ective scheme of through computations with smoothing the discontinuous
coe cients was developed for the equation of thermal conductivity by
temperature in the neighborhood of the phase transformation. This scheme is applied
to numerical simulations. On the basis of ideas in [
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 di erence problems approximating (1){(3) is proved in [
field at the instant t at the mesh nodes: comp_on_mesh (59,92%)
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 Pro ler tool was used. The results presented in Fig. 3 show that the most time-costly procedures are computation of the temperature eld `comp on mesh()` and preparation of the domain parameters ` zad teplo z grunt ()` and ` zad teplo z skv ()`.
The temperature computation procedure consists of three steps of forming and solving SLAEs for each spatial direction. We can solve each SLAE independently 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 con icts.
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.
Pro ling of parallel program shows that the proportion of the parallel code is 75%.
er DATA u ttereapm tttsan XYZ(((iii,,,jjj,,,kkk)))
n teh ehi T(t,x,y,z) fo t n ta parameters ito ld au ife t p m o C for k for j n o it rce 1) forming SLAE id- 2) solving SLAE x 3) updating T for i for k n o it ce 1) forming SLAE r id- 2) solving SLAE y 3) updating T for i for j n o it ce 1) forming SLAE ird 2) solving SLAE z 3) updating T As a model problem, let consider a seasonal freezing{thawing of the upper layers of a permafrost soil. The permafrost temperature lower than the area of in uence of the seasonal changes (lower than 10 meters) is -0.7 C. The basic 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 equation (1) and conditions (2) and (3) are determined as a result of the geophysical research.
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 e ciency
Ep = Sp=p for solving the xed problem size by various numbers p of OpenMP
threads (strong scaling). It also contains the theoretical speedup calculated using
the Amdahl law [
Sp = is proportion of the serial code (25% for our program). (4)
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
Thus, the developed mathematical model and software allow one to carry out detailed numerical calculations on long-term forecasting the temperature eld changes from di erent technical systems in the near-surface layer of soil in the permafrost zone. E ective 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.
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 Program.