A New Parallel Code Based on a Simple Piecewise Parabolic Method for Numerical Modeling of Colliding Flows in Relativistic Hydrodynamics

: A new parallel code based on models of special relativistic hydrodynamics is presented for describing interacting ﬂows. A new highly accurate numerical method is considered and ver-iﬁed. A parallel implementation of the method by means of Coarray Fortran technology and its efﬁciency are described in detail. The code scalability is 92% on a cluster with Intel Xeon 6248R NKS-1P with 192 Coarray Fortran images. Different interacting relativistic ﬂows are considered as astrophysical applications.

Godunov's method, which was developed almost 70 years ago [24], is used in most software codes to simulate hydrodynamic (in particular, relativistic) flows.It has been studied in detail for convergence [25] and the behavior of entropy [26,27].It has been modified to take into account the transfer of matter [28] and extended to solve equations of relativistic hydrodynamics [29], in particular for magnetic field flows [30,31].Almost immediately after the invention of Godunov's method, works began to appear on simplified solutions of the Riemann problem (for instance, Rusanov's scheme [32]) and decreased dissipation methods (for instance, Kolgan's scheme [33]).In 2000, a scheme called MUSCL was proposed [34].It combines the advantages of the Rusanov and Kolgan schemes.The next step in decreasing the dissipation of Godunov-type numerical methods was the use of piecewise parabolic [35,36] and piecewise cubic [37] presentations of solutions.A piecewise cubic presentation of the solution was extended by the authors to HLL-type methods [38] and an operator separation method [39].The separation of operators and the use of a combination of Godunov-and Rusanov-type schemes made it possible to formulate efficient numerical methods [40,41] and solve the problem of sound point in Godunov's method with linearized discontinuity breakdown [42].There is a tendency towards increasing the order of the polynomial for presenting the physical variables to be reconstructed in Godunov's scheme [43]).Nevertheless, there are many original works on using Kolgan's approach to decrease the dissipation in solving the Riemann problem [44,45].The present paper is essentially an extension of paper [42], using a new numerical method for solving the equations of special relativistic hydrodynamics and adding the degenerate Rusanov method with velocities limited to the speed of sound [46].This approach will allow simulating relativistic flows in highly rarefied zones, which are often present in relativistic jets [47].
One of the major advantages of the MPI 3.0 standard is an efficient implementation of one-way communications, which allows using one of the most promising parallel programming models, called Partitioned Global Address Space (PGAS).Computational experiments with continuum mechanics models have shown [48][49][50] that the code obtained by the authors using Coarray Fortran is as efficient as the MPI code.At the same time, this code is much less complicated, which allows the development of new libraries [51] and language extensions [52].Note that the multidimensional decomposition of computations in the code has practically no effect on its efficiency [53,54].
The limitations of existing solutions in terms of reproducing individual relativistic gas flows, which are primarily associated with the numerical technique used, as well as the lack of codes that use all the advantages of one-way network communications implemented in MPI 3.0 and especially in Coarray Fortran, show us the need to develop a new numerical code for modeling relativistic gas flows.The implementation priority is Coarray Fortran for us, since its use significantly simplifies the code in terms of parallel implementation, which will simplify its implementation.Section 2 provides a description of a new numerical method for solving the equations of special relativistic hydrodynamics.Section 3 is devoted to a parallel implementation of the method.In Section 4, the numerical method is verified on standard tests.Section 5 presents two astrophysical applications.Section 6 provides conclusions to the paper.

Numerical Method
Let us first introduce the physical variables: ρ-density, v-velocity vector, and ppressure.For definiteness, we take the speed of light c ≡ 1.In this case the Lorentz factor Γ is defined by the formula The condition imposed on the speed of light ensures that the velocity modulus is limited by unity.This fact will be used in the construction of the numerical scheme.In the code, we will use an ideal gas model whose equation of state is written in the form of an equation for special enthalpy h: where γ is the adiabatic index.The speed of sound c s is defined by the formula To write the equations of special relativistic hydrodynamics, we introduce the following conservative variables: relativistic density D = Γρ, relativistic momentum M j = Γ 2 ρhv j , where v j are components of the velocity vector v at j = x, y, z, and total relativistic energy E = Γ 2 ρh − p.We write the system of equations in vector conservative form from classic book [55]: where δ jk is the Kronecker symbol.The conservative variables are easily calculated from the physical variables, but for the inverse transformation the following nonlinear equation has to be solved: For this, Newton's iterative method is used: where the derivative of the function f is When a given accuracy is reached, the iterative process is terminated.The value from the previous time-step is used as the initial pressure approximation.
The system of Equation ( 4) can be rewritten in vector form as follows: Then for any cell, Godunov's scheme can be written as x,i+1,k+ 1 2 ,l+ 1 2

−F *
x,i,k+ 1 2 ,l+ 1 2 where ∆ x,y,z are the grid spacings along the corresponding coordinates, τ is a time-step, and F * are the fluxes of the corresponding quantities through the cell boundaries obtained by solving the Riemann problem.Let us describe a method for constructing a numerical scheme for solving the Riemann problem.
To construct the numerical scheme, we consider the one-dimensional equations: With the operator separation scheme, the system of Equation ( 10) is divided into the Eulerian stage with the pressure forces at work: and the Lagrangian stage at which the advective transport takes place: For the systems O (Original)-( 10), E (Eulerian)- (11), and A (Lagrangian)-( 12), we will use a Rusanov-type scheme, which can be written for an arbitrary cell in a onedimensional statement in the following form: where ∆ is the grid spacing, τ is a time-step, and F O,E,A i are the solutions of Riemann problems for equations O, E, and A on the interface between cell i − 1 2 and cell i + 1 2 .We have where L denotes the left cell and R the right one.Let us solve the system at each stage of the method.
To solve system (11), we find the eigenvalues of the Jacobian.Omitting long but obvious manipulation, we obtain where , and v n is the normal velocity component.To implement Rusanov's scheme, at the Eulerian stage we use the maximum characteristic angle: The scheme for calculating the flux is Note that the characteristics are limited by the angle associated with the speed of light, that is, the condition λ ≤ 1 is satisfied.In this way we can considerably simplify the numerical scheme.
When solving the system of Equation ( 12), the eigenvalue of the Jacobian matrix is the numerical advection speed calculated from the equation where, taking into account Rusanov's scheme (14) and the form of Equation ( 12), we obtain the following simple scheme: Thus, we have obtained fluxes for the Eulerian (E) and Lagrangian (A) stages.The flux for solving the full system of equations of special relativistic hydrodynamics is the sum of these fluxes: Taking into account that the velocity is naturally limited by the speed of light, the flux in trivial form may be written as follows: To reduce the dissipation of the numerical method, we will use a piecewise parabolic reconstruction of the physical variables.For definiteness, a piecewise parabolic function of a physical variable q(x) will be constructed on a regular grid with spacing ∆ on the interval [x i , x i+1 ].For simplicity, we will use subscript i.In general form, the parabola may be written as follows: where q i is the value in the cell center, ξ = (x − x i )∆ −1 , q i = q L i − q R i , and q ; the conservation condition is Let us describe the construction of the parabola step by step.At the first step, we construct δq i = 1/2(q i+1 − q i−1 ).To avoid the extrema of the functions, we use the formula Then, we recalculate the values at the boundary using the fourth-order interpolant At the second step, we begin to construct the local parabola itself using the formula If the local parabola is nonmonotone (this may take place at discontinuities), we rearrange the values at the boundaries q L i , q R i by the formulas Thus, the boundary values satisfy the monotonicity conditions.At the third step, the parameters of the parabola are reconstructed taking into account the new values at the cell boundaries: Recall that the local parabolas are used as important parts of the Riemann problem, since there may be discontinuities at the interfaces.
At the fourth step, the parabola is finally reconstructed, taking into account the new values at the cell boundaries: The left and right values of the physical variables are recalculated using the formulas where λ L,R is the corresponding eigenvalue at each step.The Courant-Friedrichs-Lewy (CFL) condition is a necessary condition for the stability of an explicit Godunov-type scheme.From this condition, we determine the time-step.Completing the description of the numerical method, we again use the limitation of velocities by the light speed to limit the time-step.Then, the CFL condition can be formulated as follows: where CFL is the Courant-Friedrichs-Lewy number.

Parallel Code
When using the 256 3 computational grid, one time-step is calculated in approximately a minute; a complete calculation takes approximately a day.
Of course, it is necessary to develop a parallel implementation of the code based on the geometric decomposition of the computational domain with overlapping subdomains.Decomposition occurs on the images of the Coarray Fortran program.The description and study of the parallel implementation will be described below.In this section, we consider the structure of the parallel code.For simplicity, we will describe in detail a one-dimensional version of the code.However, all computational experiments and efficiency studies will be made using a three-dimensional version (see https://gitflic.ru/project/igorkulikov/rhd3dhllk, accessed on 29 May 2022).
The parallel implementation is based on a one-dimensional geometric decomposition of the computational domain with overlapping over two cells.This is due to the use of a five-point stencil for calculating the parabolas.The interprocess communications are performed by using Coarray Fortran based on one-way communications implemented in MPI-3.0 standard.Note that the interprocess communication tools are standards of the language, not of an external library.
This implies that the software of the code being developed is portable.A schematic diagram of the geometric decomposition of calculations and the one-way communications using, as an example, a one-dimensional implementation, are shown in Figure 1, where the blue color shows the cells for the boundary conditions, the green color denotes the calculation cells, and the gray color is the overlapping of subdomain cells.Note that in one-way interaction the exchange process is always initiated by the left process (image for Coarray Fortran programs).This is due to the fact that reading and writing data of a right image is always in the first cells.Therefore, the information about the local computational domain size contained in each image need not be stored in each program image.Nevertheless, this organization of data exchanges ensures data integrity when calculating the next time-step.Implementation details of the Coarray Fortran code are given in Appendix A. To study the parallel implementation, we use two nodes of an NKS-1P cluster, each with two 24-core Intel Xeon 6248R processors with Hyper-Threading support.Let us perform a preliminary study of the program speedup on a 256 3 grid with various numbers of nodes, sockets, cores, and threads, which, upon multiplication, forms the number of images.
A study of the calculation time shows that to increase the program images optimally, the following scheme should be used: increase the number of nodes, threads, processors, and cores, as shown in Table 1.To study the speedup of the software implementation, we use a 256 3 calculation grid with different numbers of Coarray Fortran program images, and to study the scalability, a 128 3 grid per one program image.
A 17-fold speedup was achieved with 32 images and 92% efficiency when using 192 Coarray Fortran images, results are shown in Figure 2.

Verification
For verification, we formulate three one-dimensional tests to consider the dynamics of a relativistic gas on the interval [0, 1] up to t = 0.4.At point x 0 = 0.5 from the initial time the gas has a discontinuity.The adiabatic index γ = 5/3.We use non-reflecting boundary conditions.In the simulation, 200 cells for all one-dimensional tests are used.In addition to comparing a numerical solution with an analytical one, we will analyze the order of convergence of the method in L 1 norm: where u(x i ) is the exact solution at point x i , u i is the numerical result, and h is the uniform grid spacing.
First, consider a one-dimensional shock wave test with the following parameters: pressure p L = 40/3, density ρ L = 10, and velocity v L = 0.The parameters of the gas to the right of the discontinuity are: pressure p R = 10 −8 , density ρ R = 1, and velocity v R = 0.The simulation results are shown in Figures 3 and 4   Figure 5 shows the behavior of the numerical solution error and the order of convergence of the method at grid refinement.The method correctly reproduces all components of the solution.The shock wave dissipates only over two cells, which is due to the use of piecewise parabolic methods.The shock wave peak in the density function is reproduced correctly.No additional dissipation is observed on the rarefaction wave.Note that there are no entropy traces in the pressure and velocity functions.Although the error curves decrease rather smoothly, for small grids there are jumps in the order of convergence curves.With significant grid refinement, the density and velocity functions converge to 0.8, whereas the pressure function converges to 0.9.
As a second test, consider a stronger shock wave with the following parameters: pressure p L = 1000, density ρ L = 1, and velocity v L = 0.The parameters of the gas to the right of the discontinuity are: pressure p R = 10 −2 , density ρ R = 1, and velocity v R = 0.The simulation results are shown in Figures 6 and 7   Figure 8 shows the behavior of the numerical solution error and the order of convergence of the method at grid refinement.The main difficulty of this test is that behind the shock wave front there is a very thin shell, which covers only two cells in the experiment.Even at this resolution, the numerical method was able to reproduce the shell location, and only 60% of the density peak was reproduced.Note that in the shell zone the solution dissipates over two cells in the direction of the shock wave and over two cells behind the shell front.At the same time, the mass of the dissipation zone corresponds well to the mass of the shell.The other parts of the solution are reproduced correctly.The difficulties in reproducing the shell impose some restrictions on the order of convergence of the density function.Note that on small grids the shell is reproduced with a small, increasing order of convergence to 0.6.At the same time, the convergence of the velocity and pressure functions is quite stable; they converge to the first-order for the pressure function since the smooth part of the rarefaction wave prevails in the solution, and to 0.8 for the velocity function due to a large shock wave amplitude.
Let us add complexity to the second test by adding a nonzero tangential initial velocity.The parameters of the gas to the left of the discontinuity are: pressure p L = 1000, density ρ L = 1, normal velocity component v L x = 0, and tangential velocity component v L y = 0.The parameters of the gas to the right of the discontinuity are: pressure p R = 10 −2 , density ρ R = 1, normal velocity component v R x = 0, and tangential velocity component v R y = 0.99.The simulation results are shown in Figures 9 and 10   Figure 11 shows the behavior of the numerical solution error and the order of convergence of the method at grid refinement.In contrast to ideal hydrodynamics, in the equations of special relativistic hydrodynamics, the tangential velocity greatly affects the normal velocity component due to the nonlinearity of the Lorentz factor.The test is much more complicated in that the tangential velocity is given in a cold gas.Therefore, the method must correctly reproduce the interaction of the velocity resulting from the shock impact of the hot gas.Note that the method reproduces well the shock wave zone with dissipation over two cells at the front and over four cells behind the shock wave shell zone, as in the previous test.Not only the location but also the density peak amplitude is reproduced correctly.The absence of shell-type solutions causes more-regular behavior of the error and the order of convergence of the method.Thus, for the pressure function, the order of convergence is one (the reason being the same as in the previous test), and for the density and velocity functions it is 0.8.
Let us take up the problem of discontinuity breakdown in a two-dimensional statement.For this, consider gas dynamics in the [−1/2; 1/2] 2 domain up to time t = 0.4.Let us divide the domain into four parts in geographical directions: NW-[−1/2; 0] × [0; Let us define the parameters of the gas in each part as follows: ρ, v x , v y , p NE = (0.1, 0, 0, 0.01), ρ, v x , v y , p NW = (0.1, 0.99, 0, 1), ρ, v x , v y , p SE = (0.1, 0, 0.99, 1), ρ, v x , v y , p SW = (0.5, 0, 0.99, 1), The adiabatic index γ = 5/3.Nonreflecting boundary conditions are set on all boundaries.In the simulation, a grid of 200 2 cells is used.The results of the simulation are presented in Figure 12.The evolution of the interaction of discontinuity decays in a two-dimensional formulation is shown.The first feature of the test is a numerical artifact due to the interaction of contact discontinuities.The second feature of the test is the interaction of shock waves and the formation of a density jump at the front.
The interaction between domains NW/SW and SE/SW form contact discontinuities with a jump in the tangential velocities.These contact discontinuities result in a numerical artifact that expands over time in domain SW.The formation of two interacting shock waves causes a density jump in domain NE, which has the same amplitude and limited dissipation.In this test, only the interaction of different waves is investigated, since the reproduction of individual waves has already been considered in the previous test problems.

Numerical Simulation
Let us consider the problem of shock wave formation for a gas at rest in a threedimensional statement.For this, consider gas dynamics in the [−1; 1] 3 domain up to time t = 0.4.Let the domain have a spherical radius r 0 = 0.5.In the internal part of the domain, the gas parameters are: pressure p 0 = 40/3 and density ρ 0 = 10.In the external part of the domain, the gas parameters are: pressure p = 10 −8 and density ρ = 1.The adiabatic index γ = 5/3.Nonreflecting boundary conditions are used.The results of the simulation are presented in Figure 13.A grid of 400 3 cells is used in the simulation.According to the results of this test, the solution is invariant with respect to rotation in the three-dimensional case.The density peak is reproduced correctly, and the shock wave is reproduced well (its analytical solution is shown by a dotted line).This problem is considered as a central explosion of a supernova without analyzing the cause of the explosion (thermonuclear reactions or a core collapse) in the model of special relativistic hydrodynamics.
As a second problem, consider interacting relativistic jets.We will simulate two interacting galactic jets with density ρ J = 10 −3 cm −3 and radius R J = 200 pc moving with Lorentz factor Γ = 10 and having relativistic Mach number M = 8.The atmosphere of the galaxy has temperature T A = 10 7 K and density ρ A = 10 −2 cm −3 .The adiabatic index γ = 5/3.Figure 14 presents the results of simulating the evolution of the galactic jets. .Schlieren density was chosen to demonstrate the interaction of discontinuities in the jet interaction problem.Until the moment of interaction, the structure of the shock wave front and the hotpoint, which is formed in the shell due to the interaction of waves, is visible.After the interaction of jets, an almost vertical shell structure is formed between two jets, while the shock wave fronts do not contain any numerical artifacts that other implementations may have in such interactions.
In both jets, a spherical shock wave propagates forward at the speed of light.Behind the shock front there is a shell with a system of interacting shock waves.Inside the shell there is a sparse cocoon area.During the interaction of the jets, the hot point zone becomes extended and covers about 1 kpc.

Discussion and Conclusions
Distinctive features of the code are as follows: 1.
To increase the order of accuracy on smooth solutions and decrease dissipation at discontinuities, a piecewise parabolic reconstruction is used for the physical variables since they are responsible for the thermodynamic state of the continuous medium [43].

2.
We use a piecewise-parabolic reconstruction of physical variables although the calculations in the scheme are based on the use of vectors of conservative variables and their flows through cell boundaries.This is due to the fact that they determine the thermodynamic state of a continuous medium.This was emphasized in [43] when studying the accuracy of the discontinuous Galerkin method on the Einfeldt problem.

3.
When using Rusanov-type methods as part of Godunov's method, no full spectral problem is solved; this is the case, for instance, when using an equation of state of the Taub-Mathews type [56].

4.
We plan to extend the model for a multiphase medium, where the equation of state is more complicated [57].In the present paper, we neglect the right-hand side of the energy equation, although this term can be easily included in our numerical scheme [58]. 5.
On the basis of the above computational experiments, we will try to avoid using the MPI technology where possible.The Coarray Fortran technology is a good alternative to develop parallel programs for distributed memory architectures.6.
In a number of problems, a significant part of the computational domain is occupied by flows with a high value of the Lorentz factor and low pressure values.In these cases, the energy equation is used in the form: The above statements are disputable.However, we believe they are appropriate when discussing the code.To conclude: in the present paper, a new code for describing flows in a model of special relativistic hydrodynamics and implemented with Coarray Fortran parallel programming technology was presented.During the development of the code, a new original numerical method was constructed with high-order accuracy on smooth solutions and low dissipation at discontinuities.The numerical method has been thoroughly verified.A code scalability of 92% is obtained on a cluster with Intel Xeon 6248R NKS-1P with 192 Coarray Fortran images.As an astrophysical application, we considered two interacting relativistic jets.

.
i r a n k = t h i s _ i m a g e ( ) 3 .
i s i z e = num_images ( )

.
Nl oc a l = N/ i s i z e 5 .
i _ s t a r t _ i n d e x = (N/ i s i z e ) * ( irank −1) 6 .
Nl oc a l = Nl oc a l + 1 8 .
i _ s t a r t _ i n d e x = i _ s t a r t _ i n d e x + i r a n k − 1 9 .
e l s e 1 0 .
i _ s t a r t _ i n d e x = i _ s t a r t _ i n d e x + mod(N, i s i z e ) 1 1 .e n d i f calculate the first approximation of the shift index i_start_index of the array part; 6.
to obtain a uniform distribution of cells over the images, in the first processes the local array size must not exceed the other ones by more than unity; 7.
add unity to the calculated local size for all first images; 8.
shift the shift index taking into account the fact that the size of the previous images is greater than the first approximation by unity; 9.
in the subsequent processes, the first approximation of the grid part size is not recalculated; 10. recalculate the shift index with the local sizes of the first images increased by unity; 11. this completes the calculation of the shift index and local size;

Figure 1 .
Figure 1.Geometric decomposition of the calculation domain and scheme of communications.

Figure 2 .
Figure 2. Speedup (left) and scalability (right) of the parallel program.A 17-fold speedup was achieved on a 256 3 grid using 32 images of the Coarray Fortran program, and 92% efficiency was achieved using 192 Coarray Fortran images when each image had a grid of 128 3 . .

Figure 4 .
Figure 4. One-dimensional shock wave problem.Exact solution (solid line), numerical solution (circles), and difference between exact and numerical solution (red circles) for physical variables: (a) density, (b) pressure, and (c) normal velocity component.

Figure 5 .
Figure 5. One-dimensional shock wave problem.Behavior of errors (a) and order of convergence (b). .

Figure 6 .
Figure 6.One-dimensional strong shock wave problem.Exact solution (solid line) and numerical solution (circles) for physical variables: (a) density, (b) pressure, and (c) normal velocity component.

Figure 7 .
Figure 7. One-dimensional strong shock wave problem.Exact solution (solid line), numerical solution (circles), and difference between exact and numerical solution (red circles) for physical variables: (a) density, (b) pressure, and (c) normal velocity component.

Figure 8 .
Figure 8. One-dimensional strong shock wave problem.Behavior of errors (a) and order of convergence (b). .

Figure 9 .
Figure 9. One-dimensional strong shock wave problem with tangential initial velocity.Exact solution (solid line) and numerical solution (circles) for physical variables: (a) density, (b) pressure, (c) normal velocity component, and tangential velocity component (d).

Figure 11 .
Figure 11.One-dimensional strong shock wave problem with tangential initial velocity.Behavior of errors (a) and order of convergence (b).

Figure 12 .
Figure12.Two-dimensional test on discontinuity breakdown.Density function distribution at t = 0.1 (a), t = 0.2 (b), t = 0.3 (c), and t = 0.4 (d).The evolution of the interaction of discontinuity decays in a two-dimensional formulation is shown.The first feature of the test is a numerical artifact due to the interaction of contact discontinuities.The second feature of the test is the interaction of shock waves and the formation of a density jump at the front.

Figure 13 .
Figure13.Density distribution in a three-dimensional shock wave at t = 0.4.This test is an extension of the first one on the disintegration of the discontinuity to the three-dimensional case and demonstrates, in addition to the correct reproduction of the shock wave front, the invariance of the solution concerning rotation and the absence of the Carbuncle effect.

Figure 14 .
Figure 14.Schlieren density in the equatorial plane in the evolution of a galactic jet at 1000 years (a) and 2000 years (b).Schlieren density was chosen to demonstrate the interaction of discontinuities in the jet interaction problem.Until the moment of interaction, the structure of the shock wave front and the hotpoint, which is formed in the shell due to the interaction of waves, is visible.After the interaction of jets, an almost vertical shell structure is formed between two jets, while the shock wave fronts do not contain any numerical artifacts that other implementations may have in such interactions.

1 2 .
a l l o c a t e ( a r r a y ( N lo c al + 4 ) [ * ] ) 1 3 .i f ( i r a n k < i s i z e ) then 1 4 .a r r a y ( 1 ) [ i r a n k +1] = a r r a y ( N l oc al +1) a r r a y ( 2 ) [ i r a n k +1] = a r r a y ( Nl oc a l +2) 1 5 .a r r a y ( N lo c al +3) = a r r a y ( 3 ) [ i r a n k +1] a r r a y ( Nl oc al +4) = a r r a y ( 4 ) [ i r a n k +1] 1 6 .e n d i f 1 7 .sync a l l 1 8 .x = ( i + i _ s t a r t _ i n d e x − 2 ) * h − h / 2 .d0 1 9 .a r r a y ( i ) = f ( x ) 2 0 .d e a l l o c a t e ( a r r a y ) Let us consider Listing A1 in detail:1.definition of dynamic array located in shared memory between program images; 2. definition of variable irank-Coarray Fortran program image number; 3. definition of variable isize-the number of images of Coarray Fortran program; 4. with known calculation grid size N, calculate the first approximation of the size of the grid Nlocal for the corresponding image of the program; 5.
12. allocate memory for array taking into account boundary conditions and/or overlapping areas; 13. perform one-way communications with initialization of exchange between the left images of the program; 14. write two right extreme values of the local calculation grid into the first two cellsoverlapping areas of the next program image; 15. write the first two cells of the local calculation grid from the next program image into two rightmost cells-overlapping areas; 16. this completes the exchanges between the overlapping areas; 17. synchronize data in all images of the Coarray Fortran program; 18. calculate the coordinates of the center of cell of size h, taking into account the shift index; 19.assign the function value to the corresponding array element;

Table 1 .
Calculation time for one time-step when using a 256 3 grid.