* fichier :  tubedeto3d2.dgibi
************************************************************************
* Section : Chimie Combustion
************************************************************************
****************************************************************
*****      PROPAGATION D UNE DETONATION DANS UN TUBE       *****
*****      MODELE COMBUSTION H2-air de PLEXUS              *****
*****      A. BECCANTINI, SFME/LTMF, 15.01.02              *****
****************************************************************
*
* The file is divided into 5 parts
*
* 1) mesh
* 2) initial conditions and gas properties
* 3) parameters used in the computation
* 4) the main part (where the computation is realised)
* 5) the post-treatment
*

   'OPTION' 'ECHO' 0 'DIME' 3
         'TRAC' 'X' ;
         
****************************************************************
****************************************************************
*****      1D mesh                                         *****
****************************************************************
****************************************************************
*
*
*  A6 ---- A5 ------------------- A4 ----------------------- A8
*   |      |                       |                         |
*   |  AR  |        ZONE1          |         ZONE2           |
*   |      |                       |                         |
*  A1 ---- A2 ------------------- A3 ----------------------- A7
*   |      |                       |                         |
*   |  L1  |          L2           |           L3            |
*   |<---->|<--------------------->|<----------------------->|
*
*  AR = activation region
*  ZONE1
*  ZONE2
*

 RAF = 1  ;
 
 NL2SL1 = 49 ;
 NL3SL1 = 100 ;
 L1 = 0.4 ;
 
 L2 = NL2SL1 '*' L1 ;
 L3 = NL3SL1 '*' L1 ;

* If  NL2SL1 = 49 and NL3SL1 = 50, then L2 = 19.6, L3 = 20.0 ;
 
 N1 = RAF ;
 N2 = NL2SL1 '*' N1 ;
 N3 = NL3SL1 '*' N1 ;
  
* For sake of simplicity, we will only consider one mesh in y-direction

 DX = L1 '/' N1 ;
 NY = 4 ;
 NZ = NY ;
 H1 = NY '*' DX ;

 A1 = 0.0 0.0 0.0 ;
 A2 = L1  0.0 0.0 ;
 A3 = (L1 '+' L2) 0.0 0.0 ;
 A4 = (L1 '+' L2) H1 0.0 ;
 A5 = L1  H1 0.0 ;
 A6 = 0.0 H1 0.0 ;
 A7 = (L1 '+' L2 '+' L3) 0.0 0.0 ;
 A8 = (L1 '+' L2 '+' L3) H1 0.0 ;

 'OPTION' 'ELEM' 'SEG2' ;
 A1A2 = A1 'DROIT' N1 A2 ;
 A2A3 = A2 'DROIT' N2 A3 ;
 A3A4 = A3 'DROIT' NY A4 ;
 A4A5 = A4 'DROIT' N2 A5 ;
 A5A2 = A5 'DROIT' NY A2 ;
 A2A5 = 'INVERSE' A5A2 ;
 A5A6 = A5 'DROIT' N1 A6 ;
 A6A1 = A6 'DROIT' NY A1 ;
 A3A7 = A3 'DROIT' N3 A7 ;
 A7A8 = A7 'DROIT' NY A8 ;
 A8A4 = A8 'DROIT' N3 A4 ;
 A4A3 = 'INVERSE' A3A4 ;

*
**** 2D mesh
*

* DOM1 is the activation region
* DOM2 is zone 1
* DOM3 is zone 2

 'OPTION' 'ELEM' 'QUA4' ;
 DOM1P = 'TRANSLATION'  A2A5 N1 ((-1.0 '*' L1) 0.0 0.0) ;
 DOM2P = 'TRANSLATION'  A2A5 N2 (L2 0.0 0.0) ;
 DOM3P = 'TRANSLATION'  A3A4 N3 (L3 0.0 0.0) ;

* DOMTOT = total region
 
 DOMTOTP = DOM1P 'ET' DOM2P 'ET' DOM3P ;
 'ELIMINATION' (1.0D-3 '/' RAF) DOMTOTP ;
 'ELIMINATION'  DOMTOTP (1.0D-3 '/' RAF) (A1A2 'ET' A2A3 'ET'
    A3A4 'ET' A4A5 'ET' A5A2 'ET' A5A6 'ET' A6A1 'ET' A3A7
    'ET' A7A8 'ET' A8A4 'ET' A4A3) ;

 'OPTION' 'ELEM' 'CUB8' ;
 DOM1 = 'VOLUME' DOM1P 'TRANSLATION' NZ (0.0 0.0 H1) ;
 DOM2 = 'VOLUME' DOM2P 'TRANSLATION' NZ (0.0 0.0 H1) ;
 DOM3 = 'VOLUME' DOM3P 'TRANSLATION' NZ (0.0 0.0 H1) ;
 
 DOMTOT = DOM1 'ET' DOM2 'ET' DOM3 ;
 'ELIMINATION' (1.0D-3 '/' RAF) DOMTOT ;

*
**** The tables 'DOMAINE'
*

 $DOMTOT = 'MODELISER' DOMTOT 'EULER' ;
 $DOM1   = 'MODELISER' DOM1 'EULER' ;
 $DOM2   = 'MODELISER' DOM2 'EULER' ;
 $DOM3   = 'MODELISER' DOM3 'EULER' ;
 
 TDOMTOT = 'DOMA' $DOMTOT 'VF' ;
 TDOM1   = 'DOMA' $DOM1   'VF' ;
 TDOM2   = 'DOMA' $DOM2   'VF' ;
 TDOM3   = 'DOMA' $DOM3   'VF' ;
 
 QDOMTOT = TDOMTOT  . 'QUAF' ;
 QDOM1   = TDOM1    . 'QUAF' ;
 QDOM2   = TDOM2    . 'QUAF' ;
 QDOM3   = TDOM3    . 'QUAF' ;

 'ELIMINATION' QDOMTOT (0.001 '/' RAF) QDOM1 ;
 'ELIMINATION' QDOMTOT (0.001 '/' RAF) QDOM2 ;
 'ELIMINATION' QDOMTOT (0.001 '/' RAF) QDOM3 ;
  
*
**** MOD1 = model (created just to display the CHAMELEMs)
*

 MOD1 =  'MODELISER' ('DOMA' $DOMTOT 'MAILLAGE') 'THERMIQUE' ;
 
*
**** Line for graphics
*

 'OPTION' 'ELEM' 'SEG2' ;
 P1 = (DX '/' 2) (DX '/' 2) (DX '/' 2) ;
 P2 = ((L1 '+' L2 '+' L3) '-' (DX '/' 2)) (DX '/' 2) (DX '/' 2) ;
 LIGEVO = (P1 'DROIT' (N1 '+' N2 '+' N3 '-' 1) P2) 'COULEUR' 'ROUG'  ;
 'ELIMINATION' LIGEVO ('DOMA' $DOMTOT 'CENTRE') (0.001 '/' RAF) ;

* 'OPTION' 'SAUVER' ('CHAINE' 'mesh.sauv_raf' RAF) ;
* 'SAUVER' RAF DOMTOT $DOMTOT $DOM1 $DOM2 $DOM3 LIGEVO MOD1 ;

****************************************************************
****************************************************************
*****      Initial conditions and gas properties.          *****
****************************************************************
****************************************************************

* 'OPTION' 'ECHO' 1 'TRAC' 'X' ;

* RAF = 1 ;
* 'OPTION' 'REST' ('CHAINE' 'mesh.sauv_raf' RAF) ;
* 'RESTITUER' ;

*
**** 0 from a numerical point of view
*

 zero = 1.0d-9 ;
 
*
**** (Non-Uniform) initial conditions.(DOM1 et DOM2)
*    We compute the conservative variables.
*
*
* Contact discontinuity between DOM2 and DOM3 (same pressure,
* different temperature)
*

 PINI = ('MANUEL' 'CHPO' (('DOMA' $DOM1 'CENTRE') 'ET'
        ('DOMA' $DOM2 'CENTRE')) 1 'SCAL' 99700 'NATU' 'DISCRET') 'ET'
        ('MANUEL' 'CHPO' (('DOMA' $DOM3 'CENTRE')) 1 'SCAL' 99700
         'NATU' 'DISCRET') ;
 TINI = ('MANUEL' 'CHPO' (('DOMA' $DOM1 'CENTRE') 'ET'
        ('DOMA' $DOM2 'CENTRE')) 1 'SCAL' 285. 'NATU' 'DISCRET') 'ET'
        ('MANUEL' 'CHPO' (('DOMA' $DOM3 'CENTRE')) 1 'SCAL' 400
         'NATU' 'DISCRET') ; 
 WINI =  'MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 3 'UX' 0.0 
         'UY' 0.0 'UZ' 0.0 ;
*
**** From initial conditions on molar fractions, we duduce
*    the data for the virtual mixture of burned/unburned
*    gases
         
*
**** Molar fractions (before combustion)
*

 XH21 = ('MANUEL' 'CHPO' (('DOMA' $DOM1 'CENTRE') 'ET'
        ('DOMA' $DOM2 'CENTRE')) 1 'H2' 0.3 'NATU' 'DISCRET') 'ET'
        ('MANUEL' 'CHPO' (('DOMA' $DOM3 'CENTRE')) 1 'H2' 0.2
         'NATU' 'DISCRET') ;
 XH2O1 = 'MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'H2O' zero
         'NATURE' 'DISCRET' ;
 XO21 = ('MANUEL' 'CHPO' (('DOMA' $DOM1 'CENTRE') 'ET'
        ('DOMA' $DOM2 'CENTRE')) 1 'O2' 0.147 'NATU' 'DISCRET') 'ET'
        ('MANUEL' 'CHPO' (('DOMA' $DOM3 'CENTRE')) 1 'O2' 0.2
         'NATU' 'DISCRET') ;
 XN = XH21 'ET' XH2O1 'ET' XO21 ;
 CHUN = ('MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'SCAL' 1.0D0) ;
 XN21 =  CHUN '-' ('PSCAL' XN CHUN  ('MOTS' 'H2' 'O2' 'H2O')
         ('MOTS' 'SCAL' 'SCAL' 'SCAL') ) ;
 XN21 = 'NOMC' 'N2' XN21 'NATURE' 'DISCRET' ;
 XN = XN 'ET' XN21 ;
 
*
***** Gas properties
*

 PGAZ = 'TABLE' ;

* Polynomial degree of specific heats
 
 PGAZ . 'NORD' = 4 ;

* Species explicitly treated in the Euler Equations

 PGAZ . 'ESPEULE' = 'MOTS' 'H2  ' 'O2  ' 'H2O ' ;

* Species non explicitly treated

 PGAZ . 'ESPNEULE' = 'N2  ';

* Single gases properties

 PGAZ .  'H2  ' = 'TABLE'  ;
 PGAZ .  'H2O ' = 'TABLE'  ;
 PGAZ .  'N2  ' = 'TABLE'  ;
 PGAZ .  'O2  ' = 'TABLE'  ;

* R (J/Kg/K)

 mH2 = 2. '*' 1.00797E-3 ;
 mo2 = 2. '*' 15.9994E-3 ;
 mH2O = mh2 '+' (0.5 '*' mo2) ;
 mN2 = 2 '*' 14.0067E-3 ;
 RGAS = 8.31441 ;
 
 PGAZ .  'H2  ' . 'R' = RGAS '/' mh2 ;
 PGAZ .  'H2O ' . 'R' = RGAS '/' mh2o ;
 PGAZ .  'N2  ' . 'R' = RGAS '/' mn2 ;
 PGAZ .  'O2  ' . 'R' = RGAS '/' mo2 ;

* Polynomials regressions coefficients

 PGAZ .  'H2  ' . 'A' = 'PROG'  9834.91866 0.54273926 0.000862203836
                               -2.37281455E-07 1.84701105E-11 ;
 PGAZ .  'H2O ' . 'A' = 'PROG' 1155.95625 0.768331151 -5.73129958E-05
                              -1.82753232E-08 2.44485692E-12 ;
 PGAZ .  'N2  ' . 'A' = 'PROG' 652.940766 0.288239099 -7.80442298E-05
                              8.78233606E-09 -3.05514485E-13 ;
 PGAZ .  'O2  ' . 'A' = 'PROG' 575.012333  0.350522002 -0.000128294865
                              2.33636971E-08 -1.53304905E-12;
                              
* Formation enthalpies at 0K (J/Kg))
*
*     h_i(0K)  = h_i(T0) '-' \int_0^{T0} cp_i(x) dx 
*              = h_i(T0) '-' (\int_0^{T0} cv_i(x) dx '+' R_i * T0)
*     

 PGAZ .  'H2  ' . 'H0K' = -4.195D6 ;
 PGAZ .  'H2O ' . 'H0K' = -1.395D7 ;
 PGAZ .  'N2  ' . 'H0K' = -2.953D5 ;
 PGAZ .  'O2  ' . 'H0K' = -2.634D5 ;


* 'OPTION' 'SAUVER' 'gas_properties.sauv' ;
* 'SAUVER' PGAZ ;

* 
**** We compute the mass fractions
*

 MMOL = ('MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'H2' Mh2
         'NATURE' 'DISCRET') 'ET'
        ('MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'O2' Mo2
         'NATURE' 'DISCRET') 'ET'
        ('MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'H2O' Mh2o
         'NATURE' 'DISCRET') 'ET'
        ('MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'N2' Mn2
         'NATURE' 'DISCRET') ;

 MTOT = 'PSCAL' MMOL XN ('MOTS' 'H2' 'O2' 'H2O' 'N2')
   ('MOTS' 'H2' 'O2' 'H2O' 'N2') ;
    
 YH21 = (XH21 '*' MH2) '/' MTOT ;
 YO21 = (XO21 '*' MO2) '/' MTOT ;
 YN21 = (XN21 '*' MN2) '/' MTOT ;
 YH2O1 = (XH2O1 '*' MH2O) '/' MTOT ;

*
**** The conservative variables
*

 RINI = (('NOMC' 'SCAL' yh21) '*' (PGAZ .'H2' . 'R')) '+'
        (('NOMC' 'SCAL' yo21) '*' (PGAZ .'O2' . 'R')) '+'
        (('NOMC' 'SCAL' yh2o1) '*' (PGAZ .'H2O' . 'R')) '+'
        (('NOMC' 'SCAL' yn21) '*' (PGAZ .'N2' . 'R'))
   ;
   
 RN = PINI '/' (RINI '*' TINI) ;

* Internal energy

 eini = 'MANUEL' 'CHPO' ('DOMA' $DOMTOT 'CENTRE') 1 'SCAL' 0.0 ;
 Tfois = Tini ;
 
 'REPETER' BL1 ((PGAZ . 'NORD') '+' 1) ;
 
    acel = (('NOMC' 'SCAL' yh21) '*' ('EXTRAIRE' (PGAZ .'H2  ' . 'A')
            &BL1)) '+'
           (('NOMC' 'SCAL' yo21) '*' ('EXTRAIRE' (PGAZ .'O2  ' . 'A')
            &BL1)) '+'
           (('NOMC' 'SCAL' yh2o1) '*' ('EXTRAIRE' (PGAZ .'H2O ' . 'A')
            &BL1)) '+'
           (('NOMC' 'SCAL' yn21) '*' ('EXTRAIRE' (PGAZ .'N2  ' . 'A')
            &BL1))  ;
            
    eini = eini '+' (acel '*' (Tfois '/' &BL1)) ;
    
    Tfois = Tfois * tini ;
    
 'FIN' BL1 ;

 RETN = RN '*' (eini '+' (0.5 '*' ('PSCAL' WINI WINI
   ('MOTS' 'UX' 'UY' 'UZ') ('MOTS' 'UX' 'UY' 'UZ')))) ;
 GN = RN '*' WINI ;
 YN = YH21 '+' YO21 '+' YH2O1 ;
 RYN = RN '*' YN ;

 'SI' VRAI ;
*   Test  
    VN P T Y GAMN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN ;
    PCEL = 'MAXIMUM' PINI ;
    TCEL = 'MAXIMUM' TINI ;
    'SI' (('MAXIMUM' ((P '-' pini) '/' PCEL) 'ABS')> 1.0D-3) ;
       'MESSAGE' 'Probleme calcul energie totale' ;
       'ERREUR' 21 ;
    'FINSI' ;
    'SI' (('MAXIMUM' ((T '-' Tini) '/' TCEL)) > 1.0D-3) ;
       'MESSAGE' 'Probleme calcul energie totale' ;
       'ERREUR' 21 ;
    'FINSI' ;
 'FINSI' ;

* 'OPTION' 'SAUV' 'ic.sauv';
* 'SAUVER' RN RYN GN RETN zero ;
* 'SAUVER' RAF DOMTOT $DOMTOT $DOM1 $DOM2 $DOM3 LIGEVO MOD1 ;


****************************************************************
****************************************************************
*****           Parameters for the computations           ******
****************************************************************
****************************************************************
 
* Upwind scheme
* METO = 'VLH' ;
 METO = 'SS' ;
 
* Iterations
* Final time
* Safety Factor fot the time step
* Second order reconstruction?

 NITER = 10000 ;
 TFINAL = 5.0D-3 ;
 SAFFAC = 0.7 ;
 LOGSO = VRAI ;

* 'OPTION' 'SAUVER' 'parameters.sauv' ;
* 'SAUVER'  METO  NITER  TFINAL  SAFFAC  LOGSO ;

****************************************************************
****************************************************************
*****           The computation                           ******
****************************************************************
****************************************************************

************************
**** PROCEDURE DETO ****
************************

 'DEBPROC'  DETO ;
 'ARGUMENT' PGAZ*'TABLE' TN*'CHPOINT' RYN*'CHPOINT' DTI*'FLOTTANT' ;

*
* dRYO2/dT  = MO2 * A * (T ** -B) EXP( - TA / T)
*             (RYH2 ** C) * (RYO2)
*             
**** Les parametres physiques
*
*    Ts = temperature de seuil
*

    TS = 'MANUEL' 'CHPO' ('EXTRAIRE' TN 'MAILLAGE') 1 'SCAL' 600 ;
    TA = 8310.;
    A = 1.1725D14;
    B = 0.91;
    CON = 2;

    H0H2  = PGAZ .  'H2  ' . 'H0K' ;
    H0O2  = PGAZ .  'O2  ' . 'H0K' ;
    H0H2O = PGAZ .  'H2O ' . 'H0K' ;

    RYH2 = 'EXCO' 'H2' RYN 'H2' 'NATU' 'DISCRET' ;

    RYO2 = 'EXCO' 'O2' RYN 'O2' 'NATU' 'DISCRET' ;

    RYH2O = 'EXCO' 'H2O' RYN 'H2O' 'NATU' 'DISCRET' ;

    RYH2 RYO2 RYH2O DELTAE = 'FLAM' 'ARRHENIU' TS A B CON TA
                            H0H2 H0O2 H0H2O DTI RYH2 RYO2 RYH2O
                            TN ;

    RY =  RYH2 'ET' RYO2 'ET' RYH2O ;

 'FINPROC' RY DELTAE ;
 
****************************
**** FIN PROCEDURE DETO ****
****************************


*
**** The table of gas properties 
*    (PGAZ)

* 'OPTI' 'REST' 'gas_properties.sauv' ;
* 'REST' ;

* 
**** The mesh
*    RAF  = refinement parameter
*    DX   = mesh dimension
*    DOM1 = activation region
*    DOM1 and DOM2 = zone 1
*    DOM3 = zone 2
*    DOMTOT = total region
*    LIGEVO (line for evolutions)
*    The  MMODEL objects
*       $DOMTOT
*       $DOM1
*       $DOM2
*       $DOM3
*    MOD1 = model object to plot 2D graphics

*
**** Initial conditions
*
*    RYN   : densities for H2, O2, H2O (kg/m^3)
*    RN    : total density (kg/m^3)
*    GN    : total momentum (kg/m^3 m/s)
*    RETN  : total energy (J/m^3)
*
*    zero  : "the numerical value of zero"
*
*

* 'OPTI' 'REST' 'ic.sauv' ;
* 'REST' ;

*
**** Parameters for the computation
*
*    METO   : upwind scheme
*    NITER  : max iterations number
*    TFINAL : final time
*    SAFFAC : safety factor for the time step
*    LOGSO  : second order reconstruction ?

*  'OPTION' 'REST' 'parameters.sauv' ;
*  'RESTITUER' ;
  

*
**** Chemical time step 
*
*    We compute d(RHY2)/dT at the von Neumann state
*
*    Molar reaction rate
*

 TVN = 2000 ;
 TA = 8310.;
 A = 1.1725D14;
 B = 0.91;
 CON = 2;

 H0H2  = PGAZ .  'H2  ' . 'H0K' ;
 H0O2  = PGAZ .  'O2  ' . 'H0K' ;
 H0H2O = PGAZ .  'H2O ' . 'H0K' ;

 RYNRED = 'REDU' RYN ('DOMA' $DOM2 'CENTRE') ; 
 RYH2 = 'MAXIMUM' ('EXCO' 'H2' RYNRED 'SCAL') ;
 RYO2 = 'MAXIMUM' ('EXCO' 'O2' RYNRED 'SCAL') ;
 
 OMEGA = A '*' (RYH2 '**' CON) '*' RYO2 '*' (TVN '**' (-1.0 '*' B))
   '*' ('EXP' (-1.0 '*' TA '*' (TVN '**' -1))) ;

* OMEGA in (mol/m3/s)
* OMEGAH2 in (kg/m3/s)

 OMEGAH2 =  OMEGA '*' 2 '*' 2E-3 ;
 DELTATC = RYH2 '/' OMEGAH2 ; 

* NB DELTATC << DELTAT_CON
 
*
**** Zone d'activation = etat AICC
*

 VN PN TN YN GAMN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN ;
 RYH2 = ('EXCO' 'H2' RYN 'H2' 'NATU' 'DISCRET') ;
 RYO2 = ('EXCO' 'O2' RYN 'O2' 'NATU' 'DISCRET') ;
 RYH2O = ('EXCO' 'H2O' RYN 'H2O' 'NATU' 'DISCRET') ;
 TN = TN '-' ('REDU' TN ('DOMA' $DOM1 'CENTRE')) '+'
      ('MANUEL' 'CHPO' ('DOMA' $DOM1  'CENTRE') 1 'SCAL' TVN) ;
 TS = 'MANUEL' 'CHPO' ('EXTRAIRE' TN 'MAILLAGE') 1 'SCAL' (0.9 '*' TVN)
    ;
 RYH2 RYO2 RYH2O DELTARE = 'FLAM' 'ARRHENIU' TS A B CON TA
                            H0H2 H0O2 H0H2O (1.0D5 '*' DELTATC) 
                            RYH2 RYO2 RYH2O TN ;

 RYN =  RYH2 'ET' RYO2 'ET' RYH2O ;
 RETN = RETN '+' DELTARE ;

 RN0 = 'COPIER' RN ;
 GN0 = 'COPIER' GN ;
 RETN0 = 'COPIER' RETN ;
 RYN0 = 'COPIER' RYN ;

*
**** Parameter for the time loop 
*

* Names of the residuum components

 LISTINCO = 'MOTS' 'RN' 'RUX' 'RUY' 'RUZ' 'RETN' 'H2' 'O2' 'H2O' ;
  
*
**** Geometrical coefficient to compute gradients
*

 GRADR CACCA COEFSCAL =  'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'NOLIMITE'
           ('MOTS' 'SCAL') RN ;
          
 GRADV CACCA COEFVECT =  'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'NOLIMITE'
           ('MOTS' 'UX' 'UY' 'UZ') GN ;

 TPS = 0.0 ;

*
**** The chemistry time step is too little
*    We take DT_CHEM > DELTATC 
*
*

 DT_CHEM = 50 '*' DELTATC ;

*
**** Temporal loop
*

 'MESSAGE' ;
 'MESSAGE' ('CHAINE' 'Methode = ' METO) ;
 'MESSAGE' ('CHAINE' 'SAFFAC      = ' SAFFAC) ;
 'MESSAGE' ;

  VN PN TN YN SN GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN RSN ;

  TNM1 = 'COPIER' TN ;

 'TEMPS' 'ZERO' ;
 'REPETER' BL1 NITER ;

*
**** Primitive variables
*
    
    VN PN TN YN GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN 
      TNM1  ;

    TNM1 = 'COPIER' TN ;

    'SI' LOGSO ;
    
       GRADR ALR = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
           ('MOTS' 'SCAL') RN  'GRADGEO' COEFSCAL ;

       GRADP ALP = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
           ('MOTS' 'SCAL') PN  'GRADGEO' COEFSCAL ;

       GRADV ALV = 'PENT' $DOMTOT 'CENTRE' 'EULEVECT' 'LIMITEUR'
          ('MOTS' 'UX' 'UY' 'UZ') VN   'GRADGEO'  COEFVECT ;
          
       GRADY ALY = 'PENT' $DOMTOT 'CENTRE' 'EULESCAL' 'LIMITEUR'
           ('MOTS' 'H2' 'O2' 'H2O') YN  'GRADGEO' COEFSCAL ;

       ROF VITF PF YF = 'PRET' 'PERFTEMP'  2 1
                            $DOMTOT  PGAZ
                            RN    GRADR  ALR
                            VN    GRADV  ALV
                            PN    GRADP  ALP
                            YN    GRADY  ALY ;

    'SINON' ;
                              
       ROF VITF PF YF  = 'PRET' 'PERFTEMP'  1 1 
                              $DOMTOT  PGAZ
                              RN   
                              VN   
                              PN   
                              YN    ;

    'FINSI' ;
                              
    RESIDU DELTAT =  'KONV' 'VF' 'PERFTEMP' 'RESI' METO
         $DOMTOT PGAZ LISTINCO  ROF VITF PF YF ;

    DT_CON = SAFFAC '*' DELTAT ;     
*
**** The time step linked to tps
*

    DTTPS = (TFINAL '-' TPS) * (1. '+' ZERO) ;
    
*
**** Total time step
*

    DTMIN = 'MINIMUM' ('PROG' DT_CON DTTPS DT_CHEM) ;
    
*
**** Increment of the variables (convection)
*

    RESIDU = DTMIN '*' RESIDU ;
    
    DRN = 'EXCO' 'RN' RESIDU 'SCAL' ;
    DGN = 'EXCO' ('MOTS' 'RUX' 'RUY' 'RUZ') RESIDU
       ('MOTS' 'UX' 'UY' 'UZ') ;
    DRETN = 'EXCO' 'RETN' RESIDU 'SCAL' ;
    DRYN  = 'EXCO' ('MOTS' 'H2' 'O2' 'H2O') RESIDU
            ('MOTS' 'H2' 'O2' 'H2O') ;
    
    
    TPS  = TPS '+' DTMIN ;
    RN   = RN '+' DRN ;
    GN   = GN '+' DGN ;
    RETN = RETN '+' DRETN ;
    RYN = RYN '+' DRYN ;
    
*
**** Increment of the variables (chemical source)
*

    VN PN TN YN GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN 
      TNM1  ;

    RYN DRETN = DETO PGAZ TN RYN DTMIN ;
    RETN = RETN '+' DRETN ;


    'SI' (((&BL1 '/' 20) '*' 20) 'EGA' &BL1) ;
        'MESSAGE' ('CHAINE' 'ITER =' &BL1 '  TPS =' TPS) ;
    'FINSI' ;

    'SI' (TPS > TFINAL) ;
       'QUITTER' BL1 ;
    'FINSI' ;

 'FIN' BL1 ;
 'TEMPS' ;

* 'OPTION' 'SAUVER' ('CHAINE' 'result.sauv_' RAF 'tps_' TPS ) ;
* 'SAUVER'  ;
 
****************************************************************
****************************************************************
*****           The post-treatment                        ******
****************************************************************
****************************************************************

* 'OPTI' 'REST' 'result.sauv_1tps_5' ;
* 'REST' ;

 GRAPH = VRAI ;
 GRAPH = FAUX ;

* 
**** The mesh
*

 'SI' GRAPH ;
    'TRACER' DOMTOT 'TITR' ('CHAINE' 'Maillage: '
       ('NBEL' DOMTOT) ' elts');
 'FINSI' ;

*
**** Initial conditions
*

  MOD1 =  'MODELISER'  ('DOMA' $DOMTOT 'MAILLAGE') 'THERMIQUE' ;
 
 'SI' GRAPH ;
    VN PN TN YN GAMMAN = 'PRIM' 'PERFTEMP'
       PGAZ RN0 GN0 RETN0 RYN0  ;
    CHM_RN  =  'KCHA'  $DOMTOT  'CHAM'  RN0  ;
    CHM_VN  =  'KCHA'  $DOMTOT  'CHAM'  VN  ;
    CHM_PN  =  'KCHA'  $DOMTOT  'CHAM'  PN  ;
    CHM_TN  =  'KCHA'  $DOMTOT  'CHAM'  TN  ;
    CHM_YN  =  'KCHA'  $DOMTOT  'CHAM'  YN  ;
    'TRAC' CHM_RN  MOD1 
       'TITR' ('CHAINE' 'rho at t=' 0.0)   ;
    'TRAC' CHM_VN  MOD1 
       'TITR' ('CHAINE' 'v at t= '  0.0)   ;
    'TRAC' CHM_PN  MOD1 
       'TITR' ('CHAINE' 'p at t= '  0.0)   ;
    'TRAC' CHM_TN  MOD1 
       'TITR' ('CHAINE' 'T at t= '  0.0)   ;
    'TRAC' CHM_YN  MOD1 
       'TITR' ('CHAINE' 'y at t= '  0.0) ;
 'FINSI' ;

*
**** Results
*

  VN PN TN YN  GAMMAN = 'PRIM' 'PERFTEMP' PGAZ RN GN RETN RYN ;
  
 'SI' GRAPH ;

    CHM_RN  =  'KCHA'  $DOMTOT  'CHAM'  RN  ;
    CHM_VN  =  'KCHA'  $DOMTOT  'CHAM'  VN  ;
    CHM_PN  =  'KCHA'  $DOMTOT  'CHAM'  PN  ;
    CHM_TN  =  'KCHA'  $DOMTOT  'CHAM'  TN  ;
    CHM_YN  =  'KCHA'  $DOMTOT  'CHAM'  YN  ;
    
    'TRAC' CHM_RN  MOD1 
       'TITR' ('CHAINE' 'rho at t=' TPS)   ;
    'TRAC' CHM_VN  MOD1 
       'TITR' ('CHAINE' 'v at t= '  TPS)   ;
    'TRAC' CHM_PN  MOD1 
       'TITR' ('CHAINE' 'p at t= '  TPS)   ;
    'TRAC' CHM_TN  MOD1 
       'TITR' ('CHAINE' 'T at t= '  TPS)   ;
    'TRAC' CHM_YN  MOD1 
       'TITR' ('CHAINE' 'y at t= '  TPS)   ;
 
 'FINSI' ;

*
*** Evolution Objects
*

 SS = 'COORDONNEE' 1 LIGEVO;

 evxx = 'EVOL' 'CHPO' ss LIGEVO ;
 lxx = 'EXTRAIRE'  evxx 'ORDO' ;

 x0 = 'MINIMUM' lxx;
 x1 = 'MAXIMUM' lxx ;

 titolo  = 'CHAINE'  ' at t = ' TPS '  Methode = ' METO
   ' SF = ' SAFFAC ;

 VNN = 'EXCO' 'UX' VN ;
 VNT = 'EXCO' 'UY' VN ;

* rho

 evro = 'EVOL'  'CHPO' RN LIGEVO ;
 lro  = 'EXTRAIRE'  evro 'ORDO' ;
 evro = 'EVOL' 'MANU' 'x' lxx 'ro' lro ;
 tro  = 'CHAINE' 'ro ' titolo ;

* u
      
 evu  = 'EVOL' 'CHPO' VNN LIGEVO ;
 lu   = 'EXTRAIRE'  evu 'ORDO' ;
 evu  = 'EVOL' 'MANU' 'x' lxx 'u' lu ;
 tu   = 'CHAINE' 'u ' titolo ;

* v
      
 evv  = 'EVOL' 'CHPO' VNT LIGEVO ;
 lv   = 'EXTRAIRE'  evv 'ORDO' ;
 evv  = 'EVOL' 'MANU' 'x' lxx 'v' lv ;
 tv   = 'CHAINE' 'v ' titolo ;
      

* p
 
 evp  = 'EVOL' 'CHPO' PN LIGEVO ;
 lp   = 'EXTRAIRE'  evp 'ORDO' ;
 evp  = 'EVOL' 'MANU' 'x' lxx 'p' lp ;
 tp   = 'CHAINE' 'p  ' titolo ;

 
* T
 
 evT  = 'EVOL' 'CHPO' TN LIGEVO ;
 lT   = 'EXTRAIRE'  evT 'ORDO' ;
 evT  = 'EVOL' 'MANU' 'x' lxx 'T' lT ;
 tT   = 'CHAINE' 'T  ' titolo ;

 'SI' GRAPH ;
    'DESSIN' evro 'TITRE' tro 'XBOR' x0 x1 'MIMA' ;
    'DESSIN' evp  'TITRE' tp 'XBOR' x0 x1  'MIMA' ;
    'DESSIN' evu 'TITRE'  tu 'XBOR' x0 x1  'MIMA' ;
    'DESSIN' evv 'TITRE'  tv 'XBOR' x0 x1  'MIMA' ;
    'DESSIN' evT  'TITRE' tT 'XBOR' x0 x1  'MIMA' ;
 'FINSI' ;


* Test

 PVN = 3.2E6 ;
 PCJ = 1.166 ;
 RHOVN = 3.7 ;
 RHOCJ = 1.4 ;

  
 PMAX = 'MAXIMUM' PN ;
 'SI' ((PMAX '>' PVN) 'OU' (PMAX '<' PCJ)) ;
    'MESSAGE' 'Pression = ???' ;
    'ERREUR' 5 ;
 'FINSI' ;
 
 RHOMAX = 'MAXIMUM' RN ;
 'SI' ((RHOMAX '>' RHOVN) 'OU' (RHOMAX '<' (0.85 '*' RHOCJ))) ;
    'MESSAGE' 'Densite = ???' ;
    'ERREUR' 5 ;
 'FINSI' ;
 
 'FIN' ;