DNN-MG: A Hybrid Neural Network/Finite Element Method with Applications
  to 3D Simulations of the Navier-Stokes Equations
We extend and analyze the deep neural network multigrid solver (DNN-MG) for
the Navier-Stokes equations in three dimensions. The idea of the method is to
augment a finite element simulation on coarse grids with fine scale information
obtained using deep neural networks.
  The neural network operates locally on small patches of grid elements. The
local approach proves to be highly efficient, since the network can be kept
(relatively) small and since it can be applied in parallel on all grid patches.
However, the main advantage of the local approach is the inherent
generalizability of the method. Since the network only processes data of small
sub-areas, it never ``sees'' the global problem and thus does not learn false
biases.
  We describe the method with a focus on the interplay between the finite
element method and deep neural networks. Further, we demonstrate with numerical
examples the excellent efficiency of the hybrid approach, which allows us to
achieve very high accuracy with a coarse grid and thus reduce the computation
time by orders of magnitude.
Nils Margenberg, Robert Jendersie, Christian Lessig, Thomas Richter
DNN-MG: A Hybrid Neural Network/Finite Element Method with Applications to 3D Simulations of the Navier-Stokes Equations

###### Abstract

We extend and analyze the deep neural network multigrid solver (DNN-MG) for the Navier-Stokes equations in three dimensions. The idea of the method is to augment of finite element simulations on coarse grids with fine scale information obtained using deep neural networks.

This network operates locally on small patches of grid elements. The local approach proves to be highly efficient, since the network can be kept (relatively) small and since it can be applied in parallel on all grid patches. However, the main advantage of the local approach is the inherent good generalizability of the method. Since the network is only ever trained on small sub-areas, it never "sees" the global problem and thus does not learn a false bias.

We describe the method with a focus on the interplay between finite element method and deep neural networks. Further, we demonstrate with numerical examples the excellent efficiency of the hybrid approach, which allows us to achieve very high accuracies on coarse grids and thus reduce the computation time by orders of magnitude.

## 1 Introduction

Accurate flow simulations remain a challenging task. The success of deep neural networks in machine translation, computer vision, and many other fields on the other hand, has lead to a growing (and renewed) interest to apply neural networks to problems in computational science and engineering. The combination of classical finite element approximation techniques with deep neural networks adds new aspects to the pure numerics-oriented approaches. At higher Reynolds numbers, accurate simulations, especially in 3D become increasingly difficult, and the classical methods reach their limits. Although the finite element method is highly efficient and established for the discretization of the Navier-Stokes equations, fundamental problems, such as the resolution of fine structures or a correct information transport between scales, are still not sufficiently solved. Even if linear solvers with optimal complexity (such as
multigrid methods) are used, nonlinear equations such as the Navier-Stokes equations show more and more local features under mesh refinement and these turn out to be a hurdle for the approximation.

We usually consider finite elements as a hierarchical method, where natural hierarchies of finite element meshes and spaces can be constructed (adaptively) to yield approximations of increasing accuracy. This hierarchical setup is then used within multigrid methods for efficient solution of the algebraic system. Here, we connect such a hierarchic finite element approach with neural networks in a hybrid setting: while coarse meshes are handled in the traditional finite element multigrid way, updates in finer meshes are learned by the neural network. They are hence used whenever a full resolution of the effects does not seem possible or efficient.

We label this approach the Deep Neural Network Multigrid Solver (DNN-MG) as it is based on hierarchies of meshes and functions spaces and combines the tools of multigrid solvers with deep neural networks (cf. Figure 1). However, the neural networks are used to predict updates to the nonlinear problem and hence, the approach can be used for an upsampling of any kind of finite element (also of finite volume or finite difference) solutions towards a representation on a finer mesh.

We demonstrate the efficiency, generalizability, and scalability of our proposed approach using 3D simulations. We first train the neural network on data from the classical flow with one circular obstacle. To analyze the generalization capability of DNN-MG, we test the network on channel flows with one or two obstacles at different Reynolds numbers. The obstacles have an elliptical cross-section with varying eccentricities. The obtained solutions as well as the lift and drag functionals demonstrate that DNN-MG obtains considerably better accuracy than a coarse multigrid solution on mesh-level \(L\) while taking less than 3% of the computation time required by a full solution on \(L+2\) levels. Therefore, DNN-MG offers a speedup by a factor of 35. DNN-MG's efficiency is evident, as it requires only double the time of a coarse mesh solution on level \(L\), yet offers substantial improvements in the overall solution quality.

The paper is organized as follows. In the Section 2, we review related work on using neural networks for the simulation of partial differential equations. Sections 3 provides a recap of the Navier-Stokes equations' solution using the geometric multigrid method, while Section 4 introduces the neural networks used later in our numerical experiments. In Section 5, we present the deep neural network multigrid solver, discussing its design in a general form, which makes it applicable to other problems. Finally, we present our numerical results in Section 6.

Figure 1: The idea of DNN-MG is to use a multigrid solver to solve the problem on the first \(L\) levels of the mesh hierarchy. The obtained solution is then prolongated to a finer level \(L+J\), where a neural network predicts corrections to the solution. These corrections are incorporated into the time evolution through the time stepping and enter through the right-hand side.
## 2 Related works

We discuss different approaches to the simulation of physical systems using deep neural networks.

The most direct approach to utilize (deep) neural networks for the solution of partial differential equations is to represent the solution of the PDE as such a network. In [35] already, a least squares formulation has been employed to minimize the residual in the training process. The idea has recently been used again in the Deep Ritz Method [18] and a variety of similar approaches have emerged in the last years and coined Physics Informed Neural Network (PINNs) in the literature [38]. For a comprehensive review of PINN and related approaches, we refer to the comprehensive review [17]. Such approaches use an existing partial differential equation and can be considered data free since no explicit training data is required. They are, however, static in that the solution to one specific problem instance is learned. PINNs show superiority for instance in very high dimensional problems [40] but they do not reach the efficiency of established methods when standard problems, e. g. from elasticity or fluid dynamics, are considered (cf. [24] for case studies on the Poisson, Allen-Cahn and Schrodinger equation). The idea of PINNs has also been extended to learning solution operators [39, 14, 34]. The main advantage of this approach is that once the solution operator is trained, it can be applied to other settings. This enables the application to inverse problems [50, 30, 13] or optimal control [44]. Inspired by classical time stepping schemes, another approach to a neural network-based simulation of (evolutionary) partial differential equations is to consider these as time series and predict the next state of the system based on the previous ones. With this approach, network architecture for sequential data can be used, in particular recurrent neural networks and their extensions such as Long-Short-Term-Memory units (LSTM) [28] and Gated Recurrent Units (GRUs) [16] that we also used for the DNN-MG Navier-Stokes solver [45, 46].

At the other end of the spectrum are approaches that are purely data-based. For example, in [19] flow fields are learned in convolutional neural networks based on the flow geometry and using techniques from image processing. Recent work [53, 9, 36] uses purely data-driven neural network models trained on historical weather measurements and these perform on par with one of the most sophisticated weather forecasting models based on partial differential equations. The authors of [22] use transformer models for the prediction of dynamical systems. In the context of dynamical systems, the combination of Data Assimilation, Uncertainty Quantification, and Machine Learning techniques has gained significant interest, we refer to [15] for an in-depth review. For a recent overview on different neural network based approaches to approximate partial differential equations and related inverse problems we refer to [51].

One of the central open questions in the current literature is to what extent and how physical constraints or existing knowledge about the dynamics should be used for the training of a neural network, see also [23]. The Deep Neural Network Multigrid Solver (DNN-MG) is a hybrid approach that combined the model based finite element method with deep neural networks that also integrate model knowledge. DNN-MG is not a PINN since we do not represent the solution by only aim for fine scale updates. Also, the general setup of DNN-MG aims for generalization, and it is closer to learning a numerical solver than to learning a specific solution.

The close incorporation of the deep neural network into the finite element method, e. g. that we learn fine-scale solution coefficients, gives rise to analytical tools that can be used in first steps of error estimation [48].

In our first work on DNN-MG [46] we introduced a method which combines a geometric multigrid method with ANNs but only considered 2D simulations. We also adapted the method In 2D to ensure divergence-freedom by construction through a network that predicts the streamfunction [45]. In later studies we showed that the corrections added by the neural network improve the guesses in the Newton method, which leads to a reduction in wall-time compared to the coarse grid solutions [47]. Similar to the ideas of DNN-MG, the authors of [42, 43] use numerical solutions of higher order to train a neural network which generates corrective forcing terms. Similarly, in [20], low-order simulations
of wave equations are enhanced using neural networks to provide corrective forcing. This approach includes improvements in the temporal resolution of the simulations. Learning corrective forcing terms shares similarities with the approach of DNN-MG, as both methods introduce corrections to the solution through the right-hand side, which can be seen as corrective forcing terms. However, DNN-MG goes further by directly correcting the solution, potentially providing advantages in computing goal quantities like drag and lift. All these methods have a resemblance to reduced order models (ROMs). Along these lines, in [4], a ROM based on adaptive finite elements and correction terms combined with Artificial Neural Networks (ANNs) is presented. In [32] the authors developed a mesh based approach to ANN based solvers, leveraging existing FEM theory. The authors of [12] propose a neural framework for optimizing weak formulations in FEM. An ANN acting as control variable is included in the weak form and the minimization of a cost function yields desirable approximations where known data or stabilization mechanisms can be easily incorporated. Ainsworth and Dong propose a framework for adaptively generating basis functions, each represented by a neural network with a single hidden layer, which may be used in a standard Galerkin method [1]. In [49] the authors augment PDEs by ANNs to represent unknown terms. The hybrid model is then discretized using FEM.

The MgNet framework introduced in [25] draws connections from convolutional neural networks (CNNs) to Multigrid methods and thereby improve the design and understanding of CNNs. The authors of [41] propose a framework for unsupervised learning of Algebraic Multigrid (AMG) prolongation operators for linear systems which are defined directly on graphs. In the Multigrid framework, selecting an appropriate smoother is crucial to ensure convergence and efficiency, but the optimal choice of a smoother is problem dependent and is in general a major challenge for many applications. In [29], the authors optimize the smoothers in Multigrid methods, which they form by CNNs.

## 3 Numerical Methods

We solve model problems of incompressible viscous flow around a cylinder in 2D and 3D domains. This leads to the solution of the Navier-Stokes equations. Our approach involves a geometric multigrid method with a cell-based Vanka-smoother as a preconditioner to a Newton-Krylov method.

We consider a domain \(\Omega\in\mathbb{R}^{d}\) with \(d\in\{2,\,3\}\) and Lipschitz-continuous boundary \(\Gamma\) and a bounded time interval \([0,\,T]\). For solving the stationary Navier-Stokes equations we seek the velocity \(\mathbf{v}:\ [0,\,T]\times\Omega\to\mathbb{R}^{d}\) and pressure \(p:\ [0,\,T]\times\Omega\to\mathbb{R}\), such that

\[\partial_{t}\mathbf{v}+(\mathbf{v}\cdot\triangledown)\mathbf{v}- \nu\Delta\mathbf{v}+\nabla p =f \text{on }[0,\,T]\times\Omega \tag{1}\] \[\triangledown\cdot\mathbf{v} =0 \text{on }[0,\,T]\times\Omega,\]

where \(\nu>0\) is the kinematic viscosity and \(f\) an external force. The initial and boundary conditions are given by

\[\mathbf{v}(0,\,\cdot) =\mathbf{v}_{0}(\cdot) \text{on }\Omega \tag{2}\] \[\mathbf{v} =\mathbf{v}^{D} \text{on }[0,\,T]\times\Gamma^{D}\] \[\nu(\mathbf{n}\cdot\triangledown)\mathbf{v}-p\mathbf{n} =0 \text{in }[0,\,T]\times\Gamma^{N},\]

where \(\mathbf{n}\) denotes the outward facing unit normal vector on the boundary \(\partial\Omega\) of the domain. The boundary \(\Gamma=\Gamma^{D}\cup\Gamma^{N}\) is the union of \(\Gamma^{D}\) with Dirichlet boundary conditions and \(\Gamma^{N}\) with Neumann type conditions.

### Notation and variational formulation of the Navier-Stokes equations

By \(L^{2}(\Omega)\) we denote the space of square integrable functions on the domain \(\Omega\subset\mathbb{R}^{d}\) with scalar product \((\cdot,\cdot)\) and by \(H^{1}(\Omega)\) those \(L^{2}(\Omega)\) functions with weak first derivative in \(L^{2}(\Omega)\). The function spaces for the
velocity and pressure are then

\[\begin{split}\mathbf{V}&:=\mathbf{v}^{D}+H_{0}^{1}( \Omega;\Gamma^{D})^{d},\quad H_{0}^{1}(\Omega;\Gamma^{D})^{d}\,:=\left\{\,\mathbf{ v}\in H^{1}(\Omega)^{d}\,:\,\mathbf{v}=0\text{ on }\Gamma^{D}\right\}\\ L&:=\left\{\,p\in L^{2}(\Omega),\text{ and, if }\Gamma^{N}= \varnothing,\,\int_{\Omega}p\,\mathrm{d}x=0\right\},\end{split} \tag{3}\]

where \(\mathbf{v}^{D}\in H^{1}(\Omega)^{d}\) is an extension of the Dirichlet data on \(\Gamma^{D}\) into the domain. We normalize the pressure to yield uniqueness, if Dirichlet data is given on the whole boundary. With these spaces, the variational formulation of (1) is given by

\[\begin{split}(\partial_{t}\mathbf{v},\,\boldsymbol{\phi})+( \mathbf{v}\cdot\triangledown\mathbf{v},\,\boldsymbol{\phi})+\nu(\triangledown \mathbf{v},\,\triangledown\boldsymbol{\phi})-(p,\,\triangledown\boldsymbol{ \phi})&=(\mathbf{f},\,\boldsymbol{\phi})\quad\forall\boldsymbol{ \phi}\in H_{0}^{1}(\Omega;\Gamma^{D})^{d}\\ (\triangledown\mathbf{v},\,\xi)&=0\qquad\quad \forall\xi\in L\\ \mathbf{v}(0,\,\cdot)&=\mathbf{v}_{0}(\cdot)\quad \text{ on }\Omega.\end{split} \tag{4}\]

Let \(\Omega_{h}\) be a quadrilateral or hexahedral finite element mesh of the domain \(\Omega\) satisfying the usual requirements on the structural and form regularity such that the standard interpolation results hold, compare [54, Section 4.2]. By \(h_{T}\) we denote the diameter of an element \(T\in\Omega_{h}\) and by \(h\) the maximum diameter of all \(T\in\Omega_{h}\).

### Semidiscretization in space

For the finite element discretization of (4) we choose equal-order continuous finite elements of degree two for velocity and pressure on hexahedral meshes. By \(\mathbb{Q}_{h}^{(r)}\) we denote the space of continuous functions which are polynomials of maximum directional degree \(r\) on each element \(T\in\Omega_{h}\). Then we define the discrete trial- and test-spaces for the discretization of (4) as \(\mathbf{v}_{h}\), \(\boldsymbol{\psi}_{h}\in\mathbf{V}_{h}=[\mathbb{Q}_{h}^{(2)}]^{d}\) and \(p_{h}\), \(\xi_{h}\in L_{h}=\mathbb{Q}_{h}^{(2)}\). Since the resulting equal order finite element pair \(\mathbf{V}_{h}\times L_{h}\) does not fulfill the inf-sup condition, we add stabilization terms of local projection type [5]. We further add convection stabilization, also based on local projections [6]. The resulting semidiscrete variational problem then reads

\[\begin{split}&(\partial_{t}\mathbf{v}_{h},\,\boldsymbol{\psi}_{h})+( \mathbf{v}_{h}\cdot\triangledown\mathbf{v}_{h},\,\boldsymbol{\psi}_{h})+\nu( \triangledown\mathbf{v}_{h},\,\triangledown\boldsymbol{\psi}_{h})-(p_{h},\, \triangledown\boldsymbol{\psi}_{h})\\ &+\sum_{T\in\Omega_{h}}(\mathbf{v}_{h}\cdot\triangledown\mathbf{ v}_{h}-\pi_{h}\mathbf{v}_{h}\cdot\triangledown\pi_{h}\mathbf{v}_{h},\,\delta_{T}( \mathbf{v}_{h}\cdot\triangledown\boldsymbol{\psi}_{h}-\pi_{h}\mathbf{v}_{h} \cdot\triangledown\pi_{h}\boldsymbol{\psi}_{h}))=(\mathbf{f},\,\boldsymbol{ \psi}_{h})\quad\forall\boldsymbol{\psi}_{h}\in\mathbf{V}_{h}\\ &\qquad\qquad\qquad(\triangledown\mathbf{v}_{h},\,\xi_{h})+ \sum_{T\in\Omega_{h}}\alpha_{T}(\triangledown(p_{h}-\pi_{h}p_{h}),\triangledown (\xi_{h}-\pi_{h}\xi_{h}))=0\qquad\qquad\forall\xi_{h}\in L_{h}\end{split} \tag{5}\]

Here, we denote by \(\pi_{h}\,:\,\mathcal{Q}_{h}^{(2)}\to\mathcal{Q}_{h}^{(1)}\) the interpolation into the space of linear finite elements and by \(\alpha_{T},\beta_{T}\) two local stabilization parameters specified in (7).

### Time discretization

For temporal discretization, the time interval \([0,\,T]\) is split into discrete time steps of uniform size

\[0=t_{0}<t_{1}<\cdots<t_{N}=T,\;k=t_{n}-t_{n-1}.\]

The generalization to a non-equidistant time discretization is straightforward and only omitted for ease of presentation. We define \(\mathbf{v}_{n}\,:=\,\mathbf{v}_{h}(t_{n})\) and \(p_{n}\,:=\,p_{h}(t_{n})\) for the fully discrete approximation of velocity and pressure at time \(t_{n}\) and apply the second order Crank-Nicolson method to (5), resulting in the fully
discrete problem

\[(\mathbb{V}\cdot\mathbf{v}_{n},\,\xi_{h})+\sum_{T\in\Omega_{h}} \alpha_{T}(\mathbb{V}(p_{n}-\pi_{h}p_{n}),\mathbb{V}(\xi_{h}-\pi_{h}\xi_{h}))=0 \forall\xi_{h}\in L_{h},\] \[\frac{1}{k}(\mathbf{v}_{n},\,\boldsymbol{\phi}_{h})\,+\frac{1}{2 }(\mathbf{v}_{n}\cdot\mathbb{V}\mathbf{v}_{n},\,\boldsymbol{\phi}_{h})+\frac{ \nu}{2}(\mathbb{V}\mathbf{v}_{n},\,\mathbb{V}\boldsymbol{\phi}_{h})-(p_{n},\, \mathbb{V}\cdot\boldsymbol{\phi}_{h})\] \[+\sum_{T\in\Omega_{h}}(\mathbf{v}_{n}\cdot\mathbb{V}\mathbf{v}_{ n}-\pi_{h}\mathbf{v}_{n}\cdot\mathbb{V}\pi_{h}\mathbf{v}_{n},\,\delta_{T}( \mathbf{v}_{n}\cdot\mathbb{V}\boldsymbol{\phi}_{h}-\pi_{h}\mathbf{v}_{n}\cdot \mathbb{V}\pi_{h}\boldsymbol{\phi}_{h}))=\frac{1}{k}(\mathbf{v}_{n-1},\, \boldsymbol{\phi}_{h})\] \[+\frac{1}{2}(\mathbf{f}_{n},\boldsymbol{\phi}_{h})+\frac{1}{2}( \mathbf{f}_{n-1},\,\boldsymbol{\phi}_{h})-\frac{1}{2}(\mathbf{v}_{n-1}\cdot \mathbb{V}\mathbf{v}_{n-1},\,\boldsymbol{\phi}_{h})-\frac{\nu}{2}(\mathbb{V} \mathbf{v}_{n-1},\,\mathbb{V}\boldsymbol{\phi}_{h})\quad\forall\boldsymbol{ \phi}_{h}\in\mathbf{V}_{h}. \tag{6}\]

The right hand side only depends on the velocity \(\mathbf{v}_{n-1}\) at the last time step \(n-1\) and we will denote it as \(\mathbf{b}_{n-1}\) in the following.

The stabilization parameters \(\alpha_{T}\) and \(\delta_{T}\) depend on the mesh Peclet number and with two parameters \(\alpha_{0}>0\) and \(\delta_{0}\geq 0\) we define them as

\[\alpha_{T}=\alpha_{0}\left(\frac{\nu}{h_{T}^{2}}+\frac{|\mathbf{v}_{h}|_{ \infty}}{h_{T}}+\frac{1}{k}\right)^{-1},\quad\delta_{T}=\delta_{0}\left(\frac {\nu}{h_{T}^{2}}+\frac{|\mathbf{v}_{h}|_{\infty}}{h_{T}}+\frac{1}{k}\right)^{ -1}, \tag{7}\]

see [11] for details. Introducing the unknown \(\mathbf{x}_{n}=(p_{n},\,\mathbf{v}_{n})\), we write a short form of the equations (6) as

\[\mathcal{A}_{h}(\mathbf{x}_{n})=F_{h}, \tag{8}\]

where \([\mathcal{A}_{h}(\mathbf{x}_{n})]_{i}\) and \([F_{h}]_{i}\) are the left and right hand sides of (6) for all test functions \((\xi_{h}^{i},\,\boldsymbol{\phi}_{h}^{i})\).

In the presented form the Crank-Nicolson time discretization is sufficiently robust for smooth initial data \(v_{0}\); we refer to [26] for small modifications with improved robustness and stability.

### Solution of the algebraic systems

Discretization in space and time (6) leads to a nonlinear system of equations with a saddle-point character. The nonlinearity problem is solved by Newton iteration, where we usually omit the dependency of the stabilization parameters (7) on the velocity when setting up the Jacobian. The Jacobian matrix within the Newton iteration is kept for several iterations and also for consecutive time steps and only re-assembled when the nonlinear convergence rates deteriorate above a certain threshold, usually set to \(0.1\) or \(0.05\). As outer solver for the linear problems GMRES iteration [31] is employed which is preconditioned with a geometric multigrid solver [7].

To cope with the saddle-point structure of the incompressible Navier-Stokes equations and to allow for efficient shared memory parallelization, a Vanka smoother is used within the geometric multigrid. Therefore, small matrices describing the local problems on each mesh element are inverted in a parallel loop with Jacobi coupling. Using quadratic finite elements on hexahedral meshes these local matrices are of size \(108\times 108\). The complete layout of the algebraic solver and its parallelization is described in [21].

## 4 Neural Networks

Various neural network architectures are in principle suitable for our approach. Having previously worked with recurrent neural networks designed for sequence prediction, we settled for feedforward neural networks for this work. The simple structure makes them faster to train and to evaluate. Furthermore, we found that with proper regularization, the feedforward neural network described below performs more consistent than the alternatives we tried.
A general feedforward neural network is a composition of multiple parameterized nonlinear functions defined by

\[F(x):=F(x;\ \mathcal{W})=(f_{L}\circ\cdots\circ f_{1})(x;\ \mathcal{W})\,, \tag{9}\]

where \(\mathcal{W}\) denotes the set of optimizable (learnable) parameters. The \(i\)-th layer of the prototypical ANN we consider here consists of a weight matrix \(W_{i}\in\mathbb{R}^{n\times n_{i+1}}\) and a nonlinear function \(\sigma:\ \mathbb{R}\to\mathbb{R}\), called the activation function, which is applied component wise, \(\sigma\big{(}Wx\big{)}_{i}=\sigma\big{(}(Wx)_{i}\big{)}\). Furthermore, we employ two techniques common in deep learning to accelerate training and to improve generalization. Batch normalization \(N_{i}:\ \mathbb{R}^{n_{i+1}}\to\mathbb{R}^{n_{i+1}}\), applied before each activation, re-scales each component by

\[N_{i}(x)=\frac{x-\mu_{i}}{s_{i}}\odot\gamma_{i}+\beta_{i},\]

where \(x\odot y\) is the Hadamard product of two matrices (or tensors) and \(\mu_{i}\in\mathbb{R}^{n_{i+1}}\) and \(s_{i}\in\mathbb{R}^{n_{i+1}}\) are respectively, mean and standard deviation of samples estimated during training. The additional parameters \(\gamma_{i}\in\mathbb{R}^{n_{i+1}}\) and \(\beta_{i}\in\mathbb{R}^{n_{i+1}}\) are a learnable, component-wise affine transformation. Skip connections are inserted where input and output dimensions match, by adding the inputs of a layer to its outputs, leading to layers \(f_{i}\) of the structure

\[f_{1}(x) =\sigma_{1}(N_{1}(W_{1}x))\] \[f_{i}(x) =\sigma_{i}(N_{i}(W_{i}x))+x,\quad i=2,...\,,\,L-1\] \[f_{L}(x) =W_{L}x.\]

## 5 The Deep Neural Network Multigrid Method

In this section, we review the Deep Neural Network Multigrid (DNN-MG) solver, that we introduced in [46]. DNN-MG leverages a deep neural network to predict the correction of a coarse mesh solution that has been prolongated onto one or multiple finer mesh levels. The objective is to obtain solutions that are more accurate than using the coarse mesh alone, while increasing computational efficiency compared to performing full multigrid computations on the fine mesh levels. We develop the DNN-MG solver in a general formulation, while maintaining the connection to the Navier-Stokes equations, for which we developed DNN-MG. We begin by providing an overview of the modifications made in DNN-MG compared to the MG method during a time step of the Navier-Stokes simulation. We continue by describing the structure, inputs and outputs of the Neural-Network in DNN-MG. In particular, the network is designed to make localized predictions over specific parts of the simulation domain.

NotationIn the context of the Geometric Multigrid we always have a hierarchy of meshes \(\Omega_{h}^{0},...\,,\ \Omega_{h}^{L+J}\), obtained by successive refinement of the initial mesh \(\Omega_{h}^{0}\). We denote the finite element space \(\mathbf{V}_{h}\) and \(L_{h}\) introduced in Section 3.2 on level \(l\in\{0,...\,,\,L+J\}\) by \(V_{h}^{l}\) and \(L_{h}^{l}\). With their nested structure, these spaces reflect the grid hierarchy. Then, let \(\Omega_{h}^{l}\) be the quadrilateral or hexahedral finite element mesh on level \(l\in\{0,...\,,\,L+J\}\) with \(N_{\text{cells}}^{l}=\text{card}(\Omega_{h}^{l})\) cells and the index set \(Z_{l}\) of all global degrees of freedom with cardinality \(N_{\text{dof}}^{l}:=\text{card}(Z_{l})\). This notation and the description in Section 5.2 is inspired by the description of the Vanka smoother in [2].
### Time stepping using DNN-MG

We present a detailed description of one time step in the simulation of the Navier-Stokes equations using the DNN-MG solver. The computations involved are summarized in Algorithm 1.

At the beginning of time step \(n\), we solve for the unknown velocity \(\mathbf{x}_{n}^{L}\) and pressure \(p_{n}^{L}\) on the coarse level \(L\) using the classical Newton-Krylov simulation (Algorithm 1, lines 4-6) as described in Section 3.4. Subsequently, we prolongate \(\mathbf{x}_{n}^{L}\) to \(\tilde{\mathbf{x}}_{n}^{L+J}\) on the next finer level \(L+J\) (Algorithm 1, line 7) where a richer function space \(L_{h}^{L+J}\times\mathbf{V}_{h}^{L+J}\) is available.

Using the prolongated solution, we calculate the residual on level \(L+J\) for the input to the neural network (Algorithm 1, line 8). The residual is calculated according to (8). We then use the neural network of the DNN-MG solver to predict the velocity update \(\mathbf{d}_{n}^{L+J}\), which represents the difference \(\mathbf{d}_{n}^{L+J}=\tilde{\mathbf{x}}_{n}^{L+J}-\tilde{\mathbf{x}}_{n}^{L+J}\) between the prolongated \(\tilde{\mathbf{x}}_{n}^{L+J}\) and the unknown ground truth solution \(\tilde{\mathbf{x}}_{n}^{L+J}\) on level \(L+J\) (Algorithm 1, line 9). The prediction is based on \(\tilde{\mathbf{x}}_{n}^{L+J}\in\mathbf{V}_{h}^{L+J}\times L_{h}^{L+J}\) and \(\mathbf{r}_{n}^{L+J}\) and further utilizes information about the local mesh structure \(\omega\in\mathbb{R}^{n_{\text{coo}}}\). This information is extracted from \(\Omega_{h}^{L+J}\).

To conclude the \(n^{\text{th}}\) time step, we compute the right-hand side \(b_{n+1}^{L+J}\) of Eq. 6 on level \(L+J\) using the corrected velocity \(\tilde{\mathbf{v}}_{n}^{L+J}+\mathbf{d}_{n}^{L+J}\) (Algorithm 1, line 10), and then restrict it to level \(L\) (Algorithm 1, line 11). This right-hand side, incorporating the neural network-based correction, becomes part of the multigrid computations in the next time step (Algorithm 1, line 4), thereby influencing the overall time evolution of the flow. This approach of constructing the right-hand side \(b_{n+1}^{L+J}\) on level \(L+J\) and subsequently restricting it is a crucial aspect of the DNN-MG solver.

We note that in the case of the Navier-Stokes equations the pressure is handled implicitly on level \(L\) in the multigrid solve and is not directly receive a correction by the neural network. However, the pressure is included in the network's input through the prolongated solution and the residual, as it is an integral part of the solution and is indirectly influenced through the corrections to the velocity.

### The Neural Network of DNN-MG

The neural network component forms the core of the DNN-MG solver. It plays a crucial role in enhancing computational efficiency, facilitating fast training procedures, and ensuring robust generalization across
iverse flow regimes beyond the training set. We achieve this through the careful design of a compact neural network architecture with a moderate number of parameters and a localized, patch-based structure. Figure 2 provides an overview of the network component, illustrating the local approach, which we introduce next.

A patch-based neural networkTo ensure computational efficiency, the DNN-MG approach adopts a neural network that operates on localized patches of the mesh. The network's input and output are constrained to a patch, which reduces the computational cost compared to a prediction over the entire domain. This local approach leads to more compact neural networks, with a smaller number of parameters.

We define a patch \(P_{L,J}^{M}\subset\Omega_{h}^{L+J}\) as a collection of cells \(K\in\Omega_{h}^{L+J}\) on level \(L+J\). The neural network operates on a single patch \(P_{L,J}^{M}\) and predicts the velocity update \(\mathbf{d}_{n,P}^{L+J}\) specifically for the degrees of freedom within that patch. Since a single neural network is used in parallel on each patch, all patches must have the same structure and consist of the same number of fine mesh elements. Formally, let \(L(P_{L,J}^{M})\times\mathbf{V}(P_{L,J}^{M})\subset L_{h}^{L+J}\times\mathbf{V} _{h}^{L+J}\) be the function space associated with the degrees of freedom in patch \(P_{L,J}^{M}\). The neural network learns a mapping

\[\begin{split}\mathcal{N}\,:\,\left(L(P_{L,J}^{M})\times\mathbf{V} (P_{L,J}^{M})\right)^{2}\times\mathbb{R}^{n_{\text{Geo}}}&\to L(P_{ L,J}^{M})\times\mathbf{V}(P_{L,J}^{M})\,,\\ (\tilde{\mathbf{x}}_{n}^{L+J}|_{P_{L,J}^{M}},\,\mathbf{r}_{n}^{L+J }|_{P_{L,J}^{M}},\,\omega_{P_{L,J}^{M}})&\mapsto\mathbf{d}_{n,P}^ {L+J}\,,\end{split} \tag{10}\]

that predicts a local velocity update \(\mathbf{d}_{n,P}^{L+J}\in\mathbf{V}(P_{L,J}^{M})\) based on local data which is also restricted to the patch \(P_{L,J}^{M}\). This prediction is repeated independently for each patch in the domain, covering the entire domain \(\Omega_{h}^{L+J}\). The predictions of shared degrees of freedom in adjacent patches are averaged to ensure consistency. We now formulate the process we just sketched in more mathematical terms.

We define a patch \(P_{L,J}^{M}\), as a collection of cells on level \(L+J\) that are formed by \(J\) successive refinements of a coarse mesh cell \(M\in\Omega_{h}^{L}\). In other words, \(P_{L,J}^{M}\) comprises the \((2^{d})^{J}\) cells obtained by successive refinement of the coarse mesh cell \(M\in\Omega_{h}^{L}\) for \(J\) levels and can be defined rigorously by

\[P_{L,J}^{M}\coloneqq\{K\in\Omega_{h}^{L+J}\mid\tilde{K}\cap M\neq\emptyset,\, M\in\Omega_{h}^{L}\}. \tag{11}\]

Additionally, we define \(\Omega_{h}^{\text{patch}}\coloneqq\left\{P_{L,J}^{M}\text{ for }M\in\Omega_{h}^{L}\right\}\). Note that \(\text{card}(\Omega_{h}^{\text{patch}})=N_{\text{cells}}^{L}\). We denote the set of global degrees of freedom associated with the patch \(P_{L,J}^{M}\) by \(Z(P_{L,J}^{M})\subset Z_{L+J}\). The cardinality of \(Z(P_{L,J}^{M})\)

Figure 2: DNN-MG adopts a local approach where the neural network works on small neighborhoods of the simulation domain called patches, and independently provides a velocity correction for each one.
is denoted by \(N_{\text{dof}}^{P}:=\text{card}(Z(P_{L,J}^{M}))\). Further, the index set \(\hat{Z}(P_{L,J}^{M}):=\{0,...\,,\,N_{\text{dof}}^{P}-1\}\) contains all local degrees of freedom on \(P_{L,J}^{M}\). For a given patch \(P_{L,J}^{M}\) and a local degree of freedom with index \(\hat{\imath}\in\hat{Z}(P_{L,J}^{M})\) the mapping

\[\text{dof}\,:\,\Omega_{h}^{L+J}\times\hat{Z}(P_{L,J}^{M})\to Z_{l}\,,\quad(P_{L,J}^{M},\,\,\hat{\imath})\mapsto i\,, \tag{12}\]

provides the uniquely defined global index \(i\in Z_{L+J}\).

We can now define the prediction of the defect by the neural network in precise terms. It consists of multiple steps:

1. _Local restriction by \(\mathcal{Q}:\,\mathbb{R}^{N_{\text{dof}}^{L+J}}\to\mathbb{R}^{N_{\text{cells} }^{L}\times N_{\text{dof}}^{P}}\)_: Using a local restriction operator, the globally defined data of size \(N_{\text{dof}}^{L+J}\) is transformed into local, patchwise defined input of size \(N_{\text{dof}}^{P}\) suitable for the DNN. This is done for each patch, resulting in a matrix of size \(\mathbb{R}^{N_{\text{cells}}^{L}\times N_{\text{dof}}^{P}}\).
2. _Patch-wise prediction by the neural network \(\mathcal{N}:\,\mathbb{R}^{(2N_{\text{dof}}^{P}+n_{\text{Geo}})}\to\mathbb{R}^ {N_{\text{dof}}^{P}}\)_: For each patch, the neural network \(\mathcal{N}\) maps the local input data of size \((2N_{\text{dof}}^{P}+n_{\text{Geo}})\) to a local output of size \(N_{\text{dof}}^{P}\) predicting the velocity update for the degrees of freedom within the patch \(P_{L,J}^{M}\).
3. _Global extension by \(\mathcal{L}:\,\mathbb{R}^{N_{\text{cells}}^{L}\times N_{\text{dof}}^{P}}\to \mathbb{R}^{N_{\text{dof}}^{L+J}}\)_: The locally defined output of the DNN is processed by the extension operator \(\mathcal{L}\), which maps the local data of size \(N_{\text{cells}}^{L}\times N_{\text{dof}}^{P}\) to a single globally defined defect vector of size \(N_{\text{dof}}^{L+J}\).

In order to define \(\mathcal{Q}\), we need the \(P_{L,J}^{M}\)-local restriction \(\mathcal{Q}:\,\mathbb{R}^{N_{\text{dof}}^{L+J}}\times\Omega_{h}^{\text{Patch} }\to\mathbb{R}^{N_{\text{dof}}^{P}}\) defined by

\[(\mathcal{Q}(\mathbf{x},\,P_{L,J}^{M}))[\hat{\imath}]=\mathbf{x}[\text{dof}(P_ {L,J}^{M},\hat{\imath})]\,,\quad\text{for}\,\,\hat{\imath}\in\hat{Z}(P_{L,J}^{ M})\,. \tag{13}\]

Then, we define \(\mathcal{Q}\) such that it maps a global vector \(\mathbf{x}\in\mathbb{R}^{N_{\text{dof}}^{L+J}}\) to a matrix containing local restrictions to a patch for all \(P_{0},...\,,\,P_{N_{\text{cells}}^{L}-1}\in\Omega_{h}^{\text{Patch}}\), stacked on top of each other:

\[\mathcal{Q}:\,\mathbb{R}^{N_{\text{dof}}^{L+J}}\to\mathbb{R}^{N_{\text{cells} }^{L}\times N_{\text{dof}}^{P}},\quad\mathbf{x}\mapsto\left(\mathcal{Q}_{P}(x,P_{0}),...\,,\,\mathcal{Q}_{P}(x,\,P_{N_{\text{cells}}^{L}-1})\right)^{\top} \tag{14}\]

After applying \(\mathcal{Q}\) to the residual \(\mathbf{r}_{n}^{L+J}\) and \(\tilde{\mathbf{x}}_{n}^{L+J}\), we obtain the patch-wise, local data required as input to the DNN. According to (14), the local data is stacked along the first dimension, forming a single batch. The network is then evaluated independently for each patch, with each row of the input representing a different patch. The output of the DNN then consists of a batch of the patchwise defect \(\mathbf{d}_{n,P}^{L+J}\). To obtain the global defect vector \(\mathbf{d}_{n}^{L+J}\), we use the extension operator \(\mathcal{L}\) to transfer this batch of local data back to the global domain. To this end, we first introduce a \(P_{L,J}^{M}\)-local extension operator \(\mathbf{L}:\,\mathbb{R}^{N_{\text{dof}}^{P}}\times\Omega_{h}^{\text{Patch}}\to \mathbb{R}^{N_{\text{dof}}^{L+J}}\) defined by

\[(\mathbf{L}(\mathbf{y},\,P_{L,J}^{M}))[\imath]=\begin{cases}\mathbf{y}[\hat{ \imath}]\,,\quad\text{if}\,\,\exists\hat{\imath}\in\hat{Z}(P_{L,J}^{M}):\,\, \imath=\text{dof}(K,\hat{\imath}),\,K\in P_{L,J}^{M}\,,\\ 0\,,\qquad\text{if}\,\,\imath\notin Z_{l}(P_{L,J}^{M})\end{cases}\,.\]

We further need a scaling vector \(\mu\in\mathbb{R}^{N_{\text{dof}}^{L+J}}\) which contains the reciprocal of the valence of a degree of freedom, i. e. the number of patches a degree of freedom is contained in. Then we define \(\mathcal{L}\) such that it maps \(N_{\text{cells}}^{L}\) local vectors to one global vector as

\[\mathcal{L}:\,\mathbb{R}^{N_{\text{dof}}^{P}\times N_{\text{cells}}^{L}}\to \mathbb{R}^{N_{\text{dof}}^{L+J}}\,,\quad\mathbf{y}\mapsto\mu\odot\sum_{P_{L,J }^{M}\in\Omega_{h}^{\text{Patch}}}\mathbf{L}(\mathbf{y}[M],\,P_{L,J}^{M})\,. \tag{15}\]
For the sake of brevity, we used \(M\) not only as an element of \(\Omega_{h}^{L}\) in (15) but also to index the rows of the collection of local vectors.

This concludes the description of all operations necessary to integrate the neural network with its local approach into the multigrid solver. In particular

\[\mathcal{L}\circ\mathcal{N}(Q(\tilde{\mathbf{x}}_{n}^{L+J}),\,Q(\mathbf{r}_{n}^ {L+J}),\,\Omega_{h}^{L+J}):\,L_{h}^{L+J}\times\mathbf{V}_{h}^{L+J}\to L_{h}^{L +J}\times\mathbf{V}_{h}^{L+J}\,,\quad\tilde{\mathbf{x}}_{n}^{L+J}\mapsto \mathbf{d}_{n}^{L+J}\,, \tag{16}\]

is a well-defined operation (cf. Algorithm 1 line 9). The mapping (10) also fits into this framework and interpreting the definition (13) in a finite element space context, the restrictions \(\tilde{\mathbf{x}}_{n}^{L+J}|_{p_{L,J}^{M}}\) and \(\mathbf{r}_{n}^{L+J}|_{p_{L,J}^{M}}\) are obtained. The network operates on each patch independently, which allow us to take advantage of parallelization through the localized nature of the problem. This shows that the local approach is a key feature in our method, which enables the network to handle a wide range of flow scenarios while maintaining computational efficiency, as we will show later in Section 6. The rest of this section consists of a detailed description of the inputs to the neural network and the training methodology.

Neural network inputsThe effectiveness of neural network prediction hinges upon the selection of appropriate network inputs that provide enough information regarding the coarse flow characteristics at level \(L\) and the local patch geometry, akin to the significance of mesh cell geometry in classical finite element simulations. A well-informed choice of these inputs holds the potential to keep the number of parameters in the neural network low, thereby reducing evaluation time during runtime and minimizing the data and time prerequisites for training.

In our implementation, we use the following inputs to the neural network:

* nonlinear residuals \(\mathbf{r}_{n}^{L+J}|_{p_{L,J}^{M}}=\left(F_{h}^{L+J}-\mathcal{A}_{h}^{L+J}( \tilde{\mathbf{x}}_{n}^{L+J})\right)|_{p_{L,J}^{M}}\in\mathbb{R}^{N_{ud}^{p}}\) of the prolongated coarse mesh solution for (6);
* the velocities and pressure \(\tilde{\mathbf{x}}_{n}^{L+J}|_{p_{L,J}^{M}}\in\mathbb{R}^{N_{ud}^{p}}\) on mesh level \(L+J\);
* the geometry of the cells, which for the 3D case consist of
* the edge lengths \(h^{c}\in\mathbb{R}^{12}\);
* the lengths of the diagonals connecting the vertices which do not share the same face \(a^{c}\in\mathbb{R}^{4}\);
* for each vertex, the average of the 3 angles between the faces \(a^{c}\in\mathbb{R}^{8}\).

The nonlinear residual plays a crucial role as it quantifies the local error or defect at each degree of freedom on level \(L+J\) and therefore carries essential information about the underlying partial differential equation. Furthermore, the velocity and pressure fields contain details about the current solution and provide knowledge about the flow within the current patch. It is important to note that the geometric information of the cells also plays a significant role in our approach, as we impose no specific restrictions on the cells, except for the standard shape regularity requirements commonly employed in finite element methods.

Training of DNN-MGThe training of DNN-MG is based on simulations of the Navier-Stokes equations for which a multigrid representation of the velocity \(\mathbf{v}_{n}^{I}\) with two levels \(L\) and \(L+J\) is available. The velocity \(\mathbf{v}_{n}^{L+J}\) thereby serves as ground truth \(\tilde{\mathbf{v}}_{n}^{L+J}\). The goal of the training is to optimize the network parameters such that the norm \(\left\|(\tilde{\mathbf{v}}_{n}^{L+J}+\mathbf{d}_{n}^{L+J})-\tilde{\mathbf{v}}_ {n}^{L+1}\right\|\) of the difference between the predicted velocity \(\tilde{\mathbf{v}}_{n}^{L+J}+\mathbf{d}_{n}^{L+J}\) (i. e. the velocity after the correction) and the ground truth \(\tilde{\mathbf{v}}_{n}^{L}\) is minimized throughout the simulation. Typically, the training data for the network comprises only a few snapshots of simulation data from simple scenarios, as we will see in the next section.
## 6 Numerical examples

In the following paragraphs we will document different numerical simulations. The test cases are based on a well established 3D Navier-Stokes benchmark problem [55, 10]. While the neural network is trained on the original benchmark configuration, we study its performance and generalization on a set of modified problems. As perturbations we consider increasing or decreasing the Reynolds number and also a substantial change to the problem geometry by introducing a second obstacle. It is important to stress that the network is trained only once and no retraining is applied for the further test cases.

### Implementation aspects

The DNN-MG method is implemented mainly using two libraries: The FEM library Gascoigne3D [8] and the machine learning library PyTorch [52]. The Newton-Krylov geometric multigrid implemented in Gascoigne3D method has been shown to be efficient in [21].

We implement the deep neural networks with Libtorch, the C++ interface of PyTorch. Using C++ for the implementation allows us to couple the numerical methods and neural networks in a performant manner without unnecessary overhead. PyTorch also supports parallelization with MPI, which is used here to distribute training on multiple GPUs.

### The 3D Benchmark Flow Around a Cylinder

We examine a variant of the three-dimensional flow benchmark presented in [55] in a setting similar to the _3D-2Z_ one, with different Reynolds numbers \(Re\in\{200,240\}\). We show the geometry in Figure 3. While the original benchmark description considers an obstacle with circular cross-section, we allow elliptical ones for testing the generalizability of the approach. We define the Reynolds number as

\[\mathrm{Re}=\frac{\bar{\mathbf{v}}\cdot L}{v}, \tag{17}\]

where \(\bar{\mathbf{v}}\) is the average inflow velocity and \(L\) the length of the major axis of the obstacle. When the cross-section of the cylinder is circular, \(L\) represents its diameter. If the cross-section is elliptical, \(L\) corresponds to the height of the obstacle, as the major axis is always parallel to the \(y\)-axis.

The time step size is chosen as \(k=0.008\) on the interval \(I=(0,8]\). The initial velocity is \(\mathbf{v}(0)=\mathbf{0}\), and the inflow boundary condition is prescribed for \(t>0\), with a smooth startup process. The flow is driven

Figure 3: Geometry of the training and test scenario with a parabolic inflow profile at \(\Gamma_{\mathrm{in}}\), do-nothing boundary conditions at the outflow boundary \(\Gamma_{\mathrm{out}}\) and no-slip conditions on the walls \(\Gamma_{\mathrm{wall}}\). The center of the obstacle is at (0.5, 0.2).
by a Dirichlet profile \(\mathbf{v}=\mathbf{v}^{D}\) at the left boundary \(\Gamma_{\text{in}}\) given by

\[\begin{split}\mathbf{v}_{in}^{3d}(t,\,y,\,z)&=\omega( t)\frac{y(H-y)(H^{2}-z^{2})}{(H/2)^{2}H^{2}}\frac{9}{8}\bar{\mathbf{v}}_{3d}\,\,\,\text{ on}\,\,\,\Gamma_{in}\,:=0\times[0,\,H],\\ \omega(t)&=\begin{cases}\frac{1}{2}-\frac{1}{2}\cos(5 \pi t),&t\leq\frac{1}{5}\\ 1&t>\frac{1}{5}.\end{cases}\end{split} \tag{18}\]

\(H=0.41\) is the height and width of the channel and \(\omega(t)\) acts, as a regularization of the startup phase. On the wall boundary \(\Gamma_{wall}\) and on the obstacle we use no-slip boundary conditions. On the outflow boundary \(\Gamma_{\text{out}}\) we use a do-nothing outflow condition [27]. The average flow rate is \(\bar{\mathbf{v}}_{3d}=1\) and the viscosity is \(v=5\times 10^{-4}\). The training data is generated considering the circular obstacle, i. e. \(L=0.1\) which results in the Reynolds number \(\text{Re}=200\). For testing, the elliptical cross-section with \(L=0.12\) is considered, leading to the slightly increased Reynolds number \(\text{Re}=240\). See Section 5.2 for details.

As in the classical benchmark description we analyze the drag and lift coefficients of the obstacle

\[J_{\text{drag}}=\int_{\Gamma}\big{(}v\forall\mathbf{v}\mathbf{n}-p\mathbf{n} \big{)}\mathbf{e}_{x}\,\mathrm{d}s,\quad J_{\text{lift}}=\int_{\Gamma}\big{(} v\forall\mathbf{v}\mathbf{n}-p\mathbf{n}\big{)}\mathbf{e}_{y}\,\mathrm{d}s. \tag{19}\]

To compute the drag and lift coefficients, we use the Babuszka-Miller trick; cf. [10, 3] and rewrite the surface integrals over the volume for obtaining super convergence of fourth order.

### Neural Network parameters

DNN-MG uses a neural network that operates on patches. In our experiments, a patch is once or twice refined cell on level \(L\). In Table 2 we list the dimensions and number of trainable parameters of the resulting neural networks. The hidden size is arbitrary, but a size of 512 or 750 has been shown to give good results, while keeping the network size relatively small. We generally use 8 hidden layers of the same size. Unless mentioned otherwise, the results presented here were obtained with a hidden size of 750. Networks with a width of 512 are considered only for comparison.

### Neural Network training

Training data are generated on grids with three successive levels of refinement, see Table 1. Here, we denote by MG(3) the coarse grid. The grid levels MG(4) and MG(5) are the training data for a network prediction over one, respectively over two grid levels.

A single simulation is sufficient since, by the patch-wise application of the network, a single simulation provides \(N_{c}\times N_{T}\) training items, where \(N_{c}\) is the number of patches and \(N_{T}\) is the length of the time series. From the experience gained in the 2D case [46], we consider a small subinterval of \(I=[0,\,8]\). Here we

\begin{table}
\begin{tabular}{l l r r} \hline \hline  & Level & \# DoF & Patch size \\ \hline Coarse & MG(3) & 390 720 & 0 \\ Reference for DNN-MG(\(3+1\)) & MG(4) & 3 003 520 & 8 \\ Reference for DNN-MG(\(3+2\)) & MG(5) & 23 546 112 & 64 \\ \hline \hline \end{tabular}
\end{table}
Table 1: Spatial meshes with increasing number of multigrid levels and degrees of freedom (DoF). For the DNN-MG setups we indicate the size of the patches (counted as number of fine mesh elements) where the network is applied.
use \(t\in[4,\,7]\) resulting in \(N_{T}=375\). The interval from [2, 4] is used as validation set. The full dataset consists of 1 TB of data.

The networks were trained using 2 GPU nodes, each with 2 Nvidia A100 GPUs and 2 Intel Xeon Platinum 8360Y CPUs. The number of MPI processes were chosen equal to the number of GPUs such that one MPI process uses one GPU and CPU. Using this setup, the training of the feedforward network takes one day, regardless of the network size. Although the data loading process is implemented efficiently and the average load of the GPUs ranges from 80% for small networks to 90% for large networks, the limiting factor in terms of performance is the memory bandwidth.

As loss function \(\mathcal{L}\) we employed a simple \(l^{2}\)-loss

\[\mathcal{L}=\sum_{n=1}^{N_{T}}\frac{1}{N_{c}}\sum_{c=1}^{N_{c}}l^{2}(\tilde{u }_{n,c}^{1+J}+d_{n,c}^{L+J},\,\tilde{u}_{n,c}^{L+J})^{2}=\sum_{i=1}^{N_{T}} \frac{1}{N_{c}}\sum_{c=1}^{N_{c}}\left|\left(\tilde{v}_{n,c}^{L+J}+d_{n,c}^{L+ J}\right)-\tilde{v}_{n,c}^{L+J}\right|_{2}^{2},\quad J\in\{1,\,2\}.\]

The parameter \(J\) is the number of levels we skip with the prediction by the DNN. We use a Tikhonov-regularization with a scaling factor \(\alpha=10^{-5}\). To optimize our model, we utilize the AdamW [37] optimizer with a maximum of 1000 epochs. Finally, we select the model with the lowest validation loss. The convergence of the training for the different network configurations is reported in Figure 4.

\begin{table}
\begin{tabular}{r r r r r} \hline \hline Predicted levels & \(N_{\text{in}}\) & \(N_{\text{hidden}}\) & \(N_{\text{out}}\) & Trainable parameters \\ \hline
1 & 240 & 512 & 81 & 2 007 552 \\
2 & 1024 & 512 & 375 & 2 559 488 \\
1 & 240 & 750 & 81 & 4 190 250 \\
2 & 1024 & 750 & 375 & 4 998 750 \\ \hline \hline \end{tabular}
\end{table}
Table 2: The input size \(N_{\text{in}}\), hidden size \(N_{\text{hidden}}\), output size \(N_{\text{out}}\) and the number of trainable parameters of the feedforward network, for different numbers of predicted levels.

Figure 4: The average training and validation loss of the feedforward network architecture with hidden size (a) 512 and (b) 750. Due to the data efficiency of DNN-MG and the large number of batches (\(\approx\) 900 000) relative to the network size, we were able to achieve a good training loss after a single epoch and there is only little improvement in the training loss beyond that point.
### Flow prediction

For testing we use the flow around an elliptic obstacle with an increased height of \(0.12\) compared to the training data (see Figure 3). The resulting finite element meshes have the same structure as in the training case, but the elements are distorted. We observe that DNN-MG is indeed able to predict high frequency fluctuations that are not visible in the coarse mesh solution. In particular in the vicinity of the obstacle the quality of the solution is strongly enhanced with distinct features in the wake being apparent in the DNN-MG simulation.

In addition to a viscosity of \(v=5\times 10^{-4}\) and \(Re=240\), we consider fluids with \(v=4\times 10^{-4}\) corresponding to a Reynolds number of \(Re=300\) and \(v=6.67\times 10^{-4}\) which corresponds to a Reynolds number of \(180\). Both of these Reynolds numbers are calculated for the elliptic obstacle with a height of \(0.12\), which is our typical test case. The network is trained only for \(Re=200\) using the obstacle with circular cross-section.

Channel flow at \(\mathbf{Re=180}\)In Figure 6 and Table 3, the drag and lift forces at Reynolds number \(180\), computed using DNN-MG and the classical MG method on different levels are presented. The results show the predictive capacity of the DNN-MG approach in reproducing the dynamic behavior of the flow. In terms of the drag coefficient, DNN-MG yields a slight improvement in terms of the magnitude. The lift coefficients on the other hand improve greatly compared to the coarse solution, as the overall flow dynamics are well reconstructed, albeit with temporal deviations. We observe this in Figure 5 too, where different flow patters emerge in the wake of the flow. Furthermore, with respect to the minimum and maximum lift values, our approach is in good agreement with the reference solution. Experience with finite element simulations usually show that the lift is the functional that is harder to approximate by far [10].

Channel flow at \(\mathbf{Re=240}\)Similarly, Figure 7 and Table 4 depict the drag and lift forces at Reynolds number \(240\), obtained through DNN-MG and the multigrid solution on the coarse and reference level. Figure 5 show a comparison of the velocity fields obtained using DNN-MG with the coarse mesh and reference solution. The local corrections improve the solution at a global level and DNN-MG is clearly able to capture and improve the flow characteristics. In particular, it has the capacity to significantly enhance the flow resolution, yielding improved solutions and reproducing intricate features that were previously unobservable at lower levels. This underscores the method's efficacy in capturing fine-scale details. We see similar behavior as in \(\mathbf{Re=240}\), indicating good generalization of DNN-MG to different Reynolds numbers, considering the network was trained on \(\mathbf{Re=200}\).
Figure 5: Comparison of the magnitude of the velocity field for the channel with a cylindrical obstacle with an elliptical cross-section at \(\mathrm{Re}=240\) and time \(t=7.5\). The results are shown for DNN-MG(\(3+2\)), along with the reference solution obtained with 5 levels (MG(5)), as well as the coarse mesh solution at level 3 (MG(3)). The DNN-MG(\(3+2\)) method utilizes a network trained on simulation data from a round obstacle. The first three images show a cross-section of the domain in the \(xy\)-plane and the others show a cross-section of the domain in the \(xz\)-plane.
\begin{table}
\begin{tabular}{l c c c c c c c c} \hline \hline Type & \(\min C_{d}\) & \(\max C_{d}\) & \(\overline{C_{d}}\) & amp \(C_{d}\) & \(\min C_{l}\) & \(\max C_{l}\) & \(\overline{C_{l}}\) & amp \(C_{l}\) \\ \hline MG(3) & 1.7571 & 1.7662 & 1.7591 & 0.0091 & 0.0083 & 0.0235 & 0.0157 & 0.0152 \\ MG(4) & 1.7606 & 1.7713 & 1.7638 & 0.0107 & -0.0016 & 0.0208 & 0.0104 & 0.0224 \\ MG(5) & 1.8472 & 1.8595 & 1.8494 & 0.0123 & -0.0178 & 0.0069 & -0.0033 & 0.0247 \\ DNN-MG(3 + 1)-512 & 1.7082 & 1.7575 & 1.7186 & 0.0493 & 0.0025 & 0.0241 & 0.0136 & 0.0216 \\ DNN-MG(3 + 2)-512 & 1.8593 & 1.9112 & 1.8742 & 0.0519 & -0.0223 & 0.0025 & -0.0119 & 0.0248 \\ DNN-MG(3 + 1)-750 & 1.7257 & 1.7652 & 1.7374 & 0.0394 & -0.0235 & 0.0045 & -0.0049 & 0.0280 \\ DNN-MG(3 + 2)-750 & 1.7657 & 1.8052 & 1.7774 & 0.0394 & -0.0235 & 0.0045 & -0.0049 & 0.0280 \\ \hline \hline \end{tabular}
\end{table}
Table 3: The min, max, mean and amplitude of drag and lift functionals in the same setting as above in Figure 6. The best results in comparison to the \(MG\)(5) solution are highlighted with color coding: 1st best, 2nd best, 3rd best.

Figure 6: The drag (left) and lift (right) functionals for the channel with one obstacle with an elliptical cross section at \(\mathrm{Re}=180\) are shown. The results are presented for three discretizations: The coarse mesh MG(3), the reference MG(5), and the DNN-MG(\(3+2\)) improved by a neural network.
\begin{table}
\begin{tabular}{l c c c c c c c c} \hline \hline Type & \(\min C_{d}\) & \(\max C_{d}\) & \(\overline{C_{d}}\) & amp \(C_{d}\) & \(\min C_{l}\) & \(\max C_{l}\) & \(\overline{C_{l}}\) & amp \(C_{l}\) \\ \hline MG(3) & 1.6810 & 1.6920 & 1.6840 & 0.0111 & 0.0081 & **0.0495** & 0.0255 & 0.0415 \\ MG(4) & 1.6783 & 1.6938 & 1.6869 & 0.0155 & -0.0204 & 0.0373 & 0.0098 & 0.0576 \\ MG(5) & 1.7724 & 1.8096 & 1.7910 & 0.0372 & -0.0303 & **0.0852** & **0.0123** & **0.1155** \\ DNN-MG(3 + 1)-512 & 1.6305 & 1.6748 & 1.6425 & 0.0444 & -0.0364 & 0.0224 & -0.0002 & 0.0588 \\ DNN-MG(3 + 2)-512 & 1.6866 & 1.7931 & 1.7537 & 0.1066 & -0.0347 & 0.0363 & 0.0089 & 0.0707 \\ DNN-MG(3 + 1)-750 & 1.6215 & 1.6664 & 1.6309 & 0.0448 & 0.0008 & 0.0304 & **0.0146** & 0.0296 \\ DNN-MG(3 + 2)-750 & **1.6966** & 1.7167 & 1.7079 & 0.0200 & -0.0345 & 0.0426 & 0.0059 & **0.0774** \\ \hline \hline \end{tabular}
\end{table}
Table 4: The min, max, mean and amplitude of drag and lift functionals in the same setting as above in Figure 7. The best results in comparison to the \(MG\)(5) solution are highlighted with color coding: 1st best, 2nd best, 3rd best.

Figure 7: The drag (left) and lift (right) functionals for the channel with one obstacle with an elliptical cross section at \(\mathrm{Re}=240\) are shown. The results are presented for three discretizations: The coarse mesh MG(3), the reference MG(5), and the DNN-MG(\(3+2\)) improved by a neural network.
\begin{table}
\begin{tabular}{l c c c c c c c c} \hline \hline Type & \(\min C_{d}\) & \(\max C_{d}\) & \(\overline{C_{d}}\) & amp \(C_{d}\) & \(\min C_{l}\) & \(\max C_{l}\) & \(\overline{C_{l}}\) & amp \(C_{l}\) \\ \hline MG(3) & 1.6291 & 1.6449 & 1.6339 & 0.0159 & 0.0081 & 0.0337 & 0.0185 & 0.0256 \\ MG(4) & 1.6125 & 1.6507 & 1.6336 & 0.0381 & -0.0592 & 0.0724 & 0.0142 & 0.1316 \\ MG(5) & 1.7375 & 1.7794 & 1.7578 & 0.0419 & -0.0527 & 0.0404 & -0.0058 & 0.0931 \\ DNN-MG(3 + 1)-512 & 1.4532 & 1.5883 & 1.4985 & 0.1351 & -0.0047 & **0.0355** & 0.0161 & 0.0401 \\ DNN-MG(3 + 2)-512 & 1.4923 & 1.7125 & 1.5776 & 0.2202 & -0.0071 & 0.0656 & 0.0266 & 0.0728 \\ DNN-MG(3 + 1)-750 & 1.5660 & 1.6347 & 1.5858 & 0.0687 & -0.0235 & 0.0714 & 0.0184 & **0.0950** \\ DNN-MG(3 + 2)-750 & **1.6490** & **1.6910** & **1.6626** & **0.0419** & -0.0655 & 0.0664 & **0.0109** & 0.1319 \\ \hline \hline \end{tabular}
\end{table}
Table 5: The min, max, mean and amplitude of drag and lift functionals in the same setting as above in Figure 8. The best results in comparison to the \(MG\)(5) solution are highlighted with color coding: 1st best, 2nd best, 3rd best.

Figure 8: The drag (left) and lift (right) functionals for the channel with one obstacle with an elliptical cross section at \(\mathrm{Re}=300\) are shown. The results are presented for three discretizations: The coarse mesh MG(3), the reference MG(5), and the DNN-MG(\(3+2\)) improved by a neural network.
Figure 9: Relative velocity errors (first row) and pressure errors (second row) for coarse solutions with and without ANN correction (MG(3) and DNN-MG(\(3+2)-750\)) compared to the reference solution on level \(L+2\) for different Reynolds numbers. As further evidence of the innovative capacity of DNN-MG we show the errors of the MG(4) solution as well, which we consistently outperform.
### Performance measurements

For the DNN-MG method to be advantageous, it must achieve two goals: Firstly, it must enhance the accuracy compared to the solution on a coarse mesh. Secondly, it must reduce the computation time compared to a solution on a fine mesh. However, the Newton-Krylov geometric multigrid method already is of optimal complexity [33, 21]. In [46] this already raised the question what advantage DNN-MG offers, and we addressed it there for the 2D case. Here we adapt the theoretical findings to the 3D case.

Complexity analysisThe Newton-Krylov geometric multigrid method has linear complexity \(\mathcal{O}(N)\) in the number of DoFs \(N\). For each global mesh refinement the number of DoFs increase by a factor of 8, i. e. \(N^{L+1}\approx 8N^{L}\). The constant hidden in \(\mathcal{O}(N)\), however, can be very significant, since on average 5 Newton steps are required in each time step and within each Newton step one has to compute on average 12 GMRES steps with one sweep of the geometric multigrid solver as preconditioning. In addition, the geometric multigrid method usually involves six pre-smoothing and six post-smoothing steps. Thus, on each mesh level \(L\), \(L-1\),..., 0, a total of approximately 700 Vanka smoothing steps must be performed. One Vanka step thereby requires the inversion of about \(\mathcal{O}(N^{L}/8)\) small matrices of size \(108\times 108\), one on each element of the mesh, resulting in approximately \(108^{3}\simeq 1\,259\,712\) operations. Since the complexity of all mesh levels sums up to about \(N^{L}+N^{L-1}+...\,N_{0}\approx N^{L}(1+8^{-1}+...+8^{-L})\approx\frac{8}{7}N ^{L}\) we can estimate the effort for the complete solution process on level \(L\) as \(5\cdot 12\cdot(6+6)\cdot 108^{3}\cdot\frac{8}{7}\approx 10^{9}N^{L}\). We thereby only counted the effort for smoothing and neglected all further matrix, vector and scalar products. If we were to solve the problem on mesh level \(L+1\), the required effort would increase by a factor of 8 to approximately \(\approx 8\cdot 10^{9}N^{L}\) operations.

On the other hand, the DNN-MG method only requires the prolongation of the solution to the next finer mesh and the evaluation of the neural network on each patch. If we again consider patches of minimal size (one patch is one element of mesh layer \(L+1\)) about \(\mathcal{O}(N^{L+1}/8)=\mathcal{O}(N^{L})\) patches must be processed. The effort for the evaluation of the network can be estimated by the number of trainable parameters with \(N_{\text{c}}\) inputs. The DNN-MG approach requires only one evaluation of the network in each

Figure 10: On the left the runtime of DNN-MG and MG is plotted over multiple levels of refinement. We see that DNN-MG has less than 50% overhead compared to a coarse mesh simulation. On the right we see the contributions of different parts of the program to the wall clock time in percent. We observe that the evaluation of the ANN takes less than 3% of the time.
time step. The number of trainable parameters is for all network models of the order of single digit millions. Hence, the effort for correcting the level \(L\) Newton-Krylov solution by the neural network on level \(L+1\) can be estimated by \(10^{7}N^{L}\), which is almost a factor of thousand lower than number of FLOPS needed to solve on \(L+1\) (\(\approx 8\cdot 10^{9}N^{L}\)).

For the prediction of two levels, the estimate for the computational cost of DNN-MG remains largely unchanged. While the size of the patches and the inputs and outputs of the networks increase, the number of patches and network evaluations don't change. Surprisingly, we don't need to significantly increase the network size and therefore give an upper bound of \(2\cdot 10^{7}N^{L}\) trainable parameters. This is even more favorable than the prediction of a single mesh level, as the numerical solution on level \(L+2\) incurs a much higher cost of approximately \(64\cdot 10^{9}N^{L}\).

Computational performance of DNN-MGIn Figure 10 (a) the runtimes of DNN-MG and a purely numerical MG solution are plotted. While for the MG solution with each added mesh level, the wall time increases by a factor of 8, the wall time of DNN-MG increases at most by a factor of 1.5. Note that the wall time of DNN-MG(\(3+2\)) and MG(\(5\)) thereby differ by a factor of 35. Although this has to be weighed against the loss in accuracy, DNN-MG is clearly capable of increasing the computational efficiency. Moreover, DNN-MG(\(3+2\)) consistently outperforms MG(\(4\)) in terms of the accuracy in previous tests. Remarkably, it still has an advantage in terms of the runtime which once more underlines the computational efficiency of DNN-MG.

As can be seen in Figure 10 (b), the time spent in the solver in DNN-MG is constant, which shrinks the share of the solver in the wall time while other parts of the program gain weight. The main increase in computation time is due to the increased cost of the assembly of the right hand side and the evaluation of the functionals, since these calculations are done on the higher levels. Note that these operations are also of order \(O(N)\), but with a much smaller constant than other operations in the program, e. g. the operations mentioned above in our theoretical considerations. This makes the DNN-MG method highly efficient for \(1\) or \(2\) level predictions. More levels are not feasible at this point due to limitations in the underlying FEM library. Further investigations in the scalability of the method when predicting more mesh-levels is part of future work. Predicting weights of higher order polynomials instead of the weights of basis functions on refined meshes could further reduce the computational burden and increase the computational efficiency.

The evaluation of the DNN takes less than 3% of the wall time, which is less than our previous results [46] in 2D. This is due to substantial improvements of our implementation: collecting the input for the neural network and processing the output is parallelized and the evaluation of the network is done on a GPU. Therefore, DNN-MG not only reduces the cost of accurate simulations, but it is also well-suited for heterogeneous computing environments.

### Generalizability

In order to ensure the practicality of the DNN-MG method, it must generalize well to similar flows beyond those seen during training. Previous results on the 3D benchmark demonstrate the network's ability to do so under small geometric perturbations and varying Reynolds number. Here, we evaluate the network's performance under more substantial changes considering a channel with two obstacles shown in Figure 11. As for the single obstacle case, we test the performance under varying Reynolds numbers with a viscosity of \(v=5\times 10^{-4}\) and \(Re=240\), \(v=4\times 10^{-4}\) (\(Re=300\)) and \(v=6.67\times 10^{-4}\) (\(Re=300\)). We adopt some settings from the single obstacle case. For instance, the boundary conditions, in particular the inflow at the left boundary \(\Gamma_{\text{in}}\) given by (18), the initial velocity \(\mathbf{v}(0)=\mathbf{0}\) and we chose the time step size \(k=0.008\) on the interval \(I=(0,8]\).
We reuse the network trained for the single-obstacle channel at Reynolds number \(Re=200\). The neural network component of DNN-MG operates on level 5 as fine resolution and level 3 as coarse resolution in this section. In order to not overburden this paper we leave out the detailed results of a single level prediction. The complexity of the finite element meshes is given in Table 6.

Channel with two obstacles at \(\mathbf{Re=180}\)In Figure 13 and Table 7, the drag and lift coefficients obtained using DNN-MG at Reynolds number 180 are presented. The results demonstrate the method's ability to accurately predict the dynamics of the flow. In terms of the drag, the improvement by DNN-MG greatly improves the magnitude of the drag and closely follows the trajectory of the reference solution up to the small-scale dynamics. The lift is very similar in that we are able to reconstruct the overall dynamics of the flow, although we do not get the same temporal behavior. We also observe this in Figure 12, through the different flow patterns. In terms of min and max of the lift we are very close to the reference solution.

Channel with two obstacles at \(\mathbf{Re=240}\)Likewise, Figure 14 and Table 8 show the drag and lift forces at Reynolds number 240, further substantiating the accuracy of DNN-MG and its ability to improve the solution. Figure 12 further demonstrates the effectiveness of DNN-MG in accurately capturing the velocity magnitudes for scenarios on which the DNN was not trained. The figure shows the capability of the method to effectively increase the resolution of the flow, to improve the solution and to add features that could not be observed on lower levels, which are present on higher levels. This is a great example of the efficacy of the localized predictions of DNN-MG, which enable this generalization. Analogous to the case with one obstacle, we observe good generalization of DNN-MG to different Reynolds numbers. The great results we found for the drag and lift coefficients and the velocity field at \(\mathbf{Re=180}\), translate similarly to \(\mathbf{Re=240}\) case.

\begin{table}
\begin{tabular}{l r r r} \hline \hline  & Level & \# DoF & Patch size \\ \hline Coarse & MG(3) & 533 148 & 0 \\ Reference for DNN-MG(3 + 1) & MG(4) & 4 097 340 & 8 \\ Reference for DNN-MG(3 + 2) & MG(5) & 32 115 324 & 64 \\ \hline \hline \end{tabular}
\end{table}
Table 6: Spatial meshes with increasing number of multigrid levels and unknowns.

Figure 11: Geometry of the training and test scenario with a parabolic inflow profile at \(\Gamma_{\text{in}}\), do-nothing boundary conditions at the outflow boundary \(\Gamma_{\text{out}}\) and no-slip conditions on the walls \(\Gamma_{\text{wall}}\). The center of the obstacle is at (0.5, 0.2).
Figure 12: Comparison of the magnitude of the velocity field for the channel with two cylindrical obstacles at time \(t=7.5\). The first obstacle has a circular cross-section and the second obstacle has an elliptical cross-section, which results in a Reynolds number of \(\mathrm{Re}=240\). The results are shown for DNN-MG(\(3+2\)), along with the reference solution obtained with 5 levels (MG(5)), as well as the coarse mesh solution at level 3 (MG(3)). The DNN-MG(\(3+2\)) method utilizes a network trained on simulation data from a round obstacle. The first three images show a cross-section of the domain in the \(xy\)-plane and the others show a cross-section of the domain in the \(xz\)-plane.
\begin{table}
\begin{tabular}{l c c c c c c c c} \hline \hline Type & \(\min C_{d}\) & \(\max C_{d}\) & \(\widehat{C_{d}}\) & amp \(C_{d}\) & \(\min C_{l}\) & \(\max C_{l}\) & \(\widehat{C_{l}}\) & amp \(C_{l}\) \\ \hline MG(3) & 0.7674 & 1.3213 & 1.2176 & 0.5539 & -0.0087 & 0.0313 & **0.0097** & 0.0400 \\ MG(4) & 0.7905 & 1.5601 & 0.7696 & 1.3747 & -0.0091 & 0.0346 & 0.0061 & 0.0404 \\ MG(5) & 1.5089 & 1.7903 & 1.6636 & 0.2814 & -0.0872 & **0.1428** & **0.0097** & 0.2299 \\ DNN-MG(\(3+1\))-512 & 1.5179 & 1.7259 & 1.6342 & 0.2079 & -0.0487 & **0.0907** & 0.0280 & 0.1394 \\ DNN-MG(\(3+2\))-512 & 1.6038 & 1.6457 & 1.6265 & 0.0419 & -0.0432 & -0.0112 & -0.0237 & 0.0320 \\ DNN-MG(\(3+1\))-750 & 1.5980 & 1.6897 & 1.6436 & 0.0917 & -0.0545 & 0.0784 & 0.0183 & 0.1330 \\ DNN-MG(\(3+2\))-750 & 1.5591 & 1.7051 & 1.6376 & 0.1460 & -**0.0591** & **0.0834** & **0.0096** & **0.1426** \\ \hline \hline \end{tabular}
\end{table}
Table 7: The min, max, mean and amplitude of drag and lift functionals in the same setting as above in Figure 13. The best results in comparison to the MG(5) solution are highlighted with color coding: 1st best, 2nd best, 3rd best.

Figure 13: The drag (left) and lift (right) functionals for the channel with two obstacles at \(\mathrm{Re}=180\) are shown. The results are presented for three discretizations: The coarse mesh MG(3), the reference MG(5), and the DNN-MG(\(3+2\)) improved by a neural network.
\begin{table}
\begin{tabular}{l c c c c c c c c} \hline \hline Type & \(\min C_{d}\) & \(\max C_{d}\) & \(\widehat{C_{d}}\) & amp \(C_{d}\) & \(\min C_{l}\) & \(\max C_{l}\) & \(\widehat{C_{l}}\) & amp \(C_{l}\) \\ \hline MG(3) & 0.7674 & 1.3213 & 1.2176 & 0.5539 & -0.0087 & 0.0313 & **0.0097** & 0.0400 \\ MG(4) & 0.7905 & 1.5601 & 0.7696 & 1.3747 & -0.0091 & 0.0346 & 0.0061 & 0.0404 \\ MG(5) & 1.5090 & 1.7907 & 1.6642 & 0.2818 & -0.1577 & 0.1693 & **0.0083** & 0.3270 \\ DNN-MG(3 + 1)-512 & 1.3311 & 1.6166 & 1.3971 & 0.2855 & -0.0505 & 0.0960 & 0.0291 & 0.1466 \\ DNN-MG(3 + 2)-512 & 1.4353 & 1.5870 & 1.5126 & 0.1517 & -0.1398 & 0.0790 & -0.0528 & 0.2188 \\ DNN-MG(3 + 1)-750 & 1.4527 & 1.6565 & 1.5547 & 0.2038 & -0.0832 & 0.1507 & 0.0219 & 0.2339 \\ DNN-MG(3 + 2)-750 & **1.4644** & **1.6823** & **1.5929** & 0.2179 & -0.1325 & **0.1684** & 0.0279 & **0.3009** \\ \hline \hline \end{tabular}
\end{table}
Table 8: The min, max, mean and amplitude of drag and lift functionals in the same setting as above in Figure 14. The best results in comparison to the MG(5) solution are highlighted with color coding: 1st best, 2nd best, 3rd best.

Figure 14: The drag (left) and lift (right) functionals for the channel with two obstacles at \(\mathrm{Re}=240\) are shown. The results are presented for three discretizations: The coarse mesh MG(3), the reference MG(5), and the DNN-MG(3 + 2) improved by a neural network.
\begin{table}
\begin{tabular}{l c c c c c c c c} \hline \hline Type & \(\min C_{d}\) & \(\max C_{d}\) & \(\widehat{C_{d}}\) & amp \(C_{d}\) & \(\min C_{l}\) & \(\max C_{l}\) & \(\widehat{C_{l}}\) & amp \(C_{l}\) \\ \hline MG(3) & 0.7674 & 1.3213 & 1.2176 & 0.5539 & -0.0087 & 0.0313 & **0.0097** & 0.0400 \\ MG(4) & 0.7905 & 1.5601 & 0.7696 & 1.3747 & -0.0091 & 0.0346 & 0.0061 & 0.0404 \\ MG(5) & 1.5092 & 1.7908 & 1.6568 & 0.2816 & -0.0871 & 0.1354 & 0.0122 & 0.2225 \\ DNN-MG(3 + 1)-512 & 1.3049 & 1.5990 & 1.5111 & 0.2941 & -0.0471 & 0.1589 & 0.0317 & 0.2060 \\ DNN-MG(3 + 2)-512 & 1.3021 & 1.5548 & 1.3596 & 0.2526 & -0.1398 & 0.0790 & -0.0528 & 0.2188 \\ DNN-MG(3 + 1)-750 & 1.4050 & 1.6214 & 1.4465 & 0.2164 & -0.0684 & 0.2384 & 0.1214 & 0.3068 \\ DNN-MG(3 + 2)-750 & **1.4387** & **1.6319** & **1.5499** & 0.1932 & -**0.0957** & **0.1343** & **0.0067** & **0.2300** \\ \hline \hline \end{tabular}
\end{table}
Table 9: The min, max, mean and amplitude of drag and lift functionals in the same setting as above in Figure 15. The best results in comparison to the MG(5) solution are highlighted with color coding: 1st best, 2nd best, 3rd best.

Figure 15: The drag (left) and lift (right) functionals for the channel with two obstacles at \(\mathrm{Re}=300\) are shown. The results are presented for three discretizations: The coarse mesh MG(3), the reference MG(5), and the DNN-MG(3 + 2) improved by a neural network.
Figure 16: Relative velocity errors (first row) and pressure errors (second row) for coarse solutions with and without ANN correction (MG(3) and DNN-MG(\(3+2)-750\)) compared to the reference solution on level \(L+2\) for different Reynolds numbers. As in the case with a single obstacle, we show the errors of the MG(4) solution as well, which we again outperform.
Channel with two obstacles at \(\mathbf{Re}=300\)For Reynolds number \(300\), Figure 15 and Table 9 present the drag and lift forces predicted by DNN-MG. The results demonstrate the method's ability to capture the complex flow characteristics at higher Reynolds numbers. We see that DNN-MG generalizes well to different geometries and Reynolds numbers. Considering the scenarios we tested so far, we determine that DNN-MG indeed offers great potential in terms of computational efficiency.

Figure 16 demonstrates the good performance of DNN-MG in terms of velocity and pressure error. The figure shows a great improvement of the accuracy over an MG(3) solution and remarkably even outperforms the MG(4) solution. This further substantiates the efficiency in terms of runtime we demonstrated in Section 6.6. Overall, DNN-MG is robust with respect to changes of the geometry and material parameters and is able to improve solutions to PDEs in terms of characteristic quantities of interest and error measures.

## 7 Conclusion

We have presented the deep neural network multigrid solver (DNN-MG) which uses a deep neural network to improve the efficiency of a geometric multigrid solver, e. g. for the simulation of the Navier-Stokes equations. Previously demonstrated for 2D simulations, we extended DNN-MG to 3D and reformulated it in a rigorous manner. Despite the increased complexity of 3D flows, the algorithm remained applicable and delivered even larger speed-ups compared to 2D, while retaining the efficiency, generalizability, and scalability.

We established the efficacy of DNN-MG for large-scale simulations in regimes where direct solvers are not feasible anymore. The overhead is small and in 3D we can accelerate high fidelity simulations by a factor of \(35\). This comes at a trade-off in terms of quality compared to the high fidelity solution. However, we were able to consistently outperform classical solutions on the intermediate levels, both in terms of solution quality and wall time. This efficiency-performance trade-off establishes DNN-MG as a highly promising approach for accelerating numerical solution methods for PDEs.

We ascertained a great generalization capability of DNN-MG due to the local approach. We only trained the network with a single flow scenario at \(\mathbf{Re}=200\) and got great results for lower and higher Reynolds numbers. Our results showed that DNN-MG substantially improves the solution accuracy, as measured by the \(l^{2}\)-errors, and reduces the errors of the drag and lift functionals across all tested scenarios. The experiments with a two obstacle flow at different Reynolds numbers were particularly successful.

Especially noteworthy is that DNN-MG is able to predict flow profiles of completely different dynamics. The grid finenesses MG(3) and MG(5) considered by us produce flows with very different character. Thus, especially in the test problem with two obstacles, there are hardly any oscillations on the coarse grid. These are correctly reproduced by DNN-MG. The local approach is thus able to identify global structures of the Navier-Stokes solution and to correct the solution in both local and temporal dynamics.

Although DNN-MG generalizes well, there are limits to this approach. In future work we want to develop an online learning approach to retrain the network adaptively based on the uncertainty of the predictions. Further, we want to incorporate physical information, by including the residual of the PDE into the loss function. We further want to investigate the application of DNN-MG to other PDEs.

### Acknowledgement

NM acknowledges support by the Helmholtz-Gesellschaft grant number HIDSS-0002 DASHH. Computational resources (HSUper) were provided by the project hpc.bw, funded by dtec.bw -- Digitalization and
Technology Research Center of the Bundeswehr. dtec.bw is funded by the European Union - NextGenerationEU. The work of TR was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - 314838170, GRK 2297 MathCoRe.
