Many researchers have made outstanding contributions in the field of traditional discretization methods of fluid flow simulation Bristeau et al. (1985); Ferziger et al. (2002); Peyret and Taylor (2012); Fletcher (2012); Toro (2013). Based on these theories, Molenkamp et al. Molenkamp (1968) calculated the numerical solution of the convection equation using various finite-difference approximations, and determined that only the Roberts–Weiss approximation convected the initial distribution correctly, but required a huge computational cost. Mikula et al. Mikula et al. (2014) proposed an inflow implicit/outflow explicit finite volume method based on finite volume space discretization and semi-implicit time discretization to solve advection equations. The basic idea is that outflows from cells are handled explicitly, and inflows are handled implicitly. The method achieved outstanding results in terms of stability and computational accuracy. Zhao et al. Zhao et al. (2019) proposed a new improved finite volume method for solving one-dimensional advection equations under the framework of the second-order finite volume method. The method first applied the scalar conservation law to the elements in the finite volume method (FVM) to ensure its conservation in time and space and to ensure advection (i.e., conservation of transport physical quantities); then the time integral values of adjacent grid boundaries are equalized; finally, the equation is established to obtain a numerical solution. Experiments showed that this method has better stability and fewer disspation than the traditional FVM and can maintain the accuracy of the solution. Akitoshi Takayasu et al. Takayasu et al. (2019) proposed a verification calculation method for one-dimensional advection equations with variable coefficients, which was based on spectral methods and semigroup theory. They mainly provided a method for verification calculation using the C ₀ semigroup on the complex sequence space l ²,which comes from the solution of the Fourier series. Experiments showed that the given strict error proved the correctness of the exact solution, and the solution has high precision and fast solution speed. Although traditional discretization method shave achieved high solution accuracy, they have the problem of high computational cost if outstanding solution accuracy is desired.

To solve the problem of high computational cost while calculating high-precision solutions in traditional discretization methods, researchers have proposed acceleration techniques to speed up the numerical discretization solution. Multigrid technology stood out among various approaches Dwyer et al. (1982); Brown (1982); Berger and Oliger (1984); Phillips and Schmidt (1984); Phillips and Schmidt (1985); Zhang (1997); Mazhukin et al. (1993); Jin et al. (2014). Among them, Brown et al. Brown (1982) used the multigrid mesh-embedding technique to solve three-dimensional transonic potential flow. They used small grids to model regions of large local gradients and large-scale grids to model regions with relatively small gradients. Their method improved the speed of solving equation discretization schemes. Phillips et al. Phillips and Schmidt (1984) proposed a multilevel multigrid method combined with a Taylor series interpolation scheme as the best discretization acceleration scheme after comparing the use of simple multigrid and multilevel multigrid methods. Based on the previous method, Phillips Phillips and Schmidt (1985) used multigrid combined with multilevel acceleration technology to realize the accelerated solution of scalar conservation equations. In addition, they proposed a fast finite difference solution to the passive scalar advection-diffusion equation. Although these acceleration methods reduced the computational cost while maintaining high accuracy, high computational cost remained a problem due to the need to retain high solution grid modeling in some complex fluid regions.

In recent years, machine learning has been used in the numerical solution of partial differential equations, which have made enormous progress. The combination of machine learning and traditional discretization methods improved the accuracy of the solution and accelerated the numerical calculation Raissi et al. (2019); Ji et al. (2021); Vinuesa and Brunton (2021); Patel et al. (2021); Eliasof et al. (2021); Cai et al. (2022). Based on these methods, O. Obiols-Sales et al. Obiols-Sales et al. (2020) proposed a coupled deep learning and physics simulation framework (CFDNet) to accelerate the convergence of Reynolds-averaged Navier–Stokes simulations. CFDNet was designed to use a single convolutional neural network at its core to predict the main physical properties of fluids, including velocity, pressure, and eddy viscosity. In this paper, CFDNet was evaluated for various use cases, and the results showed that CFDNet significantly speeded up the numerical solution and proved that CFDNet generalized well. Vadyala Shashank Reddy et al. Vadyala et al. (2022) determined the numerical solution of the one-dimensional advection equation using different finite-difference approximations and physical informatic neural networks (PINNs). They trained a neural network to solve supervised learning tasks that obeyed any given laws of physics described by general non-linear partial differential equations. The PINNs approximation was compared with other schemes through experiments, and the results showed that the prediction results obtained by the PINNs approximation were the most accurate. Pathak et al. Pathak et al. (2020) proposed a hybrid ML-PDE solver that combined machine learning and traditional solving methods of the partial differential equation. It can obtain meaningful high-resolution solution trajectories while solving system PDEs at lower resolutions. The ML part of the solver extracted spatial features by using u-net as the model structure to predict the error accumulated in the short time interval between the evolution of the coarse grid and the solution of the system at a higher resolution. The predicted error can optimize the solution generated by the coarse grid to obtain a solution close to that generated by the fine grid, enabling high-precision solutions at low accuracy. Y. Bar-Sinai Bar-Sinai et al. (2019) designed a data-driven discretization scheme using a deep-learning algorithm. They used neural networks to estimate spatial derivatives that were optimized end-to-end to best satisfy equations on low-resolution grids. The resulting numerical method was very accurate, eventually achieving the same computational accuracy as the standard finite difference method at 4 to 8 times coarser resolution than the standard finite difference method. Zhuang [38] improved the model structure and loss function based on Y. Bar-Sinai and applied it to passive scalar advection in a two-dimensional unsteady flow. They used a convolutional neural network to learn spatial discretization coefficients to calculate spatial derivatives. Then, they combined them with traditional numerical methods to calculate time derivatives to obtain the numerical solution of partial differential equations. This method achieved a high-precision solution with a low computational cost. Ranade et al., 2021 developed DiscretizationNet, a machine learning-based PDE solver that combined essential features of existing PDE solvers with ML techniques. They used a discretization-based scheme to approximate spatiotemporal partial derivatives and a CNN-based generative encoder-decoder model with PDE variables as input and output features for iteratively generating equation solutions. Although these methods addressed the problem of traditional methods, their solution accuracy was limited due to the problems of ignoring spatiotemporal characteristics and input redundancy.

If the velocity field is divergence-free, the advective form of the scalar concentration field C ( x → , t ) for a given velocity field u → ( x → , t ) is as follows Zhuang et al. (2021):

The objective of the numerical solution for the passive scalar advection in 2-D unsteady flow is to predict the concentration field distribution at each time step in the future under the influence of the randomly changing velocity field given the initial concentration field. In this paper, we predict the concentration field distribution results in the 32 time steps to demonstrate the ability of our model to make multi-step predictions. We employ a rolling forecasting scheme in which we input multiple velocity fields between t ₀ and t ₁ into the prediction model and combine the concentration field distribution at t ₀ to predict the concentration field distribution at t ₁ .Then, we input multiple velocity fields between t ₁ and t ₂ into the model and combine the concentration field distribution at t ₁, predicted by our model to predict the concentration field distribution at t ₂. According to this calculation rule, we use multiple velocity fields between t_n and t_{n+ 1} and the concentration field distribution predicted at time t_n to predict the concentration field distribution at t_{n+ 1}.By repeating this process, we can get the passive scalar advection solution at each time in the future. Therefore, the key to our multi-step prediction method is to recursively predict the concentration field distribution at a single step, i.e., the numerical solution of passive scalar advection at the next time step. We propose the time-sequence-involved space discretization neural network (TSI-SD) to predict the space discretization coefficient for the space derivative and then combine the finite volume method to calculate the numerical solution of the next time step.

The framework of the proposed method is shown in Figure 2 . This is a fusion framework of deep learning (TSI-SD) and a traditional numerical method (FVM) for end-to-end numerical solutions of passive scalar advection equations. It consists of three modules: the spatial discretization coefficient prediction module (SDCPM), the concentration template extraction module (CTEM), and the concentration solver module based on finite volume numerical format (CSM). For the set of multiple velocity fields between the time steps t_n and t_{n+ 1}, wedecompose each velocity field into two sub-velocity fields in the horizontal and vertical directions (along the x-axis and y-axis) to obtain the velocity field set in the two directions. In the next step, we build the time-sequence-involved space discretization neural network (TSI-SD) in the SDCPM. TSI-SD extracts the spatiotemporal features from the decomposed velocity field sets in the two directions separately and then fuses them to obtain the spatial discretization coefficient of each grid point. After that, we input the coefficients into the CSM and combine them with the surrounding point concentration template of each grid point obtained by the CTEM to calculate the spatial derivative. Finally, we could calculate the concentration of each grid point at the next moment t_{n+ 1},, that is, the concentration field of t_{n+ 1} by the FVM in the CSM.

The equation-solving process can be roughly described in the following three steps:

Next, we will provide details of our proposed framework for end-to-end numerical solutions of passive scalar advection equations. First, we introduce the SDCPM and the TSI-SD in the section entitled ‘Spatial Discretization Coefficient Prediction Module’. Then, we describe in detail our proposed CTEM and CSM modules in the Concentration Template Extraction Module’ and ‘Concentration Solver Module Based on Finite Volume Numerical Format’ sections, respectively. Finally, we discuss the loss function in the ‘Loss Function’ section 3.6

In this module, we design the time-sequence-involved space discretization neural network to predict the spatial discretization coefficients, and the prediction function is

where U s the set of multiple two-dimensional velocity fields between t nd t+1, of which size is nW s the weight of our neural network. The time interval between the velocity fields is 1 n

We decompose U into velocity field groups U_x in the horizontal direction (along the x-axis) and U_y in the vertical direction (along the y-axis),

Then, we extract the spatiotemporal features separately for the decomposed velocity field sets in the two directions. Taking U_x as an example, we input the velocity field at each time step from the velocity field set { u t x , u t + 1 n x , … , u t + 1 − 1 n x } as the spatial feature of each time step into the different conv-lstm structural unit,

and S_k is regarded as the spatial feature at time t + 1 k . A CONV-LSTM structural unit contains convolution operations and long-short-term memory unit processing operations. The calculation steps can be written in the following form,

where i_k is the input gate, which is used to calculate how much information of the current state to retain. f_x is the forget gate, and its function is to calculate how much of the output information of the previous moment is discarded. g_x is theinformation extracted from the current state. C_{k -1} is the information of the previous moment. C_k is the final state at the current moment, calculated by f_k , c_{k -1}, g_k , and i_k . o_ki the output gate, which is used to calculate how much information needs to be output (to the cell at the next moment). h_k is the final output information of the state, which is calculated by o_k and c_k . w _xi ,w _hi ,w _xf ,w _hf ,w _xg ,w _hg ,w _xg ,w _ho ,b _i ,b _f ,b _g ,and b _o are the weights designed in our neural network, and these weights will be updated during the model training process.

After the information processing and transmission of n conv-lstm structural units, the information h t + 1 − 1 n output by the last unit is obtained. The final spatiotemporal fusion information I_x of the horizontal velocity field is calculated by using the output information.

In the same way, we obtain the final spatiotemporal fusion information I _y f the vertical velocity field.

After obtaining the spatiotemporal fusion information I _x and I _y in two directions, it is necessary to re-fuse the spatiotemporal features in the horizontal and vertical directions on I _x and I _y concat(), is a feature merging operation that integrates two features in a new dimension. After the feature merging operation, convolution is performed on the merged features to process the spatial information of the merged spatiotemporal features. Finally, the spatial discretization coefficient matrix α is obtained.

The dimension of the α matrix is (s,s,template_size*2) , where s is the side length of the input two-dimensional velocity field, and (s, s) is the dimension of the two-dimensional velocity field. template_size is the number of weights required for each grid point. We divide α into the grid upper boundary space discretization coefficient α_up and the grid right boundary space discretization coefficient α_right with dimensions (s, s, template_size).

This module and the next module follow the numerical solution part of the traditional advection equation adopted by Zhuang et al. Zhuang et al. (2021), and adopt the spatial derivative of the classical Euler algorithm.

In the previous part, we calculated the spatial discretization coefficient templates α_up and α_right , the dimensions of which are (s, s, template_size). Therefore, we need to find the surrounding grid point concentration templates C_up and C_right corresponding to the position of the coefficient template, the dimensions of which are both (s, s, template_size), which indicates that the number of surrounding grid point concentrations required for each point in the two-dimensional space field is template_size. As shown in Figure 3 , we input the two-dimensional concentration field C_t at time t, and its dimension is (s, s). We model the upper and right boundaries of each point in the two-dimensional matrix and obtain the concentration values of m*n grid points around it as the grid point concentration template, where

m and n are the length and width of the two-dimensional grid point concentration template. Finally, we obtain C_up and C_right with dimensions (s, s, template_size).

In this module, we first calculate the upper boundary concentration C _{up_edge} and the right boundary concentration C_{up_edge}

SUM() is the defined summation of the last dimension of the matrix, i.e., after the matrix of (s, s, template_size) is obtained through the dot product operation, the last dimension is summed to obtain the boundary concentration C_edge with dimension size (s, s) .The lower boundary concentration C_{lower_edge} and the left boundary concentration C_{left_edge} of the grid point can be directly obtained from the upper boundary concentration of the adjacent grid below its position and the right boundary concentration of the adjacent grid to the left of its position. Then, we can obtain the boundary velocity u_edge bythe same method as the calculation of the concentration boundary and boundary flux via C_edge u_edge. After obtaining the flux at the four boundaries of the grid, the traditional finite volume method is used to calculate the time derivative to obtain the concentration field distribution at the next time step, as shown in Figure 2 .

The format of the mean absolute error (MAE) used to train our model is as follows:

where C_{t +1} is the concentration field at time t+1 redicted by our model, and C ^ t + 1 is the high-precision numerical solution at 16*16 low resolution grids. The numerical solution is calculated using 128*128 high resolution grids by the second-order Vanleer format and then transformed to the solution at 16*16 low resolution grids by the dimensionality reduction method Zhuang et al. (2021).

We used the theory of divergence-free velocity field described by Saad and Sutherland (2016) to generate a divergence-free random velocity field set with the resolution of 128*128. Then, the set was divided into two parts of divergence-free random velocity field sets: the training part and the test part. These two parts were completely different to ensure the generalization of the model.

For the training set, we generated a variety of random initial concentration fields and used the second-order Vanleer numerical format to calculate the numerical solution of the equation, i.e., the concentration fields at multiple time steps with the resolution of 128*128 based on the set of divergence-free random velocity fields in the training part. It is worth noting that for the C_{t +1)} to be generated, our model needs to input the velocity field i at time t , t + 1 n , … , t + 1 − 1 n . Therefore, we set the time step length of the velocity field to be smaller than the concentration field in the generation process to ensure that the velocity field in the time interval from t to t+1 could be generated. Next, we sampled both the velocity field and the concentration field at intervals to obtain a high-precision velocity field and concentration field with a resolution of 16*16 using the dimensionality reduction method Zhuang et al. (2021). Each training sample included an input part and an output part. The input part was the velocity field and the concentration field C_t at time t at multiple time steps in the time interval from t to t+1,and the output part was the concentration field C_{t +1}. The test set generation process was consistent with the training set, but it was necessary to ensure that the random initial concentration field generated in the test set was different from the training set.

The initial and boundary conditions for the velocity and concentration fields were set as follows. The size of the two-dimensional velocity field and the two-dimensional concentration field were both [0,1] × [0,1]. Our velocity field was a divergence-free random velocity field, and the magnitude of the velocity was limited between -1 and 1. The concentration field used periodic boundary conditions, and its initial condition is to

set the concentration value range between 0 and 1. The calculation process is shown in formulas (23)-(27).

The C(x, y) is as shown in Figure 3 . C represents the concentration value.

In this part, four SOTA numerical solution methods for passive scalar advection in a two-dimensional unsteady flow were selected as our baseline: (1) traditional solvers based on 16*16 resolution grids using the second-order Vanleer discretization format(Vanleer 16*16) Lin et al. (1994); (2) traditional solvers based on 32*32 resolution grids using the second-order Vanleer discretization format (Vanleer 32*32) Lin et al. (1994); (3) traditional solvers based on 64*64 resolution grids using the second-order Vanleer discretization format (Vanleer 64*64) Lin et al. (1994); and (4) a hybrid solver based on a CNN and the finite volume method (CNN+FVM) Zhuang et al. (2021).

We first compared our TSI-SD method with traditional solvers, in which TSI-SD uses a 16*16 low-resolution grid. As shown in Figure 4 , the TSI-SD method maintained the smallest prediction error over 32 time steps, which demonstrates that our method achieved a higher solution accuracy than the traditional method at a resolution of 4× lower than the traditional method.

Then, we compared TSI-SD with the CNN-FVM solver trained based on the previous spatial discretization scheme Zhuang et al. (2021). The CNN-FVM method is currently one of the most outstanding methods for solving partial differential equations in deep learning. It has been proven to achieve very good prediction and solution results in various partial differential equations, such as Burgers’ equation Bar-Sinai et al. (2019), and the advection equation Zhuang et al. (2021), Additionally, the method has been proven effective at solving complex Navier–Stokes equations Kochkov et al. (2021),and results are as accurate as baseline solvers, with 8–10× finer resolution in each spatial dimension, resulting in 40- to 80-fold computational speedups. The original CNN-FVM solver has a prediction error of 0.0043, a single-step prediction time of 0.2712s, and a single-sample training time of 4 ms per round during the training process. Our single-step solver had an error of 0.0029, a single-step prediction time of 0.2474s, and a single-sample training time of 2ms per round during training. Our single-step error was 32.56% lower than the previous method, and the iterative prediction error after 32 steps was greatly reduced. As shown in Figure 5 , our method also outperformed the CNN-FVM solver in continuous prediction results within 32 time steps.

The reason why our solver outperformed the CNN-FVM solver in training time, prediction time, and prediction accuracy is as follows. In the spatial discretization coefficent prediction part, the inputs of the CNN solver’s prediction deep-learning model are the concentration field with (batch_size,1,grid,_size,grid_size) and the two velocity field (along the x-axis and y-axis) at a time step with (batch_size,2,grid,_size,grid_size). The input to our TSI-SD was the horizontal velocity fields along the x-axis at two time steps with (batch_size,2,grid,_size,grid_size) and the vertical velocity fields along the y-axis at two time steps with (batch_size,2,grid,_size,grid_size), so our input size was larger than the previous input size. However, in the model part, the CNN-FVM solver used a five-layer convolutional neural network to process the data collected by the concentration field and the velocity field with (batch_size,3,grid,_size,grid_size) .We used the structure of a 1-layer

convolutional neural network to process the horizontal and vertical velocity fields respectively, and then a one-layer convolutional neural network was used to process the integrated features. After inference analysis, our model parameters were fewer than the original model parameters, which resulted in a shorter training time and prediction time in our model compared with the original model training time. This was also confirmed by a saved parameter file size comparison.

Finally, we demonstrated the evolution prediction effect of an initial concentration field under different models after 32 iterations. As shown in Figure 6 , the third row shows the prediction effect of our model. The first row is our high-precision numerical solution generated using a 128*128 high-resolution grid. The second row shows how we use the averaging operation to obtain a high-precision numerical solution at a low resolution of 16*16, which is used as the ground truth of our model. The fourth and fifth rows are the results obtained using the second-order Vanleer 16*16 and CNN-FVM solvers. Figure 5 shows that our model is better than the CNN-FVM and traditional second-order Vanleer 16*16 solvers. In Figure 6 , C represents the concentration value.

In this part, we used different sets of time steps as the time series information input to TSI-SD, so that our model could extract different time features to predict the spatial discretization coefficient. The best prediction result represents the velocity fields at the selected time steps that have the greatest influence on the coefficients. Figure 7 shows that when the set of fine velocity fieldd { u t , u t + 1 n , … , u t + 1 − 1 n } was selected to replace velocity field set {u _t−n+1,…,u _t−1,u _t } to predict u_{t +1} could reduce the prediction error of the model. The experimental result demonstrates that the set of fine velocity fields extracts spatiotemporal features more effectively. That is because the time interval of the velocity field set we chose was close to the time of the predicted concentration field, so the correlation between the velocity field set and the predicted concentration field was strong. The model could learn the spatiotemporal influence of the velocity field on the concentration field from this set of velocity fields, which could accurately predict the spatial discretization coefficient. Meanwhile, the prediction error of the velocity field using { u t , u t + 1 2 } is the best, and experiments demonstrated that it involves lower computational cost; therefore, so we finally choose the

velocity field of { u t , u t + 1 2 } as the velocity field input of our final model. We think that for the 16*16 lower resolution grid, the model learned the time-space correlation between the velocity field set and the concentration field well through the analysis of the velocity fields at two times through a large amount of training data, which is also consistent with the experimental results as shown. In future studies, we will conduct more experiments on higher-resolution grids to obtain the optimal number of time steps after increasing computing power.

In this section, we carried out an experiment to prove the excellent performance of our model under a constant velocity field and a two-dimensional deforming flow velocity field. We generated the concentration under a constant velocity field, and the two-dimensional deformation flow concentration field under the velocity field:

The predicted performance is shown in Figure 8 and Figure 9 . C represents the concentration value. Our model achieved outstanding prediction results in the iterations of 32 time steps. However, at the same time,

there are also the following problems: even under a simple constant velocity field, the prediction effect will become worse and worse with the long-term iteration due to the accumulation of errors predicted by the model at each time step. We will try to fix this in the future.

In this part, we verified the advantage of only taking the velocity field as the input feature on our model. A contrast model that adds the concentration field as feature input was designed to prove our inference. The contrast model was identical to ours except that the concentration field features were fused with the spatiotemporal features extracted from the horizontal and vertical velocity fields (along the x-axis and y-axis) in the fusion module. Figure 10 shows that the iteration errors on 32 time steps of our model are lower than those of the contrast model. Therefore, we proved that the input of the concentration field information was redundant and verified our conclusion: the spatial discretization coefficients are strongly correlated with the velocity field at multiple time steps before, while the concentration field information becomes redundant when predicting the coefficients. In other words, the change in the velocity field is the main factor for the change in the concentration field. Our model extracts effective spatiotemporal features from the velocity field set to learn the influence of the change of the velocity field set on the change of the concentration field, which is very helpful for predicting the spatial discretization coefficient.

In this part, we proved that the spatial discretization coefficients predicted by our model have upwind properties on a constant velocity field. A two-dimensional velocity field U ₁ with a horizontal velocity field (along the x-axis) of +1 and a vertical velocity field (along the y-axis) of +1, and a two-dimensional velocity field U ₂ with a horizontal velocity field of -1 and a vertical velocity field of -1, were designed to prove our model’s upwind properties on a constant velocity field. Under the two velocity fields, the visualization process of the concentration coefficients of the upper and right boundaries of grid points A and B was completed.

As shown in Figure 11 , C represents the concentration value and Coefficient represents the coefficient value. For the upper boundary, the concentration on the right boundary of the constant velocity field is mainly determined by the concentration of the two adjacent grids. When the horizontal speed is +1 (i.e., the direction is to the right), the grid coefficient on the left of the right boundary of grid A is greater than the grid coefficient on the right; when the horizontal speed is −1 (i.e., the direction is to the left), the grid coefficient on the left of the right boundary of grid A is smaller than the grid coefficient on the right. For the right boundary, the concentration on the upper boundary of the constant velocity field is also mainly determined by the concentration of the two adjacent grids. When the vertical speed is +1 (i.e., the direction is downward), the grid coefficient above the lower boundary of grid A is greater than the grid coefficient below; when the horizontal speed is −1 (i.e., the direction is upward), the grid coefficient above the lower boundary of grid A is smaller than the grid coefficient below.

The concentration coefficient of another spatial grid point B is almost the same as that exhibited by A. Therefore, our grid coefficient has nothing to do with the distribution of the concentration field, but only with the distribution of the velocity field. The concentration field distributions at point A and point B are completely inconsistent, but under the same velocity field, the predicted spatial discretization coefficient distributions are basically the same, which proves that there is no significant correlation between the concentration field distribution and the spatial discretization coefficient.