# Numerical Calculations of the GRP Scheme for Nonconservative Ideal Fluid Mechanics Equations Relying on Parallel Calculation of Multifluid Grids

Authors: Huijing Zhan; Mingze Wu
Journal: Advances in Multimedia
 (2021)

---

## Abstract

In order to solve the numerical method of nonconservative ideal hydrodynamics equations, the viscous perturbation technique for solving nonconservative hydrodynamics equations is improved and tested by solving the Riemann problem. The calculation of nonconservative ideal fluid mechanics is based on the GRP format. This article aims at the calculation method of nonconservative ideal fluid mechanics in the GRP format. Riemann and the corresponding periodic vortex are processed. The multifluid network processing method in the article is compared with the current method. The result can prove that this format can be used to solve the nonconservative ideal fluid dynamics equation of multiple values in the GRP format group, its computing power is strong, and the result of the solution is accurate.


---

## Body

## 1. Introduction

Relying on the multifluid grid parallel computing method can indicate the process of computing the GRP. The algorithm can be used to complete the traditional computing tasks that can be completed on a large computer or even a supercomputer on a PC or HPC composed of multiple GRPs. Using this, the method can effectively reduce the cost [1–3], and at the same time, it can provide a new path and method for the long-term development of CFD as a computationally intensive discipline. However, it can be seen that the numerical value obtained by using this type of scheme to solve the nonconservative equations of the GRP scheme is more unreasonable than that obtained by using the conservative GRP scheme to solve the nonconservative ideal hydrodynamics [4–6]. Up to now, there is no theory to solve the problem of convergence of the numerical solution [7, 8].Aiming at the problem of low efficiency in the numerical calculation of nonconservative ideal hydrodynamics equations, this paper discretizes the values by using parallel calculation algorithm; that is, it is transformed into the GRP format and unstructured grid storage method. The combination of interpolation and reconstruction is used to ensure that the accuracy of the space is improved to the second stage. The calculation of the above format is completed on several values in the GRP format of the equation system by using the multifluid grid parallel calculation algorithm. Finally, this article shows, through the example analysis, the accuracy and practicability of the calculation method of the text.


## 2. Parallel Computing Method

There is a multi-index calculation system composed ofn calculated objects u1,u2,…,un and m indicators x1,x2,…,xm; xij=xjxii=1,2,…,n;j=1,2,…,m is the calculated data matrix (decision matrix) of the observed value of the calculated object ui on the indicator xj, which can be expressed as [9](1)A=xijn×m=x11x12⋯x1mx21x22⋯x2m⋯⋯⋯⋯xn1xn2⋯xnm.Among them,m,n≥3. The data in A are the normalized data after preprocessing.Parallel computing method process is described in this article as general transformation.(2)yi=fxi1,xi2,…,xin,i∈N.Among them,f is the positive transformation function, yi is the comprehensive calculated value of the calculated object ui, and u1,u2,…,un is sorted according to the y1,y2,…,yn value from large to small to complete the comparison of the advantages and disadvantages of u1,u2,…,un.Hypothesis 1.
Any calculated object has the dual goal of “opening up the gap between competitors” and “developing its own special warfare” and comprehensively highlights its parallelism.
Hypothesis1 is a quantitative description of the idea of autonomous parallel computing.Definition 1.
Setαij and βij, respectively, as the amount of column parallelism and row parallelism of the calculated object uii∈N on the indicator xjj∈M, and they satisfy(3)αij=1n−1∑k≠ixij−xkj,i∈N,j∈M,k∈N,βij=1m−1∑p≠jxij−ip,i∈N,j∈M,p∈M.
Ifλij=μαij+ηβij,i∈N,j∈M, is used, then λij is called the autonomous parallel amount of the calculated object uii∈N with respect to the indicator xjj∈M, where μ is the competitive target coefficient, η is the development target coefficient, μ,η∈0,1, and μ+η=1.
The column parallel quantityαiji∈N,j∈M reflects the strength difference between the jth index of the calculated object ui and the overall n−1 counted objects, and the row parallel quantity βij reflects the overall jth index of the calculated object ui and other m−1 items.


### 2.1. Parallelization of Boundary Conditions and Residual Value Calculation

In order to make the parallel processing meet the GRP calculation of boundary conditions and residuals, it is necessary to enable the GRP algorithm to carry out parallel operation and obtain the corresponding data kernel thread on GRP according to the boundary conditions and residual operator kernel function [12]. The code written on GRP can be parallelized and calculated to realize the upper boundary conditions of GRP and the kernel corresponding to the residual value operator function. At the same time, the loop body in each function corresponds to the thread corresponding to the kernel on the GRP (Figure 1).Figure 1
The boundary condition calculation thread mapping process.


### 2.2. Parallel Solving of Linear Systems

Through the numerical analysis of the GRP scheme for nonconservative ideal hydrodynamics equations and the particle swarm optimization in the search space, the GRP scheme for nonconservative ideal hydrodynamics equations is composed of several numerical calculation components, which are parallelized by the multifluid grid parallel computing method. Therefore, the GRP format of nonconservative ideal fluid mechanics equations is used to define several numerical calculations using diversity and accuracy. Comprehensive evaluation is used to complete the nonconservative ideal fluid mechanics equations’ GRP format. Several numerical calculation options are passed in the nonconservative ideal fluid mechanics equations. The analysis results of different numerical calculations in the GRP format group are mainly obtained by multifluid grid parallel calculation by forming a group of particles in the search space of nonconservative ideal fluid mechanics equations. This uses a multifluid grid parallel calculation method to parallelize it, and the repeated calculation process is as follows [13]:(1)
Parallel solution (4) as the initial solution for postintegration.(4)DΔQin0=−Rin.(2)
Calculate several inner iterative steps to solveΔQink, and get(5)DΔQink+∑j∈NiLUΔQJNK−1=−Rin1≤k≤kmax.(3)
Set the final interpretation of thenth time step to be the solution of the final integration.(6)ΔQin=ΔQinkmax.The calculation results are used in some numerical calculations of the GRP format of the nonconservative ideal fluid mechanics equations, and some numerical calculations and analyses of the GRP format of the nonconservative ideal fluid mechanics equations are obtained [14].


### 2.3. Ways to Improve Parallel Efficiency

The small units of GRP multiline processes and GRP threads are operated. If the threads are in a thread wave, then the branch flows are operated in the serial form. As a result, operating efficiency is reduced. Because the number of each direction controller in the calculation space and the number of faces and the size of the thread speed are not equal, first, calculate whether it is within its range. It can be seen that the branch flow is more prominent in the structural network [15]. In this article, the unstructured grid method is used to store the multidimensional calculation units (control volume, surface, and node) of the structured grid into one dimension, and the unit needs to be expanded at the end to avoid the branch flow, as shown in Figure 2. The number of extended units cannot exceed the size of the maximum thread speed. In this way, the number of branch flows can be zero. For unstructured grids, only extension units are added.Figure 2
Thread mapping of structural grid unit calculation.


## 2.1. Parallelization of Boundary Conditions and Residual Value Calculation

In order to make the parallel processing meet the GRP calculation of boundary conditions and residuals, it is necessary to enable the GRP algorithm to carry out parallel operation and obtain the corresponding data kernel thread on GRP according to the boundary conditions and residual operator kernel function [12]. The code written on GRP can be parallelized and calculated to realize the upper boundary conditions of GRP and the kernel corresponding to the residual value operator function. At the same time, the loop body in each function corresponds to the thread corresponding to the kernel on the GRP (Figure 1).Figure 1
The boundary condition calculation thread mapping process.


## 2.2. Parallel Solving of Linear Systems

Through the numerical analysis of the GRP scheme for nonconservative ideal hydrodynamics equations and the particle swarm optimization in the search space, the GRP scheme for nonconservative ideal hydrodynamics equations is composed of several numerical calculation components, which are parallelized by the multifluid grid parallel computing method. Therefore, the GRP format of nonconservative ideal fluid mechanics equations is used to define several numerical calculations using diversity and accuracy. Comprehensive evaluation is used to complete the nonconservative ideal fluid mechanics equations’ GRP format. Several numerical calculation options are passed in the nonconservative ideal fluid mechanics equations. The analysis results of different numerical calculations in the GRP format group are mainly obtained by multifluid grid parallel calculation by forming a group of particles in the search space of nonconservative ideal fluid mechanics equations. This uses a multifluid grid parallel calculation method to parallelize it, and the repeated calculation process is as follows [13]:(1)
Parallel solution (4) as the initial solution for postintegration.(4)DΔQin0=−Rin.(2)
Calculate several inner iterative steps to solveΔQink, and get(5)DΔQink+∑j∈NiLUΔQJNK−1=−Rin1≤k≤kmax.(3)
Set the final interpretation of thenth time step to be the solution of the final integration.(6)ΔQin=ΔQinkmax.The calculation results are used in some numerical calculations of the GRP format of the nonconservative ideal fluid mechanics equations, and some numerical calculations and analyses of the GRP format of the nonconservative ideal fluid mechanics equations are obtained [14].


## 2.3. Ways to Improve Parallel Efficiency

The small units of GRP multiline processes and GRP threads are operated. If the threads are in a thread wave, then the branch flows are operated in the serial form. As a result, operating efficiency is reduced. Because the number of each direction controller in the calculation space and the number of faces and the size of the thread speed are not equal, first, calculate whether it is within its range. It can be seen that the branch flow is more prominent in the structural network [15]. In this article, the unstructured grid method is used to store the multidimensional calculation units (control volume, surface, and node) of the structured grid into one dimension, and the unit needs to be expanded at the end to avoid the branch flow, as shown in Figure 2. The number of extended units cannot exceed the size of the maximum thread speed. In this way, the number of branch flows can be zero. For unstructured grids, only extension units are added.Figure 2
Thread mapping of structural grid unit calculation.


## 3. Numerical Method of the GRP Scheme for Solving Nonconservative Hyperbolic Equations

This paper takes a one-dimensional problem as an example to propose two types of new numerical methods based on nonconservative ideal fluid mechanics for solving nonconservative equations of GRP schemes and introduces the existing GRP schemes for solving equations.Consider the uniform spatial grid division with the step length ofΔx, the time step length is denoted as Δx, and the approximate value tn of the unknown function at the center point of the grid Ii=xi−1/2,xi+1/2 at the time of known xi is assumed to be Wi, and the unknown function needs to be calculated at the time of tn+1 approximate value.Method 1.
LetW=ρ,ρu,ρeT; the equation group of the GRP format can be equivalently written as(7)∂W∂t+AW∂W∂x=0.
Here,(8)AW=010−u22uγ−1−γeuγeu.
The characteristic value ofA(W) is(9)λ1=u−c,λ2=u,λ3=u+c.
Among them, the speed of soundc=γp/ρ; the corresponding matrix R(W) composed of eigenvectors in the right column and the matrix L(W) composed of eigenvectors in the left row are, respectively,(10)RW=111u−cuu+cγe0γe,LW=u2c−12c12γe10−1γe−u2c12c12γe.
The main steps ofMethod 1 are as follows:(1)
CalculateRi+1/2=RWi+1/2 and Li+1/2=LWi+1/2 by formula (12); among them, Wi+1/2=Wi+Wi+1/2(2)
Letv¯j=Li+1/2Wf,j=i−2,…,i+3, regard v¯j as the grid average of a certain function v, use the 5th-order WENO reconstruction technique to reconstruct the approximate value v^i+1/2± of the function on the ith grid boundary, and then inversely transform it to the physical space; that is, W^i+1/2±=Ri+1/2v^i+1/2±.(3)
CalculateAi=AWi by formula (10), and let Ai±=Ai±αiI, where I is a 3 × 3 identity matrix and αi=max1≤l≤3,i−2≤j≤i+3λji(4)
Letbi±=Ai±W^i+1/2±−W^i−1/2± form a semidiscrete format dWit/dt=−1/Δxbi++bi−(5)
Use the third-order TVD Runge–Kutta method to solve the semidiscrete problem and advance one time step
Traditional algorithm:(11)W=ρ,ρu,ρeT,Q=0,0,uT,F˜W=ρu,ρu2+p,ρeuT.
Then, the GRP format of the equation group can be expressed in a form similar to (7):(12)∂W∂t+∂F˜W∂x=p∂Q∂x.
The Jacobian matrix of fluxF˜W is(13)A˜W=∂F˜W∂W=110−u22uγ−1−eueu.
The characteristic value ofA˜W is(14)λ˜1=u−c˜,λ˜2=u,λ˜3=u+c.
Among them,c˜=p/ρ; the corresponding matrix R˜W composed of eigenvectors in the right column and the matrix L˜W composed of eigenvectors in the left row are, respectively,(15)R˜W=111u−c˜uu+c˜e0e,L˜W=u2c˜−12c˜12e10−1e−u2c˜12c˜12e.
The main steps of the traditional algorithm are as follows:(1)
Calculated by formula (18), Ri+1/2=R˜Wi+1/2 and Li+1/2=L˜Wi+1/2; among them, Wi+1/2=Wi+Wi+1/2.(2)
Local feature decomposition, that is, letVj=Li+1/2Wj and Gj=Li+1/2F˜Wj,j=i−2,…,i+3.(3)
The local Lax–Friedrichs flow function is split; that is, letGj±=Gj±αiVj,j=i−1,…,i+3. Here, αi=max1≤l≤3,i−2≤j≤i+3λ˜jl.(4)
Makev¯j=Gj+,j=i−2,…,i+3, regard v¯j as the grid average value of a certain function v, use the 5th-order WENO reconstruction technique to reconstruct the approximate value v^i+1/2− of the function on the ith grid boundary, and make G^i+1/2+=v^i+1/2−, and then inversely transform it into the physical space. That is, F^i+1/2+=Ri+1/2G^i+1/2+. Similar use of Gj−j=i−2,…,i+3 representations can be regarded as F^i+1/2− representations.(5)
Due toujj=i−2,…,i+3, reconstruct the approximate value on the ith grid boundary using the 5th-order WENO reconstruction technique u^i+1/2−; then, Q^i+1/2=0,0,u^i+1/2−T.(6)
MakeF^i+1/2=F^i+1/2++F^i+1/2−, and the semidiscrete format(16)dWitdt=−1ΔxF^i+1/2−F^i−1/2+piQ^i+1/2−Q^i−1/2.(7)
Use the third-order TVD Runge–Kutta method to solve the semidiscrete problem and advance one time step.
Direct discretion of nonconservative systems will lead to incorrect shock wave velocity and false jumps through shock wave propagation. For this reason, this paper uses viscous perturbation technology to add a correction term to the right end of the GRP format (9) and (14) of the equation system. Then, the above two methods are used to solve; that is, Method 1 is used to solve the GRP format of the equation group.(17)∂W∂t=AW∂W∂x=Δt2D.
The GRP format of equations is solved with traditional algorithms.(18)∂W∂t+∂F˜W∂x=p∂Q∂x+Δt2D.
Among them,(19)D=00ρuxuxa1+γ−1uxρea2/c2Δt/Δx−ρutut.
Here,a1=u+c+u−c, a2=u+c−u−c, and time derivative term ut=−uux−γ−1ρex/ρ. The spatial derivative term is approximated by the central difference quotient. D used in this paper is similar to that in [6], which is derived from the first-order upwind style, but from the numerical experiments in the next section, it can be seen that it can also be used for the above two types of nonconservative ideals.
In order to compare with the above two new numerical methods mentioned in this article, this article also introduces the GRP format mentioned in [5].
GRP format: the definition matrix(20)signAW=RWsignΛLW.
The diagonal elements of the diagonal matrixsignΛ are(21)signλl=1,λl>0,0,λl=0,l=1,2,3,−1,λl<0.
The main steps to solve the GRP scheme of the nonconservative system of equations are as follows:(1)
Calculated by formula (12), RWi+1/2 and LWi+1/2. Calculated by formula (11), AWi+1/2. The eigenvalues of λ1, λ2, and λ3 are calculated from equations (20) and (21). signAi+1/2=signAWi+1/2, which protects Wi+1/2=Wi+Wi+1/2.(2)
CalculateF^i+1/2=F˜i+1−Fi˜/Δx+pi+pi+1/2Qi+1−Q/Δx. The expressions of F˜ and Q here are shown in equation (15).(3)
ByWin+1=Wi−Δt/2I−signAi+1/2F^i+1/2+I+signAi−1/2F^i−1/2, advance a time step to get the approximate value of tn+1 to time Win+1.


## 4. Test and Result Analysis

In this section, through numerical tests, the numerical methods of the above three nonpreserved equations in the GRP form and the nonconservative ideal fluid mechanics in the GRP form of solution equation (4) are tested and compared.Example 1.
(Sod’s Riemann problem). The initial conditions are(22)ρ,u,p=1,0,1,forx<0.5,0.125,0,0.1,forx>0.5.
When both sides of the equation meet the boundary conditions of the emergency tributary, in the examples cited in this article, they are all set to beγ = 1.5. First, divide the interval [0, 1] into 400 grids, and use nonideal fluid dynamics to solve the nonconservative ideal fluid dynamics equations in turn. Use the algorithm cited in this article to parallelize it, and then divide it into 1200 grids at equal intervals of [10, 1] = division, mainly in the way of multifluid grid parallel computing to solve the GRP format operation time. It is basically consistent with the method used and then carried out according to the 400 grids to obtain the calculation time. The calculation time needs t = 0.16. The numerical calculation results that can be obtained by this calculation can be referred to in Figure 3. The corresponding shock wave in the figure needs to use a proportion of about 3. This problem has a certain shock wave strength. It can be seen from Figure 3 that reasonable numerical results can be obtained using the algorithm in this paper. The shock wave density value obtained in the GRP format is slightly larger than the true value.Figure 3
Sod’s Riemann problem. (a)N = 400, the numerical results of the nonconservative ideal fluid dynamics equation in GRP form (4) are solved by nonconservative ideal fluid dynamics. (b) N = 400, Method 1 is the numerical result of solving nonconservative ideal hydrodynamic equations in the GRP form. (c) N = 400, the traditional algorithm solves the numerical results of the nonconservative ideal fluid mechanics equations in the GRP form. (d) N = 1200, numerical results of nonconservative ideal hydrodynamic equations in the GRP form.
(a)(b)(c)(d)Example 2.
(Lax’s Riemann problem). The initial conditions are(23)ρ,u,p=0.445,0.698,3.528,forx<0.5,0.5,0,0.571,forx>0.5.
Both sides of the equation divide the [0, 1] interval into 400 grids under emergency tributary boundary conditions, and nonideal fluid mechanics is used in turn. It can be solved by using the parallel calculation method of multifluid grids. Then, in 1200, solve (8) in the GRP format on the subgrid. The shock wave left and right pressure ratio is about 4.5; this problem has a strong shock wave intensity. The optimal numerical results are obtained by parallel calculation using nonconservative ideal hydrodynamics. The multifluid grid parallel calculation method solves the nonconservative ideal fluid mechanics equations in the GRP form, and the numerical results are relatively close. Compared with the GRP format of the nonconservative ideal fluid dynamics solution to the nonconservative ideal fluid dynamics equation, the calculation of the shock wave is slightly inferior, and the numerical result obtained by the GRP format is viscous, and the density value on the left side of the shock wave has a larger error.
Under adiabatic conditions, the preserved form of the quadratic compressibility ideal fluid mechanics equation group in the GRP form is(24)∂ρ∂t+∂ρu∂x+∂ρv∂y=0,∂ρu∂t+∂ρu2∂x+∂ρuv∂y+∂p∂x=0,∂ρu∂t+∂ρur∂x+∂ρv2∂y+∂p∂y=0,∂ρΕ∂t+∂ρE+pu∂x+∂ρE+pv∂y=0.
Among them,E=e+u2+v2/2, e=p/γ−1ρ. Here, u and v represent the velocity components of the fluid along the work axis and F-axis, respectively, and the other quantities have the same meaning as above. GRP format (24) of the system of equations is written in the nonconservative form as(25)∂ρ∂t+∂ρu∂x+∂ρv∂y=0,∂ρu∂t+∂ρu2∂x+∂ρuv∂y+∂p∂x=0,∂ρu∂t+∂ρuv∂x+∂ρv2∂y+∂p∂y=0,∂ρe∂t+∂ρeu∂x+∂ρev∂y+p∂u∂x+∂v∂y=0.
According to the above two examples of one-dimensional operations, the optimal solution cannot be obtained by solving the nonconservative ideal hydrodynamic equation. Conserved ideal fluid mechanics equation GRP format (24) uses traditional algorithms to solve nonpreserved equation GRP format (25) and displays the numerical results.Example 3.
The initial conditions for the two-dimensional periodic eddy current problem are(26)Tx,y=1−ε2γ−18π2γe1−r2,Sx,y=1,ux,y=1−εy−52πe1−r2/2,vx,y=1−εx−52πe1−r2/2.
Among them,x,y∈Ω:=x,y|0≤x,y≤10. The corresponding directions of x and y all meet the periodic boundary conditions, r2=x−52+y−52, and ε represents the strength of the vortex. T represents the absolute temperature, S represents the thermodynamic entropy, and they satisfy T = P/ρ and S=P/ρr. This paper assumes ε = 5 at 81 × 81 which is solved using a multifluid grid parallel computing algorithm. When t = 10 and x = 5, the density will slightly change along the Y-axis direction, as shown in Figure 4(a). Refine the grid to 161 × 161, and use the traditional algorithm to solve (25). The overall error of density along the Y-axis is shown in Figure 4(b). It can be seen from the figure that the traditional algorithm has the same high accuracy as the nonconservative ideal fluid mechanics at the extreme point, but the viscosity correction term added in this paper causes a certain error. Figure 4(b) shows that the error will increase with the density of the grid and reduce.Figure 4
Two-dimensional periodic vortex problem. (a) 812 uniform grid, x = 5, and the solid line represents the true solution of the density. (b) The error of the density obtained by the traditional algorithm on the consistent grid of 812 and 1612, and x = 5.
(a)(b)Example 4.
Regarding the two-dimensional Riemann problem (refer to [7]), the initial conditions are(27)ρuvpT=1.5,0,1.5,0.8<x<1,0.8<y<1,0.5323,1.206,0,0.3,0<x<0.8,0.8<y<1,0.138,1.206,1.206,0.029,0<x<0.8,0.8<y<0.8,0.5323,0,1.206,0.3,0.8<x<1,0<y<0.8.
Here,x,y∈Ω:=x,y|0≤x≤1,0≤y≤1x, y. The direction adopts tight tributary boundary conditions, which is the problem of a more complicated flow structure caused by shock wave interaction. Solve by nonideal fluid dynamics (24) on a uniform spatial grid of 200 × 200, use the traditional algorithm to solve (25), and calculate the density distribution at time t = 0.8 with reference to Figure 5. It can be seen from the figure that the two results are basically the same.Figure 5
Two-dimensional Riemann problem. (a) Use WEN05 to solve the density contours of GRP format (24) of the conservative equations. (b) The traditional algorithm solves the contour value with density of 20.
(a)(b)


## 5. Conclusions

This paper proposes a numerical method for solving the GRP scheme of the nonconservative ideal fluid mechanics equations. This numerical calculation method can show that several numerical calculation methods of the GRP scheme of the nonconservative ideal fluid mechanics equations can effectively calculate the values. The method is stable in operation. When the amount of calculation is required, it has no direct relationship with the calculation method used. Finally, an example analysis shows that the method has high calculation efficiency, good accuracy, and certain use value.



---

*Source: 5118005-2021-12-09.xml*