 # Efficient computation of bifurcation diagrams with a deflated approach to reduced basis spectral element method

The majority of the most common physical phenomena can be described using partial differential equations (PDEs), however, they are very often characterized by strong nonlinearities. Such features lead to the coexistence of multiple solutions studied by the bifurcation theory. Unfortunately, in practical scenarios, one has to exploit numerical methods to compute the solutions of systems of PDEs, even if the classical techniques are usually able to compute only a single solution for any value of a parameter when more branches exist. In this work we implemented an elaborated deflated continuation method, that relies on the spectral element method (SEM) and on the reduced basis (RB) one, to efficiently compute bifurcation diagrams with more parameters and more bifurcation points. The deflated continuation method can be obtained combining the classical continuation method and the deflation one: the former is used to entirely track each known branch of the diagram, while the latter is exploited to discover the new ones. Finally, when more than one parameter is considered, the efficiency of the computation is ensured by the fact that the diagrams can be computed during the online phase while, during the offline one, one only has to compute one-dimensional diagrams. In this work, after a more detailed description of the method, we will show the results that can be obtained using it to compute a bifurcation diagram associated to a problem governed by the Navier-Stokes (NS) equations.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1. Introduction and motivation

Usually, when one wants to numerically compute a bifurcation diagrams, one has to combine many numerical methods in order to obtain it. However, since numerous solutions have to be computed and the involved system has to be solved several times, the computational cost of the task may be prohibitive in practical scenarios. In order to decrease such a cost, we decided to implement a technique based on four different methods.
Firstly, we rely on the reduced basis method [rozza_rom], this is important because, after a very expensive offline phase, the computation of a solution in the online one is, in a repetitive computational environment, very efficient. In fact, the former is used to generate a low-dimensional space defined as a combination of some of the most important solutions obtained during such a phase, while, during the latter, the solutions are sought in , and the affine decomposition of the operators can be exploited to further speed up the computation [supremizer]. This way, it is possible to efficiently compute any solution and to discretize the entire bifurcation diagram only in the online phase, significantly reducing its associated computational cost. Secondly, we decided to use the spectral element method (SEM) [sem]

to compute the solutions in the offline phase (such solutions are called snapshots or full order solutions). This is important because the solutions obtained with such a discretization method are characterized by a lower number of degrees of freedom when compared to their counterpart computed with the standard finite element method (FEM). Moreover, the computational cost of the offline phase can be further decreased using the static condensation method (also known as Schur complement

[golub]) to efficiently solve the resulting linear system [nektar_online]. Finally, we implemented an elaborated deflated continuation method to compute the snapshots in the offline phase and the complete bifurcation diagram in the online one. Such a technique is composed by two different parts: the continuation method [dijkstra] and the deflation one [classic_deflation]. The former is used to follow a known branch of the diagram, while we used the latter to compute the first solutions on unknown branches. The main idea behind the continuation method is to exploit the iterative solver and the continuity of each branch to obtain a solution very similar to the previous one and, therefore, belonging to the same branch of the latter. On the other hand, the deflation prevents the iterative solver from converging to known solutions and, this way, if it converges, it will converge to yet unknown fields.
The described approach is very general and can be used in several scenarios; among them, we decided to use it to compute a bifurcation diagram in a computational fluid dynamics (CFD) framework. In particular, the present work is strongly related to [local_rom], [martin_geom], [martin_curve_walls] and [martin0]. In fact, a possible future application of this work could be the mitral valve regurgitation [cardio]. This is a cardiac disease characterized by an inverse blood flux from the left ventricle to the left atrium. The main tool to detect and analyze such a disease is the echocardiography, however, when the blood flow undergoes the Coanda effect [tritton], it is very complex to quantify the flow rate. Therefore, it can be useful to use the direct simulation of the flow to properly analyze the images obtained via echocardiography.
We thus considered a two-dimensional channel with a narrow inlet, that represented, in a very simplified way, the mitral valve and the left atrium. Such a domain is shown, with the mesh used in the offline phase, in figure 1. Figure 1. Domain Ω and mesh used in the offline phase to compute the full order solutions

Here the left vertical wall is the inlet, the right one is the outlet and the remaining ones represent the heart walls. It can easily be noted that the mesh is very coarse, in fact only 19 elements are involved. However, thanks to the SEM, it is possible to compute very accurate solutions even with meshes similar to this one, in fact, it relies on high-order ansatz functions inside each element that allow an exponential decay of the error [hp_conv]. On the contrary, if one had computed the same solutions with the FEM, a much finer mesh would have been required because of the algebraic convergence that characterizes such a method [fem_error_b]. Moreover, since the convergence is much faster, it is possible to reach the same level of accuracy with a significantly lower number of degrees of freedom, thus further increasing the efficiency of the method, both in the assembly of the matrices and when the associated systems is solved.
Finally, it is important to remark that the results presented in section 6 could be obtained using the deflated continuation method without the reduced basis approach and with an arbitrary discretization technique. However, the computational cost of the process would be prohibitive because numerous solutions have to be computed to build a bifurcation diagram; furthermore, such a number exponentially increases with the number of varying parameters, when similar grids are used along all the directions in parameter space. On the contrary, exploiting the RB method, it is possible to compute few one-dimensional bifurcation diagrams during the offline phase and to reconstruct the other dimensions only in the online one, when each solution can be computed much more efficiently. However, it is convenient to exploit high-order methods as the SEM to speed up the expensive offline phase, in fact, one has to compute many full order solutions in order to obtain a reduced space able to capture all the branches in the online phase.
The paper is organized as follows: in section 2 we will focus on the formulation of the problem of interest, deriving the weak formulation from the strong equations and presenting the two linearization techniques [burger_notes] that will be used. Then, in section 3, the SEM will be described, with a particular focus on the static condensation method. The latter is a technique that can be used to significantly speed up the computation of the solution of the involved linear system exploiting the fact that two kinds of modes are present. In section 4, we will describe the RB method, the proper orthogonal decomposition (POD) [pod_tol_error] used to construct the reduced space and the affine decomposition [rozza_rom] exploited to ensure the efficiency of the method. Subsequently, in section 5, the continuation method and the deflation one will be discussed with a particular focus on their implementation. The results obtained with the described methods will be presented in section 6, that will be divided as follows: in the first part we will focus on a bifurcation diagram with a single parameter to prove that it can be computed during the online phase whereas, in the second part, we will highlight the efficiency of the method considering an additional parameter. Finally, we will talk about some of the future perspectives of this work in section 7.
We thus highlight that the advantages of the proposed method are strongly related to the way in which the four techniques are connected. In fact, even if the continuation and the deflation techniques allow one to compute a bifurcation diagram, they are very expensive and, without a reduced order model, the computational cost of the process may be prohibitive. Therefore, we decided to rely on the SEM and on the computation of a limited number of snapshots to perform the offline phase as efficiently as possible and to exploit the deflated continuation both in the offline and in the online phase. This way we could effectively obtain all the required solutions without exploiting any prior knowledge about the structure of the solution manifold.

We remark that the SEM is based on the open source software

Nektar++ version 4.4.0 (see [nektar]), while the reduced order model and the deflated continuation method described in this work have been implemented in ITHACA-SEM

## 2. Problem formulation

Let us consider the steady and incompressible Navier-Stokes equations [ns] in the domain :

 (1) {\uu⋅∇\uu−νΔ\uu+∇p=0in Ω,∇⋅\uu=0in Ω,

where is the velocity, the pressure normalized over a constant density and the kinematic viscosity assumed constant. In order to better characterize the flow regime, it is important to observe that its main features can be summarized in a non-dimensional quantity named Reynolds number [tritton] and defined as follows:

 Re=ULν,

where is a characteristic velocity of the flow, a characteristic length of the domain and a characteristic viscosity. System (1) can be obtained assuming that the Reynolds number is moderate and that the flow is steady.
It is important to observe that, when the equations in system (1) are normalized, they can be written in terms of non-dimensional quantities as:

 (2) {\uu⋅∇\uu−1ReΔ\uu+∇p=0in Ω,∇⋅\uu=0in Ω.

Since the structure of the mass balance equation does not change, while the only non-dimensional parameter in the momentum balance equation is the Reynolds number, one can conclude that it is the only meaningful one. Therefore, all the features present in the bifurcation diagrams that will be shown can be discussed in terms of . However, we decided to use as parameters the viscosity and a multiplicative factor on the Dirichlet inlet boundary condition, which we will denote by . It will be possible, in Section 6, to observe that the strict relation between the involved parameters and the Reynolds number allows one to explain the fact that the bifurcation points can be grouped, according to their nature, in specific and predictable curves.
To consider a well-posed problem, we supplemented system (1) with proper boundary conditions: a stress free boundary condition on the velocity at the outlet, a no-slip Dirichlet boundary condition on the physical walls and the following Dirichlet boundary conditions at the inlet

 \uu=[uv]=[20s(5−y)(y−2.5)0],x=0, y∈(2.5,5),

where is the second parameter that we included in the model and a parabolic profile has been imposed to consider a more realistic condition. It should be observed that the dimension of the domain, described in figure 1, and the constant in front of the parabolic profile in the inlet boundary condition are used to obtain values of in the desired range when and .
Furthermore, we introduce the variational formulation [segal], required by the Spectral element method (SEM) to obtain a discrete solution of the Navier-Stokes system. This will also be fundamental for the applicability of the Reduced basis method (RB). When deriving the weak formulation, one has to set appropriate functional spaces:

 (H10,∂ΩD(Ω))2={\ww∈(H1(Ω))2∣\ww=0 on ∂ΩD},
 L20(Ω)={q∈L2(Ω)∣∫Ωqdx=0},

where is the portion of the boundary where Dirichlet boundary conditions are imposed.
Moreover, one multiplies the equations in system (1) by the appropriate test functions, integrates over the entire domain and integrates by parts the integrals associated to and obtaining the following weak formulation: find and such that

 (3) ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩ν∫Ω∇\uu⋅∇\vv+∫Ω(\uu⋅∇\uu)\vv−∫Ωp∇⋅\vv=0∀\vv∈(H10,∂ΩD(Ω))2,∫Ωq∇⋅\uu=0∀q∈L20(Ω).

Denoting as and as , and introducing the following bilinear and trilinear forms:

 a(\vv,\ww)=ν∫Ω∇\vv⋅∇\ww ∀\vv,\ww∈V, b(\ww,q)=−∫Ω(∇⋅\ww)q ∀\ww∈V, ∀q∈Q, c(\uu,\vv,w)=∫Ω(\uu⋅∇\vv)\ww ∀\uu,\vv,\ww∈V,

problem (3) can be expressed, in a more compact way, as follows: find and such that

 (4) {a(\uu,\vv)+c(\uu,\uu,\vv)+b(\vv,p)=0∀\vv∈V,b(\uu,q)=0∀q∈Q.

Such a notation will be useful in section 3.1 to describe the static condensation method and we are going to exploit it to easily describe the two most common linearization techniques [burger_notes].
The first one, named Oseen iteration, relies on the fact that, when the iterative solver is converging, two subsequent approximations are very similar, therefore, the nonlinear term can be approximated as follows:

 c(\uuk+1,\uuk+1,\vv)≈c(\uuk,\uuk+1,\vv),

where is the approximation obtained in the last iteration, while is the unknown one. Such a linearization technique is very common because it can be easily implemented and the associated iterative solver is very stable, but the convergence is only linear. On the other hand, one can exploit the Newton method to obtain a quadratic convergence, although a more accurate initial guess is required. To derive the latter, one expresses the unknown approximation as:

 \uuk+1=\uuk+δ\uu,

where is the variation between the unknown solution and the last computed one . This way, it is possible to approximate the nonlinear term as follows:

 c(\uuk+1,\uuk+1,\vv)≈c(\uuk+1,\uuk,\vv)+c(\uuk,\uuk+1,\vv)−c(\uuk,\uuk,\vv).

In this work we exploited both techniques in order to increase the effectiveness of the deflation method (see Section 5.2). In fact, the Oseen iteration ensures a too slow convergence but, using only the Newton method, the iterative solver often diverges when solving the deflated problem.

## 3. The spectral element method

In this section the main features of the SEM [canuto_spettrali] will be described; since we are interested in the efficiency of the method, the main focus will be on the static condensation method. This is a technique that can be used to significantly reduce the computational cost to get the solution of the obtained linear system.
Let us consider the following Galerkin formulation [canuto_spettrali2], derived from problem (4): find and such that:

 (5) {a(\uu,\vv)+c(\uu,\uu,\vv)+b(\vv,p)=0∀\vv∈Vδ,b(\uu,q)=0∀q∈Qδ,

where and are, respectively, two finite dimensional approximations of and . It is thus possible to consider the bases and associated to the two components of the velocity and associated with the pressure. Therefore, discrete velocity and pressure can be expressed as follows:

 u=Nu∑i=1uiϕui,v=Nv∑i=1viϕvi,p=Np∑i=1piϕpi,

where , and are scalar coefficients that characterize the velocity and pressure fields. Moreover, in the SEM the basis functions are high-order polynomials inside the associated elements [hp_fem]. In particular, in this work we decided to use the stable pair , i.e. the velocity is represented by a polynomial of order while the pressure by one of order inside each element (see [p-p-21] and [p-p-22] for a more detailed explanation of the approach). The high order of the polynomials implies two main consequences: firstly, since several degrees of freedom are associated to each element, it is possible to obtain accurate solutions even with very coarse meshes like the one showed in figure 1. In general, this is convenient because generating fine meshes as the ones required by the FEM is an expensive operation [mesh_refinement]. Secondly, this way it is possible to ensure the exponential convergence of the method instead of the algebraic one that characterizes the FEM [hp_refinement] when the solution is smooth enough. This is important because the same level of accuracy can be obtained with a significantly lower number of degrees of freedom and, therefore, the assembly of the linear system and the computation of the solution are much more efficient. Furthermore, exploiting the high order polynomials, it is possible to rely on more accurate differentiation and integration formulas [gauss_quadrature].

### 3.1. The static condensation method

The static condensation method is a technique that allows one to solve a linear system much more efficiently, exploiting a specific structure of the associated matrix. It should be noted that, even if such a technique may be used to exploit a domain decomposition approach in a reduced basis framework [sc_rb], in this work we used it only in the offline phase to increase the efficiency of the SEM and, therefore, we did not link it with the RB method. Moreover, we remark that we will outline the method as described in [nektar_online].
To obtain the required structure, while solving the Navier-Stokes equations with the SEM, one has to split the velocity degrees of freedom into different groups. The first one contains the interior modes, i.e. all the basis functions with support inside a single element, while the second one contains all the remaining basis functions, that will be denoted as boundary modes [sem]. It is crucial to observe that two interior modes associated to two different elements are orthogonal to each other because the measure of the intersection of their supports is zero. Such a property can be exploited to properly sort all the degrees of freedom to obtain block matrices with very small blocks associated to the elements. Denoting as

the vector of the velocity degrees of freedom associated to the interior modes, as

the vector associated to the boundary ones and as the one associated to the pressure modes, it is possible to expand the linear system associated to problem (5) as follows:

 (6) ⎡⎢ ⎢⎣A−DTbndB−Dbnd0−Dint˜BT−DTintC⎤⎥ ⎥⎦⎡⎢⎣\uubndp\uuint⎤⎥⎦=⎡⎢⎣fbnd0fint⎤⎥⎦.

It can be observed that the elements of the submatrix represent the relations between pairs of boundary modes, while the ones of the submatrices and are, respectively, associated to boundary-interior and interior-interior pairs. Finally, the interaction between the pressure and the velocity is summarized in and . The elements of such matrices can be computed as follows:

 A[i][j]=c(ϕjbnd,\uuk,ϕibnd)+c(\uuk,ϕjbnd,ϕibnd)+a(ϕjbnd,ϕibnd), ∀i,j=1,…,Nbnd B[i][j]=c(ϕjint,\uuk,ϕibnd)+c(\uuk,ϕjint,ϕibnd)+a(ϕjbnd,ϕiint), ∀i=1,…,Nbnd, ∀j=1,…,Nint ˜BT[i][j]=c(ϕjbnd,\uuk,ϕiint)+c(\uuk,ϕjbnd,ϕiint)+a(ϕjbnd,ϕiint), ∀i=1,…,Nint, ∀j=1,…,Nbnd C[i][j]=c(ϕjint,\uuk,ϕiint)+c(\uuk,ϕjint,ϕiint)+a(ϕjint,ϕiint), ∀i,j=1,…,Nint Dbnd[i][j]=b(ϕjbnd,ϕip), ∀i=1,…,Np, ∀j=1,…,Nbnd Dint[i][j]=b(ϕjint,ϕip), ∀i=1,…,Np, ∀j=1,…,Nint fbnd[i]=f(ϕibnd)+c(\uuk,\uuk,ϕibnd), ∀i=1,…,Nbnd fint[i]=f(ϕiint)+c(\uuk,\uuk,ϕiint), ∀i=1,…,Nint

where , and are, respectively, the number degrees of freedom associated to the velocity boundary modes, to the velocity interior ones and to the pressure ones. Moreover, the terms and are associated to the external forces that act on the system. However, since we assumed their absence, these two terms can be neglected.
It should be observed that such expressions hold when the Newton method is employed, however, when using the Oseen one, the first terms in the expansions of , , and and the last ones in and have to be discarded. Moreover, it should be noted that the submatrix is block diagonal and, therefore, it is easy to invert. Denoting

as the identity matrix, we can premultiply the previous system by the matrix:

 K=⎡⎢⎣I0−BC−10IDintC−100I⎤⎥⎦,

obtaining the following one:

 (7) ⎡⎢ ⎢ ⎢ ⎢⎣A−BC−1˜BT−DTbnd+BC−1DTint0−Dbnd+DintC−1˜BT−DintC−1DTint0˜BT−DTintC⎤⎥ ⎥ ⎥ ⎥⎦⎡⎢⎣\uubndp\uuint⎤⎥⎦=⎡⎢ ⎢⎣fbnd−BC−1fintDintC−1fintfint⎤⎥ ⎥⎦.

It is important to observe that, this way, the third equation has been decoupled from the other ones and the associated unknowns can be easily obtained after having solved the remaining block. Let us focus on the first block system involving and , that can be written, simplifying the notation and considering a new set of matrices, as:

 (8) [^A^B^C^D][b^p]=[^fbnd^fp],

where is a vector that contains both and the mean pressure (or a degree of freedom associated to it) while accounts for the remaining pressure coefficients. A second level of static condensation can be obtained repeating the previous steps, indeed it is possible to modify equation (8) into the following one:

 (9) [^A−^B^D−1^C0^C^D][b^p]=⎡⎣^fbnd−^B^D−1^fp^fp⎤⎦.

This way, one can solve the first row of equation (9) to obtain , then one can substitute it into its second row to obtain and, finally, these quantities can be used with the third row of problem (7) to compute . The computational cost of such a process is significantly lower than the direct computation of the solution of problem (6), in fact is the only unknown vector that has to be computed solving a linear system.

## 4. The reduced basis method

In this section the RB method [rozza_rom] will be briefly described. Such a technique can be used to efficiently compute the solution of a system of PDEs and, therefore, it is often used in optimization, real time queries, optimal control, design and uncertainty quantification. In this work, since we are interested in the discretization of bifurcation diagrams with many parameters, it has been used to compute the numerous required solutions. It should be observed that the same result could be obtained without exploiting the RB method, but the computational cost would be prohibitive in this case. Let us consider the following abstract parametric problem: given , find such that

 (10) a(u(μ),v;μ)=f(v;μ)∀v∈V,

where is a vector of scalar parameters, is a symmetric, coercive, bilinear and continuous operator for any parameter in the parameter space . Analogously, is a linear and continuous operator for any . It is then possible to consider the discrete version of problem (10), that reads as follows: given , find such that

 (11) a(uδ(μ),v;μ)=f(v;μ)∀v∈Vδ,

where is a finite dimensional space that depends on the chosen discretization method. We remark that similar discrete and parametric problems have been deeply investigated, we thus refer to [canpar1], [canpar2] and [canpar3] for a more detailed analysis.
In this work the employed discretization method is the SEM, we can thus obtain accurate solutions with a space smaller than the one used in the FEM. However, in order to significantly increase the efficiency of the solver, one would like to use spaces described by only few basis functions. Defining as the finite-dimensional space used in the online phase of the reduced basis method, and assuming that , one can exploit the same formulation of problem (11) to define a reduced problem that can be solved much more efficiently: given , find :

 (12) a(uδ(μ),v;μ)=f(v;μ)∀v∈Vrb.

It is thus necessary to be able to generate a small space capable to accurately discretize the continuous solutions of problem (10). To do so, one can exploit problem (11) and different numerical techniques. In this work we only used the Proper Orthogonal Decomposition (POD) [rozza_rom] to construct such a reduced space because its downsides are balanced by the advantages of the deflated continuation method [tau_deflation] that will be described in Section 5. In fact, two of the main disadvantages of the POD method are related to the high computational cost of the offline phase (indeed many solutions have to be computed) and to the fact that the parameter space has to be properly sampled. On the other hand, the deflated continuation method allows one to efficiently compute the required solutions and to automatically select only proper values of the parameter.

### 4.1. The proper orthogonal decomposition

The POD can be used to generate a reduced space that is optimal in the norm, in fact is the -dimensional space that minimizes the following quantity [pod_tol_error]:

 √1M∑μ∈PMinfvrb∈Vrb∥uδ(μ)−vrb∥2V,

where is a finite sampling of of dimension . In order to construct one considers a symmetric and linear operator defined as follows:

 C(vδ)=1MM∑i=1(vδ,ψi)Vψi,vδ∈VM,

where and .

Then, one computes the eigenvalue-eigenfunction pairs

s.t. for any defined as:

 (C(ξi),ψj)V=λi(ξi,ψj)V,1≤j≤M.

Finally, it is possible to sort the eigenvalues in descending order and the associated eigenfunctions accordingly, and generate the reduced space with the first eigenfunctions. This way, it is possible to prove that the error obtained approximating the solutions of with the ones in is given by:

 (13)  ⎷1MM∑i=1∥ψi−PNrb(ψi)∥2V= ⎷M∑i=Nrb+1λi,

where is a projection operator over defined as

 PNrb(v)=Nrb∑i=1(v,ξi)Vξi.

Therefore, from the discrete point of view, in order to obtain a reduced space with the desired properties, one has to compute the eigenvalues and the eigenvectors of the correlation matrix of the snapshots and has to express the latter as linear combination of full order basis functions. One of the most efficient method to perform such a task is the Singular Value Decomposition (SVD)

[fast_svd].

### 4.2. The affine decomposition

In order to efficiently solve problem (12), it is important to seek the solution in a lower dimensional space, as the one that can be obtained through the POD method, whose accuracy is ensured thanks to property (13). However, in order to benefit of this low dimensionality, by splitting the computation in the two phases, the model has to fulfill some additional hypothesis. In particular, in order to have an online phase which is independent from the number of degrees of freedom of the approximation, namely , one can exploit the so-called affine decomposition [supremizer2].
The main idea is to precompute several matrices and vectors in the offline phase that will be used for every instance of a new parameter, during the online one, to rapidly assemble the linear system. This is important because, this way, such an assembly does not scale with the dimension of and all the operations of the online phase only depend on . Let us denote the linear system that has to be solved in the online phase as

 Arb(μ)urb=frb(μ).

Using this notation, the affine decomposition assumption reads as:

 Arb(μ)=Qa∑q=1θqa(μ)Arbq,
 frb(μ)=Qf∑q=1θqf(μ)frbq.

It is important to observe that and are expressed as linear combination of matrices or vectors that do not depend on the parameter and, therefore, its contribution is entirely included in the coefficients. This way it is possible to precompute for any and for any in the offline phase and, subsequently, assemble and

computing only the scalar coefficients. However, it should be noted that it is not always possible to directly exploit such a structure but, when required, it can be approximated with the so called Empirical interpolation method (EIM)

[eim].

## 5. Numerical computation of bifurcation diagrams

Let us consider the following nonlinear parametric equation:

 (14) L(\uu;μ)=0,

where the function belongs to the functional space , the parameter belongs to and is a nonlinear operator. Since, due to the nonlinearity, several solutions can exist for the same value of the parameter, it is important to summarize them in a single diagram in order to highlight the main properties of the solutions manifold associated to the parameter space. Such diagrams are called bifurcation diagrams [kielhofer], the information is often summarized by means of a scalar output function, while the parameter is represented on the other axis. It should be noted that, if , more than one axis are required to properly represent the parameter and it could be useful to rely on more advanced visualization techniques to show the diagram [bif_diagr_nd].
In this work we decided to compute the bifurcation diagrams with a combination of the continuation method and the deflation method. The main idea behind the coupling of this two techniques is that, on the one hand, one can entirely discretize a specific branch of the diagram given a solution belonging to it while, on the other hand, the deflation method is used to compute the first solutions of the new branches that are required by the continuation method. Such an approach, where the two described techniques are alternated in order to discover and follow each branch of the diagram, is called deflated continuation and is described in [tau_deflation]. However, we implemented a more elaborated version of it, in fact we decided to use two different versions of the continuation method and to pair a novel approach to the deflation one in order to increase the efficiency and the effectiveness of the deflated continuation method. Such modifications will be better explained in the following sections.

### 5.1. Continuation method

In order to simplify the notation in the discussion, let us assume that, in problem (14), there exists a solution for any value of the parameter and that . The first assumption is useful to consider well-posed problems, while the second one simplifies the description of the method. However, the fact that is not restrictive, indeed one can consider and let vary only one component at a time. This way the discussed approach can be extended to problems characterized by more parameters.
Since problem (14) is nonlinear, an iterative solver is required to compute a solution. However, in order to ensure its convergence, it needs an initial guess close enough to the sought solution, such a guess can be obtained with the continuation method. The aim of the technique, in fact, is to compute a proper initial guess in order to allow the solver to converge to a solution on the same branch of the last one. This way, any arbitrary branch can be entirely reconstructed repeating several times the following procedure [seydel].
Initially, one assumes to know solutions () on a branch of the diagram, they are associated to the parameter values and one wants to compute the solution associated to the value . Given the input , the output of the continuation method is a function close to the unknown solution . Using as an initial guess, the iterative solver can efficiently converge to . Then, the process can be repeated to compute exploiting and the previous solutions, and so on.
In this work we implemented two different versions of the continuation method. The first one, that we will denote as simple continuation, is the simplest way to exploit the already computed information to obtain an initial guess and is characterized by . On the other hand, the pseudo-arclength continuation [keller] is a more advanced technique and exploits the last two solutions to compute the subsequent initial guess.
Let us first consider the simple continuation. Exploiting the continuity of each branch, one can assume that, if the step size is very small, the solutions and will be very similar to each other. This way, can be considered a good approximation of and used as initial guess to compute it. The advantages of such an approach are that it is inexpensive and requires a single solution, however, properly setting the quantity

is complex, even if it is possible to rely on prior knowledge or heuristics.

Unfortunately, such a choice is very important, in fact a too wide step may make the solver diverge or converge to solutions on other branches, because of the significant difference between and . On the other hand, a too short step would ensure the convergence to the correct solution but the computational cost would dramatically increase [allgower].
In order to properly set the step size and to improve the accuracy of the obtained initial guess, one can choose to rely on the pseudo-arclength continuation method. In such a technique the next value of the parameter is considered as an unknown and an alternative parametrization of the branch, characterized by the curve arclength , is taken into account. To derive the system that has to be solved, as explained in [keller], let us consider the following normalization equation:

 (15) N(\uu,μ;ΔSi)≐˙\uuTi(\uu−\uui)+˙μi(μ−μi)−ΔSi=0,

where is a point on a regular portion of the branch and is the unit tangent to the curve in such a point. Equation (15) characterizes the plane orthogonal to the vector such that the distance between and its projection on the plane is . Moreover, if the line described by is a good approximation of , the orthogonal projection of on the plane is very similar to the sought solution . A representation of such a process is outlined in figure 2.

Consequently, this projection can be used as a good initial guess. In order to compute it, one can solve the linear system:

 (16) [Li\uuLiμ˙\uui˙μi][Δ\uuiΔμi]=[LiΔSi],

where the subscripts are associated to the derivation operation, the superscripts represent the point where the function is evaluated, i.e. and , and the following notation has been used:

 Δ\uui=˜\uui+1−\uui,Δμi=μi+1−μi.

This system can be obtained by linearizing the following one, obtained combining equation (14) with equation (15), with the Newton method:

 (17) {L(\uu;μ)=0,˙\uuTi(\uu−\uui)+˙μi(μ−μi)−ΔSi=0.

Furthermore, since the quantities and are, in general, not available, they have to be approximated. We decided to use the following approximations:

 ˙\uu≃\uui−\uui−1ΔSi−1,˙μ≃μi−μi−1ΔSi−1,

even if several alternatives exist.
The main advantage of such a version of the continuation method is that the subsequent value of the parameter is automatically chosen. This way, the steps are wider in very smooth regions, while they can be much shorter near the singularities. This is important because, when one wants to compute a bifurcation diagram, there are regions where the solution varies very rapidly, and regions where two solutions are very similar even if they are associated to two values of the parameters very far from each other. We decided to iteratively modify to further improve the effectiveness of the method, however, good results can be obtained also fixing it after some experimental observations.
In general, it is possible to observe that such regions are, respectively, close to and far from a bifurcation point. Moreover, the pseudo-arclength continuation is more accurate than the simple one because it relies on a branch linearization. However, sometimes such a technique cannot be used. This issue can arise in two different scenarios. Firstly, when one wants to compute the second solution, only a single solution is available and, therefore, the technique cannot be applied. Secondly, right after a bifurcation point, the number of solutions varies from one iterations to the next one. This implies that it is not possible to exploit two solutions on the same branch to compute the initial guess. In such scenarios, we decided to use the simple continuation with a step size proportional to the last value of after a bifurcation point and with a very short step after the computation of the first solution.
Finally, we highlight that, since the structure of the matrix associated to problem (16) is different from the one of problem (6), we decided to implement a bordering algorithm [keller] to solve the former with the static condensation method.

### 5.2. Deflation method

As discussed in the previous section, the continuation method allows one to follow each branch of the diagram, however, it requires a first solution on the branch to properly work. In this work such solutions have been computed in two different ways, in fact we decided to compute the very first solution using the zero initial guess because of the lack of prior knowledge, while we used the deflation method to obtain the first solutions on unknown branches. Such a technique has been initially developed to compute multiple roots of a polynomial and, in order to explain it, it is convenient to discuss, as in [classic_deflation], a simplified scenario first. Let us consider the scalar polynomial

 (18) p(x)=c0m−1∏j=0(x−xj),

where is a scaling factor and are distinct roots. Moreover, let us assume to be able to numerically compute just a single solution for any arbitrary polynomial with a numerical method. This way, one can easily compute root but one cannot obtain the other ones. However, to overcome such an issue, it is possible to consider the deflated polynomial

 (19) p1(x)=c01(x−xi)m−1∏j=0(x−xj).

It should be observed that the polynomial in equation (19) is characterized by the same roots of the one in (18) except for the -th one. It is therefore possible to exploit the available algorithm to obtain a root of , consequently obtaining the second one of . Such a process can be repeated in order to obtain all the existing roots of , however, it should be remembered that, when such roots are computed with a numerical method, the obtained results are approximations of the exact ones and, therefore, the deflated polynomial is still characterized by all the roots of the previous one. Anyway, if the approximation is accurate enough, the numerically deflated polynomial is very similar to the analytical deflated one and the difference is significant only in a very small neighbourhood of the deflated root. It is thus possible to exploit such a technique to compute distinct roots that are far from each other.
The described approach can be then generalized to be applied to systems of PDEs, however, before discussing the generalized method, it is convenient to introduce the concept of deflation operator. In equation (19) the function is called deflation operator and it is responsible for removing a root from the polynomial; such an operator can be generalized in the following way (see [classic_deflation]). Let us denote by and two Banach spaces and by an open subset of the additional Banach space . Moreover, let be a Fréchet differentiable operator and be its Fréchet derivative. Then, let be an invertible linear operator for each and for each . If, for any Fréchet differentiable operator for which the following properties hold:

 L(\ww)=0,L′(\ww) is nonsingular,

and for any arbitrary sequence such that , the following inequality holds

 (20) liminfi→∞∥M(\uui;\ww)L(\uui)∥Z>0,

then is a deflation operator. In order to generalize equation (19), one considers the deflation operator and the following deflated system:

 (21) G(\uu;μ)≐M(\uu;\ww)L(\uu;μ)=0,

that is characterized by the same solutions of except for .
Denoting by the identity operator, it is possible to consider the most straightforward generalization of the deflation operator introduced in equation (19):

 M(\uu;\ww)=I∥\uu−\ww∥U.

It can be observed that, since such an operator tends to zero when is very far from , the iterative solver may converge to unphysical solutions if it is directly employed because the exact residual is multiplied by a factor that tends to zero. Even if several alternatives exist, we preferred to implement a very simple deflation operator, that can be expressed as:

 M(\uu;\ww)=I+I∥\uu−\ww∥pU,

with in the offline phase and in the online; such quantities have been fixed via experimental observations.
The main advantage of the deflation method consists in the ability to discover unknown branches without any prior knowledge, however, if other branches exist, one cannot be sure that they will be found with such a technique. In fact, if a branch is too far from any known solution, the solver may diverge before reaching the region of attraction of any solution in . Therefore, it is advisable to fix a meaningful maximum number of iterations for the iterative solver when the deflated system is solved. Such a threshold should be high enough in order to let enough time to the iterative solver, however, if too many iterations are available, a lot of computational resources will be wasted when new branches cannot be found.
We decided to fix such a quantity equal to 150 in the offline phase and to 300 in the online one (because each iteration is significantly less expensive). Consequently, the deflation method is the bottleneck of the deflated continuation when applied at each step because the linear system has to be assembled and solved many times. Therefore, if one knows the position of the bifurcation points, it is advisable to use the deflation only in those regions.
Moreover, as described in [tau_deflation], it is possible to increase the efficiency of the deflation observing that one does not need to assemble and solve the linear systems associated to the deflated problem. In order to describe a more efficient approach, let us assume that one wants to solve system (21) with the Newton method, that is characterized by the following iteration:

 JGΔ\uuG=−G,

where is the associated Jacobian. In the same way, we will refer to the corresponding iteration on the undeflated system as:

 (22) JLΔ\uuL=−L.

Assuming that and exploiting the Sherman-Morrison formula [golub], one can derive the following relation:

 Δ\uuG =−J−1GG =−(MJL+LM′T)−1(ML) =−(M−1J−1L−M−1J−1LLM′TM−1J−1L1+M′TM−1J−1LL)(ML) =(1−M−1M′TJ−1LL1+M−1M′TJ−1LL)(−J−1LL).

One can observe that the quantity that multiplies the term is a scalar quantity, therefore we will refer to it as

 τ≐1−M−1M′TJ−1LL1+M−1M′TJ−1LL=11+M−1M′TJ−1LL.

Finally, if one is able to efficiently compute a solution of the undeflated system, it is possible to exploit the same algorithm to obtain a solution of the deflated one. In fact, one simply has to multiply each Newton iteration for , that can be computed, after having solved system (22), very efficiently with a scalar product and scalar operations. Its cost is thus linear with the number of degrees of freedom. On the other hand, if one wants to directly solve the deflated system, one has to explicitly construct the matrix , that is full, and solve the associated linear system without exploiting the sparsity of the original matrix .

#### 5.2.1. An approach to improve the deflation method

In the previous section we explained how to efficiently deflate a system modifying the residual of the iterative solver associated to the undeflated problem instead of constructing and solving the deflated one. However, we analyzed the values assumed by the new scalar factor and we observed that they were always very low (very often its absolute value was lower than . This phenomenon implied serious consequences because, using the formula

 Δ\uuG=τΔ\uuL,

one observes that when . Consequently, the iterative solver would have required too many iterations and the computational cost would have been prohibitive. We thus analyzed the behaviour of and its relation with the solutions exploiting the fact that, since , it cannot change the direction of the undeflated residual. However, even if the direction cannot differ between and , their orientation can because can be positive or negative. Therefore, we decided to modify the values of while maintaining its sign to avoid losing changes in the orientation.
Firstly, we decided to fix two lower bounds in order to prevent to assume values too close to zero through the following formula:

 (23) τ=⎧⎪⎨⎪⎩τ−tif τ>τ−t,τ+tif τ<τ+t,τotherwise .

In this work we decided to use and , however, such values are very problem dependent. Secondly, we multiplied by a scaling factor in order to avoid using the thresholds too often. It is important to observe that we initially fixed and then we multiplied it by 1.75 each time switched from negative to positive. In fact, this change of sign implies that the deflation was preventing the solver to converge to a solution () but, when the current iteration was too far from it, the solver was attracted again by (). Instead, we observed that bigger values of could help the solver to better escape from the region of attraction of a solution and, therefore, we slightly increased it each time such an escape failed.

 τ=⎧⎪⎨⎪⎩τ−tif cτ>τ−t,τ+tif cτ<τ+t,cτotherwise .

Finally, we decided to use the Newton iteration when the current iteration was close to known solutions, while we exploited the Oseen one otherwise to improve the stability of the solver.

## 6. Numerical results

In this section we will discuss the results that can be obtained using a combination of the described techniques. However, since one of the hypothesis of the reduced basis method is the presence of a smooth solution manifold with just a single solution associated to each parameter value that, in this work, does not hold because of the pitchfork bifurcation points, we will first prove that such an approach is able to accurately discretize a bifurcation diagram with a single parameter. For instance, in figure 3, a bifurcation diagram where such hypothesis do not hold is shown. Subsequently, we will move on to the description of a bifurcation diagram with two parameters, where we will be able to appreciate the efficiency of the computation. Figure 3. Bifurcation diagram with multiple solutions associated to the same parameter value

In order to better understand the different existing solutions, it is convenient to analyze the ones in figure 4, where the horizontal velocity is highlighted by the colour gradient. It is important to observe that more than one solution is associated to the values 0.6 and 0.3 of the viscosity. Moreover, some of them are axisymmetric while other ones are not, however, due to the symmetry of the domain and of the boundary conditions, it is always possible to reflect a solution over the horizontal axis of symmetry to obtain another solution. Such a phenomenon can be observed, for instance, in figures sol2 and sol3. This is important because the function that will be used to compute the bifurcation diagram is the following one:

 f(\uu)≐\mypm∫Ω∥\uu−R(\uu)∥2,

where is the velocity, is an operator that reflects a solution over the horizontal symmetry axis and the sign is positive one if the jet hugs the upper wall or negative one otherwise. Therefore, will be equal to zero if a solution is perfectly symmetrical, while its absolute value will increase with the asymmetry of the velocity field. Exploiting such an interpretation, one can immediately conclude that the bifurcation diagrams will be symmetrical because each solution can be mirrored over the symmetry axis obtaining another solution such that .

### 6.1. Results with a single varying parameter

As previously written, in this section we will discuss the results that can be obtained with only a single parameter. Therefore, we will not consider the efficiency of the computation but we will focus on the application of the described techniques to accurately discretize a bifurcation diagram. Let us remember that, in order to obtain a bifurcation diagram during the online phase, we performed the following steps. Firstly, the SEM (see section 3) has been used to compute the snapshots involved in the reduced space generation. In order to choose the proper values of the parameter and obtain multiple solutions for any , we respectively exploited the continuation and the deflation methods (see section 5). Subsequently, the snapshots are grouped to construct the reduced space with the POD (see section 4) and, finally, the continuation and the deflation are used again in the online phase to efficiently reconstruct the diagram.
The first diagrams that we want to show represent the decay of the eigenvalues of the correlation matrix obtained in the POD method. In figure svdecay1 the decay associated to the reference setting is shown; it consists in 24 snapshots distributed on three different branches of a bifurcation diagram with a single bifurcation point. It can be observed that the decay is exponential even if there is a singular point (for instance, similar behaviours have been proved in [exp_decay] and in [kolmogorov] without considering singular points). Moreover, it should be noted that the decay is very similar when the continuation steps are smaller and, therefore, all the snapshots are closer to the bifurcation point (figure svdecay2). On the other hand, the decay remains exponential but it is faster or slightly different if, respectively, the snapshots belong to the same branch or if we increase their number from 24 to 100. Such a result is important because an exponential decay of the eigenvalues implies, thanks to relation (13), that the approximation error of the reduced spaces exponentially decreases with respect to its own dimension.

The next figure shows the entire one-dimensional bifurcation diagram on which we will discuss (figure 6), it is computed during the offline phase and each point corresponds to a snapshot. Figure 6. Bifurcation diagram computed during the offline phase, each point is associated to a snapshot and each branch is characterized by a different colour and a different marker

It can be observed that it includes two bifurcation points and five different branches. In order to compute it, we decided to use as first solution the one obtained for and, then, we decreased the viscosity computing the solutions with the deflated continuation method. On the other hand, we highlight that, if one is not interested in the entire diagram but only in the position of the bifurcation points, it is convenient to perform an eigenvalues analysis as discussed in [pichi_bif] and in [eigen_study].
Since, as discussed in the introduction of this section, the output functional can be considered as a measure of the asymmetry of the solutions, one can observe that the solutions in figures sol11, sol12 and sol13 are associated to the middle branch. Moreover, the ones in figures sol2 and sol7 are two characteristic velocity fields of the first and of the second part of the upper branch that is born from the first bifurcation point (the one on the right, since we decrease the viscosity).
Consequently, the solution in figures sol3 and sol6 are associated to the lower branch because they are the mirrored solutions of the latter and, finally, the fields in figures sol4 and sol5 are representative of the solutions on the last two branches. Obtaining such a diagram is very expensive because many different solutions have to be computed with the full order solver. In fact it has been obtained using 224 snapshots with 7372 degrees of freedom and, to obtain them, each step of the iterative solver spent almost 0.67 seconds. We remark that with 7372 degrees of freedom we were able to compute accurate solutions thanks to the SEM. However, if one uses other discretization techniques or requires a very high accuracy, the computational cost can significantly increase. For instance, in [local_rom], a very similar model has been computed with the FEM exploiting 90876 degrees of freedom to obtain the desired accuracy.
It is important to note that, even if the computation of this diagram is very expensive, such a result can be obtained during the online phase much more efficiently. The obtained diagram is shown in figure 7, it can be seen that it is possible to reconstruct both the bifurcation points and all the branches. Moreover, since in this simulation the solutions were associated to 37 degrees of freedom, the iterative solver approximately spent only seconds for every iteration with a speedup of about 4 order of magnitude. The statistics in table 1 can be used to observe that, if a bifurcation diagram is discretized with at least 770 solutions, then the described approach is more efficient than the direct computation of a diagram without the RB method. Finally, we remark that these quantities have been obtained exploiting a prior knowledge obtained from previous experiments about the position of the bifurcation points to use the deflation only when new branches can be found. However, without such a knowledge one should use the deflation at each step of the continuation, extremely increasing the computational cost of the process and making the described approach even more convenient. Unfortunately, because of the data structures used in Nektar++ to perform static condensation, some of the computations of the online phase depended on and, therefore, the time required to compute, on average, a solution given the previous ones significantly increased ruining the speedup obtained in a single iteration. In fact the online solver could generally converge in 4 iterations, but the average time obtained dividing the total required time needed to get the complete diagram by the number of computed solutions is much higher than 4 times the time required to perform a single step of the iterative solver.

In order to quantify the accuracy of the obtained result, we decided to perform an empirical error analysis. Thus we reprojected the reduced solutions on the full order space and used them as initial guesses for the iterative solver. This way we compared the obtained full order solutions with the reduced ones and we computed the associated relative error that is shown in figure 8. In order to properly interpret the diagram, it is important to remark that the tolerance used by the offline iterative solver was and, therefore, the relative error is almost always very close to such a quantity. However, coherently with the fact that one of the hypothesis of the reduced basis method regards the presence of a smooth solution manifold, the error increases in the neighbourhoods of the bifurcation points.

Moreover, it is possible to observe the behaviours of the average and maximum relative errors in table 2. The first column is associated to a neighbourhood of the bifurcation point near , the second one to the portion of the domain between the singular points, the third to a neighbourhood of the other bifurcation point and, finally, the fourth to the entire diagram. It can be noted that, even if the bifurcation points strongly affect the accuracy, the global average error is only one order of magnitude above the solver threshold.

It is interesting to note that the error exponentially decreases with respect to the dimension of the reduced space, both in smooth regions and close to the singular points. This phenomenon can be observed in figure 9, where the error is associated to solutions in a small neighbourhood of the first bifurcation point. The fact that the maximum error decreases almost as the average one implies that the accuracy of the less accurate solution exhibits the same behaviour as the other ones. Figure 9. Relative errors over the dimension of the reduced space

### 6.2. Results with two parameters

Finally, in this section we will address the problem of efficiency. In the previous section we showed that the eigenvalues decay of the POD matrix is the expected one and that, exploiting the continuation and the deflation method, it is possible to construct a bifurcation diagram in the offline phase and to reconstruct it accurately and efficiently during the online one.
However, the diagram discretized during the online phase loses part of its importance because the same branches have to be discretized also in the offline one. Otherwise, if one wants to obtain it only in the online one without exploiting the continuation and the deflation offline, the reduced space would be too small and the secondary branches could not be found online. For instance, let us assume that during the offline phase one computes all the symmetric solutions, i.e. the ones such that , then the reduced space is able to generate only symmetric solutions and, therefore, the other branches can not be discovered.
Furthermore, if one is interested in computing bifurcation diagrams with more than one parameter, the number of solutions required to discretize it would significantly increase and the computational cost of the high-order simulation would become prohibitive. We remark that the coupling of the aforementioned techniques allows us to efficiently reconstruct a representation of the -dimensional manifold induced by the solutions, that would be infeasible without reduction strategies. Therefore, we decided to slightly change the described approach, computing only two one-dimensional bifurcation diagrams during the offline phase, and refining the grid associated to the second parameter only in the online phase. Such an approach can be generalized, when even more parameters are involved, by computing a small set of one-dimensional diagrams during the offline phase and generating all the other dimensions of the complete diagram only online.
As previously written, the second parameter in this work is a scaling of the inlet Dirichlet boundary condition. The bifurcation diagram that can be obtained with the described approach is shown in figure 10. It can be observed that it is possible to entirely reconstruct the two-dimensionality of the diagram with both the bifurcation points, whose position changes according to the value of . The obtained diagram is coherent with the one in figures 6 and 7 and with the nature of the two involved parameters. In fact the inlet boundary conditions is strongly related to the Reynolds number () and, since such a non-dimensional quantity is the one responsible for the bifurcation, one could expect that the bifurcation points are moved along straight lines in the -direction because multiplying by a scalar is equivalent, in terms of , to dividing by . Moreover, this property explains the fact that, when increases from 0.8 to 1, the bifurcation phenomena are anticipated, i.e. the bifurcation points are associated to higher viscosities. Figure 10. Bifurcation diagram efficiently obtained in the online phase with two parameters and two bifurcations. The colour gradient remarks the value f(\uu) in each point in order to allow the reader to more easily interpret the diagram

We wanted to consider a second parameter that could add new information to the one-dimensional diagram but, since these two parameters are so deeply related to , we decided to analyze the POD eigenvalues decay to observe if a behaviour different than the one observed in the one-dimensional case could be found. A significantly stronger decay would have mean that and contained the same information and, therefore, the same solutions that can be obtained varying the viscosity could be obtained properly rescaling the ones obtained varying .
We analyzed three different scenarios. Firstly, we considered snapshots belonging to a diagram with 2 values of and 15 of , that is very similar to the actual offline phase and it is reported in figure svdecay5. In such a setting, the observed decay follows the behaviour of the reference one-dimensional setting. Secondly, in figure svdecay6 we used only 1 value of and 30 of obtaining 30 solutions on the same branch. Consequently, this decay agrees to the one in figure svdecay3. Finally, in figure svdecay7 we decided to use a mixed approach considering 5 values of and 6 of , the obtained decay is slightly faster than the first one, because of the relation between the parameters, but it is still exponential.