# 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).

## 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'


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]);

% We do this twice for security


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 where is the base flow, is the eigenmode, and .

[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.


## 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 16-Jul-2021 00:08:48
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....


## 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 = 100;
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 16-Jul-2021 00:32:56
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 100 -dt 0.0185 -iout 10  -i0 5000 -t0 100
-> Found time series spaning from iter = 5001 to 5100 ; time = 100.019 to 101.85
-> Found 11 Snapshots in this folder. Importing them....


## 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

z =

'/PRODCOM/Singularity/freefem++/v4.7-20201117/install/bin/FreeFem++-mpi - version 4.7 (Tue Nov 17 21:25:09 CET 2020 - git v4.7) 64bits
Usage: /PRODCOM/Singularity/freefem++/v4.7-20201117/install/bin/FreeFem++-mpi [FreeFEM arguments] filename [script arguments]
FreeFEM arguments:
-f:     [filename]  script file name
-v:     [verbosity] level of FreeFEM output (0 - 1000000)
-nw:                no graphics
-wg:                with graphics
-ne:                no edp script output
-cd:                change directory to script directory
-cdtmp:             change directory to /tmp
-jc:                just compile
-ns:                same as -ne
-nowait:            do not wait graphics at the end
-nc:                without console (MS Windows only)
-log:               with console (MS Windows only)
-wait:              wait graphics at the end
-fglut: [filename]  redirect graphics in file
-glut:  [command]   use custom glut
-gff:   [command]   use custom glut (with space quoting)
-check_plugin [filename]        just try if the plugin is correct
-?:                 show help

without default ffglut: ffglut

FreeFEM website: https://freefem.org/
FreeFEM documentation: https://doc.freefem.org/
FreeFEM forum: https://community.freefem.org/
FreeFEM modules: https://modules.freefem.org/

@article{FreeFEM,
AUTHOR = {Hecht, F.},
TITLE = {New development in FreeFem++},
JOURNAL = {J. Numer. Math.},
FJOURNAL = {Journal of Numerical Mathematics},
VOLUME = {20}, YEAR = {2012},
NUMBER = {3-4}, PAGES = {251--265},
ISSN = {1570-2820},
MRCLASS = {65Y15},
MRNUMBER = {3043640},
URL = {https://freefem.org/}
}
'

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         | 16-juil.-2021 00:07:42 | 2192
2     | FFMESH_000002.msh | adapted         | 16-juil.-2021 00:08:02 | 2132
3     | FFMESH_000003.msh | adapted         | 16-juil.-2021 00:08:06 | 3838
4     | FFMESH_000004.msh | mirror          | 16-juil.-2021 00:08:18 | 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     | 16-juil.-2021 00:08:20 | 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     | 16-juil.-2021 00:07:48 | FFMESH_000001.msh | 1
2     | FFDATA_000002.txt | BaseFlow     | 16-juil.-2021 00:07:53 | FFMESH_000001.msh | 10
3     | FFDATA_000003.txt | BaseFlow     | 16-juil.-2021 00:07:59 | FFMESH_000001.msh | 60
4     | FFDATA_000004.txt | BaseFlow     | 16-juil.-2021 00:08:40 | 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    | 16-juil.-2021 00:08:47 | ./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     | 16-juil.-2021 00:07:48 | 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,

Index | Name              | Type         | Date                 | Mesh file
1     | FFDATA_000001.txt | MASK         | 16-juil.-2021 00:08:01 | FFMESH_000001.msh
2     | FFDATA_000002.txt | BaseFlow     | 16-juil.-2021 00:08:03 | FFMESH_000002.msh
3     | FFDATA_000003.txt | MASK         | 16-juil.-2021 00:08:05 | FFMESH_000002.msh
4     | FFDATA_000004.txt | BaseFlow     | 16-juil.-2021 00:08:07 | FFMESH_000003.msh
5     | FFDATA_000005.txt | Addition     | 16-juil.-2021 00:08:48 | FFMESH_000004.msh

#################################################################################
.... Checking for simulation results in folder ./CYLINDER/TimeStepping
#################################################################################
.... CONTENT OF  ./CYLINDER/TimeStepping/RUN000001/ :
Simulation completed

Time-stepping simulation started 16-Jul-2021 00:08:48
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     | 16-juil.-2021 00:08:49 | 60         | 0          | 0          | 0.64359    | 0.0025001
Snapshot_00000010.txt    | DNSField     | 16-juil.-2021 00:08:52 | 60         | 0.2        | 10         | 0.64378    | 0.0024967
(skipping intermediate)    |
Snapshot_00004990.txt    | DNSField     | 16-juil.-2021 00:26:04 | 60         | 99.8       | 4990       | 0.69101    | -0.059569
Snapshot_00005000.txt    | DNSField     | 16-juil.-2021 00:26:06 | 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)"
#################################################################################
.... CONTENT OF  ./CYLINDER/TimeStepping/RUN000002/ :
Simulation completed

Time-stepping simulation started 16-Jul-2021 00:32:56
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 100 -dt 0.0185 -iout 10  -i0 5000 -t0 100
-> Found time series spaning from iter = 5001 to 5100 ; time = 100.019 to 101.85
-> Found 11 Snapshots in this folder.
#################################################################################

.... SNAPSHOTS FOUND IN DIRECTORY ./CYLINDER/TimeStepping/RUN000002/
(couples of .txt/.ff2m files )

Total : 11 files. Displaying only first and last ones
Name                     | Type         | Date                 | Re         | t          | it         | Fx         | Fy
Snapshot_00005000.txt    | DNSField     | 16-juil.-2021 00:32:57 | 60         | 100        | 5000       | 0.69117    | -0.059683
Snapshot_00005010.txt    | DNSField     | 16-juil.-2021 00:33:00 | 60         | 100.185    | 5010       | 0.69129    | -0.058293
(skipping intermediate)    |
Snapshot_00005090.txt    | DNSField     | 16-juil.-2021 00:33:16 | 60         | 101.665    | 5090       | 0.691      | -0.0050703
Snapshot_00005100.txt    | DNSField     | 16-juil.-2021 00:33:18 | 60         | 101.85     | 5100       | 0.69093    | 0.0043671

To import one snapshot, use SF_Load. For instance to import last one : "Snap = SF_Load('TimeStepping',2,5100)"
#################################################################################

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
datatype
datastoragemode
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+'));