{'fraction_non_alphanumeric': 0.030013022805051637, 'fraction_numerical': 0.009570247433295981, 'mean_word_length': 4.583361515049036, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Structure Preservation for the Deep Neural Network Multigrid Solver

The simulation of partial differential equations is a central subject of numerical analysis and an indispensable tool in science, engineering and related fields. Existing approaches, such as finite elements, provide (highly) efficient tools but deep neural network-based techniques emerged in the last few years as an alternative with very promising results. We investigate the combination of both approaches for the approximation of the Navier-Stokes equations and to what extent structural properties such as divergence freedom can and should be respected. Our work is based on DNN-MG, a deep neural network multigrid technique, that we introduced recently and which uses a neural network to represent fine grid fluctuations not resolved by a geometric multigrid finite element solver. Although DNN-MG provides solutions with very good accuracy and is computationally highly efficient, we noticed that the neural network-based corrections substantially violate the divergence freedom of the velocity vector field. In this contribution, we discuss these findings and analyze three approaches to address the problem: a penalty term to encourage divergence freedom of the network output; a penalty term for the corrected velocity field; and a network that learns the stream function, i.e. the scalar potential of the divergence free velocity vector field and which hence yields by construction divergence free corrections. Our experimental results show that the third approach based on the stream function outperforms the other two and not only improves the divergence freedom but in particular also the overall fidelity of the simulation.


Introduction
The tremendous progress on deep learning in the last decade for applications such as image recognition and machine translation [26] has motivated their use also for the simulation of physical systems. Despite promising results that have already been obtained on this, many fundamental questions remain open. One of them is how physical constraints, such as energy and momentum conservation or divergence freedom, can be integrated into a neural network-based simulation-and if an explicit consideration of this is necessary or meaningful.
In the literature, two principle approaches to enforce physical constraints have been proposed. In the first one, they are imposed by construction, e.g. through an appropriate network design or by representing the physical quantities so that their properties and relationship have to be satisfied. Examples are symplectic neural network [36,12] that impose symplecticity by construction or the representation of fluid velocity fields in a divergence free basis with learning being performed on the basis function coefficients [47,52]. The second and more common approach to physical constraints in neural networks is to encourage that these hold but not strictly enforce them. Such a "soft" formulation is typically accomplished using a penalty term in the loss function but can also be realized, e.g., with a Bayesian formulation with the the physical constraints (including an entire partial differential equation) as priors. The approach also naturally accommodates realworld data where properties such as energy conservation almost never perfectly hold (e.g. because a system is not closed). However, to what extent a constraint is satisfied after training is difficult to control through a penalty term and its addition to the loss function increases the training complexity and can hinder convergence.
In this work, we explore how divergence freedom of the fluid velocity vector field can be ensured when a neural network is used in the simulation of the incompressible Navier-Stokes equations. In particular, we build on the recently introduced Deep Neural Network Multigrid Solver (DNN-MG) [34] that tightly integrates a neural network into a multigrid finite element solver to replace the computations on one or multiple finest mesh levels by a neural network-based prediction. In [34], the efficiency of DNN-MG for the simulation of the incompressible, instationary Navier-Stokes equations was demonstrated. However, the predicted velocity vector field is not divergence free. We hence study in this paper three different approaches to encourage or enforce divergence freedom within the DNN-multigrid formulation. First, we modify the network's loss function and include a penalty term to limit the divergence of the predicted velocity field. Second, we penalize the divergence of the corrected velocity field, which has the potential to also reduce the divergence residual of the finite element velocity field. And third, we design a neural network that learns the stream function, i.e. the potential associated with the velocity vector field, and hence provides by construction divergence free corrections.
Our results demonstrate that encouraging or enforcing divergence freedom improves the overall quality of the simulations, measured, e.g., by the drag and lift functionals, although the three approaches we developed lead only to rather modest reductions in the divergence. For the two approaches based on a penalty term, this results from a delicate balance between prediction accuracy and divergence freedom that needs to be attained but that it is difficult to satisfy in practice. The approach based on the stream function provides overall the best results and achieves lift and drag functionals close to a fine reference solution. At the same time, the corresponding network is more efficient to train. We attribute this improved performance to a formulation that better respects the intrinsic structure of the partial differential equations.
The remainder of the paper is structured as follows. After discussing related work, we recall the DNN-MG solver in Sec. 3. In Sec. 4 we describe our approaches to ensure divergence freedom within DNN-MG and present and discuss our numerical results. We conclude in Sec. 5 where we also provide an outlook for possible directions for future work.

Related Work
Recently, there has been an increasing interest to employ (deep) neural networks for the simulation of physical systems. In the following, we will discuss existing approaches with an emphasis on those that consider partial differential equations and that are most pertinent to our own work.
Neural network-based techniques for partial differential equations can, broadly speaking, be classified into three categories. The first one, e.g. the DeepRitz method by E and Yu [16] or the physics-informed neural networks by Raissi, Karniadakis and coworkers [41,32], use neural networks in place of a basis or frame to represent solutions. Training then amounts to fitting the network representation to the initial or boundary value problem. This has hence to be repeated for each problem instance but the approach yields typically a mesh free solution that scales well to high dimensions [16] and the methodology applies to a wide range of partial differential equations. The second approach takes inspiration from classical time stepping schemes and considers the solution of the PDE as the problem of predicting the next state from the previous ones. The well developed network architectures for sequential data can hence be employed, e.g. recurrent neural networks (RNN) with Long-Short-Term-Memory units (LSTM) [22] or Gated Recurrent Units (GRUs) [13], or temporal convolutional networks, cf. [1]. Examples of this approach are [23,49,53] and the work by Han, Jentzen, and E [20] who use a formulation of semi-linear, parabolic PDEs as backward stochastic differential equations and then realize each time step of these as a layer in a deep neural network. Since we combine a classical finite element based simulation, which relies on a time stepping scheme, with a neural network, our work also falls into this category. The third class of neural network-based techniques for partial differential equations uses these to represent solution operators. They hence use the initial condition as well as parameters on which the solution depends as input to the neural network and it then maps these to the solution. For example, in [7] the solution operator for parametric partial differential equations is learned based on a discretization obtained via model reduction and [28,27] learn mappings between infinite dimensional spaces that can be solution operators. Lu, Jin, Karniadakis [31] learn nonlinear operators from data, based on the result that any such operator can be approximated by a neural network [11].
For flow problems, Raissi et al. [42] demonstrated that the velocity and pressure of the Navier-Stokes equations can be learned to good approximation only from observing passively advected tracer particles. Nabian and Meidani [38], Yang and Perdikaris [54], and Raissi and Karniadakis [40] exploited that not only observations are available but also a known analytic model. These authors hence also use the deviation of the network prediction from the model as an additional penalty term in the training loss. Kasim et al. [25] recently demonstrated neural network-based simulations for a broad range of applications, including fluid dynamics, by also optimizing the network architecture itself during training. Eichinger, Heinlein and Klawonn [33] use convolutional neural networks and techniques from image processing to learn flow patterns for the Navier-Stokes flow around objects of different shape.
Next to the above direct approaches to the simulation of partial differential equations using neural networks, there have been different attempts to integrate these into existing numerical formulations. As in our work, the objective is then to combine the benefits of neural networks and classical methods and obtain techniques that would be difficult or impossible with either approach alone. For elasticity, neural network-based representations for constitutive relations were learned in [50,6,44] and then used in classical simulations, e.g. based on finite elements. Wiewel et al. [53] presented a simulation of the Navier-Stokes equation where an LSTM-based deep neural network is used to predict the pressure correction within a classical finite element simulation. In [30] a convolutional neural network is interpreted as a finite difference discretization of a PDE and this was also proposed in [45]. Bar-Sinai, Hoyer, Hickeney and Brenner [2] learn general sub-grid scale models, used, e.g., with finite difference or finite volume discretization. In the context of large eddy-type simulations, the learning of closure models has also been considered, e.g. in [29]. Wan, Dodov, Lessig, Dijkstra, and Sapsis [52] combine a wavelet-based discretization and perform learning in the coefficient space where spacefrequency correlations are resolved but, so far, without simulation. Recently, Stevens and Colonius [49] enriched finite difference and finite volume simulations with a deep neural network to resolve sharp discontinuities (as in Burger's equations) and to deal with PDEs with chaotic behavior.
To ensure that neural network-based predictions respect physical invariants, such as energy conservation or divergence freedom, network architectures tailored towards physical simulations have been proposed. Next to an implicit enforcement as in Physics Informed Neural Network (PINN) [32], various authors proposed architecture where desired properties hold by construction. Finizi, Wang and Wilson [18] and Cranmer et al. [15], e.g., learn the dynamics of classical mechanical systems by directly learning the Hamiltonian or Lagrangian functions. In related work, [36,12,24] developed neural networks that incorporate the symplectic structure of Hamiltonian mechanics and they demonstrate that this improves generalization and prediction accuracy. Based on the work by Bruna and Mallat [10], Bietti and Mairal [8] proposed neural network architectures that are invariant under group transformation, including diffeomorphisms. A neural network that is directly defined on a Lie group, and hence also enjoys group equivariance, was proposed by Bekker [5].

Deep Neural Network Multigrid Solver
In this section we summarize the Deep Neural Network Multigrid Solver (DNN-MG) introduced in [34]. For a detailed description the reader is referred to the original paper.

Finite Element Discretization of the Incompressible Navier-Stokes equations
We want to solve the incompressible, instationary Navier-Stokes equations given by where Reynolds number, and f an external force. The initial and boundary conditions are where n denotes the outward facing unit normal on the boundary ∂Ω of the domain. On the outflow boundary Γ N we consider the do-nothing outflow condition [21], which is well established to model artificial boundaries; see the discussion in [9] on possible shortcomings and a variation that is able to guarantee uniqueness of a solution. We discretize Eq. 1 using a weak finite element formulation with where φ h and ξ h are test functions and W (r) h is the space of continuous functions which are polynomials of degree r on each mesh element T ∈ Ω h in the mesh domain Ω h . The resulting equal order finite element pair V h × L h does not fulfill the inf-sup condition. We hence use stabilization terms of local projection type [3] with parameter α T = α 0 · Re · h 2 T and projection π h : W into the space of linear polynomials.
With the second order Crank-Nicolson method for time discretization, the solution of Eq. 1 subject to Eq. 2 at each time step amounts to determining the state x n = (v n , p n ) such that for all test functions φ i h and ξ i h . Eq. 3 is a large nonlinear system of algebraic equations that is solved by Newton's method based on the initial guess x for l = 1, 2, . . . Here we denote by A (x (l−1) ) the Jacobian of A at x (l−1) , which for the problem at hand can be computed analytically, cf. [43,Sec. 4.4.2]. Each Newton step requires the solution of the linear system in Eq. 5, where the system matrix A (x (l−1) ) is sparse but non-symmetric and not definite due to the saddle point structure of the underlying Navier-Stokes equation. To approximate the solution to Eq. 5 with optimal robustness, we employ the generalized minimal residual method (GMRES) introduced by Saad [46]. The convergence is accelerated through a preconditioner that is realized by a single sweep of a geometric multigrid solver.
The geometric multigrid method is based on a hierarchical approximation of a linear system on a sequence of finite element spaces V 0 ⊂ V 1 ⊂ · · · ⊂ V L defined over a hierarchy of meshes Ω 0 , . . . , Ω L = Ω h . Instead of treating the linear system on its domain, given by the finest mesh level L (as one would do in a traditional solver) the high frequencies there are smoothed and the remaining errors are treated on lower levels. This is accomplished by computing the residual on level L and restricting it to the next coarse level L − 1. There, the process is repeated and this is performed recursively until a coarsest level has been reached where a direct solver is employed. The updates from the coarser levels are then prolongated back to level L and the coarse-to-fine and fine-tocoarse iteration is repeated until a prescribed error tolerance has been met. The mesh transfer from fine to coarse is accomplished with L 2 -projections, known as restrictions, and from coarse to fine with interpolations, known as prolongations. The smoothing that preceded the restrictions is realized with a simple smoothing operator S(A l , b l , x l ) that yields an approximate solution of the linear system l b l and that aims to quickly reduce all high frequency components of the residual b l − A l x l . In our implementation we use a simple iteration of Vanka type [51], which allows for easy parallelization and gives very good performance with less than 5 pre-and post-smoothing steps [17]. while not converged do Newton's method in Eq. 5 3: Set up rhs of Eq. 3 for next time step Restriction of rhs to level L 10: end for

Deep Neural Network Multigrid Solver
The solution of the Newton iteration in Eq. 5 using GMRES and the geometric multi-grid method is highly efficient. However, since one has multiple GMRES steps for every step of the Newton method and for each one also one up-and down-sweep of the multigrid method, there is still a substantial amount of computations. The majority is thereby required for the finest mesh level (typically about three quarters in the multigrid method) so that a refinement of the mesh, or, equivalently, the use of an additional mesh level L+1, incurs a substantial increase in the computational costs. The Deep Neural Network Multigrid Solver (DNN-MG) addresses this bottleneck and replaces computations on one or multiple finest mesh levels with a neural network-based prediction. An overview of DNN-MG is provided in Algo. 1.
To ensure that the DNN-MG solver is efficient and flexible, its neural network is build around a few critical design principles: 1. The neural network operates patch-wise, i.e. locally on small neighborhoods of the mesh domain Ω L+k . This ensures the generalizability of DNN-MG to different flow regimes and meshes.
2. The neural network has memory. This ensures that complex flow behavior can be predicted and that corrections are coherent in time.
3. The neural network uses the residual of the nonlinear problem in Eq. 3 on level L + k as input. This provides rich information about the required correction.
Together, these design principles lead to a surprisingly compact neural network architecture that at the same time generalizes well to different flow regimes and domains.
In particular, the implementation in [34] has just 8634 trainable parameters with an architecture based on Gated Recurrent Units. It can hence be trained in a few hours with a small number of example flows as training set.
As shown in Algo. 1, the neural network-based correction of DNN-MG is applied at the end of every time step after the Newton solver has computed an updated velocity v L n on level L, cf. Eq. 5. To obtain the correction, v L n is first prolongated to level L + 1, yieldingṽ L+1 n = P(v L n ). The neural network is then employed individually for each patch P i (a mesh element on level L in the most simple case as in [34]) and predicts a velocity correction d L+1 n,i . The inputs to the neural network are thereby also entirely local and include the residual of Eq. 3 on level L + 1, the prolongated velocityṽ L+1 n , geometric properties such as the patch's aspect ratio, and possibly other quantities such as the Péclet number over the patch. With the correction d L+1 n , a provisional right hand side b L+1 n+1 = Rhs(ṽ L+1 n + d L+1 n , f n , f n+1 ) of Eq. 5 is formed on level L + 1 and then restricted to level L, i.e. we compute b L n+1 = R(b L+1 n+1 ). Restricting the right hand side instead of the corrected velocityṽ L+1 n + d L+1 n is thereby essential to propagate information back to level L since, ideally, R(P(v L n ) + d L+1 n ) = v L n . The corrected b L n+1 is then used in the next time step in the Newton solve, which is again improved by a neural network-based correction at the end of the time step.
The neural network of DNN-MG is trained using a high fidelity finite element solution obtained on the fine mesh level L + 1 and with the loss function Here, N is the number of time-steps that make up the training data and, since DNN-MG operates strictly local over patches, for each one the loss accumulates the local residuals in the second sum. Due to this strictly local structure of DNN-MG, a very large data set is generated by a single resolved simulation. In [34] it was demonstrated that the DNN-MG algorithm substantially improves lift and drag functionals for the classical channel flow around a cylinder and generalized well to flow regimes not seen during training, e.g. a flow without an obstacle or a flow in an L-shaped domain.

Enforcing Divergence Freedom in DNN-MG
In our first study [34], we have shown that DNN-MG is able to increase the approximation quality at very little extra cost. The error in the velocity vector field but also in derived quantities such as the drag and lift functionals of the flow around an obstacle are reduced considerably with an increase in the computation time that is orders of magnitude below those for adding an additional mesh level. We have, however, also observed that the network is not good at ensuring a solenoidal velocity vector field. Although the lossfunction aims at minimizing the distance to the high fidelity solution, which weakly satisfies divergence freeness to a certain degree, the neural network minimizes Eq. 6 with a vector field with a substantial divergence.

Violation of Divergence Freedom in Existing Formulation
The quadratic equal-order finite element approach (v h , p h ) ∈ V h × L h that provides the basis for the DNN-MG algorithm in our current implementation is stabilized with local projections [3], cf. Eq. 4. The discrete divergence equation thus reads and since ∇ · V h ⊂ L h only a perturbed divergence condition is satisfied

Approaches for Enforcing Divergence Freedom
In the following we will discuss three approaches for enforcing divergence freedom. Numerical results will be presented afterwards in Sec. 4.3.

A Penalty Approach for Divergence Free Learning
The first approach to enforce divergence freedom of the neural network output of the DNN-MG method is through a modification of the loss function by a squared divergence as penalty term. One then has where is the original loss function in Eq. 6 and d L+1 = N (v L ,ṽ L+1 ) is the correction by the network based on the coarse solution v L n and its prolongatioñ v L+1 . The parameter γ > 0 is used to control the weighting of the loss function. Eq. 8 penalizes the divergence of the correction and, through this, the neural network should at least not further distort the divergence of the velocity field.
An obvious extension of this idea is to penalize the divergence of the corrected solution v L +d L+1 . This is, in principle, advantageous since it could also correct for the divergence residual in the finite element formulation, cf. Eq. 7. We thus introduce the second, modified loss function

A Strictly Divergence Free Network Architecture
Instead of including a penalty term in the loss function to ensure divergence freedom in a weak sense, our third approach modifies the neural network architecture so that only solenoidal vector fields are returned. We accomplish this by learning the stream function ψ, which can be seen as a scalar potential for the divergence free velocity vector field [14], For an easy and efficient integration into the finite element framework, we are looking for divergence free velocity corrections d L+1 by the neural network that are in the velocity h . This holds when we learn stream functions ψ h ∈ L L+1 h that are in the pressure space L L+1 h and the construction is then also again local on each patch.
For simplicity we assume that a patch P i corresponds exactly to one element on the fine mesh Ω L+1 and that it is a quadrilateral parallel to the axes on an integer lattice.
The 9 local scalar basis functions of the second order finite elements can then be defined in the simple monomial basis, and a stream function ψ h ∈ L L+1 h is described by the 8 basis function coefficients s P i = (s 2 , . . . , s 9 ) for this basis (the stream function has no harmonic part and ξ 0,0 can thus be disregarded). Since we work with L h = W (2) h , the velocity v ψ h associated with ψ h lies in L h | P , as desired, and it can be represented in the eight velocity basis functions derived from the ξ α,β (x, y) = x α y β , and the basis function coefficients for the stream function and the derived velocity are identical.
To be able to work within the usual finite element theory, we relate the non-standard velocity basis functions v ψ i to pointwise values using the Lagrange basis that satisfies for the 9 Lagrange points δ ij . Each of the derived velocity functions v ψ i has a unique representation in this basis with coefficients η x ij , η y ij for i = 2, . . . , 9 and j = 1, . . . , 9. The pointwise values of the velocity vector field can thus be reconstructed as v ψ  Figure 2: Geometry of the training scenario with a parabolic inflow profile Γ in , donothing boundary conditions at the outflow boundary Γ out and no-slip conditions on the walls Γ wall . The center of the obstacle is at (0.3, 0.15). For the test scenario the obstacle is shifted such that the center is at (0.3, 0.25).
With the above, we can train for the eight stream function coefficients s P i = (s 2 , . . . , s 9 ) that through Eq. 11 and Eq. 13 correspond to a divergence free velocity correction

Numerical Evaluation
To evaluate the efficacy of the three approaches developed in Sec. 4.2 to improve the divergence freedom of the velocity vector field, we implemented these using PyTorch [39] and Gascoigne 3D [4] based on the framework developed in [34]. For the penalty-based approaches with the loss functions L p 1 and L p 2 the network and other aspects of the training process remain unchanged compared to [34] where the network had 8634 parameters. The network architecture we use for the stream function-based approach with L ψ is also similar to the one in [34] but the convolutional layers are replaced by fully connected ones of size 8 × 8. This results in a network with 9000 parameters.
The flow we considered in our experiments is a modification of the classical benchmark of a flow around a cylinder [48] where the obstacle is now a square, see Fig. 2. The modification simplified the implementation, since we were able to use uniform meshes, but does not change the qualitative behavior of interest to us and in particular one still has a non-stationary wake behind the obstacle as in the original benchmark [48]. The singularities at the corners of the obstacles have been excluded in the evaluation below.
Divergence freedom We first compare the divergence ∇ · v(t) L 2 (Ω L ) obtained with the three approaches (DNN-MG-L i ) as well as a coarse multigrid solution (MG), the original DNN-MG method proposed in [34] (DNN-MG), and, as reference, a fine multigrid solution on level L + 1 (MG-fine). The results are presented in Fig. 1. They show that the approaches from Sec. 4.2 are not able to significantly reduce the divergence. While DNN-MG-L ψ based on the strictly divergence free output performs best, the reduction is still modest. Ultimately, however, this is not surprising since and the perfectly divergence free correction does not address the divergence residual in v L . The second approach, using L p 2 , was designed to address this but no reduction could be observed in practice. We investigated the disappointing performance of the modified loss functions L p 1 and L p 2 that included a penalty term and found that it stems from an unstable balance between the minimization of the divergence and the fulfillment of the finite element Galerkin equation: If the weighting factor γ for the divergence is chosen close to zero then the network behaves like the classical DNN-MG approach. If we choose larger values for γ > 0, however, then the network is not able to sufficiently minimize the discrepancy to the high fidelity Galerkin solution and overall the approximation quality suffers. We tested γ in the range 10 −1 , · · · , 10 −10 but did not find a value where the two error terms in the loss functions L p 1 and L p 2 were sufficiently balanced. For γ > 10 −7 we, in fact, observed that the penalty term causes the output to converge towards zero, which is a divergence free velocity vector field but not the desired correction. The effect was most pronounced for DNN-MG-L p 2 . Fig. 3 shows the drag and lift functionals of the obstacle for the different approaches. In Table 1 we also show the corresponding minimum and maximum values, the mean values, and the amplitudes of the oscillation as well as the frequencies of the periodic flow pattern. In the lower table of Table 1 we report the relative errors with respect to the high fidelity solution MG (fine).

Drag and Lift functionals
The results show that the non-smooth corners of the square obstacle have a negative effect on the performance of the original DNN-MG approach since it no longer provides the improvement in the functional values observed in [34] for a smooth obstacle. Among the approaches proposed in the present work, DNN-MG-L ψ based on learning the stream function provides the most promising results. It is able to reduce the relative error in all functional outputs and almost matches the values of MG (fine) in Fig. 3. Furthermore, the training time for the neural network is in this case only ≈ 75% in comparison to the other DNN-MG approaches.
For DNN-MG-L ψ we observe a deviation of the oscillation frequency that is larger than for the other approaches. We believe that this defect stems from a complicated interplay of temporal dissipation and spatial discretization. For a discussion we refer to [34] and also to [35] where the same effect is observed for a parallel time stepping scheme. In the context of the deep neural network multigrid approach it must be investigated if this problem can be resolved by training with high fidelity solutions that come from discretizations that are finer in space and in time, i.e. by working with a finer time step. the fine reference solution MG (fine). To obtain meaningful comparisons, the outputs are synchronized such that we consider snapshots starting from a maximum of the drag functional, cf. Fig. 3. The results verify that all approaches proposed in Sec. 4.2 are able to reduce the error in the velocity field compared to MG and DNN-MG.

Discussion
The results presented in this section are surprising in that the loss functions of Sec. 4.2 designed to reduce the divergence improve the overall quality of the obtained results, e.g. in terms of the lift and drag functionals, but do not lead to significantly better divergence freedom. More work is hence required to obtain hybrid finite element / neural network-based simulations that improve this important property. For the loss functions with the penalty factor it should in particular be reconsidered how these can be formulated so that their two constitutive terms do not compete with each other as we observed in our experiments. The efficiency of DNN-MG-L ψ based on the stream function is not entirely surprising to us. The network there only has to learn a scalar function instead of two components of a vector field that, furthermore, satisfy a nontrivial constraint through the divergence freedom. The stream function formulation also encodes the structure of the dynamics in a different way than the usual velocity one, cf. [37], and this leads often to simulations  that preserve structural properties such as energy or vorticity. In the finite element context this was already observed by Fix [19]. We believe that this partly also explains the reduced training time for DNN-MG-L ψ .

Conclusion
We have presented three approaches to enforce divergence freedom within the hybrid finite element / neural network DNN-MG method [34] for the simulation of the Navier-Stokes equations. The first two of these use a penalty term that encourages divergence freedom while the third one guarantees a divergence free correction by learning the stream function. The penalty term-based approaches did not provide significant improvements compared to DNN-MG for both the divergence and the drag and lift functionals. The third approach based on the stream function yielded significantly better overall results but also only a rather modest improvement for the divergence. More work is hence required in the future to obtain hybrid finite element / neural network simulations that respect the divergence freedom of the Navier-Stokes equations. For the penalty-based approaches it should be investigated how the two terms in the  loss function can be balanced. DNN-MG-L 2 is thereby particularly interesting since it allows, in principle, to provide an overall reduction by correction v L towards a lower divergence. For DNN-MG-L ψ it would be interesting to investigate the efficiency when also the velocity-stream function formulation is used for the simulation.
In future work, we also want to investigate the considerably more challenging problem of a hybrid finite element / neural network simulation for the Navier-Stokes equations in 3D. How other physical invariants, e.g. energy or Kelvin's circulation theorem, can be incorporated into hybrid simulations and what benefits this provides should also be investigated in more detail. Stability, consistency and convergence are, furthermore, pertinent open questions that should be addressed. schaft (DFG, German Research Foundation) -Project-ID 422037413 -TRR 287.