User talk:Ferrandm

Hypotheses

The hypotheses to simplify the Navier-Stokes equations are :

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

${\rho }{\dfrac {DE}{Dt}}={\dfrac {\partial {{\rho }E}}{\partial {t}}}+{\dfrac {\partial {{\rho }u_{i}E}}{\partial {x_{i}}}}={\dot {\omega _{T}}}+{\dfrac {\partial }{\partial {x_{i}}}}({\lambda }{\dfrac {\partial {T}}{\partial {x_{i}}}})-{\dfrac {\partial }{\partial {x_{i}}}}({\rho }{\sum }_{k=1}^{N}h_{s,k}Y_{k}V_{k,i})+{\dfrac {\partial {\sigma _{ij}}u_{i}}{\partial {x_{j}}}}+{\dot {Q}}+{\rho }{\sum }_{k-1}^{N}Y_{k}f_{k,i}(u_{i}+V_{k,i})$

The term ${\dfrac {\partial }{\partial {x_{i}}}}({\rho }{\sum }_{k=1}^{N}h_{s,k}Y_{k}V_{k,i})$ can be ignored with regard to ${\dot {\omega _{T}}}$, which is the heat release due to combustion.

This release can be defined by this equations : ${\dot {\omega _{T}}}=-\sum _{k=1}^{N}\Delta h_{f,k}^{0}{\dot {\omega _{k}}}=\nu _{F}M_{F}{\dot {\omega }}Q_{m}$, where $Q_{m}$ is the specific heat of reaction, $\Delta h_{f,k}^{0}$ is the standard enthalpy of formation and ${\dot {\omega _{k}}}$ is the reaction rate.

${\rho }{\sum }_{k=1}^{N}h_{s,k}Y_{k}V_{k,i}$ is the power produced by volume forces $f_{k}$ on species k.

${\dot {Q}}$ is the heat source term which is null here since there is no electric spark, laser or radiative flux.

$Y_{k}$ is the mass fraction of the species k.

The tensor $\sigma _{ij}$ is composed of the tensor of constraints and the tensor of pressure : $\sigma _{ij}=\tau _{ij}-P\delta _{ij}$

Thanks to these hypotheses, the total energy balance becomes :

${\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}$.

With $S_{5}$ which corresponds to the heat release due to combustion, ${\dot {\omega _{T}}}$.

• The term $S_{i}^{*}$ is not used. Thanks to it we can define other forces like electromagnetical force or weight.

• $S_{k}$ is either a fuel or a combustive source term.

$S_{6}=\nu _{F}M_{F}{\dot {\omega }}$ if it is a fuel source term or $S_{7}=\nu _{O}M_{O}{\dot {\omega }}$ if it is an oxidant source term.

$M_{F}$ and $M_{O}$ are respectively the molar mass of fuel and oxidant.

• ${\dot {\omega }}$ can be determined by the Arrhenius law :

${\dot {\omega }}=B_{1}(\phi )\rho Y_{F}Y_{O}e^{\dfrac {\beta }{\alpha }}e^{\dfrac {-T_{a}}{T}},B_{1}(\phi )=cst$

• 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, $\lambda$ is given by :

$\lambda ={\dfrac {\mu C_{p}}{P_{r}}}$

Where $P_{r}$ is the Prandtl number :

$P_{r}={\dfrac {\nu }{D_{th}}}$

This number compares the momemtum distribution with the heat distribution.

• The Lewis number for the species k, $Le_{k}$, is :

$Le_{k}={\dfrac {\lambda }{\rho C_{P}D_{k}}}$ $={\dfrac {D_{th}}{D_{k}}}$

Where $D_{k}$ is the species diffusion.

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

• The Schmidt number, $S_{ck}$, is :

$S_{ck}=P_{r}Le_{k}$

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

$D_{th}=D_{k}$ and $D_{k}={\dfrac {D}{\rho }}$

$P_{r}={\dfrac {\mu }{D}}$

Therefore, the diffusion coefficient can be defined as :

$D={\dfrac {\mu }{S_{ck}}}$

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

${\vec {\phi }}=-D.{\vec {\nabla }}C$

${\vec {\phi }}$ 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. Consequently as the gas is viscous, the compressible effects are not dominating. What involves that the bulk viscosity is useless.

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

$C_{P}={\dfrac {7}{2}}r$

$C_{V}={\dfrac {5}{2}}r$

• The gas respects the law of perfect gases :

$P=\rho rT$

• The combustion is an irreversible transformation whose the creation of entropy is compensated with an entropy given up 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 :

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

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

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

$\alpha _{s}=0.76$

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

$\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}$

The Kronecker symbol, $\delta _{ij}$, egals 1 if i=j, 0 else.

• The acoustic Reynolds number is defined such as :

$R_{c}={\dfrac {cL}{\nu }}$

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

Closure equations

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

• $P=\rho rT$

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

     $\alpha _{s}=0.76$

• $q_{j}=\lambda {\dfrac {\partial T}{\partial x_{j}}}$

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

• $\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}$

Parameters

Makefile

The makefile allows to compile routines and subroutines.

In order to compile, the command "make" will be typped in the terminal.

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

The configuration of compilation should be modified in function of the cluster considered (Antares, Jade, ...) :

• FC indicates what is the used compiler. For instance : if the language is fortran90 and MPI, it will be written FC = mpif90.
• A part of the code is written in C, it is about 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.

-p is an option of compilation. It is used in order to have gprof command.

Since various computer can be used, the compilation have to be conditional. That is to say that several configurations of compilation are available in function of the computer.

The common files of the executable for calculation and the common files of the executable 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 then 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.sh

This file concerns the creation and the copies of files in various directories during the code execution.

It is possible to activate or to deactivate the function "echo" thanks to the command "set -x" or "set +x".

It is necessary to remember the case dimension with the variable "NDIM".

The variable "NAME" gives the name of the executable which calculates the case. This executable will be launched thanks to the function "mpirun".

The total number of processors corresponds to NPX x NPY x NPZ.

For instance, if the name of the executable is Mezo3D and if the number of processors by direction is 8, then NAME = Mezo3D and the executable will be launched like this : mpirun -n 512 Mezo3D.

The file fort.10.sh exports the variables NX, NY, NZ, NPX, NPY, NPZ, dir and NP so that they can be used by other shell files.

Creating directories requires to use the subroutine generate(). Thus the file "launch-run" must be compilated with the option -g.

The subroutine generate() 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".

It is possible to submit all the jobs to the computer thanks to the subroutine submit(). So the file launch-run must be compilated with the option -s. A job is submit with the command "llsubmit job.ll".

The end of this file reminds to the user what is the dimension of the case and how many processors are used.

Vtk_bin_write

This file concerns the distribution of processors and the choice of the calculator.

If the calculator is ATON, then the value entered would be 1 and so, it would be 0 for the calculator VARGAS :

#define ATON 1
#define VARGAS 0

The distribution of processors can be defined in function of

• Coordinates :

 for (k = *dnz - 1; k < *fnz; k++)
 {  for (j = *dny - 1; j < *fny; j++)
    {  for (i = *dnx - 1; i < *fnx; i++)
       {  cbuf[0] = *(coorx + i);
          cbuf[1] = *(coory + j);
          cbuf[2] = *(coorz + k);
          FWRITE(cbuf, sizeof(double), 3, fptr);
       }
    }
 }

• Vectors : in that case it will depend on the dimension of the problem, NDIM. For instance, if NDIM = 3 :

else //Ndim=3
 {  for (k = *dnz - 1; k < *fnz; k++)
    {
       for (j = *dny - 1; j < *fny; j++)
       {
          for (i = *dnx - 1; i < *fnx; i++)
          {
             cbuf[0] = *(vectx + i + j*(*nnx) + k*(*nnx)*(*nny));
             cbuf[1] = *(vecty + i + j*(*nnx) + k*(*nnx)*(*nny));
             cbuf[2] = *(vectz + i + j*(*nnx) + k*(*nnx)*(*nny));
             FWRITE(cbuf, sizeof(double), 3, fptr);
          }
       }
    }
 }

• Scalar : this function will define the id number of the processor :

for (k = *dnz - 1; k < *fnz; k++)
    {
       for (j = *dny - 1; j < *fny; j++)
       {
          for (i = *dnx - 1; i < *fnx; i++)
          {
             buf = *(scalar + i + j*(*nnx) + k*(*nnx)*(*nny));
             FWRITE(&buf, sizeof(double),1, fptr);
           }
       }
   }

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 has to be created like this :

• Working name : locating more easily the calculations.

@ job_name = $NAME-$dir 

• 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. 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 really necessary memory by process has to get closer as much as possible to the wanted memory.

@ 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 send in case of overtaking of elapsed time.

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

Test case

Profiling

Profiling is a way to optimize the code skill. It consists of 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 CPU consummate time as well as the routines which are called from them.
• Flat profile : list of the routines in the decreasing order of the CPU consummate time and the number of calls.

Analysis of the flat profile

Description

• % time : percentage of total CPU time consumed by the routine.
• Cumalative 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 CPU time consumed. Thanks to the cumulative seconds, we can deduce that the total CPU time consumed egals 0.25 second. The number of seconds represented by the subroutine derivative2 on slice is about 0.05 second. This subroutine is invoked 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 does it spend.

The routines are indexed in decreasing order by comsumption of CPU time.

This report details the percentage of the total time that was spent in a function and its children. The total of this percentages is not equal to 100 %. For instance the % run time 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, of additional cost and of 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.

It is also necessary 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 the 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. Else, 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 way uses the subroutines MPI_ISSEND() and MPI_IRECV().

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 till 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 :

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

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

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

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

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

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

$Q^{i}$ and $Q_{1}^{i}$ are Np x Nvar matrix, with Np, number of points and Nvar, number of variables.