Generative prediction of flow field based on the diffusion model

Abstract

We propose a geometry-to-flow diffusion model that utilizes the input of obstacle shape to predict a flow field past the obstacle. The model is based on a learnable Markov transition kernel to recover the data distribution from the Gaussian distribution. The Markov process is conditioned on the obstacle geometry, estimating the noise to be removed at each step, implemented via a U-Net. A cross-attention mechanism incorporates the geometry as a prompt. We train the geometry-to-flow diffusion model using a dataset of flows past simple obstacles, including the circle, ellipse, rectangle, and triangle. For comparison, the CNN model is trained using the same dataset. Tests are carried out on flows past obstacles with simple and complex geometries, representing interpolation and extrapolation on the geometry condition, respectively. In the test set, challenging scenarios include a cross and characters ‘PKU’. Generated flow fields show that the geometry-to-flow diffusion model is superior to the CNN model in predicting instantaneous flow fields and handling complex geometries. Quantitative analysis of the model accuracy and divergence in the fields demonstrate the high robustness of the diffusion model, indicating that the diffusion model learns physical laws implicitly.

1 Introduction

Understanding the interaction between fluid and solid objects is crucial for optimizing various engineering systems (Shelley & Zhang, 2011). Computational fluid dynamics (CFD) plays a vital role in revealing the complex fluid-solid interactions (Tong et al., 2021) and enhancing design efficiency. However, CFD faces several challenges, including the expertise required to generate high-quality computational meshes, especially for complex geometries (Perot, 2011), and the substantial computational resources needed for high-fidelity simulations (Kern et al., 2024). These challenges have motivated the exploration of CFD-accelerating approaches.

Machine learning (ML) has emerged as a promising solution to address the limitations of traditional CFD methods (Kochkov et al., 2021). By leveraging the data-driven models, ML has demonstrated the potential to significantly accelerate simulations (Vinuesa & Brunton, 2022). Various ML approaches have been applied to CFD successively, including super-resolution of flow fields (Fukami et al., 2019), estimating forces of a plate in flows (Tong et al., 2022), physics-informed neural network (PINN) for flow simulation (Chu et al., 2022), velocity interpolation in multiphase flows (Hu et al., 2024). Despite these advancements, many existing ML methods struggle with extrapolation (Wu et al., 2020). The deficit of ML methods on generalization to out-of-distribution scenarios poses a challenge for their application in CFD (Fukami et al., 2024), where extrapolation on flow parameters and solid geometry is often expected.

Generative artificial intelligence has recently revolutionized many fields, particularly large language models (Biever, 2023) and image generation (Pan et al., 2023). Generative models aim to learn the underlying probability distribution of a given dataset. Once trained, they generate outputs by sampling from the distribution. They have been applied to scientific research, such as the material design (Fu et al., 2023) and active control (Kim et al., 2024). Examples in CFD include the super-resolution of turbulence field via a generative adversarial network (Kim et al., 2021) and flow field prediction using a variational autoencoder (Ranade et al., 2021). However, most existing works still focus on interpolation within the conditions listed in the training set, lacking further examination of out-of-distribution generalization.

The diffusion model (Yang et al., 2023) has emerged as a powerful probabilistic method, with predominant formulations including the denoising diffusion probabilistic model (DDPM) (Sohl-Dickstein et al., 2015) and score-based generative models (Song & Ermon, 2020). The diffusion model progressively perturbs data distribution to the Gaussian by injecting noise, then learns to reverse this process for sample generation. Shan et al. (2024) raised a physics-informed diffusion model for super-resolution. Qiu et al. (2024) trained physics-informed diffusion model to output time evolution of flow past a cylinder and inside veins. The ability of the diffusion model to learn distributions and its potential for extrapolation makes it a promising approach for accelerating CFD and exploring complex flow phenomena around bodies with diverse geometries.

In this work, we mainly focus on the two-dimensional (2D) flow past an obstacle. We develop a geometry-to-flow (G2F) diffusion model in which the obstacle geometry serves as a prompt to guide the gradual denoising process, generating the flow field around obstacles. The G2F diffusion model is trained with elementary geometries. Evaluations using simple and complex geometries validate the model in interpolation and extrapolation tasks. Our results demonstrate the model’s capability to generate flow fields for a wide range of obstacle shapes, including those outside the training distribution.

2 Geometry prompt to flow field output

2.1 DDPM

We use the obstacle geometry as a prompt to generate a flow past the obstacle as output through the DDPM. Figure 1 provides a brief overview of this G2F diffusion model. There are two diffusion processes in the DDPM – the forward diffusion perturbs data to noise, and the reverse diffusion converts noise back to data (Sohl-Dickstein et al., 2015). The DDPM learns a probabilistic representation of the underlying data, thereby facilitating the generation of high-quality samples that resemble the target distribution.

For data $\boldsymbol{z}_{0}$ sampled from a distribution $q(\boldsymbol{z}_{0})$, the forward diffusion utilizes a Markov transition kernel

$$q\left(\boldsymbol{z}_{t}|\boldsymbol{z}_{t-1}\right)=\mathcal{N}\left(\boldsymbol{z}_{t};\sqrt{1-\beta_{t}}\boldsymbol{z}_{t-1},\beta_{t}\boldsymbol{I}\right),\quad t=1,\dots,T$$ (1)

to inject noise, generating a sequence of noisy samples, $\boldsymbol{z}_{1},\dots,\boldsymbol{z}_{T}$ in total $T$ steps, where $\beta_{t}\in\left(0,1\right)$ is the diffusion rate, $\mathcal{N}$ denotes the Gaussian distribution, and $\boldsymbol{I}$ is the identity matrix. For $T\rightarrow\infty$, $q\left(\boldsymbol{z}_{T}\right)\approx\mathcal{N}\left(\boldsymbol{0},\boldsymbol{I}\right)$, allowing to approximate $\boldsymbol{z}_{T}$ as a Gaussian vector.

Figure 1: Schematic for the reverse diffusion process in the G2F diffusion model. It generates a sample by removing noises in a series of steps through a U-Net. The geometry prompt is injected into the U-Net via a cross-attention mechanism, which takes the geometry as the prompt to control the generation process.

To generate samples $\boldsymbol{z}_{0}$, DDPM conducts a reverse diffusion process on the Gaussian vector $\boldsymbol{z}_{T}$ using a learnable Markov transition kernel

$$p_{\boldsymbol{\theta}}\left(\boldsymbol{z}_{t-1}|\boldsymbol{z}_{t};\boldsymbol{c}\right)=\mathcal{N}\left[\boldsymbol{z}_{t-1};\boldsymbol{\mu}_{\boldsymbol{\theta}}\left(\boldsymbol{z}_{t},t\right),\boldsymbol{\Sigma}_{\boldsymbol{\theta}}\left(\boldsymbol{z}_{t},t\right),\boldsymbol{c}\right],$$ (2)

where $\boldsymbol{\theta}$ denotes model parameters, $\boldsymbol{c}$ denotes the condition, the mean $\boldsymbol{\mu}_{\boldsymbol{\theta}}$ and variance $\boldsymbol{\Sigma}_{\boldsymbol{\theta}}$ are parameterized by the neural network (NN). The parameters $\boldsymbol{\theta}$ are trained to make the reverse diffusion process a good approximation of the forward one, minimizing the Kullback-Leibler divergence (Kullback & Leibler, 1951) between $q\left(\boldsymbol{z}_{0},\boldsymbol{z}_{1},\dots,\boldsymbol{z}_{T}\right)$ and $p_{\boldsymbol{\theta}}\left(\boldsymbol{z}_{0},\boldsymbol{z}_{1},\dots,\boldsymbol{z}_{T};\boldsymbol{c}\right)$. To optimize training, Ho et al. (2020) proposed a loss function based on the mean-square-error (MSE) of the predicted noise,

$$\mathcal{L}=\mathbb{E}_{t,\boldsymbol{z}_{0},\boldsymbol{\epsilon}}\left[\left\|\boldsymbol{\epsilon}-\boldsymbol{\epsilon}_{\boldsymbol{\theta}}\left(\boldsymbol{z}_{t},t;\boldsymbol{c}\right)\right\|^{2}\right],$$ (3)

where $\mathbb{E}$ denotes the expectation, and $\boldsymbol{\epsilon}_{\boldsymbol{\theta}}\left(\boldsymbol{z}_{t},t;\boldsymbol{c}\right)$ is a NN predicting the noise vector $\boldsymbol{\epsilon}\sim\mathcal{N}\left(\boldsymbol{0},\boldsymbol{I}\right)$.

2.2 Geometry prompt

In the G2F diffusion model sketched in figure 1, we introduce the geometry prompt as the condition $\boldsymbol{c}$ in (2), and use a U-Net (Ronneberger et al., 2015) with a cross-attention mechanism (Rombach et al., 2022) to realize $\boldsymbol{\epsilon}_{\boldsymbol{\theta}}(z_{t},t;\boldsymbol{c})$.

The U-Net consists of a shrinking path, an expanding path, and some skip connections between the two. The shrinking path gradually compresses a low-dimensional, high-resolution field to a high-dimensional, low-resolution field. It contains five layers, with sizes of (64, 64, 256), (32, 32, 512), (16, 16, 1024), (8, 8, 1024), and (4, 4, 2048), followed by an average pool layer to obtain encoded data with size of (1, 1, 2048). The expanding path follows the reverse process.

To incorporate the geometry prompt, we first encode the geometry information into a vector of 0 and 1, representing the fluid and solid, respectively. The prompt is encoded by two NNs, having the same sizes as the last and the second-to-last layers in the expanding path. Outputs from these two NNs are then combined into the last and the second-to-last layers in the expanding path. The cross-attention (Rombach et al., 2022) in the last two layers balances the interpolation and extrapolation in our practice.

Note that DDPM generates results from random noises, which means that DDPM’s results differ even with the same prompt. Here we focus on the geometry prompt, so we only refer the generated flow fields to a geometry, without specifying a particular time.

2.3 Datasets

We evaluate the G2F diffusion model based on 2D flows past various obstacles. The dataset includes a range of obstacle geometries, as illustrated in figure 2. The G2F diffusion model is trained using the flow field around simple geometries, including the circle, ellipse, rectangle, and triangle. They are often used as building blocks to construct complex geometries. To assess the model’s ability to extrapolate to different geometries, the test set includes flows around obstacles not present in the training set. This includes several cases with parallelogram obstacles. Furthermore, the model is tested on a cross and characters ‘PKU’ with convex shapes.

Figure 2: Datasets for evaluating the G2F diffusion model. Elementary obstacle geometries including the circle, ellipse, rectangle, and triangle are used for training; Complex obstacle geometries including the parallelogram, cross, and characters ‘PKU’ are employed for testing. Each sample $\boldsymbol{z}_{0}=\left[u,v,p,\omega\right]$ includes velocity, pressure, and vorticity.

The geometry complexity is characterized by the roundness $r=4\pi S/C^{2}$, where $S$ and $C$ are the area and circumference, respectively. The roundness is bounded by $r=1$ for the circle and $r\rightarrow 0$ for extremely complex geometry. The cross in the test set has $r=0.321$, smaller than all shapes in the training set, and the characters ‘PKU’ have $r\approx 0$.

Focusing on the geometry condition, we evaluate the model using flow field data with the same free-stream velocity $U_{\infty}=1$, viscosity $\mu=1\times 10^{-2}$, and density $\rho=1$. The circle radius is 1, and the characteristic length of the obstacles ranges from 1 (a rectangle or ellipse with its longitudinal axis aligned in the streamwise direction) to 4 (a rectangle or ellipse with its longitudinal axis aligned in the spanwise direction). Correspondingly, the Reynolds number $\Rey$ varies from 100 to 400. We solve the incompressible Navier-Stokes equations to obtain the flow field around various obstacles as the ground truth (GT) using OpenFOAM (Greenshields, 2022). The body-fitted meshes with 55000 cells are employed for a computational domain of $\left(-20,60\right)\times\left(-20,20\right)$, where the geometric centroid of the obstacle is positioned at the spatial location $(x,y)=(0,0)$. Each flow-field snapshot has the velocity $(u,v)$, pressure $p$, and vorticity $\omega$. Note that we interpolate the flow field to a $128\times 128$ uniform grid over a view box of $\left(-2,10\right)\times\left(-4,4\right)$ to obtain $z_{0}=[u,v,p,\omega]$ for the model. For each case, 50 snapshots are sampled with the time step of $1$.

In the implementation, we empirically set $T=400$ and $\beta_{t}=1\times 10^{-4}+5\times 10^{-5}t$. The parameters $\boldsymbol{\theta}$ in the NN are optimized to minimize the loss function in (3). The parameters are initialized by the Kaiming method (He et al., 2015) and updated using the Adam optimizer (Kingma & Ba, 2015). The initial learning rate is 0.001, subsequently reduced to one-tenth of its preceding value after every 2000 epochs.

In comparison, we trained a CNN model using the same U-Net and dataset. The CNN model loads geometry prompts at the entrance of the shrinking path and outputs flow fields at the exit of the expanding path, similar to Shirzadi et al. (2022). The CNN parameters are optimized to minimize the MSE of outputs against the GT. We specifically chose to compare with CNN as it utilizes the same U-Net architecture as the G2F diffusion model, allowing for a direct comparison of the two approaches. The physics-informed model, such as PINN (Karniadakis et al., 2021), is not compared here because the presented diffusion model does not incorporate explicit physical constraints. The physics-informed diffusion model is discussed in supplementary material and will be further studied within the same network structure in future work.

3 Results

We evaluate the G2F diffusion model against the obstacles in the training and test sets. Given the random nature of the model, where a specific time is not designated, a snapshot in the GT solution is selected for comparison with the model output. We measure the similarity between the model output and GT using the $L_{1}$ error

$$L_{1}\left(\boldsymbol{z}_{0}^{M}\right)=\left\|\boldsymbol{z}_{0}^{M}-\boldsymbol{z}_{0}^{\mathrm{GT}}\right\|,$$ (4)

where the superscript $M$ denotes the model output. Then the GT snapshot with the minimal $L_{1}\left(\boldsymbol{z}_{0}^{M}\right)$ is used for comparison below.

3.1 Training set

We first present the flow past a cylinder as a standard ML-enhanced CFD case. Figure 3 compares the G2F diffusion and CNN models against the GT. At $\Rey=200$, the flow is featured by periodic vortex shedding.

\begin{overpic}[height=126.47249pt]{pic/round/u.pdf} \put(-8.0,95.0){(a)} \end{overpic}

\begin{overpic}[height=126.47249pt]{pic/round/u-plot.pdf} \put(0.0,95.0){(b)} \end{overpic}

\begin{overpic}[height=126.47249pt]{pic/round/p.pdf} \put(-8.0,95.0){(c)} \end{overpic}

\begin{overpic}[height=126.47249pt]{pic/round/p-plot.pdf} \put(0.0,95.0){(d)} \end{overpic}

Figure 3: Flow past a cylinder generated by the G2F diffusion model and CNN, compared with the GT: (a) contour of $u$, (b) profiles of $u$ in the wake at $x=0$, 2, and 6, (c) contour of pressure $p$, and (d) pressure on the obstacle surface. The vertical dashed lines in the upper left panel mark the stations at $x=0$, 2, and 6.

Figures 3(a) and (c) compare the instantaneous contours of the streamwise velocity $u$ and pressure $p$, respectively. Both models perform well in the stagnation region. The G2F diffusion model captures the von Kármán vortex street in the wake. The model output well approximates the alternating shear layers and low-pressure zones, although it exhibits minor discrepancies in values. In contrast, the CNN model fails to generate the asymmetric flow structures, displaying a ‘time-averaged’ flow field. This is because the diffusion model learns the distribution $q\left(\boldsymbol{z}_{0}\right)$ of the data $\boldsymbol{z}_{0}$, where a sample from the distribution approximates one snapshot of the flow field. Meanwhile, the CNN is trained to minimize the distance between its output and all snapshots, resulting in an ‘averaged’ result.

We plot the profiles of $u$ at $x=0$, 2, and 6 in figure 3(b). The G2F diffusion model agrees with the GT on the shift of velocity peaks and valleys along the development of the wake. However, quantitative errors remain in the model output. Note that the G2F diffusion model returns negative inside the obstacle, where $u$ should be 0. Due to averaging over all snapshots, the CNN model gives symmetric $u$ profiles at all three stations.

We further compare the pressures on the obstacle from different models in figure 3(d). The pressure on the obstacle surface is plotted against the angle in the polar coordinate, as illustrated in figure 3(c). The G2F diffusion model follows the GT, capturing the peak at the stagnation point and the subsequent asymmetric pressure drop on both sides. It also accurately generates the pressure rises at $\varphi=\pi/2$ and $3\pi/2$, indicating a good prediction on the boundary layer separation. Although the CNN gives reasonable results on the stagnation point and the local minimum position, its solution includes strong unphysical oscillations.

Overall, the G2F and CNN models exhibit distinct behaviours in predicting the flow field. The G2F diffusion model produces results that qualitatively match the instantaneous flow field, while the CNN model provides a time-averaged flow field. Similar results are observed in other geometries. This difference arises from the goals of the two deep-learning NNs. The G2F diffusion model is trained to minimize the KL divergence, equivalent to maximizing the likelihood. Sampling from the distribution $p_{\boldsymbol{\theta}}(\boldsymbol{z}_{0};\boldsymbol{c})$ then returns data close to one snapshot in training data. The CNN model struggles with asymmetry because it minimizes the sum of errors between its output and all snapshots, leading to an averaged result.

3.2 Test set

We employ obstacle geometries absent from the training set to test the extrapolation capability of the G2F diffusion model, including parallelograms, a cross, and the characters ‘PKU’. The parallelograms exhibit similarities to the simple geometries in the training set. Furthermore, the cross and characters ‘PKU’ are out-of-distribution cases, characterized by combinations of elementary geometries and convex shapes.

Figure 4 compares the velocity and pressure fields generated by the G2F diffusion and CNN models around the cross with the GT. The G2F diffusion model generates the instantaneous flow field with shedding vortices and clear boundary layers. Stagnation zones on the upstream side of the cross are predicted correctly, implied by the high-pressure values on the obstacle surface. Meanwhile, the CNN model returns obvious oscillations around the obstacle, and the fields are averaged in the wake. Moreover, the CNN model incorrectly predicts two negative pressure regions near the front arm of the cross, implying that it treats the cross arms as standalone objects in the flow field. The G2F diffusion model addresses this challenge effectively, demonstrating the extrapolation capability of generative models.

\begin{overpic}[height=126.47249pt]{pic/cross/u.pdf} \put(-8.0,95.0){(a)} \end{overpic}

\begin{overpic}[height=126.47249pt]{pic/cross/u-plot.pdf} \put(0.0,95.0){(b)} \end{overpic}

\begin{overpic}[height=126.47249pt]{pic/cross/p.pdf} \put(-8.0,95.0){(c)} \end{overpic}

\begin{overpic}[height=126.47249pt]{pic/cross/p-plot.pdf} \put(0.0,95.0){(d)} \end{overpic}

Figure 4: Flow past a cross generated by the G2F diffusion model and CNN, compared with the GT: (a) contour of $u$, (b) profiles of $u$ in the wake at $x=$ 0, 2, and 6, (c) contour of pressure $p$, and (d) pressure on a unit circle crossing the four arms. The vertical dashed lines in the upper left panel mark the stations at $x=0$, 2, and 6.

Due to the two cavities downstream of the cross, there exists a local minimum following the upper edge of the cross, as shown around $y=1.5$ in the $u$ profile at $x=2$. This is not found for simple geometries in the training set. The G2F diffusion model indicates a weak shear layer at this position, although it cannot accurately reproduce this structure. In addition, the G2F diffusion model gives a narrower wake flow compared with the GT, due to two possible reasons. First, the cross-attention mechanism makes the geometry more influential on the flow field around the obstacle (Vaswani et al., 2017), while the far field is less affected by the geometry. Second, a lack of obstacle geometries with large widths in the $y$-direction in the training set makes it difficult for the G2F diffusion model to predict a wide wake flow.

Generating flow fields around the characters ‘PKU’ is challenging due to the sharp corners, convex curves, and separated components, none of which are included in the training set. Given the difficulty in simulating this flow with a simple mesh, we present the generated streamwise velocity $u$ and pressure $p$ by the two models in figure 5, without the GT. The G2F diffusion model demonstrates a remarkable extrapolation capability, capturing the major flow features, despite some artefacts such as the spurious flows through the characters ‘P’ and ‘U’ in the velocity field. Again, the CNN model introduces noise at character boundaries and provides time-averaged-like fields. Model performance on other obstacle geometries is provided in supplementary material.

Figure 5: Flow past characters ‘PKU’ generated by the G2F diffusion model and CNN, on contours of $u$ and $p$.

All geometries’ results are summarized in figure 6. The black and red symbols represent the geometries in the training and test sets, respectively. The horizontal axis measures the roundness $r$ of geometries. Note that all geometries in the training set have $r\geq 0.5$, while several geometries in the test set have $r<0.5$, indicating out-of-distribution cases.

\begin{overpic}[height=151.76854pt]{pic/extra/uv.pdf} \put(6.0,95.0){(a)} \end{overpic}

\begin{overpic}[height=151.76854pt]{pic/extra/div_e.pdf} \put(5.0,95.0){(b)} \end{overpic}

Figure 6: Comparison of results for different obstacle geometries: (a) normalized $L_{1}$ error between ML methods and GT and (b) divergence. It shows that the G2F diffusion model outperforms the CNN in understanding physical laws, and achieves good performances in extrapolation tasks.

We compare the accuracy of the two ML models in figure 6(a) using the $L_{1}$ error of $\left[u,v\right]$ normalized by the difference of the maximum and minimum values in the GT, $L_{1}([u,v])/([u,v]_{\max}^{\mathrm{GT}}-[u,v]_{\min}^{\mathrm{GT}})$. The pressure and vorticity are not included, considering the differences in the order of magnitude. Most of the G2F diffusion model’s results are better than CNN’s, especially in the test set, indicating that the diffusion model is more accurate and robust. We also notice that CNN performs well for some slender obstacles, as these cases have weaker vortex shedding. Consequently, the wake flows in these cases are close to time-averaged results by the CNN model.

We further evaluate the model realizability via the divergence $\mathcal{D}_{\mathrm{div}}=\left\|\partial_{x}u+\partial_{y}v\right\|$ averaged over the domain, which quantifies the discrepancy from the incompressible flow. Results in figure 6(b) demonstrate that $\mathcal{D}_{\mathrm{div}}$ increases overall with the decrease of $r$. For all cases, the G2F diffusion model yields smaller $\mathcal{D}_{\mathrm{div}}$ than the CNN model, and the increase rate of $\mathcal{D}_{\mathrm{div}}$ with the decrease in $r$ is smaller. This confirms the robustness of the G2F diffusion model across various geometries, also implying that incorporating a divergence-free constraint may not significantly contribute towards better performance. The physics-informed diffusion model is discussed in supplementary material. Additionally, flows with $\Rey$ closer to the training-set average (around 200) yield better results, suggesting that the data at a broader range of $\Rey$ should be incorporated in the training set to enhance the model performance.

4 Conclusion

We propose a G2F diffusion model for predicting flows past an obstacle with various geometries. The model employs a U-Net architecture with a cross-attention mechanism to incorporate geometry information as a prompt, guiding the reverse diffusion process to generate flow fields. A dataset including 2D flows with different geometries is established for training and testing. The G2F diffusion model is trained using flow fields for simple obstacle geometries and evaluated on both interpolation and extrapolation tasks.

We assess the G2F diffusion model using simple and complex obstacle geometries present and absent in the training set, respectively. The model results are compared with those from a CNN model and the GT. The results demonstrate that the diffusion model outperforms the CNN model, particularly in generating instantaneous flow fields and handling out-of-distribution geometries. In particular, the diffusion model successfully captures essential flow features while maintaining physical consistency, even in challenging scenarios. It reproduces instantaneous flow fields by learning data distribution, while the CNN model yields time-averaged data. Furthermore, the diffusion model’s ability to learn physical laws is evident from its consistent performance across various geometries, as quantified by the divergence.

Beyond the current preliminary study on the application of the diffusion model in CFD, there are several issues that deserve further investigation. In addition to the obstacle geometry, $\Rey$ is another key parameter, emphasizing the importance of diverse training data for improving the model’s generalization capabilities. Future research includes expanding the dataset to cover a wider range of Reynolds numbers and geometries, exploiting the scaling law, and enhancing the model’s ability to handle diverse flow conditions.

While results indicate that the diffusion model has the potential to learn the underlying physics implicitly, incorporating physics information could accelerate the sampling process, as discussed in supplementary material. Since the conditioned diffusion model suffers from the diversity problem (Gu et al., 2024), further exploration of condition implementation approaches is necessary to address this challenge. Additionally, investigating the application of the diffusion model in generating time-dependent flow field data can provide insights for understanding and predicting complex fluid dynamics.

[Funding]Numerical simulations were performed on the TH-2A supercomputer in Guangzhou, China. This work has been supported in part by the National Natural Science Foundation of China (Nos. 11925201, 52306126, 11988102, and 92270203), the National Key R&D Program of China (Grant No. 2020YFE0204200), and the Xplore Prize.

[Declaration of interests]The authors report no conflict of interest.

[Author ORCIDs]Jiajun Hu, https://orcid.org/0009-0005-5192-876X; Zhen Lu, https://orcid.org/0000-0002-6729-3771; Yue Yang, https://orcid.org/0000-0001-9969-7431

[Author contributions]Y.Y., Z.L, and J.H. designed research. J.H. and Z.L. preformed research. All the authors discussed the results and wrote the manuscript. All the authors have given approval for the manuscript.