Revision as of 13:06, 7 September 2012 (view source) Bossenne (Talk | contribs) (→‎Hypotheses) ← Older edit		Revision as of 13:07, 7 September 2012 (view source) Bossenne (Talk | contribs) (→‎Data Structures) Newer edit →
Line 91:		Line 91:
	* Optimized non-blocking MPI communications		* Optimized non-blocking MPI communications
−	* Homogeneous mesh refinement	+	* Homogeneous mesh refinement in order to be more precise
	* VTK post-treatment possible		* VTK post-treatment possible

---

Revision as of 13:07, 7 September 2012

Work in progress, come back later ! -- Partie publique

Presentation

H-Allegro is a parallel, high accuracy, fully compressible, reactive, direct Navier Stokes solver. It is the MPI-parallelized evolution of former solver ALLEGRO. H-Allegro has mainly been developed during the PhD work of Eric Albin and Marianne Sjostrand-Cuif, by Yves D'Angelo (PhD Director) and also by the PARALGO company.

Models

The model is defined thanks to a one step Arrhenius' law ^[1] and partially premixed gas for chemistry, the Schmidt number is assumed as constant for the mixing and the viscosity is described by a power law.

Numerics

• Spatial : the case is solved thanks to a 6th-order finite difference scheme and an original staggered scheme developped by E. Albin ^[2].

• Temporal : the resolution is carried out by a 3rd order explicit time integration (Runge-Kutta3) ^[3].

• Boundary conditions : the 3D-NSCBC processing takes into account the inputs and outputs waves.

Equations and Hypotheses

Equations used in H-Allegro

The set of equations solved by H-Allegro are the fully compressible reactive Navier-Stokes equations in conservative form. They write :

• ${\displaystyle {\dfrac {\partial {\rho }}{\partial {t}}}+{\dfrac {\partial {{\rho }U_{j}}}{\partial {x_{j}}}}=0}$

• ${\displaystyle {\dfrac {\partial {{\rho }U_{i}}}{\partial {t}}}+{\dfrac {\partial {\rho }U_{i}U_{j}}{\partial {x_{j}}}}+{\dfrac {\partial {P}}{\partial {x_{i}}}}={\dfrac {\partial {\tau _{ij}}}{\partial {x_{j}}}}}$

• ${\displaystyle {\dfrac {\partial {{\rho }E}}{\partial {t}}}+{\dfrac {\partial {(P+{\rho }E)}U_{j}}{\partial {x_{j}}}}={\dfrac {\partial {q_{j}}}{\partial {x_{j}}}}+{\dfrac {\partial {\tau _{ij}U_{i}}}{\partial {x_{j}}}}+S_{5}}$

• ${\displaystyle {\dfrac {\partial {{\rho }Y_{k}}}{\partial {t}}}+{\dfrac {\partial {{\rho }Y_{k}U_{j}}}{\partial {x_{j}}}}={\dfrac {\partial {q_{j}^{k}}}{\partial {x_{j}}}}+S_{k}}$

${\displaystyle \displaystyle \rho }$ is the density, ${\displaystyle \displaystyle U_{j}}$ is the speed, ${\displaystyle \displaystyle P}$ is the pressure, ${\displaystyle \displaystyle \tau _{ij}}$ is the viscosity tensor, ${\displaystyle \displaystyle E}$ is the mass energy, ${\displaystyle \displaystyle q_{j}}$ is the heat flow by gas thermal conduction, ${\displaystyle S_{5}={\dot {\omega }}_{T}}$ is the heat release due to combustion, ${\displaystyle \displaystyle Y_{k}}$ is the mass fracion of the species k and ${\displaystyle \displaystyle S_{k}}$ is the source term for species k.

${\displaystyle S_{6}={\nu }_{F}M_{F}{\dot {\omega }}}$ is the fuel source term, while ${\displaystyle S_{7}={\nu }_{O}M_{O}{\dot {\omega }}}$ is the oxidant source term.

Closure equations for the code

• ${\displaystyle \rho E=\rho C_{V}T+{\dfrac {1}{2}}\rho U_{i}^{2}}$

• ${\displaystyle \displaystyle P=\rho rT}$

• ${\displaystyle \mu =\mu _{0}\left({\dfrac {T}{T_{0}}}\right)^{\alpha _{s}}}$           ${\displaystyle \displaystyle \alpha _{s}=0.76}$

• ${\displaystyle q_{j}=\lambda {\dfrac {\partial T}{\partial x_{j}}}}$           ${\displaystyle \lambda ={\dfrac {\mu {C_{p}}}{Pr}}}$

• ${\displaystyle q_{j}^{k}=D{\dfrac {\partial Y_{k}}{\partial x_{j}}}}$           ${\displaystyle D={\dfrac {\mu }{Sc_{k}}}}$

• ${\displaystyle \tau _{ij}=\mu \left({\dfrac {\partial U_{i}}{\partial x_{j}}}+{\dfrac {\partial U_{j}}{\partial x_{i}}}\right)-{\dfrac {2}{3}}\mu {\dfrac {\partial U_{k}}{\partial x_{k}}}\delta _{ij}}$

${\displaystyle \displaystyle T}$ is the temperature, ${\displaystyle \displaystyle \mu }$ is the dynamic vicosity, ${\displaystyle \displaystyle \lambda }$ is the thermal conductivity, ${\displaystyle \displaystyle C_{p}}$ is the heat capacity at constant pressure, ${\displaystyle \displaystyle Pr}$ is the Prandtl number, ${\displaystyle \displaystyle D}$ is the diffusion coefficient.

Hypotheses

The hypotheses to simplify the Navier-Stokes equations are :

• The total energy balance according to Poinsot-Veynante is described below, if ${\displaystyle E=e_{s}+{\dfrac {1}{2}}u_{i}^{2}}$.

• The heat flows due to mass fraction gradients and the species molecular distribution due to temperature gradients (Soret and Dufour effects) are neglected ^[4].

• The Laplace coefficient of perfect gas is considered constant :

${\displaystyle \gamma ={\dfrac {C_{P}}{C_{V}}}=1.4}$

The heat capacity at constant pressure and the heat capacity at constant volume are defined as follows :

${\displaystyle C_{P}={\dfrac {\gamma {R}}{\gamma {-1}}}}$ ${\displaystyle C_{P}={\dfrac {R}{\gamma {-1}}}}$

As ${\displaystyle \displaystyle \gamma }$ is constant, they are also constant.

• In order to take the acoustic waves into account, the acoustic Reynolds number is defined such as :

${\displaystyle Re_{ac}={\dfrac {cL}{\nu }}}$

c is the speed of sound :${\displaystyle 340m.s^{-1}}$, L is the characteristic length and ${\displaystyle \nu }$ is the kinematic viscosity of air :${\displaystyle 1.45e^{-5}m.s^{-1}}$.

Data Structures

• 1D, 2D, 3D structured solver

• Optimized non-blocking MPI communications

• Homogeneous mesh refinement in order to be more precise

• VTK post-treatment possible

• Binary files

• Unformatted Fortran90 format data files

Software Engineering

• Fortran90 with modules

• Keyword-based input file

Gallery

H-Allegro

Test image

Performances

Courbe ! (Accélération en WS)

Partie privée -- User's manual

Parameters

Makefile

The makefile allows to compile routines and subroutines. In order to compile, the command "make" will be written in the terminal.

The default target of this file is "all : :".

Since various computers can be used, the compilation has to be conditional. That is to say that several configurations of compilation are available in function of the considered cluster (Antares, Jade, ...) :

• FC indicates which compiler is used. For instance : if the language is fortran90 and MPI, it will be written FC = mpif90.
• A part of the code is written in C, namely the file Vtk_bin_write.c. So the compiler for the C language has to be defined thanks to CFLAGS.
• FFLAGS defines the optimization options for the Fortran compiler. For instance :

FFLAGS = -O3 -ipo -p.

"O3" is an optimization option of level 3. From this level, semantics of the routine can be modified. The optimization will be more aggressive concerning the memory access and the buckles management. The compilation time will be more important than for a lower level and the code size will be also bigger.

"ipo" is an interprocedural optimization option. It authorizes the inlining of present subroutines in seperated source files.

"p" is a profiling option.

The common files of the executable for calculation and for postprocessing have to be defined.

For example :

all : : Mezo3D df6_pospro
COMMON_OBJS += communs.o
COMMON_OBJS += lectus.o

Then Allegro files are compiled and an executable is created :

Mezo3D : $(ALLEGRO_OBJS)
$(QUIET_LINK)$(FC) $(FFLAGS) -o $@ $(filter %.o,$ ˆ)

The same thing is done for the postprocessing files :

df6_pospro : $(POSTPRO_OBJS)
$(QUIET_LINK)$(FC) $(FFLAGS) -o $@ $(filter %.o,$ˆ)

An option named clean is available. This one allows to erase files with extensions .o, .mod, .d. The executable files, the solution files and the gmon.out files are also deleted.

It is possible to change this option if necessary.

Fort.10

This file is the one where the parameters of the mesh, the chemistry... are set.

GRID

This part is dedicated to the parameters of the mesh.

• Ndim : spatial dimension of the simulation (1 for a 1D simulation, 2 for a 2D one, 3 for a 3D one).
• Nx, Ny, Nz : number of points respectively in the directions ${\displaystyle {\vec {x}}}$, ${\displaystyle {\vec {y}}}$ and ${\displaystyle {\vec {z}}}$.
• Npx, Npy, Npz : number of processors allocated to the resolution, respectively in the directions ${\displaystyle {\vec {x}}}$, ${\displaystyle {\vec {y}}}$ and ${\displaystyle {\vec {z}}}$. The total number of points in the i direction is ${\displaystyle Ni.Npi}$, and the total number of processors used for the calculation is ${\displaystyle Npx.Npy.Npz}$.

Then the geometry of the box where the simulation will take place is defined.

• Xmax, Ymax, Zmax : maximal lenghts of the area, respectively in the directions ${\displaystyle {\vec {x}}}$, ${\displaystyle {\vec {y}}}$ and ${\displaystyle {\vec {z}}}$.
• IGridx, IGridy, IGridz : 1 if the grid is uniform in the considered direction, 0 otherwise.
• Alphax, Alphy, Alphaz : stretching ratios. (NO MORE AVAILABLE IN THE STAGGERED VERSION OF THE CODE)
• Betax, Betay, Betaz : stretching rates. (NO MORE AVAILABLE IN THE STAGGERED VERSION OF THE CODE)
• XP0, YP0, ZP0 : stretching positions. (NO MORE AVAILABLE IN THE STAGGERED VERSION OF THE CODE)

Finally, the kind of boundary conditions is defined with IType. 3 choices are given :

• Periodic (0),
• General (1),
• Symmetrical (2).

(NO MORE AVAILABLE IN THE STAGGERED VERSION OF THE CODE)

NSCBC

NSCBC means Navier-Stokes Characteristic Boundary Conditions. With these boundary conditions, the different types of input and output are taken into account. They have to be defined on each face of the box.

• Periodic (0)
• Non-reflecting outflow (1)
• Non-reflecting inflow (2)
• Hard inflow (3)
• Wall (4)

PARAM

The needed constants are defined here.

• Nsc : number of chemical species. For example, Nsc = 3 if there are three species in the calculation (for example, fuel, oxidant and a neutral species).

Remark: If we considere a reactive case, Nsc=1 will refer to the fuel, Nsc=2 will refer to the oxydizer.

• Re_ac : acoustic Reynolds of the flow.

${\displaystyle Re_{ac}={\dfrac {\rho {cL}}{\mu }}}$

c : sound speed ; L : characteristic length. At this stage the user should be clear in which unity he defines his geometry (1 mm, 1 cm ?), the density and the viscosity of the fluid at the reference temperature, as well as the sound speed. exemple: ${\displaystyle c=340m.s^{-1}}$, ${\displaystyle L=1mm=0.001m}$, ${\displaystyle \rho =1,2kg.m^{-3}}$, ${\displaystyle \mu =1,85.10^{-5}kg.m^{-1}.s^{-1}}$ lead to ${\displaystyle Re_{ac}=22050}$

• Gamma : ratio between heat capacity at constant pressure and heat capacity at constant volume.

${\displaystyle \gamma ={\dfrac {C_{p}}{C_{v}}}}$

• Pr : Prandtl number of the flow.

${\displaystyle Pr={\dfrac {\nu }{D_{th}}}}$

Where ${\displaystyle D_{th}}$ is thermal difffusivity ratio.

• Sc : Schmidt number of each species. Species 1 is the fuel, species 2 air and species 3 a neutral one.

REACTION

In this part, the parameters affecting the chemistry are defined.

• I_react : if 1, the chemistry will be resolved, if 0, it won't.
• CFL_react : CFL (Courant-Friedrichs-Lewy condition) of the reaction.

${\displaystyle CFL={\dfrac {u\Delta {t}}{\Delta {x}}}}$

• Alpha : thermal expansion ratio.

${\displaystyle \alpha ={\dfrac {T_{b}-T_{u}}{T_{b}}}}$

${\displaystyle T_{b}}$ : temperature of the burnt gases ; ${\displaystyle T_{u}}$ : temperature of the unburnt gases.

• Beta : Zeldovitch number.

${\displaystyle \beta =\alpha {\dfrac {T_{a}}{T_{b}}}}$

${\displaystyle T_{a}}$ : temperature of the activation barrier so that the reaction can happen.

• Damla_factor : Damköhler factor which defines the operational conditions of a reaction.

${\displaystyle Da=Kt1\beta .exp\left(-{\dfrac {\beta }{\alpha }}\right)}$

• NuF*MF and NuO*MO : stoichiometric coefficient respectively for the fuel and the oxidant.

${\displaystyle \displaystyle Y=\nu {M}}$

M : molar mass of the species.

• YF_O and YO_O : Mass fraction of fuel and oxidizer in the reference flows of fuel and of oxydizer.

Example : air and fuel are injected separately, so YF_O=1 and YO_O=0.23.

FILES

In this section, the back-up files are managed (name, creation frequency...). These files are storing raw data.

• PREFIX : prefix of the filenames.
• I_fresh_run : if 1, a new calculation is started. If 0, the code restarts a calculation at a specified point.
• I_read_unit : the code will continue the calculations starting from the file with this number.
• I_save_unit : number of the first new file while calculating.
• I_print : number of iterations that have to be done to print a message.
• CFL_no : settles the equation CFL.
• I_iter_save : number of iterations that have to be done to create a data file.
• Time_save : time to save the elapsed time.
• N_iter_stop : maximal number of iterations to stop the calculation.
• Time_end : maximal time to stop the calculation.

POST

This section deals with the post-treatment parameters.

• PREFIXpost : prefix of the files that will be subjected to the post-treatment.
• Ichoice : format of the output data. 1 to use Paraview, 2 to use Gnuplot.
• Xstart : first read data file that will be post-treated.
• Xend : last read data file.
• Xinc : increment to jump files.
• Ystart : first output file to be written.

At this stage some cleaning of the code would be useful for futur users.

• NormDir : plane normal direction (1,2,3). (WARNING in all_post NormDir may be read as REDIM in all_post.f90, it is a bug that hasn't been fixed yet, the best is to let all the following values at 0)
• PlaneCoord : position of the plane on XP2(NormDir), between 1 and ${\displaystyle Ni.Npi}$ for the direction i.
• REDIM : if 0, the results will be adimensional, if 1 they will be dimensional.
• Lo* : reference length (${\displaystyle 1m}$).
• Co* : reference sound speed (${\displaystyle 1m.s^{-1}}$).
• Rho_o* : reference density (${\displaystyle 1kg/m^{3}}$).
• To* : reference temperature (${\displaystyle 1K}$).

SCALAR

This part is made to let the user customize the code, by entering parameters that he could refer to using the subroutine LECTUS with the tag scalar as it is made for example in the file all_post.f90. The value of the scalars are given here.

• Uo : speed of the flow.
• Cv : heat capacity at constant volume.
• Rv, Xo, Yo : parameters for particular cases.

Fort.10.sh

Fort.10.sh is a file written in shell language in purpose to ease the creation of initial directories for the submission of new jobs.

As it has been first thought, it works with the file launch-run.sh. In launch-run.sh the variables NX, NY, NZ (that define the mesh size per processor) and the variables NPX, NPY, NPZ (that define the number of processors in each direction) are defined to be exported in fort.10.sh. In this way, several calculations can be launched by changing only a few parameters.

Initialization example :

sh launch-run.sh -g

The option -g calls the subroutine generate(). This subroutine creates a directory named by the variable "dir". In this directory, two other directories will be created. A directory "IN" will contain the input files while a directory "OUT" will contain the output files. Then the template of the executable "NAME" is redirected towards a file named fort.10. The file fort.10 will be located in the directory "IN". This step is realized like this :

sh $Name.template > $dir/IN/fort.10.

When the file fort.10 is created, the subroutine copies the executable calculating the case, the executable dealing the results as well as the file fort.11 in the directory "IN". Then the template of the job is redirected towards the file named "job".

Vtk_bin_write

This file concerns the writting of output files on several processors. The user may have to change flags in the file depending on the cluster he wishes to use. Two configurations are available and cover all the clusters used with the code. For example, it is possible to launch calculations either on ATON or on VARGAS. In that case, it is necessary to put the value 1 for the computer which is going to be used. For instance, if the used computer is ATON, then the value entered will be 1 and so, it will be 0 for the computer VARGAS :

#define ATON 1
#define VARGAS 0

Otherwise, if the used computer is VARGAS, then the value entered will be 1 and it will be 0 for the computer ATON :

#define ATON 0
#define VARGAS 1

Be careful, some options will be possibly modified in function of the cluster.

The distribution of processors can be defined in function of :

• Coordinates.
• Vectors : in that case, it will depend on the dimension of the problem, NDIM.
• Scalar : this function will define the id number of the processor.

Job

The job file is used if the calculation of the problem is made on a computer which is not local. The extension of this file is .ll because it is a LoadLeveler file.

A job for Antares has to be created like this :

• Working name : locating easier the calculations.

@ job_name = NAME 

• Standard exit file :

@ output = \(job_name).o\(jobid)

• Error exit file :

@ error = \(job_name).e\(jobid)

• Working type :

@ job_type = MPICH

• Total number of wanted processus :

@ total_tasks = NP

• Maximal elapsed time for the job :

@ wall_clock_limit = 24 :00 :00

• Maximal memory by MPI process : if a process exceeds the specified limit, then the calculation is killed (it will be written "out of memory" in the error exit file). Conversely, if this parameter is overestimated, the waiting time will be longer. It is preferable to execute calculations with a data limit subordinate or equal in 3gb. The wanted memory has to get the closest as possible to the real necessary memory by process.

@ data_limit = 550mb

• Initial directory to send : this directory contains all the data for the calculation.

@ cri_initialdir = PWD/dir/IN

• Final directory for the results : the computer has to send the results to this directory.

@ cri_finaldir = PWD/dir/OUT

• User mail : an e-mail will be sent in case of overtaking of elapsed time as well as when the calculation has ended.

@ notify_user = name@adress.fr

• Buffering working directory :

cd $LOCAL_WORK_DIR

The arborescence of the file must be defined below this line. Generally, it corresponds to the arborescence of the final directory.

The result files are moved in a directory named Spool before being copied in the final directory.

Launch-run.sh

This file is a tool and is not essential to run the code, but can save a lot of time in case of a parameters study. For exemple it may be interesting to have a way to change the parameters automatically in order to reduce the possibility of mistakes.

This file automates the changes to be done on the quantities of mesh points or processors in each direction. Instead of giving the value for Nx, Ny, Nz, Npx, Npy and Npz in the control file, a new file (template) is created from the original one. In this template, the value are replaced by a notation : in the direction i, the number of points is $NI and the number of processors is $NPI. The whole text is inserted in the command echo" ". The beginning of the file is shown below.

echo"
====================================================
ALLEGRO1-2-3D CONTROL FILE
====================================================
# GRID #

$NDIM                Ndim
$NX                  Nx
$NY                  Ny
$NZ                  Nz
$NPX                 Npx
$NPY                 Npy
$NPZ                 Npz

...
"

The value that will be used are then written in the launch-run.sh. This file has to begin with #!/bin/sh in order to be executable for the shell. Each case of calculation is written on a sole line.

#!/bin/sh

NDIM=3 NX=40 NY=40 NZ=40 NPX=1 NPY=1 NPZ=1 ./fort.10.sh $@
NDIM=3 NX=40 NY=40 NZ=40 NPX=2 NPY=1 NPZ=1 ./fort.10.sh $@

The template will pass on the data to the fort.10.sh, which in turn will create the fort.10 file according the value given by the launch-run.sh.

In this example, the first calculation is in 3D, uses a 40x40x40 mesh and works on one processor. The second one is also in 3D, but the size of the mesh is 80x40x40 and works on 2 processors.

When doing a list of cases, there are 2 possibilities to run it :

• the code runs once each case one after the other (as shown above)
• or the code runs only one case by one. In this case, only one line is not a comment, and the calculation has to be launched manually after every change.

All_init.f90

This file allows the user to initialize the problem, that is to say to give the initial conditions. For example a pressure peak can be defined, the density or the temperature can be fixed and so on.

In order to describe a function (for instance a pressure function) on all the points of the scalar grid S, 2 buckles are necessary to make vary the indexes I and J in the directions ${\displaystyle {\vec {x}}}$ and ${\displaystyle {\vec {y}}}$ (if it is a 2D case, otherwise for a 3D case, 3 buckles will be used by adding the index K). These indexes give knot coordinates. Since the indexes are appropriate for processors and not for physical problem, they can not be used to code a case, that is why variables XP2(I,1) and XP2(J,2) are used. I indicates the position in the direction x, J indicates the position in the direction y, 1 indicates the first direction and 2, the second one.

The speed is described by WK1 and WK2 if the case is two-dimensional, adding WK3 if the case is 3D. In a test case with a pressure peak, we consider these components as null.

The conservative variable Q(I,J,K,5) which corresponds to the energy (the number 5 indicates that is energy), is always defined as :

${\displaystyle Q(I,J,K,5)={\dfrac {P_{atm}}{\gamma -1}}+{\dfrac {1}{2}}{\dfrac {WK1^{2}+WK2^{2}+WK3^{2}}{Q(I,J,K,4)}}}$

The conservative variable Q(I,J,K,4) corresponds to the density.

The routine interpolate1 is called in order to pass a scalar from the grid S to the grid of one momentum component. Thanks to this routine, speeds are expressed on V1, V2 and V3 grids.

Organization of H-Allegro

The routines and subroutines used by the H-Allegro solver are presented below. They are shared out among 5 folders, namely :

1. DER
2. GRID
3. POST
4. SOLV
5. UTIL

Since the code uses various grids, it is necessary to define the changes that operations are doing. In other words, the differentiations or the interpolations include a grid change, so they have to be defined for each grid. The index 1 in the names of the subroutines stands for a change from the scalar grid to a momentum component one (called ${\displaystyle V_{i}}$), while the index 2 stands for the opposite change, going from the momentum component grid (${\displaystyle V_{i}}$) to the scalar one (S).

For example, if the derivative of ${\displaystyle \rho }$ (a scalar) is wanted on the scalar grid, there are two possibilities to proceed :

• the scalar is differentiated with Diff1 (which changes the grid from S to ${\displaystyle V_{i}}$) and then interpolated with Interp2 (changes the grid from ${\displaystyle V_{i}}$ to S).
• the scalar is interpolated with Interp1 (changes the grid from S to ${\displaystyle V_{i}}$) and then differentiated with Diff2 (changes the grid from ${\displaystyle V_{i}}$ to S).

In both cases, ${\displaystyle {\dfrac {\partial {\rho }}{\partial {x}}}}$ will be obtained on the S grid as wanted.

DER/DF6

The folder PADE is not described because the Pade scheme is not used in the MPI version of H-Allegro.

The number at the end of the names indicates that the calculation includes a change of grid. 1 means going from the S grid to the ${\displaystyle V_{i}}$ one, while 2 means the contrary, going from the ${\displaystyle V_{i}}$ grid to the S one.

• der.f90  : the subroutines to do the various interpolations and differentiations.
  • INIT_DER : initialization.
  • DERIVATIVE1, DERIVATIVE2: derivative of a function along the chosen axis.
  • INTERPOLATE1, INTERPOLATE2 : interpolation of a function along the chosen axis.
  • INTERP_SURF1, INTERP_SURF2 : interpolation of a function along the axis Iaxe, on a face with a normal Iface (Iface is perpendicular to Iaxe).
  • INTERP_LINE1, INTERP_LINE2 : interpolation of a function along the axis Iaxe, on a line with a normal Iface (Iface is perpendicular to Iaxe).

• der_sub.f90  : the subroutines to do the differentiations and interpolations, taking the boundary conditions into account.
  • echange : exchanges data with the nearby processors.
  • echange_line : exchanges data with the processors on another line.
  • echange_surf : exchanges data with the processors on another surface.
  • derivative1_on_slice, derivative2_on_slice : differentiation when the boundary conditions are symmetric.
  • derivative1_periodic_on_left, derivative2_periodic_on_left : differentiation when the boundary conditions are periodic on the left.
  • derivative1_periodic_on_right, derivative2_periodic_on_right : differentiation when the boundary conditions are periodic on the right.
  • derivative1_general_on_left, derivative2_general_on_left : differentiation when the boundary conditions are general on the left.
  • derivative1_general_on_right, derivative2_general_on_right : differentiation when the boundary conditions are general on the right.
  • interpolate1_on_slice, interpolate2_on_slice : interpolation when the boundary conditions are symmetric.
  • interpolate1_periodic_on_left, interpolate2_periodic_on_left : interpolation when the boundary conditions are periodic on the left.
  • interpolate1_periodic_on_right, interpolate2_periodic_on_right : interpolation when the boundary conditions are periodic on the right.
  • interpolate1_general_on_left, interpolate2_general_on_left : interpolation when the boundary conditions are general on the left.
  • interpolate1_general_on_right, interpolate2_general_on_right : interpolation when the boundary conditions are general on the right.

GRID

• all_dimg.f90  : initialization of the dimensions and the various variables.

• grid_bounds.f90  : manages the reading of the domain faces on the different grids.
  • INIT_BOUNDS : determines the presence of faces, edges and corners NSCBC3D.
  • COPY_FACE_from_S, COPY_FACE_from_Vi, COPY_FACE_from_Vij : copies faces according the grid.
  • SUB_COPY_FACE : algorithm to copy the data of a face.
  • PASTE_FACE_to_S, PASTE_FACE_to_Vi : pastes according the grid.
  • SUB_PASTE_FACE : algorithm to paste the data of a face.
  • ADD_FACE_to_S, ADD_FACE_to_Vi : adds a face according the grid.
  • SUB_ADD_FACE : algorithm to add a face.

• grid_dim.f90  : creation of the grid.
  • LECTUS_GRID : initializes the grid with the parameters, the number of points and processors along each directions, checks that the indicated number of processors is correct.
  • GETINDEX : indexes the points.

• grid.f90  : initialization of the grid.
  • INIT_GRID : initializes the grid, taking into account the overlaps and calculates the size of a cell.

• plot.f90  : prepares the use of Gnuplot.
  • SAVE_GNUPLOT : creates files with the coordinates of the point and the corresponding values.
  • SAVE_GNUPLOT_FIELDS : creates files with the fields at different times.
  • SAVE_GNUPLOT_1D : does the same work as SAVE_GNUPLOT for 1D simulation.

• post_plot.f90 : creates files for the post-processing.
  • PARAVIEW_BINARY : manages the creation of files according the number of dimensions for Paraview.
  • OUT_GNUPLOT : manages the creation of files according the number of dimensions for Gnuplot.

• vtk_bin_write.c  : manages the distribution of processors and the choice of the calculator.

POST

• all_post.f90  : manages the post-processing, with the initialization and the creation of files.
  • LECTUS_POST : reads the data.
  • INIT : initialization of the variables on the various grids.
  • GETVECT : resizing or not.
  • ALLOCATE_FIELDS : allocation of the fields according to the names.
  • COMPUTE_FIELDS : resizes if needed and calculates the fields.
  • READ_FILES : depending on the rank of the processor, opens a new file and reads in another one the time and the values of the scalars.
  • ROTATIONNEL : calculates the curl along the number of dimensions.
  • PROJECT : projection of a vector component onto the S grid.
  • DILUTION : calculation of the dilution (the thermal diffusivity ratio varies according to the space).

• all_post_cut.f90  : all_post_f90 without some subroutines (GETVECT, ALLOCATE_FIELDS, COMPUTE_FIELDS, ROTATIONNEL, PROJECT and DILUTION), but has another one.
  • SAVE_T : creates files keeping the values of the temperature.

SOLV

• all_chemistry.f90  : manages the chemistry of the code.
  • LECTUS_REACT : initializes the parameters of the reaction.
  • ARRHENIUS : calculates Arrhenius law (${\displaystyle {\dot {\omega }}}$).
  • ENERGY_SOURCE : calculates the RHS of the energy source term.
  • SCALAR_SOURCE : calculates the RHS of the scalar source term.
  • REACTION_RATE : calculates the reaction rate.
  • CHECK_CFL_REACT : calculates the time-step of the reaction.
  • FLAME_SPEED : calculates the propagation speed of the flame.

• all_solver.f90  : resolution of the case.
  • NEGATIVE : spotting of the bugs.
  • CONVECTION1 : calculates the convection, the result is on the S grid.
  • CONVECTION2 : calculates the convection, the result is on the ${\displaystyle V_{i}}$ grid.
  • CONTINUITE : calculates ${\displaystyle {\dfrac {\partial {\rho {U_{j}}}}{\partial {x_{j}}}}}$ (continuity equation) on the S grid.
  • ENERGY : calculates ${\displaystyle {\dfrac {\partial {\rho {E}}}{\partial {t}}}}$ (conservation of energy) on the S grid.
  • MOMENTUM : calculates ${\displaystyle {\dfrac {\partial {\rho {U_{i}}}}{\partial {t}}}}$ (conservation of momentum) and the stress tensor ${\displaystyle \tau _{ij}}$.
  • SCALAR : calculates the scalars and the diffusion terms.
  • RHD_INIT : initializes the variables (speed, pressure, temperature...).
  • RHD : calculation of the variables (energy, momentum...), the reaction rate and the boundary conditions.
  • TIME_STEPPING : calculates the time advancement (solver Runge-Kutta, order 3).
  • CHECK_CFL : calculates the time-step.

• all_dim.f90  : declares the dimensions of the functions send_left, recv_left...

• all_equ.f90  : calculates some variables.
  • RHOETHERM : calculates ${\displaystyle \rho {C_{v}T}=\rho {E}-{\dfrac {1}{2}}\rho {U_{i}^{2}}}$.
  • SOUND : calculates the sound speed ${\displaystyle c={\sqrt {\rho {C_{v}T}}}}$.
  • TEMP_PRES : calculates the pressure and the temperature.
  • ENTROPIE : calculates the entropy.
  • VISCOSITY : calculates the viscosity.

• all_nscbc.f90  : manages the boundary conditions.
  • LECTUS_NSCBC : allocates the type of boundary conditions to each face.
  • INIT_BC : initialization and allocation of tables for the faces.
  • COPY_BOUNDS_S, COPY_BOUNDS_Vi, COPY_BOUNDS_Vij : copy of a value in the chosen grid on the faces, edges and corners of the S grid.
  • OUTFLOW_3D : correctives terms for the DQM to the faces (correction of normal convection terms and viscosity terms).
  • SOFT_INFLOW_3D : change of the DQM on a face (adding the transverse terms of convection).
  • HARD_INFLOW_3D : change of the DQM on a face (values of the speed and energy forced, calculation of the transverse terms, viscosity and convection terms).
  • NOSLIP_WALL_3D : change of the DQM on a face (transverse and normal values of ${\displaystyle {\dfrac {\partial {\rho {U}}}{\partial {t}}}}$ forced to 0, calculation of the transverse terms, viscosity and convection terms and ${\displaystyle {\dfrac {\partial {\rho {E}}}{\partial {t}}}}$).
  • FACE_PERIODIC : change of the DQM (correction of the normal convection terms).
  • NSCBC : implementation of the changes for the faces.

• all_param.f90  : reading of some parameters and calculations.
  • LECTUS_PARAM : reading of the parameters (Sc, ${\displaystyle Re_{ac}}$, ${\displaystyle \gamma }$, Pr) and calculation of ${\displaystyle {\dfrac {r}{C_{v}}}}$, the atmospheric pressure and the perfect gas constant.
  • LECTUS_FILES : reading of the parameters for the output and the time-stepping.

UTIL

• lectus.f90
  • LECTUS_MARK : reads and compares marks.

• minmax.f90  : search of minimum and maximum.
  • GETMIN : finds the minimum of a function and the corresponding point.
  • GETMAX : finds the maximum of a function and the corresponding point.
  • GETERRMAX : finds the maximal error.
  • GETMINMAX : finds the maximum and the minimum of a function and the corresponding points.

• mtfort.f90  : generates pseudorandom real numbers.

• random.f90  : generates various pseudorandom numbers.

• thi.f90  : generates an Isotropic Homogeneous Turbulence (THI) for a colocated grid.
  • LECTUS_PREFIX : reads the prefix of the files.
  • LECTUS_THI : reads the parameters used for the THI.
  • ADD_BOUNDS2D : adds boundaries in 2D.
  • ADD_BOUNDS3D : adds boundaries in 3D.
  • SAVE_FIELDS : creates files with the values of the fields.

Conditioning

After some calculations, the distribution of the cores may be changed in order to continue the calculation or to do the post-treatment :

• increase the cores number with the help of the routine Divise,
• decrease the cores number with the help of the routine Reunite.

Divise

The output files of the previous calculation will be used, so it is advisable to copy them into the "Divise" folder and work with the copy. In fact, the only file to change is create2.sh.

• The first step is to make sure that the division according the chosen direction will give a number of points which is an integer (use even numbers).

For example, the previous configuration was :

NX=80  NY=50  NZ=40
NPX=1  NPY=2  NPZ=2

You may want to divise according the ${\displaystyle {\vec {x}}}$ direction, and the new configuration will be :

NX=40  NY=50  NZ=40
NPX=2  NPY=2  NPZ=2

But if in the previous configuration NX=85, the division is impossible.

• Specify the back-up times of the output files, which corresponds to the last part of the suffix (e.g. for TESTCF0001_0007, the back-up time is 7).

for T in 'seq -f %02g ${\displaystyle n}$ ${\displaystyle m}$'; do

${\displaystyle n}$ is the first back-up time to consider and ${\displaystyle m}$ is the last back-up time.

• Specify the input and output data. The index e stands for the input while the index s stands for the output. Nvar corresponds to the variables number : if there are no chemistry, Nvar=5 (${\displaystyle \rho }$, E, ${\displaystyle U_{i}}$), otherwise Nvar=5+Nsc. The name of the input files is defined in the fort.10.

Using the example above with a division along the ${\displaystyle {\vec {x}}}$ axis :

NXe=80 NYe=50 NZe=40 NPXe=1 NPYe=2 NPZe=2 NXs=40 NYs=50 NZs=40 NPXs=2 NPYs=2 NPZs=2 Nvar=8 Tps_save=$T ENTER=TESTCF EXIT=OUTPUT sh./create_xyz -x $@
gfortran filereading.f90 divise_x.f90 -o divX.x
./divX.x
rm divX.x
rm divise_x.f90
rm TESTCF*

If the division is done on another axis, the option for the execution of create_xyz.sh should be changed, as well as the name of divise_x.f90 and divX.x. When several divisions are planned, the code lines above are copied and changed as wanted. The routine will be executed with the command "sh create2.sh".

• Before running the new calculation, the file fort.10 has to be changed in order to have the correct data (GRID and FILES parts).

The FILES part for the example is as follows :

# FILES #
OUTPUT      PREFIXpost
1           I_fresh_run
${\displaystyle m}$          I_read_unit
${\displaystyle m+1}$       I_save_unit

Reunite

The output files of the previous calculation will be used, so it is advisable to copy them into the "Reunite" folder and work with the copy. The only file to change is create.sh.

• Specify the back-up times of the output files.

for T in 'seq -f %02g ${\displaystyle n}$ ${\displaystyle m}$

${\displaystyle n}$ is the first back-up time and ${\displaystyle m}$ is the last one.

• Specify the input data. The output are not described, so take care when doing several concatenations, the output of the first concatenation is only written as the input of the second one.

Npx=2 Npy=2 Npz=2 Nx=40 Ny=50 Nz=40 Nvar=8 OUTPUT=RESULT INPUT=TESTCF T=$@ sh reunix.sh > RX.f90 $@
gfortran RX.f90 filereading.f90 -o X.x
./X.x
rm X.x
rm RX.f90
rm TESTCF*

In this example, the configuration before the concatenation is :

NX=40  NY=50  NZ=40
NPX=2  NPY=2  NPZ=2

And after the concatenation, there is only one core for the ${\displaystyle {\vec {x}}}$ direction :

NX=80  NY=50  NZ=40
NPX=1  NPY=2  NPZ=2

If the division is done on another axis, the name of reunix_x.f90 and RX.x should be changed. When several concatenations are planned, the code lines above are copied and changed as wanted. The routine is executed with the command "sh create.sh".

• Before running the new calculation, the file fort.10 has to be changed according the new configuration.

Specific case

In the case of calculations with the cluster Antares, there may be some compatibility troubles. The option of compilation -frecord-marker can be used to solve these troubles. There are two possibilities :

• -frecord-marker=4 this option has to be written in the makefile before the compilation with the other options (FFLAGS). The output files are then in the correct format.
• -frecord-marker=8 this option is written in the file create2.sh (for the division) or create.sh (concatenation) as follows :

gfortran filereading.f90 divise_x.f90 -frecord-marker=8 -o divX.x
gfortran RX.f90 filereading.f90 -frecord-marker=8 -o X.x

Post-treatment tools

The post-treatment is managed by the file all_post.f90.

• Parameters : variables definition (${\displaystyle \gamma }$, ${\displaystyle \gamma -1}$, ${\displaystyle P_{atm}}$, XNuF_MF, XNuO_MO, YF_O, YO_O).
• Dimensions and allocations.
• Init Grid Derivative :
  • Average calculation :${\displaystyle Average=\sum {\dfrac {moyenne(IP,N)}{compteur}}}$               With ${\displaystyle compteur=Xend-Xstart}$
  • Y and id writing.
  • Vtk (Ichoice = 1, Paraview), Gnu (Ichoice = 2, gnuplot) files creation or slices creation, vtk files (Ichoice = 3, Paraview). Be careful, this last option is only available for a three-dimensional case.
  • Calculated fields writing.
• Lectus post :
  • Reading of the quantities number to be calculated (cf fort.10, fort.11).
  • Lines number counting.
  • Reading of the variables names to be calculated.
  • Reading of the files names to be read.
  • Post-processing data reading (Prefix post, Ichoice, Xstart).
• Init :
  • Speed field initialization on grid Vi.
  • Speed gradient initialization on grid S or Vij.
  • Speed divergence initialization on grid S.
  • Temperature and pressure initialization on grid S.
  • Temperature gradient initialization on grid Vj.
  • Viscosity initialization on grid S. Be careful, it is dynamic viscosity, ${\displaystyle \mu }$.
• Getvect :
  • Calculation of WQ(IP,Icmp), with WQ varying from 1 to Ncmp :

if REDIM = 0 ${\displaystyle \displaystyle WQ(IP,Icmp)=Q(IP,Icmp)}$, else ${\displaystyle \displaystyle WQ(IP,Icmp)=Q(IP,Icmp)xcoeff_{redim}}$

• Allocate Fields :
  • Nfields calculation in function of the variables : ${\displaystyle Nfields=\sum _{1}^{NIn}Ndim}$, or ${\displaystyle Nfields=\sum _{1}^{NIn}Nsc}$, or ${\displaystyle Nfields=\sum _{1}^{NIn}Nfields+1}$
  • This subroutine defines variables which can be post-processed. Check up that variables correpond to what you want. The available variables are : Pres, Temp, Rho, Rho U, Rho E, Y, U, We, Mach, Mach T, Mach P, Vel, SL, Grad T, Grad Z, Grad Y1, Div U, Rot U, e, S, FP, Z, C, Visc, Dilut, Q criteria, tak, c, Epsi.
• Compute Fields :
  • Calculation of curls, vectors, scalars.
  • Calculation of WQ(IP,Ifield).
• Read Fields :
  • Vtk files construction.
• Vorticity :
  • The curl calculation in function of Ndim : RotX, RotY, RotZ.
• Critere Q :
  • ${\displaystyle Q\_criteria={\dfrac {1}{2}}(Q_{ij}Q_{ij}-S_{ij}S_{ij})}$
• Takeno (cf : Thesis of Marianne Sjostrand)
  • Takeno index varies from -1 to 1. It is defined as : ${\displaystyle tak={\dfrac {{\dfrac {\partial Y_{F}}{\partial x_{i}}}{\dfrac {\partial Y_{O}}{\partial x_{i}}}}{|{\dfrac {\partial Y_{F}}{\partial x_{i}}}||{\dfrac {\partial Y_{O}}{\partial x_{i}}}|}}}$
• Grad Grad :
  • The index of Takeno is defined as : ${\displaystyle tak=\sum {\dfrac {\partial c}{\partial x_{i}}}{\dfrac {\partial c}{\partial x_{i}}}}$
• Grad Croise :
  • The index of Takeno is defined as : ${\displaystyle tak=\sum {\dfrac {\partial c}{\partial x_{i}}}{\dfrac {\partial z}{\partial x_{i}}}}$
• Project :
  • Calculation of projection : ${\displaystyle S=U_{i}{\dfrac {\dfrac {DF}{Dx_{i}}}{|{\dfrac {DF}{Dx_{i}}}|}}}$
• Dilution :
  • Calculation of the coefficient D : ${\displaystyle D=\sum Q(IP,Ivar)}$
• Dissipation :
  • Calculation with the field of immediate speeds.

${\displaystyle \epsilon =\tau _{ij}{\dfrac {\partial U_{i}}{\partial x_{j}}}}$              ${\displaystyle \tau _{ij}={\dfrac {\mu }{e}}({\dfrac {\partial U_{i}}{\partial x_{j}}}+{\dfrac {\partial U_{j}}{\partial x_{i}}}-{\dfrac {2}{3}}{\dfrac {\partial U_{k}}{\partial x_{k}}})}$

Test cases

Scalability Test

Folder : Poiseuille.

It is important to keep in mind that the results of the test will depend on the case and on the used computer. For this test, a simplified 3D Poiseuille flow is used, with simple physics (no speed, constant density...). In order to make a test case of the code on the cluster Antares, it is enough to inform these files : job.ll.template and launch-run.sh.

Job.ll.template

• Modify your e-mail address in order to receive information on the calculation progress.
• Possibly, modify the parameters "wall_clock_limit" and "data_limit" to gain places in the waiting line or to increase the memory attributed by process.

Jobs submission

• Compile files :

make

• Create files for every test case :

sh launch-run.sh -g

• Submit jobs to a cluster :

sh launch-run.sh -s

The terminal asks if the job "XXXXX" is to be submitted. Just enter "y". If an error occurs saying the file is not executable, make it executable with the command :

chmod 777 "filename"

Strong Scaling

The goal of this part is to define the optimal size of the blocs. The total size of the mesh is fixed, for example 80x80x80, and the load per core will vary when the cores number changes. The different configurations of the cores distribution are presented below.

Total core number	Cores number (${\displaystyle {\vec {x}}}$)	Cores number (${\displaystyle {\vec {y}}}$)	Cores number (${\displaystyle {\vec {z}}}$)
1	1	1	1
2	2	1	1
	1	2	1
	1	1	2
4	2	2	1
	2	1	2
	1	2	2
8	2	2	2
16	4	2	2
	2	4	2
	2	2	4
32	4	4	2
	4	2	4
	2	4	4
64	4	4	4
128	8	4	4
	4	8	4
	4	4	8
256	8	8	4
	8	4	8
	4	8	8
512	8	8	8

For example, if the code runs on 64 cores, which means 4x4x4 cores, the size of the bloc for a core is 20x20x20.

Speed-up

The acceleration is defined by : ${\displaystyle acc(n)={\dfrac {T_{ref}}{T_{n}}}}$             (see Scalability for more information)

Graphique !

Parallel efficiency

The parallel efficiency is defined by : ${\displaystyle eff(n)={\dfrac {acc(n)}{\dfrac {n}{n_{ref}}}}}$             (see Scalability for more information)

Graphique !

Optimal size of blocs

The optimal size of blocs is found with the help of the reduced efficiency : ${\displaystyle E=t.{\dfrac {n}{N_{noeuds}N_{ite}}}}$             (see Optimum size of blocs for more information)

It is plotted against the total number of nodes per core. The curve is supposed to have a plateau corresponding to the optimal size of the blocs.

Graphique !

In this case, the plateau, which is actually quite small, is observed for a bloc size between 20x20x20 and 40x40x40.

Weak Scaling

Launch-run.sh

Having determined thanks to the strong scaling an optimal slice for the block size, it is possible to obtain the optimal block size with the weak scaling. Thus tests are realized for a block size of 40x40x40, 32x32x32 and 20x20x20 on various combinations of cores. You can choose your own block size.

Distribution of cores :

Total core number	Cores number (${\displaystyle {\vec {x}}}$)	Cores number (${\displaystyle {\vec {y}}}$)	Cores number (${\displaystyle {\vec {z}}}$)
1	1	1	1
2	1	1	2
	1	2	1
	2	1	1
4	1	1	4
	2	2	1
	2	1	2
8	1	1	8
	2	2	2
	4	2	1
	2	4	1
	2	1	4
16	2	1	8
	1	2	8
32	2	2	8
	4	1	8
	1	4	8
64	4	2	8
	2	4	8
128	8	2	8
	2	8	8
256	8	4	8
	4	8	8
512	8	8	8

Example of a test line :

NX=32 NY=32 NZ=32 NPX=8 NPY=8 NPZ=8 ./fort.10.sh $@

That means that the case is tested with a block size of ${\displaystyle 32^{3}}$ on 8x8x8 (512) cores.

Data Processing

When the jobs are submitted and the calculation is realized, the computer sends files results to the final directory which was defined in the template of the job. Execution times are in files log.

Parallel efficiency calculation The objective is to get the closest as possible to the theoretical parallel efficiency whose value is 1. To obtain the parallel efficiency of the studied test case, the following equation has to be used :

${\displaystyle eff(n)={\dfrac {T_{ref}}{T_{n}}}}$             (see Scalability for more information)

figure

The parallel efficiency is drawn in function of the cores number. The principle of the weak scaling is to keep a constant load by core, so the total size of the mesh increases. This explains why the parallel efficiency does not increase and stabilizes from a certain cores number.

Speed-up calculation

The objective is to obtain a speed-up which keeps the shape of the theoretical speed-up curve. The equation allowing to determine this speed-up is the following one :

${\displaystyle acc(n)={\dfrac {n}{n_{ref}}}{\dfrac {T_{ref}}{T}}}$             (see Scalability for more information)

figure

Determination of the optimal block size

An optimal zone for the block size is determined thanks to the strong scaling. Tests with the weak scaling are realized on various block sizes belonging to this slice. It is then enough to observe the speed-up and efficiency curves and to determine the block size which will give the best efficiency with the best speed-up.

Example :

After analysing the strong scaling results, it was deducted that the optimal block size is between 20x20x20 and 40x40x40.

	20x20x20	32x32x32	40x40x40
Efficiency	0.30	0.49	0.55
Speed-up (512 cores)	19.38	31.47	35.13

figure acc 20x20x20

figure acc 32x32x32

figure acc 40x40x40

In our case, the optimal block size is 40x40x40 since it has a better efficiency and has a better speed-up for a number of points higher than the other block sizes.

Cylindrical bidimensional flame

Folder : FLAM2Dcyl.

This test case corresponds to a turbulent flame in 2D, with a cylindrical shape. The size of the box is 20cmx20cm and the mesh is 3200x3200.

The growing flame is introduced in an isotrop homogeneous turbulence (THI) which is decreasing and initialized on a non-periodic domain.

Spherical tridimensional flame

Folder : FLAM3Dsph.

This test case corresponds to a turbulent flame in 3D, with a spherical shape. The size of the box is 3cmx3cmx3cm and the mesh is 480x480x480.

As the cylindrical bidimensional flame case, the growing flame is introduced in an isotrop homogeneous turbulence (THI) which is decreasing and initialized on a non-periodic domain.

Tridimensional meso-combustor

Folder : MEZO3D.

This case consists of a fuel-air mixing in a meso-combustor in three dimensions. There are two inputs (fuel and air), and an output in the wall.

• all_chemistry.f90 : this file allows to adjust the speed of a flame according to the local equivalence ratio for the laminated cases. If the case doesn't use laminated chemistry then the function XK_CH4_Air(Phi) must be deleted :

WK1(IP) = EXP(TaM/Tem(IP)+Ta_over_Tb)*(Q(IP,IF)*Q(IP,IO)/Q(IP,4))*XK_CH4_Air(Phi)

• all_nscbcmezopropre.f90 : Navier Stokes characteristic boundary conditions of faces. The Neumann conditions are applied to the wall :

${\displaystyle {\dfrac {\partial Y_{F}}{\partial x_{1}}}=0}$

3 points are used to described the domain : P, PP and PPP. P is on the face, PP and PPP are in the volume. For these points, the variable Y is considered and not ${\displaystyle \rho Y}$, so it is necessary to divise ${\displaystyle \rho Y}$ by ${\displaystyle \rho }$.

• communs.f90 : coordinates of the fuel and air entries centers, this file contains variables common to several files.
• all_init_propre.f90 : this file is used to start a calculation. Twall and Twall2 allow to create a fresh gas heart surrounded by a warm gas. A linear regression is used to describe temperature closed to the wall if ${\displaystyle T\leq TT}$. The fuel and air entries are defined (radius, speed of fluids, temperature, initial mass fractions).
• all_init_restart.f90 : this file is used to take back a calculation from a output file. In this file, the fuel and air entries are not defined. If you want to run the case from a save file, you have to enter the value 0 for the variable I_fresh_run (cf fort.10) and to indicate the number of the chosen save file. If you want to run the case from an output file and to modify some properties (for instance, the wall temperature), you have to enter the value 2 for the variable I_fresh_run
• all_solver.f90 : this file is used to solve the case (for example : error treatment, continuity, convection). A patch of viscosity must be added.
• fort.11 : list of variables which are post-processed.
• allegro.f90 : initializes the calculation and solves the case.

Derivative routines test

Folder : Test1_sin.

With this test case, it is possible to verify the derivative routines for a simple sinus.

NSCBC conditions test

Folder : Test4_Vortex2Dnscbc.

This case tests the NSCBC (Navier-Stokes Characteristic Boundary Conditions) in 2D or 3D. These conditions can be changed in the file vortex2D.template in the "NSCBC" part. As reminder, there are 5 boundary conditions :

• periodic
• non-reflecting outflow
• non-reflecting inflow
• hard inflow
• wall.

Flame speed

Folder : VITESSE_DE_FLAMME.

This case tests the flame speed along the local equivalence ratio. The chosen values of the local equivalence ratio are given as a parameter (Phi) in the file "launch-runs.sh".

Profiling and Scalability

Optimum size of blocs

In order to avoid the loss of time, it is advised to use the optimum size of blocs and the optimum cores number, which depend on the studied case as well as the cluster.

The optimum size of blocs is determined with the help of the reduced efficiency :

${\displaystyle E=T_{n}.{\dfrac {n}{N_{noeuds}N_{ite}}}}$

${\displaystyle T_{n}}$ is the execution time for ${\displaystyle N_{ite}}$ and a total number of nodes ${\displaystyle N_{noeuds}}$.

The values are plotted against the total number of nodes per core ${\displaystyle \textstyle \left({\dfrac {N_{noeuds}}{n}}\right)}$. This curve is supposed to have a plateau corresponding to the optimum size of the blocs. If it is not the case, the test has to be done again with smaller sizes of blocs. Once the optimum size of blocs found, it will be used in the following tests and calculations.

To be more precise, the Weak Scaling method is used to compare the evolutions of speed-up and parallel efficiency according the cores number. The configuration with less difference with the theory will be favoured. The optimum cores number is then defined for the chosen configuration using the evolution of the speed-up : when the real curve stops following the theorical one, the optimum number has been exceeded.

Scalability

Principle

When a code works on many processors, the calculation domain is divided and the code may lose performances because of the communications between the processors. The scalability test ^[5] checks the behaviour of the code when the number of cores increases, especially that the performances are quite stable.

Two tests are possible :

• Weak Scaling, the work load per core remains constant. This means that when the number of cores increases, the size of the cells decreases.
• Strong Scaling, the size of the problem remains constant whatever the number of cores working. The work load decreases when the number of cores increases.

Some conditions have to be followed when testing the code. The maximal number of cores ${\displaystyle n_{max}}$ is defined by the biggest planed configuration. From this, a number of cores for the reference execution can be specified :

${\displaystyle n_{ref}\leq {\dfrac {n_{max}}{4}}}$

Moreover, it is advised to use up to 3 intermediate values in order to visualize the evolution of the speed-up and the parallel efficiency.

• Weak Scaling :
  • Speed-up ${\displaystyle acc(n)={\dfrac {n}{n_{ref}}}{\dfrac {T_{ref}}{T_{n}}}}$
  • Parallel efficiency ${\displaystyle eff(n)={\dfrac {T_{ref}}{T_{n}}}}$

• Strong Scaling :
  • Speed-up ${\displaystyle acc(n)={\dfrac {T_{ref}}{T_{n}}}}$
  • Parallel efficiency ${\displaystyle eff(n)={\dfrac {acc(n)}{\dfrac {n}{n_{ref}}}}}$

In both cases, the theoretical speed-up is defined by : ${\displaystyle acc_{theo}(n)={\dfrac {n}{n_{ref}}}}$

The value of the theoretical parallel efficiency is 1.

Since the goal of this test is also to find out the maximal number of cores used without loss of perfomances, it is not necessary to run the code with a high number of iterations for the calculation. Actually, 10 iterations are enough to test.

Execution time

In order to determine the speed-up and the parallel acceleration, the execution time is required. The box below shows the way to have the execution time for the first processor, which orders the others, printed in a new file.

character*13:: FILEMPI
real(8) :: cpu_time

cpu_time = MPI_WTime()
# Time Stepping

temps_ex = MPI_Wtime() - cpu_time
if(rang_dans_monde.EQ.0) then
  FILEMPI='infotemps.txt'
  open(unit=1234567,file=FILEMPI,form='formatted',status='replace')
  write(1234567,*) temps_ex
  close(unit=1234567,status='keep')
  endif
# End

This way works for a calculation done on a local cluster, but it has to be changed for a calculation on Antares for example :

if(rang_dans_monde.EQ.0) then
  temps_ex = MPI_Wtime() - cpu_time
  print*, 'ALLEGRO End_RN, OK', temps_ex
  print*, '=========================================='
  endif

The job.ll.template is then slightly modified :

# Execution d'un programme MPI
time mpirun.Impi ./$NAME > log

The data usually printed on the screen (including the execution time) will be available in the "log" file.

Profiling

Profiling is a way to optimize the code skill. It consists in the determination of the parts which use most resources.

Thanks to the command "gprof", it is possible to create a report of execution. The executable which calculates the case must be compilated with the option -pg, then it must be executed in order to obtain the report, gmon.out.

The file gmon.out is divided into 2 parts :

• Call-graph profile : list of the routines in decreasing order of the consumed CPU time as well as the routines which are called from them.
• Flat profile : list of the routines in the decreasing order of the consumed CPU time and the number of calls.

Analysis of the flat profile

Description

• % time : percentage of total CPU time consumed by the routine.
• Cumulative seconds : partial sum of CPU time.
• Self seconds : CPU time of a routine.
• Calls : number of calls of a routine.
• Self ms/call : average of CPU time consumed by call to the routine, exclusive time.
• Total ms/call : total time (exclusive + inclusive) consumed by call to the routine.

Analysis

Having to test a calculation on various combinations of processors, it seems clearly that the most consumer routines are :

• der_sub_MOD_interpolate1_on_slice
• der_sub_MOD_derivative2_on_slice
• all_solver_MOD_momentum
• der_sub_MOD_interpolate2_on_slice

For instance for a calculation on 4x1x1 processors, the flat profile will be the following :

figure

It indicates that the subroutine derivative2_on_slice represents 20 % of total consumed CPU time. Thanks to the cumulative seconds, we can deduce that the total consumed CPU time egals 0.25 second. The CPU time of the subroutine derivative2_on_slice is about 0.05 second. This subroutine is called 545 times. We can say that the average number of milliseconds spent in this subroutine and its children per call is 0.09 ms.

Remark : If the calculation is run on 2x2x1 instead of 4x1x1, the number of calls of the subroutines derivative and interpolate periodic are divided by 2, while the number of calls of the subroutines derivative and interpolate on slice increases.

Analysis of the call-graph profile

This part is very interesting because it allows to know which routine calls another one and how long it spends. The routines are indexed by decreasing order of CPU time comsumption.

This report details the percentage of the total time that was spent in a function and its children. The total of these percentages is not equal to 100%. For instance the run time percentage for the main egals 100% because it is a routine which takes a lot of time, but the routine solving the time stepping has a run time of 97.3%. This report gives also the amount of time that was propagated from the function into the parent (self), the amount of time that was propagated from the children into the function (children), the number of calls from a routine to another.

Figure

For instance, if we consider the subroutine der_sub_MOD_echange (5), we can see that the parent functions are der_MOD_interpolate2 (16), der_MOD_derivative1 (15), der_MOD_derivative2 (11), der_MOD_interpolate1 (8). This subroutine is called 1014 times by the routines der_MOD_interpolate2 (93 calls), der_MOD_derivative1 (174 calls), der_MOD_derivative2 (327 calls) and der_MOD_interpolate1 (420 calls). The time spent in this function is about 0.01 second.

The children functions of this subroutine are der_sub_MOD_interpolate1_on_slice (7), der_sub_MOD_derivative2_on_slice (9), der_sub_MOD_derivative1_on_slice (17), der_sub_MOD_interpolate2_on_slice (18), der_sub_MOD_interpolate1_periodic_on_right (32) and so on.

The total amount of time propagated into this subroutine by its children is 0.13 second.

We can see that the subroutine echange calls 545 times the subroutine derivative2_on_slice. This last one is called only by the subroutine echange since it is indicated 545/545.

The subroutine echange spends 56 % of the total time in run time.

Optimization of communication time

The run time consists of the calculation time and the communication time. Profiling aims also to minimize this communication time which can be decomposed in the time of latency, additional cost and transfer.

The additional cost is the time of preparation of a message.

How to optimize the communication time ?

The first step is to cover the communications by the calculations.

Then, it is necessary to avoid copying out messages in buffering memory spaces and to minimize additional costs led by the repetitive calls to the subroutines of communication.

Finally, it is better to send an only big message rather than several small messages.

How to know the execution time of a routine ?

This information is given by the subroutine of MPI communication : MPI_WTIME().

This subroutine informs the user of the run time as it is called. Therefore, if this subroutine is called at the beginning and at the end of a routine, the run time can be deducted.

What way is it preferable to use for sending and the reception of messages ?

Standard sending : If the MPI makes a buffering copy of the message then the sending ends when this one is ended. Otherwise, the sending ends when the reception of the message is ended.

Blocking synchronous sending : This sending way uses the subroutines MPI_SSEND() and MPI_RECV().

The blocking synchronous sending allows to avoid buffering copies of messages, in this way it also avoids additional costs.

Thanks to it, a factor 2 is gained with regard to the standard mode.

Non-blocking synchronous sending  : This sending way is chosen in H-Allegro. It uses the subroutines MPI_SSEND() and MPI_RECV().

The non-blocking way allows to cover the communications by the calculations. The sending of a message ends when the reception was posted and when the reading of the message is ended.

The sending is coupled with the reception thus a factor 3 is gained with regard to the standard sending.

So that this mode can work correctly, it is necessary to add MPI_WAIT() after MPI_ISSEND() and MPI_IRECV(). It allows to synchronize the process until the end of the request.

Model of calculation

Time discretization

A third order Runge-Kutta method is used. This one allows to improve the precision of time integration and the stability.

The equations used in H-Allegro are defined below :

${\displaystyle {\dfrac {Q_{1}^{1}-Q_{1}}{\Delta t/3}}+{\dfrac {8}{5}}RHS(Q_{1})=0}$

${\displaystyle Q^{2}=Q^{1}-{\dfrac {\Delta t}{4}}RHS(Q_{1})}$

${\displaystyle {\dfrac {Q_{1}^{2}-Q^{2}}{\Delta t/2}}+{\dfrac {5}{6}}RHS(Q_{1}^{1})}$

${\displaystyle {\dfrac {Q^{3}-Q^{2}}{\Delta t}}+{\dfrac {3}{4}}RHS(Q_{1}^{2})}$

${\displaystyle \displaystyle Q^{i}}$ is the matrix of gross conservative data at the ith step and at the time t.

${\displaystyle Q_{1}^{i}}$ is the same matrix, at the time ${\displaystyle t+\Delta t}$.

${\displaystyle \displaystyle Q^{i}}$ and ${\displaystyle Q_{1}^{i}}$ are Np x Nvar matrices, with Np, number of points and Nvar, number of variables. Those matrices store the conservative raw variables.

Spatial discretization

Grid arrangements

H-Allegro uses a hybrid-colocated-staggered grid which allows a good accuracy and robustness of the code, while having an unambiguous definition of the boundary conditions.

For example, in 2D and for a rectangular structured grid, there are 4 grids, namely :

• a cell-centered grid, S, dedicated to the scalars ${\displaystyle \rho }$, ${\displaystyle \rho {E}}$ and ${\displaystyle \rho {Y}}$
• different face-centered grids corresponding to momentum (or vector) components. In this case, 2 grids exist : ${\displaystyle V_{i}}$ for ${\displaystyle \rho {U_{i}}}$ and ${\displaystyle V_{j}}$ for ${\displaystyle \rho {U_{j}}}$.
• a ${\displaystyle \tau _{ij}}$ grid dedicated to the boundary conditions and that minimises the number of interpolations.

Finite differences

The finite differences method was discovered by Taylor Brook. Adapted by Lagrange to the resolution of differential equations, it makes it possible to interpolate or differentiate a discrete function, with a chosen accuracy.

${\displaystyle f_{n+b_{k}}=f_{n}+b_{k}\Delta {x}\left({\dfrac {\partial {f}}{\partial {x}}}\right)_{n}+{\dfrac {(b_{k}\Delta {x})^{2}}{2!}}\left({\dfrac {{\partial }^{2}{f}}{\partial {x}^{2}}}\right)_{n}+...+{\dfrac {(b_{k}\Delta {x})^{m}}{m!}}\left({\dfrac {{\partial }^{m}{f}}{\partial {x}^{m}}}\right)_{n}+O(\Delta {x}^{m+1})}$

${\displaystyle b_{k}}$ is the step.

Actually, the accuracy of the scheme depends on the Stencil, which means it depends on the number of points taken into account in order to calculate the differential. A Stencil of 2 means that 2 points of each side are considered.

Taking into account the chosen numerical scheme and a symmetrical Stencil, the differential of a function at the point n is written :

${\displaystyle {\dot {f_{n}}}=\sum _{k=1}^{N}{{\dot {\beta }}_{k}}.{\dfrac {f_{n+b_{k}}}{{\Delta }x}}+O({\Delta }x^{N-1})}$

Where ${\displaystyle {\dot {\beta _{k}}}={\bar {\beta _{k}}}\left({\dfrac {1}{b_{k}}}-\sum _{l=1}^{N+1}{\dfrac {\bar {\beta _{l}}}{b_{l}}}\right)}$ and ${\displaystyle {\bar {\beta }}_{k}=\prod _{j=1etj{\neq }k}^{N}{\dfrac {b_{j}}{n_{j}-b_{k}}}}$

N : value of the Stencil. ${\displaystyle {\dot {\beta _{k}}}}$ : balancing ratio, giving more importance to the points near the considered one.

Formulae for the hybrid-colocated-staggered scheme :

• First point

${\displaystyle {\dot {f_{1}}}=-{\dfrac {46}{15}}{\dfrac {f_{1}}{\Delta {x}}}+{\dfrac {15}{4}}{\dfrac {f_{3/2}}{\Delta {x}}}-{\dfrac {5}{6}}{\dfrac {f_{5/2}}{\Delta {x}}}+{\dfrac {3}{20}}{\dfrac {f_{7/2}}{\Delta {x}}}}$

• General case :

${\displaystyle {\dot {f_{n}}}={\dfrac {75}{64}}{\dfrac {f_{n+1/2}-f_{n-1/2}}{\Delta {x}}}-{\dfrac {25}{384}}{\dfrac {f_{n+3/2}-f_{n-3/2}}{\Delta {x}}}+{\dfrac {3}{640}}{\dfrac {f_{n+5/2}-f_{n-5/2}}{\Delta {x}}}}$

Interpolation and differentiation operations

In order to calculate the various variables, the code uses interpolations and differentiations. Multidimensional arrangements are considered as combinations of monodimensional arrangements, which means that operators defined in 1D will be also used in 3D to do interpolation and differentiation operations.

• The interpolation is specially used in order to change the grid (for example, needing the value of the scalar to obtain the value of a component). There are 2 interpolations possible :
  • changing from the scalar grid to a momentum component (or ${\displaystyle \tau _{ij}}$) one, using the subroutine Interp1,
  • changing from a momentum component (or ${\displaystyle \tau _{ij}}$) grid to the scalar one, using the subroutine Interp2.

• The differentiation is used to differentiate, but it also changes the grid.
  • Diff1 differentiates and changes the grid from the scalar one to a momentum component one,
  • Diff2 differentiates and changes a momentum component grid to the scalar one.

Notes and References

1. ↑ Van Kalmthout, p.50 Theoretical and Numerical combustion - T. Poinsot and D. Veynante, edition of 2001
2. ↑ Using staggered grids with characteristic boundary conditions when solving compressible reactive Navier-Stokes equations - E. Albin, Y. D'Angelo and L. Vervisch
3. ↑ Minimal storage time-advancement schemes for spectral methods - A. Wray
4. ↑ Multicomponent Transport Algorithms - A. Ern and V. Giovangigli
5. ↑ http://www.idris.fr/su/Scalaire/babel/passage\_production.html