{'fraction_non_alphanumeric': 0.03193938231344715, 'fraction_numerical': 0.012353731056972952, 'mean_word_length': 4.362168278132072, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A neural network multigrid solver for the Navier-Stokes equations

We present the deep neural network multigrid solver (DNN-MG) that we develop for the instationary Navier-Stokes equations. DNN-MG improves computational efficiency using a judicious combination of a geometric multigrid solver and a recurrent neural network with memory. DNN-MG uses the multi-grid method to classically solve on coarse levels while the neural network corrects interpolated solutions on fine ones, thus avoiding the increasingly expensive computations that would have to be performed there. This results in a reduction in computation time through DNN-MG's highly compact neural network. The compactness results from its design for local patches and the available coarse multigrid solutions that provides a"guide"for the corrections. A compact neural network with a small number of parameters also reduces training time and data. Furthermore, the network's locality facilitates generalizability and allows one to use DNN-MG trained on one mesh domain also on different ones. We demonstrate the efficacy of DNN-MG for variations of the 2D laminar flow around an obstacle. For these, our method significantly improves the solutions as well as lift and drag functionals while requiring only about half the computation time of a full multigrid solution. We also show that DNN-MG trained for the configuration with one obstacle can be generalized to other time dependent problems that can be solved efficiently using a geometric multigrid method.


Introduction
The last decade has seen enormous progress on the use of deep neural networks for applications in machine translation, computer vision, and many other fields, e.g. [33,47]. These advances lead to a growing (and renewed) interest to apply neural networks also to problems in computational science and engineering, including for the simulation of partial differential equations. There, the existence of efficient numerical methods raises the questions how neural networks can be combined with, e.g. multigrid finite element methods, to leverage the benefits the different methodologies provide.
Towards this end, we propose the deep neural network multigrid solver (DNN-MG) that we develop in this paper for the simulation of the instationary Navier-Stokes equations. DNN-MG combines a geometric multigrid solver with a recurrent neural network to replace the computations on one or multiple finest mesh levels with a network-based correction. For this, the Navier-Stokes equations are first solved on the first L levels of the multigrid mesh hierarchy using a classical multigrid solver. The obtained solution is then interpolated (or prolongated) to a finer level L + K where the network predicts a correction of the velocity vector field towards the unknown ground truth. The right hand side of the problem is then formed on the fine level L + K and its restriction to level L is used in the classical multigrid solve in the next time step. Through this, the neural network-based corrections feed back to the coarse levels and affect the overall time evolution of the flow while the costs of the multigrid computations on fine levels are avoided.
The key to the efficiency of the DNN-MG solver is the use of a highly compact neural network with a small number of parameters. The compactness is achieved by the network's locality, i.e. that it operates independently on patches of the mesh domain, e.g. a single mesh element on level L or a collection of few adjacent elements. It is also enabled by the availability of the coarse multigrid solution that provides significant information about the flow behavior and thus reduces the state that needs to be carried in the network. A local, patch-based network also reduces the required amount of training data and time, since even a single flow contains already a multitude of patches and different behaviors. The latter also facilitates the networks' ability to generalize to flows not seen during training. To ensure that the neural network can model complex flow behavior, it contains a Gated Recurrent Unit (GRU) with memory that incorporates a flow's past into a prediction, similar to how many classical time integration schemes use multiple past time steps for their update. The GRU's memory and the use of the network's prediction in the right hand side at time step n + 1 facilitates a stable and temporally consistent time integration.
We demonstrate the efficacy of the DNN-MG solver for a 2D laminar flow in a channel with obstacles with varying elliptical eccentricities. We first consider the classical flow with one obstacle and train with different eccentricities than used for testing. 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 level L while requiring only about half the computation time of a full solution on level L + 1. For a prediction using level L + 2 even better results are obtained with only a very small additional overhead compared to DNN-MG on level L+1. To analyze DNN-MG's ability to generalize to flow regime's not seen during training, we also consider the channel flow with no or with two obstacles while using the neural network trained for the one obstacle case. In both instances, DNN-MG is able to provide high fidelity solutions and for the two obstacle flow it again reduces the error of the lift and drag functionals substantially.
The remainder of the paper is structured as follows. In the next section we review related work that uses neural networks for the simulation of partial differential equations. In Sec. 3 we then recall the solution of the Navier-Stokes equations using the geometric multigrid method before providing an introduction to recurrent neural networks in Sec. 4. In Sec. 5 we introduce the deep neural network multigrid solver, discuss its design and also sketch its generalization to other problems. Our numerical results are presented in Sec. 6.

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 in particular the Navier-Stokes equations.
For learning, partial differential equations are often considered as sequential time series. Neural network architectures developed for these data can then be employed. The most common architecture are recurrent neural networks (RNN), e.g. with Long-Short-Term-Memory units (LSTM) [26] to avoid the vanishing gradient problem. Because of the prevalence of sequential data, LSTMs led to many variations over the years. One are Gated Recurrent Units (GRUs) which, while not as powerful as LSTMs, provide often similar performance and are easier to train [13]. An alternative to LSTM-type architectures are temporal convolutional neural networks (TCNs) where the temporal dependence is modeled using a causal (i.e. one-sided) convolution in the time variable. Bai et al. [1] demonstrated that this can provide substantial improvements in prediction accuracy over recurrent neural networks with memory. A different approach to modeling temporal dependencies was developed by Voelker et al. [50] who represented the time domain in Legendre polynomials and learned basis function coefficients. Through this, they obtained a neural network model with continuous time. A related idea are the neural ordinary differential equations architectures by Chen et al. [11] where the depth of the network is continuous.
Approaches for a direct description of the time evolution of partial differential equations using neural networks go back more than 20 years, e.g. [32]. In recent work, Raissi et al. [42] demonstrated that the velocity field 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 [56], and Raissi and Karniadakis [41] 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. [30] recently demonstrated neural network-based simulations for a broad range of applications, including fluid dynamics, by also optimizing the network architecture itself during the training process. Eichinger, Heinlein and Klawonn [34] 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 neural network-based approaches to the simulation of partial differential equations, there have been different attempts to integrate neural networks into existing numerical techniques. As in our work, the objective is 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 [48,7,44]. These were then used in classical simulations, e.g. based on finite elements. Wiewel et al. [54] presented a simulation of the Navier-Stokes equation where an LSTM-based deep neural network is used to predict the pressure correction. Xi et al. [55] and Werhahn et al. [53] use neural networks to upsample a coarse simulation and obtain physically plausible high frequency details. An approach based on the splitting method for high dimensional PDEs was developed by Beck et al. [3] where one obtains a set of smaller learning problems that are easier to treat than a single large one. Wan and Sapsis [51] describe the motion of finite size particles in fluid flows by combining analytic arguments with a neural network that models the discrepancy between the (highly) idealized analytic model and real world observations.
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 in the literature. The authors of [37,12,28] developed neural networks that incorporate the symplectic structure of Hamiltonian mechanics and they demonstrate that this improves generalization and accuracy. In a related line of research, ideas from nonlinear analysis and stability theory were exploited. Erichson et al. [16], for example, improve the accuracy of the learning process by ensuring that the learned system has the same Lyapunov stability as the system of interest. Similarly, Gin et al. [20] train a neural network to provide a coordinate transformation that linearizes a nonlinear PDE, in the spirit of the Cole-Hopf transform for Burgers equation, but using a neural network to find and represent the transform. For this, they use a Koopman operator representation. Related to this is the work Zhang et al. [57] where a neural network is used to improve the dynamically orthogonal (DO) and bi-orthogonal (BO) methods for stochastic pdes.
An alternative line of research aims to learn surrogate neural network representations of physical systems from the governing (stochastic) partial differential equations without recourse to simulations or measurements [19,58,29]. This avoids the training data generation that is often a bottleneck with the approaches above and allows one to quantify the uncertainty of the prediction. Karumuri et al. [29], e.g., obtain their loss function through a variational formulation of the stochastic PDEs they are considering. Related to this is the Deep Ritz Method by E and Yu [15] where the ansatz of the Ritz method is used but with neural networks for the function approximation. Qin et al. [40] show that residual networks can be seen as exact one-step time stepping methods and these provide building blocks for multi-step methods.
The performance of (deep) neural networks has also been analyzed theoretically. Already in the 1980s and 1990s, the first results demonstrated the potential of neural networks as universal function approximators [14,27] that do not suffer from the curse of dimensionality [2]. More recent results on the approximation properties of deep neural networks focused in particular on the required size of a network to achieve certain approximation bound, e.g. [8]. Mallat [35] provided insight into the effectiveness of deep neural networks using tools from harmonic analysis and pointed out the importance of invariants for their understanding. For partial differential equations, currently only few results exist. Kutyniok et al. [31] studied the theoretical power of forward, deep neural networks for reduced order models and they show that the solution map, which yields the basis coefficients of the reduced basis for given values of the parameters, can be efficiently approximated using neural networks and in a manner essentially independent of the dimension. Grohs et al. [21] study the spacetime error for PDEs when Euler steps are realized using neural networks. To our knowledge, for partial differential equations no theoretical results for recurrent neural networks with memory exist.

Finite element discretization of the Navier-Stokes equations
In the following, we provide an overview of the solution of the Navier-Stokes equations using the geometric multigrid method.
We consider a domain Ω ∈ R d with d ∈ {2, 3} and Lipschitz-continuous boundary Γ and a bounded time interval [0, T ]. The solution of the instationary Navier-Stokes equations then seeks the velocity v : [0, T ] × Ω → R d and pressure p : where Re > 0 is the Reynolds number and f an external force. The initial and boundary conditions are given by where n denotes the outward facing unit normal vector on the boundary ∂Ω of the domain. The boundary Γ = Γ D ∪ Γ N is split into subsets Γ D with Dirichlet boundary conditions and Γ N with Neumann type conditions.

Variational formulation of the Navier-Stokes equations
Eq. 1 can be discretized with finite elements. We briefly sketch the variational formulation in the following. By L 2 (Ω) we denote the space of square integrable functions on the domain Ω ⊂ R d with scalar product (·, ·) and by H 1 (Ω) those L 2 (Ω) functions with weak first derivative in L 2 (Ω). The function spaces for the velocity and pressure are then where v D ∈ H 1 (Ω) d is an extension of the Dirichlet data on Γ D into the domain. Given Dirichlet data on the complete boundary, the pressure is normalized to yield uniqueness. With these spaces, the variational formulation of Eq. 1 is given by Let Ω h be a quadrilateral or hexahedral finite element mesh of the domain Ω satisfying the usual requirements on the structural and form regularity such that the standard interpolation results hold, compare [43,Section 4.2]. For the finite element discretization of Eq. 4, we choose W (r) h as the space of continuous functions which are polynomials of degree r on each element T ∈ Ω h . If the finite element mesh contains elements of generalized quadrilateral or hexahedral type, e.g. to resolve curved boundaries, we use isoparametric finite elements, see [43, Section 4.2.1] for details and the specific realization chosen for this work. The discrete trial-and test-spaces for the discretization of Eq. 4 we throughout the paper be v h , φ h ∈ V h = [W (2) h ] d and p h , ξ h ∈ L h = W (2) h . Since the resulting equal order finite element pair V h × L h does not fulfill the inf-sup condition, we add additional stabilization terms of local projection type [4]. The resulting semidiscrete variational problem then reads where π h : W denotes the local projection operator into the space of linear polynomials. The local stabilization parameter α T is chosen as where we usually use α 0 = 0.1 and where h T is the local mesh size of element T . We refer to [4,9] for details. Since the numerical experiments described in Section 6 consider moderate Reynolds numbers only, no stabilization of the convective term is required.

Time discretization
For temporal discretization, the time interval [0, T ] is split into discrete time steps of uniform size 0 = t 0 < t 1 < · · · < 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 v n := 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 Eq. 5, resulting in the fully discrete problem The right hand side only depends on the the velocity v n−1 at the last time step n − 1 and we will denote it as b n−1 in the following. The Crank-Nicolson time discretization in Eq. 6 corresponds to a nonlinear system of equations that can be solved, e.g., with the Newton-Krylov method. In the presented form it is sufficiently robust for smooth initial data v 0 ; we refer to [23] for small modifications with improved robustness and stability.

Newton-Krylov solution
In each time step, Eq. 6 amounts to a large nonlinear system of algebraic equations. By introducing the unknown x = (v n , p n ), it can be written as for all test functions φ i h and ξ i h . The nonlinear problem is solved by Newton's method based on the initial guess x (0) = (v n−1 , p n−1 ) and the iteration l = 1, 2, . . .
Algorithm 1 The geometric multigrid method for the solution of the linear system A L x L = b L given on a finest level L. If all high frequency errors of S are smoothed at a constant rate, the method achieves optimal complexity O(n). The multigrid solver is initiated on the finest mesh level L and then used recursively: if l − 1 = 0 then Coarse-grid solution 6:

end if
10: x l ← s l + P(l, c l−1 ) Prolongation 11: return s l 13: end procedure where we denote by A (x (l−1) ) the Jacobian of A at x (l−1) , i.e. the matrix of first partial derivatives. For the finite element discretization of the Navier-Stokes equations this is easily computed analytically, cf. [43,Sec. 4.4.2].
Each Newton step requires the solution of a linear system of equations, 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 with optimal robustness, we employ the generalized minimal residual method (GMRES) by Saad [45]. The convergence is thereby accelerated through a geometric multigrid solver that is used as a preconditioner with a single sweep performed for every GMRES step. We usually require about 10 linear solves for each nonlinear Newton iteration.

The geometric multigrid method
Multigrid methods are a class of efficient algorithms for the solution of linear systems that can achieve the optimal complexity of O(n), where n is the number of degrees of freedom (in our case proportional to the number of mesh vertices). We employ the geometric multigrid method as a preconditioner for the GMRES iterations. It is summarized in Algo. 1.
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 solving the linear system on the finest mesh level L, as one would do with traditional solvers, geometric multigrid methods smooth high frequency errors with a simple iterative method and treat the remaining errors on lower levels. Usually, smoothing is applied at the beginning (Algo. 1, line 2) and at the end (Algo. 1, line 11) of each iteration. In between, the remaining residual is computed (Algo. 1, line 3) and restricted to the next coarse level L − 1 (Algo. 1, line 4). On this coarser level, the multigrid iteration is then called recursively (Algo. 1, line 8) until the coarsest level L − 1 = 0 has been reached where a direct solver is employed (Algo. 1, line 6). The update from the coarse mesh L − 1 is then prolongated back to level L (Algo. 1, line 10).
The mesh transfer from fine to coarse levels is accomplished with L 2 -projections, known as restrictions, and from coarse to fine ones with interpolations, known as prolongations. The smoothing that preceded the restrictions uses a smoothing operator S(A l , b l , x l ). It is usually a simple iteration that yields an approximation to the solution of the linear system, i.e.
and that aims to quickly reduce all high frequency components of the residual b l − A l x l . If this happens on all levels l at a constant rate, then the geometric multigrid method achieves the optimal complexity O(n). In our implementation we use a simple iteration of Vanka type [49], which allows for easy parallelization and gives very good performance with less than 5 pre-and post-smoothing steps [18,17]. The idea of the Vanka type smoother is to exactly solve small subproblems and to replace the inverse of where R T is the restriction to those nodes that belong to one mesh element T ∈ Ω l and R T T its transpose. Considering piecewise quadratic finite elements, 9 nodes are attached to each element such that the local matrices to be inverted have the dimension R T A l R T T ∈ R 27×27 , since each node comprises the scalar pressure and two velocity components.
For a more detailed discussion of the geometric multigrid method we refer to [22,10] and the realization in Gascoigne 3D [5] used throughout this work is described in [6].

Recurrent neural networks with memory
Different types of neural networks exist. The ones originally introduced under the name are today known as feed-forward neural networks and used for tasks such as image recognition or data mining. Recurrent neural networks are an extension where an activation variable a k n ∈ R p k allows to propagate information in a discrete time variable indexed by n, such as those occurring in time series (e.g. the stock market) or sequential data (e.g. text). An extension of this model uses network nodes with memory. These allow to more effectively model long term dependencies and also to avoid the vanishing gradient problem that otherwise occurs with recursive neural networks [25]. In our work, we use Gated Recurrent Units (GRUs) [13] that are somewhat simpler than the better known Long-Short-Term-Memory (LSTM) units [52] but are therefore also easier to train. At while not converged do Newton-GMRES method for Eq. 6 3: Set up rhs of Eq. 6 for next time step Restriction of rhs to level L 10: end for layer k in the network, a GRU is given by where denotes the element-wise product and the W (·) k are weight matrices and bias vectors determined by training. The so called update gate vector z k n ∈ R q in Eq. 11 controls to what extend h k n−1 is carried over to the output h k n at the current time n and the reset gate vector r k n ∈ R q controls the contribution of h k n−1 to the nonlinearity of the cell. Together, z k n ∈ R q and r k n ∈ R q thus control the memory of a GRU cell, i.e. to what extent the past hidden output h k n−1 contributes to the current one h k n .

A deep neural network multigrid solver
In this section, we develop the deep neural network multigrid solver (DNN-MG). It uses a recurrent neural network to predict the correction of a coarse mesh solution that has been prolongated (or interpolated) onto one or multiple mesh levels. Through this, we can obtain solutions that are more accurate than those obtained with the coarse mesh only while being computationally more efficient than performing the full (multigrid) computations on the fine mesh level(s). We will develop the DNN-MG solver in the context of the Navier-Stokes equations and return to the general formulation at the end of the section. To simplify the exposition, we will also assume that the neural network operates on only one finer level L + 1. However, bridging more levels is possible as we demonstrate in Sec. 6.6. Figure 1: The neural network part of DNN-MG predicts a correction d L+1 n of the prolongated velocity P(v L n ) towards the unknown ground truthv L+1 n , which traditionally would be obtained using the multigrid solver for level L + 1. The correction d L+1 n is incorporated into the time integration by using the improved solutionṽ L+1 n + d L+1 n to build a right hand side b L+1 n+1 and restricting it to level L. The restriction b L n+1 is then used in the next time step in the multigrid solver. Temporal coherence of the correction (and hence the overall flow) is achieved through the memory of the GRU in the neural network.

Time stepping using DNN-MG
We begin by detailing one time step of the simulation of the Navier-Stokes equations using the DNN-MG solver. The computations are summarized in Algo. 2 and a conceptual depiction is provided in Fig. 1.
At the beginning of time step n, we first solve for the unknown velocity v L n and pressure p L n on the coarse level L using the classical Newton-Krylow simulation (Algo 2, lines 2-5) as described in Sec The training inputs are velocity fields on level L and L + 1, with the latter one serving as reference solution. From these as well as the meshes Ω L and Ω L+1 on the respective levels we derive the input to the neural network that operates patch-wise (cf. Fig. 3). At application time, we use a classical multigrid solution on level L and apply the neural network part of DNN-MG to obtained a corrected solution P(v L n )+d L+1 n on level L+1. With it, a right hand side b L+1 n is computed that is then restricted to level L where it is used in the multigrid solver in the next time step.

The Neural Network of DNN-MG
At the heart of the DNN-MG solver is its neural network component. Care in its design is required to ensure that it improves the efficiency of the runtime computations, is easy to train, and generalizes well to flow regimes not seen during training. These objectives will be attained through a very compact neural network with a small number of parameters and a local, patch-based design. An overview of the network component is provided in Fig. 2. behavior on the entire simulation domain, necessitating substantial training time and data. The domain would also be "baked" into the network, hampering generalization.
To avoid these issues, DNN-MG uses a neural network that operates patch-wise over small neighborhoods of the mesh. In its most local variant each patch consists of the unknowns that belong to one single mesh element only, cf. Fig. 3; a more global approach is easily realized by considering larger patches consisting of multiple adjacent elements. Our network takes as input thus only information from one patch and predicts the velocity update d L+1 n only for the mesh nodes in the patch, with the prediction done independently for all patches to cover the entire simulation domain (duplicate predictions for adjacent patches are averaged). The locality is one key to the compactness of the neural network and hence also to its efficient evaluation. Furthermore, since a single neural network is used for all patches in the domain, training with a small number of flows exposes the network to a large number of different flow behaviors. This facilitates the network's ability to generalize to flows not seen during training and reduces training time and data.
Neural network inputs The efficacy of the neural network prediction relies critically on suitable network inputs that provide rich information about the coarse flow behavior on level L as well as the local patch geometry (similar to how the geometry of mesh cells play an important role in classical finite element simulations). In particular, with carefully selected inputs one can potentially reduce the number of parameters in the neural network, which reduces the evaluation time at runtime and the amount of data and time required for training.
When the most localized patches directly supported by the mesh are used (as we do in most of our numerical experiments) each patch consists of exactly one mesh element on level L with N L = 4 nodes and thus with N L+1 = 9 nodes on level L + 1. Our inputs to the neural network are then: • nonlinear residuals r L+1 n = f L+1 − A L+1 (P(x L n )) ∈ R 3 N L+1 of the prolongated coarse mesh solution for Eq. 7; • the velocities v L n ∈ R 2N L on mesh level L; • Péclet numbers Pe L ∈ R N L ; • geometry of the cells, in particular the edge lengths h c ∈ R 4 , the aspect ratios a c ∈ R 4 (of two neighboring and two opposite edges each), and the angles between the edges α c ∈ R 4 .
We hence have in total 51 inputs to the neural network. The most important input is the nonlinear residual r L+1 n = f L+1 − A L+1 (P(x L n )), cf. Eq. 7, that provides a measure of the local error (or defect) at every mesh node on level L + 1 and thereby incorporates information about the partial differential equation. We also use the Péclet number Pe L = qṽ n L /ν where ν is the fluid's viscosity and q a characteristic length scale, which is the patch size in our case. The Péclet number describes the ratio between transport and diffusion on each node on level L and it hence measures which part of the Navier-Stokes equations dominates the local behavior of the flow.
Neural network architecture To predict a defect d L+1 n that is consistent with the fluid flow, we use recurrent neural networks with memory. Through this, the past simulation states . . . , v n−2 , v n−1 affect the predicted update d n at time step n.
The specific network architecture we use for DNN-MG is shown in Fig. 4. At its heart is a GRU unit with memory, cf. Sec 4. It could in principle be replaced by another cell with memory, such as an LSTM unit or a TCN, but since the GRU unit performed well in our experiments and is, through its simpler design, easier to train we used it in the present work. The two convolutional blocks in the network are designed to learn local dependencies, with one learning the interactions of the two velocity components at each vertex and the other the spatial dependence of each component across a patch. Each of these blocks consist of two layers, where the second one is just the transposed version of the first layer so that input and output sizes are equal. The final convolutional layer reduce the dimensionality of the data from the multiple filters that form the preceding convolutional layers to the output d n of the network.
Training of DNN-MG The training of DNN-MG is based on simulations of the Navier-Stokes equations for which a multigrid representation of the velocity v l n with two levels L and L + 1 is available. The velocity v L+1 n thereby serves as ground truthv L+1 n . The training then finds network weight (in the GRU unit, the dense and the convolutional layers) such that the norm (ṽ L+1 n + d L+1 n ) −v L+1 n of the difference between the predicted velocityṽ L+1 n + d L+1 n (i.e. the velocity after correction, Algo. 2, line 9) and the ground truthv L n is minimized over the course of a simulation.

A general formulation of DNN-MG
After the detailed description of DNN-MG for the solution of the Navier-Stokes equations we will in the following show how it can be applied to a general PDE in variational form u ∈ V : is a semi-linear form (linear in the second argument) and V is a Banach space of test functions for the problem. Given a suitable Galerkin formulation in the discrete subspace V h ⊂ V the abstract Newton iteration for the solution of the problem A h (x h ) = f h can be formulated exactly as for the Navier-Stokes equations, compare Eq. 7 and Eq. 8. The Newton-Krylov approach based on the preconditioned GMRES iteration and a geometric multigrid solver as preconditioner is highly robust and it can be applied to a variety of PDE models (as is done in Gascoigne [5]). The DNN-MG solver in Algorithm 2 can hence in principle be applied to various PDEs.

Algorithmic Complexity of DNN-MG
The Newton-Krylov geometric multigrid method can already achieve optimal complexity, which raises the question what advantage DNN-MG has to offer. Before addressing the question from a practical point of view in the next section, we sketch in the following an answer from a theoretical point of view. The Newton-Krylov geometric multigrid framework achieves linear complexity O(N ) in the number N of degrees of freedom, which roughly quadruples with each global mesh refinement, i.e. N L+1 ≈ 4N L . The constant hidden in 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 10 GMRES steps with one sweep of the geometric multigrid solver as preconditioning. Furthermore, in the geometric multigrid typically 5 pre-smoothing and 5 post-smoothing steps have to be determined. Hence, a total of about 500 Vanka smoothing steps must be performed on each of the mesh layers L, L − 1, . . . , 0. One Vanka step thereby requires the inversion of about O(N l /4) small matrices of size 27 × 27, one on each element of the mesh, cf. Eq. 9, resulting in approximately 27 3 ≈ 20 000 operations. Since the complexity of all mesh levels sums up to about N L + N L−1 + . . . N 0 ≈ N L (1 + 4 −1 + · · · + 4 −L ) ≈ 4 3 N L we can estimate the effort for the complete solution process on level L as 5 · 10 · (5 + 5) · 27 3 · 4 3 ≈ 10 7 N L . We thereby only counted the effort for smoothing and neglected all further matrix, vector and scalar products. Solving the problem on an additional mesh level L+1 would require about four times this effort, i.e. ≈ 4 · 10 7 N L operations.
On the other hand, the DNN-MG method only requires the prolongation of the solution to the next finer mesh and there, as main effort, 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 O(N L+1 /4) = 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 c inputs. The DNN-MG approach asks for only one such evalution of the network in each time step. The specific network model used for the numerical test cases comprises 8634 trainable parameters. Hence, the effort for correcting the level L Newton-Krylov solution by the neural network on level L + 1 can be estimated by 8634N L , which is negligible in comparison to the effort of solving on level L + 1 (≈ 40 7 N L ). In this section, we present numerical results obtained with the deep neural network multigrid solver (DNN-MG) for the Navier-Stokes equations. We consider a classical benchmark problem describing the laminar flow around a cylinder, cf. [46]. In contrast to the original benchmark configuration, we parametrize the geometry with an ellipsoid with varying aspect ratio as obstacle. This allows us to easily generate different but similar training data and also provides a basic test for the generalization ability of our technique. We also report on the runtime performance of DNN-MG and study its ability to generalize to flow regimes that differ from the training data. We will work with DNN-MG with the classical simulation operating up to level L and the neural network on level unknowns multigrid levels coarse 2124 3 (L) fine 8088 4 (L + 1) reference 31536 5 (L + 2) Table 1: Spatial meshes with increasing number of multigrid levels and unknowns.

Numerical examples
L + 1 except for Sec. 6.6 where we will consider prediction across two levels and with larger patches.

Setup
We consider a variation of the well-established "laminar flow around a cylinder" introduced by Schäfer and Turek [46], in particular the unsteady 2D-2 testcase where the laminar flow around a circular obstacle at Reynolds number Re = 100 leads to a stable periodic flow pattern. The geometry of the flow domain is shown in Fig 5. The flow is driven by a Dirichlet profile v = v D at the left boundary Γ in given by where H = 0.41 is the height of the channel and the function ω(t) regularizes the startup phase of the flow. On the wall boundary Γ wall and on the obstacle we prescribe no-slip boundary conditions v = 0 and on the outflow boundary Γ out we use a do-nothing outflow condition [24]. The average flow rate is v avg = 4 3 which, through the different obstacles used for training and testing, results in the Reynolds numbers Re train = 166 and Re test = 200, respectively.
In all numerical examples the temporal step size is k = 0.01. In space, we consider the sequence of meshes in Table 1. Finer levels result from a uniform refinement of the "coarse" mesh conforming to the obstacle (see [5] for details on conforming boundary handling).

Neural network parameters
As discussed in Sec. 5, DNN-MG uses a local neural network that operates over patches of the mesh domain. For our numerical experiments, we defined a patch to be the refinement of a coarse mesh cell on level L, i. e. we used the most localized patches that are naturally supported by the mesh structure. This resulted in a network (cf. Figure 4) with only 8634 trainable parameters. With none of the layers having a bias and N in = 51 (the input length), N GRU = 32 (the size of the GRU state), N out = 18 (the output length), there are 3(32 · 32 + 32 · 51) = 7968 parameters in the GRU and 32 · 18 = 576 in the fully connected layer. Each of the convolutional blocks that learns vertex interactions has, furthermore, 36 parameters and the last convolutional layer has 18 parameters. Note that the length of the feature vector (input length) only depends on the number of elements that make up a patch but not on the number of elements or patches in the mesh. This is due to our local patch-based design and in particular, the network can be used independently of the number of elements, refinement level of the mesh, the domain Ω and the overall flow situation. However, changing these parameters too much could affect the performance of the network. This is closely linked to the generalizability, which we demonstrate in Section 6.7.
Each convolutional block learns 4 filters, in the first layer the number of filters is expanded to 4 and in the following transposed layer it is reduced to 1. In the last layer the outputs of the convolutional blocks are concatenated along the filter dimension, which then has to be reduced to a length of 1. This is achieved by a a convolution with length 1 (summing along the filter dimension and weighting the output), resulting in 18 trainable parameters.

Neural network training
As training data for DNN-MG we used three (L + 1)-level multi-grid simulations of the flow around a cylinder, each with a slightly different elliptical obstacle with height 0.1. The small number of simulations was sufficient since, by the patch-wise application of the network, a single simulation provides N c ×N T training items, where N c is the number of patches and N T is the length of the time series. By the periodicity of the flow it also sufficed to consider a small number of full periods. We hence used t ∈ [1,6] resulting in N T = 500 (we used t > 1 to avoid the start-up phase of the simulation where the quasi-periodic flow forms).
As loss function L we employed a simple l 2 -loss Tikhonov-regularization was used with a scaling factor α = 10 −4 . Training was performed with the ADAM optimizer with a limit of 1000 epochs. At the end, the model with the lowest validation loss was saved. This resulted in a training time of 4h 23m on an Intel Xeon Gold 6254.

Flow prediction
For testing we used the flow around an elliptic obstacle with an increased height of 0.12 compared to the training data (see Fig. 5). Fig. 6 shows the velocity magnitude at time Re ∇v − pI n · e 2 ds, respectively, where e 1 = (1, 0) T and e 2 = (0, 1) T . The results in Fig. 8 and Table 2 demonstrate that DNN-MG is able to substantially correct both the drag and lift functionals. Maximum and minimum values but also the frequency of oscillation of the augmented simulation are significantly closer to the fine mesh one that serves as reference. Fig. 10, finally, shows the relative velocity and pressure errors  Table 1 reveals a frequency shift which is depending on the spatial finite element discretization. We show the startup phase [0, 3] and the interval [9,10], where the flow is fully developed.
where · 2 denotes the Euclidean norm of the discrete solution. DNN-MG is able to robustly reduce the velocity error from 20% − 27% to about 10%. While the relative pressure error is also reduced by a factor of 2 on average, some oscillatory peaks remain. These stem from the shift in periodicity as discussed in the following remark.
Remark. Fig. 7 shows the lift functional for classical multigrid simulations with three different resolutions both in the startup phase t ∈ [0, 3] and for t ∈ [9,10] where the flow profile is fully developed. It can be seen that on coarser meshes it takes longer for the periodic flow profile to develop (Fig. 7, left) and the period of oscillations is getting shorter with increasing mesh resolution (Fig. 7, right). This shift in frequency and phase as a function of mesh resolution prevents a direct comparison of results obtained on different levels in a multigrid-type algorithm such as DNN-MG. In [36], a shifted variant of the Crank-Nicolson scheme is discussed to remedy the problem. Since this cannot be used directly for DNN-MG, we instead compare the functional values in shifted intervals [t * , t * + 1], where t * is chosen as the first peak in the lift functional after time t = 9s, identified separately for the coarse mesh solution, DNN-MG and the fine mesh one on level L + 1. The results reported in Fig. 8, Fig. 10, Fig. 14 and Table 2, Table 3 are those for the temporally adjusted flow in [t * , t * + 1].

Timings
The use of the DNN-MG method in applications is only meaningful when it not only improves the accuracy compared to a coarse mesh solution (as demonstrated in the last subsection) but also reduces the computation time compared to solving on a finer mesh. As shown in Figure 11, DNN-MG saves 55 % of the wall clock time required for a fine mesh solution. It hence indeed improves the computational efficiency. Compared to solving on a coarse mesh, the runtime increases by 37 %, but with an improved accuracy for DNN-MG. Fig. 11 also shows that the evaluation of the neural network requires

Prediction across 2 levels
So far we only considered the prediction of one level, i. e. a correction on level L w. r. t. level L + 1. In this section we investigate the prediction of a correction using level L + 2.
The DNN-MG algorithm does not change conceptually in this case. Prolongation and restriction, however, now transfer between L and L + 2. One patch corresponds again to one element on the mesh of level L but this now includes 4 × 4 equally subdivided elements of the mesh on level L + 2. A small modification to the neural network is necessary to accommodate the larger input and output sizes of 99 and 50, respectively. We also increase its internal size to 60 960 trainable parameters by setting N GRU = 96 to account for the larger number of degrees of freedom. Figure 9 and Table 2 show the obtained values for the drag and lift functionals. When comparing to the drag and lift for the L + 1 case in Fig. 8, it must be noted that different data is used for training in the two cases. A comparison should thus consider whether the  Figure 11: Comparison of the wall-clock times of the different mehods DNN-MG method is able to learn the respective data (i.e. L + 1 and L + 2, respectively) and not whether the total error can be reduced. Figure 9 and Table 2 show that with the larger mesh and newtork in the L + 2 case one obtains predictions for the drag and in particular the lift that are of very high quality and significantly improved compared to the predictions for L + 1. For L + 2, minimum and maximum of the lift have an error of less than 1% whereas the lift error was about 15% in the L + 1 case. Further, the enhanced frequency and amplitude of the lift functional is especially notable.
The computation time of DNN-MG with a two level prediction is very similar to the one level case. The additional prolongation and restriction computations do not add noticeable workload. Copying the larger input and output data and evaluating the larger network also adds less than 1% to the runtime compared to DNN-MG (L + 1). However, we need about 50% more time to generate the training data for the two level predictions and the training times increases by ≈ 5%.

Generalizability
For the practicality of the DNN-MG method it is important that a neural network trained on one flow performs well also on similar ones, i.e. that it generalizes to flows beyond those seen during training. The results in Sec. 6.4 already demonstrated that the network is able to do so under small perturbations of the geometry of the obstacle. In the following, we consider three more substantial changes to the flow: a channel without an obstacle, a channel with two obstacles, and a flow in an L-shaped domain, Unless stated otherwise, we reuse the network trained as discussed in Sec. 6.3 for the channel with a single obstacle. In this section we again use level L as coarse resolution and L + 1 as fine one on which the neural network component of DNN-MG operates.
Channel flow without obstacle In our first test case we considered the stationary, 1dimensional flow that develops in the channel of Fig. 5 when the obstacle close to the inflow is removed. The stationary flow is challenging for DNN-MG in our setup since it is trained on oscillatory flow data only. Fig. 12 shows the DNN-MG solution that has been obtained. We observe a slight deformation of the velocity field (which should be strictly constant in the horizontal direction) but an otherwise stable and in particular stationary flow. The network is thus able to generalize to this entirely different flow regime. We conjecture that the network's ability to do so results from our network design that learns individual patches. Hence, even for the non-stationary flow around a cylinder that served as training data one has at least locally also relatively stationary flows in the training data.
Channel with two obstacles The second generalization test was also an extension of the initial benchmark problem in Sec. 6.4 where we replaced the one obstacle with two. These are centered at (0.2, 0.18) and (0.6, 0.24), respectively, with their minor and major axes being 0.9 and 1.1. Since the incoming flow is already disturbed by the first obstacle before reaching the second one, the configuration differs substantially from the training scenario with only one. To represent the modified geometry of this test case we employ a finite element mesh with 9, 741 degrees of freedom on level L and 37, 917 for the fine mesh of level L + 1, The magnitude of the velocity field obtained with DNN-MG as well as multigrid solutions on levels L and L + 1 are shown in Figure 13. All fields look similar, although visually DNN-MG is closer to the reference than the coarse solution in the vicinity of the two obstacles and differs more in the far field on the right. To quantitatively assess the performance of DNN-MG we again used the drag and lift functionals, which we now evaluated for the second obstacle. Fig. 14 and Table 3 show that DNN-MG is able to substantially improve the estimated amplitude, despite being trained with only one obstacle. Moreover, the results were obtained using a network trained on a different mesh, which was possible since our neural network operates on local patches that are agnostic to the global mesh structure. The frequency in Fig. 14 Table 3: Drag and the lift functionals for the generalization of DNN-MG trained on the flow around one cylinder for the channel with two obstacles. As reference we show the classical multi-grid solutions on the coarse and the fine meshes. L-shaped domain As a final test case we considered an L-shaped domain without obstacle as shown in Fig. 15. Using DNN-MG with the trained neural network from Sec. 6.4 fails in this case since the training data contains only flows in the positive xdirection (the predicted solution then has a unphysical drift in the x-direction also in the lower arm of the L-shape). To alleviate this problem, we retrained DNN-MG with an augmented data set where the existing data was also used with the x-and y-axes swapped. Through this, no new and expensive simulations on level L + 1 were necessary, and the training converged faster since the optimization could start from an already trained network. In Figure 15 the magnitude of the resulting velocity fields for DNN-MG and multigrid solutions on levels L and L + 1 are depicted. As can be seen there, DNN-MG does not provide a substantial improvement over the coarse multigrid solution. In our opinion, this is not surprising given that the training data does not contain a flow past a sharp edge. Nonetheless, the DNN-MG simulation remains stable and the neural network correction does not deteriorate the solution quality of the coarse simulation on level L.

Conclusion
We have presented the deep neural network multigrid solver (DNN-MG) that uses a recurrent neural network to improve the efficiency of a geometric multigrid solver, e.g. for the simulation of the Navier-Stokes equations. The neural network is integrated into the multigrid hierarchy and corrects on fine levels the prolongated solutions (instead of the smoothing operations used classically), affecting the overall time evolution through Figure 15: Comparison between simulations with L + 1 multigrid levels, DNN-MG, L multigrid levels (from left to right).
the right hand side of the problem. Central to DNN-MG is the compactness of its neural network. This is enabled through the network's locality, namely that it operates independently on small patches of the domain, and the coarse multigrid solution that provides a "guide" for the network's corrections. The compactness is vital for the computational efficiency of DNN-MG and also reduces the training data and time that are required. The locality furthermore facilitates generalizability and even allows one to use a network trained on one mesh domain with an entirely different one. The use of recurrent neural networks with memory thereby enables DNN-MG to account for time dependencies in a flow and obtain predicted corrections that are temporally coherent. We analyzed the efficacy of DNN-MG for the classical Navier-Stokes benchmark problem of a laminar flow around an obstacle. Our results show that DNN-MG is able to substantially improve the accuracy of solutions and provide significantly reduced errors for the drag and lift functionals. It thereby requires only approximately half the computation time of a multigrid solution at the same resolution and only a modest increase compared to a coarse one. At the same time, DNN-MG is able to generalize substantially beyond the flow used for training and with the neural network trained on the one obstacle flow we were able to obtain high fidelity solutions also for channel flows with no or two obstacles.
The experiment with the flow in the L-shaped domain, where we do not obtain improvements over a coarse multigrid solution, shows the limits of the current DNN-MG implementation. We believe that a more complex neural network and more variability in the training data will be able to alleviate the limitation. We want to study this in future work. There, we also want to generalize our current implementation to 3D flows. Although DNN-MG is already considering the physics of the underlying problem by using the nonlinear residual on level L + 1, it is worth exploring how additional model knowledge such as the divergence freedom of the predicted velocity fields can be integrated into DNN-MG. In future work, we also want to consider the application of DNN-MG to other partial differential equations, cf. from nonlinear elasticity, cf. Sec. 5.3.