LINEAR Stability ANALYSIS of the wake of a cylinder with STABFEM

This scripts demonstrates how to use StabFem to study the instability in the wake of a 2D cylinder in incompressible flow using LINEAR STABILITY ANALYSIS (LSA)

The script reproduces all figures of section 3 in Fabre et al (ASME-AMR, 2019)

The script performs the following calculations :

  1. Generation of an adapted mesh
  2. Base-flow properties for Re = [2-40]
  3. Stability curves St(Re) and sigma(Re) for Re = [40-100]

The raw source code is here.

Note : the present version of this script uses the current StabFem syntax which has been improved since the publication of the paper. You can see the original script used in the paper here but the syntax employed there may not work any more.

( Updated july 2021, StabFem 3.8 )

Contents

First we set a few global variables needed by StabFem

close all;
addpath('../../SOURCES_MATLAB');
SF_Start('verbosity',2);
SF_DataBase('create','./WORK/');
NOTICE  - Initialization already done

Explanation of these starting options : 'verbosity' controls the ammount of messages produced by StabFem. Recommended values are : '2' -> only important messages, no FreeFem output (recommended for published scripts such as this one) '4' -> more messages plus FreeFem output (recommended for ongoing work) SF_DataBase creates the DataBase directory where result files are stored

##### CHAPTER 1 : COMPUTING THE MESH WITH ADAPTMESH PROCEDURE

We create the initial mesh using the driver SF_Mesh.m which launches the program Lshape_Mesh.edp and imports the generated mesh.

ffmesh = SF_Mesh('Mesh_Cylinder.edp','Options',{'Xmin',-40,'Xmax',80,'Ymax',40});

We initially compute a base flow for a low value of the Reynolds number.

bf=SF_BaseFlow(ffmesh,'solver','Newton_2D_simple.edp','Re',1);
WARNING -  Folder ./WORK/MESHBF not found : creating it !

Here SF_BaseFlow.m is the generic driver to be used to interface FreeFem++ solvers computing a "base flow" (i.e. a steady solution of Navier-Stokes equations) using Newton iteration.

In the present example, we use the solver Newton_2D_simple.edp which is a simplified version, made as short as possible for pedagogical reasons, but which should work only for the present case of a cylinder. In this program, the only argument to be provided through the 'Options' argument is the Reynolds number.

Instead, you may use the generic solver Newton_2D.edp which recognizes a much larger number of options and can be customized to one's needs through macros (see StabFem manual).

Note that the solver is specified here by options "'solver','Newton_2D_simple.edp'" when first calling the driver. This solver will be used by default for all subsequent calls to SF_BaseFlow and so it will not be useful to specify it in the sequel.

We now do a first mesh adaptation.

bf=SF_Adapt(bf,'Hmax',5);

Here we use the SF_Adapt.m function to adapt the mesh to the base flow. This function (which is a wrapper for the FreeFem++ program AdaptMesh.edp) allows to perform mesh adaptation using up to 8 fields with various data type (scalar P1 or P2 data ; composite data for instance [P2xP2xP1] ,...) Have a look at these programs to understand how it works !

We then progressively increase the Reynolds number up to 60, and readapt.

bf=SF_BaseFlow(bf,'Re',10);
bf=SF_BaseFlow(bf,'Re',60);
bf=SF_Adapt(bf,'Hmax',5);
bf=SF_BaseFlow(bf,'Re',60);

Mind that we need to recompute the baseflow after adaptation.

We now compute an eigenvalue/eigenmode pair using shift-invert method (the "shift" is a guess for the eigenvalue ; the performance of the algorithm strongly relies on the appropriate choice of this parameter !) (Note that here we compute the adjoint eigenmode which will be used for mesh adaptation ; see paper AMR, 2018 to learn more on mesh adaptation strategies)

[ev,em] = SF_Stability(bf,'solver','Stab_2D.edp','shift',0.04+0.76i,'nev',1,'Options',{'Symmetry','A','type','A'})
ev =

   0.0471 - 0.7447i


em = 

  struct with fields:

                mesh: [1x1 struct]
            filename: './WORK/EIGENMODES//FFDATA_000001.txt'
     DataDescription: 'Eigenmode for a 2D-incompressible problem '
            datatype: 'EigenModeA'
     datastoragemode: 'CxP2P2P1'
     datadescriptors: 'ux,uy,p'
        meshfilename: './WORK/MESHES/FFMESH_000003.msh'
    baseflowfilename: './WORK/BASEFLOWS/FFDATA_000004.txt'
              lambda: 0.0471 - 0.7447i
            INDEXING: [1x1 struct]
                 sym: -1
               shift: 0.0400 + 0.7600i
                iter: 1
                vort: [1540x1 double]
                  Fy: 0.5000 + 0.0000i
                  Fx: 0
             AEnergy: 0.7448
                  ux: [17520x1 double]
                  uy: [17520x1 double]
                   p: [8760x1 double]
                type: 'D'

Here SF_Stability.m is the generic driver to be used for FreeFem++ linear stability programs compute either a single eigenvalue (with shift-invert with nev=1) or a collection of eigenvalues (with Arnoldi/SLEPC with nev>1). In the present case (2D incompressible problem), the relevant program is Stab_2D.edp.

As with SF_BaseFlow, the solver specified at first call of the driver will be used by default in the sequel so it is not necessary to specify it again.

NB 'Options' is a cell-array list of options to be passed to the solvers and detected by getARGV. In the present solver we have:

1/ 'Symmetry','A' to specify that we are looking for antisymmetric modes (this is the default value and we will not specify it in the next calls)

2/ 'Type','A' to specify that we want to compute the Adjoint modes (default value is 'D' for direct modes)

We then readapt the mesh to both the base flow and the computed eigenmode.

bf = SF_Adapt(bf,em,'Hmax',5);
bf = SF_BaseFlow(bf); % must recompute after adapt

We do this twice, and at second time we also use a 'MASK' which will force the grid size to be at most 0.15 in the wake region defined by -2<x<5 and 0<y<2.

[ev,em] = SF_Stability(bf,'shift',0.04+0.76i,'nev',1,'Options',{'type','A'});
MASK = SF_Mask(bf,[-2 5 0 2 0.15]);
bf=SF_Adapt(bf,em,MASK,'Hmax',5);
bf = SF_BaseFlow(bf);

CHAPTER 1b : DRAW FIGURES for mesh and base flow

figureformat='png';
system('mkdir FIGURES');

%  plot the mesh (full size)
figure();SF_Plot(bf,'mesh','xlim',[-40 80],'ylim',[0 40]);
title('Initial mesh (full size)');
box on; set(gca,'FontSize', 18); saveas(gca,'FIGURES/Cylinder_Mesh_Full',figureformat);

%  plot the mesh (zoom)
figure();SF_Plot(bf,'mesh','xlim',[-1.5 4.5],'ylim',[0 3]);
box on; set(gca,'FontSize', 18); saveas(gca,'FIGURES/Cylinder_Mesh',figureformat);

%   plot the base flow for Re = 60
figure();
SF_Plot(bf,'vort','xlim',[-1.5 4.5],'ylim',[0 3],'cbtitle','\omega_z','colormap','redblue','colorrange',[-2 2]);
ls
box on;  set(gca,'FontSize', 18); saveas(gca,'FIGURES/Cylinder_BaseFlowRe60',figureformat);
mkdir: impossible de créer le répertoire «FIGURES»: Le fichier existe
autorun.m
CYLINDER_LINEAR.m
CYLINDER_LINEAR.oldcache
CYLINDER_LINEAR_SFLAUNCH.m
CYLINDER_MESHTESTS.m
CYLINDER_NONLINEAR.m
DEMO_DataBase_Management.m
DEMO_DataBase_Management.oldcache
DEMO_NEWOPTIONSFORNEWTON_CYLINDER.m
DEMO_NEWOPTIONSFORNEWTON_CYLINDER.oldcache
ErrorCatching.m
FIGURES
freefem++.pref
html
literature_data
Mesh_Cylinder.edp
Newton_2D_simple.edp
Param_Adaptmesh.idp
Param_Mapping.edp
POSTPROCESS_LINEAR_CYLINDER.m
POSTPROCESS_LINEAR_CYLINDER.oldcache
SCRIPT_CYLINDER_ALLFIGURES.m
SCRIPT_CYLINDER_Short.m
SF_AutoInclude.idp
SF_Custom.idp
SF_Geom.edp
SmartMesh_Cylinder.m
WORK

Chapter 2 : eigenmode, adjoint, sensitivity

We now recompute the eigenmode along with the adjoint and structural sensitivity on the new mesh

[ev,em] = SF_Stability(bf,'shift',0.04+0.76i,'nev',1);
[ev,emA] = SF_Stability(bf,'shift',0.04+0.76i,'nev',1,'Options',{'type','A'});
emS = SF_Sensitivity(bf,em,emA,'solver','Sensitivity2D.edp');
Here is how to plot the eigenmode, adjoint and structural sensitivity
figure();SF_Plot(em,'ux','xlim',[-2 8],'ylim',[0 5],'colormap','redblue','colorrange','cropcentered',...
    'boundary','on','bdlabels',2,'bdcolors','k');
box on; set(gca,'FontSize', 18); saveas(gca,'FIGURES/Cylinder_EigenModeRe60_AdaptS',figureformat);  %

figure();SF_Plot(emA,'ux','xlim',[-2 8],'ylim',[0 5],'colormap','redblue','colorrange','cropcentered','boundary','on','bdlabels',2,'bdcolors','k');
box on;  set(gca,'FontSize', 18); saveas(gca,'FIGURES/Cylinder_EigenModeAdjRe60',figureformat);

figure();SF_Plot(emS,'S','xlim',[-2 4],'ylim',[0 3],'colormap','ice','boundary','on','bdlabels',2,'bdcolors','k');
hold on; SF_Plot(bf,'psi','contour','only','clevels', [0 0],'CColor','b','CStyle','monochrome','ColorBar','off','xlim',[-2 4],'ylim',[0 3],...
'colormap','ice','colorrange',[min(real(emS.S)), max(real(emS.S))]);

box on;
pause(0.1);

Chapter 2b : Demonstration of "spectrum explorator" mode

plotoptions = {'ux','xlim',[-2 8],'ylim',[0 5],'colormap','redblue','colorrange','cropcentered',...
                'boundary','on','bdlabels',2,'bdcolors','k'};
ev = SF_Stability(bf,'nev',20,'shift',0.74i,'PlotSpectrum',true,'PlotModeOptions',plotoptions);
figure(100); box on; pos = get(gcf,'Position'); pos(4)=pos(3)*.75;set(gcf,'Position',pos); % resize aspect ratio
set(gca,'FontSize', 14);
saveas(gca,'FIGURES/SpectrumExplorator',figureformat);

CHAPTER 3 : Loop over Re for Base Flow and eigenvalue computation

The objective is to reproduce results in figures 3b and 5 from the AMR 2019 paper : drag of base flow and eigenvalue as function of Reynolds. We will do both base-flow and eigenvalue computations is a single loop and will do the post-processing in a subsequent step using Database explorator.

Notes :

- We will do the eigenvalue computation in continuation mode (option "'shift','cont'), meaning that the "guess" is extrapolated from previous calculations This method is extremely fast but requires an initialization with a first shift before entering the loop.

- The threshold tracker mode is enabled (option 'Threshold','on'); accordingly during the loop if a threshold is encountered (i.e., the real part of lambda crosses zero), its characteristics (here, the corresponding Re) will be directly written in a subfolder 'THRESHOLD' of the database.

- Another possibility is to do the computations in separate loops and generate arrays of data to be plotted within the loops. This method was the one initially used in the AMR 2019 paper, see original script SCRIPT_CYLINDER_ALLFIGURES.m we recommend to use the present method based on data base instead.

Initialisation of loop

Re_BF = [40 : 2.5: 80]; % range of Reynolds to explore.
bf = SF_BaseFlow(bf,'Re',Re_BF(1)); % first
[ev,em] = SF_Stability(bf,'shift',-.03+.72i,'nev',1); % initial shift has been determinated previously

The loop

    for Re = Re_BF
        bf = SF_BaseFlow(bf,'Re',Re);
        ev = SF_Stability(bf,'nev',1,'shift','cont','threshold','on');
    end

After the loop we generate a summary of the database:

sfs = SF_Status
################################################################################
...  SUMMARY OF YOUR DATABASE FOLDER :    ./WORK/
#################################################################################
 
.... CONTENT OF DIRECTORY ./WORK/MESHES :
     (list of meshes previously created/adapted ; couples of .msh/.ff2m files )
 
Index | Name              | generation mode | Date                 | Nv       | H         
1     | FFMESH_000001.msh | initial         | 08-sept.-2021 01:45:57 | 2192     | 40        
2     | FFMESH_000002.msh | adapted         | 08-sept.-2021 01:46:04 | 825      | 40        
3     | FFMESH_000003.msh | adapted         | 08-sept.-2021 01:46:10 | 1540     | 40        
4     | FFMESH_000004.msh | adapted         | 08-sept.-2021 01:46:18 | 3480     | 40        
5     | FFMESH_000005.msh | adapted         | 08-sept.-2021 01:46:33 | 4550     | 40        
 
#################################################################################
#################################################################################
 
.... CONTENT OF DIRECTORY ./WORK/BASEFLOWS
     (couples of .txt/.ff2m files )
 
 Index | Name              | Type         | Date                 | Mesh file         | Re         | Fx        
 1     | FFDATA_000001.txt | BaseFlow     | 08-sept.-2021 01:46:03 | FFMESH_000001.msh | 1          | 5.5507    
 2     | FFDATA_000002.txt | BaseFlow     | 08-sept.-2021 01:46:07 | FFMESH_000002.msh | 10         | 1.3936    
 3     | FFDATA_000003.txt | BaseFlow     | 08-sept.-2021 01:46:10 | FFMESH_000002.msh | 60         | 0.64084   
 4     | FFDATA_000004.txt | BaseFlow     | 08-sept.-2021 01:46:14 | FFMESH_000003.msh | 60         | 0.64317   
 5     | FFDATA_000005.txt | BaseFlow     | 08-sept.-2021 01:46:25 | FFMESH_000004.msh | 60         | 0.64396   
 6     | FFDATA_000006.txt | BaseFlow     | 08-sept.-2021 01:46:42 | FFMESH_000005.msh | 60         | 0.6441    
 7     | FFDATA_000007.txt | BaseFlow     | 08-sept.-2021 01:47:35 | FFMESH_000005.msh | 40         | 0.75419   
 8     | FFDATA_000008.txt | BaseFlow     | 08-sept.-2021 01:47:44 | FFMESH_000005.msh | 40         | 0.75419   
 9     | FFDATA_000009.txt | BaseFlow     | 08-sept.-2021 01:47:56 | FFMESH_000005.msh | 42.5       | 0.73625   
 10    | FFDATA_000010.txt | BaseFlow     | 08-sept.-2021 01:48:08 | FFMESH_000005.msh | 45         | 0.71985   
 11    | FFDATA_000011.txt | BaseFlow     | 08-sept.-2021 01:48:20 | FFMESH_000005.msh | 47.5       | 0.70477   
 12    | FFDATA_000012.txt | BaseFlow     | 08-sept.-2021 01:48:32 | FFMESH_000005.msh | 50         | 0.69084   
 13    | FFDATA_000013.txt | BaseFlow     | 08-sept.-2021 01:48:44 | FFMESH_000005.msh | 52.5       | 0.67792   
 14    | FFDATA_000014.txt | BaseFlow     | 08-sept.-2021 01:48:56 | FFMESH_000005.msh | 55         | 0.66589   
 15    | FFDATA_000015.txt | BaseFlow     | 08-sept.-2021 01:49:08 | FFMESH_000005.msh | 57.5       | 0.65464   
 16    | FFDATA_000016.txt | BaseFlow     | 08-sept.-2021 01:49:20 | FFMESH_000005.msh | 60         | 0.6441    
 17    | FFDATA_000017.txt | BaseFlow     | 08-sept.-2021 01:49:32 | FFMESH_000005.msh | 62.5       | 0.63419   
 18    | FFDATA_000018.txt | BaseFlow     | 08-sept.-2021 01:49:44 | FFMESH_000005.msh | 65         | 0.62485   
 19    | FFDATA_000019.txt | BaseFlow     | 08-sept.-2021 01:49:57 | FFMESH_000005.msh | 67.5       | 0.61603   
 20    | FFDATA_000020.txt | BaseFlow     | 08-sept.-2021 01:50:09 | FFMESH_000005.msh | 70         | 0.60767   
 21    | FFDATA_000021.txt | BaseFlow     | 08-sept.-2021 01:50:21 | FFMESH_000005.msh | 72.5       | 0.59974   
 22    | FFDATA_000022.txt | BaseFlow     | 08-sept.-2021 01:50:33 | FFMESH_000005.msh | 75         | 0.5922    
 23    | FFDATA_000023.txt | BaseFlow     | 08-sept.-2021 01:50:45 | FFMESH_000005.msh | 77.5       | 0.58502   
 24    | FFDATA_000024.txt | BaseFlow     | 08-sept.-2021 01:50:58 | FFMESH_000005.msh | 80         | 0.57816   
 
#################################################################################
 
.... CONTENT OF DIRECTORY ./WORK/EIGENMODES
     (couples of .txt/.ff2m files )
 
 Index | Name              | Type         | Date                 | BaseFlow file     | lambda           | sym       
 1     | FFDATA_000001.txt | EigenModeA   | 08-sept.-2021 01:46:16 | ./WORK/BASEFLOWS/FFDATA_000004.txt | 0.047057-0.74471i | -1        
 2     | FFDATA_000002.txt | EigenModeA   | 08-sept.-2021 01:46:28 | ./WORK/BASEFLOWS/FFDATA_000005.txt | 0.047093-0.74496i | -1        
 3     | FFDATA_000003.txt | EigenMode    | 08-sept.-2021 01:46:51 | ./WORK/BASEFLOWS/FFDATA_000006.txt | 0.047102+0.74488i | -1        
 4     | FFDATA_000004.txt | EigenModeA   | 08-sept.-2021 01:46:56 | ./WORK/BASEFLOWS/FFDATA_000006.txt | 0.047102-0.74488i | -1        
 5     | FFDATA_000005.txt | EigenMode    | 08-sept.-2021 01:47:17 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.047013+0.74317i | -1        
 6     | FFDATA_000006.txt | EigenMode    | 08-sept.-2021 01:47:18 | ./WORK/BASEFLOWS/FFDATA_000006.txt | 0.047102+0.74488i | -1        
 7     | FFDATA_000007.txt | EigenMode    | 08-sept.-2021 01:47:18 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.047217+0.7478i | -1        
 8     | FFDATA_000008.txt | EigenMode    | 08-sept.-2021 01:47:18 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.049829+0.73609i | -1        
 9     | FFDATA_000009.txt | EigenMode    | 08-sept.-2021 01:47:18 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.049031+0.72546i | -1        
 10    | FFDATA_000010.txt | EigenMode    | 08-sept.-2021 01:47:18 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.05433+0.74503i | -1        
 11    | FFDATA_000011.txt | EigenMode    | 08-sept.-2021 01:47:18 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.052625+0.75474i | -1        
 12    | FFDATA_000012.txt | EigenMode    | 08-sept.-2021 01:47:18 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.054777+0.73205i | -1        
 13    | FFDATA_000013.txt | EigenMode    | 08-sept.-2021 01:47:19 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.051487+0.76121i | -1        
 14    | FFDATA_000014.txt | EigenMode    | 08-sept.-2021 01:47:19 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.051291+0.71741i | -1        
 15    | FFDATA_000015.txt | EigenMode    | 08-sept.-2021 01:47:19 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.056443+0.73421i | -1        
 16    | FFDATA_000016.txt | EigenMode    | 08-sept.-2021 01:47:19 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.050788+0.71357i | -1        
 17    | FFDATA_000017.txt | EigenMode    | 08-sept.-2021 01:47:19 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.050714+0.76884i | -1        
 18    | FFDATA_000018.txt | EigenMode    | 08-sept.-2021 01:47:19 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.054759+0.76128i | -1        
 19    | FFDATA_000019.txt | EigenMode    | 08-sept.-2021 01:47:19 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.050221+0.77174i | -1        
 20    | FFDATA_000020.txt | EigenMode    | 08-sept.-2021 01:47:20 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.048911+0.69954i | -1        
 21    | FFDATA_000021.txt | EigenMode    | 08-sept.-2021 01:47:20 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.04566+0.69467i | -1        
 22    | FFDATA_000022.txt | EigenMode    | 08-sept.-2021 01:47:20 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.056405+0.70905i | -1        
 23    | FFDATA_000023.txt | EigenMode    | 08-sept.-2021 01:47:20 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.064952+0.73754i | -1        
 24    | FFDATA_000024.txt | EigenMode    | 08-sept.-2021 01:47:20 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.064872+0.72878i | -1        
 25    | FFDATA_000025.txt | EigenMode    | 08-sept.-2021 01:47:39 | ./WORK/BASEFLOWS/FFDATA_000007.txt | -0.030424+0.72009i | -1        
 
 
 You have previously computed 42 eigenvalues in this session
 The full list will be available as field  'EIGENVALUES' of this function's result
 
 
 In these computations you have identified 1 neutral ceigenvalues
 The full list will be available as field  'THRESHOLD' of this function's result
 
#################################################################################
 
.... CONTENT OF DIRECTORY ./WORK/MESHBF
     (couples of .txt/.ff2m files )
 
 Index | Name              | Type         | Date                 | Mesh file         | Re         | Fx        
 1     | FFDATA_000001.txt | BaseFlow     | 08-sept.-2021 01:46:03 | FFMESH_000001.msh | 1          | 5.5507    
 
#################################################################################
 
.... CONTENT OF DIRECTORY ./WORK/MISC
     (couples of .txt/.ff2m files )
   NB this directory may contain various secondary files produced by StabFem,
      such as flowfields projected after adaptation but not recomputed, adaptation masks, etc... 
 
 Index | Name              | Type         | Date                 | (no bf or mesh indication) | Re         | eigenvalue
 1     | FFDATA_000001.txt | BaseFlow     | 08-sept.-2021 01:46:05 | (unknown)        | 1          | ???       
 2     | FFDATA_000002.txt | BaseFlow     | 08-sept.-2021 01:46:11 | (unknown)        | 60         | ???       
 3     | FFDATA_000003.txt | BaseFlow     | 08-sept.-2021 01:46:20 | (unknown)        | 60         | ???       
 4     | FFDATA_000004.txt | MASK         | 08-sept.-2021 01:46:31 | (unknown)       
 5     | FFDATA_000005.txt | BaseFlow     | 08-sept.-2021 01:46:36 | (unknown)        | 60         | ???       
 6     | FFDATA_000006.txt | Sensitivity  | 08-sept.-2021 01:46:58 | (unknown)        | 60         | 0.047102  
 
#################################################################################

sfs = 

  struct with fields:

         MESHES: [1x5 struct]
      BASEFLOWS: [1x24 struct]
     EIGENMODES: [1x25 struct]
    EIGENVALUES: [1x1 struct]
     THRESHOLDS: [1x1 struct]
         MESHBF: [1x1 struct]
           MISC: [1x6 struct]

CHAPTER 4a : Figures for BASE FLOW PROPERTIES

The object "sfs" has fields 'BASEFLOWS' giving a summary of the corresponding folder. This object is an array of structures with fields Re,Fx in accordance with the metadata of the baseflow declared in the code.

BFS = sfs.BASEFLOWS
BFS = 

  1x24 struct array with fields:

    filename
    Re
    Fx
    meshfilename
    H

We use these arrays to plot Fx as function of Re:

figure(22);hold off;
plot([BFS(8:end).Re],[BFS(8:end).Fx],'b+-','LineWidth',2);
xlabel('Re');ylabel('Fx'); box on; set(gca,'FontSize', 18);
saveas(gca,'FIGURES/Cylinder_Fx_baseflow',figureformat);

Notes - It is mandatory to use square brackets here. The reason is that BFS is a array of structures, so that BFS.Re is a collection of values not an array. [EVS.Re] turns this into a regular array.

- We exclude the baseflows numbered 1 to 7 which were generated before entering the loop

- We could also have done the same thing using object 'sfs.BASEFLOWS' which also has fields 'Re', and 'Fx'. However doing so, data for Re=60 will appear several times due to the initial calculations when adapting the mesh.

Chapter 4B : plot of eigenvalues

Here is how to plot these data, using again postprocessing from data base.

Thanks to metadata inheritance (see manual, sections 4.2.2 and 4.2.3), the object 'EIGENVALUES' has fields 'lambda' (eigenvalue), as well as 'Re' and 'Fx' (metadata inherited from the base flow).

EVS = sfs.EIGENVALUES
EVS = 

  struct with fields:

    lambda: [42x1 double]
       sym: [42x1 double]
     shift: [42x1 double]
     isadj: [42x1 double]
        Re: [42x1 double]
        Fx: [42x1 double]
         H: [42x1 double]

We use these arrays to plot eigenvalue as function of Re. Note that we exclude the 25 first values which were computed before enteering the loop.

Note that contrary to object BFS used previously, here EVS is a structure of arrays not an array of structures so brackets are not mandarory.

Note also that thanks to the threshold tracker mode, the characheristics of the detected threshold are in a field 'THRESHOLDS' of the database. We will plot the corresponding point on the curves with a red circle.

figure(20);
plot(EVS.Re(26:end),real(EVS.lambda(26:end)),'b+-',sfs.THRESHOLDS.Re,0,'ro');
xlabel('Re');ylabel('$\sigma$','Interpreter','latex');
box on; set(gca,'FontSize', 18); saveas(gca,'FIGURES/Cylinder_Sigma_Re',figureformat);

figure(21);hold off;
plot(EVS.Re(26:end),imag(EVS.lambda(26:end))/(2*pi),'b+-',sfs.THRESHOLDS.Re,imag(sfs.THRESHOLDS.lambda)/(2*pi),'ro');
xlabel('Re');ylabel('St');
box on; set(gca,'FontSize', 18); saveas(gca,'FIGURES/Cylinder_Strouhal_Re',figureformat);
pause(0.1);

APPENDIX : HOW TO USE STABFEM WITHOUT MATLAB/OCTAVE

If you don't have (or don't like) the matlab/octave environment, you can perfectly use the FreeFem++ part of the softare directly in a bash terminal ! For instace here is how to perform the mesh adaptation and eigenmode computation (basically equivalent to the first part of the present program) in bash mode.

Note that each time you launch a StabFem script, a file ".stabfem_log.bash" containing all system commands launched by the drivers is produced. It is highly recommended to have a look at this file !

%{
> FreeFem++ -v 0 Mesh_Cylinder.edp -Xmin 40 -Xmax 80 -Ymax 40
> FreeFem++ -v 0 Newton_2D_simple.edp -Re 10
> cp WORK/BaseFlow.txt WORK/BaseFlow_guess.txt
> FreeFem++ -v 0 Newton_2D_simple.edp -Re 60
> cp WORK/BaseFlow.txt WORK/FlowFieldToAdapt1.txt
> FreeFem++ -v 0 ../../SOURCES_FREEFEM/AdaptMesh.edp
   $$ ENTERING ADAPTMESH.edp
   $$ Enter nfields (number of fields to be used for adaptation) >> 1
   $$ Enter storage mode of .txt file number 1 ? (string like ReP2P2P1, etc...) >> ReP2P2P1
   $$ Enter number of additional real scalars associated to flowfield in file number 0 >> 1
   (...)
> cp WORK/FlowFieldToAdapt1.txt WORK/BaseFlow_guess.txt
> cp WORK/mesh_adapt.msh WORK/mesh.msh
> FreeFem++ -v 0 Newton_2D_simple.edp -Re 60
> FreeFem++ -v 0 ../../SOURCES_FREEFEM/Stab2D.edp -shift_r 0.04 -shift_i 0.74 -nev 1
> cp WORK/BaseFlow.txt WORK/FlowFieldToAdapt1.txt
> cp WORK/Sensitivity.txt WORK/FlowFieldToAdapt2.txt
> FreeFem++ -v 0 ../../SOURCES_FREEFEM/AdaptMesh.edp
   $$ Entering ADAPTMESH.edp
   $$ Enter nfields (number of fields to be used for adaptation) >> 2
   $$ Enter storage mode of .txt file number 1 ? (string like ReP2P2P1, etc...) >> ReP2P2P1
   $$ Enter number of additional real scalars associated to flowfield in file number 0 >> ?1
   $$ Enter storage mode of .txt file number 2 ? (string like ReP2P2P1, etc...) >> ReP2
   $$ Enter number of additional real scalars associated to flowfield in file number 1 >> ?0
   (...)
> cp WORK/FlowFieldToAdapt1.txt WORK/BaseFlow_guess.txt
> cp WORK/mesh_adapt.msh WORK/mesh.msh
> FreeFem++ -v 0 Newton_2D_simple.edp -Re 60
> FreeFem++ -v 0 ../../SOURCES_FREEFEM/Stab2D.edp -shift_r 0.04 -shift_i 0.74
%}

%[[PUBLISH]]
%[[WORK]]



%
fclose(fopen('CYLINDER_LINEAR.success','w+'));