=Paper=
{{Paper
|id=Vol-2274/paper-01
|storemode=property
|title=Simulation of Thermal Processes in Permafrost: Parallel Implementation on Multicore CPU
|pdfUrl=https://ceur-ws.org/Vol-2274/paper-01.pdf
|volume=Vol-2274
|authors=Elena N. Akimova,Mikhail Yu. Filimonov,Vladimir E. Misilov,Nataliia A. Vaganova
}}
==Simulation of Thermal Processes in Permafrost: Parallel Implementation on Multicore CPU==
https://ceur-ws.org/Vol-2274/paper-01.pdf
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)