Demonstration of how to do a DNS using StabFem

In this tutorial script we show how to perform a DNS of the wake of a cylinder for Re = 60, allowing to observe the development and saturation of the well-known Von Karman vortex street. The DNS solver is based on Uzawa method using Cahouet preconditionner. The FreeFem++ code TimeStepper_2D.edp is adapted from the one used in Jallas, Marquet & Fabre (PRE, 2017).

Contents

Chapter 0 : initialization

addpath('../../SOURCES_MATLAB');
SF_Start('verbosity',4);
if (SF_core_getopt('isoctave'))
  warning(' This script may take a very long time if using octave because reading .txt files and generating movie is very slow')
end



% NB change verbosity to 4 to follow the simulation while it's running
%  or 2 to supress output from freefem++

SF_DataBase('open','./CYLINDER/');

NOTICE  - Initialization already done
WARNING - Detected verbosity = 4 in Publish mode ; not recommended ! switching to 2 instead

We use the same procedure as in CYLINDER_LINEAR.m. First we build an initial mesh (on a half-domain) and progressively increase the Reynolds :

    ffmesh = SF_Mesh('Mesh_Cylinder.edp','Params',[-40 80 40],'cleanworkdir',true);
    bf=SF_BaseFlow(ffmesh,'solver','Newton_2D.edp','Re',1,'Symmetry','S');
    bf=SF_BaseFlow(bf,'Re',10);
    bf=SF_BaseFlow(bf,'Re',60);

WARNING -  Your program uses depreciated syntax 'Params' to pass parameters 
WARNING -  It is advised to switch to new method using 'Options', and modify your Freefem solvers to detect parameters using 'getARGV'  
WARNING -  Folder ./CYLINDER/MESHBF not found : creating it !

Then we readapt the mesh to fit the base flow. We use an "adaptation mask" generated by SF_Mask.m to enforce the grid size to be at least 0.15 in a window -2<x<10; 0<y<2

    Mask = SF_Mask(bf.mesh,[-2 10 0 2 .1]);
    bf = SF_Adapt(bf,Mask,'Hmax',5);

% We do this twice for security
    Mask = SF_Mask(bf.mesh,[-2 10 0 2 .1]);
    bf = SF_Adapt(bf,Mask,'Hmax',5);

If you think this mesh is not refined enough, a this stage you can use SF_Split.m to divide all triangles in 4 !

%bf = SF_Split(bf);

We now have a convenient mesh defined on a half-domain. We then convert it into a full mesh using SF_Mirror.m

bf = SF_Mirror(bf);
bf = SF_BaseFlow(bf,'solver','Newton_2D.edp','Re',60,'Symmetry','N'); % recompute after mirror.

You now have a robust mesh well fitted to DNS. You may also try other adaptation strategies (adaptation to eigenmode, etc..). See other examples in the project !

Note that the initially generated mesh is not perfectly symmetric. In some cases (e.g. if you want to do DNS very close to the instability threshold Re ~ 47) it is desirable to have a perfectly symmetric mesh. A procedure to do that is to first generate a half-mesh, adapt it and "mirror" it. An example on how to do this is available here : SCRIPT_DNS_WITHSYMMETRICMESH.m.m

Chapter 2 : Generation of a starting point for DNS using stability results

We want to initialise the DNS with an initial condition $u = u_b + \epsilon u_1$ where $u_b$ is the base flow, $u_1$ is the eigenmode, and $\epsilon = 0.01$.

[ev,em] = SF_Stability(bf,'solver','Stab_2D.edp','shift',0.04+0.76i,'nev',1); % compute the eigenmode.

startfield = SF_Add({bf,em},'coefs',[1,0.01],'NameParams',['Re'],'Params',[bf.Re]); % creates startfield = bf+0.01*em

WARNING - The number of additional scalars 1 does not match the number of parameters entered by the user 1 The default action is to choose the structure of the first field. Therefore the first field must contain the user input fields. 
WARNING -  in SF_Add : could not process with auxiliary field vort
WARNING -  in SF_Add : could not process with auxiliary field psi

Chapter 3 : Launch a DNS

We do 5000 time steps and produce snapshots each 10 time steps. We use the driver SF_TS_Launch.m which launches the Freefem++ solver TimeStepper_2D.edp

Niter = 5000;   % total number of time steps
iout = 10;   % will generate snapshots each iout time steps
dt = 0.02;   % time step
DNSOptions = {'Re',60,'Niter',Niter,'dt',dt};

folder = SF_TS_Launch('TimeStepper_2D.edp','Init',startfield,'Options',DNSOptions,'wait',true);

The drivers returns "folder", which corresponds to the name of the directory where results are stored. To import the data in this folder we use SF_TS_Check.m as follows.

[DNSstats,DNSfields] = SF_TS_Check(folder); % Import statistics and snapshots (may be very long with octave...)

#################################################################################
.... CONTENT OF  ./CYLINDER/TimeStepping/RUN000001/ : 
 Simulation completed

 Time-stepping simulation started 23-Mar-2023 02:33:25
 Solver                      : TimeStepper_2D.edp
 Mesh file name              : ./CYLINDER/MESHES/FFMESH_000004.msh
 Initial condition file name : ./CYLINDER/MISC/FFDATA_000005.txt
 Argument string             :  -Re 60 -Niter 5000 -dt 0.02  -i0 0 -t0 0
 -> Found time series spaning from iter = 1 to 5000 ; time = 0.02 to 100
 -> Found 501 Snapshots in this folder. Importing them....
 Simulation completed

Chapter 4 : plotting the results and generating a movie

Here is how to generate a movie

h = figure;
mkdir('html');
filename = 'html/SCRIPT_DNS_EXAMPLE_movie.gif';
for i=1:length(DNSfields)
    SF_Plot(DNSfields(i),'vort','xlim',[-2 10],'ylim',[-3 3 ],'colorrange',[-3 3],...
        'title',['t = ',num2str(DNSfields(i).t)],'boundary','on','bdlabels',2);
    pause(0.1);
    frame = getframe(h);
    im = frame2im(frame);
    [imind,cm] = rgb2ind(im,256);
    if (i==1)
       imwrite(imind,cm,filename,'gif', 'Loopcount',inf,'DelayTime',0.1);
       set(gca,'nextplot','replacechildren');
    else
       imwrite(imind,cm,filename,'gif','WriteMode','append','DelayTime',0.1);
    end
end

Warning: Directory already exists.

Here is the movie

We now plot the lift force as function of time

figure(15);
subplot(2,1,1);
plot(DNSstats.t,DNSstats.Fy);title('Lift force as function of time');
xlabel('t');ylabel('Fy');
subplot(2,1,2);
plot(DNSstats.t,DNSstats.Fx);title('Drag force as function of time');
xlabel('t');ylabel('Fx');

Now we plot the pressure and vorticity along the surface

alpha = linspace(-pi,pi,201);
Xsurf = .501*cos(alpha);
Ysurf = .501*sin(alpha);
Psurf = SF_ExtractData(DNSfields(end),'p',Xsurf,Ysurf);
Omegasurf = SF_ExtractData(DNSfields(end),'vort',Xsurf,Ysurf);
figure(16);
subplot(2,1,1);title('Pressure along the surface P(r=a,\theta) at final time step')
plot(alpha,Psurf);
xlabel('\theta');ylabel('P(r=a,\theta)');
subplot(2,1,2);title('Vorticity along the surface \omega(r=a,\theta) at final time step')
plot(alpha,Omegasurf);
xlabel('\theta');ylabel('\omega_z(r=a,\theta)');

Note: To improve runtime build MEX function fftri2gridfast() from fftri2gridfast.c
Note: To improve runtime build MEX function fftri2gridfast() from fftri2gridfast.c

CHAPTER 5 : how to do a DNS starting from a previous snapshot

Suppose we want to restart the previous DNS and do another 100 time steps

dt = 7.4/400;  Niter = 400;
DNSOptions = {'Re',60,'Niter',Niter,'dt',dt,'iout',10};
StartField = SF_Load('TimeStepping',1,5000); % this will be the starting point

folder = SF_TS_Launch('TimeStepper_2D.edp','Init',StartField,'Options',DNSOptions,'wait',true);
[DNSstats,DNSfields] = SF_TS_Check(folder);

%
% In this last computation we have used 400 steps to span a full oscillation cycle, with period T = 7.4
% This will be used to generate a short movie as follows :

h = figure;
filename = 'html/SCRIPT_DNS_EXAMPLE_movie_short.gif';
for i=1:length(DNSfields)
    SF_Plot(DNSfields(i),'vort','xlim',[-2 10],'ylim',[-3 3 ],'colorrange',[-3 3],...
        'title',['t = ',num2str((i-1)*dt)],'boundary','on','bdlabels',2);
    pause(0.1);
    frame = getframe(h);
    im = frame2im(frame);
    [imind,cm] = rgb2ind(im,256);
    if (i==1)
       imwrite(imind,cm,filename,'gif', 'Loopcount',inf,'DelayTime',0.1);
       set(gca,'nextplot','replacechildren');
    else
       imwrite(imind,cm,filename,'gif','WriteMode','append','DelayTime',0.1);
    end
end

#################################################################################
.... CONTENT OF  ./CYLINDER/TimeStepping/RUN000002/ : 
 Simulation completed

 Time-stepping simulation started 23-Mar-2023 03:21:16
 Solver                      : TimeStepper_2D.edp
 Mesh file name              : ./CYLINDER/MESHES/FFMESH_000004.msh
 Initial condition file name : ./CYLINDER/TimeStepping/RUN000001/Snapshot_00005000.txt
 Argument string             :  -Re 60 -Niter 400 -dt 0.0185 -iout 10  -i0 5000 -t0 100
 -> Found time series spaning from iter = 5001 to 5400 ; time = 100.019 to 107.4
 -> Found 41 Snapshots in this folder. Importing them....
 Simulation completed

Chapter 6 : Data-Mining and post-processing

As for baseflows and eigenmodes in global stability analysis, StabFem stores a copy of all data generated by DNS in a subfolder DNSFLOWS of the working directory. To see the content of this directory, use again the SF_Status.m data-mining wizard :

clear all; close all;
SF_Start('verbosity',2);
SF_DataBase('open','./CYLINDER/');

sf = SF_Status('ALL')

 Warning  : SF_core_start was not previously lanched ; launching right now
WARNING - librairy ARPACK NOT available
################################################################################
...  SUMMARY OF YOUR DATABASE FOLDER :    ./CYLINDER/
#################################################################################
 
.... CONTENT OF DIRECTORY ./CYLINDER/MESHES :
     (list of meshes previously created/adapted ; couples of .msh/.ff2m files )
 
Index | Name              | generation mode | Date                 | Nv      
1     | FFMESH_000001.msh | initial         | 23-mars-2023 02:31:18 | 2192    
2     | FFMESH_000002.msh | adapted         | 23-mars-2023 02:32:01 | 2132    
3     | FFMESH_000003.msh | adapted         | 23-mars-2023 02:32:09 | 3838    
4     | FFMESH_000004.msh | mirror          | 23-mars-2023 02:32:19 | 7526    
#################################################################################
 
     (list of base flows associated to newly computed meshes ; couples of .txt/.ff2m files )
     NB : these base flows are most likely simply projected and not recomputed, it is not recommended to use them
 
 Index | Name              | Type         | Date                 | Mesh file         | Re        
 1     | FFDATA_000009.txt | BaseFlow     | 23-mars-2023 02:32:23 | FFMESH_000004.msh | 60        
 
#################################################################################
#################################################################################
 
.... CONTENT OF DIRECTORY ./CYLINDER/BASEFLOWS
     (couples of .txt/.ff2m files )
 
 Index | Name              | Type         | Date                 | Mesh file         | Re        
 1     | FFDATA_000001.txt | BaseFlow     | 23-mars-2023 02:31:32 | FFMESH_000001.msh | 1         
 2     | FFDATA_000002.txt | BaseFlow     | 23-mars-2023 02:31:43 | FFMESH_000001.msh | 10        
 3     | FFDATA_000003.txt | BaseFlow     | 23-mars-2023 02:31:54 | FFMESH_000001.msh | 60        
 4     | FFDATA_000004.txt | BaseFlow     | 23-mars-2023 02:33:07 | FFMESH_000004.msh | 60        
 
#################################################################################
 
.... CONTENT OF DIRECTORY ./CYLINDER/EIGENMODES
     (couples of .txt/.ff2m files )
 
 Index | Name              | Type         | Date                 | BaseFlow file     | eigenvalue       | sym       
 1     | FFDATA_000001.txt | EigenMode    | 23-mars-2023 02:33:23 | ./CYLINDER/BASEFLOWS/FFDATA_000004.txt | 0.046916+0.74495i | 0         
 
 
 You have previously computed 1 eigenvalues in this session
 The full list will be available as field  'EIGENVALUES' of this function's result
 
#################################################################################
 
.... CONTENT OF DIRECTORY ./CYLINDER/MESHBF
     (couples of .txt/.ff2m files )
 
 Index | Name              | Type         | Date                 | Mesh file         | Re        
 1     | FFDATA_000001.txt | BaseFlow     | 23-mars-2023 02:31:32 | FFMESH_000001.msh | 1         
 
#################################################################################
 
.... CONTENT OF DIRECTORY ./CYLINDER/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                 | Mesh file        
 1     | FFDATA_000001.txt | MASK         | 23-mars-2023 02:31:59 | FFMESH_000001.msh
 2     | FFDATA_000002.txt | BaseFlow     | 23-mars-2023 02:32:03 | FFMESH_000002.msh
 3     | FFDATA_000003.txt | MASK         | 23-mars-2023 02:32:06 | FFMESH_000002.msh
 4     | FFDATA_000004.txt | BaseFlow     | 23-mars-2023 02:32:12 | FFMESH_000003.msh
 5     | FFDATA_000005.txt | Addition     | 23-mars-2023 02:33:24 | FFMESH_000004.msh
 
#################################################################################
.... Checking for simulation results in folder ./CYLINDER/TimeStepping
#################################################################################
.... CONTENT OF  ./CYLINDER/TimeStepping/RUN000001/ : 
 Simulation completed

 Time-stepping simulation started 23-Mar-2023 02:33:25
 Solver                      : TimeStepper_2D.edp
 Mesh file name              : ./CYLINDER/MESHES/FFMESH_000004.msh
 Initial condition file name : ./CYLINDER/MISC/FFDATA_000005.txt
 Argument string             :  -Re 60 -Niter 5000 -dt 0.02  -i0 0 -t0 0
 -> Found time series spaning from iter = 1 to 5000 ; time = 0.02 to 100
 -> Found 501 Snapshots in this folder.
#################################################################################
 
.... SNAPSHOTS FOUND IN DIRECTORY ./CYLINDER/TimeStepping/RUN000001/
     (couples of .txt/.ff2m files )
 
Total : 501 files. Displaying only first and last ones
Name                     | Type         | Date                 | Re         | t          | it         | Fx         | Fy        
Snapshot_00000000.txt    | DNSField     | 23-mars-2023 02:33:27 | 60         | 0          | 0          | 0.64359    | 0.0025001 
Snapshot_00000010.txt    | DNSField     | 23-mars-2023 02:33:34 | 60         | 0.2        | 10         | 0.64378    | 0.0024967 
(skipping intermediate)    |
Snapshot_00004990.txt    | DNSField     | 23-mars-2023 03:11:18 | 60         | 99.8       | 4990       | 0.69101    | -0.059569 
Snapshot_00005000.txt    | DNSField     | 23-mars-2023 03:11:21 | 60         | 100        | 5000       | 0.69117    | -0.059683 
 
 To import one snapshot, use SF_Load. For instance to import last one : "Snap = SF_Load('TimeStepping',1,5000)"
 Simulation completed
 
 
#################################################################################
.... CONTENT OF  ./CYLINDER/TimeStepping/RUN000002/ : 
 Simulation completed

 Time-stepping simulation started 23-Mar-2023 03:21:16
 Solver                      : TimeStepper_2D.edp
 Mesh file name              : ./CYLINDER/MESHES/FFMESH_000004.msh
 Initial condition file name : ./CYLINDER/TimeStepping/RUN000001/Snapshot_00005000.txt
 Argument string             :  -Re 60 -Niter 400 -dt 0.0185 -iout 10  -i0 5000 -t0 100
 -> Found time series spaning from iter = 5001 to 5400 ; time = 100.019 to 107.4
 -> Found 41 Snapshots in this folder.
#################################################################################
 
.... SNAPSHOTS FOUND IN DIRECTORY ./CYLINDER/TimeStepping/RUN000002/
     (couples of .txt/.ff2m files )
 
Total : 41 files. Displaying only first and last ones
Name                     | Type         | Date                 | Re         | t          | it         | Fx         | Fy        
Snapshot_00005000.txt    | DNSField     | 23-mars-2023 03:21:18 | 60         | 100        | 5000       | 0.69117    | -0.059683 
Snapshot_00005010.txt    | DNSField     | 23-mars-2023 03:21:24 | 60         | 100.185    | 5010       | 0.69129    | -0.058293 
(skipping intermediate)    |
Snapshot_00005390.txt    | DNSField     | 23-mars-2023 03:24:27 | 60         | 107.215    | 5390       | 0.69248    | -0.061114 
Snapshot_00005400.txt    | DNSField     | 23-mars-2023 03:24:31 | 60         | 107.4      | 5400       | 0.69262    | -0.061285 
 
 To import one snapshot, use SF_Load. For instance to import last one : "Snap = SF_Load('TimeStepping',2,5400)"
 Simulation completed
 
 
 
#################################################################################

sf = 

  struct with fields:

          MESHES: [1x4 struct]
       BASEFLOWS: [1x4 struct]
      EIGENMODES: [1x1 struct]
     EIGENVALUES: [1x1 struct]
          MESHBF: [1x1 struct]
            MISC: [1x5 struct]
    TimeStepping: [1x2 struct]

We see here that the software has stored 500 snapshots in folder "RUN00001' and 10 snapshots in folder 'RUN00002'. They can be imported with SF_Load. For instance the snapshot corresponding to iteration = 40000, of run 1 :

Snap = SF_Load('TimeStepping',1,4000);
figure;
SF_Plot(Snap,'vort','xlim',[-2 10],'ylim',[-3 3 ],'colorrange',[-3 3],...
       'title','DNS Snapshot #40','boundary','on','bdlabels',2 );

The statistics are available as field 'TimeStepping' of the object returned by SF_Stat. Let us see the contents of this object:

sf.TimeStepping

ans = 

  1x2 struct array with fields:

    filename
    DataDescription
    datatype
    datastoragemode
    datadescriptors
    t
    Fx
    Fy
    Energy
    cfl
    iter
    status

If we want to plot the drag/lift phase portrait corresponding to the whole dataset (with dots) and for the last 400 points (with red line)

figure;plot(sf.TimeStepping(1).Fx,sf.TimeStepping(1).Fy,':');hold on;
plot(sf.TimeStepping(2).Fx,sf.TimeStepping(2).Fy,'r-');
xlabel('Fx (Drag)');ylabel('Fy (Lift)');

% [[PUBLISH]] (this tag is to enable automatic publication as html ; don't touch this line)

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


Published with MATLAB® R2022b