Improving Pseudo-Time Stepping Convergence for CFD Simulations With
  Neural Networks
Computational fluid dynamics (CFD) simulations of viscous fluids described by
the Navier-Stokes equations are considered. Depending on the Reynolds number of
the flow, the Navier-Stokes equations may exhibit a highly nonlinear behavior.
The system of nonlinear equations resulting from the discretization of the
Navier-Stokes equations can be solved using nonlinear iteration methods, such
as Newton's method. However, fast quadratic convergence is typically only
obtained in a local neighborhood of the solution, and for many configurations,
the classical Newton iteration does not converge at all. In such cases,
so-called globalization techniques may help to improve convergence.
  In this paper, pseudo-transient continuation is employed in order to improve
nonlinear convergence. The classical algorithm is enhanced by a neural network
model that is trained to predict a local pseudo-time step. Generalization of
the novel approach is facilitated by predicting the local pseudo-time step
separately on each element using only local information on a patch of adjacent
elements as input. Numerical results for standard benchmark problems, including
flow through a backward facing step geometry and Couette flow, show the
performance of the machine learning-enhanced globalization approach; as the
software for the simulations, the CFD module of COMSOL Multiphysics is
employed.
Anouk Zandbergen, Tycho van Noorden, Alexander Heinlein
# Improving Pseudo-Time Stepping Convergence for CFD Simulations With Neural Networks

###### Abstract

Computational fluid dynamics (CFD) simulations of viscous fluids described by the Navier-Stokes equations are considered. Depending on the Reynolds number of the flow, the Navier-Stokes equations may exhibit a highly nonlinear behavior. The system of nonlinear equations resulting from the discretization of the Navier-Stokes equations can be solved using nonlinear iteration methods, such as Newton's method. However, fast quadratic convergence is typically only obtained in a local neighborhood of the solution, and for many configurations, the classical Newton iteration does not converge at all. In such cases, so-called globalization techniques may help to improve convergence.

In this paper, pseudo-transient continuation is employed in order to improve nonlinear convergence. The classical algorithm is enhanced by a neural network model that is trained to predict a local pseudo-time step. Generalization of the novel approach is facilitated by predicting the local pseudo-time step separately on each element using only local information on a patch of adjacent elements as input. Numerical results for standard benchmark problems, including flow through a backward facing step geometry and Couette flow, show the performance of the machine learning-enhanced globalization approach; as the software for the simulations, the CFD module of COMSOL Multiphysics(r) is employed.

keywords: Computational fluid dynamics, Newton's method, pseudo-time stepping, neural network, hybrid algorithm Msc: [2020] 35Q30, 65H10, 65N12, 65N22, 68T07 +
Footnote †: journal:

## 1 Introduction

Computational fluid dynamics (CFD) simulations are highly relevant for a wide range of applications, including the fields of aerospace, environmental, and biological engineering as
well as weather predictions and medicine. Numerical simulations of Newtonian fluids, that is, fluid with a linear correlation of the local viscous stresses and strain rate, involves the solution of the Navier-Stokes equations. Depending on the flow regime, the Navier-Stokes equations exhibit strong nonlinearities, which pose challenges to the numerical nonlinear solver employed.

In this work, we consider two-dimensional stationary CFD simulations using a mixed finite element method (FEM) as the discretization. The nonlinear system of equations resulting from discretizing the steady-state Navier-Stokes equations is then solved using Newton's method, that is, via solving a sequence of linearized problems. Whereas Newton's method provides quadratic convergence near the solution, the iteration may only converge slowly or even diverge farther away from the solution. In order to improve global convergence of the method, globalization techniques have been developed; see [1] for an overview of globalization techniques for the Navier-Stokes equations. In [1], the authors categorize globalization techniques into the classes of _backtracking methods_ and _trust region methods_; cf. [2; 3]. Other examples of methods to improve the robustness of the nonlinear solver include homotopy, continuation, pseudo transient continuation, mesh sequencing methods, or nonlinear preconditioning methods; cf. [4; 5; 6; 7; 8; 9; 10; 11; 12; 13].

The combination of scientific computing and machine learning, which is also denoted as scientific machine learning (SciML), is a new field [14], which has gained a lot of attention over the past few years. One popular application area for SciML techniques are CFD simulations; see, for instance, the review paper [15]. Successful examples for the use of SciML techniques in the context of CFD are, for instance, discretization approaches using neural networks (NNs) [16; 17; 18], surrogate models [19; 20; 21; 22], and model discovery [23; 24; 25]. In [26], the temporal evolution of dynamic systems is learned via long short-term memory (LSTM) NNs. Approaches for the enhancement of CFD simulations via machine learning include the generation of optimized meshes [27], the improvement of the resolution of the simulation results [28; 29; 30; 31], or the detection of troubled cells in simulation meshes due to Gibbs oscillations arising from the use of higher-order methods [32]. We note that, in [29; 30; 31; 32], the network inputs and output are specifically chosen from local patches of elements to achieve generalizability of the machine learning model; we will apply a similar approach in this work.

The use of machine learning models for the enhancement of numerical solvers has also been investigated in several works. In the context of linear solvers, [33] and [34] have employed features generated by the proper orthogonal decomposition (POD) to develop a modified conjugate gradient (CG) solver or preconditioning techniques, respectively. The enhancement of the CG method by a nonlinear convolutional NN-based preconditioner has been investigate in [35]. In [36; 37], NNs have been employed to improve domain decomposition methods, and in [38; 39; 40; 41; 42], components of multigrid methods are selected or generated using machine learning models. In the early works [43; 44; 45; 46; 47], suitable numerical solution algorithms are being selected using machine learning techniques. Related to that, [48] investigated the prediction of a suitable combination of preconditioner and iterative method for solving a given sparse linear system.

Fewer works have considered the improvement of the nonlinear solver convergence using
machine learning. In [49], the nonlinear convergence of numerical solvers for nonlinear partial differential equations (PDE) is improved using POD features in a nonlinear preconditioning approach, and the proposed method is tested specifically for CFD problems. Moreover, in [50], dimensionless numbers and simulation properties are used as input to a random forest model to predict a relaxation parameter to accelerate a convergence of the Picard iteration for multiphase porous media flow. To the best of the authors' knowledge, there is currently no work that directly addresses the improvement of the nonlinear convergence of Newton's method for CFD simulations using NNs. The only exception is the master thesis of the first author [51], which laid the foundation for this paper.

In particular, we enhance the pseudo-transient continuation [8], or pseudo-time stepping, approach for improving the robustness of the nonlinear convergence for the steady-state Navier-Stokes equations using an NN model. Similar to previous approaches [29; 30; 31; 32], the model is trained to predict the local pseudo-time step size for each mesh element using only local input features; this aims at ensuring generalizability of our approach. We test the performance and generalizability of our approach for different geometries and flow configurations using CFD simulations with the software package COMSOL Multiphysics(r); as a reference, we compare our approach against two default pseudo-time step control mechanisms implemented in the FEM software package COMSOL Multiphysics(r). As the machine learning framework, we employ Matlab, which enables us to access data from COMSOL via an available interface between the two packages.

The paper is organized as follows. First, the Navier-Stokes equations and pseudo-time stepping are briefly recalled and discussed in section 2 respectively section 3. In section 4 the NN for the local element pseudo-time step prediction is described. The numerical results are given in section 5. Lastly, the conclusions are presented in section 6.

## 2 Navier-Stokes equations

The Navier-Stokes equations are a system of partial differential equations (PDEs) that describes the flow of Newtonian fluids, that is, fluids with a linear correlation of the viscous stresses and the local strain rate. In particular, we consider incompressible fluids, that is, fluids with a constant density; cf. [52]. The stationary Navier-Stokes equations for incompressible flow read

\[\begin{array}{rcl}\rho(\mathbf{u}\cdot\nabla)\mathbf{u}-\mu\nabla^{2}\mathbf{ u}&=&-\nabla p+\mathbf{f},\\ \rho\nabla\cdot\mathbf{u}&=&0,\end{array} \tag{1}\]

describing the conservation of momentum and mass, respectively. In these equations, \(\mathbf{u}\) is the flow velocity, \(\mu\) is the dynamic viscosity of the fluid, \(\rho\) is the density, which is assumed to be constant, \(p\) is the pressure, and \(\mathbf{f}\) is the body force acting on the flow [52]; in all our experiments, the body force is assumed to be zero.

Due to the convective term \(\rho(\mathbf{u}\cdot\nabla)\mathbf{u}\), eq. (1) is nonlinear, and the numerical treatment of this nonlinearity is the subject of this research. Very strong nonlinearities are generally related to turbulent flow. However, already in the regime of laminar flow, the nonlinearities may become severe enough to cause nonlinear solvers to diverge; hence, we will focus on
the laminar case. Flow remains laminar as long as the Reynolds number is below a critical value in the order of \(10^{5}\); see, e.g., [53]. The Reynolds number \(\mathrm{Re}\) is given by

\[\mathrm{Re}=\frac{\rho UL}{\mu}. \tag{2}\]

Here, \(U\) and \(L\) are the typical velocity and length scale. The Reynolds number describes the ratio between inertial and viscous forces; cf. [54]. At low Reynolds numbers, the flow remains laminar as the viscous forces are able to damp out disturbances in the flow. At high Reynolds numbers, the inertial forces become more important resulting in nonlinear interactions to grow, which causes turbulence.

In order to solve the stationary Navier-Stokes equations in eq. (1) on a computational domain \(\Omega\subset\mathbb{R}^{2}\) numerically, we consider a triangulation of \(\Omega\) into linear triangles. Then, we discretize eq. (1) using piecewise linear finite elements for both the velocity \(\mathbf{u}\) and the pressure \(p\) and stabilize using Galerkin least-squares (GLS) streamline diffusion and Hughes-Mallet type shock-capturing [55; 56]. This yields a discrete but nonlinear system of equations

\[N(\mathbf{u}_{h})+B^{\top}p_{h} =0, \tag{3}\] \[B\mathbf{u}_{h} =0,\]

where we denote the discrete velocity and pressure fields as \(\mathbf{u}_{h}\) respectively \(p_{h}\). The linear operators \(B\) and \(B^{\top}\) correspond to the divergence and gradient, respectively, and the nonlinear operator \(N\) corresponds to the convection-diffusion of the velocity field.

By combining the fields \(\mathbf{u}_{h}\) and \(p_{h}\) into a single field \(\mathbf{v}\), we can simply formulate the discrete nonlinear problem eq. (3) as

\[F(\mathbf{v})=0.\]

We solve this system using Newton's method, that is, by solving a sequence of linear systems of the form

\[-F^{\prime}(\mathbf{v}^{n})\Delta\mathbf{v}^{n}=F(\mathbf{v}^{n}),\]

where \(n\) is the iteration index of the Newton iteration and \(F^{\prime}(\mathbf{v}^{n})\) is the Jacobian at the linearization point \(\mathbf{v}^{n}\). Close to the solution, Newton's method converges with quadratic rate. However, for an arbitrary initial guess, the method might also diverge. As we will see in our numerical results, this may happen for a relatively wide range of cases for the stationary Navier-Stokes equations. In order to make the nonlinear iteration more robust, globalization techniques can be employed; see, e.g., [1] for an overview of globalization techniques for Newton's method for the Navier-Stokes equations. Here, we will consider the pseudo-time stepping approach [8], which we will discuss in more detail in the next section.

## 3 Pseudo-time stepping

Solving stationary nonlinear equations using Newton's method requires an initial guess close enough to the root. Globalization techniques are designed to give a result for a wider
range of initial guesses, but can stagnate at local minima [8]. Pseudo-time stepping, or pseudo-transient continuation, is an alternative to Newton's method for computing stationary solutions of time-dependent partial differential equations. The idea is to obtain a solution more robustly by solving the time-dependent equation instead of the stationary one. The idea behind the method is to frame the problem in a time-depending setting,

\[\frac{\partial\mathbf{v}}{\partial t}=-F(\mathbf{v}), \tag{4}\]

with \(F(\mathbf{v})=0\) the system of nonlinear equations of which a solution must be determined. Now the algorithm for pseudo-time stepping can be described as the numerical integration with a variable time step method of the initial value problem

\[\frac{\partial\mathbf{v}}{\partial t}=-F(\mathbf{v}),\quad\mathbf{v}(0)= \mathbf{v}_{0}; \tag{5}\]

cf. [8]. The idea is that this time integration converges to a steady state for a wider range of initial values \(\mathbf{v}_{0}\) than the standard Newton method. In addition, the method tries to increase the time step as \(F(\mathbf{v})\) approaches \(0\), such that when the iterations get closer to a steady state, the convergence becomes (close to) quadratic.

The iterations of the pseudo-time stepping algorithm are given by

\[\mathbf{v}^{n+1}=\mathbf{v}^{n}-\left((\Delta t^{n})^{-1}I+F^{\prime}(\mathbf{ v}^{n})\right)^{-1}F(\mathbf{v}^{n}), \tag{6}\]

with \(F^{\prime}(\mathbf{v}^{n})\) the Jacobian. To understand how (6) is obtained from (5) consider an Euler backward step starting from \(\mathbf{v}^{n}\) for equation (5):

\[\mathbf{z}^{n+1}=\mathbf{v}^{n}-\Delta t^{n}F(\mathbf{z}^{n+1}). \tag{7}\]

Note that \(\mathbf{z}^{n+1}\) is a root of

\[G(\xi):=\xi+\Delta t^{n}F(\xi)-\mathbf{v}^{n}, \tag{8}\]

as a function of \(\xi\). The first Newton iterate for solving \(G(\xi)=0\) with initial guess \(\xi_{0}=\mathbf{v}^{n}\) gives

\[\xi_{1} =\mathbf{v}^{n}-\left(I+\Delta t^{n}F^{\prime}(\mathbf{v}^{n}) \right)^{-1}(\mathbf{v}^{n}+\Delta t^{n}F(\mathbf{v}^{n})-\mathbf{v}^{n}),\] \[=\mathbf{v}^{n}-\left((\Delta t^{n})^{-1}+F^{\prime}(\mathbf{v}^{ n})\right)^{-1}F(\mathbf{v}^{n}).\]

From this corrector iteration, we have obtained the formula in (6). The pseudo code of the method is the given by Algorithm 1.

### Applying pseudo-time stepping to the weak form of the Navier-Stokes equations

The weak form of the Navier-Stokes equations is given by: find velocity and pressure fields \((\mathbf{u},p)\in X\) such that

\[\begin{array}{rcl}\int_{\Omega}(\rho\partial_{t}\mathbf{u}+\rho( \mathbf{u}\cdot\nabla)\mathbf{u})\cdot\phi-(\mu\nabla\mathbf{u}\cdot\nabla \phi-p\nabla\cdot\phi)&=&\int_{\Omega}\mathbf{f}\cdot\phi+\int_{ \Gamma_{N}}\mathbf{h}\cdot\phi,\\ \int_{\Omega}\rho\nabla\cdot\mathbf{u}\psi&=&0,\end{array} \tag{9}\]
for all test functions \((\phi,\psi)\in Y\), with \(X\) and \(Y\) suitable function spaces, e.g., \(X=Y=H_{0}^{1}(\Omega)^{2}\times L^{2}(\Omega)\). Here \(\mathbf{f}\) is a forcing term, and \(\mathbf{h}\) encodes the boundary conditions.

To apply pseudo-time stepping to this weak form, we first discretize in space using finite elements by replacing the function spaces \(X\) and \(Y\) by finite dimensional subspaces \(X_{h}\subset X\) and \(Y_{h}\subset Y\), consisting of continuous functions that are piecewise linear on a finite element mesh of the domain \(\Omega\) with maximum mesh size \(h\). Then, in order to apply pseudo-time stepping to the stationary form of the spatially discretized system, the same conceptual idea as described in the previous section is used, which can be summarized by first writing the time-discretized equations for one Euler backward step, and then linearizing this system for the new time step around the current time step. The Euler backward formula with time step \(\Delta t^{n}\) for the discretized weak form of (9) gives the following (discretized) nonlinear problem: find \((\mathbf{u}_{h},p_{h})\in X_{h}\) such that

\[\sum_{e\in E}\Biggl{(}\int_{\Omega^{e}}(\rho\frac{\mathbf{u}_{h} ^{n+1}-\mathbf{u}_{h}^{n}}{\Delta t^{n}}+\rho(\mathbf{u}_{h}^{n+1}\cdot\nabla) \mathbf{u}_{h}^{n+1})\cdot\phi_{j}-(\mu\nabla\mathbf{u}_{h}^{n+1}\cdot\nabla \phi_{j}-p_{h}^{n+1}\nabla\cdot\phi_{j})+\] \[-\!\!\int_{\Omega^{e}}\!\!\mathbf{f}\cdot\phi_{j}-\int_{\Gamma_{ N}^{e}}\!\!\mathbf{h}\cdot\phi_{j}\Biggr{)} =0, \tag{10}\] \[\sum_{e\in E}\Biggl{(}\int_{\Omega^{e}}\!\!\rho\nabla\cdot\mathbf{ u}_{h}^{n+1}\psi_{j}\Biggr{)} =0,\]

for all test functions \((\phi_{j},\psi_{j})\in Y_{h}\), where the summation is over the set of all mesh elements \(E=\{e_{1},e_{2},\ldots,e_{m}\}\), with the number of elements \(m\).

Linearizing this system around \((\mathbf{u}_{h}^{n},p_{h}^{n})\) results in

\[\sum_{e\in E}\Biggl{(}\int_{\Omega^{e}}(\rho\frac{\mathbf{u}_{h }^{n+1}-\mathbf{u}_{h}^{n}}{\Delta t^{n}}+\rho(\mathbf{u}_{h}^{n}\cdot\nabla) \mathbf{u}_{h}^{n+1}+\rho(\mathbf{u}_{h}^{n+1}\cdot\nabla)\mathbf{u}_{h}^{n}- \rho(\mathbf{u}_{h}^{n}\cdot\nabla)\mathbf{u}_{h}^{n})\phi_{j}+\] \[-\!\!\int_{\Omega^{e}}(\mu\nabla\mathbf{u}_{h}^{n+1}\cdot\nabla \phi_{j}-p_{h}^{n+1}\nabla\cdot\phi_{j})-\int_{\Omega^{e}}\!\!\mathbf{f}\cdot \phi_{j}-\int_{\Gamma_{N}^{e}}\!\mathbf{h}\cdot\phi_{j}\Biggr{)} =0, \tag{11}\] \[\sum_{e\in E}\Biggl{(}\int_{\Omega^{e}}\!\!\rho\nabla\cdot\mathbf{ u}_{h}^{n+1}\psi_{j}\Biggr{)} =0.\]
Solving this linear system for \((\mathbf{u}_{h}^{n+1},p_{h}^{n+1})\) constitutes one pseudo-time step. We write eq. (10) in matrix format

\[\frac{1}{\Delta t^{n}}M(\mathbf{v}^{n+1}-\mathbf{v}^{n})+F(\mathbf{v}^{n+1})=0, \tag{12}\]

where the coefficients of \((\mathbf{u}_{h}^{n},p_{h}^{n})\) w.r.t. the chosen basis of \(X_{h}\) are again combined into a single vector \(\mathbf{v}^{n}\), and where \(M=\begin{bmatrix}M_{\mathbf{u}}&\mathbf{0}\\ \mathbf{0}&\mathbf{0}\end{bmatrix}\) is a singular mass matrix with the velocity mass matrix \(M_{\mathbf{u}}\). Then, we can also write the solution of eq. (11) in the more compact form

\[\mathbf{v}^{n+1}=\mathbf{v}^{n}-(\frac{1}{\Delta t^{n}}M+F^{\prime}(\mathbf{v} ^{n}))^{-1}F(\mathbf{v}^{n}). \tag{13}\]

### Choice of the pseudo-time step

The idea is that the pseudo-time step \(\Delta t^{n}\) is initially small, such that its influence in (6) is large. When the iteration approaches convergence, the magnitude of the pseudo-time steps becomes small, such that (6) performs more like the standard Newton iteration. This has the benefit that, in the beginning, the method depends less on the initial guess, and in the end, it may converge quadratically. There are several possibilities to define a pseudo-time step that fits this property.

In addition to that, it is possible to use a pseudo-time step that uses local, i.e. mesh element dependent, information. This amounts to replacing the global pseudo-time step \(\Delta t^{n}\) in (11) by a local, mesh element dependent, pseudo-time step \(\Delta t^{n}_{e}\). In this case the mass matrix \(M\) in (13) will depend on the vector of local pseudo-time steps \(\overline{\Delta t}^{n}=(\Delta t^{n}_{e_{1}},\ldots,\Delta t^{n}_{e_{m}})\), so that (13) becomes

\[\mathbf{v}^{n+1}=\mathbf{v}^{n}-(M(\overline{\Delta t}^{n})+F^{\prime}( \mathbf{v}^{n}))^{-1}F(\mathbf{v}^{n}). \tag{14}\]

In COMSOL, for instance, the following local pseudo-time step is used:

\[\Delta t^{n}_{e}=\text{CFL}(n)\frac{h_{e}}{\|\mathbf{u}_{h}^{n}\|_{e}}, \tag{15}\]

with \(h_{e}\) the size of mesh element \(e\), \(\|\mathbf{u}_{h}^{n}\|_{e}\) the Euclidean norm of the velocity field \(\mathbf{u}_{h}^{n}\) in the mesh element \(e\), and \(\text{CFL}(n)\) a global Courant-Friedrichs-Lewy (CFL) number [54]. This global CFL number may depend on the iteration count \(n\), for instance,

\[\text{CFL}_{iter}(n)=\begin{cases}1.3^{\min(n,9)}&1\leq n\leq 20,\\ 1.3^{9}+9\cdot 1.3^{\min(n-20,9)}&20<n\leq 40,\\ 1.3^{9}+9\cdot 1.3^{9}+90\cdot 1.3^{\min(n-40,9)}&n>40;\end{cases} \tag{16}\]

see [54, p. 1108]. Alternatively, it can be given by a controller based on the nonlinear error estimate \(e_{n}\) for iteration \(n\), the given target error estimate "tol", and control parameters \(k_{P}\), \(k_{I}\), and \(k_{D}\), which are positive constants:

\[\text{CFL}_{e}(n)=\left(\frac{e_{n-2}}{e_{n-1}}\right)^{kp}\left(\frac{\text{ tol}}{e_{n-1}}\right)^{k_{I}}\left(\frac{e_{n-2}/e_{n-1}}{e_{n-3}/e_{n-2}} \right)^{k_{D}}\text{CFL}_{e}(n-1), \tag{17}\]
these options are available for controlling the global CFL number in COMSOL.

For both cases, in the initial iterations of pseudo-time stepping, the method starts with a small global CFL number and gradually increases it as the solution reaches convergence [54]. In eq. (16) this is because the CFL number increases each iteration, also shown in Figure 1. For eq. (17), when the solution is near convergence, the error estimates \(e_{n}\) will be smaller and thus the CFL number increases.

In this paper we investigate replacing these two pseudo-time step control algorithms by the use of an NN to control the local, mesh element dependent, pseudo-time step \(\Delta t_{e}^{n}\). In order to make the method easily generalizable to different meshes, the NN is provided with local information from the corresponding mesh element and its adjacent mesh elements. One could think of local information about the solution, residuals, mesh and the cell Reynolds number. For example, elements with high local residuals could then have smaller local pseudo-time steps compared to elements with low residuals, which could accelerate convergence in areas which already have low residuals. So, treating each mesh element individually using local information can accelerate convergence for the whole simulation.

## 4 Neural network for the local pseudo-time step prediction

In order to predict an element-wise local pseudo-time step for the pseudo-time stepping that yields fast and robust Newton convergence, we employ a data-based approach based on machine learning techniques; more specifically, we use artificial neural networks (ANNs). The most basic form of an ANN is a multilayer perceptron (MLP) and is given as a composed function

\[\mathcal{NN}\left(x\right)=W_{n+1}\,\,f_{n}\circ\cdots\circ f_{1}\left(x \right),\]

with

\[f_{i}\left(x\right)=\sigma\left(W_{i}x+b_{i}\right). \tag{18}\]

Here, the \(W_{i}\in\mathbb{R}^{n_{i}\times n_{i-1}}\) and \(b_{i}\in\mathbb{R}^{n_{i}}\) are the so-called weight matrices and bias vectors, which represent the linear parts of the NN function \(\mathcal{NN}\); see Figure 6 for an exemplary net

Figure 1: A plot of CFL\({}_{iter}\) depending on the nonlinear iteration count, as defined in eq. (16).
work architecture of an MLP. Given a specific network architecture, that is, the dimensions \(\left(n_{i}\right)_{i}\), the MLP is determined by the coefficients of the weight matrices and bias vectors, also denoted as the network parameters. The MLP becomes nonlinear due to composition with the activation function \(\sigma\). Here, we will employ the rectified linear unit (ReLU) function \(\sigma\left(x\right)=\max\left(0,x\right)\). MLPs, and NNs in general, are well-suited for approximating nonlinear functions; see, e.g., [57]. In order fit an NN to input data \(X=\left\{x_{i}\right\}\) and corresponding output data \(Y=\left\{y_{i}\right\}\), a loss function \(\mathcal{L}\) is minimized with respect to the network parameters:

\[\arg\min_{W_{i},b_{i}}\mathcal{L}\left(\mathcal{NN}\left(X\right),Y\right)\]

The loss function penalizes deviation of the network function from the reference output data \(Y\). In order to optimize the NN parameters (a.k.a. network training), we employ a stochastic gradient descent (SGD) method using adaptive moment estimation (Adam) [58]; the gradients are computed using the backpropagation algorithm [59]. For a more detailed introduction to deep learning and NNs; see, e.g., [60].

In this section, we will discuss how the local pseudo-time steps used as reference data in the root mean squared error (RMSE) loss function are computed as well as the composition of our training data set and the network architecture.

### Loss and optimal local pseudo-time step

Our goal is to predict local pseudo-time steps that enhance the pseudo-time stepping method in accelerating convergence of a given flow problem. In particular, we try to predict the vector of local pseudo-time steps \(\overline{\Delta t}_{opt}\) such that the next iterate is as close as possible to the solution of the nonlinear problem. We will denote this as the _the optimal local pseudo-time step_ for an iterate \(\mathbf{v}^{n}\), \(\overline{\Delta t}_{opt}(\mathbf{v}^{n})\), defined by

\[\overline{\Delta t}_{opt}(\mathbf{v}^{n}):=\arg\min_{\overline{\Delta t}}|| \mathbf{v}^{n}-\left(M(\overline{\Delta t})+F^{\prime}(\mathbf{v}^{n})\right)^ {-1}F(\mathbf{v}^{n})-\mathbf{v}^{*}||, \tag{19}\]

where \(\mathbf{v}^{*}\) satisfies \(F(\mathbf{v}^{*})=0\). The norm that is used here is the \(L^{2}\) norm of the reconstructed velocity components.

Figure 2: The local pseudo-time step optimized by SNOPT; the difference between the obtained solution and the converged solution was minimized up to a value of 3.1\(\cdot\)10\({}^{-5}\).
Since the optimal local pseudo-time step cannot be computed explicitly, we compute an approximation via an optimization procedure, using the sparse nonlinear optimizer (SNOPT) software package [61; 62]. This optimization is performed within COMSOL; see [51] for more details. For the example of a back-step geometry with laminar flow, we obtain the optimal local pseudo-time step as shown in Figure 2 by using this two step procedure.

Once we have computed the local optimal pseudo-time step for all elements in the data set, we optimize the squared error of the network prediction against these reference values. The whole data set consists of a large number of optimal local pseudo-time steps \(\Delta t_{e_{i},opt}(\mathbf{v}^{n})\) computed from various problem configurations, nonlinear iterates \(\mathbf{v}^{n}\), and elements \(e_{i}\) computed using the optimization in COMSOL. To define the complete loss, let \(\Delta t^{(i)}_{pred}\) and \(\Delta t^{(i)}_{opt}\) be the network prediction and optimized pseudo-time steps, respectively, for the \(i\)th data point in a data set with \(n\) data points. The full RMSE loss, is then given as:

\[\text{loss}=\sqrt{\frac{1}{n}\sum_{i=1}^{n}\Big{(}\Delta t^{(i)}_{pred}-\Delta t ^{(i)}_{opt}\Big{)}^{2}}. \tag{20}\]

### Neural network input

NNs are universal function approximators and can theoretically approximate any continuous nonlinear function up to arbitrary precision; see, for instance, [57; 63; 64]. Therefore, we expect that it is possible to predict the optimal local pseudo-time step, given sufficient input data. On the other hand, we want to limit the complexity of the input data and train a model with strong generalization properties. Therefore, similar to other approaches for NN-enhanced simulations [29; 30; 32], we employ local input data to predict the local optimal pseudo-time step. This means that we do not include any information about the specific global boundary value problem but only data that is given on a local patch of elements around the element of interest, that is, where the local optimal pseudo-time step is to be predicted. As a result, the trained NN model is applicable to any considered laminar flow problem. An additional benefit of using only local data is that the dimensionality and complexity of the data is being limited. Note that, as also noted in [29], due to the local

Figure 3: Data point which contains information of a patch of four elements \(e_{j,1}\) to \(e_{j,4}\) with their vertices \(v_{j,1}\) to \(v_{j,6}\); see Table 1 for all input features for the NN model on each element patch.
data approach, a single simulation already generates a rich data set: each simulation yields data on a large set of elements over a number of Newton iterations. However, to generate a training data set with sufficient variation to allow for generalization of the model, the training data may be sampled from different simulations with varying flow conditions, induced by different boundary conditions and geometries as well as mesh refinement levels.

Let us now discuss the specific data included in the model input. As mentioned before, we employ data from a patch of elements

\[P_{j}=\bigcup_{i=1}^{4}e_{j,i}.\]

In particular, the patch consists of the element for which we want to compute the local optimal pseudo-time step, \(e_{j,1}\), as well as the (up to) three adjacent elements \(e_{j,2}\), \(e_{j,3}\), and \(e_{j,4}\); cf. figure 3. Specifically, we denote two elements as adjacent if they are connected via an element edge. On each patch, we sample the data listed in table 1 with their respective location. The nodal information is explicitly available from the mixed finite element discretization used in our simulations; cf. the description of our simulation setup in section 2.

Let us briefly motivate our choices. Firstly, we include the element size, which is implied by the element edge lengths, and \(\mathbf{u}\) since they are also employed in the computation of the pseudo-time step approach in eq. (15). We complement this \(\mathbf{u}\) by \(p\) to provide a complete description of the finite element solution on the local patch. Moreover, the local information describing the convergence of the Newton iteration is supplied in the form of two different

\begin{table}
\begin{tabular}{|l|l|l|} \hline
**variable** & **unit** & **location** \\ \hline \hline element edge length & m & - \\ \hline \(u\) & \(\nicefrac{{\mathrm{m}}}{{\mathrm{s}}}\) & \(v_{j,1}\), \(v_{j,2}\), \(v_{j,3}\) \\ \hline \(v\) & \(\nicefrac{{\mathrm{m}}}{{\mathrm{s}}}\) & \(v_{j,1}\), \(v_{j,2}\), \(v_{j,3}\) \\ \hline \(p\) & Pa & \(v_{j,1}\), \(v_{j,2}\), \(v_{j,3}\) \\ \hline \(R_{u}\) & & \(v_{j,1}\), \(v_{j,2}\), \(v_{j,3}\) \\ \hline \(R_{v}\) & & \(v_{j,1}\), \(v_{j,2}\), \(v_{j,3}\) \\ \hline \(R_{p}\) & & \(v_{j,1}\), \(v_{j,2}\), \(v_{j,3}\) \\ \hline \(r_{u}\) & \(\nicefrac{{\mathrm{N}}}{{\mathrm{m}^{3}}}\) & \(v_{j,1}\), \(v_{j,2}\), \(v_{j,3}\) \\ \hline \(r_{v}\) & \(\nicefrac{{\mathrm{N}}}{{\mathrm{m}^{3}}}\) & \(v_{j,1}\), \(v_{j,2}\), \(v_{j,3}\) \\ \hline \(r_{p}\) & \(\nicefrac{{\mathrm{kg}}}{{\mathrm{(m^{3}\cdot s)}}}\) & \(v_{j,1}\), \(v_{j,2}\), \(v_{j,3}\) \\ \hline cell Reynolds number & & \(e_{j,1}\) \\ \hline \end{tabular}
\end{table}
Table 1: Variables used for input for the NN with their given unit and location on one element in the patch, in this case the central element. One data point contains information of four elements, with a non-structured order of the adjacent elements. The two different types of residuals \(r_{\alpha}\) and \(R_{\alpha}\), \(\alpha=u\), \(v\), \(p\) are given in eqs. (21) and (22).
types of residuals, i.e., the residual of the discretized system of nonlinear equations

\[F(\mathbf{v})=\begin{pmatrix}R_{u}\\ R_{v}\\ R_{p}\end{pmatrix}=\begin{pmatrix}N(u,v)+B^{\top}p\\ B(u,v)\end{pmatrix}=\begin{pmatrix}N(\mathbf{u})+B^{\top}p\\ B(\mathbf{u})\end{pmatrix} \tag{21}\]

and the residual obtained by substituting the approximate solution into the system of PDEs eq. (1):

\[r_{u} =\rho(u\partial_{x}u+v\partial_{y}u)-\mu(\partial_{x}^{2}u+ \partial_{y}^{2}u)+\partial_{x}p+f_{x},\] \[r_{v} =\rho(u\partial_{x}v+v\partial_{y}v)-\mu(\partial_{x}^{2}v+ \partial_{y}^{2}v)+\partial_{y}p+f_{y}, \tag{22}\] \[r_{p} =\rho(\partial_{x}u+\partial_{y}v).\]

The local Reynolds number, as given in equation eq. (2), provides additional information about the local nonlinearity of the flow; it is also used in papers [29; 30; 50] as input for the machine learning model. The input is normalized to have mean 0 and standard deviation by computing the z-score defined as

\[z(\sigma)=\frac{x-\mu}{\sigma},\]

with \(x\) the evaluated input data, \(\mu\) the mean of the input, and \(\sigma\) the standard deviation of the input [65]. The the mean and standard deviation for centering and scaling the data are saved and later on used to normalize input data from other simulations.

Note that, for practical reasons, in our implementation in COMSOL, we collect all input information, as listed in table 1, element by element from each patch. This means that vertex- and edge-based information is duplicated where the elements in the patch touch. As a result, we obtain a total of \(n_{0}=124\) input features for our network model. In practice, one may want to omit the redundant data, however, this is currently not straightforward based on the data structures in our interface from Matlab to COMSOL since the elements in each patch cannot be easily retrieved in a consistent ordering. Moreover, the redundant information implicitly encodes connectivity information, which could alternatively be encoded by a consistent ordering of the elements within the patch. Our numerical results in section 5 indicate that this handling of the input data does not prevent our model to learn from the data. However, in future work, we plan to investigate if using unique input data with a consistent ordering may have a positive effect on our model.

Finally, elements which are adjacent to the boundary of the computational domain \(\Omega\), may only have one or two neighboring elements, resulting in patches consisting of two or three elements only. In order to encode this in the data, we set any data of the vertices and edges outside the domain to zero; other approaches for encoding this information are possible and will, again, be considered in future work.

### Generation of training data

We train our model based on simulation data from several boundary value problems defined on back-step geometries; see Figure 2 for an exemplary back-step geometry. The
boundary conditions for the back-step flow are defined as follows: We impose a parabolic inflow velocity profile at the inlet \(\Gamma_{in}\) and a constant pressure at the outlet \(\Gamma_{out}\) on the right side of the geometry. For the other parts of the boundary, \(\Gamma_{wall}\), a no-slip boundary condition is prescribed.

As discussed in section 4.2, due to the locality of the input features, we expect that it suffices to use training data for a single type of geometry, as long as there is enough variation in the data with respect to the individual patches. In particular, we will consider two different base cases of back-step geometries, which we denote as B1 and B2. For both cases, we consider meshes with different levels of refinement and inflow velocities. This also results in flow fields with different Reynolds numbers. Moreover, we use data from different Newton iterations since, as can be seen in Figure 8, the iterate changes quite drastically between different Newton iterations. Additionally, we consider scaled variants of B1 and B2, which we denote as B1S and B2S, respectively. These are obtained by scaling the geometric dimensions of B1 and B2 by a factor of 0.1 and increasing the inflow velocity by a factor of 10; cf. Table 2 for the dimensions of the four different geometries. As a result the Reynolds number remains the same, and we would expect the nonlinear convergence to be similar. By adding the scaled versions to the training configurations, we try to force the model to learn the dependence of the nonlinearity on the Reynolds number. We summarize all configurations considered for the generation of training data in Table 3.

In total, we obtain 79 256 data points from all configurations listed in Table 3. To arrive at a balanced distribution of the training data set, we perform the sampling such that we obtain roughly the same number of data points for each element size in the data set; that is,

\begin{table}
\begin{tabular}{|l|c|c|c|c|} \hline
**back-step geometry** & **dim. inflow tunnel (m)** & **dim. outflow tunnel (m)** & **range inflow velocity (m/s)** & **range maximum element size (m)** \\ \hline \hline
**B1** & 0.05 \(\times\) 0.25 & 0.12 \(\times\) 1.15 & 0.001 – 0.015 & 0.0106 – 0.0256 \\ \hline
**B1S** & 0.005 \(\times\) 0.025 & 0.012 \(\times\) 0.115 & 0.01 – 0.15 & 0.00106 – 0.00256 \\ \hline
**B2** & 0.08 \(\times\) 0.25 & 0.22 \(\times\) 1.15 & 0.001 – 0.01 & 0.0126 – 0.0206 \\ \hline
**B2S** & 0.008 \(\times\) 0.025 & 0.022 \(\times\) 0.115 & 0.01 – 0.1 & 0.00126 – 0.00206 \\ \hline \end{tabular}
\end{table}
Table 2: Dimensions of the back-step geometries considered for the training and test data; see Figure 4 for the base geometry.

Figure 4: Back-step geometries employed: B1, B1S, B2, and and B2S vary in the width of the inlet \(\Gamma_{in}\) (red), where a parabolic velocity profile is prescribed, and the length of the blue wall segment; cf. Table 2. At the outlet \(\Gamma_{out}\), we prescribe a constant pressure boundary condition, and at all other parts of the boundary, which are denoted as \(\Gamma_{wall}\), a no-slip boundary condition is enforced. For an exemplary flow field, see Figure 8.
we select exactly \(3\,500\) data points from all configurations with the same maximum element size. In total, we can categorize \(14\) different simulations based on the maximum element size, which gives a data set of \(49\,000\) data points. To make it easier to applied batch learning for the NN, \(1\,000\) data points are discarded at random. So, in total there are \(48\,000\) data points, of which \(70\%\) are used for training our network model, \(15\%\) are used for testing and \(15\%\) are used for validation. In section 5, we will first discuss how this model performs on configurations with back-step geometry in section 5.1. Then, in sections 5.2 and 5.3, we

\begin{table}
\begin{tabular}{|l|r|r|r|r|} \hline
**back-step** & **number of** & **max. element** & **inflow velocity** & **\# iterations** \\ \cline{2-5}
**geometry** & **elements** & **size (m)** & **(m/s)** & **\# iterations** \\ \hline \hline \multirow{6}{*}{**B1**} & \multirow{6}{*}{\(3\,908\)} & \multirow{6}{*}{\(0.0206\)} & \(0.001\) & \(1\), \(10\) \\  & & & \(0.005\) & \(10\) \\  & & & \(0.01\) & \(10\) \\ \cline{2-5}  & & & \(0.003\) & \(1\), \(10\) \\ \cline{2-5}  & \multirow{3}{*}{\(3\,908\)} & \(0.0206\) & \(0.006\) & \(10\) \\  & & & \(0.009\) & \(10\) \\ \cline{2-5}  & & & \(0.002\) & \(1\), \(10\) \\  & & & \(0.007\) & \(10\) \\  & & & \(0.008\) & \(2\), \(4\) \\  & & & \(0.015\) & \(10\) \\ \cline{2-5}  & \multirow{3}{*}{\(13\,828\)} & \multirow{3}{*}{\(0.0106\)} & \(0.004\) & \(2\), \(10\) \\  & & & \(0.008\) & \(2\), \(10\) \\ \hline \multirow{6}{*}{**B1S** (scaled B1)} & \multirow{6}{*}{\(3\,908\)} & \multirow{6}{*}{\(0.00206\)} & \(0.01\) & \(1\), \(10\) \\  & & & \(0.05\) & \(10\) \\ \cline{2-5}  & & & \(0.1\) & \(10\) \\ \cline{2-5}  & \multirow{3}{*}{\(3\,558\)} & \multirow{3}{*}{\(0.00266\)} & \(0.03\) & \(1\), \(10\) \\  & & & \(0.06\) & \(10\) \\ \cline{2-5}  & & & \(0.09\) & \(10\) \\ \cline{2-5}  & \multirow{3}{*}{\(13\,828\)} & \multirow{3}{*}{\(0.00106\)} & \(0.02\) & \(1\), \(10\) \\  & & & \(0.07\) & \(10\) \\ \cline{2-5}  & & & \(0.08\) & \(2\), \(4\) \\ \cline{2-5}  & & & \(0.15\) & \(10\) \\ \cline{2-5}  & \multirow{3}{*}{\(13\,828\)} & \multirow{3}{*}{\(0.00106\)} & \(0.04\) & \(2\), \(10\) \\ \cline{2-5}  & & & \(0.00106\) & \(0.08\) & \(2\), \(10\) \\ \cline{2-5}  & & \(5\,808\) & \(0.0156\) & \(0.005\) & \(1\), \(10\) \\ \cline{2-5}  & \multirow{3}{*}{\(4\,134\)} & \(0.0186\) & \(0.05\) & \(1\), \(10\) \\ \cline{2-5}  & & \(4\,134\) & \(0.00186\) & \(0.1\) & \(2\), \(10\) \\ \cline{2-5}  & & \(4\,134\) & \(0.0186\) & \(0.01\) & \(2\), \(10\) \\ \cline{2-5}  & \multirow{3}{*}{\(8\,790\)} & \(0.00126\) & \(0.013\) & \(3\), \(10\) \\ \cline{2-5}  & & \(5\,808\) & \(0.00156\) & \(0.05\) & \(1\), \(10\) \\ \cline{2-5}  & & \(4\,134\) & \(0.00186\) & \(0.1\) & \(2\), \(10\) \\ \cline{2-5}  & \multirow{3}{*}{\(8\,790\)} & \(0.00126\) & \(0.13\) & \(3\), \(10\) \\ \cline{2-5}  & & & \(0.00126\) & \(0.13\) & \(3\), \(10\) \\ \hline \end{tabular}
\end{table}
Table 3: Different combinations of mesh sizes and flow velocities of the back-step simulations to create data. The number of obtained data points is given in column 2, which are sampled from the nonlinear iteration steps given in the last column. The dimensions of the four back-steps can be found in Table 2.
will also discuss the generalization to other configurations, that is, to Couette flow and flow around an obstacle, respectively.

### Network architecture and training

The architecture of the network model employed in this work is depicted in Figure 6: the model consists of an input layer with 124 neurons, two hidden layers with 16 neurons each, and an output layer with the neuron representing prediction of the local optimal pseudo-time step; the composition of the input data has been discussed in section 4.2.

We determined the network architecture as follows: We first performed a grid search with 6-fold cross validation on the search space given in Table 4 to determine three different network architectures with an acceptable average validation loss. Note that the loss is defined based on the prediction of the local optimal pseudo-time step; that means a low validation loss does not necessarily result in a low solution residual when the network is used. From the grid search, we determined one network with small, one with medium, and one with high networks capacity, that is, networks with:

* 2 hidden layers and 16 neurons per layer,
* 3 hidden layers and 64 neurons per layer,

\begin{table}
\begin{tabular}{|l|l|} \hline
**hyper parameter** & **search space** \\ \hline \hline \# hidden layers & \(\{2,3,4,5\}\) \\ \# neurons per hidden layer & \(\{2^{4},2^{5},2^{6},2^{7},2^{8}\}\) \\ \hline \end{tabular}
\end{table}
Table 4: Hyper parameters and search space for the grid search.

\begin{table}
\begin{tabular}{|l|l|l|} \hline
**geometry** & **inflow/wall** & **max. element size (m)** \\
**type** & **velocity (\(\nicefrac{{m}}{{s}}\))** & \\ \hline \hline
**B1** & \(0.001,0.004,0.007\), & \(0.0106:0.005:0.0256\) \\  & \(0.01,0.012,0.015\) & \\ \hline
**B1S** & \(0.01,0.04,0.07\), & \(0.00106:0.0005:0.00256\) \\  & \(0.1,0.12,0.15\) & \\ \hline
**B2** & \(0.001:0.003:0.01\) & \(0.0126,0.0156,0.0186\), \\  & & \(0.0206\) \\ \hline
**B2S** & \(0.01:0.03:0.1\) & \(0.00126,0.00156,0.00186\), \\  & & \(0.00206\) \\ \hline
**C** & \(0.01,0.03,0.04\), & \(0.014:0.002:0.022\) \\  & \(0.05,0.07,0.1\) & \\ \hline
**CS** & \(0.01,0.03,0.05\) & \(0.0028:0.0004:0.0044\) \\ \hline \end{tabular}
\end{table}
Table 5: Combinations of velocities and mesh sizes for simulations used to test the NN performance, with the back-step configurations B1–B2S in Table 2 and the Couette flow configurations C1 and C2 in Table 6.
* 4 hidden layers and 256 neurons per layer.

Then, in order to find the best model for accelerating the Newton convergence, we compare them when applied to several back-step and Couette simulations, further explained in section 5.1 and section 5.2. The combinations of the inflow or wall velocities with the maximum element size used for these simulations are given in Table 5.

For each network, the required number of nonlinear iterations per simulation was compared. Surprisingly, the largest network now performed much worse, while the smallest network performed best, which could be explained by an overfitting behavior.

Each network is trained based on the different back-step simulations given in Table 2. As mentioned before, the local optimal pseudo-time steps obtained using the SNOPT optimizer are used as the reference. In particular, the network output is compared to these optimized local pseudo-time steps via the RMSE loss function eq. (20). As the optimizer, we use

Figure 5: Training plot of network trained on the optimized pseudo-time step targets.

Figure 6: Exemplary NN with 124 inputs, 2 hidden layers with 16 neurons and an output layer. The arrows are associated with weights in the weight matrices \(W_{i}\); cf. eq. (18).
stochastic gradient descent with adaptive moments (Adam) [58] with an initial learning rate of \(0.001\). As the stopping criterion and regularization, we employ early stopping with a patience of \(150\) epochs.

## 5 Numerical results

In this section, we discuss numerical results for our NN-enhanced pseudo-time stepping approach, which we have introduced in the previous section 4. As discussed in section 4.3, we have trained the model only for back-step geometries. Our model is designed such that only local features are used as input in order to enable generalizability. First, we test our model on the _back-step flow_ training cases, with varying mesh sizes and inflow velocities within the range of the training data; cf. section 5.1. In addition to that, we investigate the performance on flipped and rotated cases. Then, in order to further test the generalization properties of the model, we test our model on _Couette cylinder flows_ in section 5.2 as well as for _flow around an obstacle_ cases in section 5.3. For the Couette cases, we vary the geometry, the mesh size, and the velocity of the inner rotating wall, and for the flow around an obstacle, back-step and Couette cylinder geometries with different kinds of obstacles are

\begin{table}
\begin{tabular}{|l|r|r|r|r|} \hline
**Couette** & **radius outer** & **radius inner** & **range wall** & **range maximum** \\
**geometry** & **circle (m)** & **circle (m)** & **velocity (m/s)** & **element size (m)** \\ \hline \hline
**C** & 0.4 & 0.2 & 0.01 – 0.1 & 0.014 – 0.022 \\ \hline
**CS** & 0.08 & 0.04 & 0.01 – 0.05 & 0.0028 – 0.0044 \\ \hline \end{tabular}
\end{table}
Table 6: Geometrical dimensions of the two different Couette cylinders with corresponding ranges of velocities for the boundary conditions and maximum element sizes.

Figure 7: Convergence plot of a B1 simulation for all three strategies for the local pseudo-time step under consideration: using the NN model as well as the strategies defined in eqs. (16) and (17).
considered.

In particular, we compare the convergence of our hybrid approach with the convergence of two classical strategies for the choice of the local pseudo-time step described in section 3.2, that is, the strategy denoted by \(\text{CFL}_{iter}\), as given by eqs. (15) and (16), and the strategy denoted by \(\text{CFL}_{e}\), as given by eqs. (15) and (17). Furthermore, in the appendix, the convergence of the network is also compared with two variants of Newton's method, that is, the standard Newton method with a constant damping factor and a variant with an adaptive choice of the damping factor; see appendix A for more details. The \(L^{2}\) norm of the "reconstructed" residual \(\tilde{R}\)

\[\|\tilde{R}\|_{L^{2}(\Omega)}=\sqrt{\int_{\Omega}\tilde{R}_{u}^{2}+\tilde{R}_{ v}^{2}+\tilde{R}_{p}^{2}\,dx} \tag{23}\]

is used as the stopping criterion. Here, \((\tilde{R}_{u},\tilde{R}_{v})\) and \(\tilde{R}_{p}\) are the residual velocity and pressure fields in \(X_{h}\) reconstructed by using the components of the residual vector \(R=(R_{u},R_{v},R_{p})\), as defined in eq. (21), as coefficients w.r.t. the chosen basis of \(X_{h}\).

Note that COMSOL uses a scaled version of the residual eq. (23) in the stopping criterion, such that the tolerance may differ depending on the problem configuration. However, when reporting the convergence results, we make sure that the same tolerance for the unscaled residual eq. (23) is used to determine convergence for the different approaches.

All simulations have been performed using COMSOL Multiphysics(r) version 6.1.

### Back-step geometries

In this subsection, we discuss the performance of our model on the same four back-step geometries which we also used for the training; cf. section 4.3 and table 2. In particular, we

Figure 8: Iterations of B1 with inflow velocity \(0.001\,\text{m/s}\) and maximum element size \(0.0156\,\text{m}\) using NN, with velocities (left) and corresponding distribution of the predicted local pseudo-time step prediction (right).
investigate if the network can accelerate the convergence when varying inflow velocity and mesh size on those configurations; cf. Table 5. As a result of those variations, a total of 80

Figure 9: The results of the back-steps for each method; the simulations are ordered from smallest required number of nonlinear iterations to largest separately for each approach. The horizontal lines indicate the average nonlinear iteration counts for the different approaches.
different simulations for the back-step geometries are analyzed.

Exemplary, we plot the convergence history for the back-step geometry B1 with inflow velocity \(0.001\,\nicefrac{{m}}{{s}}\) and a maximum element size of \(0.0156\,\mathrm{m}\) in Figure 7; in particular, we compare the local optimal pseudo-time step predicted by the network model, which we denote as NN, against the \(\mathrm{CFL}_{iter}\) and \(\mathrm{CFL}_{e}\) approaches based on the iteration count and the evolution of the residual norm. It can be observed that the network speeds up the convergence significantly compared with the two reference approaches; moreover, whereas the residual norm oscillates strongly for \(\mathrm{CFL}_{iter}\), the convergence is almost monotonous for the network and \(\mathrm{CFL}_{e}\) approaches. The qualitative behavior observed in Figure 7, where the classical choices for the CFL number yield slow convergence or even strong oscillations in the residual, is not an exception; similar convergence behavior can also be observed for Couette flow and flow around an obstacle configurations; cf. sections 5.2 and 5.3. Furthermore, in Figure 8, we depict four iterates during the Newton iteration using the network approach as well as the corresponding predicted local optimal pseudo-time step. It can be observed that both the iterate and predicted local optimal pseudo-time step still vary during the first iterations and both become stationary in the later iterations.

In Figure 9, we compare the performance with respect to the numbers of Newton iterations resulting from the different choices NN, \(\mathrm{CFL}_{iter}\), and \(\mathrm{CFL}_{e}\) for the pseudo-time step on different back-step configurations; cf. section 4.3 and table 2. Furthermore, as a first investigation of the generalizability, we additionally consider cases resulting from mirroring (BM) or anti-clockwise rotating by 90 degrees (BR) the B1 cases.
see appendix A for the results for the Newton methods. This can also be seen in the average number of iterations for the configurations B1, B2, BM, and BR; see Figures 9 and A.17. For B1S, the network performs equally well as \(\text{CFL}_{iter}\).

In particular, for cases with lower velocities (B1 and B2), the network performs best. For cases with higher velocities (B1S and B2S) but the same Reynolds numbers, the network model performs clearly worse compared to the low velocity cases. Moreover, it could be observed that, for B1S and B2S, convergence was generally less robust; in particular, we observed a relatively high ratio of cases for which none of the pseudo-time step approaches converged. Since the nonlinearity should be mostly determined by the Reynolds number rather than the magnitude of the velocities, we would like to investigate in future research whether a scaling/normalization of the network inputs can improve the performance. The same behavior can also be observed for the Newton methods, although these methods converge in more cases; cf. appendix A.

Even though the mirrored and rotated configurations, of course, yield the same but mirrored respectively rotated results, we observe slight variations in the performance of the network model. This shows a slight overfitting with respect to the B1 configurations as shown exemplarily in Figure 8, where, for instance, the velocities in positive \(x\) direction are dominant; this clearly changes when the geometry is mirrored or rotated. However, the performance of the network model is only slightly worse compared to the B1 cases, which shows that the overfitting is not severe with respect to the prediction of suitable local pseudo-time steps. In the future, we may consider group invariant networks [66] to prevent these geometric overfitting effects.

Figure 11: Convergence plot of a C simulation for all three strategies for the local pseudo-time step under consideration: using the NN model as well as the strategies defined in eqs. (16) and (17).
### Couette

Next, we discuss the generalization of our hybrid globalization approach to a different type of geometries, that is, Couette cylinders. We consider two different geometries and a range of configurations resulting from varying the wall rotation velocity and refinement of the mesh; cf. table 6. In total, we obtain 45 different configurations of Couette simulations, and in table 5, the corresponding boundary conditions and maximum element sizes for the back-step configurations previously used for training the NN model can be found.

Again, we present the iterates and corresponding predicted local pseudo-time steps during the Newton iteration for one exemplary case with C geometry; see figure 10. As for the back-step configuration in figure 8, we observe that both the solution and local pseudo-time step prediction approach a stationary distribution. For this configuration, we observe a significant speed-up of Newton convergence when using the local pseudo-time step prediction by our NN, which has only been trained on back-step geometries; cf. figure 11.

Figures 12 and A.18 confirm that the use of our NN-based approach is beneficial for most of the Couette flow cases. First of all, we notice that, for the Couette flow simulations considered, all choices for the local pseudo-time step yield convergence. However, the NN-based approach accelerates convergence in \(70\,\mathrm{\char 37}\) of the C cases and even \(93.3\,\mathrm{\char 37}\) of the CS cases. The fact that the network seems to perform better for CS than for C cases, compared with the \(\mathrm{CFL}_{iter}\) and \(\mathrm{CFL}_{e}\) choices, is remarkable because the CS cases contain higher velocities and smaller mesh element sizes compared to C. For the back-step geometries, we observed the opposite; cf. figure 9, where the network performed worse on B1S and B2S than B1 and B2 configurations.

Also when taking the two variants of Newton's method into account in the results in appendix A, the NN approach remains competitive. In particular, the NN is on average close to the best method (for the C configurations) or as good as the best method

Figure 12: The results of the Couette flow for each method; the simulations are ordered from smallest required number of nonlinear iterations to largest separately for each approach. The horizontal lines indicate the average nonlinear iteration counts for the different approaches.
(for the CS configurations).

### Flow around an obstacle

As the final type of test examples, we consider Couette and back-step flow around an obstacle. In particular, we place an obstacle inside the computational domain of the back-step B1 and Couette C configurations; cf. Table 7, where three different obstacle geometries are considered: a circle with radius of \(0.03\,\mathrm{m}\), an ellipse with a-semiaxis \(0.03\,\mathrm{m}\) and b-semiaxis \(0.02\,\mathrm{m}\), and an ellipse with a-semiaxis \(0.02\,\mathrm{m}\) and b-semiaxis \(0.03\,\mathrm{m}\). We denote the resulting configurations as back-step with obstacle (BO) and Couette cylinder with obstacle (CO); see table 7 for more details on the configurations and figure 13 for exemplary flow fields for the three different obstacle geometries in Couette flow. In total, we obtain 30 BO and 45 CO configurations.

The results of the simulations with obstacles are given in Figures 16 and A.18. Qualitatively, the NN approach compares to \(\mathrm{CFL}_{iter}\) and \(\mathrm{CFL}_{e}\) similarly as before. It performs better on \(83.3\,\mathrm{\char 37}\) of the BO cases and on more than \(73.3\,\mathrm{\char 37}\) of the CO cases; it can also be observed that the CO cases are more difficult, as in more than \(11\,\mathrm{\char

\begin{table}
\begin{tabular}{|l|l|l|l|l|} \hline
**geometry type** & **x (m)** & **y (m)** & **inflow/wall velocity (m/s)** & **max. element size (m)** \\ \hline \hline
**BO** & \(0.37\) & \(0.035\!:\!0.005\!:\!0.08\) & \(0.008\) & \(0.0126\) \\ \hline
**BO** & \(0.3,1.1\) & \(0.04\) & \(0.008\) & \(0.0126\) \\ \hline
**CO** & \(0.1\) & \(0.4\) & \(0.014\!:\!0.002\!:\!0.022\) & \(0.01,0.03,0.05\) \\ \hline \end{tabular}
\end{table}
Table 7: Combinations of object center coordinates, velocities and mesh. Here, BO are B1 back-steps in table 2 with obstacles: a circle with radius \(0.03\,\mathrm{m}\), an ellipse with a-semiaxis \(0.03\,\mathrm{m}\) and b-semiaxis \(0.02\,\mathrm{m}\), and an ellipse with a-semiaxis \(0.02\,\mathrm{m}\) and b-semiaxis \(0.03\,\mathrm{m}\). The Couette flow has geometry of C1 in table 6 with obstacles: a circle with radius \(0.04\,\mathrm{m}\), an ellipse with a-semiaxis \(0.03\,\mathrm{m}\) and b-semiaxis \(0.02\,\mathrm{m}\), and an ellipse with a-semiaxis \(0.02\,\mathrm{m}\) and b-semiaxis \(0.03\,\mathrm{m}\).

Figure 13: Velocity solutions for Couette flow with different kinds of obstacles. The wall rotation is \(0.07\,\mathrm{\SIUnitSymbolMicro m}\) and the maximum element size is \(0.02\,\mathrm{m}\). The center of the obstacle is x = \(0.1\,\mathrm{m}\) and y = \(0.4\,\mathrm{m}\).
approaches converges. Also the average number of iteration is best for the NN approach for the BO and CO configurations. For one of the considered simulations of CO with an elliptic obstacle, we again plot iterates and corresponding distributions of the local pseudo-time step in figure 14. In figure 15, we plot the residual history over the Newton iteration, showing the acceleration from the NN approach.

Finally, also compared with the Newton methods in appendix A, the NN performs clearly best, that is, both in terms of robustness and in terms of acceleration of convergence.

## 6 Conclusion

We have introduced a novel NN-based approach to improve the pseudo-time stepping algorithm. Instead of using local pseudo-time steps based on predefined global CFL numbers, optimal local pseudo-time step predictions of the network are used to compute the local pseudo-time step. The network uses local information in order to make these local pseudo-time step predictions and to achieve generalizability. As a result, it can accelerate the convergence for the simulations on which it has been trained as well as for simulations which it has not seen before. In all considered simulations, the network was able to perform better or equally well compared to the standard CFL number based pseudo-time stepping strategies in most cases; the same is true in comparison with the variants of Newton's method.

In future work, we would like to extend the method to turbulent and three-dimensional

Figure 14: Iterations of CO with wall velocity 0.04 m/s and maximum element size of 0.022 m using NN, with velocity (above) and local pseudo-time step prediction (below). The obstacle is an ellipse with a-semiaxis 0.03 m and b-semiaxis 0.05 m, with center x = 0.1 m and y = 0.4 m.
flow cases. Furthermore, we will try to improve the method presented here, for example, by generating a more diverse training data set, combining the heuristic strategies for choosing the pseudo-time step with our network approach the default CFL number based strategies and the network, or removing redundancies in the input features.

## Appendix A Comparison against variants of Newton's method

In Figures 17 and 18, we show results for all configurations considered, comparing our novel NN approach with four nonlinear solvers implemented in COMSOL, that is, the two default CFL number strategies defined in eqs. (16) and (17) as well as two variants of Newton's method. As variants of Newton's method, we employ the standard Newton method with damping factor of one and an adaptive choice of the damping factor, depending on the intermediate iterates. We denote the two approaches as Newton constant (NC) and automatic Newton (AN), respectively; cf. [54, p. 1627]. Different from Figures 9, 12, and 16, the ordering of the simulations for the different methods is the same for each approach. The simulations are first organized into different blocks, and within each block the ordering is based on the number of nonlinear iterations needed for the NN approach. The blocks are as follows: before the first vertical dashed line, the NN performed best; between the first and second vertical dashed line, the NN performed best with the same number of iterations as at least one other method; between the second and third vertical dashed line, the network performed worse than at least one other method; after the third vertical line, none of the methods converged.

The results are mostly in alignment with the previous results as the NC and AN methods often perform similarly to one of the classical CFL strategies considered. In particular, when

Figure 15: Convergence plot of a CO simulation with an ellipse with a-semiaxis 0.03 m and b-semiaxis 0.05 m, with center x = 0.1 m and y = 0.4 m for the local pseudo-time step under consideration: using the NN model as well as the strategies defined in eqs. (16) and (17).
inspecting the average iteration count, both NC and AN perform quite similarly to CFL\({}_{iter}\) for B1, B1S, BM, BR, CS, and BO configurations; in all those cases the NN approach is clearly competitive. For the CO and B2 cases, AN performs better than NC, while the NN approach performs best. Only for the B2S configuration, where the NN approach generally performed worst, AC seems to perform clearly best; however, all approaches struggle for the B2S cases. Finally, for the C configurations, all methods perform similarly well, except for NC, which fails in most of the cases.

In summary, adding the NC and AN approaches to the picture does not change the qualitative comparison much; only for the B2S cases, AN performs clearly better than the NN approach. In future work, this should be taken into account for further improving the NN model.
