The numerical modeling of the collapse of molecular cloud on adaptive nested mesh

In paper, a numerical modeling of the collapse of molecular cloud is presented. The numerical method for solving the equations of hydrodynamics is based on the extension of HLL method for using of nested adaptive mesh technologies. In case of nested meshes, if neighboring cells have the same size as the one considered cell, the Rieman problem is trivial. If the neighboring cell is larger than a considered cell, then a uniform distribution of the parameters over the large cell is assumed. Then a new subcell is allocated in a larger cell corresponding to a size of the cell considered and being adjacent to it, after which the necessary data for a computing. If a neighboring cell is less than a considered cell, then values around a small cell are averaged. The energy behavior of the rotating molecular cloud is in quantitative agreement with the results of other codes before the moment of collapse, and after the moment of collapse. It should be noted that collapse occurs at exactly the same moment in time.


Introduction
The processes of collapse of molecular clouds and stars are being actively studied today in connection with the appearance of a large number of observational data. The phenomenon of collapse takes place both at the initial stage of the evolution of the star, and at the final stage of a supernova explosion with a collapsing core. Mathematical modeling plays a more important role in the theoretical study of such processes. One of the main problems in such problems is the ratio of scales: when the size of the star is 10 9 meters the collapsed core can reach a size of 10 4 meter. Even with the use of powerful supercomputers, it is impossible to adequately model such processes.
One of approaches to solving the collapse problems is adaptive grids, which allow the local grinding of the grid. Here is an overview of the most known codes based on this approach.
ART [1]. The code is based on the technology of adaptive grid trees. Thus, the classical adaptive mesh (AMR) is represented in the form of a regular (octet) tree structure, which naturally lies on the classical data structures, and therefore all classical tree parallel algorithms are applicable to them. To solve the Poisson equation, the relaxation method is used when the Poisson equation is reformulated as: to solve the last equation is used the method of successive upper relaxation. The MUSCL scheme is used to solve gas-dynamic equations.  [2]. The program complex is based on the MUSCL and PPM (piecewise-parabolic method) methods using the AMR approach. To solve the Poisson equation, a multigrid method is used, and at each stage of the cycle the Gauss-Seidel method is used. An important feature of the code from the point of view of computational methods is the use of an approach based on the construction of a three-dimensional piecewise parabolic function [3], which is important for the invariance of the numerical method. In addition, more "smooth" symmetric constraints were used [4,5].
ENZO [6]. At the heart of the software complex is the solution of the equations of magnetic gas dynamics with allowance for cosmological expansion. To simulate a collisionless component, the N-body model is used. The code includes a large number of subgrid processes: primordial chemical kinetics, cooling / heating functions, radiation transfer, as well as star formation processes and supernova explosion effects. Several solvors are used to solve hydrodynamic equations: PPM (implemented only for gas dynamics equations), MUSCL and finite-difference method. To solve the Poisson equation, an algorithm based on a fast Fourier transform is used. The so-called structured adaptive grid is used in the complex, the basic idea is to minimize the difference between the computational grid between neighboring cells. Such a structure allows the use of regular trees, where the subdomain is divided no more than twice, which increases the efficiency of using such computed grids.
GAMER [7]. The code implements the solution of the equation of gas dynamics using the AMR approach on graphic accelerators. To solve the equations of gas dynamics, the TVD approach is used, the combination of the method based on the fast Fourier transform and the method of successive upper relaxation is used to solve the Poisson equation. Perhaps the main feature of this complex is the implementation of the AMR approach on graphics cards. So the regular grid structure is naturally projected onto the GPU architecture, while the tree structure requires special approaches. This approach consists in using octets to define a grid, which is projected onto a separate stream of the graphics card. The main problem here is the formation of dummy cells for an octet, which is about 63 % of the time, but this procedure can be implemented independently for each octet.
RAMSES [8]. The code implements a numerical solution of the equations of gravitational gas dynamics using the AMR approach, based on splitting into octets. To solve the Poisson equation, a combination of the method based on the fast Fourier transform and the Gauss-Seidel method is used. To approximate the Poisson equation, a simple five-point finite-difference approximation is used, which was replaced by a more efficient 19-point approximation, which was implemented as an extension of the RAMSES code to the case of nonclassical gravity (MOND) [9].
However, adaptive meshes are based on grinding the grid no more than twice in neighboring cells, which is due to the peculiarities of software implementations and numerical methods. This approach actually does not allow reproducing multiscale processes.
A powerful alternative to adaptive grids is the nested grid approach, when in the cell a complete task with a sufficiently detailed grid is actually solved. In fact, the use of nested grids allows to reproduce any currents of the multi-scale type. The development of such an approach will allow, for example, to eliminate the gap problem between modeling the evolution of galaxies and the process of star formation, and also to reproduce in sufficient detail the turbulent flows and processes of nuclear collapse. The original approach to creating a computational model of the hydrodynamics of astrophysical objects based on the technology of nested grids will be presented in the article. In this article, we are not concerned with the issues of parallel implementation, although it is clear that the use of nested grids requires the balancing of processor load.
In the second section, a numerical scheme for solving the equations of gas dynamics and the Poisson equation will be described in detail, and an algorithm for rebuilding the embedded grid based on the mass criterion will also be presented. The third section will be devoted to modeling the process of cloud collapse. The conclusion is given in the fourth section.

The mathematical model
In a few papers in the field of computational fluid dynamics, instead of a non-conservative equation for internal energy, an equation for entropy is used, written in a conservative form [10,11,12,13]. The undoubted advantage of this approach is the guarantee of non-decreasing entropy and the effect of negative pressure, as well as the possibility of correctly describing supersonic flows. There are also some limitations on the use of this approach (see discussion at work [14]).

∂ ∂t
To solve equations (1) and (2)  the following subsections, a more detailed description of the numerical schemes will be given.

Solution of the equations of gravitational hydrodynamics
To solve the hydrodynamic equations, a modification of the original numerical method based on a combination of the method of operator separation, the Godunov method and the HLL scheme was used. This method combines all the advantages of the above methods and has a high degree of parallelism. The details of the numerical scheme are described in [15]. The main idea of the method is to use the vector notation of the equations of hydrodynamics: where v -is a vector of conservative variables. For equation (3) the numerical scheme in one direction is written in the form: where F -solution of the Riemann problem. We omit the derivation of the numerical scheme on the basis of the use of conjugate equations and the method of separating operators and formulate the final form of the solution of the Riemann problem: For discretization using nested meshes, we introduce in the three-dimensional domain of the solution a uniform cubic root grid with the coordinates of the cell centers where h -is the step of the root grid, I max , K max , L max -is the number of cells in the directions x, y, z. In the implementation, for convenience of organizing calculations and without losing the generality of the code I max = K max = L max = N was used. In the cell with the number (i, k, l) we introduce an embedded cubic grid with the coordinates of the cell centers   3). In the case of equal cell sizes (figure 3 (in the middle)) the solution of the Riemann problem is analogous to the solution of Riemann problems on internal interfaces of a nested grid and is trivial. In the case where the cell of an adjacent nested grid is larger than the one considered ( figure 3 (left)), then a uniform distribution of the quantities over the blue cell is assumed, and even for the interface between the reduced blue cell and the pink cell, the Riemann problem is solved. In the case where the considered pink cell borders on several cells of a neighboring embedded grid (figure 3 (right)), then using the assumption of a uniform distribution of hydrodynamic quantities over the pink cell, Riemann problems are solved on all interfaces, averaged. The detailed arrangement of quantities for the organization of calculations is shown in the figure 2. Also, we do not dwell on the solution of the Poisson equation in the first stage (a detailed description of the solution method can be found in [14]).

Solution of the Poisson equation
The Successive Over-Relaxation (SOR) method is an iterative process of finding the potential on an embedded grid for given initial and boundary conditions obtained from solving the Poisson equation on the root grid. The iterative process is organized as follows.
(i) On the internal nodes of the nested grid a numerical scheme for solving the Poisson equation is written: (ii) At each iteration a discrepancy is defined: where ω -relaxation parameter. (iv) When a small residual is reached or when the specified number of iterations is exceeded A similar approach to solving the Poisson equation has proved itself in a number of programming codes, for example, in the code GAMER [7].

Selecting the time step
Step by time τ is calculated from the Courant condition where CF L -Courant number, c = √ γp/ρ -sound speed, h min -the size of the minimal cell of the nested grid.

Regularization of the solution
At the final stage of solving the hydrodynamic equations, the solution is corrected. In the case of a gas-vacuum boundary using the formula: (6) in the rest of the region there is an adjustment, which guarantees nondecreasing entropy: Such a modification provides a detailed balance of energies and guarantees nondecreasing entropy.

Algorithm for rebuilding the grid
Restructuring of the grid takes place according to the criterion of the cell mass of the root grid. The size of the nested grid is determined from the condition: where C 1,2 -scaling constants, which are selected depending on the task (the characteristic density) and the requirements for a minimum resolution. So, for the dimensionless density ρ in the problems of cloud collapse, the equation (8) is rewritten as: that is, the unit density value corresponds to a nested grid of size , increasing the same density by an order of magnitude leads to a rearrangement of the grid with [8 × 8 × 8] on [32 × 32 × 32]. Verification of the necessity of grid reorganization is carried out at each step in time, therefore the grid changes no more than twice. The grid reconstruction scheme is shown in the figure 4. Figure 4 shows the schemes of projections of conservative quantities (density, angular momentum, entropy density and total energy). After the reorganization of the grid, the nonconservative (primitive in the case of relativistic hydrodynamics) are calculated from conservative variables.

Simulation of the collapse of the cloud
We will simulate a gas cloud bounded by a sphere of radius R 0 = 100 parsecs, with mass M g = 10 7 M ⊙ , with density distribution ρ (r) ≃ 1/r and temperature T ≈ 2000 . The adiabatic exponent corresponds to hydrogen γ = 5/3. The speed of sound is c ≈ 3.8 km/s. As the size values, we choose the following values: L 0 = 100 parsec, ρ 0 = 1.2 · 10 −18 / 3 , v 0 = 21 km/s. Then, in dimensionless quantities, the problem is posed as follows: ρ = 1.0 -the density of the gas cloud at the center, p = 2 × 10 −2 -is the pressure in the gas cloud at the center, γ = 5/3 -is the adiabatic exponent, [0; 3.2] 3 -is the computational domain. The behavior of various types of energies using the original approach described in this article, and using AstroPhi [16] code is given in the figure 5. We note that the behavior of the energies is qualitatively close enough for both realizations. However, we note that the simulation shown in Figure 5 with the AstroPhi code was about 20 hours using 32 cores of the resources of the Siberian Supercomputer Center (SCCC), while code based on nested grids allowed to obtain the result in two hours using a single core resources of SSCC. Nested adaptive grids are a powerful tool for modeling hydrodynamic flows that allow you to adapt the computational process to the features of the solution. At the same time, the organization of calculations is quite complex, although it is obvious that the use of the power of many classical processors, graphics accelerators and accelerators Intel Xeon Phi will further improve code performance based on nested grids. In the feature papers, the questions of parallel implementation will be described in detail.

Discussion
In a description of the paper are not three important topics which, although not the subject of this paper has been affected, still require a description.
(i) Using adaptive nested mesh implies the possibility of multiscale modeling. Such an approach can, for example, be used in the problems of star formation in evolution of galaxies. In this case, the evolution or collision of galaxies can be modeled on the root mesh, and the cells of root mesh with a characteristic size and time, corresponding to the evolution of the interstellar medium (ISM), can be used to directly model of the star formation as an individual problem. A similar approach can be used to simulate ISM turbulence with the ability to reproduce protostellar and protoplanetary systems. (ii) Now in the code, multi-core architectures (for example, for Intel Xeon Phi) was developed, which is develop by parallel implementation of iterations of main computational cycles. The main interest is the development for the massively parallel architectures. For such an implementation would use an approach based on the splitting of two groups of processes: processes, that perform calculations on the root mesh, and processes, that dynamically appointed for the nested mesh. The method for a definition such groups requires a priori information about the problem and is adjusted in the course of computing. The distribution of such computations is planned to be described in next papers. (iii) In the paper, the first order accuracy method was described. Undoubtedly, the use of highorder methods will provide an even better numerical solution. It's also obvious, that the piecewise-parabolic method on local stencil (PPML) requires a modification for adjacent cells of different sizes to create local parabolas. We need to change two fundamental points in method: to construct a parabola from the average values on the boundary cells (now in PPML method uses a non-compact template), thus avoiding problems with different size cells, and integrate the parabola within the cell, defining the value, arriving at the boundary between cells (similar to the work [13]), which also requires information only of local cell.
The description of such a method is planned to be described in next papers.

Conclusion
The original approach to the solution of the equations of gravitational hydrodynamics based on the technology of nested grids is presented in the article. The organization of a numerical solution of the equations of gas dynamics and Poisson equations is described in detail. The algorithm of grinding and enlargement of nested meshes on the basis of the mass criterion is given. The process of cloud collapse is modeled. The behavior of various types of energies is shown in the collapse of the cloud, the simulation results are compared using nested grid technology and using the regular detailed grid implemented in AstroPhi code.