## Chapter 3Stability Analysis

After having computed the base flow it is now possible to calculate the eigenvalues and the eigenmodes of the linearised Navier-Stokes equations. Two different algorithms can be used to solve the equations:

• the Velocity Correction Scheme (`VelocityCorrectionScheme`) and
• the Coupled Linearised Navier-Stokes algorithm (`CoupledLinearisedNS`).

We will consider both cases, highlighting the similarities and differences of these two methods. In this tutorial we will use the Implicitly Restarted Arnoldi Method (IRAM), which is implemented in the open-source library ARPACK and the modified Arnoldi algorithm1 that is also available in Nektar++.

### 3.1 Velocity Correction Scheme

First, we will compute the leading eigenvalues and eigenvectors using the velocity correction scheme method. In the `\$NEKTUTORIAL/stability` folder there is a file called `Channel-VCS.xml`. This file is similar to `Channel-Base.xml`, but contains additional instructions to perform the direct stability analysis.

Note: The entire `GEOMETRY` section, and `EXPANSIONS` section must be identical to that used to compute the base flow.

Task: 3.1 Configure the following additional `SOLVERINFO` options which are related to the stability analysis.

1. set the `EvolutionOperator` to `Direct` in order to activate the forward linearised Navier-Stokes system.
2. set the `Driver` to `Arpack` in order to use the ARPACK eigenvalue analysis.
3. Instruct ARPACK to converge onto specific eigenvalues through the solver property `ArpackProblemType`. In particular, set `ArpackProblemType` to `LargestMag` to get the eigenvalues with the largest magnitude (that determines the stability of the flow).

Note: It is also possible to select the eigenvalue with the largest real part by setting `ArpackProblemType` to (`LargestReal)` or with the largest imaginary part by setting `ArpackProblemType` to (`LargestImag`).

Task: 3.2 Set the parameters for the IRAM algorithm.

• `kdim=16`: dimension of Krylov-space,
• `nvec=2`: number of requested eigenvalues,
• `nits=500`: number of maximum allowed iterations,
• `evtol=1e-6`: accepted tolerance on the eigenvalues and it determines the stopping criterion of the method.

Task: 3.3 Configure the two `FUNCTION` called `InitialConditions` and `BaseFlow`.

1. A restart file is provided to accelerate communications. Set the `InitialConditions` function to be read from `Channel-VCS.rst`. The solution will then converge after 16 iterations after it has populated the Krylov subspace.

Note: The restart file is a field file (same format as `.fld` files) that contains the eigenmode of the system.

Note: Since the simulations often take hundreds of iterations to converge, we will not initialise the IRAM method with a random vector during this tutorial. Normally, a random vector would be used by setting the SolverInfo option `InitialVector` to `Random`.

2. The base flow file (`Channel-Base.fld`), computed in the previous section, should be copied into the `Channel/Stability` folder and renamed `Channel-VCS.bse`. Now specify a function called `BaseFlow` which reads this file.

Task: 3.4 Run the solver to perform the analysis

`\$NEK/IncNavierStokesSolver Channel-VCS.xml`

At the end of the simulation, the terminal screen should look like this:

Iteration 16, output: 0, ido=99
Converged in 16 iterations
Converged Eigenvalues: 2
Magnitude   Angle       Growth      Frequency
EV: 0 1.00112     0.124946    0.0022353   0.249892
Writing: "Channel-al_eig_0.fld"
EV: 1 1.00112     -0.124946   0.0022353   -0.249892
Writing: "Channel-al_eig_1.fld"
L 2 error (variable u) : 0.0367941
L inf error (variable u) : 0.0678149
L 2 error (variable v) : 0.0276887
L inf error (variable v) : 0.0649249
L 2 error (variable p) : 0.00512347
L inf error (variable p) : 0.00135455

The eigenvalues are computed in the exponential form Me where M = |λ| is the magnitude, while θ = arctan(λi∕λr) is the phase:

 (3.1)

It is interesting to consider more general quantities that do not depend on the time length chosen for each iteration T. For this purpose we consider the growth rate σ = ln(M)∕T and the frequency ω = θ∕T.

Figure 3.1 shows the profile of the computed eigenmode. The eigenmodes associated with the computed eigenvalues are stored in the files `Channel_VCS_eig_0.fld` and `Channel_VCS_eig_1.fld`. It is possible to convert this file into VTK format in the same way as previously done for the base flow.

Task: 3.5 Verify that for the channel flow case :

 σ = 2.2353 × 10-3 ω = ±2.49892 × 10-1

and that the eigenmodes match those given in figures 3.1.

This values are in accordance with the literature, in fact in Canuto et al., 1988 suggests 2.23497 × 10-3 and 2.4989154 × 10-1 for growth and frequency, respectively.

Tip: Note that Nektar++ implements also the modified Arnoldi algorithm. You can try to use it for this test case by setting `Driver` to `ModifiedArnoldi`. You can now try to re-run the simulation and verify that the modified Arnoldi algorithm provides a results that is consistent with the previous computation obtained with Arpack.

### 3.2 Coupled Linearised Navier-Stokes algorithm

Note: Remember to use the files provided in the folder `Stability/Coupled` for this case.

It is possible to perform the same stability analysis using a different method based on the Coupled Linearised Navier-Stokes algorithm. This method requires the solution of the full velocity-pressure system, meaning that the velocity matrix system and the pressure system are coupled, in contrast to the velocity correction scheme/splitting schemes.

Inside the folder `\$NEKTUTORIAL/stability` there is a file called `Channel-Coupled.xml` that contains all the necessary parameters that should be defined. In this case we will specify the base flow through an analytical expression. Even in this case, the geometry, the type and number of modes are the the same of the previous simulations.

Task: 3.6 Edit the file `Channel-Coupled.xml`:

Note: As before the bits to be completed are identified by …in this file.

• Set the `SolverType` property to `CoupledLinearisedNS` in order to solve the linearised Navier-Stokes equations using Nektar + +’s coupled solver.
• the `EQTYPE` must be set to `SteadyLinearisedNS` and the `Driver` to `Arpack`.
• Set the `InitialVector` property to `Random` to initialise the IRAM with a random initial vector. In this case the function `InitialConditions` will be ignored.
• To compute the eigenvalues with the largest magnitude, specify `LargestMag` in the property `ArpackProblemType`.

It is important to note that the use of the coupled solver requires that only the velocity component variables are specified, while the pressure is implicitly evaluated.

Task: 3.7 Continue modifying `Channel-Coupled.xml`:

• It is necessary to set up the base flow. For the `SteadyLinearisedNS` coupled solver, this is defined through a function called `AdvectionVelocity`. The u component must be set up to 1 - y2, while the v-component to zero.

For the coupled solver, it is also necessary to define the following additional tag outside of the `CONDITIONS` tag:

1<FORCING>
2   <FORCE TYPE="StabilityCoupledLNS">
3   </FORCE>
4</FORCING>

This has already been set up in the XML file. This is necessary to tell Nektar++ to use the previous solution as the right hand side vector for each Arnoldi iteration.

Task: 3.8 Now run the solver to compute the eigenvalues

`\$NEK/IncNavierStokesSolver Channel-Coupled.xml`

The terminal screen should look like this:

=======================================================================
Solver Type: Coupled Linearised NS
=======================================================================
Arnoldi solver type    : Arpack
Arpack problem type    : LM
Single Fourier mode    : false
Beta set to Zero       : false
Shift (Real,Imag)      : 0,0
Krylov-space dimension : 64
Number of vectors      : 4
Max iterations         : 500
Eigenvalue tolerance   : 1e-06
======================================================
Initial Conditions:
- Field u: 0 (default)
- Field v: 0 (default)
Matrix Setup Costs: 0.565916
Multilevel condensation: 0.098134
Inital vector       : random
Iteration 0, output: 0, ido=-1
Writing: "Channel-Coupled.fld"
Iteration 20, output: 0, ido=1
Writing: "Channel-Coupled.fld"
Iteration 40, output: 0, ido=1
Writing: "Channel-Coupled.fld"
Iteration 60, output: 0, ido=1
Writing: "Channel-Coupled.fld"
Iteration 65, output: 0, ido=99

Converged in 65 iterations
Converged Eigenvalues: 4
Real        Imaginary
EV:  0  -0.000328987            -0
Writing: "Channel-Coupled_eig_0.fld"
EV:  1   -0.00131595            -0
Writing: "Channel-Coupled_eig_1.fld"
EV:  2   -0.00296088            -0
Writing: "Channel-Coupled_eig_2.fld"
EV:  3   -0.00526379            -0
Writing: "Channel-Coupled_eig_3.fld"
L 2 error (variable u) : 2.58891
L inf error (variable u) : 1.00401
L 2 error (variable v) : 0.00276107
L inf error (variable v) : 0.0033678

Using the Stokes algorithm, we are computing the leading eigenvalue of the inverse of the operator L-1. Therefore the eigenvalues of L are the inverse of the computed values2 . However, it is interesting to note that these values are different from those calculated with the Velocity Correction Scheme, producing an apparent inconsistency. However, this can be explained considering that the largest eigenvalues associated to the operator L correspond to the ones that are clustered near the origin of the complex plane if we consider the spectrum of L-1. Therefore, eigenvalues with a smaller magnitude may be present but are not associated with the largest-magnitude eigenvalue of operator L. One solution is to consider a large Krylov dimension specified by kdim and the number of eigenvalues to test using nvec. This will however take more iterations. Another alternative is to use shifting but in this case it will make a real problem into a complex one (we shall show an example later). Finally, another alternative is to search for the eigenvalue with a different criterion, for example, the largest imaginary part.

Task: 3.9 Set up the Solver Info tag `ArpackProblemType` to `LargestImag` and run the simulation again.

=======================================================================
Solver Type: Coupled Linearised NS
=======================================================================
Arnoldi solver type    : Arpack
Arpack problem type    : LI
Single Fourier mode    : false
Beta set to Zero       : false
Shift (Real,Imag)      : 0,0
Krylov-space dimension : 64
Number of vectors      : 4
Max iterations         : 500
Eigenvalue tolerance   : 1e-06
======================================================
Initial Conditions:
- Field u: 0 (default)
- Field v: 0 (default)
Matrix Setup Costs: 0.557085
Multilevel condensation: 0.101482
Inital vector       : random
Iteration 0, output: 0, ido=-1
Writing: "Channel-Coupled.fld"
Iteration 20, output: 0, ido=1
Writing: "Channel-Coupled.fld"
Iteration 40, output: 0, ido=1
Writing: "Channel-Coupled.fld"
Iteration 60, output: 0, ido=1
Writing: "Channel-Coupled.fld"
Iteration 65, output: 0, ido=99

Converged in 65 iterations
Converged Eigenvalues: 4
Real        Imaginary
EV:  0    0.00223509      0.249891
Writing: "Channel-Coupled_eig_0.fld"
EV:  1    0.00223509     -0.249891
Writing: "Channel-Coupled_eig_1.fld"
EV:  2    -0.0542748      0.300562
Writing: "Channel-Coupled_eig_2.fld"
EV:  3    -0.0542748     -0.300562
Writing: "Channel-Coupled_eig_3.fld"
L 2 error (variable u) : 2.58891
L inf error (variable u) : 1.00401
L 2 error (variable v) : 0.00276107
L inf error (variable v) : 0.0033678

In this case, it is easy to to see that the eigenvalues of the evolution operator L are the same ones computed in the previous section with the time-stepping approach (apart from round-off errors). It is interesting to note that this method converges much quicker that the time-stepping algorithm. However, building the coupled matrix that allows us to solve the problem can take a non-negligible computational time for more complex cases.

### 3.3 Three-dimensional Channel flow

Now that we have presented the various stability-analysis tools present in Nektar++, we conclude showing the capabilities of the code in three spatial dimensions. In the folder
`\$NEKTUTORIAL/stability3D` there are the files that are required for the stability analysis - note that we do not show the geometry and the base flow generation (we will be using the exact solution) since we have already presented these features in the previous tutorials.

The case considered is similar to the channel flow presented in section 2. However, in this case the Reynolds number is set to 10000. In order to run a three-dimensional simulation, we can either run the full 3D solver by creating a 3D geometry or use a 2D geometry and specify the use of a Fourier expansion in the third direction. The last method is also known as 3D homogenous 1D approach. Here we will present this approach.

Specifically, we use a 2D geometry and we add the various parameters necessary to use the Fourier expansion. Note that in the 2D plane we will use a `MODIFIED` expansion basis with `NUMMODES=11`.

Task: 3.10 In the file `\$NEKTUTORIAL/stability3D/PPF_R10000_3D.xml`, make the following changes:

• Add a `SOLVERINFO` tag called `HOMOGENEOUS` and set it to `1D`.
• Add two additional `SOLVERINFO` tags called `ModeType` and `BetaZero` and set them to `SingleMode` and `True`, respectively.
• Add two `PARAMETERS` called `HomModesZ` and `LZ` and set them to `2` and `1`, respectively.
• Add two other `PARAMETERS` called `realShift` and `imagShift` and set them to `0.003` and `0.2`, respectively.

Now run the solver - the terminal screen should look like this:

=======================================================================
Solver Type: Coupled Linearised NS
=======================================================================
Arnoldi solver type    : Modified Arnoldi
Single Fourier mode    : true
Beta set to Zero       : true (overrides LHom)
Shift (Real,Imag)      : 0.003,0.2
Krylov-space dimension : 64
Number of vectors      : 2
Max iterations         : 500
Eigenvalue tolerance   : 1e-06
======================================================
Initial Conditions:
- Field u: 0 (default)
- Field v: 0 (default)
- Field w: 0 (default)
Writing: "PPF_R10000_3D_0.chk"
Matrix Setup Costs: 1.97987
Multilevel condensation: 0.427631
Inital vector       : random
Iteration: 0
Iteration: 1 (residual : 4.89954)
Iteration: 2 (residual : 3.64295)
Iteration: 3 (residual : 2.54314)
....
Iteration: 20 (residual : 1.35156e-05)
Iteration: 21 (residual : 1.64786e-06)
Iteration: 22 (residual : 1.92473e-07)
Writing: "PPF_R10000_3D.fld"
L 2 error (variable u) : 3.01846
L inf error (variable u) : 2.25716
L 2 error (variable v) : 1.8469
L inf error (variable v) : 0.985775
L 2 error (variable w) : 5.97653e-06
L inf error (variable w) : 1.2139e-05
EV:  0      0.518448      -26.6405    0.00373022      0.162477
Writing: "PPF_R10000_3D_eig_0.fld"
EV:  1      0.518448       26.6405    0.00373022      0.237523
Writing: "PPF_R10000_3D_eig_1.fld"
Warning: Level 0 assertion violation
Complex Shift applied. Need to implement Ritz re-evaluation of eigenvalue.
Only one half of complex value will be correct

Now convert the two files containing the eigenvectors and visualise them in Paraview or VisIt - the solution should look like the one below:

Task: 3.11 The complete input file `\$NEKTUTORIAL/stability3D/PPF_R15000_3D.xml` has been provided to show a full 3D unstable eigenmode where β is not zero. Run this file and see that you obtain the eigenvalue 0.00248682 ±-0.158348i

Task: 3.12 You can now see what the difference when not using an imaginary shifting. Set the parameters imagShift=0, kdim=384 and nvec=196.

This should take O(500) iterations to complete and hidden in the list of eigenvalues should be the unstable values 0.00248662 ± 0.158347i. They were eigevalues 132 and 133 in my run.

### 3.4 Solutions

Completed solutions to the tutorials are available in the `completed` directory.

This completes the tutorial.