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 :
- Generation of an adapted mesh
- Base-flow properties for Re = [2-40]
- 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
- ##### CHAPTER 1 : COMPUTING THE MESH WITH ADAPTMESH PROCEDURE
- CHAPTER 1b : DRAW FIGURES for mesh and base flow
- Chapter 2 : eigenmode, adjoint, sensitivity
- Chapter 2b : Demonstration of "spectrum explorator" mode
- CHAPTER 3 : Loop over Re for Base Flow and eigenvalue computation
- CHAPTER 4a : Figures for BASE FLOW PROPERTIES
- Chapter 4B : plot of eigenvalues
- APPENDIX : HOW TO USE STABFEM WITHOUT MATLAB/OCTAVE
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);
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: []
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
ux: []
uy: []
p: []
status: 'unloaded'
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 Newton_2D_simple.edp CYLINDER_LINEAR.m Param_Adaptmesh.idp CYLINDER_LINEAR_SFLAUNCH.m Param_Mapping.edp CYLINDER_MESHTESTS.m POSTPROCESS_LINEAR_CYLINDER.m CYLINDER_NONLINEAR.m SCRIPT_CYLINDER_ALLFIGURES.m DEMO_DataBase_Management.m SCRIPT_CYLINDER_Short.m DEMO_NEWOPTIONSFORNEWTON_CYLINDER.m SF_AutoInclude.idp ErrorCatching.m SF_Custom.idp FIGURES SF_Geom.edp freefem++.pref SmartMesh_Cylinder.m html WORK literature_data WORK_DEBUG Mesh_Cylinder.edp
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 | 23-mars-2023 01:44:46 | 2192 | 40
2 | FFMESH_000002.msh | adapted | 23-mars-2023 01:45:00 | 825 | 40
3 | FFMESH_000003.msh | adapted | 23-mars-2023 01:45:12 | 1540 | 40
4 | FFMESH_000004.msh | adapted | 23-mars-2023 01:45:26 | 3480 | 40
5 | FFMESH_000005.msh | adapted | 23-mars-2023 01:45:55 | 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 | 23-mars-2023 01:44:58 | FFMESH_000001.msh | 1 | 5.5507
2 | FFDATA_000002.txt | BaseFlow | 23-mars-2023 01:45:05 | FFMESH_000002.msh | 10 | 1.3936
3 | FFDATA_000003.txt | BaseFlow | 23-mars-2023 01:45:10 | FFMESH_000002.msh | 60 | 0.64084
4 | FFDATA_000004.txt | BaseFlow | 23-mars-2023 01:45:19 | FFMESH_000003.msh | 60 | 0.64317
5 | FFDATA_000005.txt | BaseFlow | 23-mars-2023 01:45:41 | FFMESH_000004.msh | 60 | 0.64396
6 | FFDATA_000006.txt | BaseFlow | 23-mars-2023 01:46:14 | FFMESH_000005.msh | 60 | 0.6441
7 | FFDATA_000007.txt | BaseFlow | 23-mars-2023 01:48:18 | FFMESH_000005.msh | 40 | 0.75419
8 | FFDATA_000008.txt | BaseFlow | 23-mars-2023 01:48:38 | FFMESH_000005.msh | 40 | 0.75419
9 | FFDATA_000009.txt | BaseFlow | 23-mars-2023 01:49:01 | FFMESH_000005.msh | 42.5 | 0.73625
10 | FFDATA_000010.txt | BaseFlow | 23-mars-2023 01:49:23 | FFMESH_000005.msh | 45 | 0.71985
11 | FFDATA_000011.txt | BaseFlow | 23-mars-2023 01:49:45 | FFMESH_000005.msh | 47.5 | 0.70477
12 | FFDATA_000012.txt | BaseFlow | 23-mars-2023 01:50:09 | FFMESH_000005.msh | 50 | 0.69084
13 | FFDATA_000013.txt | BaseFlow | 23-mars-2023 01:50:32 | FFMESH_000005.msh | 52.5 | 0.67792
14 | FFDATA_000014.txt | BaseFlow | 23-mars-2023 01:50:58 | FFMESH_000005.msh | 55 | 0.66589
15 | FFDATA_000015.txt | BaseFlow | 23-mars-2023 01:51:22 | FFMESH_000005.msh | 57.5 | 0.65464
16 | FFDATA_000016.txt | BaseFlow | 23-mars-2023 01:51:48 | FFMESH_000005.msh | 60 | 0.6441
17 | FFDATA_000017.txt | BaseFlow | 23-mars-2023 01:52:15 | FFMESH_000005.msh | 62.5 | 0.63419
18 | FFDATA_000018.txt | BaseFlow | 23-mars-2023 01:52:39 | FFMESH_000005.msh | 65 | 0.62485
19 | FFDATA_000019.txt | BaseFlow | 23-mars-2023 01:53:03 | FFMESH_000005.msh | 67.5 | 0.61603
20 | FFDATA_000020.txt | BaseFlow | 23-mars-2023 01:53:32 | FFMESH_000005.msh | 70 | 0.60767
21 | FFDATA_000021.txt | BaseFlow | 23-mars-2023 01:53:58 | FFMESH_000005.msh | 72.5 | 0.59974
22 | FFDATA_000022.txt | BaseFlow | 23-mars-2023 01:54:26 | FFMESH_000005.msh | 75 | 0.5922
23 | FFDATA_000023.txt | BaseFlow | 23-mars-2023 01:54:54 | FFMESH_000005.msh | 77.5 | 0.58502
24 | FFDATA_000024.txt | BaseFlow | 23-mars-2023 01:55:21 | 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 | 23-mars-2023 01:45:24 | ./WORK/BASEFLOWS/FFDATA_000004.txt | 0.047057-0.74471i | -1
2 | FFDATA_000002.txt | EigenModeA | 23-mars-2023 01:45:47 | ./WORK/BASEFLOWS/FFDATA_000005.txt | 0.047093-0.74496i | -1
3 | FFDATA_000003.txt | EigenMode | 23-mars-2023 01:46:39 | ./WORK/BASEFLOWS/FFDATA_000006.txt | 0.047102+0.74488i | -1
4 | FFDATA_000004.txt | EigenModeA | 23-mars-2023 01:46:48 | ./WORK/BASEFLOWS/FFDATA_000006.txt | 0.047102-0.74488i | -1
5 | FFDATA_000005.txt | EigenMode | 23-mars-2023 01:47:43 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.047013+0.74317i | -1
6 | FFDATA_000006.txt | EigenMode | 23-mars-2023 01:47:44 | ./WORK/BASEFLOWS/FFDATA_000006.txt | 0.047102+0.74488i | -1
7 | FFDATA_000007.txt | EigenMode | 23-mars-2023 01:47:44 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.047217+0.7478i | -1
8 | FFDATA_000008.txt | EigenMode | 23-mars-2023 01:47:44 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.049829+0.73609i | -1
9 | FFDATA_000009.txt | EigenMode | 23-mars-2023 01:47:44 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.049031+0.72546i | -1
10 | FFDATA_000010.txt | EigenMode | 23-mars-2023 01:47:45 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.05433+0.74503i | -1
11 | FFDATA_000011.txt | EigenMode | 23-mars-2023 01:47:45 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.052625+0.75474i | -1
12 | FFDATA_000012.txt | EigenMode | 23-mars-2023 01:47:45 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.054777+0.73205i | -1
13 | FFDATA_000013.txt | EigenMode | 23-mars-2023 01:47:45 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.051487+0.76121i | -1
14 | FFDATA_000014.txt | EigenMode | 23-mars-2023 01:47:46 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.051291+0.71741i | -1
15 | FFDATA_000015.txt | EigenMode | 23-mars-2023 01:47:46 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.056443+0.73421i | -1
16 | FFDATA_000016.txt | EigenMode | 23-mars-2023 01:47:46 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.050788+0.71357i | -1
17 | FFDATA_000017.txt | EigenMode | 23-mars-2023 01:47:46 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.050714+0.76884i | -1
18 | FFDATA_000018.txt | EigenMode | 23-mars-2023 01:47:46 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.054759+0.76128i | -1
19 | FFDATA_000019.txt | EigenMode | 23-mars-2023 01:47:47 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.050221+0.77174i | -1
20 | FFDATA_000020.txt | EigenMode | 23-mars-2023 01:47:47 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.048911+0.69954i | -1
21 | FFDATA_000021.txt | EigenMode | 23-mars-2023 01:47:47 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.04566+0.69467i | -1
22 | FFDATA_000022.txt | EigenMode | 23-mars-2023 01:47:48 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.056405+0.70905i | -1
23 | FFDATA_000023.txt | EigenMode | 23-mars-2023 01:47:48 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.064952+0.73754i | -1
24 | FFDATA_000024.txt | EigenMode | 23-mars-2023 01:47:48 | ./WORK/BASEFLOWS/FFDATA_000006.txt | -0.064872+0.72878i | -1
25 | FFDATA_000025.txt | EigenMode | 23-mars-2023 01:48:27 | ./WORK/BASEFLOWS/FFDATA_000007.txt | -0.030424+0.72009i | -1
You have previously computed 25 eigenvalues in this session
The full list will be available as field 'EIGENVALUES' 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 | 23-mars-2023 01:44:58 | 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 | 23-mars-2023 01:45:01 | (unknown) | 1 | ???
2 | FFDATA_000002.txt | BaseFlow | 23-mars-2023 01:45:14 | (unknown) | 60 | ???
3 | FFDATA_000003.txt | BaseFlow | 23-mars-2023 01:45:29 | (unknown) | 60 | ???
4 | FFDATA_000004.txt | MASK | 23-mars-2023 01:45:52 | (unknown)
5 | FFDATA_000005.txt | BaseFlow | 23-mars-2023 01:46:00 | (unknown) | 60 | ???
6 | FFDATA_000006.txt | Sensitivity | 23-mars-2023 01:46:53 | (unknown) | 60 | 0.047102
#################################################################################
sfs =
struct with fields:
MESHES: [1x5 struct]
BASEFLOWS: [1x24 struct]
EIGENMODES: [1x25 struct]
EIGENVALUES: [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:
mesh: []
lambda: [25x1 double]
sym: [25x1 double]
shift: [25x1 double]
isadj: [25x1 double]
Re: [25x1 double]
Fx: [25x1 double]
H: [25x1 double]
status: 'loaded'
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);
Unrecognized field name "THRESHOLDS". Error in CYLINDER_LINEAR (line 321) plot(EVS.Re(26:end),real(EVS.lambda(26:end)),'b+-',sfs.THRESHOLDS.Re,0,'ro');
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+'));