Mathematical models of migration processes that take into account non-local effects caused by media's fractal properties often have an integro-differential nature. Numerical methods for solving problems for such models have a higher order of computational complexity compared to the corresponding classical methods. Therefore, for their effective practical application, the usage of high-performance computational techniques, particularly for shared memory systems, is critical. In this regard, here we study the efficiency of using a multi-threaded implementation of the TFQMR iterative algorithm for solving linear systems that arise after the discretization of the initial-boundary value problems for the non-isothermal fractionaldifferential model of moisture transport in combination with the dynamic change of time step length based on the convergence characteristics of the TFQMR algorithm. The considered computation procedure is aimed at increasing the simulation speed without explicit consideration of the features of problem being solved. The conducted studies showed that the consideration of the temperature field, which is described by an integer-order differential model, leads to a decrease in the maximum acceleration of numerical scheme's multithreaded implementation. It also leads to 8-10-times increase in simulation time due to the need to reduce time step length in accordance with different speeds of heat and mass transport processes. At the same time, the procedure for dynamical change of time step length allows performing adaptive solution of the problem without user intervention.
When taking into account non-local effects caused by medium’s fractal properties in mathematical models, it is often necessary to switch from differential to integro-differential models, particularly the models that contain the so-called fractional derivatives [1]. Numerical methods for solving problems for this class of models, such as the finite difference method [2] or spectral methods [3], have a higher order of computational complexity compared to the corresponding classical methods that are applied to integer-order differential models. Therefore, for their effective practical application, the development of high-performance computational techniques [4, 5], particularly for computational systems with shared memory, is highly important.
Two factors are critical for the performance of computational algorithms and the accuracy of the obtained solutions when using finite-difference approximation of differential equations of integer and fractional order. The first factor is the performance of algorithms for solving systems of linear algebraic equations obtained after discretization and the second one is an approach for choosing discretization steps with respect to the time and space variables that can be either fixed or dynamically changeable.
Noting that the convergence of iteration algorithms for linear systems’ solution can depend on condition numbers of systems’ matrices, here we study the efficiency of a multi-threaded implementation of the iterative TFQMR algorithm [6] for solving linear systems in combination with a dynamic time step length change based on the convergence characteristics of this algorithm. Such approach does not take into account the properties of the underlying model and the behavior of its execution time along with its influence on solutions' accuracy should be studied for each particular model or class of models to which it is applied.
The class of models we refer to in this study in the models of moisture transport. They form the basis for simulation of physical, chemical, and biological processes in soils and the results of such simulations and predictions can serve as inputs for making management decisions in agriculture and land reclamation.
Among the areas of moisture transport modeling application, for which the accuracy of forecasts and the speed of their obtainment are critical, we can highlight irrigated agriculture. In this field, in order to support decision-making [7], models based on Richards differential equation [8] are widely used. Among large literature dedicated to the increase of the accuracy of modeling based on such an equation, primarily by involving additional factors that influence the movement of moisture, we can single out two directions: the use of fractional-order models [9-12] to take into account the fractality of soils [9, 13] and the consideration of the influence of such factors as non-isothermality [14] or salinity [15] on soils’ hydro-physical properties. The combination of these two directions is relevant for further development of mathematical basis and software tools for modeling migration processes in soils.
In our previous study [12] the multi-threaded TFQMR with dynamic time step change procedure was applied to the case of isothermal fractional-differential moisture transport equation. The performance of the parallel computational scheme on different CPUs was studied without taking into account the behavior of dynamic time step change procedure. In the continuation of the studies presented in [12] the research question answered in this paper is how making the underlying model more complex combining fractional and integer-order equations will influence the acceleration of computation and will the dynamic change of time step in such case influence the accuracy of numerical method used for discretization. 2.
We consider a multi-threaded implementation of the TFQMR algorithm [6] using the OpenMP
framework that is based on the following [12]: all computations are performed in one parallel block;
all summations are implemented in parallel using the reduction directive; scalar variables are updated
in the single environment. The algorithm for solving a linear system Ax = b, dim A = N N (vectors
and matrices are here and further denoted in bold) considering as a black box the procedure of
multiplication of the left-hand side matrix on vectors can be denoted as follows ( x0 is the initial value
of solution vector):
(
(
d, xm = xm-1 + d using parallel for; (d5) calculate rv = (b − Axm ) (b − Axm ) / N using parallel for with reduction directive, calculating square root in the single environment;
(d6) algorithm execution stops if rv 1 where 1 is the accuracy threshold; (e) calculate 1 = r w using parallel for with reduction directive; (f) set = 1 / , = 1 in the single environment; (g) set y0 = w + y1 using parallel for; (h) set v Ay0 + ( Ay1 + v ) using parallel for.
We apply the TFQMR algorithm in the case when an initial-boundary value problem for a differential or integro-differential model is discretized by a finite-difference technique and is solved using a time-stepping scheme. Thus, an initial value x0 is set to a vector generated from problem’s solution on the previous step or from initial conditions.
On each time step, after solving the linear system we propose to adjust the length of time step according to the technique described in [12]. It is based on the hypothesis about the existence of a correlation between time step length and the condition number of the matrix together with the correlation of the condition number with the number of iterations of the solution algorithm.
Thus, the step is multiplied by the fixed value (which was assumed to be equal to 1.25) when the number of iterations of the TFQMR algorithm exceeds the specified maximum value (which was assumed to be equal to 30). The solution at the corresponding time step after the decrease of its length is repeated. If the number of iterations is less than 1/3 of the maximum value, the length of the next time step increases similarly.
The influence of temperature in the fractional-differential with respect to spatial variables generalization of the Richards equation [9, 10, 12] can be described according to [14] in the form of a dependency between the hydraulic conductivity k and the temperature T :
k (H ,T ) = k (H ,T0 )ea(T −T0 ) C(h) = Dx (kx (H ,T )
) + Dz (kz (H ,T ) where T0 and a are the model’s parameters. Adding the heat transport equation and such a dependency to the model described in [12], we obtain the following non-isothermal fractional-order moisture transport model: h H H z
) − S, x H z t C
T T t H x
2T = x2 + 2zT2 − cv vx Tx + vz Tz ,
0 x Lx , 0 z Lz , t 0, 0 < , < 1,
vx = kx (H ,T )
, vz = kz (H ,T )
, Dx H =
1
2 (Dx,l H + Dx,r H ),
(
The boundary conditions for equation (
The water retention curve of the soil is described according to van Genuchten model [17] in the form (h) = 0 + 1 − 0 1 + (10 h )n 1−1/n and the dependency between the hydraulic conductivity and full moisture potential is described according to Averyanov model [18] in the form
(H − z) − 0 k (H ) = k f
1 − 0 where k f is the saturated hydraulic conductivity, is the fixed exponent. second order condition
Regarding the heat transport equation (
1 t Ta (t) = Tn + (Td − Tn )1 + sin − 6 / 12 .
2 3600 The initial condition has the form T (x, 0) = Ta (0) .
As an example of practical application of the model (
The numerical solution of the initial-boundary value problem for the model (
j−1
= (xi = ihx , zk = khz ,t j = l=0 l ) : i = 1, m, k = 1, n, j = 0,1, 2,... (
Aˆ1j,i−,1k = Bˆ1j,i−,1k = 1
− − h2 x 1
− + h2 x 4 4 1 4 1 4 c k ( Hi,jk−1,Ti,jk−1 ) v c k ( Hi,jk−1,Ti,jk−1 ) v B1j,i−,1k =
Dx (k (Hij+−11,k ) + k (Hi,jk−1)), B2j,−i,1k =
D (k (Hi,jk−+11) + k (Hi,jk−1)),
z ( Hij+−11,k − Hij−−11,k ) , Aˆ2j,−i,1k = ( Hij+−11,k − Hij−−11,k ) , Bˆ2j,−i,1k =
Dx = h−1−
x (2 − )
, Dz = C(hi,jk−1) 1 4 1 4 1
− − h2 z 1
− + h2 z h−1−
z (2 − ) , Rˆ j−1 = − i,k , C c k ( Hi,jk−1,Ti,jk−1 ) v c k ( Hi,jk−1,Ti,jk−1 ) v 4 4 ( Hi,jk−+11 − Hi,jk−−11 ) , ( Hi,jk−+11 − Hi,jk−−11 ) , 1 T − 2
2 hx
1
+ 2 ,
hz
where hx = Lx / m, hz = Lz / n are the steps with respect to the spatial variables, j , j = 0,1, 2... are the
time steps’ lengths. Here and further, approximation of the function h and, similarly, other functions,
on the grid (
On the grid (
Bˆ1j,i−,1k + Ti,jk +1 Bˆ2j,−i,1k − Ti,jk Rˆi,jk−1 = ˆ ij,−k1 A1j,i−,1k =
D (k (Hij−−11,k ) + k (Hi,jk−1)), A2,i,k
j−1 = x
D (k (Hi,jk−−11) + k (Hi,jk−1)),
z Ri,jk−1 = A1j,i−,1k + A2j,−i,1k + B1j,i−,1k + B2j,−i,1k +
ij,−k1 = −hij−−11,k A1j,i−,1k − hi,jk−−11 A2j,−i,1k − hij+−11,k B1j,i−,1k − hi,jk−+11 B2j,−i,1k + +( A1j,i−,1k + A2j,−i,1k + B1,i,k 2,i,k j−1 + B j−1 −
C(hi,jk−1)
) hi,jk−1 − Sij,k − ij,−k1 + Dxhx (k (Hij,k−1+1) − k (Hij,k−1)), ˆ ij,−k1 = Ti,jk−1 CT .
Here j−1 describes the discretization of the "non-local" part of the fractional derivative and has i,k the form j−1 = i,k 2 l=2,li 1 m−1 ( i − l + 1 1− − i − l 1− ) h j−1 + xx,l,k 2 l=2,lk 1 n−1 ( k − l + 1 1− − k − l 1− ) hzjz−,i1,l , h j−1 xx,i,k
= 2Dx ( hij−−11,k A1j,i−,1k + hij+−11,k B1j,i−,1k − hi,jk−1 ( A1j,i−,1k + B1j,i−,1k )), h j−1 zz,i,k
= 2Dz (hi,jk−−11 A2j,−i,1k + hi,jk−+11B2j,−i,1k − hi,jk−1( A2j,−i,1k + B2j,−i,1k ) + 0.5hx (k (Hi,jk−+11) − k (Hi,jk−1)). 3.2.
Before the start of the solution procedure on one time step we perform parallel calculation of the values of hxjx−,1l,k , h j−1 zz,i,l
with subsequent parallel calculation of the values of A1j,i−,1k , A2j,−i,1k , B1j,i−,1k , B2j,−i,1k ,
Ri,jk−1 , ij,−k1 , and ˆ ij,−k1 . Calculation of j−1 here is performed using the previously calculated values of
i,k
hxjx−,1l,k , h j−1 . Multiplication on the left-hand side matrix of the system (
We estimate the working time of the parallel computing scheme assuming that for the time t1 (P)
spent on executing barrier synchronizations and non-parallelized blocks the estimate t1 (P) = k2 P
holds. We also assume that in the software implementation there are additional sequential operations
applied to the entire finite-difference grid. Then, the total working time can be estimated as
(
Ts (n, m, P) = Nit k1 P + k2P + k3 n m (
The initial data described in [16] were used to perform the following computational experiments. A single-layer soil model with a filtration coefficient of 15 cm/day was considered. The simulation domain was of 10 m in length and 1 m in depth. Irrigation was simulated when the average moisture content of the 0.5 m layer of soil decreased to 95% of field capacity and continued until it reached 100% of field capacity. Evapotranspiration was considered equal of 5.1 mm/day. Within the simulation domain, we model the placement of 15 drip pipelines at the depth of 0.2 m directly under 15 rows of plants.
The values of the model’s (
The dynamics of the average moisture content of the root layer for the cases of the classical
isothermal model containing only Equation (
In the case of the dynamics of the average moisture content in the root layer, the solutions obtained using the fractional-differential models simulate its slower drying compared to the classical models. Consideration of the influence of soil temperature on hydraulic conductivity, on the contrary, leads to faster drying. Thus, the considered factors have an opposite in direction, but identical in behavior, effect on the integral indicator - the dynamics of average moisture content (Fig. 1).
In the case of the changes in moisture content along the depth of soil massif (Fig. 2), their influence has different behavior. In comparison with the solutions obtained by the classical isothermal model, the consideration of the factor of spatial non-locality leads to the solutions with a decreased moisture content inside the zone moistened by emitters (0.08<z<0.3 in Fig.2) and an increased moisture content at the borders of these zones (z<0.08 and 0.3<z<0.5). Moreover, a greater influence is observed in the deeper layers of soil. Taking into account the non-isothermal nature of the process leads to an increase in the rate of moisture transport and an increase in moisture content in the entire moistened zone depending on pressure gradient. The new non-isothermal fractional-differential model describes moisture distributions that are a superposition of these influences.
0,25 r e y a lt o reo 0,245 h t n it n tnoe 0,24 c e itro 0,235 u s m e g a r veA 0,23 0,225
Classical isothermal model Fractional-differential isothermal model Fractional-differential non-isothermal model 0 10000 20000 30000 40000 50000 60000 70000 80000 Time, s
The considered problem was solved with fixed time step lengths ( = 1 , 5, 20, 30, 40, and 60) as well as with the application of the procedure for their dynamic change with the maximum number of TFQMR iterations equal to 10 and 30. The execution of the iterative procedure ended when residual becomes less than 10-12. Calculations were carried out until t=7000 on 50x500 cells grid in 4 threads. Irrigation was simulated in the whole considered time period to exclude the influence of discontinuities of the model’s coefficients. The solution at = 1 was used as a base for comparison. As an accuracy criterion we used the average absolute difference between the solutions at a fixed moment of time (here - t=7000) in the form (h) = 1 nm n m (hi(j1) − hij ) i=1 j=1 2 where hi(j1) is the water head in the cell (i,j) calculated by the numerical method with = 1 , hij is the water head in the cell (i,j) for the numerical solution h , the accuracy of which is assessed.
Accuracy estimates in the performed computational experiments linearly depended on the time step length (Fig. 3, for the procedure of time step’s dynamic change here and in Fig. 4 its average value was used) in accordance with theoretical expectations arising from the used first-order approximation of the derivative with respect to the time variable. The time spent on obtaining the solution of the problem at one time step (Fig. 4) also increased linearly with the increase in the length of the step. In this case it can be explained by the linear increase of the number of TFQMR iterations required to solve the linear systems. Let us note that an increase in the number of iterations for the TFQMR algorithm indicates, according to the results of multiple studies (see, in particular, [20]), an increase in the condition number of the matrices.
When the procedure of the dynamic change of time step with a restriction on 30 TFQMR iterations was used, the average step length was equal to 27.5 with the execution time of 6.14 s. With a restriction on 10 iterations, the average step length was equal to 9 and the execution time was equal to 14.95 s. The considered solutions’ accuracy estimate when using the average step length corresponded to the linear dependency on the step length.
Thus, the procedure of the dynamic change of time step has no influence on the order of accuracy of the numerical method. It increases the accuracy with the decrease of the restriction on the number of iterations of the TFQMR algorithm.
6,00E-06 5,00E-06 0 0 0 7 t=4,00E-06 r o f e c ren3,00E-06 e iff d te2,00E-06 u l o s b ea1,00E-06 g a r e vA0,00E+00 0,016 0,015 s ,p 0,014 e t se 0,013 itm 0,012 e fon 0,011 o e 0,01 m itn 0,009 tuo 0,008 i c xeE 0,007 0,006 0 10 20
30
The total time of modeling procedure execution on the 100x1000 cells grid up to the moment of time t = 86400 for the classical and fractional-differential models in isothermal and non-isothermal cases when running on AMD Ryzen 3 5300U CPU are given in Table 1.
For all models, the running time when executed in 8 threads was less than the time observed in the case of execution in 4 threads. This reflects the fact that the used CPU has 4 physical cores. The maximum speed-up of calculations was 1.37 for the integer-order models, 2.02 for the fractionaldifferential model without taking into account the influence of temperature, and 1.72 for the nonisothermal fractional-differential model. The higher acceleration in the case of the fractional-order models is explained by the need to calculate the value of fractional derivative’s "non-local" part present in the right-hand side of the linear systems, in particular, by the fact that this block of calculations is parallelized by data without additional synchronization or reduction operations.
Linear system’s size in the case of the non-isothermal models is twice as large as in the case of the models that do not take the temperature factor into account. Also, different speeds of the heat and mass transport processes cause an increase in the condition number of linear systems obtained when 1 2 4 8 n 100 150 200 250 300 106 106 80 93 22 152 432 2044 4633 16 101 288 1183 2504 discretizing the non-isothermal models. This, in turn, leads to a decrease in the length of time steps when applying the procedure for their dynamic change. These two factors lead to a significant increase in execution times (~10 times in the case of the classical model and ~8 times in the case of the fractional-order model). The acceleration, which here is influenced by matrices structure, does not change significantly for the classical model and decreases in the case of the fractional-differential model, since, with a doubling of the size of the matrix, the time spent on the calculation of the “nonlocal” part of the fractional derivative remains unchanged.
Total execution times of the simulation procedure up to the moment of time t = 1000 for
n = 100,150, 200, 250,300, m = 10n and the case of the fractional-differential non-isothermal model
(
Using the least squares minimization technique, we obtained the coefficients’ values for the
performance model (
Let us note that the average value of the time step here decreases linearly with the decrease in the length of the step with respect to the spatial variables.
The development of the fractional-differential model of moisture transport in the direction of taking into account the influence of temperature on the relevant processes allows describing nonclassical distributions of the water head field, in particular, when applied to the modeling of subsurface drip irrigation. Compared to both integer- and fractional-order isothermal models, the new model describes situations of increased moisture content at the bottom of zones moistened by the emitters, which is caused by temperature field gradient.
The use of parallel algorithms is important for fractional-differential models as in this case the order of computational complexity increases in comparison with integer-order models. However, the complexity of the considered non-isothermal model affects the characteristics of parallel algorithms for solving the respective initial-boundary value problem.
The studied computational procedure, aimed at increasing the simulation speed without considering the features of the underlying model, consists of a combination of a multi-threaded parallelized version of the TFQMR algorithm and a procedure for the dynamic change of time step length based on the convergence characteristics of the TFQMR algorithm. The conducted studies showed that the consideration of temperature field’s influence leads to 8-10 times increase in simulation time due to the need to reduce time step lengths in order to ensure the specified accuracy of linear systems’ solution with a given restriction on the number of TFQMR iterations. The use of the integer-order heat transport equation together with the fractional-differential equation of moisture transport leads to a decrease in the maximum acceleration of the multi-threaded implementation of the numerical scheme.
The procedure for the dynamic change of time step length, the application of which does not affect the order of accuracy of the numerical method used for the discretization of the initial-boundary problems, together with the obtained estimates of multi-threaded algorithm’s performance, make possible further development of a computationally efficient decision support system for modeling moisture transport under abnormal conditions.
[12] V. Bohaienko, A. Gladky, Multithreading performance simulating fractional-order moisture transport on AMD EPYC, Journal of Numerical and Applied Mathematics 2 (2022) 174-182. doi:10.17721/2706-9699.2022.2.20 [13] Y. Pachepsky, D. Timlin, Water transport in soils as in fractal media, Journal of Hydrology 204 (1988) 98-107. doi:10.1016/S0022-1694(97)00110-8 [14] H. Stoffregen, G. Wessolek, M. Renger, Effect of temperature on hydraulic conductivity, in: M. van Genuchten et al. (eds.), Characterization and Measurement of the Hydraulic Properties of Unsaturated Porous Media. Proceedings, International Workshop on Characterization and Measurement of Hydraulic Properties of Unsaturated Porous Media, Riverside, California, 1999, pp. 497–506. [15] L.J. Chen, Q. Feng, F.R. Li, Ch.Sh. Li, Simulation of soil water and salt transfer under mulched furrow irrigation with saline water, Geoderma 241-242 (2015) 87-96. doi:10.1016/j.geoderma.2014.11.007 [16] M. Romaschenko, V. Bohaienko, A. Bilobrova, Two-dimensional mathematical modeling of the soil water regime under drip irrigation (in Ukrainian), Bulletin of Agrarian Science, 99 4 (2021) 59-66. doi:10.31073/agrovisnyk202104-08 [17] M. van Genuchten, A closed-form equation for predicting the hydraulic conductivity of unsaturated soils, Soil Sci. Soc. Am. J. 44 (1980) 886-900. doi:10.2136/sssaj1980.03615995004400050002x [18] S. Averyanov, Filtration from canals and its influence on the groundwater regime (in
Russian), Kolos, Moscow, 1982. [19] A. Samarskii, The theory of difference schemes, CRC Press, New York, NY, 2001. [20] C.H. Ahn, W.C. Chew, J.S. Zhao, E. Michielssen, Numerical Study of Approximate Inverse Preconditioner for Two-Dimensional Engine Inlet Problems, Electromagnetics, 19 2 (1999) 131-146. doi:10.1080/02726349908908631