A physically consistent AI-based SPH emulator for computational ﬂ uid dynamics

: The integration of arti ﬁ cial intelligence (AI) into computational ﬂ uid dynamics (CFD) has signi ﬁ cantly expanded the scope of ﬂ uid modeling, allowing enhanced analysis capabilities and improved simulation performance. While Eulerian methods already bene ﬁ t extensively from AI, notably in reliable weather prediction, the application of AI to Lagrangian methods remains less consolidated. Smoothed particle hydrodynamics (SPH) is a Lagrangian mesh-less numerical method for CFD with well-established advantages for the simulation of highly dynamic free-sur-face ﬂ ows. Here, we explore an application of AI to SPH simulations, utilizing an arti ﬁ cial neural network (ANN) to estimate hydrodynamic forces between particle pairs, learning from SPH-simulated results. A model of this nature, which emulates the mathematical representation of physics, is termed an emulator . We examine the physical signi ﬁ cance of the emulator, presenting its applications in benchmark tests, assessing its faithfulness to traditional SPH simulations, and highlighting its ability to generalize and simulate test cases with varying levels of complexity beyond its training data.


Introduction
There is a growing interest in leveraging artificial intelligence (AI) to push forward the frontiers of computational fluid dynamics (CFD) [1], enabling new forms of analysis [2][3][4] and enhancing simulation performances [5,6].A common approach in this direction is the use of emulators, where AI models are trained over simulated data to learn the behavior of traditional CFD models and replace them, offering better performances [7,8].The use of AI-based emulators opens to new levels of analysis [9,10], allowing a better description of physical phenomena [11,12] and an efficient simulation of large-scale problems [8,13,14].However, the increasing widespread use of AI-based models raises the issue of rigorously validating the outputs, necessitating a deeper analysis of their reliability.Eulerian formulations can already provide high-fidelity and reliable results, with clear examples being the Deep Emulator Network SEarch (DENSE) model for weather predictions [8] or MeshGraphNet for fluid dynamics simulations [15].On the other hand, Lagrangian applications, especially in their mesh-free variants, are less consolidated, requiring a significant part of the design effort already to handle the absence of pre-established connections between nodes [16].This aspect is also reflected in the generated results.While some of these approaches are based on simplified CFD models with little physical meaning, many others only show qualitative comparisons to the reference CFD models, with quantitative analysis typically being out of their scope.It is, therefore, important to explore the realworld relevance and accuracy of mesh-free Lagrangian emulators compared to traditional equation-based methods.We will focus on smoothed particle hydrodynamics (SPH) [17], a mesh-free numerical method that allows to obtain high-fidelity simulations of complex fluids [18], specifically for free surface and multi-phase flows subject to large deformations and interacting with complex geometries [19].
Among the major groups of approaches that can be found in the literature for Lagrangian emulators, one is based on the use of graph neural networks (GNNs), where the graph is a direct way to compensate for the lack of a spatial structure in the dataset.GNNs are built so that particles are treated as the nodes and the interactions between them are represented as the edges.Specifically for the context of fluid dynamics, Lagrangian fluid graph networks [20] have been developed [21].Despite this approach can effectively cope with unstructured datasets, it affects the mesh-free nature of the reference SPH method.It introduces the need for re-meshing after large deformations, a computationally expensive job that goes against some of the possible applications of emulators.
Another major kind of approach for Lagrangian emulators uses convolutional neural networks (CNNs).These take into account the spatial distribution of the Lagrangian nodes by performing spatial convolutions over restricted neighborhoods using predefined masking functions.Originally, CNNs are designed to work with structured grids of data, as would be the pixels of an image or the nodes of an Eulerian numerical method.Therefore, their application to Lagrangian mesh-free methods is not straightforward and needs specific designs [22,23].These approaches allow a high-level emulation and can replace the process of iterating over the neighboring particles, one of the fundamental tasks in SPH simulations, altogether.However, these models show a strong dependence on the spatial features of the training data (e.g., global position of the particles), which can affect the quality of the results and the generalization capabilities.This is the case when some particles' configurations are not well represented in the training dataset.
A third kind of approach treats the emulation locally, at a particle level or particle-pairs interactions.In this case, the spatial distribution of the particles is handled with iterations over neighbors, as in the original SPH approach.Ladickỳ et al. [7] presents an application of AI to speed-up SPH simulations, where a random forest (RF) model learns and reproduces the behavior of each particle over time steps longer than those utilized in simulation.The RF takes as input features the SPH viscous and pressure forces (computed online with the classical SPH equations), as well as some particle state physical quantities, and returns an update of the particles state relative to the time step presented at training that can be longer than those used for the reference simulation.The emulator delivers good speed-up of the method with real-time performances.However, one of the main limits of RF models for regression tasks is the inability to extrapolate output values outside of the data ranges observed during training, which is a limitation to the consistency of the method and to generalization possibilities.A different solution is presented by Alexiadis [16], where an artificial neural network (ANN) is used to estimate the force between particle pairs.The inputs to the network are the relative particle distance, relative velocity, and the two particles densities.Because the network works on particle pairs, in this approach, the global spatial distribution of the particles in the training dataset is not as substantial as it is for CNN.This is a case where the concept of emulator is applied in a strong sense, as the network plays the same role of the reference equations.In the work by Alexiadis [16], different CFD methods are emulated (SPH, molecular dynamics, and discrete element method), showing good fidelity of the results to those of the reference models and good generalization capabilities.The paper also shows an example of model inversion, as the ANN is trained to estimate particle forces while only knowing total forces per particle.However, the study mainly covers the technical aspects of the development and testing of the emulator, whereas the adopted reference models only take into account little physical meaning (they are so-called "vanilla" models).This does not make possible a direct applicability of the emulator for practical use.
Here, we present an SPH emulator based on ANN, developed for physically consistent simulations, where the reference model is in line with SPH formulations adopted in practical applications.This study emphasizes the design aspects that have a strong impact on the physical fidelity of the reference SPH model.The design of the ANN recalls the idea proposed in the study by Alexiadis [16] and extends the overall structure of the emulator to handle the more complex and realistic physics.The presented design enforces symmetry in the predicted forces by only using particle-pair properties as inputs to the network.Specifically, we feed the network the product of the particle densities as a unique feature, as opposed to giving them individually.This also lightens the computational management of the network, which has a reduced number of inputs.To test the generalization capabilities, we train the emulator over a test case with intermediate complexity, represented by a dam break flow, and then we reproduce test cases with both lower and higher complexities.
In the following, we will begin by introducing the SPH model that we use as a reference for the emulator.In Section 3, we will discuss the design and training of the emulator, and Section 4 shows how the emulator reproduces the same problem used for training and how it generalizes to new cases.Finally, we will conclude with remarks in Section 5.

SPH reference model
SPH [17] is a particle-based Lagrangian mesh-free numerical method with growing applications to fluid dynamics.The Lagrangian nature of the method allows an accurate treatment of flow interfaces, including free surfaces and interactions with complex geometries [24].The absence of any predefined interconnection between nodes (i.e., meshes) allows efficient management of highly dynamic flows undergoing strong mixing [25].This model has shown to be very flexible in its application to diverse scientific and technological fields.As a part of its flexibility, SPH comes in many formulations [26,27], including numerical corrections [19,28,29], which can be chosen according to the characteristics of the simulated problem and the level of abstraction that one wants to adopt with respect to the simulated physics.

SPH formulation
The SPH model that we will take as reference is based on a quite general, weakly compressible formulation for inviscid fluids (as water is typically approximated), which is used in many validated applications [19,26,28,30,31].We will consider the Navier-Stokes equations of momentum conservation, to model the dynamics of the particles, and mass conservation, to model the evolution of particles density.Our model will be developed in 2D.In this formulation, the SPH discretization of the equation for mass conservation is written as follows: where the second term on the right-hand side stands for the density diffusion correction, presented in [32], and it is written as follows: (2) with = ε 0.01 h .The 2D equations of momentum conservation are written as follows: where .81 T is the gravity vector and Π ij is the artificial viscosity term, which helps reducing numerical high-frequency components of the velocity field, defined as follows [17]: In order to reproduce an equivalent kinematic viscosity, ν, we define the artificial viscosity coefficient α as in [17], using the relationship The pressure P is obtained from the density using Cole's equation of state [33]: where ρ 0 is the reference density of the fluid and γ is the polytropic constant that we take equal to 1.The speed of sound c 0 is taken 20 times higher than the hydrostatic velocity, i.e., , with h as the maximum height of the fluid [34].
We use the fifth order Wendland smoothing kernel [35], which has been observed to be beneficial for freesurface simulations [36].This kernel has a smoothing radius h 2 , where h is the smoothing length, defined using the smoothing factor as = ∕ α h p Δ s , and we take = α 1.33 s .

Integration scheme
We adopt a second-order predictor-corrector integration scheme described by the following steps: (1) Compute accelerations and density derivatives at instant n: and densities (predictor): (3) Compute corrected accelerations and density derivatives: (a) ( ) ) Compute new positions, velocities, and densities (corrector): The time step t Δ is computed for each particle i, and it is required to fulfill CFL-like (Courant-Friedrichs-Lewy) stability conditions [37] determined by the acceleration magnitude and speed of sound: The resulting overall time step for the model is chosen as the minimum particle time step.
A physically consistent AI-based SPH emulator for computational fluid dynamics  3

Boundary model
We adopt the dynamic boundary model [38,39], for which some parallel layers of extra particles, referred to as boundary particles, are placed along analytical boundaries.These layers are spaced at intervals of p Δ from one another, originating at the analytical boundary position.Fluid particles are thus positioned, commencing from a distance of p Δ away from the analytical boundary.These extra layers of particles should be enough to complete the support of the smoothing kernel when it is centered on the outmost fluid particle.Boundary particles are assigned a density, which is treated equally to fluid particles.The velocity of these particles is set to zero (unless a moving boundary needs to be modeled).

AI-based emulators for CFD
An AI-based CFD (AI-CFD) emulator is a model that combines the principles of CFD simulations with AI techniques to enhance fluid simulations, either in terms of performances or generated information.The AI algorithms can join an equation-based CFD model or replace parts of it and emulate the respective behavior, learning from a dataset of simulated input-output values.Incidentally, emulators are trained to mimic CFD models with specified mathematical laws, which represent the emulated entity.This is in contrast to generic data-driven approaches, where the data can be experimental and the physical laws unknown [40,41].
Applications of AI can be divided into two main tasks [42]: classification and regression.The first task aims at dividing the data into classes, or clusters, connected to each other by intraclasses similarities, which measure the homogeneity of a cluster, and interclasses dissimilarities between different clusters, which measure the distance between two clusters [43][44][45][46][47]. Therefore, the output of the model is a discrete number linked to the class membership.For the regression task, the purpose is to evince the relationships between input and output data and to make it applicable to new datasets [48].In this case, outputs are given as continuous values.For emulators, regression applications of AI are generally used, as the aim is to generate continuous simulation quantities compatible with training data.Examples of AI architectures that can be used in Lagrangian emulators for regression tasks are RFs [7] and ANNs [16,49].The first one is known for the inference speed, and the second one for better generalization ability.The two architectures are largely used for nonlinear regression tasks.

ANN-based SPH emulator
We here design an emulator for the weakly compressible SPH scheme shown in Section 2.1, where the momentum equation is replaced by an ANN, thus emulating the forces between particle pairs.
The ANN that we use is a multilayer perceptron (MLP).This is a fully connected class of feedforward ANN composed of layers of nodes (usually called neurons) with weighted interconnections.Each node computes its output from the inputs by using so-called activation functions [50].A network is typically composed of at least three layers: the first one is called input layer, the last one is called output layer, and any other layers in between are called hidden layers.The number of nodes in the input and output layers determines, respectively, the number of inputs and outputs of the network.The number of hidden layers and their nodes determines the complexity of the network and, therefore, the capability to catch details of the dataset.It is demonstrable that the standard multilayer feed-forward networks, with at least one sufficiently wide hidden layer of neurons, can approximate any kind of functions with any accuracy, and so they can be considered universal approximators [51].The process of extrapolating the relationships between input and output from a representative dataset is called training, during which the weights of the connections are set.
While the outputs of a network are clearly defined by the quantity that one is looking for, in our case, the particles interaction forces, the set of input variables provided to the ANN must be designed in order to contain all the information needed to infer input-output relationships within the dataset.This design process is called features extraction, and it is proper of AI approaches classified as machine learning techniques [42,52], as is the one presented here.These are in contrast to deep learning techniques, where features are automatically extracted from the provided input dataset by the model during the training phase [50,53].
For the design of an emulator, the relationships between inputs and output data are typically known in analytical form (the set of equations of the reference numerical model), which is a great advantage for the process of feature extraction.Looking at our equations of interest, Eq. (3), which in turn recalls Eqs (4) and ( 6), it is clear that the inputs to the network will include particle positions, velocities, and densities.Still, the design of the MLP inputs requires additional effort, primarily because velocities and positions are 2D vectors, while input variables of the MLP should be scalars.We, therefore, derive scalar quantities from the variables of interest, combining them in a way suggested by how these variables figure within the reference equations.We use three input variables, which are the relative distance between the particles, the normal component of the relative velocity, and the product of the two densities.While the first two quantities are identical to those used in [16], individual densities were used in that case as two independent inputs.Our design choice keeps the number of inputs smaller, with consequent computational advantages during training and execution of the network.Furthermore, our inputs only refer to particle-pairs quantities, forcing the network to provide symmetric interaction forces.As for the particle forces estimated by the network, they need to be vectors with the component of acceleration in the two spatial dimensions.However, considering Eq. ( 3), these vectors are, by construction, oriented as the vector of relative position between the interacting particles, which is known.Therefore, the network is designed to produce in output the module of the acceleration vector, which we then multiply by the unit vector of the relative position.The structure of the network was defined empirically during training, including three hidden layers with the following number of neurons: 9, 27, and 9.We note that this network is larger than the one used by Alexiadis [16], likely because of the higher complexity of the physics represented in the dataset.The activation function adopted for the nodes is the exponential linear unit (ELU) function [54].
As a measure to reduce the execution time of the emulator, during runtime, we substitute the use of the ANN with a lookup table that maps its input-output relationship [55], as done in the study by Alexiadis [16].This table is constructed once at the conclusion of the training phase, by sweeping the inputs of the ANN across a predefined set of values and storing the corresponding outputs.During runtime, the output values of the ANN are reconstructed from this table using linear interpolation.Although this process introduces some errors in the estimation of the output, it significantly reduces the emulator's runtime and maintains its independence from the complexity of the MLP.

Training the emulator
The dataset for the training is generated by simulating a dam break flow, with the SPH formulation presented in Section 2.1.Other than being a classical benchmark case for SPH, a dam break presents a wide range of flow conditions, evolving from the initially violent flow of the transient, toward a static steady state.Therefore, the training dataset will be representative for the different flow  6)) and considering a hydrostatic pressure profile.
In order to incorporate meaningful phases of the flow into the training dataset, we simulate 50.0 s of evolution, at the end of which the flow reaches the state shown in Figure 1.We sample the whole particle set, including boundary particles, every 0.5 s.This results into a dataset including 3.42 million of particle samples for the intermediate spatial resolution.About 80% of the total dataset is used for the training phase, while the remaining part is used for the validation phase.We found that including boundary particles within the dataset improved the quality of the results in terms of particle disorder.This is likely due to a better representativity of surface particles.
Despite the output of the MLP represents particle interaction forces, following [16], we write the loss function L i taking into account the total forces acting on particles.Specifically, the network is trained to minimize a total loss function, composed by the sum over all the neighbor lists for all the particles and time steps, namely of all the single L i , which consider the difference between the total force F i target given by the SPH model and the sum of the particle pair interaction forces f ij ANN predicted by the MLP: In this way, the sum of the estimated forces has to approximate the sum of the computed force, without knowing them individually.This approach allows to reduce the size of the training dataset, as it overcomes the need to store individual samples of the forces, in favor of storing only total forces.Total forces are computed and added to the training set for boundary particles.
Training is run for 10,000 epochs, with a learning rate of • = − α 5 10 4 [50].The Adam optimizer is used, typically reaching a final training and validation loss in the order of − 10 3 .We note that emulating the momentum equation means that we also emulate the computation of the kernel and the equation of state.Therefore, the parameters adopted in these expressions, like the p Δ , the smoothing length, or the speed of sound, are learnt during training and will not need to be specified for the emulator.On the other hand, this means that these values are fixed for the emulator, and a new model (by means of a new training) needs to be created when one needs to use different values.

Results and discussion
To test the emulator, we initially replicate the dam break problem used to generate the training dataset.Figure 2 shows a comparison of three different phases of the dam break flow, and a good agreement can be appreciated between the results obtained with the equation-based SPH model (hereafter simulations) and the AI-based SPH emulator (hereafter emulations).Some discrepancies can be seen at the level of particle configurations, especially in the regions where the details of the flow are comparable to the particle size or where the fluid behaves in a highly chaotic way (see the splashes shown in the middle and bottom frames of Figure 2).In both cases, the local flow is highly sensitive to particle positions and velocities, and even small discrepancies can generate a different content of flow features in the resulting simulations.We observed these phenomena by also reproducing the experiments, which include the three steps, i.e., equation-based SPH simulation, training, and AIbased SPH emulation, with different spatial resolutions (one lower resolution with = p Δ 1.0 m and one higher reso- lution with = p Δ 0.5 m).For both lower resolutions, where the particle size becomes more important within the domain, and higher resolutions, where more sensitive flow details are generated, different forms of discrepancies emerged during the comparison of SPH-emulator results.Appendix 1 shows more details about the results obtained with the three resolutions.However, in most practical applications, this level of comparison can be disregarded, as the discrepancies can be considered compatible with simulation tolerances.On the other hand, a robust matching can be observed in terms of major flow features, like the initial advancement of the fluid front, and the interaction with the walls, including the formation of a swirl in proximity of the left wall, visible at the bottom of Figure 2.
To quantitatively compare these results, Figure 3 shows the position of the waterfront over time in the two cases, measured as the position of the rightmost particle of fluid over time.Taking into account the scale of the fluid domain and the adopted spatial resolution, the two fronts advance with a strong matching, with a little divergence developing after 2 s of evolution.By looking at the particles configurations in Figure 4, we can see that the discrepancy in the advancement of the front is generated by a slightly different particle configuration at the front of the flows.The difference is compatible with the scale of a particle size and, therefore, is comparable to the simulation tolerance.On the other hand, a very strong agreement can be observed for the overall profile of the flow.
For the emulator to be used in a scientific context, it is also important to assess the ability to generalize and correctly predict the behavior of the fluid in conditions that have not been presented during training.To this aim, we simulate two other problems that involve the same physics of the dam break but introduce different levels of complexity.We study a flow at lower complexity by examining the behavior of a still volume of fluid in a tank.This is a first generalization test of the emulator, which, being trained over fluid in motion, should also be able to deal with fluid at rest.We used a 90 m by 15 m domain of fluid, with the same physical and simulation properties of the training fluid, and we used the emulator previously trained over the dam break to reproduce 200.0 s of evolution.We then simulated this problem with the original equation-based SPH model and compared the results.Although in this case the individual particles assume a different configuration, the overall behavior of the fluid was consistent with the expectations for the test case and with the respective SPH-based simulation.
To test the ability of the emulator to deal with flows with higher complexity than the training case, we include the presence of waterfalls, obtained by starting the flow from an elevated shell, again using the model trained over the dam break.Figure 5 shows the evolution of the flow, simulated with both the emulator and the reference SPH model.Similar to what observed for the dam break, some differences arise over simulated time, due to the discrepancies between the computed and estimated accelerations that gather during the simulation.However, the sequence of major flow dynamics is preserved.
As a remark about generalization capabilities, we recall the peculiarity of using an emulator that reproduces particle interaction quantities rather than directly modeling the spatial behavior of the fluid.In this case, the emulator learns to reproduce the numerical relationships between the variables in the reference model.Therefore, when the emulator is presented with a different problem to simulate than what was seen during training, as long as the values of the variables are in ranges where they hold the same physics learnt from the reference simulation, the emulator should be able to produce physically meaningful results.From the generalization tests we saw previously, even though some differences appeared in particle distribution with respect to the reference simulations, the fluid consistently exhibited a realistic behavior across all examined cases.
In addition, for the generalization tests, we reproduce the experiments (SPH simulation, training, and emulation) with the two different spatial resolutions ( = p Δ 1.0 m and = p Δ 0.5 m), obtaining results consistent with the respec- tive SPH-based simulations.The only limitation that we have encountered occurs when increasing flow complexity and spatial resolution, with some instabilities happening at the impact of jets or drops.However, these limitations are inherently conservative, as they primarily influence the numerical execution of the simulation and do not lead to physically unrealistic results.We have noted that training  the network with input SPH quantities sampled during the predictor phase of the integration (see Section 2.2) and SPH reference outputs sampled during the corrector phase (see Section 2.2), improves the stability of the emulator.This could be explained considering the similarity of this approach to an implicit integration scheme.Although this procedure introduces a deviation from the reference integration scheme, we deem this deviation to be negligible with respect to the achieved benefit in stability.

Conclusions
We presented an AI-CFD emulator based on the use of an ANN to emulate the SPH momentum equation and produce forces between particles starting from quantities representative of the particle state.
This emulator is meant to take into account the physical meaningfulness of the results, so it was developed by referring to SPH formulations adopted in real-world applications.We, therefore, adopt a second-order predictor-corrector integration scheme, where the density is integrated by means of a continuity equation.We adopted the dynamic boundary and studied the effect of their presence within the training set, finding that the results of the emulator were more disordered when not including the boundary particles data.Moreover, we formed all the input features from pair-related quantities, enforcing reciprocity in the interaction, which is a requirement for energy conservation.
Simulations of a dam break problem showed that the results of the emulator are faithful to those of the equationbased simulation.To validate the model, from a qualitative point of view, we analyzed the presence of major flow features in the generated results, and from a quantitative point of view, we studied the position over time of a dam break front, finding good levels of matching between the results.
We then studied the ability of the emulator to generalize for conditions never seen during training phase.We used the model trained over dam break data to different settings with more or less complexity with respect to the reference case.We have compared the emulations with the corresponding simulations, showing similar levels of agreement to what obtained for the validation.In addition, we also have changed the spatial resolution of the problems and repeated the experiments, to verify the robustness of this approach also for changes of this type.Even if particles distribution is not reproduced exactly, due to the differences in the accelerations estimation, these tests show how this emulator is able to reproduce the main features of the fluid, with a faithful appearance of jets, vortexes, and droplets.
Therefore, the emulator here presented is capable of faithfully reproducing traditional SPH simulations and generalizing over test cases with varying levels of complexity.
conditions.The fluid that we model has density = ∕ ρ 1.00 kg m 3 , artificial speed of sound = c 560.29 m/s 0 , and artificial viscosity coefficient = α 0.01.The tank is 90 m long, and the fluid is initially confined within a rectangular shape of width 30 m and height 40 m, as shown in Figure1.We will consider three spatial resolutions, for which the whole particle set is discretized by using a spatial step of = p m, resulting in 6,681, 3,423, and 2,151 particles, respectively.The initial density is set by using an inversion of the state equation (Eq.(

Figure 1 :
Figure 1: Dam break test case used to generate training data.On the left is the initial simulation setup with particles colored by density [ ∕ kg m 3 ].On the right is the final frame of the dataset with particles colored by velocity magnitude [m/s].

Figure 2 :
Figure 2: Simulation (left) and emulation (right) of a dam break with an artificial viscous term.The time frames per line are related to 1.5, 8.5, and 15.0 s of simulation.Particles are colored by velocity magnitude (m/s), with a minimum of 0.0 m/s and a maximum of 29.0 m/s.

Figure 5 :
Figure 5: Simulation (left) and emulation (right) of the waterfalls problem with artificial viscous term.The time-frames per line are related to 7.5, 22.0, and 30.0 s of simulation.Particles are colored by velocity magnitude (m/s), with a minimum of 0.0 m/s and a maximum of 29.0 m/s .