Télécharger 15wedge.dgibi

Retour à la liste

* fichier :  15wedge.dgibi
*
******************************************************************
* CALCUL DE L'ECOULEMENT SUPERSONIQUE STATIONNAIRE DANS UN CANAL *
* AVEC RAMPE INCLINEE A 15deg.                                   *
* POWELL AND DEZEEUW, AIAA-91-1542-CP                            *
* FORMULATION VF COMPRESSIBLE EXPLICITE/IMPLICITE                *
*                                                                *
* H. PAILLERE/P. GALON TTMF AOUT 1997                            *
*                                                                *
* MODIF BECCANTINI MARS 97 (reconstruction lineaire exacte)      *
* MODIF BECCANTINI MARS 01 (implicite)                           *
* MODIF BECCANTINI AUGUST 2021 (quelque sintaxe a changé         *
*                                                                *
******************************************************************
*
************************************************************************
************************************************************************
************************************************************************
************************************************************************
***************** 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' ;
             SMB = 'CHAN' 'ATTRIBUT' SMB 'NATURE' 'DISCRET' ;
*
********* 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 . 'CHINSS') 'INSS' ;
 
 RJACO1 = 'KONV' 'VF' 'PERFMONO' 'CLIM'  'JACOCONS'
        MDOMINT MLIGG LISTCONS LISTP
        (TABLIM . 'RN') (TABLIM . 'VN') (TABLIM . 'PN')
        (TABLIM . 'GAMN') (TABLIM . 'CHINSS') 'INSS' ;
 
 RCHLIM2 RCHRES2 = 'KONV' 'VF' 'PERFMONO' 'CLIM'  'RESI'
        MDOMINT MLIGD LISTCONS LISTP
        (TABLIM . 'RN') (TABLIM . 'VN') (TABLIM . 'PN')
        (TABLIM . 'GAMN') (TABLIM . 'CHOUTSS') 'OUTSS' ;
 
 RJACO2 = 'KONV' 'VF' 'PERFMONO' 'CLIM'  'JACOCONS'
        MDOMINT MLIGD LISTCONS LISTP
        (TABLIM . 'RN') (TABLIM . 'VN') (TABLIM . 'PN')
        (TABLIM . 'GAMN') (TABLIM . 'CHOUTSS') 'OUTSS' ;
 
 
 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 ********************************
************************************************************************
************************************************************************
*
 
 TYEL = 'QUA4' ;
 
 'OPTION'  'DIME' 2  'ELEM' TYEL  'ISOV' 'SULI'
    'ECHO' 1 'TRAC' 'X' ;
 
 
 GRAPH = VRAI ;
 GRAPH = FAUX ;
 
******************
******************
******************
**** MAILLAGE ****
******************
******************
******************
******************
*
*
 RAF = 16 ;
 
 NX1 = 1 '*' RAF ;
 NX2 = 1 '*' RAF ;
 NX3 = 4 '*' RAF ;
 NY = 2 '*' RAF ;
 DY = 2.0 '/' NY ;
 
 H1 = 0.268 ;
 
 A1 = (1. '-' DY) 0.0 ;
 A2 = 1.0 0.0 ;
 A3 = 2.0 0.0 ;
 A4 = 3.0 H1 ;
 A5 = 7.0 H1 ;
 A6 = (7.0 '+' DY) H1 ;
 
 B1 = (1.0 '-' DY) 2.0 ;
 B2 = 1.0 2.0 ;
 B3 = 2.0 2.0 ;
 B4 = 3.0 2.0 ;
 B5 = 7.0 2.0 ;
 B6 = (7.0 '+' DY) 2.0 ;
 
 BAS1 = A2 'DROI' NX1 A3 ;
 BAS2 = A3 'DROI' NX2 A4 ;
 BAS3 = A4 'DROI' NX3 A5 ;
 HAU3 = B5 'DROI' NX3 B4 ;
 HAU2 = B4 'DROI' NX2 B3 ;
 HAU1 = B3 'DROI' NX1 B2 ;
 
 COTEG = B2 'DROI' NY A2 ;
 COTED = A5 'DROI' NY B5 ;
 
 DOMINT = 'DALLER' (BAS1 ET BAS2 ET BAS3) COTED
     (HAU3 ET HAU2 ET HAU1) COTEG 'PLAN' ;
 
 LIGG = COTEG ;
 LIGD = COTED ;
 
 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 (DY '/' 100.)  QLIGG ;
 'ELIMINATION' QDOMINT (DY '/' 100.)  QLIGD ;
 
 'SI' GRAPH ;
    'TRACER' (DOMINT 'ET' (LIGG 'COULEUR' 'ROUG')
       'ET' (LIGD 'COULEUR' 'BLEU'))
       'TITRE' 'Domaine total' ;
 'FINSI' ;
 
 'SI' FAUX ;
    'OPTION' 'SAUVER' ('CHAINE' './D_sauv/'
    'mail' RAF '.sauv') ;
    'SAUV' ;
    'FIN' ;
 'FINSI' ;
*
***************************************************************
***************************************************************
***************************************************************
************** Initial conditions *****************************
***************************************************************
***************************************************************
***************************************************************
*
 'SI' FAUX ;
    RAF = 16 ;
    FICHER = './D_sauv/' ;
    'OPTION' 'RESTITUER' ('CHAINE' FICHER 'mail' RAF '.sauv') ;
    'RESTITUER' ;
 'FINSI' ;
*
* 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 ;
*
 
 gamscal = 1.4 ;
 GAMN = 'MANUEL' 'CHPO' (TDOMINT . 'CENTRE') 1 'SCAL' gamscal;
 
 Minf = 2.0D0 ;
 alpha = 0.D0 ;
 RAIR = 288. ;
 
 'SI' ('NEG' ALPHA 0.0) ;
   'MESSAGE' 'Attention aux C.L.' ;
   'ERREUR' 21 ;
 'FINSI' ;
 
 rhoinf = 1.0D0 ;
 rhouinf = 'COS' alpha ;
 rhovinf = 'SIN' alpha ;
 rhoEtinf = 2.0D0 '+' (gamscal '*' Minf '*' Minf
      '*' (gamscal '-' 1.0D0)) ;
 rhoEtinf = rhoEtinf '/' (2.D0 '*' gamscal '*' Minf
      '*' Minf '*' (gamscal '-' 1.0D0)) ;
 
 RN0 = 'MANUEL' 'CHPO' (TDOMINT . 'CENTRE') 1 'SCAL' rhoinf ;
 GN0 = 'MANUEL' 'CHPO' (TDOMINT . 'CENTRE') 2 'UX' rhouinf
   'UY' rhovinf ;
 RETN0 = 'MANUEL' 'CHPO' (TDOMINT . 'CENTRE') 1 'SCAL' rhoEtinf ;
 
 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 ;
 pinf = 'MAXIMUM' PN0 ;
 
 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;
 TN0 = PN0 '/' (Rair * RN0) ;
 
*
**** Plot of IC
*
 
 'SI' GRAPH ;
 
    CHM_RN   =  'KCHA' MDOMINT 'CHAM' RN0 ;
    CHM_PN   =  'KCHA' MDOMINT 'CHAM' PN0 ;
    CHM_VN  =  'KCHA' MDOMINT 'CHAM'  VN0 ;
 
    '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' 'MACHN at t='  0.0);
 
 'FINSI' ;
 
 'SI' FAUX ;
*
    'OPTION' 'SAUVER'
       ('CHAINE' FICHER  'ic'  RAF '.sauv') ;
     'SAUV' ;
 'FINSI' ;
*
*
 'SI' FAUX ;
     FICHER = './D_sauv/' ;
    'OPTION' 'RESTITUER' ('CHAINE' FICHER 'ic16.sauv') ;
    'RESTITUER' ;
 'FINSI' ;
*
*****************************************************
*****************************************************
*****************************************************
**************** La table RV   **********************
*****************************************************
*****************************************************
*****************************************************
*
 RV = 'TABLE' ;
 RV . 'ANOM' = FAUX ;
* Constant time step
 CFLREF = 50. ;
 UREF = rhouinf '/' rhoinf ;
 RV . 'DTIMP' = CFLREF * DY '/' UREF ;
*
 RV . 'ORDTPS' = 1 ;
*
 RV . 'UN' = UNCONS ;
 RV . 'LISTCONS' = LISTCONS ;
 RV . 'ERROMOTS' = NOMRET ;
*
 RV . 'CLIM' = FAUX ;
*
 RV . 'CONS' = FAUX ;
*
 RV . 'TFINAL' = 0. ;
 RV . 'NITMAEX' = 100000 ;
*
 RV . 'NITMAIN' = 10 ;
 RV . 'NITMIIN' = 10 ;
 RV . 'EPSINT' = 1.0D-3 ;
*
* 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' = FAUX ;
* RV . 'CO1' = 'MANUEL' 'CHPO' ('DOMA' MDOMINT 'CENTRE')
*              1 'SCAL' (1.0D2) ;
* RV . 'CO1MAX' = 'MANUEL' 'CHPO' ('DOMA' MDOMINT 'CENTRE')
*              1 'SCAL' (1.0D2) ;
* RV . 'CO1MIN' = 'MANUEL' 'CHPO' ('DOMA' MDOMINT 'CENTRE')
*              1 'SCAL' (1.0D0) ;
* 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' = 'VLH' ;
* PKONTAB . 'METHODE' = 'ROELM' ;
* PKONTAB . 'METHODE' = 'HLLC' ;
*
 PKONTAB . 'IMPL' =  VRAI ;
 'SI' (PKONTAB . 'IMPL') ;
*    PKONTAB . 'TYPEJACO' = 'AUSMPLUS' ;
    PKONTAB . 'TYPEJACO' = 'VLH' ;
 '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' = 'LIMITEUR' ;
    PKONTAB . 'GRADRN' = MCHRCON ;
    PKONTAB . 'GRADPN' = MCHRCON ;
    PKONTAB . 'GRADVN' = MCHVCON ;
 'FINSI' ;
* 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' . 'CHINSS' = 'MANUEL' 'CHPO' ('DOMA' MLIGG
 'CENTRE') 4 'RN' rhoinf 'UX' (rhouinf '/' rhoinf)
 'UY' (rhovinf '/' rhoinf) 'PN' pinf ;
 PKONTAB . 'TABLIM' . 'CHOUTSS' = CHPVID ;
*
***********************************
* 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' = 5 ;
 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' = 5 ;
  MATTAB . 'OUBMAT' = 0 ;
  MATTAB . 'ILUTPPIV' = 0.1 ;
  MATTAB . 'ILUTPMBL' = 0 ;
  MATTAB . 'ILUTLFIL' = 4. ;
  MATTAB . 'ILUTDTOL' = 0. ;
  MATTAB . 'NITMAX' =  3000 ;
  MATTAB . 'RESID' = 1.D-5 ;
*
**** 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
*
  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 ;
 
***************************************************************
***************************************************************
***************************************************************
************** Plot of the solutions **************************
***************************************************************
***************************************************************
***************************************************************
*
 '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 en KONV' ;
    'ERREUR' 5 ;
 'FINSI' ;
 
*
*** Test
*
 
 PCEN1 = (TDOMINT . 'CENTRE') 'POIN' 'PROC' (2.5 (H1 '/' 2.0)) ;
 'SI' GRAPH ;
    'TRACER' (DOMINT 'ET' PCEN1) 'TITRE'
     'Test 1: solution derrier le premier choc' ;
 'FINSI' ;
 ro1 = 'EXTRAIRE' RN 'SCAL' PCEN1 ;
 
 erro = ((ro1 '-' 1.7289) 'ABS') '/' 1.7289 ;
 'SI' (erro > 1.0D-2) ;
    'ERREUR' 5 ;
 'FINSI' ;
 
 p1 = 'EXTRAIRE' PN 'SCAL' PCEN1 ;
 
 erro = ((p1 '-' 0.3919) 'ABS') '/' 0.3919 ;
 'SI' (erro > 5.0D-2) ;
    'ERREUR' 5 ;
 'FINSI' ;
 
 m1 = 'EXTRAIRE' MACHN 'SCAL' PCEN1 ;
 
 erro = ((m1 '-' 1.44572) 'ABS') '/' 1.44572 ;
 'SI' (erro > 5.0D-2) ;
    'ERREUR' 5 ;
 'FINSI' ;
 
* Derrier le mach stem on a le mach min
* MMIN = 0.577350
* Il faut avoir un maillage assez fin pour l'avoir
 
 MMIN = 'MINIMUM' MACHN ;
 'MESSAGE' ('CHAINE' 'MMIN_ex =' 0.577350) ;
 'MESSAGE' ('CHAINE' 'MMIN =' MMIN) ;
 
 'SI' (MMIN > 1.0) ;
    'ERREUR' 5 ;
 'FINSI' ;
 'FIN' ;