Télécharger sine_bumpBM.dgibi

Retour à la liste

* fichier :  sine_bumpBM.dgibi
******************************************************************
* CALCUL DE L'ECOULEMENT SUBSONIQUE STATIONNAIRE DANS UN CANAL   *
* AVEC SINE-SHAPED BUMP                                          *
* FORMULATION VF COMPRESSIBLE EXPLICITE/IMPLICIT                 *
*                                                                *
* H. PAILLERE/P. GALON TTMF AOUT 1997                            *
*                                                                *
* MODIF BECCANTINI MARS 97                                       *
*                                                                *
*       BECCANTINI NOVEMBRE  98                                  *
* MODIF BECCANTINI AUGUST 2021                                   *
*                                                                *
*  Low Mach AUSM+ scheme is used for computation                 *
******************************************************************
 
 GRAPH = VRAI ;
 COMPLET = VRAI ;
 
 GRAPH = FAUX ;
 COMPLET = FAUX ; 
 
 'SI' complet ;
   RAF = 8 ;
   NITER = 20 ;
 'SINON' ;
   RAF = 5  ;
   NITER = 20 ;
 'FINSI' ;
 
 
 TYEL = 'QUA4' ;
 
 'OPTION'  'DIME' 2  'ELEM' TYEL  'ISOV' 'SULI'
   'ECHO' 1 'TRAC' 'X';
 
*
*** C.L. et initiales
*
 RO_INF = 1.4D0 ;
 P_INF  = 1.0D0 ;
 U_INF  = 0.0001D0 ;
 GAMSCAL    = 1.4D0;
 
************************************************************************
************************************************************************
***************** PARTIE PROCEDURES ************************************
************************************************************************
*
*
* Put it in utilproc
* Procedures will finish in FIN PARTIE PROCEDURES
*
 
***** $$$$ EXEXCH
*
*****************************************************
*****************************************************
** PROCEDURE EXEX POUR FORMULATION VF COMPRESSIBLE **
*****************************************************
*****************************************************
*
* We check wether the table RV is complete
*
*            
 'DEBPROC'  EXEXCH ;
 'ARGUMENT' RV*TABLE ;
 'MESSAGE' ;
 'MESSAGE' 'PROCEDURE EXEXCH' ;
*
**** Initialisation d'un CHPOINT 'ET' d'une MATRIK vide
*
 CHPVID MATVID = 'KOPS' 'MATRIK' ;
 RESPRO CHPVID MATVID ;
*
**** Les inconnues
*
 'SI' ('NON' ('EXISTE' RV  'UN')) ;
    'MESSAGE' ;
    'MESSAGE' 'UN = ???' ;
    'ERREUR' 21 ;
 'FINSI' ;
*
**** We check the compatibility between the variables and their names
*
 'SI' ('NON' ('EXISTE' RV  'LISTCONS')) ;
    'MESSAGE' ;
    'MESSAGE' 'LISTCONS = ???' ;
    'ERREUR' 21 ;
 'FINSI' ;
 'MESSAGE' ;
 'MESSAGE' 'Conservative variables check' ;
 UNCELL = 'EXCO' (RV . 'LISTCONS') (RV . 'UN') (RV . 'LISTCONS') ;
 ERRO = 'MAXIMUM' (UNCELL '-' (RV . 'UN')) 'ABS' ;
 ERRO = ERRO '/' (('MAXIMUM' UNCELL 'ABS') '+' 1.0D-15) ;
 'SI' (ERRO > 1.0D-6) ;
    'MESSAGE' ;
    'MESSAGE' 'UN = ???' ;
    'MESSAGE' 'LISTCONS = ???' ;
    'ERREUR' 21 ;
 'SINON' ;
    'MESSAGE' 'OK' ;
 'FINSI' ;
*
**** On which variables the error has to be computed?
*
 'SI' ('NON' ('EXISTE' RV  'ERROMOTS')) ;
    'MESSAGE' ;
    'MESSAGE' 'ERROMOTS = ???' ;
    'ERREUR' 21 ;
 'FINSI' ;
*
**** Definition of the table containing some results
*
 'SI' ('NON'  ('EXISTE' RV  'RESULTATS')) ;
   RV . 'RESULTATS' = 'TABLE' ;
 'FINSI' ;
* 
**** Existence d'une procédure pour imposer le conditions aux limites
*    à chaque iteration (interne ou externe)
*  
 'SI' ('NON' ('EXISTE' RV 'CLIM')) ;
    'MESSAGE' ;
    'MESSAGE' 'CLIM = ???' ;
    'ERREUR' 21 ;
 'FINSI' ;
* Par securité, on les impose tout de suite
 'SI' (RV . 'CLIM') ;
    'SI' ('NON' ('EXISTE' RV  'MAIFAN')) ;
       'MESSAGE' ;
       'MESSAGE' 'MAIFAN = ???' ;
       'ERREUR' 21 ;
    'FINSI' ;
    RV . 'UN' =  PROLIM RV ;
* DUCLIM à imposer dans l'inversion matricielle:
* pas d'increment sur les cellules fantomes
    DUCLIM = 0.0 '*' ('REDU' (RV . 'UN') (RV . 'MAIFAN')) ;
 'SINON' ;
    DUCLIM = 'COPIER' CHPVID ;
 'FINSI' ;
 
 RESPRO DUCLIM ;
*
* 
****** Mise a jour de la matrice à inverser pendant les iterations internes?
*
* RV . 'MATUPDAT' 
* RV . 'MUPINT'
*
 'SI' ('NON' ('EXISTE' RV 'MATUPDAT')) ;
    'MESSAGE' ;
    'MESSAGE' 'MATUPDAT = ???' ;
    'ERREUR' 21 ;
 'FINSI' ;
 'SI' (RV . 'MATUPDAT') ;
   'MESSAGE' ;
   'MESSAGE'
      'We always update the matrix to inverse'  ;
   LOGMAOLD = FAUX ;
 'SINON' ;
   'MESSAGE' ;
   'MESSAGE'
          'We don t always update the matrix to inverse' ;
   LOGMAOLD = VRAI ;
 'FINSI' ;
*
 RESPRO LOGMAOLD ; 
*
**** Iterations externes et/ou temps final
*    La table 'PASDETPS'
*
 'SI' ('NON'
    (('EXISTE' RV 'NITMAEX') 'OU' ('EXISTE' RV 'TFINAL'))) ;
    'MESSAGE' ;
    'MESSAGE'  'NITMAEX = ???' ;
    'MESSAGE'     'ou' ;
    'MESSAGE'  'TFINAL = ???' ;
    'ERREUR' 21 ;
 'FINSI' ;
*
**** Gas model
*
 'SI' ('NON' ('EXISTE' RV 'PGAZ')) ;
    'MESSAGE' ;
    'MESSAGE'  'PGAZ = ???' ;
    'ERREUR' 21 ;
 'FINSI' ;
*
**** BAS MACH
*
 'SI' ('NON' ('EXISTE' RV 'BASMACH')) ;
    'MESSAGE' ;
    'MESSAGE'  'BASMACH = ???' ;
    'ERREUR' 21 ;
 'SINON' ;
   'SI' (RV . 'BASMACH') ;
      'SI' (('NON' ('EXISTE' RV 'CO1')) 'OU'
            ('NON' ('EXISTE' RV 'CO2')))  ;
            'MESSAGE' ;
            'MESSAGE'  'CO1 = ???' ;
            'MESSAGE'  'CO2 = ???' ;
            'ERREUR' 21 ;
       'FINSI' ;
   'FINSI' ;
 'FINSI' ;   
*
*** La table PASDETPS
*    
 'SI' ('NON' ('EXISTE' RV 'PASDETPS')) ;
    'MESSAGE' ;
    'MESSAGE' 'We create the table PASDETPS' ;
    RV . 'PASDETPS' = 'TABLE' ;
 'SINON' ;
    'MESSAGE' ;
    'MESSAGE' 'We use the table PASDETPS that already exists' ;
 'FINSI' ;
  KTPS = RV . 'PASDETPS' ;
*
 RESPRO KTPS ;  
* Initialisation éventuelle de la table 
 'SI' (('NON' ('EXISTE' KTPS 'NUPASDT')) 'OU'
       ('NON' ('EXISTE' KTPS 'TPS')) 'OU'
       ('NON' ('EXISTE' KTPS 'TPSM'))) ;
    'MESSAGE' ;
    'MESSAGE' 'Table PASDETPS' ;
    'MESSAGE' 'NUPASDT = 0' ;
    'MESSAGE' 'TPS = 0.0' ;
    'MESSAGE' ;
    KTPS . 'NUPASDT' = 0 ;
    KTPS . 'TPSM' = 0.0D0 ;
    KTPS . 'TPS'  = 0.0D0 ;
 'FINSI' ;
* 
* KTPS contient
*
* KTPS . 'NUPASDT' = numero de pas de TPS actuel dans les itérations
* internes
*
* KTPS . 'TPSM'    = le TPS aprés (KTPS . 'NUPASDT'  '-' 2) itérations
* KTPS . 'TPS'     = le TPS aprés (KTPS . 'NUPASDT'  '-' 1) itérations
* KTPS . 'TPSP'    = le TPS aprés (KTPS . 'NUPASDT') itérations, i.e. à
* la fin de l'itération actuelle
*
*
*  
*         N.B. 'XINIT' is used in the iterative methods
*         We have decided to keep it 0
 RV . 'MATINV' . 'XINIT'  = 0.0 '*' (RV . 'UN') ;
 'FINPROC' ;
 
 
**** $$$$ EXEXIM
*
*****************************************************
*****************************************************
** PROCEDURE EXEX POUR FORMULATION VF COMPRESSIBLE **
*****************************************************
*****************************************************
*
* RV  . 'UN'      = les variables conservatives
*
* RV .  'ANOM'    = LOGIQUE (anomalie detectée ?)
*
* RV  . 'DTIMP'   = pas de tmps imposé (facultatif)
*
* RV  . 'LISTCONS'= noms des variables conservatives
*
* RV  . 'ERROMOTS'= noms des variables sur lesquelles on calcule Linf
*
* RV . 'CLIM'     = logique qui me dit si existe une procédure pour le
*                   calcul de conditions limites (procédure PROLIM)
* 
* RV . 'MAIFAN'   = (a définir si RV . 'CLIM')
*                   les cellules fantômes
*
* RV . 'TFINAL'   =  le temps final
*
* RV . 'NITMAEX'  =  le nombre d'itération externes
*
* N.B. Si RV . 'TFINAL' et RV . 'NITMAEX' sont les deux spécifiés, on
*      s'arrête quand un des deux critères est satisfait
*
* RV . 'MATUPDAT' = si VRAI, on ne met pas toujours à jour la matrice
*                   à inverser
*
* RV . 'MUPEXT'   = si (RV . 'MATUPDAT'), indice de mise à jours de la
*                   matrice jacobienne pendant les itérations externes
*
* RV . 'MUPINT'   = si (RV . 'MATUPDAT'), indice de mise à jours de la
*                   matrice jacobienne pendant les itérations internes
*
* RV . 'MUPINI'   = si (RV . 'MATUPDAT'), nombre initial de pas de temps 
*                   pendant lesquels on met à jour la matrice à chaque
*                   itération externe
*
* RV . 'MUPLIN'   = si (RV . 'MATUPDAT') et on utilise un solveur
*                   itératif pour la résolution du système linéaire,
*                   si le nombre d'itérations linéaires est > que
*                   RV . 'MUPLIN', alors on met à jour le
*                   preconditionneur
*
* RV . 'MATACC'   = si (RV . 'MATUPDAT'), logique qui nous dit si on
*                   veut utiliser la méthode de Broyden.
*
* RV . 'LISTOPER' = liste des opérateurs (ou des procédures) qui 
*                   interviennent dans le calcul (chaque opérateur a une 
*                   table associée), qui s'appelle &NOMPER ou
*                   & = position de l'opérateur dedans cette liste
*                   NOMPER = noms de l'opérateur ou de la procédure
*
* RV . 'NITMAIN' =  (à définir dans le cas d'un calcul implicite)
*                   le nombre max d'itération internes
* RV . 'NITMIIN' =  (à définir dans le cas d'un calcul implicite)
*                   le nombre min d'itération internes
*
* RV . 'EPSINT'   = (à définir dans le cas d'un calcul implicite)
*                   l'erreur pour le critère de convergence sur les
*                   itérations internes
*
* RV . 'MATHINV'  = (a définir sans la cas d'un calcul implicite)
*                   table de SOUSTYPE 'TYPINV' pour l'inversion de
*                   MATRIK (pour l'opérateur 'KRES')
*
* RV . 'RESULTATS' = table qui contient des resultats
*                    'NITERLIN' = nombre de iterations lineaires
*                    dans le solveur iteratif
*                    'ERROLIN'  = residu de l'erreur du solveur
*                    iteratif
*                    'ITERIN' = iteration interne (0, 1, ...)
*
*************************************************************************
* MODIF
*************************************************************************
* Message en anglais!
*            
 'DEBPROC'  EXEXIM ;
 'ARGUMENT' RV*TABLE ;
 'MESSAGE' ;
 'MESSAGE' 'PROCEDURE EXEXIM' ;
*
**** Initialisation and controls into exexch
*
 
 CHPVID MATVID DUCLIM LOGMAOLD KTPS = EXEXCH RV ;
 
*************************************************************************
*************************************************************************
*************************************************************************
*************************************************************************
**************   ITÉRATIONS EXTERNES ************************************
*************************************************************************
*************************************************************************
*************************************************************************
*************************************************************************
 LOGEXP = VRAI ;
* LOGEXP = variable logique qui me dit si on est en explicite ;
 LOGQIE = FAUX ;
*
**** Boucle qui s'arrête quand LOGQIE = VRAI ; i.e.
*
*     (KTPS . NUPASDT) = (RV . 'NITMAEX')
*      ou   
*     (KTPS . 'TPS') = (RV . 'TFINAL');
*
*
 'SI' ('EXISTE' RV 'DTIMP') ;
    ALPDT = RV . 'DTIMP' ;
 'SINON' ;
    ALPDT = 0.0 ;
 'FINSI' ;
 NITEX = 0 ;
 'REPETER'  BLEX  -1  ;
    NITEX = NITEX '+' 1 ;
*
****** Iteration interne
*
    RV . 'RESULTATS' . 'ITERIN' = 0 ; 
*
****** Mise à jour de la matrice à inverser
*
    'SI' ((NITEX 'EGA' 1) 'OU' (RV . 'ANOM')) ;
*         A la premiere iteration externe on doit calculer la matrice
       LOGUPEX = VRAI ;
    'SINON' ;
       'SI' LOGMAOLD ;
           MUPEXT = RV . 'MUPEXT' ;
*          On met a jour la matrice toutes les MUPEXT iterations externes
           LOGUPEX = ((NITEX '/' MUPEXT) '*' MUPEXT) 'EGA' NITEX ;
       'FINSI' ;
    'FINSI' ;
*
****** Trop d'iterations externes?
*       
    KTPS . 'NUPASDT' = KTPS . 'NUPASDT' '+' 1 ;
    'SI' ('EXISTE' RV 'NITMAEX') ;
       'SI' ( (KTPS . 'NUPASDT') '>EG' (RV . 'NITMAEX')) ;
          LOGQIE = VRAI ;
       'FINSI' ;
    'FINSI' ;
*
****** Impression
*
    'MESSAGE' ;
    'MESSAGE' ('CHAINE' 'PASDETPS = ' (KTPS . 'NUPASDT')
              '  TPS = ' (KTPS . 'TPS') '  DTM1 = ' ALPDT) ;
*
*** ALGORITHM A RESOUDRE
*
*   F(U^{n+1}) = 0
*
*   if (RV . INST) then
*
*     F(U^{n+1}) = -1. * (U^{n+1} - U^{n}) '/' (\delta t) 
*               '+' \sum_k  Res_k(U^{n+1}) '+' \sum_kexpl Res_kexpl(U^{n})
*
*   else
*
*     F(U^{n+1}) = -\sum_k Res_k(U^{n+1})
*
*   endif
*
*    
*   Avec une methode de type Newton-Raphson
*
*   a) U^{n+1,0} = U^{n}
*
*   b) for l=0,1,2,...
*
*         MAT(U^{n+1,l}) \delta U^{n+1,l} = -F(U^{n+1,l})
*
*         U^{n+1,l+1}  = U^{n+1,l} '+' \delta U^{n+1,l} 
*
*         if || \delta U^{n+1,l} ||_{\inf} < \epsilon goto c)
*
*   c)
*
*         U^{n+1} = U^{n+1,l+1}
*
*      endif
*
*** LES VARIABLES
*
*   RESIMP = \sum_k Res_k(U^{n+1,l})
*   RESEXP = \sum_kexpl Res_kexpl(U^{n})
*   MATASS = - \sum_k JACR_k(U^{n+1,l})
*   UNM = conserved variables at t^{n}
*   UNM = conserved variables at t^{n-1}
*   LREEDT = LISTREEL de (\delta t)_k t.q.
*            \delta t = min (\delta t)_k 
* 
    RESIMP = 'COPIER' CHPVID ;
    RESEXP = 'COPIER' CHPVID ;
    'SI' ('EXISTE' RV  'UNM2') ;
       RV . 'UNM3'    = 'COPIER' (RV . 'UNM2') ;
    'FINSI' ;
    'SI' ('EXISTE' RV  'UNM') ;
       RV . 'UNM2'    = 'COPIER' (RV . 'UNM') ;
    'FINSI' ;
    RV . 'UNM'  = 'COPIER' (RV . 'UN') ;
    MATASS = 'KOPS' 'MULT' 0.0  MATVID ;
    LREEDT = 'PROG' ;
*    
***********************************************
********* Boucle sur les operateurs ***********
***********************************************
*** On calcule: LREEDT
*               RESIMP
*               RESEXP
*              (MATASS)
*
    'REPETER'   BLOP  ('DIME'  (RV . 'LISTOPER')) ;
       NOMPER  = 'EXTRAIRE'  &BLOP  (RV . 'LISTOPER') ;
       NOTABLE = 'MOT'  ('TEXTE'  ('CHAINE'  &BLOP  NOMPER) ) ;
       'SI' (RV . NOTABLE . 'IMPL') ;
          'SI' ('EGA' NITEX 1) ;
             'MESSAGE' ;
             'MESSAGE' ('CHAINE' ('MOT' NOMPER) ': implicit') ;
          'FINSI' ;
*         JACO    = objet de type MATRIK (éventuellement vide)
*         RESIDU  =     "         RESIDU          "           
*         ALPHADT =     "         LISTREEL        "
          LOGEXP = FAUX ;
          JACO RESIDU ALPHADT = ('TEXTE'  NOMPER) (RV . NOTABLE) ;
          'SI' (RV . NOTABLE . 'ANOM') ;
             RV . 'ANOM' = VRAI ;
             RV . NOTABLE . 'ANOM' = FAUX ;
             'QUITTER' BLOP ;
          'SINON' ;
             MATASS = MATASS 'ET' ('KOPS' 'MULT' -1.0D0 JACO ) ;
             RESIMP = RESIMP 'ET' RESIDU ;
          'FINSI' ;
       'SINON' ;
          'SI' ('EGA' NITEX 1) ;
             'MESSAGE' ;
             'MESSAGE' ('CHAINE' ('MOT' NOMPER) ': explicit') ;
          'FINSI' ;
          RESIDU ALPHADT = ('TEXTE'  NOMPER) (RV . NOTABLE) ;
          'SI' (RV . NOTABLE . 'ANOM') ;
             RV . 'ANOM' = VRAI ;
             'QUITTER' BLOP ;
          'SINON' ;
             RESEXP = RESEXP 'ET' RESIDU ;
          'FINSI' ;
       'FINSI' ;
       LREEDT = LREEDT 'ET' ALPHADT ;
    'FIN'  BLOP ;
***********************************************
***** Fin  de la boucle sur les operateurs ****
***********************************************
*
******* Anomalie detectée (CTRL+S 9999)
*       On arrete ici
*
    'SI' ('NON' (RV . 'ANOM')) ;
*
******* On controlle la compatibilité E/S (si NITEX = 1) 
*
       'SI' ('EGA' NITEX 1) ;
*
**** Dans le cas implicite, on verifie l'existence
*    des parametres pour les itérations internes
*
          'SI' ('NON' LOGEXP) ;
             'SI' ('NON' (('EXISTE' RV  'NITMAIN') 'ET'
                ('EXISTE' RV  'EPSINT') 'ET' ('EXISTE' RV  'NITMIIN')));
                'MESSAGE' ;
                'MESSAGE' 'NITMAIN = ??? ' ;
                'MESSAGE' 'NITMIIN = ??? ' ;
                'MESSAGE' 'EPSINT  = ??? ' ;
                'ERREUR' 21 ;
             'FINSI' ;
          'FINSI' ;
       'FINSI' ;
*
******* Fin contrôle compatibilité E/S
*
*** Mise a jour de la table RV . 'PASDETPS' au debu du calcul
*
       'SI' ('EXISTE' RV 'DTIMP') ;
          ALPDT = RV . 'DTIMP' ;
       'SINON' ;
          ALPDT = 'MINIMUM' LREEDT ;
       'FINSI' ;
       ALPDT1 =  (RV . 'TFINAL') '-' (KTPS . 'TPS') ;
       'SI' ((ALPDT 'EGA' ALPDT1 (ALPDT '*' 1.0D-6)) 'OU'
             (ALPDT > ALPDT1)) ;
          KTPS . 'TPSP' =  (RV . 'TFINAL') ;
          ALPDT = ALPDT1 ;
          LOGQIE = VRAI ;
*         Dans ce cas, il faut metre a jour la matrice à inverser
*         car elle peut etre tres differnte pas rapport à celle
*         calculé precedentment
          LOGMAOLD = FAUX ;          
       'SINO' ;
          KTPS . 'TPSP' = (KTPS . 'TPS') '+' ALPDT ;
       'FINSI' ;
*
       'SI' LOGEXP ;
*************************************************************************
********** Cas explicite pour les variables conservatives ***************
*************************************************************************
          'SI' ((NITEX 'EGA' 1) 'OU' ((RV . 'ORDTPS') 'EGA' 1)) ;
             RV . 'UN' = (RV . 'UNM') '+' (ALPDT '*' RESEXP) ;
          'SINON' ;
             RV . 'UN' = ((4. '/' 3.) '*' (RV . 'UNM')) '+'
                         ((-1. '/' 3.) '*' (RV . 'UNM2')) '+'
                         (((2. * ALPDT) '/' 3.) '*' RESEXP) ;
          'FINSI' ;
          'SI' (RV . 'CLIM') ;
*             'MESSAGE' 'PROLIM after an explicit iteration' ;
             RV . 'UN' = PROLIM RV ;
          'FINSI' ;
       'SINON' ;
*************************************************************************
********** Cas implicite pour les variables conservatives ***************
*************************************************************************
*
****** Initialisation of RV . 'RESULTATS' . 'ERRONLIN' ;
*
          RV . 'RESULTATS' . 'ERRONLIN' = 'PROG' ;
*
*************************************************************************
************** Les iterations internes **********************************
*************************************************************************
*
          'REPETER' BLINT (RV . 'NITMAIN') ;
*
             RV . 'RESULTATS' . 'ITERIN' = &BLINT ;          
*
****** Implicite
*
             'SI' ((NITEX 'EGA' 1) 'OU' ((RV . 'ORDTPS') 'EGA' 1)) ;
*
*************** Implicit Euler
*             
                MAT1 = 'KOPS'  'MULT' (1. '/' ALPDT)
                        (RV . 'MATIDE') ;
                SMB = ((RV . 'UNM') '-' (RV . 'UN')) '/' ALPDT ;
             'SINON' ;
*
**************** BDF2
*             
                 MAT1 = 'KOPS'  'MULT' (1.5 '/' ALPDT)
                    (RV . 'MATIDE') ;
                 SMB = ((4. '*' (RV . 'UNM'))
                    '+' (-3. '*' (RV . 'UN'))
                    '+' (-1. '*' (RV . 'UNM2'))) '/' (2. '*' ALPDT) ;
             'FINSI' ;
*
********* Bas Mach
*
             'SI' (RV . 'BASMACH') ;
                MAT1 = MAT1 'ET' (RV . 'GAMSDTAU') ;
             'FINSI' ;
*          
             MATTOT = MAT1 'ET' MATASS ;
             RESTOT = SMB 'ET' RESIMP 'ET' RESEXP ;
*
********* LOGOLD = VRAI
*         Si methode directe, on utilise la meme parametrisation LU
*         Si methode iterative, on utilise le meme preconditionneur
*
             LOGOLD = VRAI ;
             'SI' LOGMAOLD ;
                'SI' (LOGUPEX 'ET' (&BLINT 'EGA' 1)) ;
*               LOGUPEX = VRAI -> At this external iteration we can
*               update the matrix to inverse
                   'OUBLIER' MATOLD ;
                   'MENAGE' ;
                   MATOLD = MATTOT ;
                   LOGOLD = FAUX ;
                   'MESSAGE' ;
                   'MESSAGE' 'We update the matrix' ;
                'SINON' ;
*               We don't care about LOGUPEX if (&BLINT > 1)
*               In this case we update the matrix each MUPINT-th
*               iteration
                'SI'
                    (((&BLINT  '/' (RV . 'MUPINT')) '*' (RV . 'MUPINT'))
                      'EGA' &BLINT) ;
                     'OUBLIER' MATOLD ;
                     'MENAGE' ;
                      MATOLD = MATTOT ;
                      LOGOLD = FAUX ;
                      'MESSAGE' ;
                      'MESSAGE' 'We update the matrix' ;
                   'FINSI' ;
                'FINSI' ;
             'SINON' ;
*
*            We always update the matrix to inverse
*
                'OUBLIER' MATOLD ;
                'MENAGE' ;
                LOGOLD = FAUX ;
                MATOLD = MATTOT ;
             'FINSI' ;
*
             'SI' ('EXISTE' (RV . 'MATINV')  'CONVINV') ;
                'SI' LOGOLD ;
                   NLIT = ('DIME' (RV . 'MATINV' . 'CONVINV')) ;
                   'SI' (NLIT > RV . 'MUPLIN') ;
                      'MESSAGE' ;
                      'MESSAGE' 'We update the matrix' ;
                      'OUBLIER' MATOLD ;
                      'MENAGE' ;
                      MATOLD = MATTOT ;
                   'FINSI' ;
                'FINSI' ;
             'FINSI' ;
             RV . 'MATINV' . 'MATASS' = MATOLD ;
             RV . 'MATINV' . 'MAPREC' = MATOLD ;
*
             'SI' LOGOLD ;
*
*            LOGOLD vrai si on utilise la meme
*            parametrisation LU ou le meme preconditionneur
*            pour calculer la solution du systeme lineaire
*          
                'SI' ((RV . 'MATINV' . 'TYPINV') 'EGA' 1) ;
*               Meme parametrisation LU que MATOLD
                   DELTAU = 'KRES' MATOLD 
                            'TYPI' (RV . 'MATINV')
                            'CLIM' DUCLIM
                            'SMBR' RESTOT
                            'IMPR' 0 ;
                'SINON' ;
*                  Meme preconditionneur que MATOLD
                   DELTAU = 'KRES' MATTOT
                            'TYPI' (RV . 'MATINV')
                            'CLIM' DUCLIM
                            'SMBR' RESTOT
                            'IMPR' 0 ;
                'FINSI' ;
             'SINON' ;          
                DELTAU = 'KRES' MATOLD 
                         'TYPI' (RV . 'MATINV')
                         'CLIM' DUCLIM
                         'SMBR' RESTOT
                         'IMPR' 0 ;
             'FINSI' ;
 
             RV . 'UN' = (RV . 'UN') '+' DELTAU ;
 
             'SI' (RV . 'CLIM') ;
                UNCELL = RV . 'UN' ;
*            'PROLIM after the matrix inversion' ;
                RV . 'UN' = PROLIM RV ;
*            We redefine  DELTAU to compute the error
                DELTAU = ((RV . 'UN') '-' UNCELL) '+' DELTAU ;
             'FINSI' ;
*
********* Boucle sur les operateurs implicites pour calculer RESIMP
*         et MATASS de RV . 'UN'
 
             RESIMP = 'COPIER' CHPVID ;
             MATASS = 'KOPS' 'MULT' 0.0 MATVID ;
             'REPETER'   BLOP  ('DIME'  (RV . 'LISTOPER')) ;
                NOMPER  = 'EXTRAIRE'  &BLOP  (RV . 'LISTOPER') ;
                NOTABLE = 'MOT'  ('TEXTE'  ('CHAINE'  &BLOP  NOMPER) ) ;
                'SI' (RV . NOTABLE . 'IMPL') ;
                   JACO RESIDU ALPHADT =
                      ('TEXTE'  NOMPER) (RV . NOTABLE) ;
                   'SI' (RV . NOTABLE . 'ANOM') ;
                      RV . 'ANOM' = VRAI ;
                      'QUITTER' BLINT ;
                   'SINON' ;
                      MATASS = MATASS 'ET' ('KOPS' 'MULT' 
                             -1.0D0 JACO ) ;
                      RESIMP = RESIMP 'ET' RESIDU ;
                   'FINSI' ;
                'FINSI' ;
             'FIN'  BLOP ;
*
********* Test de convergence
*
             ERRO = 'MAXIMUM' DELTAU 'ABS' 'AVEC' (RV . 'ERROMOTS')  ;
             'SI' (&blint 'EGA' 1) ;
                ERRO0 = ERRO ;
             'FINSI' ;
             RV . 'RESULTATS' . 'ERRONLIN' = (RV . 'RESULTATS'
             . 'ERRONLIN') 'ET' ('PROG' ERRO) ;
*
********* Impression
*                    
             'MESSAGE' ;
             'MESSAGE' ('CHAINE' 'PASDETPS = ' (KTPS . 'NUPASDT')
                   '  TPS = ' (KTPS . 'TPS') '  DT = ' ALPDT) ;
             'MESSAGE' ('CHAINE' 'ITERIN = ' &BLINT) ;
             'REPETER' BLERRO ('DIME' (RV . 'LISTCONS')) ;
                MOTERR = 'EXTRAIRE' (RV . 'LISTCONS') &BLERRO ;
                ERRBLE = 'MAXIMUM' DELTAU 'ABS' 'AVEC'
                             ('MOTS' MOTERR) ;
                'MESSAGE' (CHAIN 'ERREUR on ' MOTERR ' = ' ERRBLE) ;
             'FIN' BLERRO ;
             'SI' (('EXISTE' (RV . 'MATINV') 'CONVINV')) ;
                NLIT = ('DIME' (RV . 'MATINV' . 'CONVINV')) ;
                'MESSAGE' ('CHAINE' 'Linear iterations = ' NLIT );
                RV . 'RESULTATS' . 'NITERLIN' = NLIT ;
                REIT = 'EXTR' (RV . 'MATINV' . 'CONVINV') NLIT ;
                'MESSAGE' ('CHAINE' 'Linear residuum = ' REIT );
                RV . 'RESULTATS' . 'ERROLIN' = REIT ;
             'FINSI' ;
*
************ Criteres de sortie
*
             'SI' (&BLINT > (RV . 'NITMIIN')) ;
                'SI' ('EXISTE' RV 'EPSREL') ;
                   'SI' (ERRO < ((RV . 'EPSREL') '*' ERRO0)) ;
                      'QUITTER' BLINT ;
                   'FINSI' ;
                'FINSI' ;
                'SI' (ERRO '<' (RV . 'EPSINT')) ;
                   'QUITTER' BLINT ;
                'FINSI' ;
             'FINSI' ;
*             'SI' (ERRO > (1.0D3 '*' ERRO0)) ;
*                RV . 'ANOM' = VRAI ;
*                'QUITTER' BLINT ;
*             'FINSI' ;
*             'SI' (&BLINT '>EG' (RV . 'NITMAIN')) ;
*                RV . 'ANOM' = VRAI ;
*                'QUITTER' BLINT ;
*             'FINSI' ;
          'FIN' BLINT ;
*************************************************************************
**************** Fin boucle iterations internes *************************
*************************************************************************
***********************************************
************* Fin cas implicite ***************
***********************************************
       'FINSI' ;
*   ('SI' (RV . 'ANOM') ;)
    'FINSI' ;
*
    'SI' (RV . 'ANOM') ;
       LOGQIE = FAUX ;
       NITEX = NITEX '-' 1 ;
       KTPS . 'NUPASDT' = KTPS . 'NUPASDT' '-' 1 ;
       KTPS . 'TPSP' = KTPS . 'TPS' ;
       RV . 'UN' = 'COPIER' RV . 'UNM' ;
       'SI' ('EXISTE' RV  'UNM2') ;
          RV . 'UNM' = 'COPIER' (RV . 'UNM2') ;
       'FINSI' ;
    'SINON' ;
*
*** Mise a jour de la table RV . 'PASDETPS' à la fin du calcul
*
       KTPS . 'TPSM' =  KTPS . 'TPS' ;
       KTPS . 'TPS'  =  KTPS . 'TPSP' ;
    'FINSI' ;
*
**** Dernier commande de BLEX
*    
    'SI' LOGQIE ;
       'QUITTER' BLEX ;
    'FINSI' ;
 
 'FIN' BLEX ;
 
 'FINPROC' ;
 
*****************************************************
*****************************************************
** FIN PROCEDURE EXEX                              **
*****************************************************
*****************************************************
 
*****$$$$ PKON
 
*****************************************************
*****************************************************
** PROCÉDURE PKON                                  **
*****************************************************
*****************************************************
* Il faut définir:
*
* *KKONV . 'EQEX'     = table générale, qui contient toutes les 
*                       informations sur le calcul qu'on va faire.
*                       Dans ce table, on ne prend que:
*                       - la table (KKONV . 'EQEX' . 'PASDETPS')
*                       - les inconnues du problème
*                         (KKONV . 'EQEX' . 'UN')
*                       - le maillage fantome
*                         (KKONV . 'EQEX' . 'MAIFAN')
*                       - les vitesses de cut-off pour le bas Mach
*                         (KKONV . 'EQEX' . 'CO1')
*                         (KKONV . 'EQEX' . 'CO2')
*
* *KKONV . 'GAZ'      = le modelé de gaz qu'on considère
*                       - si KKONV . 'GAZ' = 'PERFMONO', alors
*                         on considère un gaz parfait mono-espèce
*                         polytropique
*
* *KKONV . 'MODELE'    = objet modele
*
* *KKONV . 'LISTCONS' = les noms des variables conservatives, i.e.
*                       densité, q.d.m., energie totale par unité de
*                       volume
*
* *KKONV . 'IMPL'     = calcul implicite ou non ?
*
* *KKONV . 'METHODE'  = méthode pour le calcul du flux convectif
*
* *KKONV . 'TYPEJACO' = 'VLH'
*                       'AUSMPLUS'
*                       'AUSMPLM'
*                       (à donner si (KKONV . 'IMPL'))
*
* *KKONV . 'ORDREESP' = ordre en espace (1 ou 2) ;
*
* *KKONV . 'GRADRN'   coeff. pour calculer le gradient de la densité
* *KKONV . 'GRADVN'   coeff. pour calculer le gradient de la vitesse
* *KKONV . 'GRADPN'   coeff. pour calculer le gradient de la pression
* *KKONV . 'VLIM'     vitesses au bord imposées
*                     (à donner si (KKONV . 'ORDREESP') = 2)
*
* *KKONV . 'LIMITEUR' = limiteur utilisé
*                       (à donner si (KKONV . 'ORDREESP') = 2)
*
* *KKONV . 'NFROZLIM' = entier
*                       On gele les limiteurs quand
*                       RV . 'PASDETPS' . 'NUPASDT' = KKONV . 'NFROZLIM'
*                       On le met dedans
*                       *KKONV . 'FROZALR'
*                       *KKONV . 'FROZALP'
*                       *KKONV . 'FROZALV'
*
* *KKONV . 'ALPHA'    = coefficient de securité pour le quel on 
*                       multiplie le pas de tps determiné par un 
*                       condition de type CFL.
*                       Il peut etre: 
*                       - un flottant
*                       - une mot que veut 'INF'
*
* *KKON . 'FACELIM' = maillage de centres de face ou on ne calcule pas
*                     le flux convectif et le jacobien
*
* *KKON . 'CLIM'    = LOGIQUE
*                     si vrai on doit definir une procedure PKOCLI
*
* *KKON . 'DT'      = output. Donne les pas de temps plus petit
*                     (r_elem/(u+c)) 
*
 'DEBPROC'  PKON ;
 'ARGUMENT' KKONV * TABLE ;
*
 CHPVID MATVID = 'KOPS' 'MATRIK' ; 
*    
 METO = KKONV . 'METHODE' ;
*
 RV = KKONV . 'EQEX' ;
*
 LOGIMP = KKONV . 'IMPL';
 LOGEXP = 'NON' LOGIMP ;
* 
 ORDESP = KKONV . 'ORDREESP' ;
*
 'SI' (('NEG' (KKONV . 'ALPHA') 'INF') 'ET'
    ('NEG' ('TYPE' (KKONV . 'ALPHA')) FLOTTANT)) ;
    'MESSAGE' 'PKON . ALPHA = ???' ;
    'ERREUR' 21 ;
 'FINSI' ;
*
 UN = 'REDU' (RV . 'UN') ('DOMA' (KKONV . 'MODELE') 'CENTRE') ;
*
 'SI' ('NEG' (KKONV . 'GAZ') 'PERFMONO') ;
*
******** EULER, monoespece, "calorically perfect" (cv = constant)
*
    'MESSAGE' ;
    'MESSAGE' 'PKON . GAZ = ???' ;
    'ERREUR' 21 ;
 'FINSI' ;
*
**** Type de jacobien (Bas Mach ?)
*
 'SI' LOGIMP ;
    ITJACO = 0 ;
    'SI' (('EGA' (KKONV . 'TYPEJACO') 'AUSMPLM') 'OU'
          ('EGA' (KKONV . 'TYPEJACO') 'RUSANOLM') 'OU'
          ('EGA' (KKONV . 'TYPEJACO') 'ROELM') 'OU'
          ('EGA' (KKONV . 'TYPEJACO') 'HLLCLM')) ;
       ITJACO = -1 ;
    'FINSI' ;
 'FINSI' ;
* 
 ITMETO = 0 ;
 'SI' (('EGA' METO 'AUSMPLM') 'OU'
       ('EGA' METO 'RUSANOLM') 'OU'
       ('EGA' METO 'ROELM') 'OU'
       ('EGA' METO 'HLLCLM')) ;
    ITMETO = -1 ;
 'FINSI' ;
*
***** Les variables conservatives
*
 MOT1 = 'EXTRAIRE' (KKONV . 'LISTCONS') 1 ;
 'SI' ('EGA' ('VALE' 'DIME') 2) ;
    NOMMOM = 'EXTRAIRE' (KKONV . 'LISTCONS')
       ('LECT' 2 3 ) ;
    NOMVEL = 'MOTS' 'UX' 'UY' ;
       MOT2 = 'EXTRAIRE' (KKONV . 'LISTCONS') 4 ;
    NOMGP = 'MOTS' 'P1DX' 'P1DY' ;
 'SINON' ;
    NOMMOM = 'EXTRAIRE' (KKONV . 'LISTCONS')
      ('LECT' 2 3 4) ;
    NOMVEL = 'MOTS' 'UX' 'UY' 'UZ' ;
    MOT2 = 'EXTRAIRE' (KKONV . 'LISTCONS') 5 ;
    NOMGP = 'MOTS' 'P1DX' 'P1DY' 'P1DZ' ;
 'FINSI' ;
 RN = 'EXCO' MOT1 UN 'SCAL' ;
 GN = 'EXCO' NOMMOM UN NOMVEL ;
 RETN = 'EXCO' MOT2 UN 'SCAL' ;
*
***** On calcule les variables primitive
*
 GAMN = (RV . 'PGAZ' . 'GAMN') ;
 VN PN = 'PRIM' 'PERFMONO' RN GN RETN  GAMN ;
*  VN PN = PPRIM RN GN RETN  GAMN ;
*  'SI' ((('MINIMUM' PN) < 0) 'OU' (('MINIMUM' RN) < 0)) ;
*     'MESSAGE' 'Negative density or pressure' ;
*      KKONV . 'ANOM' = VRAI ;
*     RN = ('ABS' RN) '+' (('MAXIMUM' RN 'ABS') '*' 1.001) ;
*     PN = ('ABS' PN) '+' (('MAXIMUM' PN 'ABS') '*' 1.001) ;
*  'FINSI' ;
*
**** Conditions aux limits ???
*
 'SI' (KKONV . 'CLIM') ;
 
*     TABLIM = 'TABLE' definie dans le programme principal ;
 
     KKONV . 'TABLIM' . 'RN' = RN ;
     KKONV . 'TABLIM' . 'VN' = VN ;
     KKONV . 'TABLIM' . 'PN' = PN ;
     KKONV . 'TABLIM' . 'GAMN' = GAMN ;
 
     CHPLIM RESLIM JACLIM = PKOLIM (KKONV . 'TABLIM') ;
     MAILLIM = ('EXTRAIRE' CHPLIM 'MAILLAGE') 'ET' (KKONV . 'FACELIM') ;
 'SINON' ;
     MAILLIM = KKONV . 'FACELIM' ;
     RESLIM JACLIM = 'KOPS' 'MATRIK' ;
     CHPLIM = 'COPIER' RESLIM ;
 'FINSI' ;
  MAILLIM = 'CHANGER' 'POI1' MAILLIM ;
*
**** La gravite
*
*
  RESGRA = 'FIMP' 'VF' 'GRAVMONO' 'RESI' (KKONV . 'LISTCONS')
               RN GN (KKONV . 'GRAVITE') ;
*               
  JACGRA = 'FIMP' 'VF' 'GRAVMONO' 'JACOCONS' (KKONV . 'LISTCONS')
               RN GN (KKONV . 'GRAVITE') ;
*
* Fin contribution gravité
*
 
*
***** On calcule les variables aux faces
*
 ORDTPS = 1 ;
 
*
 'SI' (ORDESP 'EGA'  1) ;
*       
********* Ordre 1 en espace
*
    ROF VITF PF GAMF  =  'PRET' 'PERFMONO' 1 ORDTPS
             (KKONV . 'MODELE') RN VN PN GAMN ;
 'SINON' ;
*       
********* Ordre 2 en espace
*
    'SI'  (KKONV . 'CLIM') ;
       GRADV ALV = 'PENT' (KKONV . 'MODELE') 'CENTRE' 'EULEVECT'
             (KKONV . 'LIMITEUR')  NOMVEL VN  'CLIM'
             ((KKONV . 'VLIM') '+' ('EXCO' NOMVEL CHPLIM))
             'GRADGEO' (KKONV . 'GRADVN')  ;
       GRADR ALR = 'PENT' (KKONV . 'MODELE') 'CENTRE' 'EULESCAL'
             (KKONV . 'LIMITEUR')  ('MOTS' 'SCAL') RN  
             'CLIM' ('EXCO' 'RN' CHPLIM) 'GRADGEO' (KKONV . 'GRADRN') ;
*
       GRADP ALP = 'PENT' (KKONV . 'MODELE') 'CENTRE' 'EULESCAL'
             (KKONV . 'LIMITEUR') ('MOTS' 'SCAL') PN  
             'CLIM' ('EXCO' 'PN' CHPLIM) 'GRADGEO' (KKONV . 'GRADPN') ;
*             
    'SINON' ;
       GRADV ALV = 'PENT' (KKONV . 'MODELE') 'CENTRE' 'EULEVECT'
              (KKONV . 'LIMITEUR')  NOMVEL VN  'CLIM'
              (KKONV . 'VLIM') 
              'GRADGEO' (KKONV . 'GRADVN')  ;
       GRADR ALR = 'PENT' (KKONV . 'MODELE') 'CENTRE' 'EULESCAL'
             (KKONV . 'LIMITEUR') ('MOTS' 'SCAL')  RN  
             'GRADGEO' (KKONV . 'GRADRN') ;
*
       GRADP ALP = 'PENT' (KKONV . 'MODELE') 'CENTRE' 'EULESCAL'
             (KKONV . 'LIMITEUR') ('MOTS' 'SCAL') PN  
             'GRADGEO' (KKONV . 'GRADPN') ;
    'FINSI' ;
*
****** Limiters = 0 on ghostcells
*                
    'SI' (('EXISTE' RV 'MAIFAN')) ;
       ALRLEV =  ('REDU' ALR (RV  . 'MAIFAN')) ;
       ALPLEV =  ('REDU' ALP (RV  . 'MAIFAN')) ;
       ALVLEV =  ('REDU' ALV (RV  . 'MAIFAN')) ;
       CELL = ALR ;
       ALR = ALR '-' ALRLEV ;
       'DETR' CELL ;
       CELL = ALP ;
       ALP = ALP '-' ALPLEV ;
       'DETR' CELL ;
       CELL = ALV ;
       ALV = ALV '-' ALVLEV ;
       'DETR' CELL ;
    'FINSI' ;
*
****** Frozen Limiters
*
*      External iterations
*
    'SI' ('EXISTE' KKONV 'NFROZLIM') ;
       ITEREX = RV . 'PASDETPS' . 'NUPASDT' ;
       'SI' (ITEREX 'EGA' (KKONV . 'NFROZLIM')) ;
          'MESSAGE' ;
          'MESSAGE' 'External iterations' ;
          'MESSAGE' 'We froze the limiters' ;
          KKONV . 'FROZALR' = 'COPIER' ALR ;
          KKONV . 'FROZALP' = 'COPIER' ALP ;
          KKONV . 'FROZALV' = 'COPIER' ALV ;
       'SINON' ;
          'SI' (ITEREX > (KKONV . 'NFROZLIM')) ;
*
*            Min
*          
             DAR = (ALR '-' (KKONV . 'FROZALR')) 'ABS' ;
             DAP = (ALP '-' (KKONV . 'FROZALP')) 'ABS' ;
             DAV = (ALV '-' (KKONV . 'FROZALV')) 'ABS' ;
             ALR = (ALR '+' (KKONV . 'FROZALR')) '-' DAR ;
             ALP = (ALP '+' (KKONV . 'FROZALP')) '-' DAP ;
             ALV = (ALV '+' (KKONV . 'FROZALV')) '-' DAV ;
             ALR = ALR '/' 2.0 ;
             ALV = ALV '/' 2.0 ;
             ALP = ALP '/' 2.0 ;
             KKONV . 'FROZALR' = 'COPIER' ALR ;
             KKONV . 'FROZALP' = 'COPIER' ALP ;
             KKONV . 'FROZALV' = 'COPIER' ALV ;
          'FINSI' ;
       'FINSI' ;
    'FINSI' ;
*
*      Internal iterations
*
    'SI' ('EXISTE' KKONV 'NFROZLII') ;
       ITERIN = RV . 'RESULTATS' . 'ITERIN' ;
       'SI' (ITERIN 'EGA' (KKONV . 'NFROZLII')) ;
          'MESSAGE' ;
          'MESSAGE' 'Internal iterations' ;
          'MESSAGE' 'We froze the limiters' ;
          KKONV . 'FROZALR' = 'COPIER' ALR ;
          KKONV . 'FROZALP' = 'COPIER' ALP ;
          KKONV . 'FROZALV' = 'COPIER' ALV ;
       'SINON' ;
          'SI' (ITERIN > (KKONV . 'NFROZLII')) ;
             DAR = (ALR '-' (KKONV . 'FROZALR')) 'ABS' ;
             DAP = (ALP '-' (KKONV . 'FROZALP')) 'ABS' ;
             DAV = (ALV '-' (KKONV . 'FROZALV')) 'ABS' ;
             ALR = (ALR '+' (KKONV . 'FROZALR')) '-' DAR ;
             ALP = (ALP '+' (KKONV . 'FROZALP')) '-' DAP ;
             ALV = (ALV '+' (KKONV . 'FROZALV')) '-' DAV ;
             ALR = ALR '/' 2.0 ;
             ALV = ALV '/' 2.0 ;
             ALP = ALP '/' 2.0 ;
             KKONV . 'FROZALR' = 'COPIER' ALR ;
             KKONV . 'FROZALP' = 'COPIER' ALP ;
             KKONV . 'FROZALV' = 'COPIER' ALV ;
          'FINSI' ;
       'FINSI' ;
    'FINSI' ;
*
 
*
********* Ordre 1 en temps
*
     ROF VITF PF GAMF   =  'PRET' 'PERFMONO'  ORDESP ORDTPS
                           (KKONV . 'MODELE')
                           RN  GRADR  ALR
                           VN  GRADV  ALV
                           PN  GRADP  ALP
                           GAMN ;
 'FINSI' ;
 
*
 'SI' (ITMETO 'EGA' (-1)) ;
     RESIDU DELTAT  =  'KONV' 'VF' 'PERFMONO' 'RESI' METO
                       (KKONV . 'LISTCONS') (KKONV . 'MODELE')
                       ROF VITF PF GAMF
                       MAILLIM (RV . 'CO1')
                       (RV . 'CO2') ;
 'SINON' ;
     RESIDU DELTAT  =  'KONV' 'VF' 'PERFMONO' 'RESI' METO
                       (KKONV . 'LISTCONS') (KKONV . 'MODELE')
                       ROF VITF PF GAMF
                       MAILLIM ;
 'FINSI' ;
* 
 'SI' ('EGA' ('TYPE' (KKONV . 'ALPHA')) FLOTTANT) ;
    ALPDT = 'PROG' ((KKONV . 'ALPHA') '*' DELTAT) ;
 'SINON' ;
    ALPDT = 'PROG' ;
 'FINSI' ;
 
 KKONV . 'DT' = DELTAT ;
 
*
**** Low Mach
*
 'SI' (RV . 'BASMACH') ;
     MATBM = 'KONV' 'VF' 'PERFMONO' 'GAMMCONS' (KKONV . 'LISTCONS')
             ('DOMA' (KKONV . 'MODELE') 'CENTRE')
             ('DOMA' (KKONV . 'MODELE') 'DIAMIN')
              RN VN PN GAMN (RV . 'CO1')
                  (RV . 'CO2') ;
     RV . 'GAMSDTAU' = 'KOPS' 'MULT'
         (1. '/' (KKONV . 'CFLDTAU')) MATBM ;
 'FINSI' ;
*
 'SI' LOGIMP ;
   'SI' (ITJACO 'EGA' 0) ;
       JACO =  'KONV' 'VF' 'PERFMONO' 'JACOCONS'  (KKONV . 'MODELE')
                  (KKONV . 'LISTCONS')
                  MAILLIM (KKONV . 'TYPEJACO')
                  RN VN PN GAMN ;
   'SINON' ;
          JACO =  'KONV' 'VF' 'PERFMONO' 'JACOCONS'  (KKONV . 'MODELE')
                  (KKONV . 'LISTCONS') 
                  MAILLIM  (KKONV . 'TYPEJACO')
                  RN VN PN GAMN (RV . 'CO1')
                  (RV . 'CO2') ;
   'FINSI' ;
   'RESPRO' (JACO 'ET' JACLIM 'ET' JACGRA)
            (RESIDU '+' RESLIM '+' RESGRA)
            ALPDT ;
 'SINON' ;
   'RESPRO' (RESIDU '+' RESLIM '+' RESGRA) ALPDT ;
 'FINSI' ;
*    
**** On detrui les choses qui ne servent plus
*
 'DETR' UN ;
 'OUBL' UN ;
 'DETR' RN ;
 'DETR' GN ; 
 'DETR' RETN ;
 'DETR' VN ;
 'DETR' PN ;
 'OUBL' RN ; 
 'OUBL' GN ;
 'OUBL' RETN ;
 'OUBL'  VN ;
 'OUBL' PN ;
*
**** Les MCHAML faces
*
 'DETR' ROF ;
 'DETR' VITF ;
 'DETR' PF ; 
 'DETR' GAMF ;
 'OUBL' ROF ;
 'OUBL' VITF ;
 'OUBL' PF ;
 'OUBL' GAMF ;
*
**** Les pentes et les limiteurs
*
 'SI' (ORDESP 'EGA' 2);
*
       'DETR' GRADR ;
       'DETR' GRADP ;
       'DETR' GRADV ;
       'DETR' ALR ;
       'DETR' ALP ;
       'DETR' ALV;
*
       'OUBL' GRADR ;
       'OUBL' GRADP ;
       'OUBL' GRADV ;
       'OUBL' ALR ;
       'OUBL' ALP ;
       'OUBL' ALV;
*
 'FINSI' ;
 
 'FINPROC' ;
*****************************************************
*****************************************************
** FIN PROCEDURE PKON                              **
*****************************************************
*****************************************************
 
 
*********************************************************************
**** Procedure PKOLIM ***********************************************
*********************************************************************
 
 
'DEBPROC' PKOLIM ;
'ARGUMENT' TABLIM*'TABLE' ;
 
 LISTP = 'MOTS' 'RN' 'UX' 'UY' 'PN' ;
 
 RCHLIM1 RCHRES1 = 'KONV' 'VF' 'PERFMONO' 'CLIM'  'RESI'
        MDOMINT MLIGG LISTCONS LISTP
        (TABLIM . 'RN') (TABLIM . 'VN') (TABLIM . 'PN')
        (TABLIM . 'GAMN') (TABLIM . 'CHINRI') 'INRI' ;
 
 RJACO1 = 'KONV' 'VF' 'PERFMONO' 'CLIM'  'JACOCONS'
        MDOMINT MLIGG LISTCONS LISTP
        (TABLIM . 'RN') (TABLIM . 'VN') (TABLIM . 'PN')
        (TABLIM . 'GAMN') (TABLIM . 'CHINRI') 'INRI' ;
 
 RCHLIM2 RCHRES2 = 'KONV' 'VF' 'PERFMONO' 'CLIM'  'RESI'
        MDOMINT MLIGD LISTCONS LISTP
        (TABLIM . 'RN') (TABLIM . 'VN') (TABLIM . 'PN')
        (TABLIM . 'GAMN') (TABLIM . 'CHOUTR') 'OUTRI'  ;
 
 RJACO2 = 'KONV' 'VF' 'PERFMONO' 'CLIM'  'JACOCONS'
        MDOMINT MLIGD LISTCONS LISTP
        (TABLIM . 'RN') (TABLIM . 'VN') (TABLIM . 'PN')
        (TABLIM . 'GAMN') (TABLIM . 'CHOUTR') 'OUTRI' ;
 
 
 CHPLIM = RCHLIM1 '+' RCHLIM2 ;
 RESLIM = RCHRES1 '+' RCHRES2 ;
 JACLIM = RJACO1 'ET' RJACO2 ;
 
 'RESPRO' CHPLIM RESLIM JACLIM ;
 'FINPROC' ;
 
*********************************************************************
**** Procedure CALC *************************************************
*********************************************************************
*
* Personal procedure
*
 'DEBP' CALC ;
 'ARGU' RVX*'TABLE' ;
*
 RV  = RVX . 'EQEX' ;
 UN = RV . 'UN' ;
 MCALC = RVX . 'MODELE' ;
*
* During the external iteration (except the first; in which
* UNM = UN of the previous step )
*
 'SI' ('NEG' 0 (RV . 'RESULTATS' . 'ITERIN')) ;
 
    RN = 'REDU' ('DOMA' MCALC 'CENTRE')
                ('EXCO' UN NOMDEN ('MOTS' 'SCAL')) ;
    GN = 'REDU' ('DOMA' MCALC 'CENTRE')
             ('EXCO' NOMMOM UN NOMVEL) ;
    RETN = 'REDU' ('DOMA' MCALC 'CENTRE')
                ('EXCO' NOMRET UN ('MOTS' 'SCAL'))  ;
    VN PN = 'PRIM' 'PERFMONO' RN GN RETN (RV . 'PGAZ' . 'GAMN') ;
*    
    'SI' ('EGA' (RVX . 'COMPT') 0) ;
       RVX . 'IT' = 'PROG' ;
       RVX . 'ER' = 'PROG' ;
    'SINON' ;
       CHPVOL = 'DOMA' MCALC  'VOLUME' ;
       VOLTOT = 'XTY' CHPVOL ('MANUEL' CHPO
      ('DOMA' MCALC 'CENTRE') 1 'SCAL' 1.) ('MOTS' 'SCAL')
       ('MOTS' 'SCAL');
       ERRO = 'ABS' (PN '-' (RVX . 'PN0')) ;
       ERRO =  'XTY' ERRO CHPVOL ('MOTS' 'SCAL')
       ('MOTS' 'SCAL') ;
       ERRO = ERRO '/' VOLTOT ;
       RVX . 'IT' = (RVX . 'IT') 'ET' ('PROG'
          (RVX . 'COMPT')) ;
       RVX . 'ER' = (RVX . 'ER') 'ET' ('PROG' ERRO) ;
    'FINSI' ;
    RVX . 'PN0' = PN ;
    RVX . 'COMPT' = (RVX . 'COMPT') '+' 1 ;
 'FINSI' ;
*
* We change the cut-off during the internal iterations
*
 'SI' FAUX ;
 'SI' ((RV . 'RESULTATS' . 'ITERIN') > 0) ;
     CUMAX = RV . 'CO1MAX' ;
     CUMIN = RV . 'CO1MIN' ;
     DCU = (CUMIN '-' CUMAX) '/' 5  ;
     CU = CUMAX '+' (DCU * (RV . 'RESULTATS' . 'ITERIN')) ;
*    Take the max between CU and CUMIN
     A = CUMIN '+' CU ;
     B = 'ABS' (CUMIN '-' CU) ;
     CU = 0.5 '*' (A '+' B) ;
*
     RV . 'CO1' = CU ;
     RV . 'CO2' = CU ;
 'FINSI' ;
 'FINSI' ; 
*  
 IALPDT = 'PROG' ;
 IRESU IJACO ='KOPS' 'MATRIK' ;
*
 'MENAGE' ;
*
 'SI' (RVX . 'IMPL') ;
     'RESPRO' IJACO IRESU IALPDT ;
 'SINON' ;
     'RESPRO' IRESU IALPDT ;
 'FINSI' ;
*
 'FINP' ;
*
 
************************************************************************
************************************************************************
***************** FIN PARTIE PROCEDURES ********************************
************************************************************************
************************************************************************
************************************************************************
 
************************************************************************
************************************************************************
**************************** MESH **************************************
************************************************************************
************************************************************************
 
 NY = 5 '*' RAF ;
 NX1 = 4 '*' RAF ;
 NX2 = 2 '*' NX1 ;
 NX3 = NX1  ;
 NX = (NX1 '+' NX2 '+' NX3)  ;
 DX = (4.0 '/' NX)  ;
 
 A0 = -2.0 0.0 ;
 A1 = -1.0 0.0 ;
 A2 = 1.0 0.0 ;
 A3 = 2.0 0.0 ;
 A4 = 2.0 1.0 ;
 A5 = -2.0 1.0 ;
 
*
**** LIGB
*
 
 LIGB1 = A0 'DROIT' NX1 A1 ;
 
* LIGB2 (On utilise un propriete de 'ET' ; si 'ET' change ?)
 
 xcel = ('COORDONNEE' 1 A1) '+' DX ;
 ycel = 0.1 '*' ( 1.0 '+' ('COS' (180 '*' xcel)));
 ACEL = xcel ycel ;
 LIGB2 = A1 'DROIT' 1 ACEL ;
 'REPETER' BL1 (NX2 '-' 2) ;
    ACEL0 = ACEL ;
    xcel = xcel '+' DX ;
    ycel = 0.1 '*' ( 1.0 '+' ('COS' (180 '*' xcel)));
    ACEL = xcel ycel ;
    LIGB2 = LIGB2 'ET' (ACEL0 'DROIT' 1 ACEL) ;
 'FIN' BL1;
 LIGB2 = LIGB2 'ET' (ACEL 'DROIT' 1 A2) ;
 
 
 LIGB3 = A2 'DROIT' NX3 A3 ;
 
 LIGB = LIGB1 'ET' LIGB2 'ET' LIGB3 ;
 
*
**** LIGH
*
 
 LIGH = A4 'DROIT' NX A5 ;
 
*
**** DOMINT
*
 
 DOMINT =  LIGH 'REGLER' NY ('INVERSE' LIGB) ;
 LIGCON = 'CONTOUR' DOMINT ;
 
*
*** LIGG
*
 
 LIGG = LIGCON 'ELEM' 'COMP' A5  A0 ;
 
*
**** LIGD
*
 
 LIGD = LIGCON 'ELEM' 'COMP' A3  A4 ;
 
*
**** MODEL OBJECTS
*
 
 MDOMINT  = 'MODELISER' DOMINT 'EULER' ;
 MLIGG    = 'MODELISER' LIGG   'EULER' ;
 MLIGD    = 'MODELISER' LIGD   'EULER' ;
 
*
**** Creation of DOMAINE tables via the MODEL object
*
 
 TDOMINT = 'DOMA' MDOMINT 'VF' ;
 TLIGG   = 'DOMA' MLIGG   'VF' ;
 TLIGD   = 'DOMA' MLIGD   'VF' ;
 
 QDOMINT = TDOMINT . 'QUAF' ;
 QLIGD = TLIGD . 'QUAF' ;
 QLIGG = TLIGG . 'QUAF' ;
 
 'ELIMINATION' QDOMINT (DX '/' 100.)  QLIGG ;
 'ELIMINATION' QDOMINT (DX '/' 100.)  QLIGD ;
 
 'SI' GRAPH ;
    'TRACER' (DOMINT 'ET' (LIGG 'COULEUR' 'ROUG')
       'ET' (LIGD 'COULEUR' 'BLEU'))
       'TITRE' 'Domaine total' ;
 'FINSI' ;
 
 
*******************************************
****** LIGNE de Post-traitement ***********
*******************************************
 
 LIGB = LIGCON 'ELEM' 'COMPRIS' A0 A3 ;
 POIN0 =  LIGB 'POIN' 1;
 X1 Y1 = 'COORDONNEE' POIN0 ;
 'REPETER' BLLIM (('NBNO' LIGB) '-' 1) ;
    X0 = X1 ;
    Y0 = Y1 ;
    POIN0 = LIGB 'POIN' (&BLLIM '+' 1) ;
    X1 Y1 = 'COORDONNEE' POIN0 ;
    XFAC = (X0 '+' X1) '/' 2 ;
    YFAC = (Y0 '+' Y1) '/' 2 ;
    PFAC = (TDOMINT . 'FACE') 'POIN' 'PROC' (XFAC YFAC);
    GEOFAC1 = (TDOMINT . 'FACEL') 'ELEM' 'APPUYE'
             'LARGEMENT' PFAC ;
*
**** Tranformation en POI1
*
    GEO1POI1 = 'CHANGER' 'POI1' GEOFAC1 ;
    PCEL11 = GEO1POI1 'POIN'  1 ;
    PCEL12 = GEO1POI1 'POIN'  2 ;
*
**** Il faur verifier que PFAC = PCEL12 = PCEL22
*    ('NBEL' GEO1POI1) = ('NBEL' GEO2POI1) = 2
*
   'SI'  ( ('NBEL' GEO1POI1) 'NEG' 2);
       'MESSAGE' ;
       'MESSAGE'
          'Probleme dans la creation de la ligne pour le post.';
       'MESSAGE' ;
       'ERREUR' 21 ;
   'FINSI' ;
   'SI'  ( PCEL12 'NEG' PFAC);
       'MESSAGE' ;
       'MESSAGE'
          'Probleme dans la creation de la ligne pour le post.';
       'MESSAGE' ;
       'ERREUR' 21 ;
   'FINSI' ;
*
*** Creation d'un maillage SEG2
*
   'SI' (&BLLIM 'EGA' 1);
       PCEN0 = PCEL11 ;
    'SINON' ;
      'SI' (&BLLIM 'EGA' 2);
          LIGPOST = 'MANUEL' 'SEG2' PCEN0 PCEL11 ;
          PCEN0 = PCEL11 ;
      'SINO' ;
          LIGPOST = LIGPOST 'ET' ('MANUEL' 'SEG2' PCEN0 PCEL11) ;
          PCEN0 = PCEL11 ;
      'FINSI' ;
    'FINSI' ;
 'FIN' BLLIM ;
 
 'SI' GRAPH ;
    'TRACER' (DOMINT 'ET' (LIGPOST 'COULEUR' 'VERT'))
       'TITRE' 'LIGPOST' ;
 'FINSI' ;
 
***************************************************************
***************************************************************
***************************************************************
************** Initial conditions *****************************
***************************************************************
***************************************************************
***************************************************************
*
*** C.L. et initiales
*
 
 C_INF = GAMSCAL * P_INF '/' RO_INF ;
 C_INF = C_INF '**' 0.5 ;
 M_INF = U_INF '/' C_INF ;
 
* 
* Names of conserved variables and others
*
 NOMDEN = 'MOTS' 'RN' ;
 NOMMOM = 'MOTS' 'RUX' 'RUY' ;
 NOMRET = 'MOTS' 'RETN' ;
 NOMVEL = 'MOTS' 'UX' 'UY' ;
* NOMPRE = 'MOTS' 'NOMPRE' ;
 LISTCONS = NOMDEN 'ET' NOMMOM 'ET' NOMRET ;
*
 
 GAMN = 'MANUEL' 'CHPO' (TDOMINT . 'CENTRE') 1 'SCAL' gamscal;
 
 
 rek_inf = 0.5D0 * u_inf * u_inf * ro_inf ;
 rei_inf = P_INF '/' (gamscal '-' 1.0D0) ;
 
 RN0 = 'MANUEL' 'CHPO' (TDOMINT . 'CENTRE') 1 'SCAL' ro_inf ;
 GN0 = 'MANUEL' 'CHPO' (TDOMINT . 'CENTRE') 2 'UX' (ro_inf * u_inf)
   'UY' 0.0 ;
 RETN0 = 'MANUEL' 'CHPO' (TDOMINT . 'CENTRE') 1 'SCAL'
    (rei_inf '+' rek_inf)  ; 
 
 UNCONS = ('NOMC' 'RN' RN0 'NATU' 'DISCRET') 'ET' ('NOMC'
    ('MOTS' 'UX' 'UY') GN0 ('MOTS' 'RUX' 'RUY') 'NATU' 'DISCRET') 'ET'
    ('NOMC' 'RETN' RETN0 'NATU' 'DISCRET') ; 
 
 VN0 PN0 = 'PRIM' 'PERFMONO' RN0 GN0 RETN0 GAMN ;
 ERRO = ('MAXIMUM' (PN0 '-' p_inf) 'ABS') / p_inf ;
 'SI' (ERRO > 1.0D-16) ;
    'ERREUR' 21 ;
 'FINSI' ;
 
 aa = 'EXCO' 'UX' VN0  ;
 ERRO = ('MAXIMUM' (aa '-' u_inf) 'ABS') / u_inf ;
 'SI' (ERRO > 1.0D-12) ;
    'ERREUR' 21 ;
 'FINSI' ;
 
 VN2 = 'PSCAL' VN0 VN0 ('MOTS' 'UX' 'UY') ('MOTS' 'UX' 'UY') ;
 C2 = GAMN '*' ( PN0 '/' RN0) ;
 CN0 = C2 '**' 0.5 ;
 
 MACH2 = VN2 '/' C2;
 MACHN0 = MACH2 '**' 0.5;
 
*
**** Plot of IC
*
 
 'SI' GRAPH ;
 
    CHM_RN   =  'KCHA' MDOMINT 'CHAM' RN0 ;
    CHM_PN   =  'KCHA' MDOMINT 'CHAM' PN0 ;
    CHM_VN  =  'KCHA' MDOMINT 'CHAM'  VN0 ;
    CHM_MACH  =  'KCHA' MDOMINT 'CHAM' MACHN0 ;
 
    'TRACER'  CHM_RN  MDOMINT
       'TITR'  ('CHAINE'  'RN at t='  0.0);
    'TRACER'  CHM_PN  MDOMINT
       'TITR'  ('CHAINE'  'PN at t='  0.0);
    'TRACER'  CHM_VN MDOMINT
       'TITR'  ('CHAINE' 'VN at t='  0.0);
    'TRACER'  CHM_MACH MDOMINT
       'TITR'  ('CHAINE' 'MACH at t='  0.0);
 
 'FINSI' ;
 
 
*****************************************************
*****************************************************
*****************************************************
**************** La table RV   **********************
*****************************************************
*****************************************************
*****************************************************
*
 RV = 'TABLE' ;
 RV . 'ANOM' = FAUX ;
* Constant time step
 CFLREF = 500.0 ;
 UREF = u_inf ;
 UREF = c_inf ;
 RV . 'DTIMP' = CFLREF * DX '/' UREF ;
*
 RV . 'ORDTPS' = 1 ;
*
 RV . 'UN' = UNCONS ;
 RV . 'LISTCONS' = LISTCONS ;
 RV . 'ERROMOTS' = NOMRET ;
*
 RV . 'CLIM' = FAUX ;
*
 RV . 'CONS' = FAUX ;
*
 RV . 'TFINAL' = 1.0E6 ;
 RV . 'NITMAEX' = NITER ;
*
 RV . 'NITMAIN' = 6 ;
 RV . 'NITMIIN' = 6 ;
 RV . 'EPSINT' = 1.0D-16 ;
*
* RV . 'FEXT' = 20 ;
* RV . 'FINT' = 1 ;
*
 RV . 'LISTOPER' = 'MOTS'  'CALC' 'PKON'  ;
* 
 CALCTAB = 'TABLE' ;
 RV . '1CALC' = CALCTAB ;
 PKONTAB = 'TABLE'   ;
 RV . '2PKON' = PKONTAB ;
*
 RV . 'MATIDE' = 'KOPS' 'MATIDE' LISTCONS ('DOMA' MDOMINT 'CENTRE')
   'MATRIK'  ;
*
**** Parametres de la procedure PROLIM
*
* Gas model
* Gas properties
 RV . 'PGAZ' = 'TABLE' ;
 RV . 'PGAZ' . 'R' = Rair ;
 RV . 'PGAZ' . 'GAMN' = gamn ;
* RV . 'PGAZ' . 'MU' = 0.0 ;
* RV . 'PGAZ' . 'LAMBDA' = 0.0 ;
* Bas Mach (PKON)
 RV . 'BASMACH' = VRAI ;
 RV . 'CO1' = 'MANUEL' 'CHPO' ('DOMA' MDOMINT 'CENTRE')
              1 'SCAL' (10 * U_INF) ;
 RV . 'CO1MAX' = 'MANUEL' 'CHPO' ('DOMA' MDOMINT 'CENTRE')
              1 'SCAL' (10 * U_INF) ;
 RV . 'CO1MIN' = 'MANUEL' 'CHPO' ('DOMA' MDOMINT 'CENTRE')
              1 'SCAL' (10 * U_INF) ;
 RV . 'CO2' = 'COPIER' ( RV . 'CO1') ;
* RV . 'CO2' can be Modified into PDIF 
*****************
** CALCUL *******
*****************
* Personal procedure
* 
 CALCTAB . 'ANOM' = FAUX ;
 CALCTAB . 'EQEX' = RV ;
 CALCTAB . 'MODELE' = MDOMINT ;
* We call this procedure during the external iterations -> 
 CALCTAB . 'IMPL' = VRAI ;
 CALCTAB . 'COMPT' = 0 ;
*
*****************
* PKON **********
*****************
 CHPVID MATVID  = 'KOPS' 'MATRIK' ;
 PKONTAB . 'CFLDTAU' = 10000. ;
 PKONTAB . 'ANOM' = FAUX ;
 PKONTAB . 'EQEX' = RV ;
 PKONTAB . 'GAZ' = 'PERFMONO' ;
 PKONTAB . 'MODELE' = MDOMINT ;
 PKONTAB . 'LISTCONS' = LISTCONS ;
* PKONTAB . 'METHODE' = 'CENTERED' ;
* PKONTAB . 'METHODE' = 'RUSANOLM' ;
 PKONTAB . 'METHODE' = 'AUSMPLM' ;
* PKONTAB . 'METHODE' = 'AUSMPLUS' ;
* PKONTAB . 'METHODE' = 'VLH' ;
* PKONTAB . 'METHODE' = 'ROELM' ;
* PKONTAB . 'METHODE' = 'ROE' ;
* PKONTAB . 'METHODE' = 'HLLCLM' ;
*
 PKONTAB . 'IMPL' =  VRAI ;
 'SI' (PKONTAB . 'IMPL') ;
*    PKONTAB . 'TYPEJACO' = 'AUSMPLUS' ;
*    PKONTAB . 'TYPEJACO' = 'VLH' ;
*    PKONTAB . 'TYPEJACO' = 'RUSANOLM' ;
    PKONTAB . 'TYPEJACO' = 'AUSMPLM' ;
 'FINSI' ;
*
 PKONTAB . 'ORDREESP' = 2 ;
 'SI' ((PKONTAB . 'ORDREESP') 'EGA' 2) ;
    PKONTAB . 'VLIM' = CHPVID ;
    CHPENLIM = 'MANUEL' 'CHPO' (('DOMA' MLIGG 'CENTRE') 'ET'
       ('DOMA' MLIGD 'CENTRE')) 1 'SCAL' 0.0 ;
    CHPE2LIM = 'MANUEL' 'CHPO' (('DOMA' MLIGG 'CENTRE') 'ET'
      ('DOMA' MLIGD 'CENTRE')) 2
       'UX' 1.0 'UY' 0.0 ;
    TOTO TITI MCHRCON = 'PENTE' MDOMINT 'CENTRE' 'EULESCAL'
             'LIMITEUR' ('MOTS' 'SCAL') RN0 'CLIM' CHPENLIM ;
    TOTO TITI MCHVCON = 'PENTE' MDOMINT 'CENTRE' 'EULEVECT'
             'LIMITEUR' ('MOTS' 'UX' 'UY') GN0 'CLIM' CHPE2LIM ;
    PKONTAB . 'LIMITEUR' = 'NOLIMITE' ;
    PKONTAB . 'GRADRN' = MCHRCON ;
    PKONTAB . 'GRADPN' = MCHRCON ;
    PKONTAB . 'GRADVN' = MCHVCON ;
 'FINSI' ;
 PKONTAB . 'NFROZLIM' = 50 ;
* Gravité
* PKONTAB . 'CORRGRA' = FAUX ;
 PKONTAB . 'GRAVITE' =  'MANUEL' 'CHPO' ('DOMA' MDOMINT 'CENTRE') 2
   'UX' 0.0 'UY' 0.0 ;
*
 MAIVID = 'DIFF' ('DOMA' MLIGG 'CENTRE')
   ('DOMA' MLIGG 'CENTRE') ;
 PKONTAB . 'FACELIM' = MAIVID ;
 PKONTAB . 'ALPHA' = 'INF' ;
 PKONTAB . 'CLIM' = VRAI ;
 PKONTAB . 'TABLIM' = TABLE ;
*
* Normalisation du profile parabolique, pour obtenir que
* le flux rentrante vaut MINJ
*
 PKONTAB . 'TABLIM' . 'CHINRI' = 'MANUEL' 'CHPO' ('DOMA' MLIGG
 'CENTRE') 4 'RN' ro_inf 'UX' U_INF
 'UY' 0.0 'PN' p_inf ;
 PKONTAB . 'TABLIM' . 'CHOUTR' = 'MANUEL' 'CHPO' ('DOMA' MLIGD
 'CENTRE') 4 'RN' ro_inf 'UX' U_INF
 'UY' 0.0 'PN' p_inf ;
*
*********************************** 
* Inversion de la matrice *********
***********************************
*
* RV . 'MATUPDAT' updating
* RV . 'MUPEXT'   external iterations updating
* RV . 'MUPINT'   internal iterations updating
* RV . 'MUPLIN'   We update the matrix if in the previous
*                 internal iteration the number of linear
*                 iterations to solve the system were bigger
*                 than (RV . 'MUPLIN')
 RV . 'MATUPDAT' = FAUX ;
 'SI' ('NON' ( RV . 'MATUPDAT')) ;
    RV . 'MUPEXT' = 1 ;
    RV . 'MUPINT' = 50 ;
    RV . 'MUPLIN' = 4000 ;
 'FINSI' ;
 RV . 'MATINV' = 'TABLE' 'METHINV' ;
 MATTAB = RV . 'MATINV' ;
* MATTAB .  'TYPINV' = 1 : methode exact
* MATTAB .  'TYPINV' = 5 ; GMRES
 MATTAB .  'TYPINV' = 1 ;
 MATTAB .  'IMPINV' = 0 ;
* 
* Matrice pour assurer que la matrice à inverser est correctement
* assemblé
*
*  MATTAB  . 'MATASS'  definie en EXEXIM
*  MATTAB  . 'MAPREC'         " 
*  MATTAB  . 'XINIT'          "   
*
 
* Methode de numerotation de DDL
  MATTAB . 'TYRENU' = 'RIEN' ;
  MATTAB . 'PCMLAG' = 'APR2' ;
  MATTAB . 'GMRESTRT' = 500 ;
  MATTAB . 'SCALING' = 1 ;
* ILU 3
* ILUT (dual truncation) 5
* ILUTP 9
* Dans le cas ILUT il faut definir
* ILUTLFIL
* ILUTDTOL
  MATTAB . 'PRECOND' = 7 ;
  MATTAB . 'OUBMAT' = 0 ;
  MATTAB . 'ILUTPPIV' = 0.1 ;
  MATTAB . 'ILUTPMBL' = 0 ;
  MATTAB . 'ILUTLFIL' = 4. ;
  MATTAB . 'ILUTDTOL' = 0. ;
  MATTAB . 'NITMAX' =  3000 ;
  MATTAB . 'RESID' = 1.D-10 ;
*
**** Save results into a file
*
 SI FAUX ;
    FICHER1 = ('CHAINE' FICHER 'main' RAF
           (PKONTAB . 'METHODE')
           'OE' (PKONTAB . 'ORDREESP')
           'OT' (RV . 'ORDTPS') '.sauv') ;
    'OPTION' 'SAUVER' FICHER1 ;
    'MESSAGE' 'We save results into' ;
    'MESSAGE' FICHER ;
    'MESSAGE' ;
 'FINSI' ;
*
**** Exexcution EXEXIM
*
 'SI' FAUX ;
   LISTTPS = ('PROG' 100.0) ;
   'REPETER' BL1 ('DIME' LISTTPS) ;
      RV . 'TFINAL' = 'EXTRAIRE' LISTTPS &BL1 ;
      'TEMPS' 'ZERO' ;
         EXEXIM RV ;
      'TEMPS' 'IMPR' ;
*
**** To save memory storage
*
      'OUBLIER' MATTAB MATASS ;
      'OUBLIER' MATTAB MAPREC ; 
      'MENAGE' ;
      'SI' FAUX ;
         'SAUVER' 'LABEL' ('CHAINE' &BL1) ;
      'FINSI' ;
    'FIN' BL1 ;
 SINON ;
*   'TEMPS' 'ZERO' ;
    EXEXIM RV ;
*   'TEMPS' 'IMPR' ;
 'FINSI' ;
*
***************************************************************
***************************************************************
***************************************************************
************** Plot of the solutions **************************
***************************************************************
***************************************************************
***************************************************************
*
* GRAPH = VRAI ;
*
 'SI' FAUX ;
     FICHER =
     './D_sauv/main16VLHOE2OT1.sauv' ;
     'OPTION' 'RESTITUER' FICHER ;
     'RESTITUER' ;
 'FINSI' ;
*
* LOGGNP = logical. If VRAI, files for gnuplot are created
*
 LOGGNP = VRAI ;
 VALECH = 1 ; 
 
 'OPTION' 'ECHO' VALECH 'TRAC' 'X' ;
 
 TPS = RV . 'PASDETPS' . 'TPS' ;
 METO = PKONTAB . 'METHODE' ;
 ORD_ESP = PKONTAB . 'ORDREESP' ;
 ORD_TPS = RV . 'ORDTPS' ;
 
 NTOTEL = 'NBNO' ('DOMA' MDOMINT 'CENTRE') ;
 
*
**** Les variables conservatives
*
 
 GAMN = 'REDU' GAMN ('DOMA' MDOMINT 'CENTRE') ;
 RN = 'EXCO' NOMDEN (RV . 'UN') ('MOTS' 'SCAL') ;
 GN = 'EXCO' NOMMOM (RV . 'UN') ('MOTS' 'UX' 'UY') ;
 RETN = 'EXCO' NOMRET (RV . 'UN') ('MOTS' 'SCAL') ;
 
*
**** Les variables primitives
*
 
 VN PN =   'PRIM' 'PERFMONO' 
           RN GN RETN GAMN ;
 
 CSON2 = (GAMN '*' PN) '/' RN ;
 VN2 = 'PSCAL' VN VN ('MOTS' 'UX' 'UY') ('MOTS' 'UX' 'UY') ;
 MACHN2 = VN2 '/' CSON2 ;
 MACHN = MACHN2 '**' 0.5 ;
 
 TITOLO =  ('CHAINE' METO ', OE = ' ORD_ESP ', OT =' ORD_TPS 
    ', NBEL = ' NTOTEL ' ,tps =' TPS) ;
 
 EVO1 = 'EVOL' 'MANU' 'it' (CALCTAB . 'IT') 'err'
    ((CALCTAB . 'ER') '+' 1.0D-16) ;
 
*
*** GRAPHIQUE DES SOLUTIONS
*
 
 'SI' GRAPH ;
 
    'TRACER' DOMINT 'TITRE' ('CHAINE' 'Maillage') ;
 
    CHM_RN   =  'KCHA' MDOMINT 'CHAM' RN ;
    CHM_PN   =  'KCHA' MDOMINT 'CHAM' PN ;
    CHM_VN  =  'KCHA' MDOMINT 'CHAM'  VN ;
    CHM_MACH =   'KCHA' MDOMINT 'CHAM' MACHN ;
 
    'TRACER'  CHM_RN  MDOMINT ('CONTOUR' DOMINT)
       'TITR'  ('CHAINE' ' RN : ' TITOLO) ;
    RNV   =  'ELNO' TDOMINT ('REDU' RN TDOMINT . 'CENTRE') ;
    'OPTION' 'ISOV' 'LIGN' ;
    'TRACER' DOMINT RNV ('CONTOUR' DOMINT)
      'TITRE' ('CHAINE' 'RN : ' TITOLO) 15 ;
    'OPTION' 'ISOV' 'SULI' ;
    'TRACER'  CHM_PN  MDOMINT  ('CONTOUR' DOMINT)
       'TITR'  ('CHAINE' 'PN : ' TITOLO) ;
    PNV   =  'ELNO' TDOMINT ('REDU' PN TDOMINT . 'CENTRE') ;
    'OPTION' 'ISOV' 'LIGN' ;
    'TRACER' DOMINT PNV ('CONTOUR' DOMINT)
      'TITRE' ('CHAINE' 'PN : ' TITOLO) 15 ;
    'OPTION' 'ISOV' 'SULI' ;
    'TRACER'  CHM_MACH MDOMINT  ('CONTOUR' DOMINT)
       'TITR'  ('CHAINE' 'MACHN : ' TITOLO);
    MACHNV   =  'ELNO' TDOMINT ('REDU' MACHN TDOMINT . 'CENTRE') ;
    'OPTION' 'ISOV' 'LIGN' ;
    'TRACER' DOMINT MACHNV ('CONTOUR' DOMINT)
      'TITRE' ('CHAINE' 'Mach : ' TITOLO) 15 ;
    'OPTION' 'ISOV' 'SULI' ;
    'TRACER'  CHM_VN  MDOMINT  ('CONTOUR' DOMINT)
       'TITR'  ('CHAINE'  'VN : ' TITOLO) ;
    VECN = 'VECTEUR' VN (8. '/' RAF) 'UX' 'UY' 'JAUNE' ;
    'TRACER' DOMINT VECN ('CONTOUR' DOMINT)
       'TITRE' ('CHAINE' 'VN : ' TITOLO) ;
*
    'OPTION' 'ISOV' 'SULI' ;
    'DESSIN' EVO1 LOGY 'TITRE' 'Error' ;
*
 'FINSI' ;
 
*
**** Test de convergence
*
 
 AA = 'EXTRAIRE' (CALCTAB . 'ER')  ('DIME' (CALCTAB . 'ER')) ;
 
 'SI'  (AA > 1.0D-8) ;
    'MESSAGE' 'Probleme' ;
    'ERREUR' 5 ;
 'FINSI' ;
*
 
 XLIGPOST = 'COORDONNEE' 1 LIGPOST ;
 EVCOOR = 'EVOL' 'CHPO' XLIGPOST 'SCAL' LIGPOST ;
 LISTX = 'EXTRAIRE' EVCOOR 'ORDO' ;
 
 EVMACH =  'EVOL' 'CHPO' MACHN 'SCAL' LIGPOST ;
 LISTMACH = 'EXTRAIRE' EVMACH 'ORDO' ;
 EVMACH = 'EVOL' 'MANU' 'x' LISTX 'Mach' LISTMACH ;
 'SI' GRAPH ;
   'DESSIN' EVMACH 'MIMA' 'TITRE' TITOLO ;
 'FINSI' ;
 
 MMAX = 'MAXIMUM' MACHN ;
 MMIN = 'MINIMUM' MACHN ;
 
 MMAX_REF = 1.47 * M_INF ;
 MMIN_REF = 0.863 * M_INF ;
 
 ERRO = 'ABS' (MMAX_REF '-' MMAX)  '/' M_INF ;
 'SI'  (ERRO > 1.0D-2) ;
    'MESSAGE' 'Probleme' ;
    'ERREUR' 5 ;
 'FINSI' ;
*
 
 ERRO = 'ABS' (MMIN_REF '-' MMIN)  '/' M_INF ;
 'SI'  (ERRO > 1.0D-2) ;
    'MESSAGE' 'Probleme' ;
    'ERREUR' 5 ;
 'FINSI' ;
 
 'FIN' ;