CC BY
8
Discrete & Continuous Models 2019 27 (1) 21—32
. „ . http://journals.rudn.ru/miph
& Applied Computational Science
Computational modeling and simulation
Research article
PACS 07.05.Tp, 02.60.Pn, 02.70.Bf DOI: 10.22363/2658-4670-2019-27-1-21-32
Parallel algorithm for numerical solution of heat equation in complex cylindrical domain
Alexander Ayriyan1, Jan Busa Jr.1'2
1 Laboratory of Information Technologies, JINR, Dubna, Russia 2 Institute of Experimental Physics, Slovak Academy of Sciences, Kosice, Slovakia
(received: April 6, 2019; accepted: October 21, 2019)
In this article we present a parallel algorithm for simulation of the heat conduction process inside the so-called pulse cryogenic cell. This simulation is important for designing the device for portion injection of working gases into ionization chamber of ion source. The simulation is based on the numerical solving of the quasilinear heat equation with periodic source in a multilayered cylindrical domain. For numerical solution the Alternating Direction Implicit (ADI) method is used. Due to the non-linearity of the heat equation the simple-iteration method has been applied. In order to ensure convergence of the iteration process, the adaptive time-step has been implemented. The parallelization of the calculation has been realized with shared memory application programming interface OpenMP and the performance of the parallel algorithm is in agreement with the case studies in literature.
Key words and phrases: quasilinear heat equation, multilayer cylindrical geometrical structure, pulse periodic source, parallel algorithm, thermal gates
1. Introduction
The purpose of this work is to develop algorithms for simulation of the heat conduction process inside the so-called pulse cryogenic cell [1,2]. Such simulations are important for designing the cell that implements "the thermal gates" of a portion injection of working gases into the ionization chamber of a multiply charged ion source [3]. While reliable operation of mechanical valves for pulsed injection of gaseous mixtures in the millisecond range at cryogenic temperatures is practically impossible, the use of gas temperature properties at cryogenic temperatures can be a real alternative. Indeed, the vapor pressures of various gases have strong dependency on the temperature [4] in the interval between temperatures of liquid helium 4.2 K and liquid nitrogen
© Ayriyan A., Busa Jr. J., 2019
This work is licensed under a Creative Commons Attribution 4.0 International License http://creativecommons.org/licenses/by/4.0/
78 K [5], typical temperature terminals in the cryogenic technique used as thermostats with large capacity.
The cryogenic cell is a multilayer cylinder (see Fig. 1) placed inside a vacuum chamber. The thermal process in the cylinder starts by passing an electrical current through its conductive layer. It allows heating the outer layer of the cylinder to upper desired temperature (maximal critical temperature). After switching the current off, the cylinder is let to cool down to lower desired temperature (minimal critical temperature). This process periodically repeats for a given period of time based on the requirements of the experiment. A copper core of the cylinder serves as a cooler for outer layers during and after the heating process. It is connected to a liquid helium temperature terminal. The core is separated by the electrical insulator from the conductive layer. It is made to avoid the flow of the electrical current into the core. The last layer is a thin coating which prevents molecules of working gases from binding to the conductive layer made of graphite.
Figure 1. A schematic view of the quarter of the cell slice through the axis. The bottom line is the cylinder axis (the symmetry axis, r = 0). The cooler (the copper core rod) cools the
cell by contact with the temperature terminal - 4 (liquid helium). The heater - 2 (the conductive layer) heats the cell up by the way of passing of the electrical current. The inner insulator - 1 - is needed to prevent the electrical current outflow from the heater to the cooler. The outer insulator - 3 - prevents molecules of working gases from binding to the
conductive layer
2. Initial-boundary value problem
Let us consider the heat equation describing thermal evolution in the closed
cylindrical domain Q = {(r, z) | r G [0, rmax(z)], z G [0, zmax(r)]}:
dT 1 d ( dT \ d ( dT \
p(r)c"(Tr)m = 1 dr (rA(T'r)Sir) + dZ (A(T-r) &) +
+ X(T, t, r; tper,t src / ) V /
where the thermal coefficients are the non-linear function of the temperature and they have discontinuities of the first kind along the radial direction at rm (m = 0,1, 2, 3, see Fig. 1). The source function producing the periodic
process of heating can be expressed in the form
12(r)
X r; ¿per, ¿src ) X (T) r* p(t; ¿per, ¿src )
S
C
(2)
where p(t; iper,isrc) is the periodic normalized function with parameters tper (time of period) and tsrc (time of heating less than or equal to tper); n G N0 is the index of a period; x(T, r) is temperature depended specific resistivity with discontinuities at the given values of r; I(r) is the electrical current amplitude which has a finite value I0 only in the source layer (see Fig. 1), everywhere else it is zero; SC is the cross-section of the source layer. The periodic normalized function p(t; iper,isrc) is given as a rectangular pulse one [6]:
u(t; tper,tsrc) ^ ^ \ß (t ntper) ^ (t n ¿per tsrc)]
(3)
n=0
where 0(i) is Heaviside step function [7,8], however, in order to make the processes of "turn on" and "turn off" more realistic, the model of the transient process (4) (see Fig. 2) has been implemented:
v(t; ¿per? ¿src? ¿trs? C)
E
n=0
erf( 1
t
trs
- erf ( (
L t n ¿per ¿src i ¿trs
. (4)
Here erf(t) is the error function [7,8].
nt
n-t +t
per src
Figure 2. Dashed line: the periodic step function u(t) as a combination of the Heaviside ones. Solid line: the "realistic" source function for £ = 4 and Z = 2 representing the transient model. The dot-dashed lines show the ends of the transient process at "turn on"
and "turn off"
1
2
t
Such smoothing of the processes of "turn on" and "turn off" of the source also helps to stabilize numerical simulations.
The initial condition can be set by the assumption that at t = 0 the cryogenic cell has been already cooled down by the temperature terminal:
T (r,z,t = 0) = To, (5)
where T0 = 4.2 K is the temperature of liquid helium.
The boundary conditions can be expressed after the assumption that the temperature flow across the boundary of the domain is zero except for the right side where a connection to the temperature terminal exists (see Fig. 1):
dT
— = 0 V (r,z) G in \{(r,z): z = Zmax}, (g)
T = To V (r, z) G {(r, z) : z = Zmax},
where in is the boundary of the domain n, n is the normal vector of the boundary in. This assumption is motivated by the following statements:
— the cryonics cell is installed in the vacuum chamber, therefore, there is no convective heat transfer;
— the working temperature is too low for the appearance of thermal radiation;
— we neglect the energy for evaporation of gas molecules;
— there is no temperature flow through the axis r = 0 due to the axial symmetry.
3. Discretization and numerical method
Numerical solution of Eq. (1) can be obtained by using a shifted non-uniform grid (see Fig. 3):
u = {(t,r, z) | 0 ^ t< TO, tfc+i = tk + Tfc+1, k G No;
0.5hi ^ r ^ rmax - 0.5hWj-i, ri+i = n + hi+i, i = 0,...,Nj - 1; 0.5ni ^ z ^ zmax - 0.5nMi-i, Zj+i = Zj + ^j+i, j = 0,..., Mi - 1}. (7)
r
Figure 3. The discretization of the 2D domain. The non-uniform grid is shifted in order to have points at the centers of boxes. Different layers (materials) have different step size (they
are differently colored)
z
It seems usual to make uniform steps in the subdomain corresponding one layer. The shifted grid has no points at the boundary of discontinuity. One can also use a special grid - the grid with points at the boundary of materials, in this case one has to take care of approximation of the thermal coefficients and the source function at the boundary.
The initial boundary value problem Eqs. (1)-(6) has been approximated on the grid (7) by the alternating direction implicit (ADI) schemes [9-12]:
pij CVij
pij CVij
T- ■ — T •
0.5t f--- T- ■ 0.5t
ar [ t i,j ] + az [ ti,j ]+ x ar [ ti,j ] + az [ t^ ]+ x
(8)
(9)
where Tjj is the temperature distribution on the next time layer, T^ is the temperature distribution on the half layer (in between the next and current time layers), Tj,j is the value on the current time layer and t is the time step,
Pi,j = P(Ti,j), cvij = cy (T ), Xi,j = X (T ). The spatial finite difference operators in Eqs. (8) and (9) are:
Ar j 11
\ Ti+1j Tij ri+1 Ai+1 j—h—
i+l
- ri-2 Ai-2 ,j
Ti,j Ti-1,. hi
Az [ Tij ] = -1
TT
Ai,j+2 nj+l AiJ-1
TT
(10) (11)
where i = 1,..., Nj — 1, j = 1,..., Mj — 1, h = r — ri-i, nj = zj — zj-i, hi = (hj+i + hj) /2, nj = (nj+i + nj)/2, Ajj = Am(Tjj), cyjj = cym(Tjj),
Xi,j = Xm(Tij), rj± 1/2 = (ri + ri±1)/2, Aj±i/2,j = Am(Tjj + Ti±1,j)/2, Ai,j±i/2 =
Am(Tj,j + Ti,j±i)/2.
Due to the non-linearity of Eqs. (8)-(9) (when the thermal coefficients and the source function depend on temperature) the simple-iteration method has been applied for calculation of the sought-for function on the half and next time layers. The recursive forms for (8)-(9) are expressed as follows [11-13]:
Ts+i _t
ps c s Tjj Tj,j =
pi,j cyi,j o 5t
1 1
ri hi
Ts+1 _ Ts+1
1 A?, ! /i+1,j - r,
'i+2' li+1
h
i+1
t s+1_t s+1
ij i 1 j ! A l .-;-
"1 i-2j h
+
+ Az [ Tij ]+ X?., (12)
here it is assumed that when s ^ œ, then T ? ^ T, p? ^ p, cV? ^ cV, A? ^ A, and X? ^ X.
The iteration process starts with the initial condition that Tj = Tij, and stops after fulfilling the following criteria
||Ts+1 - Ts||c = max |Ts+1 - Ts|| < e. (13)
The values of the sought-for function on next time-layer are obtained as this
rps+1 rp — - Ti,j — T i,j pi,j cVi,j o
= Ar [ ti,j ] +1 n j
Ts+1 _ Ts+1 Ts+1 _ Ts + 1
^s Ti,j+1 — Ti,j _ ^s Ti,j — Ti,j-1
i,j+1 nj+1 i,j-1 Vj
+ Xi,j, (14)
with the initial condition T^0 = Ti,j. Same as before T^ T and
^ A. This iteration process stops after fulfilling the same criteria (13). The systems of linear algebraic equations (12) and (14) are solved by the Thomas method [13-15].
In order to ensure convergence of the iteration process (12)-(14), the adaptive time-step has been implemented. If the counter of iterations s exceeds some maximal value, meaning that the process does not converge fast enough, the time-step t is divided by 2 and the iteration process is restarted (back to the Eq. (12)).
4. Parallel algorithm
Main computational complexity comes from repetitive calculation of Eqs. (12) and (14) across each time layer. A solution of the first of them is needed as a start for solving the second one and it is again used as a start for calculation of the sought-for function on the next time layer. Due to rather low complexity of solving of one system of linear equations ((12) or (14)) it is better to use parallelization based on shared memory since for distributed memory parallelization the cost of data transfer would be too high. Therefore, we opted for OpenMP [16,17].
In Fig. 4 one can see the flowchart of the algorithm. After initializing the solution - setting tprd, tsrc, I, and t0 = 0, the program repeats iterations until the requested time is reached. In one step of evolution it first initializes the estimate of the time step t as being equal ttrs/1000 if we are in the transition process or tsrc/100 elsewhere. After this it starts the iteration process in order to obtain the solution Ti,j at the time tk + t/2 using (12). If the obtained solution is precise enough, it alternates direction and continues by the way of (14) to the solution Ti,j at the time tk + t. If this solution is again precise enough, the algorithm sets the total time to tk = tk-1 + t and the actual solution Ti,j = T^j. If in any of the previous tests the number of iterations exceeds the number of maximum iterations (iter - see Fig. 4), then the timestep t is divided by 2 and the calculation returns to the beginning of the evolution step.
The whole algorithm terminates when the total time t reaches the desired value or when we realize that we have entered into the periodic process (temperature changes periodically depending on switching the current in conductive layer on and off).
Load iprd, iSrc, I] Initialize solution Tj Set k = 0,tk=0 ''
>
P p
c-i
M ¡=
ux p
c-l
Tl
p
p sr
leg
Figure 4. Schematics of the algorithm with parts running using OpenMP highlighted by the dashed line
to
5. Numerical results
The calculations have been performed on the HybriLIT computational platform [18, 19] using the Skylake processor Intel Xeon Gold 6154 [20] containing 18 CPU cores providing two threads per core (36 threads in total) under OS Scientific Linux 7.5 (Nitrogen) [21].
As an example we have studied the case when tper = 10-1, tsrc = 10-2, ttrs = 10-4, and 10 = 0.5742. These values of parameters were chosen to realize the "thermal gates" for maximal and minimal critical temperatures 42K and 37K, respectively [22]. The cell size was selected as follows: z0 = 4, zmax = 5, r* = 0.24, r* = 0.245, r* = 0.25, rmax = 0.2501. The whole domain was split into 100 parts along the z axis at the first layer (core) and to 80 in the other layers. Along the r axis, individual layers were divided (starting from the core) into 800, 200, 200, and 10 parts respectively. This discretization of the domain was chosen in order not to split the domain into too many parts and at the same time have enough information about the solution. It is noteworthy that we have much more steps in the radial direction. It is so because we expect the flow in-between the layers to be more active than the relaxation towards the terminal. Nevertheless, in our experiments we have also densified the grid to see the impact of the grid on both the precision of the solution and the calculation time. The results are given in Fig. 5, it shows the temperature fields at different moments of evolution inside the required steady periodic temperature regime.
Figure 5. Temperature fields when the periodic temperature regime is achieved. Left panel: the temperature field at the moment before turning on the source. Right panel: the temperature field at the moment when the source was just turned off
The field is shown in the whole cylindrical domain right before turning on the source (at the left panel). The field is in a state of maximum relaxation in that moment. At the right side the temperature filed is shown in the moment when source is turning off. The results demonstrate that temperature significantly changes only in outer layers, the core rode is practically not heated, it works as thermostat. The evolution of the temperature at the surface of the cryogenic cell at the z = 0 is given in Fig. 6, it is evident that is takes time to achieve the required steady periodic temperature regime.
i
it
ffi
66.5 67 67.5
> £
P E3
M
w< P
Tl
p
p^ fT
Figure 6. Temperature evolution on the surface at the £ = 0 (black line). Red colored upper dashed line represents the maximal critical temperature (42 Iv) and blue colored lower dashed line represents the minimal critical temperature (37Iv). Left panel: Temperature dependence on the time from 0 up to 70 sec. Central panel: the beginning of the evolution (the first 3.5 sec.). Right panel: the last 3.5 sec. of the evolution
when the periodic steady temperature regime is achieved
oq o
3 6 9 12 15 18 21 24 27 30 33 36 omp threads
U U0 3 6 9 12 15 18 21 24 27 30 33 36 omp threads
20
18
16
,— 14
13 12
10
(l>
6 8
H 6
4
2
0
omp threads
Figure 7. The experimental performance of the algorithm. Left panel: the calculation time in hours. Central panel: the speed up
of calculations. Right panel: the efficiency of parallelization
I.O CO
6. Summary and conclusions
We have developed the algorithm for numerical simulations of heat evolution inside the multilayered cylindrical domain with the periodic source. One can see from Fig. 6 that the required periodical temperature regime cannot be obtained immediately, it has a setup mode which is much longer than one period, thus it has to be taken into account when designing the pulse cryogenic cell. The simulations show the possibility of the realization of "thermal gates" for a particular set of parameters. The algorithm has been integrated to the hybrid algorithm MPI+OpenMP for solving the optimization problem of the heat source characteristics (tper, tsrc, and I0) of the pulse cryogenic cell [22].
The performance of the parallel algorithm (see Fig. 7) is in agreement with the case studies in literature, e.g., [23,24]. As it is shown in the picture that the saturation calculation time been achieved at 18 threads, after acting the hyper-threading the speedup stops. Thus there is no reason to involve in calculations more than this number of threads for the considered problem with the given grid size (1210 x 100).
Acknowledgments
Authors thank H. Grigorian, E. E. Donets and E. Ayryan for fruitful discussions and help, it has allowed to improve the article. The work has been partially supported by the Russian Foundation for Basic Research under the projects № 18-51-18005 and № 19-01-00645. Computations have been conducted on the basis of the HybriLIT heterogeneous computing platform (LIT, JINR).
References
1. D. E. Donets, E. E. Donets, T. Honma, K. Noda, A. Y. Ramzdorf, V. V. Salnikov, V. B. Shutov, E. D. Donets, Physics research and technology developments of electron string ion sources, Review of Scientific Instruments 83 (2012) 02A512. doi:10.1063/1.3678660.
2. K. Katagiri, A. Noda, K. Suzuki, K. Nagatsu, A. Y. Boytsov, D. E. Donets, E. D. Donets, E. E. Donets, A. Y. Ramzdorf, M. Nakao, S. Hojo, T. Wakui, K. Noda, Cryogenic molecular separation system for radioactive 11C ion acceleration, Review of Scientific Instruments 86 (12) (2015) 123303. doi:10.1063/1.4937593.
3. A. Ayriyan, J. Busa Jr., E. E. Donets, H. Grigorian, J. Pribis, Algorithm and simulation of heat conduction process for design of a thin multilayer technical device, Applied Thermal Engineering 94 (2016) 151-158. doi:10.1016/j.applthermaleng.2015.10.095.
4. R. E. Honig, H. O. Hook, Vapor pressure data for some common gases, David Sarnoff Research Center, RCA, Princeton, 1960.
5. D. B. Chelton, D. B. Mann, Cryogenic data book (1956).
6. P. Symons, Digital waveform generation, Cambridge University Press, New York, 2013.
7. M. Abramowitz, I. A. Stegun, Handbook of mathematical functions, Tenth Edition, Dover Publications, New York, 1972.
8. G. A. Korn, T. M. Korn, Mathematical handbook for scientists and engineers, Revised Edition, Dover Publications, New York, 2013.
9. D. W. Peaceman, H. H. Rachford Jr., The numerical solution of parabolic and elliptic differential equations, Journal of the Society for Industrial and Applied Mathematics 3 (1955) 28-41.
10. N. N. Yanenko, Fractional step methods for solution of multidimensional problems of mathematical physics [Metod drobnykh shagov resheniya mnogomernykh zadach], Nauka, Moscow, 1967, in Russian.
11. A. A. Samarskii, P. N. Vabishchevich (Eds.), Additive difference schemes for problems of mathematical physics [Additivnyye skhemy dlya zadach matematicheskoy fiziki], Nauka, Moscow, 1999, in Russian.
12. A. A. Samarsky, The theory of difference schemes, Marcel Dekker Inc., New York, 2001.
13. N. N. Kalitkin, Numerical methods [Chislennyye metody], Second Edition, BHV-Peterburg, Saint Petersburg, 2013, in Russian.
14. L. H. Thomas, Elliptic problems in linear differential equations over a network, Tech. rep., Columbia University, New York (1949).
15. A. A. Samarskii, P. N. Vabishchevich, Computational heat transfer, Vol. 1, John Wiley & Sons Ltd., Chichester, England, 1995.
16. L. Dagum, R. Menon, OpenMP: an industry standard API for shared-memory programming, Computational Science & Engineering, IEEE 5 (1) (1998) 46-55.
17. OpenMP Architecture Review Board, OpenMP application programming interface version 5.0 (November 2018).
URL https://www.openmp.org/wp-content/uploads/OpenMP-API-Specification-5.0.pdf
18. HybriLIT group, Heterogeneous platform "HybriLIT" [cited 14.01.2019]. URL http://hlit.jinr.ru/en/
19. G. Adam, M. Bashashin, D. Belyakov, M. Kirakosyan, M. Matveev, D. Podgainy, T. Sapozhnikova, O. Streltsova, S. Torosyan, M. Vala, L. Val-ova, A. Vorontsov, T. Zaikina, E. Zemlyanaya, M. Zuev, IT-ecosystem of the HybriLIT heterogeneous platform for high-performance computing and training of IT-specialists, in: V. Korenkov, A. Nechaevskiy, T. Zaikina, E. Mazhitova (Eds.), CEUR Workshop Proceedings, Vol. 2267, 2018, pp. 638-644.
20. Intel corporation, Intel Xeon Gold 6154 Processor [cited 14.01.2019]. URL https://www.intel.com/content/www/us/en/products/ processors/xeon/scalable/gold-processors/gold-6154.html
21. Fermilab, CERN, Scientific Linux [cited 14.01.2019]. URL http://scientificlinux.org
22. A. Ayriyan, J. B. Jr., H. Grigorian, E. E. Donets, Solving the optimization problem for designing a pulsed cryogenic cell, Physics of Particles and Nuclei Letters 16 (5) (2019 (Accepted to print)) 300-309. doi:10.1134/S1547477119030026.
23. V. A. Tokareva, O. I. Streltsova, M. I. Zuev, Parallel realizations of locally one-dimensional difference schemes for solving the initial-boundary value problems for the multidimensional heat equation, in: A. V. Friesen, V. Chudoba (Eds.), Proceedings of the XX International Scientific Conference of Young Scientists and Specialists (AYSS-2016), 2016, pp. 142-147.
24. M. D. Schuster, The heat equation: high-performance scientific computing case study, Computing in Science & Engineering 20 (5) (2018) 114-127. doi:10.1109/MCSE.2018.05329820.
For citation:
A. Ayriyan, J. Busa Jr., Parallel algorithm for numerical solution of heat equation in complex cylindrical domain, Discrete and Continuous Models and Applied Computational Science 27 (1) (2019) 21-32. DOI: 10.22363/26584670-2019-27-1-21-32.
Information about the authors:
Alexander S. Ayriyan (Russian Federation) — researcher of the Laboratory of Information Technologies, Joint Institute for Nuclear Research (e-mail: [email protected], phone: +7(496)216-35-98, ORCID: https://orcid.org/0000-0002-5464-4392, ResearcherID: O-5259-2015, Scopus Author ID: 24528624400) Jan Busa Jr. (Russian Federation) — PhD in Mathematics, Senior Researcher of the Laboratory of Information Technologies, Joint Institute for Nuclear Research (e-mail: [email protected], phone: +7(496)216-35-98, ORCID: https://orcid.org/0000-0003-3807-0373, ResearcherID: AAB-6463-2019, Scopus Author ID: 57188571341)
