{'fraction_non_alphanumeric': 0.03196275520373587, 'fraction_numerical': 0.012400694780603093, 'mean_word_length': 4.598517455762793, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 99, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Performance and accuracy assessments of an incompressible fluid solver coupled with a deep Convolutional Neural Network

The resolution of the Poisson equation is usually one of the most computationally intensive steps for incompressible fluid solvers. Lately, Deep Learning, and especially Convolutional Neural Networks (CNN), has been introduced to solve this equation, leading to significant inference time reduction at the cost of a lack of guarantee on the accuracy of the solution. This drawback might lead to inaccuracies and potentially unstable simulations. It also makes impossible a fair assessment of the CNN speedup, for instance, when changing the network architecture, since evaluated at different error levels. To circumvent this issue, a hybrid strategy is developed, which couples a CNN with a traditional iterative solver to ensure a user-defined accuracy level. The CNN hybrid method is tested on two flow cases, consisting of a variable-density plume with and without obstacles, demostrating remarkable generalization capabilities, ensuring both the accuracy and stability of the simulations. The error distribution of the predictions using several network architectures is further investigated. Results show that the threshold of the hybrid strategy defined as the mean divergence of the velocity field is ensuring a consistent physical behavior of the CNN-based hybrid computational strategy. This strategy allows a systematic evaluation of the CNN performance at the same accuracy level for various network architectures. In particular, the importance of incorporating multiple scales in the network architecture is demonstrated, since improving both the accuracy and the inference performance compared with feedforward CNN architectures, as these networks can provide solutions 1 10-25 faster than traditional iterative solvers.

to the numerical setup: a small error on the gradient estimation may lead to very large differences in the solution. As a result, the CFD community has produced an extreme effort over the last decades to improve numerical schemes [1] towards high-order discretization with high accuracy, most often at a price of a very high CPU cost. Such simulations have been made possible by massively parallel computing (High Performance Computing, HPC) combined with a constant amelioration of the computational hardware. More recently, the development and improvement of Graphical Processing Units (GPUs) has also enabled to accelerate calculation codes, and has opened the path to Machine Learning (ML) techniques that were not yet developed in the CFD community.
The theoretical foundations of ML models were developed in the mid-XX-th century. In 1958 Rosenblatt [2] developed an Artificial Neural Network (Perceptron), while the backpropagation was introduced in 1960 by Kelley [3] for control applications (although the first successfully application was obtained by Rumelhart et al. [4] in 1986). Convolutional Neural Networks (CNNs) appeared slightly later, when LeCun et al. [5] created a network that extracted spatial features from 2D images. The performance of CNNs on image classification problems [6] and the development of more efficient GPUs then sky-rocketed the use of ML in various domains. Mostly driven by the computer vision community, open-source frameworks were developed, which led to an even further spread of Neural Networks (NNs). This recent data-driven wave finally touched the CFD community, who made substantial efforts in the past 5 years to apply ML to flow problems.
ML techniques are optimization algorithms that extract patterns and structures from data, "generalizing" the valuable information found on its training data. When applied to Fluid dynamics, this usually translates into the creation of Reduced-order Models (ROM), a simple representation of the information contained in the dataset. These simpler representations can be used in two different ways in CFD, as reviewed by Brunton et al. [7]: either discovering unknown models, or improving existing models. Discovering unknown models from data has been a long-standing problem in the scientific community. Throughout the history of humanity, people built models to describe the behavior of observed phenomena, in order to make future predictions. Nowadays, with the advance of technology, a time has reached where a lot of data of different processes are available, yet being still under-exploited. This has promoted Machine Learning, with the aim to benefit from these data, by extracting features to develop reduced models. For instance, Brunton et al. [8], followed by Rudy et al. [9], established governing equations for the chaotic Lorenz system or for the vortex shedding behind a cylinder using a sparse dynamics model and a sparse regression model. Another example is the use of Machine Learning for problems where the governing equations are well known, but too expensive to resolve, as in turbulent flows, for example Zhang et al. [10] used neural networks to predict subgrid-scale correction factors on turbulent flows. Fukami et al. [11] applied CNNs to reconstruct turbulent flows that respect the energy spectrum, while other works [12,13,14] used Generative Adversarial Networks (GANs) to reconstruct turbulent fields. Similarly, CNNs have also been used in turbulent combustion to predict the sub-grid flame wrinkling [15].
When the resolution of PDEs is computationally expensive, Machine Learning techniques can be used to develop surrogate models that capture the essence of the studied phenomena with a reduced computational cost. Classical reduced-order models, [16,17,18] include proper orthogonal decomposition (POD) or dynamic mode decomposition (DMD). These methods can be combined with Machine Learning techniques, such as autoencoders [19] or long shortterm memory (LSTM) neural networks [20] in order to get satisfactory results on problems like the 2D Navier-Stokes equations. The first studies introducing ML to solve PDEs date back to 1990, with the work of Lee et al. [21] followed by Dissanayake et al. [22] in 1994. Dissanayake's team proposed a multilayer feedforward network to solve a linear Poisson equation and a thermal conduction problem on a non-linear heat generation test case. Neural Networks can be described as 'universal approximator' functions [23], which enable to easily differentiate their inputs with respect to their output. That way, they constituted a novel approach where the classical finite-element or finite-difference discretization problem was substituted by an unconstrained optimization problem. Following the reduced-order-modeling philosophy, their objective was to obtain an easy-to-implement solver that could rapidly give accurate results. Lagaris et al. [24] further exploited the idea, adding both Neumann and Dirichlet boundary conditions. Gonzalez-Garcia et al. [25] developed a different approach, where they modeled the right-hand side of a PDE with a neural network, in order to turn the problem into a more simple Ordinary Differential Equation. They showed that the network enables to generalize from experimental data which might be incomplete, or noisy, testing on the Kuramoto-Shivanisky equation. All these works profit from the differentiable nature of neural networks. However, these networks were rapidly limited by the available computational capability and were restricted to a single particular PDE, with not much margin for generalization.

Configuration and methodology 2.1 Plume Case
Buoyancy forces are induced by density variations in the presence of gravity. The study of buoyancy-driven flows has been a long-standing problem for the fluid mechanics community [40], as they correspond to a wide variety of remarkable geophysical flows, from volcanic eruptions to avalanches. The Earth's climate and ocean's behavior are strongly dependent of buoyancy, as the ocean's currents and tides depend on the creation and distribution of water density variations. Buoyant coastal currents for example, are responsible of redistributing the fresh water coming from rivers and other sources into the ocean, carrying as well sediments or pollutants which can strongly impact natural ecosystems. A thorough analysis of buoyancy-driven flows, especially focused on oceanic flows, can be found in the work of Chassignet et al. [41]. The civil engineering community also studies flows driven by buoyancy forces, which control the ventilation of buildings [42] or even entire cities, with direct effect on the propagation of pollutant particles [43]. The propagation of fires in buildings has also led to extended work with theoretical, numerical and experimental models for smoke propagation, as shown by Ahn et al. [44].
Usually the term buoyancy-driven flow refers to the situation where a lighter moving fluid is surrounded by a heavier fluid, whereas a heavier moving fluid surrounded by a lighter fluid is called a density-driven flow. This work focuses on plumes, which result from injecting a light fluid into a heavier quiescent environment. Plumes can be found in a wide variety of physical phenomena, such as the flow escaping the chimney of a locomotive steam engine or even the smoke elevating from a cigarette. These types of flows have been studied for over 80 years, after Zeldovich's [45] work on similarity analysis, and Schmidt's [46], who obtained analytical expressions for the mean velocity and temperature profile in turbulent plumes. Despite the presence of a wide variety of experimental [47,48] and theoretical [49] works in this domain, the phenomenology underlying plume physics is still not completely mastered. Numerous works have used CFD to understand this type of flows, based either on Reynolds Averaged Navier-Stokes equations [50,51] or LES [52]. Even for simplified incompressible and inviscid plumes, the interaction between the gravity source term and the momentum-driven forces results in an interesting and complex behavior. The misalignment of the density and pressure gradients causes the generation of vorticity and the development of instabilities [40]. Thus, even with a simplified setup, a variety of phenomena with different spatial scales can be found. For further information, the reader can refer to the reviews of Turner [53] and Hunt et al. [54].
The plume case has also been widely treated in the computer vision community, due to the difficulty of calculating realistic smokes in real-time, which requires efficient computational solvers, and explains why Machine Learning was first introduced in this context. One can cite the works of Chu et al. [55] with a CNN, Kim et al. [56] using GAN with a novel stream function-based loss and Tompson et al. [37] who introduced the idea of a CNN coupled to an incompressible fluid solver.

3
The two configurations studied in this work are displayed in Fig. 1. The flow is assumed inviscid (infinite Reynolds number) and adiabatic. Thus, the plume behavior is driven by only one dimensionless number, called the Richardson number (Eq. 1), which compares the buoyancy and inertia forces: where Δ / 0 is the density contrast between the light and heavy fluid, g is the gravity, a characteristic length, and the velocity of the injected lighter fluid. High Richardson numbers imply that the flow is piloted by the buoyant forces, whereas low Richardson numbers indicate that the flow is driven by momentum (also known as jets). In both cases, the domain is a square box of length , with an inflow centered in the bottom face of the domain, of radius = 0.29 . In the second configuration, a cylinder of diameter = 0.312 is added, aligned with the horizontal axis and centered at 0.625 on the vertical axis (Fig. 1). The objective is to evaluate the behavior of the proposed learning-based approach when the flow impinges objects and walls. Figure 1: Plume configuration with and without cylinder

Numerical setup
For the CFD study, a fluid solver based on the open-source solver Mantaflow [57] is used to solve the Euler equations. A variable-density incompressible formulation is employed with external forces due to gravity. To model buoyancy-driven flows under the incompressibility assumption, the Boussinesq approximation [58] is applied, which holds for small density variations compared to the mean density, i.e., Δ / 0 1. Note that the continuity equation with the Boussinesq approximation and the incompressibility constraint leads to the advection of the density variation Δ as a passive scalar, i.e., being constant along a streamline.
The Euler equations are solved in two steps: (i) the advection step and (ii) the pressure projection step. The computational domain is discretized on a uniform cartesian grid with a finite difference scheme on a Marker-and-Cell (MAC) staggered grid [59]. First, the density is convected along streamlines, using an unconditionally stable, second-order accurate Maccormack semi-Lagrangian method [60]. Once the density field is advected, the momentum equation is solved through a splitting method, where the velocity field is first advected similarly to the density field, and then modified by the external forces with a forward Euler scheme, resulting in the velocity field * . Since the pressure term in the momentum equation and the mass conservation equations have not yet been used, the velocity field * does not satisfy the incompressible kinematic constraint ∇ · * = 0. To do so, this pressure term is used as a correction term to ensure a divergence-free velocity field, so that the corrected velocity reads: Taking the divergence of both sides of Eq. (2) and imposing the continuity equation ∇ · u = 0 yields the well-known Poisson equation: This equation is solved to obtain the pressure correction, ensuring that the velocity field at the next time-step satisfies both the momentum and the mass conservation equations. This correction is applied using Eq. (2).
The Poisson equation is a second-order elliptical partial differential equation, which is well known to be computationally expensive, taking up to 80% of the computational time of the incompressible solver [39]. For this reason, this work only focuses on the resolution of this particular equation. This equation may be solved using direct methods [61,62], such as the LU factorization, which can be interpreted as a Gaussian Elimination Method (GEM). The idea behind these methods is to recast the A matrix of the linear system = into an upper triangular and a lower triangular matrix = , which enables to decompose the problem into two successive GEMs, = and = respectively. These methods can provide solutions at the machine-precision, however, they tend to be computationally expensive, especially for large computational domains since the cost of solving a single GEM is proportional to 2 3 /3 operations, where is the number of points in the mesh. Iterative methods have been proposed to avoid the direct solution of the system. For such methods, a trade-off has to be chosen between accuracy and computational cost. Despite the need of multiple iterations to achieve convergence, those methods are usually more efficient than direct approaches since their computational cost is proportional to 2 per iteration. Among the iterative methods, the Jacobi method, which was first developed in 1845 [63], by the German mathematician Carl Gustav Jacobi, decomposes the matrix A into a diagonal matrix D and the upper and lower diagonal matrices L and U. The linear system can be expressed as = − ( + ) , which is equivalent to = −1 ( − ( + ) ). An iterative process is obtained by solving the previous equation in the form +1 = −1 ( − ( + ) . As any iterative method, the number of iterations required towards convergence strongly depends on the initial guess 0 : an additional result of the present work compared with previous studies is to use the neural network prediction as initial guess to speed-up the Jacobi method. More advanced methods [64,65] exist for the Poisson equation, yet this work will focus only on the Jacobi method for the sake of stability, its easy implementation on GPUs, as well as its straightforward hybridisation with neural networks.

A deep-learning approach to solve the Poisson equation
To overcome the difficulties of iterative methods, the present study focuses on solving the Poisson equation with a Convolutional Neural Network (CNN). Particular attention is put on its ability to tackle new cases, not seen during the training phase. To further ensure both a fast and reliable prediction of the corrected velocity field, this CNN method will be hybridized with the Jacobi method, using the CNN prediction as initial guess to significantly speed-up the Jacobi approach.
Convolutional Neural Networks perform multiple linear transformations, typically = + , followed by non-linear activation functions to output the desired field. The training of the neural network consists in the optimization of the weights and biases to minimize a user-defined objective function L. CNNs differ from traditional MLP as they do not apply the weights and biases to the entire input fields. Indeed, CNNs apply the transformation locally through a kernel of limited size, typically 3 × 3, which is then 'slid' through the entire input field. This type of network is widely used for image recognition tasks, as the kernels can be interpreted as filters that extract patterns and structures from the input images. CNNs are therefore particularly adapted to fluid mechanics, since the underlying physics are mostly driven by flow structures (jets, vortices, stagnation points etc.). Additionally, compared with MLP, CNNs usually require a smaller number of tunable parameters ( and ), and are theoretically independent of the size of the input field: the same trained network can be reused on cases with different resolutions or meshes. However, a current limitation of CNNs is the need of perfectly uniform cartesian grids, whereas most of today CFD simulations are performed on non-uniform, possibly unstructured, meshes. While interpolation from non-uniform meshes onto a uniform cartesian grid is possible, recent works have efficiently adapted convolution operations to unstructured meshes [66], making CNNs even more suitable for CFD applications. However, these methods are still complex and not fully mature, so CNNs are still widely used on CFD. For the sake of simplicity, this study considers a perfectly uniform cartesian grid. The spatial resolution will however be varied to validate the convergence properties of the method, as well as to evaluate its computational performance on different grid sizes.

Network architectures
The choice of the neural network architecture is usually driven by the user's experience, and further improved by a trial and error approach. Since the network architecture is critical, tools for automatic architecture search have been also developed [67]. While discovering efficient network architectures, these automated data-driven search strategies require an extremely high computational time since multiple trainings have to be performed, as shown by Zoph et al. [68] who used in parallel 500 Nvidia p100 GPUs during 4 days to find an optimal network for scalable image classification. Architectures are even more critical in physics, and especially fluid mechanics, since the physics usually involve a large range of different scales which need to be captured accurately by the network. Towards this objective, Geng et al. [69] have employed such method to discover optimal neural networks dedicated to multi-scale analysis of seismic waves. However, despite those automatic strategies, the choice of the network architecture and its effect on accuracy and inference time is often ignored. In particular, guidelines still need to be established for physics-related tasks. Here, different types of network architectures are investigated, in order to analyze their effect on the resolution of the Poisson equation. Note that the objective of this work is not to find the optimum architecture for the resolution of the Poisson equation, but rather to provide general trends and understanding of the architecture effect on the network accuracy and performance. For all architectures tested here, the pressure field is computed using a CNN for which the inputs are the uncorrected velocity divergence ∇ · * and a boolean field describing the object geometry.
First, a simple "feed-forward" convolutional network is considered, denoted hereafter MonoScale. This network (Fig. 2) contains 395601 tunable parameters distributed in 9 layers, following a straightforward architecture. For convolution operations, a replication padding is used to ensure that all the feature maps keep the size of the original fields. This type of padding creates ghost cells on the image boundaries with a value copied from the closest image pixel. This network is the most simple strategy to tackle the Poisson problem, as the input is just passed through a series of convolutional layers, without any further modification or post-processing. * * * * * * The second architecture is the MultiScale (Fig. 3), introduced by Mathieu et al. [70] for video image prediction, and introduced to solve the Poisson equation by Ajuria et al. [39]. It has also been employed in other studies on fluid mechanics, such as in Fukami et al. [71,72] focusing on the super-resolution of turbulent flows, or applied to the propagation of acoustic waves by Alguacil et al. [73]. The idea behind this architecture is to feed the neural network with several representations of the same inputs, focusing on different scales. To do so, the original input is interpolated on coarser meshes: while the largest resolution retains all scales, smaller resolutions focus only on the large scales. The network then encodes information linked to the different spatial scales at each level, which helps the network on the generalization task, and avoids spurious high-frequency oscillations on the large features of the image [70]. Here the MultiScale architecture (Fig. 3) has three scales, of sizes 2 1/4 , 2 1/2 and 2 respectively. The first scale interpolates the original images on a quarter-cell size mesh. The resulting image is then interpolated on a half-cell size mesh to twice its size ( 2 1/2 ). The middle scale takes as an input both the interpolation of the original entries on a half-cell size mesh, and the output of the previous scale. Finally, the resulting image is interpolated to the original size 2 , and is concatenated to the input fields. It is then fed to the final scale branch. The first scale contains 5 layers, while both the intermediate and final scales contain 7 layers. Non-linear ReLU activation functions are placed after each convolution layer, except on the last two ones, to enable the network to compute positive and negative outputs.
The MonoScale and MultiScale networks are then compared with the well-known Unet architecture (Fig. 4), which was first developed by Ronneberger et al. [74] for the segmentation of biomedical images. Apart from the initial biology-related applications, this network has also been used in many other fields, including for regression tasks in fluid mechanics. For instance, Lapeyre et al. [15] employed this architecture on turbulent sub-grid scale modeling for reacting flows. The structure of the Unet is similar to the MultiScale, as it also combines information extracted Figure 3: MultiScale network, with 418640 parameters, where • corresponds to convolution with kernels of size 5×5, * to kernels of size 3×3 and to kernels of size 1×1. R corresponds to the ReLU activation function, indicates a bilinear downsampling operation, whereas corresponds to the bilinear interpolation. Each box indicates the number of feature maps, or channels ( ) present at each layer.
at different spatial scales. However, the Unet differs from MultiScale on the scale treatment: while several scales are managed in parallel in the MultiScale network and the input image is fed directly to each branch, the Unet acts as a simple feed-forward network with a decreasing-increasing resolution to encode information in a low-dimensional manifold (at the smaller scale). The image input is fed to the network only at the initial resolution. Additionally, skip connections are imposed at each scale, allowing information to by-pass the lower-resolution treatments and to avoid the vanishing gradient issue inherent to deep network architectures. *  Note that those architectures are similar to multigrid and multi-fidelity approaches developed in CFD. Indeed, the MultiScale network corresponds to multi-fidelity approaches [75], where the original problem is solved on a coarse grid, the solution is then used as an initial guess to a more refined grid [76]. The Unet on the other hand is similar to multigrid solvers [77], where the initial problem is solved with a computationally non-expensive solver, and the result is corrected with approximations computed on coarser grids.
In comparison with classical deep neural networks employed for classification tasks which may contain 10 6 (e.g. ResNets [78]) to 10 8 (e.g. VGG-16 [79]) tunable weights, the proposed architectures are relatively small, containing all about 4 10 5 parameters. Since the main goal of the present approach is to accelerate the classical iterative Poisson solvers, a limited number of parameters have been chosen. Note however that this cannot constitute a guideline since the performance of the networks will be evaluated at iso-level of accuracy: the number of parameters is, therefore, a non-trivial trade-off between accuracy and inference time, for which no rule exists in the literature. To investigate this accuracy-performance link, a smaller MultiScale containing 1.3 10 5 parameters is introduced in complement to the three previous networks. The structure is exactly the same as the MultiScale, but only the number of filters per layer is changed, resulting in a network with around 3 times fewer parameters. This network will be denoted SmallScale in the following.
Consequently, four deep neural networks will be analyzed on the two test cases introduced in Section 2.1. For each architecture, both their error levels and distributions, as well as their performances on inference time, will be assessed.

Training and Loss function
The networks introduced in Section 3.1 are trained using a procedure similar to Tompson et al. [37] and Ajuria et al. [39]. Compared with most studies on Machine Learning using a known output as target in the loss function L for a supervised training [5], here a semi-supervised learning strategy is used, where no 'ground truth' field is needed. To do so, the residual of the continuity equation is used as loss function, where the divergence of the velocity field is computed using a first-order forward finite difference scheme. Thus, training the network is equivalent to minimizing this residual, i.e., to enforce the mass conservation equation, which is the goal of solving the Poisson equation in an incompressible solver. Note that this approach can be considered as a Physics-Driven Neural Network (PDNN), demonstrating how a target physical equation can be introduced into the learning strategy. It differs from Physics-Informed Neural Networks (PINN, [28]), which combines a semi-supervised physics-driven approach with a standard supervised learning method in which a 'ground truth' target is required. The main benefit of PDNN over PINN, is the possibility of adding long-term loss effects during the training by constraining the network to produce consistent predictions in time. Figure 5: PDNN learning strategy by combining a short-term and a long-term loss. The tiled box for the CNN indicates that the network parameters are frozen (i.e. they are the same as the one used in the network at time ).
To train the network, the widely-used Mean Squared Error (MSE) metric [80] is introduced, which computes the L2 distance between two fields. As the loss function of the PDNN is the residual of the continuity equation, the MSE metric is applied directly to the divergence of the corrected velocity, which can be interpreted as the L2 distance to zero. Since the solution of the Poisson equation is the pressure correction term, the network output has to be coupled with the fluid solver to produce the desired corrected velocity and compute its divergence and L2 norm. Note that the gradient of this correction (Eq. 2) is known, so that a classical backpropagation can still be employed. However, if the network is only trained by evaluating the divergence level after a single timestep, no information is given to the network on how its error can be amplified by the non-linear advection at the next timestep. To circumvent this issue, a new term is introduced into the PDNN loss. This new term, known as long-term loss, corresponds to the divergence field of the uncorrected velocity field obtained several time steps after the initial prediction (Fig. 5). To do so, the network is trained inside the CFD solver: the prediction of the CNN at time t is fed through the advection and the CNN-based Poisson solver to compute the next n timesteps. The resulting divergence is computed (long-term loss in Fig. 5), and added to the initial divergence (short-term loss). The network weights are then updated using a gradient descent algorithm. Note that in the backpropagation algorithm, the chain rule [81] is used, yet here the advection part is not differentiated, and therefore does not appear explicitly in the gradients. The number of time steps to compute the long-term divergence term is either 4 (90% of the times) or 16 (10% of the times). As a summary, the total loss function is: where and are two hyperparameters controlling the relative importance of the short-term and long-term losses.
As a semi-supervised training procedure is considered, a complete training dataset is not necessary to train the network since no 'ground truth' target is needed. However, physical initial conditions are still possible, in particular to enforce the spatial coherence and structures relevant from a physical point of view, and to avoid overfitting. Thus, even for this PDNN, a training dataset was computed with the open-source code Mantaflow [57]. This dataset consists of closed boundary domains, with randomly initialized velocity fields and no density variations (i.e. the Richardson number is = 0 for all training examples). Random geometries are placed in arbitrary locations, and velocity divergence sources are introduced in the domain to initialize the flow motion. As a result, the dataset consists of 320 simulations on a 128 × 128 uniform cartesian grid. Each simulation is run for 64-time steps.

Hybrid Methodology
The main difficulty when building a Machine Learning-based CFD solver is to guarantee its precision and robustness, in particular for cases far from the learning database. To do so, Ajuria et al. [39] introduced a hybrid strategy, where the network error is tracked over time. Depending on this error level, the network prediction is used as an initial guess for the Jacobi method to further improve its quality. Note that a similar procedure was also proposed by Hsieh et al. [82], where an iterative solver is coupled to a deep linear network to guarantee the accuracy level of the network prediction. An overview of the hybrid method is shown in Fig. 6: after the advection step, the pressure correction term is predicted by the deep neural network architectures proposed in Section 3.1. At this point, no guarantee is given on the precision of the method, possibly leading to large errors or even numerical instabilities [39]. The hybrid strategy is then applied, based on an error associated with the divergence velocity field, denoted E ∞ . Ajuria et al. [39] proposed to use the maximum of the divergence velocity field, i.e. E ∞ = |∇ · u|, showing a significant improvement of the solver accuracy and robustness. Note however that no physical evidence of this choice is given. In particular, there is no formal proof that reducing the maximum of the velocity divergence actually leads to a more physical behavior of the simulation. One objective of this paper is to propose a choice of E that guarantees a physical behavior of the simulation, towards a fast, reliable and robust CNN-based incompressible solver. Figure 6: Sketch of the hybrid method, which is activated depending on the error E compared with the threshold value E .

Error Quantification
In order to compare the performance of the several architectures defined in Section 3.1, a detailed analysis of the network errors is required. Indeed, evaluating and comparing the time of inference of a method, in particular with neural networks, is relevant only if performed at a fixed error level. Thus, this section intends to characterize the accuracy of each architecture, so that a fair comparison of performances can be achieved in Section 5.

Preliminary results without hybrid approach
First, the four networks described in Section 3.1 are tested without the hybrid approach on a plume test case with a Richardson number = 14.8, in order to evaluate the network accuracy and its generalization capabilities (since the training dataset was obtained for = 0). As the Boussinesq approximation is used, the density variation should remain small and is set here to Δ / 0 = 0.01. The lighter fluid is injected at the inlet boundary with the velocity = 0.01 −1 . The CFD domain size is = = 512 , discretized on a uniform cartesian grid of 512 × 512 cells. The gravity is set to = 0.002 −2 , and the inlet radius is = 74, 25 , taken as the characteristic length to define the Richardson number. The divergence is normalized with the inlet velocity and the characteristic length, and time is normalized as well by the characteristic time needed for the plume to reach the top boundary at = ∞, i.e., in the pure jet configuration with no buoyancy. These normalized quantities are therefore defined as h yh x Figure 7: Sketch of the plume-cylinder configuration.h andh are the coordinates of the plume head location.
The plume head position is defined as the highest point of the plume at a time (Fig. 7). On the vertical axis, this position is notedh = ℎ / . In the case of the plume impinging the cylindrical obstacle, the plume head position along the horizontal axis is also measured, denoted byh = ℎ / (Fig. 7). The vertical positionh allows the identification and comparison of the plume velocity when rising due to both the advection by the jet and the buoyancy forces. The horizontal positionh measures the flow deviation induced by the cylinder. This measurement is performed in both the left and right sides of the CFD domain in order to quantify the asymmetry of the flow field. Analyzing the plume head's position for the several networks on the two test cases presented in Fig. 1 reveals that deep neural networks are able to solve the Poisson equation in an Euler incompressible solver with good precision, even in cases not seen during training (Fig. 8, tested at = 14.8). However, slightly different plume developments are obtained depending on the network architecture. As introduced in previous work [39], the Jacobi solver is taken as the reference for the two studied test cases. The accuracy of the Jacobi solver solution increases with the number of iterations, but as the convergence rate does not follow a linear behavior, a trade-off between the number of iterations and the desired accuracy is usually necessary. The grey zone in Fig. 8 represents the solution range obtained with between 200 and 10000 iterations. Note that the plume development is slower with higher number of iterations. While the solution deviation of the test case without obstacle remains narrow, the cylinder test case shows a larger dispersion. Due to computational cost, it was however not possible to perform more than 10000 iterations of the Jacobi solver.
For the no-cylinder case, the MonoScale network results in a faster head propagation, which falls outside the confidence interval of the Jacobi solvers, whereas all other network architectures produce an overall good physical behavior of the plume. The same trend is observed in the case with the obstacle, with a faster propagation and an early flow deviation for the MonoScale. Note that the symmetry is preserved by all networks since no difference is observed between left and right horizontal positionsh . The final horizontal length of MonoScale does not match the ones of the other networks. This is due to a wrong prediction of the flow structure and roll-up of the plume by the MonoScale network. This is highlighted in Fig. 9, showing the density iso-contours (white lines) and the percentiles of the divergence error (gray background field) at˜= 0.29 (no-cylinder case) and˜= 0.41 (case with obstacle). A significant difference of the plume shape is observed between the MonoScale and the other networks. The MonoScale network predicts larger vortices, probably due to an inaccurate pressure field prediction which negatively affects the baroclinic torque. However, if the divergence error is analyzed, the MonoScale follows the most regular distribution, whereas the rest of the networks seem to struggle in the surroundings of the cylinder. Moreover, the Unet network shows a quite irregular divergence distribution, with problems at boundaries. The complete evolution of the plume for both cases are provided in A (Figs. 19 and 20). In order to provide robust and reliable CNN-based incompressible solvers, Ajuria et al. [39] proposed to track the network error in time, and introduced a hybrid strategy to improve the network predictions if this error level becomes unacceptable. Previously, the head location has been employed for measuring the error, yet it is case-dependent and requires a reference simulation, which prevents this type of error to be used in new unknown simulations. The percentiles of the divergence error displayed in Fig. 9 highlight the locations in the flow field where the neural networks have difficulties, but it cannot be used as an absolute error level as well. Thus, the objective is to define an absolute unsupervised error E, which can be evaluated for any plume case without reference, and is consistent with the physicalbased error given by the head position. Here, two measures are considered: (i) the maximum of the divergence velocity field E ∞ = (|∇ · u|), and (ii) the mean of this divergence field, E 1 = (|∇ · u|). This type of error control is typically used in iterative algorithms, such as the Jacobi method or the Conjugate Gradient method [83]. This assumes a uniform error distribution, which implies that the maximum or mean error value is representative of the overall flow behavior. However, this hypothesis does not necessarily hold for neural networks, since their error distribution is usually non-uniform: neural networks may lead to a very high localized error, near a boundary condition for example. Figure 10 shows these two errors E 1 and E ∞ for both test cases without (a) and with (b) obstacle, for the various network architectures. Interestingly, the 'best' network is different for the two types of error. With the error based on E ∞ , the MultiScale network outperforms the other architectures, in particular for the case with cylinder, whereas the SmallScale leads to the highest error. Regarding E 1 , the Unet gives a low error in both cases, similar to the MultiScale network. However, based on this error, this is the MonoScale which is the worst network in both situations, with the highest mean divergence. These test cases reveal how critical is the choice of the absolute error E when comparing and then choosing network architectures.

Controlling the error level using the hybrid approach
Before assessing the performance of each neural network in Section 5, it is essential to ensure that the hybrid approach is able to control the physical targetsh andh . To do so, this Section intends to determine if E 1 or E ∞ is able to control the physical behavior of the simulation. In other words, a proper definition of the error E should guarantee that all networks hybridized with the same threshold E < E lead to the exact same evolution ofh andh in time.  First, a test is carried out in order to choose between E 1 and E ∞ as a proper indicator of the flow behavior, characterized byh andh . The lowest values obtained in Fig. 10 are chosen as threshold values. With these thresholds, the four different networks are evaluated and compared with the Jacobi method using a number of iterations also driven by the same threshold. Results on the divergence errors and head's position are displayed in Figs. 11 and 12 for the two test cases. First, it can be noticed that the specified thresholds are correctly followed by the networks and the Jacobi methods, since for the maximum error (Sub-figure-a in Figs. 11 and 12) and mean error (Sub-figure-b in Figs. 11 and 12), the divergence curves are superimposed. This demonstrates how the hybrid approach is able to guarantee a level of accuracy for the neural network predictions, which are otherwise unreliable. The resulting head positionsh andh , representative of the physical behavior of the simulations, are also provided (Sub-figure-c and Sub-figure-d in Figs. 11 and 12). Interestingly, it is found that imposing the same mean error to all methods yields exactly the same time evolution of the physical targets, for both cases with and without obstacles. However, imposing the maximum error as a threshold leads to different, yet close, simulation behaviors. Typically, the same trend as in Fig. 8 is observed, with a too rapid rise of the plume when predicted by the Jacobi method (★) and the MonoScale network ( ). This faster convection speed is highlighted in the second case, where the plume is impinging the cylinder obstacle at shorter times for these two methods, resulting in an early flow deviation characterized by a non-nullh position (Fig. 12-c). Note however that prescribing a lower maximum error, typically E = 0.18, yields also a similar plume evolution for all networks (B).

Analysis of the error distribution
To further understand how the error definition (E 1 or E ∞ ) is driving the overall behavior of the simulations, the time-evolution of the spatial distribution of the error is studied. To do so, at a given normalized time˜, the distribution of the divergence error is computed by a Kernel Density Estimation (KDE, [84]). The evolution of the KDE for the plume case without obstacle is displayed in Fig. 13. The same analysis on the second test case corresponding to the plume impinging the cylinder is performed in C. Results are displayed at four timesteps:˜= 0.1, 0.2, 0.29 and 0.39.

13
Both E ∞ (top) and E 1 (bottom) have been tested. The corresponding threshold values E have been set to the lowest values of the error obtained by all networks without hybrid approach (Section 4.1). When following the maximum divergence threshold (top), it can be noticed that the error distributions differ depending on the method employed. A significant difference is observed at˜= 0.1 for the Jacobi iterative method, for which the error distribution is spreading over a wide range, from 0 to 0.3. In comparison, all neural networks at this normalized time produce an error close to 0.05, resulting in a sharp unimodal density function. Note that such behavior is ideal since the error is well controlled in space, yet the threshold value (E = 1.5 at˜= 0.1, Fig. 11-a) is not directly linked to the most probable error (E ≈ 0.05 at˜= 0.1). At the same time˜= 0.1, using E 1 (Fig. 13, bottom) generates exactly the same error distribution for all approaches, including the Jacobi iterative method. This result generalizes the conclusions established from Figs. 11 and 12: when specifying the error as the spatially averaged E 1 error on the CFD domain, the hybrid approach yields exactly the same error distribution (whatever the method used or the network architecture), and consequently produces the same evolution of the plume head positionh, i.e., the same physical behavior of the simulation. Interestingly, the error distribution is similar to the case following E ∞ : the probability density function of the divergence error is unimodal, with a peak located at |∇ ·˜| ≈ 0.05. However, compared with the use of the maximum error, here the threshold value (E ≈ 0.05 at˜= 0.1, Fig. 11-b) directly corresponds to the most probable error.  As the simulation is running, a clear advantage of using E 1 can be observed at˜= 0.2. Not only it yields the same distribution error for all methods, but this distribution is close to the one obtained at previous times (e.g.˜= 0.1), revealing that the hybrid approach is consistent in time. On the contrary, the use of E ∞ leads to very different error distributions. It is worth noting that the Jacobi method (★) and the MonoScale ( ) exhibit a distribution shift towards higher error values, thus explaining their higher mean values ( Fig. 11-a). This result is consistent with the unphysical behavior of these two simulations (too fast rising of the plume), since the evolutions ofh are strongly correlated with the mean value error (Figs. 11-b and d).

< l a t e x i t s h a 1 _ b a s e 6 4 = " t 4 r F C r Z g Q g A D J a V o j v W 4 h 0 + E t e Y = " > A A A C 2 X i c j V H L S s N A F D 3 G 9 z s + d m 6 C R X B V U q n o s i i C S w V b C 6 2 U y T i 1 Q / M i m Q i 1 d O F O 3 P o D b v W H x D / Q v / D O m I J a R C c k O X P u P W f m 3 u v F v k y V 6 7 6 O W e M T k 1 P T M 7 N z 8 w u L S 8 v 2 y m o t j b K E i y q P / C i p e y w V v g x F V U n l i 3 q c C B Z 4 v j j 3 u o c 6 f n 4 t k l R G 4 Z n q x e I i Y F e h b E v O F F E t e 7 0 Z M N X h z O 8 f D V r 9 p g z b q j d o 2 Q W 3 6 J r l j I J S D g r I 1 0 l k v 6 C J S 0 T g y B B A I I Q i 7 I M h p a e B E l z E x F 2 g T 1 x C S J q 4 w A B z p M 0 o S 1 A G I 7 Z L 3 y v a N X I 2 p L 3 2 T I 2 a 0 y k + v Q k p H W y R J q K 8 h L A + z T H x z D h r 9 j f v v v H U d + v R 3 8 u 9 A m I V O s T + p R t m / l e n a 1 F o Y 9 / U I K m m 2 D C 6 O p 6 7 Z K Y r + u b O l 6 o U O c T E a X x J 8 Y Q w N 8 p h n x 2 j S U 3 t u r f M x N 9 M p m b 1 n u e 5 G d 7 1 L W n A p Z / j H A W 1 n W K p X N w 9 L R c q B / m o Z 7 C B T W z T P P d Q w T F O U C X v G z z i C c 9 W w 7 q 1 7 q z 7 z 1 R r L N e s 4 d u y H j 4 A v m u Y T A = = < / l a t e x i t >
At larger timesteps (˜= 0.29 and 0.39), the error distribution is spreading for all cases toward higher error levels (from 0 to 0.3 for most cases). Again, defining the error and threshold values using E ∞ (top) is not able to control the whole error distribution. The Unet (•) produces the most unimodal distribution with a peak at a low error of 0.05, while the Jacobi method (★) and MonoScale network ( ) generate a flat error distribution with larger errors up to 0.6. However, following E 1 (bottom) is still producing consistent error distributions, with low discrepancies between the various methods, including the MonoScale and Jacobi. A slight shift of the error distribution can be observed, with a most probable error close to 0.1 (instead of 0.05 for early times). This is due to the fact that the present test case is performed with a varying threshold value E (˜): this value is chosen as the lowest error value obtained without the hybrid approach, which is therefore not fixed in time. Figure 11-b reveals that the threshold value of the mean error is actually increasing in time, from 0.05 to early timesteps towards 0.1 at the end of the simulation (˜= 0.6). Overall, the description of the error E and threshold value E using the spatially averaged error over the CFD domain has shown a remarkable ability to control the whole distribution error in time and space, whatever the method used. For all times, the distribution is unimodal, with a most probable error in the range 0.05 − 0.1, in good agreement with the threshold value of the mean error employed for the hybrid approach. This is not surprising since the error distribution is a unimodal symmetric probability function, for which the most probable and mean values are equal. In other words, prescribing the mean error is constraining the whole error distribution, whereas following the maximum value E ∞ only constrain the distribution tail, which seems insufficient to control the physics of the system, in particular here the rising speed of the plume head. Note that similar trends and conclusions are observed with the plume impinging the cylindrical obstacle ( Fig. 26 in C), even if the error distributions are no more unimodal: the hybrid approach based on the mean error is still able to control the whole distribution error whatever the network architecture. For completeness, the whole spatial distributions of the error are displayed for the two cases (with and without obstacle) and both error definitions (E 1 and E ∞ ) in A.
Finally, this study reveals that the neural networks hybridized with a Jacobi method are able to guarantee the accuracy level of the whole simulation in time, in particular when choosing an error defined as the spatially averaged divergence error in the CFD domain. Consequently, in the following, the hybrid approach based on E 1 is chosen to ensure that all networks have a similar accuracy with the same plume head evolution in time, a pre-requisite to perform a fair comparison of the network performances.

Performance assessment
In this section, the performance of neural networks is assessed for the time of inference, i.e., the time needed to produce the pressure correction. The analysis is innovative in two ways: (i) a fair performance comparison is performed using the hybrid approach, allowing the assessment of the time of inference at a fixed error level whatever the method (Jacobi or neural network architecture), and (ii) these performances are evaluated for several grid sizes. Note that even if CNN are theoretically capable of dealing with CFD domains of arbitrary size and resolution, they have in practice difficulties which such cases, probably because of boundary effects [85]: here again, the hybrid approach is useful to guarantee the error level of the solution when changing the grid size. To do so, as suggested in Section 4.2, the mean divergence error is used to define the threshold E . All simulations were performed on the same GPU card, namely here a NVIDIA Tesla V100 with 32 Gb of memory. Note that the performance evaluation on multiple GPUs is out of the scope of the present study.
To measure the fluid solver performance, the following times are defined: • : Inference time taken by a neural network or the Jacobi solver to output one pressure field when the divergence field is inputted, i.e., not taking into account the time to correct the velocity field or the time spent on the extra Jacobi iterations of the hybrid solver.
• : Time taken to perform the entire pressure projection step. This includes the time to perform the first pressure inference, the extra Jacobi iterations in the hybrid process, as well as the correction of the velocity field.
• : Time taken to perform an entire iteration of the fluid solver, i.e., the advection and pressure projection steps.

Network Comparison
Using the hybrid strategy, a non-trivial trade-off has to be made between accuracy (usually implying more complex and deeper architectures) and fast inference time (requiring small networks). In practice, a less accurate network with a fast inference time will require more Jacobi iterations to reach the desired accuracy level, thus limiting drastically its performance. Comparing the performance of the several network architectures used in this work will provide a first insight on this trade-off, in order to establish guidelines for future developments of AI-based solvers.
Consequently, the performance assessment of the neural networks is split in two steps: (i) assess the inference time of each CNN to produce the output, and (ii) evaluate the number of Jacobi iterations needed by each network, using the CNN prediction as initial guess, to reach the target accuracy. Table 5.1 shows the inference time of each network, and reveals that is not directly associated with the number of parameters of the network. Indeed, whereas MonoScale, MultiScale and Unet have the same number of parameters, they also have very different inference times. Similarly, the SmallScale contains three times fewer parameters compared with Unet, but still requires more time to output its prediction. This is due to the specific architectures of the Unet and MultiScale networks which involve several scales: filters and their associated parameters at low resolution (capturing the large scales of the flow field) are applied on reduced inputs, therefore limiting the number of operations to perform, and thus reducing ( Fig. 14-b). Thus, the more parameters at lower scale, the faster the network is. This is confirmed by the associated number of floating-point operations (FLOP) for each network architecture, displayed in Fig. 14 for each network and several grid sizes, from 1024 to 4.2 × 10 6 cells. It reveals that the FLOP number depends linearly on the grid resolution and is directly correlated with the inference times of Tab. 1: the 5 scales of the Unet result in a low FLOP number, and therefore a reduced , comparable to the SmallScale network with three times fewer parameters. The MultiScale network also contains 3 scales, but most parameters are concentrated at the larger resolution, which is therefore dominating the FLOP number and inference time. Only a small gain is observed for this architecture compared with a standard MonoScale network. Since the several networks have different error levels, the number of Jacobi iterations to reach the desired accuracy level is decisive in the overall performance of the AI-based code. Figure 15 shows the number of Jacobi iterations needed to reach the threshold E , for both cases without (Sub-figure-a) and with (Sub-figure-b) obstacle. As expected, the classical Jacobi solver (★) and the MonoScale network ( ) need numerous iterations, typically between 300 and 400. Using a multiscale architecture (either MultiScale or SmallScale) improves significantly this situation, since only 100 to 200 iterations are needed for the same accuracy level. Note that for the case with obstacle, the MultiScale network needs almost no Jacobi iteration after the plume has impinged the cylindrical obstacle (˜> 0.3). In contrast, the Unet outperforms all networks, with almost no Jacobi iteration for all cases. However, it requires 50 − 100 iterations after the plume has impinged the cylinder (˜> 0.3), which is consistent with the previous finding that Unet has difficulties close to walls (Fig. 9). As a general conclusion, all networks outperform the classical Jacobi solver, proving that using neural networks as an initial guess to the iterative solver is relevant. However, this AI-based computational strategy is efficient only if the time needed to produce the initial guess by the CNN leads to a total time lower than using the Jacobi solver alone.
To verify this, the time to solve the Poisson equation ( ) is now compared for all methods in Fig. 16-a. As a comparison, the computational time to solve the advection equation is = 35.5 . The bar plot ( Fig. 16-b) highlights the acceleration factor = / for each network architecture, where and are the computational times of the Jacobi method and the hybridized network respectively. First, it can be noticed that all AI-based approaches are faster than the Jacobi solver alone, which proves the efficiency of the proposed computational strategy using CNN to predict the solution of the Poisson equation. As expected, the MonoScale only provides a small gain, about 1.2 − 1.3, since (i) it needs numerous iterations of the Jacobi solver to achieve the desired accuracy level, and (ii) it has the largest inference time (Tab. 1). Note that the latter is evaluated at 14.6 , which indicates that even for this network with large inference time, it is still negligible compared to the computational time of Jacobi iterations. Interestingly, the large (MultiScale) and small (SmallScale) multiscale networks achieve a comparable acceleration factor, respectively ≈ 3.7 and 2.4 without obstacle, but differ when the plume impinges the cylinder, with ≈ 7.5 and ≈ 3.1. Indeed, the larger network achieves a better accuracy, thus limiting the number of Jacobi iterations. However, the accuracy-inference time trade-off is not trivial when dealing with deep neural networks embedded in a CFD solver, as in test cases where few Jacobi iterations are needed, the SmallScale could outperform the larger MultiScale network, thanks to the smaller network inference time. Finally, the Unet architecture outperforms all other methods in the case without obstacle with an acceleration factor = 26.1. As Unet has difficulties handling boundaries, its performance drops in the case with obstacle, yet still accelerating the CFD solver by a factor = 6.6, slightly below the MultiScale network performance. Note however that for longer simulations, the MultiScale network would further outperform Unet since requiring no Jacobi iteration for times larger than˜= 0.3.
(a) (b) Figure 16: Time (a) and acceleration factor (b) for the non-cylinder test case ( ) and cylinder test case ( ), for the four studied networks: MonoScale, MultiScale, × SmallScale, and • Unet, as well as the ★ Jacobi solver.

Performance evaluation on various grid sizes
The previous performance assessment was made at a fixed resolution of 2.6 10 5 cells. Recent work [85] revealed that networks can encode spatial location and resolution on the training dataset, which can lead to poor performance when applied to a different resolution than the one used during the training phase (2.6 10 5 cells in this work). This drawback is here tackled by the hybrid approach, which ensures the desired accuracy level whatever the grid size. It is also found in this work that the present CNN architectures are capable of generalizing to different grid sizes, even if no clear explanation can be given for this specificity. As a reminder, Fig. 14 showed that the number of operations per network evolves linearly with the grid size. However, Fig. 17 shows a different behavior when looking at , with two distinct regimes: (i) the network performances follow the same behavior as in Fig. 14 when the domain size is large enough, with a linear increase with grid size (ii) remains constant for small grids. This phenomenon is related to the CPU overhead, as the time taken by the GPU kernels is inferior to the command transmission between the CPU and the GPU, which depends only on the network architecture but not the grid size. Therefore, networks with a simpler architecture result in a lower CPU overhead, and so a faster inference. It should be noted how for example the MultiScale and SmallScale networks show a similar inference time when the GPU kernels are saturated, as the architecture and layer size are exactly the same. This saturation should be taken into account for small resolution test cases, for instance, by designing dedicated simple architectures.

Analysis of the network performances
In order to better understand the link between the network architecture and its performances, the time taken for different scales is displayed in Fig. 18, for both the MultiScale (a) and Unet (b) networks when the spatial resolution is varying. This figure shows the evolution of the 3 main scales of the MultiScale network, depending on the spatial resolution. As expected, the largest scale takes the major part of the inference time for both networks, performing most of the operations due to the larger size of the domain. Moreover, it is shown that the smaller scales are more prone to the command overhead, which makes the saturated regime occurring on a larger resolution range, for instance, the scale n/4 of the MultiScale is saturated until the resolution 512 512, and the scale n/16 of the Unet until 1024 1024. Exploiting smaller scales with more parameters could open the path to new computational strategies with large computational domains divided into multiple smaller subdomains of the appropriate size, for example, adapted to each scale of the network.

Conclusions
This work focuses on the resolution of the Poisson equation with convolutional neural networks applied to incompressible fluid flows. A plume test case, with and without obstacle, is chosen since parametrized by a single dimensionless quantity: the Richardson number. The well-known MultiScale and Unet architectures are analysed, and compared to a simpler architecture as well as with a traditional Jacobi solver. In order to ensure a user-defined accuracy level, a hybrid strategy is used where the neural networks are coupled with a Jacobi solver when a threshold error criterion is not satisfied. This work shows that a threshold based on the mean divergence of the flow field ensures a consistent physical behavior of the fluid flow. The various networks are compared in both accuracy and inference time when using the hybrid strategy. It is found that the network behavior varies with the studied case. When there is no obstacle in the domain, the Unet network outperforms the other networks both in accuracy and performance. However, when an obstacle is introduced, its accuracy decreases because of difficulties close to walls, resulting in a simulation time superior to the MultiScale network. For both cases, the MultiScale and Unet provide speed-ups of order 5-7, demonstrating the potential of such a method. Moreover, it is shown that the multiple scales of the Unet and MultiScale networks result in a faster inference time, as long as the resolution is larger that 256 256 . For lower resolutions, the CPU overhead controls the inference time, enabling the simpler MonoScale network to obtain the fastest inference time. As a conclusion, this work has revealed the potential of coupling CNN with multiple scales to a CFD solver to produce fast and reliable flow predictions.

A Spatial distribution of the divergence error
To complement this study in Section 4, the divergence distribution at four timesteps is shown here. Note that these divergence percentiles show the spatial distribution of the divergence, but they do not give any information about the absolute value of the studied error. Figure 19 represents the divergence distributions at 4 timesteps for the no cylinder case where no hybrid strategy is used, i. e., the divergence distribution of the raw output of the four networks. Qualitatively, all networks behave similarly, as the main divergence sources are found on the plume head and vortices. However, when closely analysed, each network counts with its distinct behavior related to the network architecture. The MultiScale and SmallScale architectures behave similarly, showing a smooth divergence distribution, with the presence of a divergence stride close to the domain edges, probably related to the network padding, while the Unet network seems to struggle predicting the pressure field close to the walls (particularly the top and left wall). Moreover, the MonoScale network shows the smoothest divergence distribution, even if its flow topology clearly shows a faster plume propagation.  Similarly, Fig 20 represents the same divergence distribution but on the cylinder case. The divergence distribution follows a similar behavior as the one observed in Fig. 19. However, the presence of the cylinder highlights potential difficulties of networks to predict flow-wall interactions. In particular, it is shown that the Unet network particularly faces such a difficulty, as the major part of its error concentrates near the cylinder wall. Analysing the flow topology, shows that the MonoScale network creates larger structures, probably due to an incorrect prediction of the baroclinic-torque.  show the divergence distribution and density iso-contours of the four studied networks, as well as the Jacobi solver with the hybrid approach. Figures 21 and 22 correspond to the no cylinder test case with E = E ∞ (Fig. 21) and E = E 1 (Fig. 22). As the threshold level is set to the minimum value obtained by any networks without any Jacobi  Figure 20: Divergence percentiles and density iso-contours (in white) of the four studied networks for the cylinder case at timesteps˜= 0.14, 0.27, 0.41 and 0.55 when the networks prediction are not coupled with the hybrid strategy.
iterations, the Unet network is usually the one needing the least number of Jacobi iterations. Thus, the Unet's output closely matches the raw output, which shows the error concentration on the top and left wall, yet the global flow topology is consistent with other networks. Figures 23 and 24 correspond to the cylinder test case for both thresholds E = E ∞ (Fig. 23) and E = E 1 (Fig. 24). As the less accurate network, the MonoScale network is mainly controlled by the Jacobi iterations, which results in similar flow topologies between the MonoScale and the Jacobi solver. As previously mentioned, the Unet network concentrates its error around the cylinder, as it struggles to correctly predict the pressure field around objects. This results in a quite different error distribution, even when E = E 1 , as even if the Unet network needs some Jacobi iterations to reach the specified divergence threshold, the resulting error field is still controlled by the error concentrated around the object and walls. However, and as previously mentioned, even if the divergence distribution is different, the flow topology perfectly matches the resulting flow of the other networks.
So to sum up, it is clear that following the maximum error with the hybrid approach yields different plume behaviors, with a too fast rising for the MonoScale network and the Jacobi method, whereas the mean error provides consistent numerical simulations whatever the network architectures. For all cases and methods, several zones of high errors can be distinguished: (i) at the inlet where the lighter fluid is injected, (ii) above the plume head, and (iii) at the core of the vortices produced by the baroclinic torque, in particular at later times (˜= 0.29 and 0.39) for the no cylinder test case. In that context, the Unet is producing additional error patterns, with high errors close to the boundaries, either near the cylinder walls or close to the CFD domain edges. Moreover, the error distribution in space is not symmetric, compared with all other methods. Note however that percentiles do not reflect the absolute error level, and the Unet has shown a good overall accuracy for both cases, even in the presence of the cylindrical obstacle.  Figure 24: Divergence percentiles and density iso-contours (in white) of the four studied networks for the cylinder case at timesteps˜= 0.14, 0.27, 0.41 and 0.55 when the networks prediction are coupled with the hybrid strategy. The chosen threshold is E = E 1 , where E = (E 1 ) obtained by any network at a given timestep, without being coupled with the hybrid strategy.

B Sensitivity of the hybrid solver on the threshold value E
The evolution of the head position with a lower threshold is shown in Fig. 25 for the plume test case without obstacle. If the maximum threshold level is set to E = 0.37 ( Fig. 25-a), which corresponds to a divergence threshold lower than the maximum threshold that could be obtained with any network (Fig. 10), the plume does not follow a homogeneous behavior. In particular, the Jacobi method and MonoScale network produce a higher divergence level. If the threshold is further decreased by a factor 2 (E = 0.18, Fig. 25-b), all methods now produce simulations with the same evolution of the plume head position. To further understand the condition which dominates the flow behavior, the mean divergence of the velocity field at each time step is displayed with dotted lines in Fig. 25. It shows the evolution of E 1 for both threshold values E = 0.37 (a) and E = 0.18 (b). As the maximum threshold is decreased, the spatially-averaged error gets similar for all networks, resulting in the same plume evolution in time. Note that discrepancies are still observed between these errors, the Unet outperforming all other networks. C Evolution of the error distribution for the plume-cylinder case Figure. 26 shows the KDE for the plume impinging a cylinder, for the cases where the threshold E is equal to both E ∞ and E 1 for four timesteps (˜= 0.14, 0.27, 0.41 and 0.55). As previously mentioned, when E = E ∞ , the four networks, as well as the Jacobi solver, show a different behavior. This can be particularly appreciated at timestep˜= 0.27, although all the networks follow a unimodal distribution, the MonoScale network and the Jacobi solver, previously identified as the less accurate solvers, have a divergence peak around |∇ · u| = 0.2, whereas the MultiScale and SmallScale networks behave similarly with a peak around |∇ · u| = 0.1, and the Unet network outperforms the rest with a peak around |∇ · u| = 0.05. When E = E 1 , the MonoScale, MultiScale, SmallScale and Jacobi solver follow an almost identical distribution, showing a divergence peak at around |∇ · u| = 0.1. The Unet behavior should be highlighted, since its distribution has a wider peak, centered on around |∇ · u| = 0.05. This is related to the threshold level, which follows the lowest value of the mean divergence of all the studied networks. At timestep˜= 0.27, this value corresponds to the Unet network, so to achieve this threshold level, the other networks rely on Jacobi iterations which homogenize the divergence distribution, whilst the Unet's KDE corresponds to the network's output. Despite this difference, the plume flow is not affected, as macroscopically the flow structures are identical for all the cases when E = E 1 .