Algorithm to Generate Liutex Core Lines Based on Forward Liutex Magnitude Gradient Lines

Abstract

Vortex definition and identification are extremely important for the study of fluid dynamics research. Liutex is a newly proposed concept that correctly represents vortex. Liutex is a vector whose direction is the local rotation axis and whose magnitude is twice the angular speed. To identify the unique structure of a vortex, a method known as the Liutex Core Line method has been developed, which displays the rotational core axis of a vortex. However, the original method is a manual method, which is not practical for real application, and an automatic algorithm is required for practical usage. Xu et al. proposed an algorithm by selecting the best line from a group of candidate lines, which is an important progress. In this chapter, from another perspective to solve this problem, a new algorithm is introduced based on forward Liutex magnitude gradient lines. Since gradient lines have the feature that they advance to the local maximums, the route will still result in a unique line, which avoids the process to find the best line. This algorithm has achieved some success for the Lambda vortex in early boundary layer transition.

1. Introduction

Defining and visualizing vortices have been a challenge in fluid mechanics for several decades, even though the vortex is the essential part of formulating the turbulent flow that consists of countless vortices in different sizes and strengths. Turbulence generation and sustenance are still a mystery after intensive research efforts of centuries. Kaczorowski et al. [1] and Xi et al. [2] performed deep theory and experiment research on this topic, especially the Rayleigh–Bénard (RB) convection. In 2019, Zhou et al. elaborated the hydrodynamic instabilities that induced turbulent mixing in wide areas including inertial confinement fusion, supernovae, and their transition criteria. Zhou [3, 4] described in detail Rayleigh–Taylor and Richtmyer–Meshkov instability and some related models. Zhou’s analysis is systematical, comprehensive, and sophisticated for flow instability and vortex generation, covering long history and state-of-the-art advances in turbulence research, which has clearly shown guidance for further and deeper scientific research. Since vortex is the sinus and muscles of turbulence, it is of great importance to find out the definition of vortex. Human’s understanding of vortex has gone through three stages [5].

The first generation is vorticity-based methods which have several weaknesses. One severe weakness is the contamination by shear. This shear contamination can be easily seen in the near-wall laminar boundary layer where there are high levels of shear with small amounts of rotational motion. High vorticity magnitude is seen in areas of high shear and low rotation, while areas with stronger rotation can be found to have relatively lower vorticity magnitude values [6]. Though the misconception that a vortex is defined by vorticity is still being published in textbooks and research papers, we know this is not appropriate. Then, the second generation that uses eigenvalue-based vortex identification criteria was developed. The second generation of vortex identification methods mainly used the eigenvalues of the velocity gradient tensor. Among the second generation, there are some popular methods such as Q criterion [7], $∆$ criterion [8], ${\lambda }_{\mathit{ci}}$ criterion [9], ${\lambda }_2$ criterion [10], and some other methods. Although these methods made advances over the first generation, they were still contaminated by shear in different degrees. Another problem is that they are all scalar-based. These methods require a user-specified threshold to draw an iso-surface. Since the threshold was somewhat arbitrary, there must be many different thresholds used to represent the vortical structure and different thresholds may give totally different vortical structures. These scalar-based methods also prevent us from determining the rotation axis except that in the paper of ${\lambda }_{\mathit{ci}}$ [9], Zhou et al. indicated the rotation axis direction aligns the eigenvector but ${\lambda }_{\mathit{ci}}$ was still defined as a scalar and it is contaminated by shear as well due to the fact that the non-orthogonal transformation is used in the definition. To find the proper definition of vortex, Liu et al. proposed Liutex [11], which is a vector, in 2018. The direction of Liutex is the local rotation axis, and the magnitude of Liutex represents twice the angular speed.

Since the invention of Liutex, many researchers have used and tested it and reported that it is a strictly mathematical vortex definition and the currently best vortex identification method. Cuissa et al. [12] found that Liutex matches the analytical result of Lamb-Osceen vortex. Shen et al. [13] commented that Liutex is the only method that fully identifies small-scale vortices in vertical slit fishways simulation. Xu et al. found Liutex similarity, i.e., both the frequency and wavenumber spectrum of Liutex match the -5/3 law [14, 15, 16]. More reports on Liutex application can be found from [17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29] and four published books [30, 31, 32, 33]. Liutex is considered as the vortex definition in this chapter.

Since Liutex is a vector-based concept that provides the direction as well as the magnitude of the rigid rotation, we can integrate the Liutex vector passing through a point within the vortex region to obtain a special Liutex core line. The concept of Liutex core line was first proposed in [34]. In short, Liutex core lines are formed by points where local maximal Liutex magnitude is reached. When calculating Liutex core lines using computer programs, some numerical difficulties are raised. Because of the numerical errors, the points selected by Liutex core line definition are a cluster of points rather than a single point. To solve this problem, Xu et al. [35] proposed seed point selection method. Li et al. [36] further improved this method by choosing the best line among the candidate lines generated by seed points. In this chapter, a new method, from another perspective, is reported to calculate the Liutex core lines based on the idea that the forward Liutex magnitude gradient lines converge to the Liutex core lines. So, instead of drawing two-directional Liutex magnitude lines, we only draw the forward lines. This leads to one advantage that the converging lines are unique so that people do not need to select one line from several candidate lines. This method achieves preliminary success for identifying Liutex core lines of the Lambda vortex in early boundary layer transition.

This chapter is organized as follows: The definitions of Liutex and Liutex core lines are introduced in Section 2 and 3, respectively. In Section 4, the manual method and Xu’s algorithm are reviewed. Then, the new proposed algorithm is introduced in Section 5. In Section 6, the proposed automatic method is applied to find the Liutex core lines from the data obtained from a Direct Numerical Simulation (DNSUTA) [37] of flow transition in a boundary layer validated by researchers from UTA and NASA Langley. The results of the DNSUTA were also compared to other DNS results where the DNSUTA was shown to have remarkable consistency with [38]. Conclusions are made in Section 7.

2. Liutex

Liutex $\overrightarrow{R}$ is a vector such that

$$\overrightarrow{R}=R\overrightarrow{r}$$E1

where its direction $\overrightarrow{r}$ is the eigenvector of the velocity gradient tensor $\mathit{grad}\overrightarrow{v}$ that satisfies $\overrightarrow{\omega }·\overrightarrow{r}>0$ where $\overrightarrow{\omega }$ is the vorticity and the magnitude R can be calculated by the following equation:

$$R=\left(\overrightarrow{\omega }·\overrightarrow{r}\right)-\sqrt{{\left(\overrightarrow{\omega }·\overrightarrow{r}\right)}^2-4{{\lambda }_{\mathit{ci}}}^2}$$E2

$\overrightarrow{r}$ is the real eigenvector of the velocity gradient tensor matrix representing local rotation axis, and R refers to the twice angular speed.

3. Liutex core line

The Liutex core line [34] (or the vortex rotation axis) can be defined as a Liutex line passing through the points that satisfy the condition:

$$\nabla R\times \overrightarrow{r}=0,R>0$$E3

where $\nabla R$ represents the Liutex magnitude gradient vector and $\overrightarrow{r}$ represents the direction of the Liutex vector. The Liutex core lines can be colored according to the vortex strength (Liutex magnitude) and provide a unique representation of the vortex structure.

The Liutex core concept comes from the idea that Liutex core is located at the position where Liutex reaches the maximum in the plane perpendicular to the local Liutex direction as shown in Figure 1. Since it is the local maxima in the plane, the projection of the Liutex gradient in the plane is zero, and thus, its only possible direction is perpendicular to the plane, which is parallel to the local Liutex direction. So, Eq. (3) is used to define Liutex core.

Figure 1.

Liutex magnitude gradient lines distribution for the local maximal Liutex magnitude in the plane perpendicular to the Liutex direction.

4. Manual method and Xu’s algorithm to calculate Liutex core lines

After defining the Liutex core line in [34], a manual method was proposed simultaneously in the same paper based on visualization software, for instance, Tecplot, Paraview, and so forth. After loading data files into these software, slices, iso-surfaces, and streamlines can be drawn, which are used in the manual method. The steps of the manual method are listed below.

Manual method:

Step 1: Draw the iso-surface of Liutex to get a general idea of the vortex structure.

Step 2: Set a slice that exhibits the contour of Liutex magnitude.

Step 3: Determine the concentration of line of some Liutex magnitude gradient lines passing through the slice.

Step 4: Find the intersection point of the concentration line and the slice and create a Liutex line passing through the intersection point.

Use Burgers vortex as an example. The Bugers vortex is an analytical solution of the Navier-Stokes governing equations and is a typical widely used vortex model. Its velocity field can be described in the cylindrical coordinates as follows:

$$v_r=-\mathit{\alpha r}$$E4

$$v_z=2\mathit{\alpha z}$$E5

$$v_{\theta }=\frac{\mathrm{\Gamma }}{2\mathit{\pi r}}g\left(r\right)$$E6

where $\alpha >0$ and $\mathrm{\Gamma }>0$ are constants, and $g\left(r\right)$ can be expressed as

$$g\left(r\right)=1-e^{\left(-\frac{\alpha r^2}{2\nu }\right)}$$E7

Some streamlines of the Burgers vortex is shown in Figure 2. It can be seen that Burgers vortex is like an upside-down cyclone where the flow converges into the center and then stays in the center and is stretched. As stated in the step 1, the iso-surface of Liutex is shown in Figure 3, which shows there is a vortex center inside the iso-surface. Then, according to step 2 and step 3, a slice cutting the iso-surface is created, and gradient of Liutex magnitude lines is drawn, as shown in Figure 4. Clearly, all gradient lines converge to one point, and that point is the seed point. Next, we draw the Liutex line passing through the seed point, and this line represents the vortex center, and the result is shown in Figure 5.

Figure 2.

Streamlines of the Burgers vortex.

Figure 3.

Liutex iso-surface of the Burgers vortex.

Figure 4.

Liutex magnitude distribution on a slice with Liutex gradient lines.

Figure 5.

Liutex core line of the Burgers vortex.

The Burgers vortex is a simple vortex model. Then, a more complicated case is tested. We will use the direct numerical simulation of the flat plate boundary layer transition. Similarly, we can first draw the iso-surface of Liutex (shown in Figure 6) and find the vortex structures are much more complicated than the Burgers vortex. Next, we display Liutex gradient lines and find the intersection point of the lines and the slice as shown in Figure 7. Then, we draw the Liutex line passing through that point. In Figure 8, colors are used to represent the Liutex magnitude and we can see the rotation strengths.

Figure 6.

Liutex iso-surface of the flat plate boundary layer transition.

Figure 7.

Liutex iso-surface, a slice, Liutex gradient lines, and the intersection point.

Figure 8.

Liutex iso-surface and Liutex core line.

Obviously, the manual method is low in efficiency and impractical for generating all Liutex core lines for complicated vortex. So, many scientists aim to find an automatic way. Li and Xu et al. proposed an algorithm to automatically calculate Liutex core lines in [36].

Xu’s algorithm:

Step 1: Read position and velocity data.

Step 2: Compute Liutex vector.

Step 3: Search vortex core points by the criterions $\left|\nabla R\times \overrightarrow{r}\right|<\epsilon$ and ${\mathrm{\Omega }}_R>{\mathrm{\Omega }}_0$ where ${\mathrm{\Omega }}_R$ represents relative rotation strength and can be calculated by

$${\mathrm{\Omega }}_R=\frac{{\left(\overrightarrow{\omega }·\overrightarrow{r}\right)}^2}{2\left[{\left(\overrightarrow{\omega }·\overrightarrow{r}\right)}^2-2{\lambda }_{\mathit{ci}}^2+2{\lambda }_{\mathit{cr}}^2+{\lambda }_r^2\right]+\epsilon }$$E8

where $\epsilon$ is a small positive number,${\lambda }_{\mathit{cr}}$ and ${\lambda }_{\mathit{ci}}$ are the real and complex eigenvalues of the velocity gradient tensor, and ${\lambda }_r$ is the real eigenvalue of the velocity gradient tensor.

Step 4: Compute vortex core lines based on selected points in step 3.

Step 5: Group the vortex core lines.

Step 6: Delete the false and inaccurate vortex core lines by the third criterion.

Step 7: Output data of vortex core lines.

To make it simple, Xu’s algorithm indeed proposes a series of criterions on how to select the seed points. His method first selects some seed point candidates by applying some criterions. And then, it groups the seed point candidates of the same vortex together. After that, it selects the best one from the seed point candidates in each group.

5. New algorithm to calculate Liutex core line

As explained in Section 3, Liutex direction is parallel to Liutex magnitude gradient direction at vortex core positions, so theoretically, Liutex lines and Liutex gradient lines will be the same in the location where the Liutex core lines pass. However, in practice, due to numerical errors, the points that the numerical program finds will not be exactly located on the same single Liutex line, which requires selecting the best one by a specified algorithm in Xu’s paper. Another perspective to solve this problem is to make use of the property of gradient lines. Gradient lines have a good property that they can converge to the local maximums. In Figure 9, people can see that even though the seed points are much away from the vortex center, the generated Liutex gradient lines will converge to a unique line finally. So, even though selected points have numerical errors, we will still get a unique line. Based on this idea, a new algorithm is proposed as follows.

Figure 9.

Liutex iso-surface and Liutex gradient lines.

Step 1: Set the Liutex threshold and exclude the non-vortex points whose Liutex magnitudes are below the threshold.

Step 2: Find all seed points satisfying

$$\left|\nabla R\times \overrightarrow{r}\right|<\epsilon$$E9

where $\epsilon$ is a small number. Originally, the Liutex core line definition is $\nabla R\times \overrightarrow{r}=0$. Because computer has numerical errors, we set a small number $\epsilon$ rather than 0.

Step 3: Draw forward Liutex gradient lines starting from the selected seed points.

The proposed method has the following merits. First, this method does not have a high demand on the accuracy of the seed points. Xu’s method already shows that it is hard to get enough accurate seed points and to ensure a unique line, people need to pick one seed point from seed point candidates. But the proposed method can tolerate some errors on the seed point selection. Secondly, since the proposed method can tolerate some errors on the seed point selection, the parameters in this method are not case-sensitive because people do not rely on the criterions to filter out all fake seed points.

6. Test case

In this section, the proposed automatic method is applied to find the Liutex core line for the data obtained from a Direct Numerical Simulation [6] (DNSUTA) of flow transition in a boundary layer, which has been validated by researchers from UTA and NASA Langley. The mesh of this simulation includes $1920\times 128\times 241$ points in streamwise (x), spanwise (y), and wall normal (z) directions. The parameters of the flow are shown in Table 1. $M_{\infty }$ is the inflow Mach number; Re represents Reynolds number, which is defined as $\mathit{Re}=\frac{{\rho }_{\infty }{u}_{\infty }{\delta }_{\mathit{\text{in}}}}{{\mu }_{\infty }}$ using the inflow boundary layer displacement thickness as the length reference; $x_{\mathit{\text{in}}}$ is the distance between the leading edge of the flat plate and upstream boundary of the computational domain; Lx and Ly are the lengths of the computational domain along x and y directions, respectively; $L{z}_{\mathit{\text{in}}}$ means the height at inflow boundary; $T_w$ and $T_{\infty }$ represent the temperatures of wall and free stream, respectively. The illustration of the computational domain is shown in Figure 10.

$M_{\infty }$	Re	$x_{\mathit{\text{in}}}$	Lx	Ly	$L{z}_{\mathit{\text{in}}}$	$T_w$	$T_{\infty }$
0.5	1000	300.79${\delta }_{\mathit{\text{in}}}$	$798.03{\delta }_{\mathit{\text{in}}}$	$22{\delta }_{\mathit{\text{in}}}$	$40{\delta }_{\mathit{\text{in}}}$	273.15 K	273.15 K

Table 1.

Flow parameters.

Figure 10.

Computational domain.

This simulation uses the 3D compressible Navier–Stokes equations in curvilinear coordinates, which can be expressed as the following,

$$\frac{1}{\mathrm{J}}\frac{\partial \mathbf{Q}}{\mathrm{\partial t}}+\frac{\partial \left(\mathbf{E}-{\mathbf{E}}_{\mathrm{v}}\right)}{\mathrm{\partial \xi }}+\frac{\partial \left(\mathbf{F}-{\mathbf{F}}_{\mathrm{v}}\right)}{\mathrm{\partial \eta }}+\frac{\partial \left(\mathbf{H}-{\mathbf{H}}_{\mathrm{v}}\right)}{\mathrm{\partial \zeta }}=0$$E10

$$\mathit{Q}\mathbf{=}\left(\begin{array}{c}\rho \\ \mathit{\rho u}\\ \mathit{\rho v}\\ \mathit{\rho w}\\ e\end{array}\right),\mathit{E}\mathbf{=}\frac{1}{J}\left(\begin{array}{c}\mathit{\rho U}\\ \mathit{\rho uU}+p{\xi }_x\\ \mathit{\rho vU}+p{\xi }_y\\ \mathit{\rho wU}+p{\xi }_z\\ \left(e+p\right)U\end{array}\right),\mathit{F}\mathbf{=}\frac{1}{J}\left(\begin{array}{c}\mathit{\rho V}\\ \mathit{\rho uV}+p{\eta }_x\\ \mathit{\rho vV}+p{\eta }_y\\ \mathit{\rho wV}+p{\eta }_z\\ \left(e+p\right)V\end{array}\right)$$E11

$$\mathit{H}\mathbf{=}\frac{1}{J}\left(\begin{array}{c}\mathit{\rho W}\\ \mathit{\rho uW}+p{\zeta }_x\\ \mathit{\rho vW}+p{\zeta }_y\\ \mathit{\rho wW}+p{\zeta }_z\\ \left(e+p\right)W\end{array}\right),{\mathit{E}}_v\mathbf{=}\frac{1}{J}\left(\begin{array}{c}0\\ {\tau }_{\mathit{xx}}{\xi }_x+{\tau }_{\mathit{yx}}{\xi }_y+{\tau }_{\mathit{zx}}{\xi }_z\\ {\tau }_{\mathit{xy}}{\xi }_x+{\tau }_{\mathit{yy}}{\xi }_y+{\tau }_{\mathit{zy}}{\xi }_z\\ {\tau }_{\mathit{xz}}{\xi }_x+{\tau }_{\mathit{yz}}{\xi }_y+{\tau }_{\mathit{zz}}{\xi }_z\\ q_x{\xi }_x+q_y{\xi }_y+q_z{\xi }_z\end{array}\right)$$E12

$${\mathit{F}}_v\mathbf{=}\frac{1}{J}\left(\begin{array}{c}0\\ {\tau }_{\mathit{xx}}{\eta }_x+{\tau }_{\mathit{yx}}{\eta }_y+{\tau }_{\mathit{zx}}{\eta }_z\\ {\tau }_{\mathit{xy}}{\eta }_x+{\tau }_{\mathit{yy}}{\eta }_y+{\tau }_{\mathit{zy}}{\eta }_z\\ {\tau }_{\mathit{xz}}{\eta }_x+{\tau }_{\mathit{yz}}{\eta }_y+{\tau }_{\mathit{zz}}{\eta }_z\\ q_x{\eta }_x+q_y{\eta }_y+q_z{\eta }_z\end{array}\right),{\mathit{E}}_v\mathbf{=}\frac{1}{J}\left(\begin{array}{c}0\\ {\tau }_{\mathit{xx}}{\zeta }_x+{\tau }_{\mathit{yx}}{\zeta }_y+{\tau }_{\mathit{zx}}{\zeta }_z\\ {\tau }_{\mathit{xy}}{\zeta }_x+{\tau }_{\mathit{yy}}{\zeta }_y+{\tau }_{\mathit{zy}}{\zeta }_z\\ {\tau }_{\mathit{xz}}{\zeta }_x+{\tau }_{\mathit{yz}}{\zeta }_y+{\tau }_{\mathit{zz}}{\zeta }_z\\ q_x{\zeta }_x+q_y{\zeta }_y+q_z{\zeta }_z\end{array}\right)$$E13

where $J$ is the Jacobian of the coordinate transformation; ${\xi }_x,{\xi }_y,{\xi }_z,{\eta }_x,{\eta }_y,{\eta }_z,{\zeta }_x,{\zeta }_y,\mathit{\text{and}}{\zeta }_z$ are coordinate transformation metrics; $\mathit{Q}$ is the vector of conserved quantities; $\mathit{E},\mathit{F},\mathit{\text{and}}\mathit{H}$ are the inviscid flux vectors and ${\mathit{E}}_{\mathit{v}},{\mathit{F}}_{\mathit{v}},\mathit{\text{and}}{\mathit{H}}_{\mathit{v}}$ are the viscous vectors; and ${\tau }_{\mathit{xx}},{\tau }_{\mathit{yx}},{\tau }_{\mathit{zx}},{\tau }_{\mathit{xy}},{\tau }_{\mathit{yy}},{\tau }_{\mathit{zy}},{\tau }_{\mathit{xz}},{\tau }_{\mathit{yz}},\mathit{\text{and}}{\tau }_{\mathit{zz}}$ are viscous stress components.

Compact schemes are used to do the spatial discretization. It has the following form:

$${\beta }_{-}{f}_{j-2}^{\prime }+{\alpha }_{-}{f}_{j-1}^{\prime }+f_j^{\prime }+{\alpha }_{+}{f}_{j+1}^{\prime }+{\beta }_{+}{f}_{j+2}^{\prime }=\frac{1}{h}\left(b_{-}{f}_{j-2}+a_{-}{f}_{j-1}+c{f}_j+a_{+}{f}_{j+1}+b_{+}{f}_{j+2}\right)$$E14

where $f_j^{\prime }$ is the derivative at point j. For the sixth-order scheme used in this simulation,

$${\beta }_{-}=0,{\alpha }_{-}=\frac{1}{3},{\beta }_{+}=0,{\alpha }_{+}=\frac{1}{3},b_{-}=-\frac{1}{36},a_{-}=-\frac{7}{9},a_{+}=\frac{7}{9},b_{+}=\frac{1}{36}$$E15

The total variation diminishing (TVD) third order Runge–Kuntta method is used for the time discretization. The equations are:

$$Q^{\left(0\right)}=Q^{\left(n\right)}$$E16

$$Q^{\left(1\right)}=Q^{\left(0\right)}+\mathrm{\Delta }t{R}^{\left(0\right)}$$E17

$$Q^{\left(2\right)}=\frac{3}{4}{Q}^{\left(0\right)}+\frac{1}{4}{Q}^{\left(1\right)}+\frac{1}{4}\mathrm{\Delta }t{R}^{\left(1\right)}$$E18

$$Q^{\left(n+1\right)}=\frac{1}{3}{Q}^{\left(0\right)}+\frac{2}{3}{Q}^{\left(2\right)}+\frac{2}{3}\mathrm{\Delta }t{R}^2$$E19

The simulation result is shown in Figure 11. Spanwise vortex appears first followed by the Lambda vortex. We will use the Lambda vortex region to test the proposed method.

Figure 11.

DNS result of the flat plate boundary layer transition.

Applying step 1 and step 2 of the proposed method, we can get the seed points as shown in Figures 12–14. Drawing forward Lituex gradient lines starting from these seed points, we get Figures 15–17. If we draw Liutex lines from these seed points, we will get several separate lines rather than a unique line as shown in Figures 18–20. The proposed new algorithm obtains the unique vortex core line and does not need to select the best line from a group of lines generated by the seed points. Even though the obtained seed points have numerical errors, the output will still be convergent to a unique line. Figures 21–23 show Liutex magnitudes at the vortex core positions by colors. They show the rotation strength of the vortex core line is not uniform and is varied at different locations. This is what people cannot get by using the iso-surface methods as iso-surface is a surface where the magnitudes are equal to the threshold.

Figure 12.

Selected seed points by step 1 and step 2.

Figure 13.

(top view) Selected seed points by step 1 and step 2.

Figure 14.

(side view) Selected seed points by step 1 and step 2.

Figure 15.

Selected seed points and corresponding forward Liutex gradient lines.

Figure 16.

(top view) Selected seed points and corresponding forward Liutex gradient lines.

Figure 17.

(side view) Selected seed points and corresponding forward Liutex gradient lines.

Figure 18.

Liutex lines generated by the same seed points.

Figure 19.

(top view) Liutex lines generated by the same seed points.

Figure 20.

(side view) Liutex lines generated by the same seed points.

Figure 21.

Liutex iso-surface and colorful forward Liutex gradient lines.

Figure 22.

(top view) Liutex iso-surface and colorful forward Liutex gradient lines.

Figure 23.

(side view) Liutex iso-surface and colorful forward Liutex gradient lines.

7. Conclusion

This study reports a method to automatically draw Liutex core lines from a new perspective. The proposed method is based on the property of gradient lines that they are converging toward to the local maximum points. Instead of Liutex lines, Liutex gradient lines are used in the proposed method and give this method two major merits. Firstly, the uniqueness comes from the property of gradient lines. Secondly, the proposed method can tolerate some numerical errors as these errors can easily be eliminated during the process that gradient lines converge to the center. However, this is not a mature method, and some work still needs to be done in the future, for example, developing it into a universal method that works for different cases.

Acknowledgments

The project is currently supported by NSF under Grant #2300052. The authors are thankful for the support by the UTA Department of Mathematics, which houses the UTA Vortex and Turbulence Research Team. The authors are also grateful to Texas Advanced Computing Center (TACC) for providing computation time.

Data availability

The data that supports the findings of this study are available from the corresponding author upon reasonable request.