A Spectral Element Reduced Basis Method for Navier-Stokes Equations with Geometric Variations

We consider the Navier-Stokes equations in a channel with a narrowing of varying height. The model is discretized with high-order spectral element ansatz functions, resulting in 6372 degrees of freedom. The steady-state snapshot solutions define a reduced order space through a standard POD procedure. The reduced order space allows to accurately and efficiently evaluate the steady-state solutions for different geometries. In particular, we detail different aspects of implementing the reduced order model in combination with a spectral element discretization. It is shown that an expansion in element-wise local degrees of freedom can be combined with a reduced order modelling approach to enhance computational times in parametric many-query scenarios.

We consider the flow through a channel with a narrowing of variable height. A reduced order model (ROM) is computed from a few high-order SEM solves, which accurately approximates the high-order solutions for the parameter range of interest, i.e., the different narrowing heights under consideration. Since the parametric variations are affine, a mapping to a reference domain is applied without further interpolation techniques. The focus of this work is to show how to use simulations arising from the SEM solver Nektar++ [7] in a ROM context. In particular, the multilevel static condensation of the high-order solver is not applied, but the ROM projection works with the system matrices in local coordinates. See [1] for further details. This is in contrast to our previous work [21], since numerical experiments have shown that the multilevel static condensation is inefficient in a ROM context. Additionally, we consider affine geometry variations. With SEM as discretization method, we use global approximation functions for the high-order as well as reduced-order methods. The ROM techniques described in this paper are implemented in open-source project ITHACA-SEM1.
The outline of the paper is as follows. In Sec. 2, the model problem is defined and the geometric variations are introduced. Sec. 3 provides details on the spectral element discretization, while Sec. 4 describes the model reduction approach and shows the affine mapping to the reference domain. Numerical results are given in Sec. 5, while Sec. 6 summarizes the work and points out future perspectives.

Problem Formulation
Let Ω ∈ R 2 be the computational domain. Incompressible, viscous fluid motion in spatial domain Ω over a time interval (0, T) is governed by the incompressible Navier-Stokes equations with vector-valued velocity u, scalar-valued pressure p, kinematic viscosity ν and a body forcing f: Boundary and initial conditions are prescribed as with d, g and u 0 given and The Reynolds number Re, which characterizes the flow [8], depends on ν, a characteristic velocity U, and a characteristic length L: We are interested in computing the steady states, i.e., solutions where ∂u ∂t vanishes. The high-order simulations are obtained through time-advancement, while the ROM solutions are obtained with a fixed-point iteration.

Oseen-Iteration
The Oseen-iteration is a secant modulus fixed-point iteration, which in general exhibits a linear rate of convergence [11]. Given a current iterate (or initial condition) u k , the next iterate u k+1 is found by solving linear system: Iterations are typical stopped when the relative difference between iterates falls below a predefined tolerance in a suitable norm, like the L 2 (Ω) or H 1 0 (Ω) norm.

Model Description
We consider the reference computational domain shown in Fig. 1, which is decomposed into 36 triangular spectral elements. The spectral element expansion uses modal Legendre polynomials of order p = 11 for the velocity. The pressure ansatz space is chosen of order p − 2 to fulfill the inf-sup stability condition [12,13]. A parabolic inflow profile is prescribed at the inlet (i.e., x = 0) with horizontal velocity component u x (0, y) = y(3 − y) for y ∈ [0, 3]. At the outlet (i.e x = 8) we impose a stress-free boundary condition, everywhere else we prescribe a no-slip condition.
The height of the narrowing in the reference configuration is µ = 1, from y = 1 to y = 2. See Fig. 1. Parameter µ is considered variable in the interval µ ∈ [0.1, 2.9]. The narrowing is shrunken or expanded as to maintain the geometry symmetric about line y = 1.5. Fig. 2, Fig. 3, and Fig. 4 show the velocity components close to the the steady state for µ = 1, 0.1, 2.9, respectively. The viscosity is kept constant to ν = 1. For these simulations, the Reynolds number (6) is between 5 and 10, with maximum velocity in the narrowing as characteristic velocity U and the height of the narrowing characteristic length L. For larger Reynolds numbers (about 30), a supercritical pitchfork bifurcation occurs giving rise to the so-called Coanda effect [15,22,21], which is not subject of the current study. Our model is similar to the model considered in [9], i.e. an expansion channel with an inflow profile of varying height. However, in [9] the computational domain itself does not change.

Spectral Element Full Order Discretization
The Navier-Stokes problem is discretized with the spectral element method. The spectral/hp element software framework used is Nektar++ in version 4.4.02. The  discretized system of size N δ to solve at each step of the Oseen-iteration for fixed µ can be written as where v bnd and v int denote velocity degrees of freedom on the boundary and in the interior of the domain, respectively, while p denotes the pressure degrees of freedom. The forcing terms on the boundary and interior are denoted by f bnd and f int , respectively. The matrix A assembles the boundary-boundary coupling, B the boundary-interior coupling,B the interior-boundary coupling, and C assembles the interior-interior coupling of elemental velocity ansatz functions. In the case of a Stokes system, it holds that B =B T , but this is not the case for the Oseen equation because of the linearized convective term. The matrices D bnd and D int assemble the pressure-velocity boundary and pressure-velocity interior contributions, respectively.
The linear system (7) is assembled in local degrees of freedom, resulting in block matrices A, B,B, C, D bnd and D int , each block corresponding to a spectral element. This allows for an efficient matrix assembly since each spectral element is independent from the others, but makes the system singular. In order to solve the system, the local degrees of freedom need to be gathered into the global degrees of freedom [1].
The high-order element solver Nektar++ uses a multilevel static condensation for the solution of linear systems like (7). Since static condensation introduces intermediate parameter-dependent matrix inversions (such as C −1 in this case) several intermediate projection spaces need to be introduced to use model order reduction [21]. This can be avoided by instead projecting the expanded system (7) directly.
Next, we will take the boundary-boundary coupling across element interfaces into account. Let M denote the rectangular matrix which gathers the local boundary degrees of freedom into global boundary degrees of freedom. Multiplication of the first row of (7) by M T M will then set the boundary-boundary coupling in local degrees of freedom: The action of the matrix in (8) on the degrees of freedom on the Dirichlet boundary is computed and added to the right hand side. Such degrees of freedom are then removed from (8). The resulting system can then be used in a projection-based ROM context [18], of high-order dimension N δ × N δ and depending on the parameter µ :

Reduced Order Model
The reduced order model (ROM) computes accurate approximations to the highorder solutions in the parameter range of interest, while greatly reducing the overall computational time. This is achieved by two ingredients. First, a few high-order solu-tions are computed and the most significant proper orthogonal decomposition (POD) modes are obtained [18]. These POD modes define the reduced order ansatz space of dimension N, in which the system is solved. Second, to reduce the computational time, an offline-online computational procedure is used. See Sec. 4.1.
The POD computes a singular value decomposition of the snapshot solutions to 99.99% of the most dominant modes [10], which define the projection matrix U ∈ R N δ ×N used to project system (9): The low order solution x N (µ) then approximates the high order solution as x(µ) ≈ Ux N (µ).

Offline-Online Decomposition
The offline-online decomposition [10] enables the computational speed-up of the ROM approach in many-query scenarios. It relies on an affine parameter dependency, such that all computations depending on the high-order model size can be moved into a parameter-independent offline phase, while having a fast input-output evaluation online.
In the example under consideration here, the parameter dependency is already affine and a mapping to the reference domain can be established without using an approximation technique such as the empirical interpolation method. Thus, there exists an affine expansion of the system matrix A(µ) in the parameter µ as The coefficients Θ i (µ) are computed from the mapping x = T k (µ)x + g k , T k ∈ R 2×2 , g k ∈ R 2 , which maps the deformed subdomainΩ k to the reference subdomain Ω k . See also [16,17]. Fig. 5 shows the reference subdomain Ω k for the problem under consideration. For each subdomainΩ k the elemental basis function evaluations are transformed to the reference domain. For each velocity basis function u = (u 1 , u 2 ), v = (v 1 , v 2 ), w = (w 1 , w 2 ) and each (scalar) pressure basis function ψ, we can write the transformation with summation convention as: with The subdomain Ω 5 (see Fig. 5) is kept constant, so that no interpolation of the inflow profile is necessary. To achieve fast reduced order solves, the offline-online decomposition expands the system matrix as in (11) and computes the parameter independent projections offline, which are stored as small-sized matrices of the order N × N. Since in an Oseen-iteration each matrix is dependent on the previous iterate, the submatrices corresponding to each basis function are assembled and then formed online using the reduced basis coordinate representation of the current iterate. This is the same procedure used for the assembly of the nonlinear term in the Navier-Stokes case [18].

Numerical Results
The accuracy of the ROM is assessed using 40 snapshots sampled uniformly over the parameter domain [0.1, 2.9] for the POD and 40 randomly chosen parameter locations to test the accuracy. Fig. 6 (left) shows the decay of the energy of the POD modes. To reach the typical threshold of 99.99% on the POD energy, it takes 9 POD modes as RB ansatz functions. Fig. 6 (right) shows the relative L 2 (Ω) approximation error of the reduced order model with respect to the full order model up to 6 digits of accuracy, evaluated at the 40 randomly chosen verification parameter locations. With 9 POD modes the maximum approximation error is less than 0.7% and the mean approximation error is less than 0.5%.
While the full-order solves were computed with Nektar++, the reduced-order computations were done in ITHACA-SEM with a separate python code. To assess the computational gain, the time for a fixed point iteration step using the full-order system is compared to the time for a fixed point iteration step of the ROM with dimension 20, both done in python. The ROM online phase reduces the computational time by a factor of over 100. The offline time is dominated by computing the snapshots and the trilinear forms used to project the advection terms. See [18] for detailed explanations.

Conclusion and Outlook
We showed that the POD reduced basis technique generates accurate reduced order models for SEM discretized models under parametric variation of the geometry. The potential of a high-order spectral element method with a reduced basis ROM is the subject of current investigations. See also [6]. Since each spectral element comprises a block in the system matrix in local coordinates, a variant of the reduced basis element method (RBEM) ( [19,20]) can be successfully applied in the future.