User talk:Bossenne

From www.coria-cfd.fr
Revision as of 14:07, 28 June 2012 by Bossenne (Talk | contribs) (3DPeriodic)

Jump to: navigation, search

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

Motivation

H-Allegro is a parallel, high accuracy, compressible, reactive, direct Navier Stokes solver (DNS). It is the MPI-parrallelized evolution of the former code ALLEGRO. From ALLEGRO, H-Allegro has mainly been developed by the PARALGO company, and two PhD students: Eric Albin and Marianne Sjostrand.


Equations and Hypotheses

Equations

Hypotheses

The hypotheses to simplify the Navier-Stokes equations are :

  • The total energy balance according to Poinsot-Veynante is described below, if  :

The term can be ignored with regard to , which is the heat release due to combustion.

This release can be defined by these equations : , where is the specific heat of reaction, is the standard enthalpy of formation and is the reaction rate.

is the power produced by volume forces on species k, considered as non-existent since there are no volume forces.

is the heat source term which is null here since there is no electric spark, laser or radiative flux.

is the mass fraction of the species k.

The tensor is composed of the tensor of constraints and the tensor of pressure :

Thanks to these hypotheses, the total energy balance becomes :

.


With corresponding to the heat release due to combustion, .

  • The term is not used. It usually defines other forces like electromagnetical forces or weight.
  • is either a fuel or a combustive source term.

if it is a fuel source term or if it is the oxidant one.

and are respectively the molar mass of fuel and oxidant.

  • can be determined by the Arrhenius law :

  • The heat flows due to mass fraction gradients and the species molecular distribution due to temperature gradients are neglicted (cf : Soret and Dufour).
  • The termal conductivity, is given by :

Where is the Prandtl number :

This number compares the momemtum distribution with the heat distribution.

  • The Lewis number for the species k, , is :

Where is the species diffusion.

We will consider the Lewis number egals to 1 in order to simplify the physics of the pre-mixture flame.

  • The Schmidt number, , is :

  • From the relation of these last three numbers, we can deduct that :

and

Therefore, the diffusion coefficient can be defined as :

  • We use the Law of Fick which is a simplified diffusion law :

is the flux of particles density.

C is the particles density.

Thus the diffusion coefficients of the various species are characterized by the number of Lewis.

  • We consider that the gas is a perfect, viscous, reagent and diatomic gas. As a consequence of the viscosity, the compressible effects are not dominating, which involves that the bulk viscosity is neglicted.
  • The gas is reagent thus the mixture of various species is not isothermal, they must be individually followed. It implies that the calorific capacities depend on the temperature and on the composition.
  • For a diatomic gas, the calorific capacities can be defined as follows :

  • The gas respects the law of perfect gases :

  • The combustion is an irreversible transformation whose creation of entropy is compensated with an entropy given by the system to the outside, because of a thermal transfer, thus we can consider that the transformation is isentropic. Consequently, the law of Laplace for the thermodynamics is applicable :

  • According to the power law, the dynamic viscosity depends only on the temperature :

  • The tensor of the constraints, by respecting the hypothesis of a Newtonian fluid, is :

The Kronecker symbol, , egals 1 if i=j, 0 else.

  • The acoustic Reynolds number is defined such as :

c is the speed of sound : 340 m., L is the characteristic length and is the kinematic viscosity of air : 1.45e-5 m²..

Closure equations

     


Models

  • Mixing : constant Schmidt number
  • Chemistry : one step Arrhenius' law
  • Viscosity : power law


Numerics

  • Spatial : 6th-order finite difference scheme
  • Temporal : 3rd order explicit time integration (Runge-Kutta3)
  • Chemical :
    • One step chemistry based on an Arrhenius' law
    • partially premixed gas
  • Boundary conditions : the 3D-NSCBC processing, applied in a referential attached to the local streamlines.


Data Structures

  • 1D, 2D, 3D structured solver
  • Optimized non-blocking MPI communications
  • Homogeneous mesh refinement
  • 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
Test image 2

Performances

Scalability Test

Partie privée -- User's manual

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 , and
  • different face-centered grids corresponding to momentum (or vector) components. In this case, 2 grids exist : for and for .
  • a 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.

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 :

avec et

N : value of the Stencil.  : balancing ratio, giving more importance to the points near the considered one.


Formulae for the hybrid-colocated-staggered scheme :

  • First point
  • General case :


Interpolation and differentiation operations

In order to calculate the various variable, 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 inerpolation 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 ) one, using the subroutine Interp1,
    • changing from a momentum component (or ) 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.

Parameters

fort.10/Mezo3D - Executable file

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 , and .
  • Npx, Npy, Npz : number of processors allocated to the resolution, respectively in the directions , and . The total number of points in the i direction is , and the total number of processors used for the calculation is .

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 , and .
  • IGridx, IGridy, IGridz : 1 if the grid is uniform in the considered direction, 0 otherwise.
  • Alphax, Alphy, Alphaz : stretching ratios.
  • Betax, Betay, Betaz : stretching rates.
  • XP0, YP0, ZP0 : stretching positions.

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

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


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 scalar, defined by the equations. For example, Nsc = 3 if the scalars are , and .
  • Re_ac : acoustic Reynolds of the flow.

c : sound speed ; L : characteristic length.

  • Gamma : ratio between heat capacity at constant pressure and heat capacity at constant volume.
  • Pr : Prandtl number of the flow.

Where 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 resolve, if 0, it won't.
  • CFL_react : CFL (Courant-Friedrichs-Lewy condition) of the reaction.
  • Alpha : thermal expansion ratio.

 : temperature of the burnt gases ;  : temperature of the unburnt gases.

  • Beta : Zeldovitch number.

 : temperature of the activation barrier so that the reaction can happen.

  • Damla_factor : Damköhler factor which defines the operational conditions of a reaction.
  • NuF*MF and NuO*MO : stoichiometric coefficient respectively for the fuel and the oxidant.

M : molar mass of the species.

  • YF_O and YO_O : stoichiometric coefficients in the reference flow of fuel and oxidant.


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 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.
  • NormDir : plane normal direction (1,2,3).
  • PlaneCoord : position of the plane on XP2(NormDir), between 1 and for the direction i.
  • REDIM : if 0, the results will be adimensional, if 1 they will be dimensional.
  • Lo* : reference length (1m).
  • Co* : reference sound speed (1m/s).
  • Rho_o* : reference density ().
  • To* : reference temperature (1K).


SCALAR

The value of the scalars are given here.

  • Uo : speed of the flow.
  • Cv : heat capacity at constant volume.
  • Rv :
  • Xo :
  • Yo :


Launch-run.sh

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 $@

In this example, the first calculation is in 3D, uses a 40x40x40 mesh and works on a 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.

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 [1] 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 is defined by the biggest planed configuration. From this, a number of cores for the reference execution can be specified :

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.

  • Speed-up :
  • Parallel efficiency :

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


Organization of H-Allegro

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

  1. 3DPeriodic le nom change, non ?
  2. Allegro


3DPeriodic

  • all_chemistry.f90  : cf all_chemistry.f90 in SOLV


  • all_init.f90  : initializes various parameters and fields.
    • LECTUS_SCAL : reads the initial radius of the flame and the filtration thickness of speed.
    • INIT_FIELDS : initialization various fields and determination of the speed of the premixed flame.
    • INIT_PROF : initializes the input conditions and the relaxation ratios.


  • all_solver.f90  : cf all_solver.f90 in SOLV


  • all_nscbcmezo.f90  : cf all_nscbc.f90 in SOLV, but this one has another subroutine.
    • Mask_mask :  ?


  • all_post.f90  : cf all_post.f90 in SOLV. There is one more subroutine is this version.
    • Q_critere :  ?


  • communs.f90  : indicates the value of some common parameters (wealth, initial speed...).
    • FTS :  ?
    • Mask_VISC :  ?


  • ouvre_file.f90  : reads the time and the scalars.

ALLEGRO

DER/DF6

The folder PADE is not described because the Pade scheme is not used.

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 one, while 2 means the contrary, going from the 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 neighbours A préciser !
    • echange_line :  ?
    • echange_surf :  ?
    • derivative1_on_slice, derivative2_on_slice : differentiation when the boundary conditions are symmetric A vérifier !
    • 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 A vérifier !
    • 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
    •  ?


  • 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

A faire

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 :  ?
    • ROTATIONNEL : calculates the curl along the number of dimensions.
    • PROJECT :  ?
    • DILUTION :  ?


  • 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 ().
    • 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 grid.
    • CONTINUITE : calculates (continuity equation) on the S grid.
    • ENERGY : calculates (conservation of energy) on the S grid.
    • MOMENTUM : calculates (conservation of momentum) and the stress tensor .
    • 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  : declarations (number of scalars and variables).
    •  ?


  • all_equ.f90  : calculates some variables.
    • RHOETHERM : calculates .
    • SOUND : calculates the sound speed .
    • 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 forced to 0, calculation of the transverse terms, viscosity and convection terms and ).
    • 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, , , Pr) and calculation of , 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 :  ?


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

Notes et Références

  1. http://www.idris.fr/su/Scalaire/babel/passage\_production.html