* fichier :  precmat.dgibi
************************************************************
**** OPERATEUR 'KOPS'                                   ****
**** Matrice de preconditionnement écoulements bas mach ****
**** divisé par le pas de temps locale.                 ****
****                                                    ****
**** A. BECCANTINI, S. V. KUDRIAKOV,  LTMF              ****
**** NOVEMBRE    2001                                   ****
************************************************************

 'OPTION'  'ECHO' 1 ;
 'OPTION'  'DIME' 2 ;
 'OPTION'  'ELEM' 'TRI3' ;
 'OPTION'  'TRAC' 'X' ;

 GRAPH = FAUX ;
* GRAPH = VRAI ;
 
*
*** MAILLAGE
*
 
 P0 = 0.0D0 0.0D0 ;
 P1 = 3.0D0 0.0D0 ;
 P2 = 3.0D0 3.0D0 ;
 P3 = 0.0D0 3.0D0 ;
 P4 = 6.0D0 0.0D0 ;
 P5 = 6.0D0 3.0D0 ;

 RAF = 1 ;
 
 N1 = 2 '*' RAF ;
 N2 = 2 '*' RAF ;
 N3 = 2 '*' RAF ;
 N4 = 2 '*' RAF ;
 N5 = 1 '*' RAF ;
 N6 = 2 '*' RAF ;
 N7 = 2 '*' RAF ;

 LINEXT1 = ((P0 'DROIT'  N1 P1) 'ET' (P1 'DROIT'  N2 P2) 'ET' 
           (P2 'DROIT'  N3 P3) 'ET' (P3 'DROIT'  N4 P0)) ;
          
 LINEXT2 = ((P1 'DROIT'  N2 P2) 'ET' (P2 'DROIT'  N5 P5) 'ET' 
           (P5 'DROIT'  N6 P4) 'ET' (P4 'DROIT'  N7 P1)) ;          

 'OPTION' 'ELEM' QUA4 ;
 DOM1 = 'SURFACE'  LINEXT1 'PLAN' ;

 'OPTION' 'ELEM' TRI3 ;
 DOM2 = 'SURFACE'  LINEXT2 'PLAN' ;

 DOMTOT = DOM1 'ET' DOM2;
 'ELIMINATION' 1D-6 DOMTOT;

 CHPBRU = ('NOMC' 'UX' ('BRUI' 'BLAN' 'UNIF' 0.0  0.4 DOMTOT)
           'NATU' 'DISCRET') 'ET'
          ('NOMC' 'UY' ('BRUI' 'BLAN' 'UNIF' 0.0  0.2 DOMTOT)
           'NATU' 'DISCRET') ;
            
 'FORME' CHPBRU ;

 'SI' GRAPH ;
    'TRACER' DOMTOT ;
 'FINSI' ;

 $DOMTOT = 'DOMA' DOMTOT ;
 DIAMIN = 'DOMA' $DOMTOT 'DIAMIN' ;

*
**** 'KONV' pour la creation de la matrice identité
*

 LMOT = 'MOTS' 'R' 'UX' 'UY' 'RET' ;

 MAT1 = 'KOPS' 'MATIDE' LMOT ($DOMTOT . 'CENTRE') 'MATRIK' ;

 RN = ('BRUI' 'BLAN' 'UNIF' 1.0  0.4 ($DOMTOT . 'CENTRE')) ;
 PN = ('BRUI' 'BLAN' 'UNIF' 1.11  0.3  ($DOMTOT . 'CENTRE')) ;
 GAMN = ('BRUI' 'BLAN' 'UNIF' 1.4  0.1 ($DOMTOT . 'CENTRE')) ;
 CN = GAMN '*' (PN '/' RN) ;
 CN = CN '**' 0.5 ;

* Supersonic case
 
 VN = ('NOMC' 'UX' ('BRUI' 'BLAN' 'UNIF' 1.5  0.25 ($DOMTOT . 'CENTRE'))
       'NATU' 'DISCRET') 'ET'
      ('NOMC' 'UY' ('BRUI' 'BLAN' 'UNIF' 1.0  0.2  ($DOMTOT . 'CENTRE'))
       'NATU' 'DISCRET') ;

 VN = VN '*' CN ('MOTS' 'UX' 'UY') ('MOTS' 'SCAL' 'SCAL')
      ('MOTS' 'UX' 'UY') ;

* En ce cas particulier, MAT2 = matrice identité '/' dt
* et DT = DIAMIN '/' (QQ '+' CN)

 UINF = 'MANUEL' 'CHPO' ($DOMTOT . 'CENTRE') 1 'SCAL' 0.0 ;
 UINF1 = 'MANUEL' 'CHPO' ($DOMTOT . 'CENTRE') 1 'SCAL' 0.0 ;

 MAT2 = 'KONV' 'VF' 'PERFMONO' 'GAMMCONS'  LMOT ($DOMTOT . 'CENTRE')
         DIAMIN RN VN PN GAMN UINF UINF1 ;

 QQ = 'PSCAL' VN VN ('MOTS' 'UX' 'UY') ('MOTS' 'UX' 'UY') ;
 QQ = QQ ** 0.5 ;
 UNSDT =  (QQ '+' CN) '/' DIAMIN ;

 UTEST = ('MANUEL' 'CHPO' ($DOMTOT . 'CENTRE') 4
       'R' 11.0 'UX' 12.0 'UY' 13.0 'RET' 14.0 'NATU' 'DISCRET') ;
 UTEST1 = UTEST '*' UNSDT ;

 ERRO = 'MAXIMUM'
   (('KOPS'  'MULT' MAT1 UTEST1) '-' ('KOPS'  'MULT' MAT2 UTEST))
   'ABS'  ;
 
 'SI' (ERRO > 1.0D-12) ;
      'MESSAGE' 'Probleme en ET or en KOPS ' ;
      'ERREUR' 5 ;
 'FINSI' ;

 
* Dans de vitesse subsonique, on calcules les elts de la matrice dans
* l'elt3
 
 DOM3  = 'ELEM' DOMTOT 'QUA4' 3 ;
 $DOM3 = 'DOMA' DOM3  'INCL' $DOMTOT ;

 RN3 = 'MANUEL' 'CHPO' ($DOM3 . 'CENTRE') 1 'SCAL' 1.11 'NATU' 'DISCRET'
 ;
 PN3 = 'MANUEL' 'CHPO' ($DOM3 . 'CENTRE') 1 'SCAL' 7.17 'NATU' 'DISCRET'
 ;
 GAMN3 = 'MANUEL' 'CHPO' ($DOM3 . 'CENTRE') 1 'SCAL' 1.4 'NATU'
 'DISCRET' ;
 CN3 = GAMN3 '*' (PN3 '/' RN3) ;
 CN3 = CN3 '**' 0.5 ;
 VN3 = 'MANUEL' 'CHPO' ($DOM3 . 'CENTRE') 2 'UX' 0.1 'UY' 0.2
        'NATU' 'DISCRET' ;
 VN3 = VN3 '*' CN3 ('MOTS' 'UX' 'UY') ('MOTS' 'SCAL' 'SCAL')
      ('MOTS' 'UX' 'UY')  'NATU' 'DISCRET' ;

 RN = (RN  '-' ('REDU' RN ($DOM3 . 'CENTRE'))) '+' RN3 ;
 PN = (PN  '-' ('REDU' PN ($DOM3 . 'CENTRE'))) '+' PN3 ;
 VN = (VN  '-' ('REDU' VN ($DOM3 . 'CENTRE'))) '+' VN3 ;
 GAMN = (GAMN  '-' ('REDU' GAMN ($DOM3 . 'CENTRE'))) '+' GAMN3 ;

 
 MAT = 'KONV' 'VF' 'PERFMONO' 'GAMMCONS'  LMOT ($DOMTOT . 'CENTRE')
       DIAMIN  RN VN PN GAMN UINF  UINF1 ;

 UNZERO = ('MANUEL' 'CHPO' ($DOMTOT . 'CENTRE') 4
       'R' 0.0 'UX' 0.0 'UY' 0.0 'RET' 0.0 'NATU' 'DISCRET') ;
 UNRN3 = ('MANUEL' 'CHPO' ($DOM3 . 'CENTRE') 1
       'R' 1.0  'NATU' 'DISCRET') 'ET' UNZERO ;
 UNVX3 = ('MANUEL' 'CHPO' ($DOM3 . 'CENTRE') 1
       'UX' 1.0  'NATU' 'DISCRET') 'ET' UNZERO ;
 UNVY3 = ('MANUEL' 'CHPO' ($DOM3 . 'CENTRE') 1
       'UY' 1.0  'NATU' 'DISCRET') 'ET' UNZERO ;
 UNRETN3 = ('MANUEL' 'CHPO' ($DOM3 . 'CENTRE') 1
       'RET' 1.0  'NATU' 'DISCRET') 'ET' UNZERO ;

*
* \Gamma(1,1) = 1.19000000000000394
* \Gamma(1,2) = -0.252726909164141866
* \Gamma(1,3) = -0.505453818328283733
* \Gamma(1,4) = 0.840406455469216551
* \Gamma(2,1) = 0.571367728421105434E-01
* \Gamma(2,2) = 0.924000000000000044
* \Gamma(2,3) = -0.151999999999999941
* \Gamma(2,4) = 0.252726909164141866
* \Gamma(3,1) = 0.114273545684221087
* \Gamma(3,2) = -0.151999999999999941
* \Gamma(3,3) = 0.696000000000000063
* \Gamma(3,4) = 0.505453818328283733
* \Gamma(4,1) = 4.33849594594603438
* \Gamma(4,2) = -5.77081405705304729
* \Gamma(4,3) = -11.5416281141060946
* \Gamma(4,4) = 20.1899999999999906
*

* We have to divide them by DT
* In this case uref = qq
*
*
 QQ = 'PSCAL' VN VN ('MOTS' 'UX' 'UY') ('MOTS' 'UX' 'UY') ;
 QQ = QQ ** 0.5 ;
 qq = 'EXTRAIRE' QQ ('POIN' 1 ($DOM3 . 'CENTRE')) 'SCAL' ;
 CN = ((GAMN * PN) '/' RN) ** 0.5 ;
 a  = 'EXTRAIRE' CN ('POIN' 1 ($DOM3 . 'CENTRE')) 'SCAL' ;
 ur = qq ;
 mref=ur '/' a ;
 mref2=mref '*' mref ;
 mach=qq '/' a ;
 mach2=mach '*' mach ;
 lambda=(1.0D0 '-' mref2) '**' 2 ;
 lambda=(lambda '*' mach2) '+' (4.0D0 '*' mref2) ;
 lambda=(lambda '**' 0.5) * a ;
 lambda=( (1.0D0 '+' mref2) '*' qq) '+' lambda ;
 lambda=lambda '*' 0.5D0 ;
 dtau= ('EXTRAIRE' DIAMIN 'SCAL'
   ('POIN' 1 ($DOM3 . 'CENTRE'))) / lambda ;

 ERRTOL = 1.0D-8 ;

* 'XTY' UNR3 donne la premiére ligne

 GAMMA11 = ('XTY' UNRN3 ('KOPS' MAT 'MULT' UNRN3) LMOT LMOT) * dtau ;
 'SI' (('ABS' (GAMMA11 '-' 1.19000000000000)) '>' ERRTOL) ;
    'MESSAGE' ;
    'MESSAGE' 'GAMMA11' ;
    'MESSAGE' ;
    'ERREUR' 5 ;
 'FINSI' ;

 GAMMA12 = 'XTY' UNRN3 ('KOPS' MAT 'MULT' UNVX3) LMOT LMOT ;
 GAMMA12 = GAMMA12 * dtau ;
 'SI' (('ABS' (GAMMA12 '-' -0.25272690916414)) '>' ERRTOL) ;
    'MESSAGE' ;
    'MESSAGE' 'GAMMA12²' ;
    'MESSAGE' ;
    'ERREUR' 5 ;
 'FINSI' ;

 GAMMA13 = 'XTY' UNRN3 ('KOPS' MAT 'MULT' UNVY3) LMOT LMOT ;
 GAMMA13 = GAMMA13 * DTAU ;
 'SI' (('ABS' (GAMMA13 '-' -0.50545381832828)) '>' ERRTOL) ;
    'MESSAGE' ;
    'MESSAGE' 'GAMMA13' ;
    'MESSAGE' ;
    'ERREUR' 5 ;
 'FINSI' ;

 GAMMA14 = 'XTY' UNRN3 ('KOPS' MAT 'MULT' UNRETN3) LMOT LMOT ;
 GAMMA14 = GAMMA14 * DTAU ;
 'SI' (('ABS' (GAMMA14 '-' 0.840406455469217)) '>' ERRTOL) ;
    'MESSAGE' ;
    'MESSAGE' 'GAMMA14' ;
    'MESSAGE' ;
    'ERREUR' 5 ;
 'FINSI' ;

* 'XTY' UVNX3 donne la deuxieme ligne

 GAMMA21 = 'XTY' UNVX3 ('KOPS' MAT 'MULT' UNRN3) LMOT LMOT ;
 GAMMA21 = GAMMA21 * DTAU ;
 'SI' (('ABS' (GAMMA21 '-' 0.571367728421105434E-01)) '>' ERRTOL) ;
    'MESSAGE' ;
    'MESSAGE' 'GAMMA21' ;
    'MESSAGE' ;
    'ERREUR' 5 ;
 'FINSI' ;

 GAMMA22 = 'XTY' UNVX3 ('KOPS' MAT 'MULT' UNVX3) LMOT LMOT ;
 GAMMA22 = GAMMA22 * DTAU ;
 'SI' (('ABS' (GAMMA22 '-'  0.924000000000000044)) '>' ERRTOL) ;
    'MESSAGE' ;
    'MESSAGE' 'GAMMA22²' ;
    'MESSAGE' ;
    'ERREUR' 5 ;
 'FINSI' ;

 GAMMA23 = 'XTY' UNVX3 ('KOPS' MAT 'MULT' UNVY3) LMOT LMOT ;
 GAMMA23 = GAMMA23 * DTAU ;
 'SI' (('ABS' (GAMMA23 '-' -0.151999999999999941)) '>' ERRTOL) ;
    'MESSAGE' ;
    'MESSAGE' 'GAMMA23' ;
    'MESSAGE' ;
    'ERREUR' 5 ;
 'FINSI' ;

 GAMMA24 = 'XTY' UNVX3 ('KOPS' MAT 'MULT' UNRETN3) LMOT LMOT ;
 GAMMA24 = GAMMA24 * DTAU ;
 'SI' (('ABS' (GAMMA24 '-'  0.252726909164141866)) '>' ERRTOL) ;
    'MESSAGE' ;
    'MESSAGE' 'GAMMA24' ;
    'MESSAGE' ;
    'ERREUR' 5 ;
 'FINSI' ;

* 'XTY' UNVY3 donne la deuxieme ligne

 GAMMA31 = 'XTY' UNVY3 ('KOPS' MAT 'MULT' UNRN3) LMOT LMOT ;
 GAMMA31 = GAMMA31 * DTAU ;
 'SI' (('ABS' (GAMMA31 '-' 0.114273545684221087)) '>' ERRTOL) ;
    'MESSAGE' ;
    'MESSAGE' 'GAMMA31' ;
    'MESSAGE' ;
    'ERREUR' 5 ;
 'FINSI' ;

 GAMMA32 = 'XTY' UNVY3 ('KOPS' MAT 'MULT' UNVX3) LMOT LMOT ;
 GAMMA32 = GAMMA32 * DTAU ;
 'SI' (('ABS' (GAMMA32 '-' -0.151999999999999941)) '>' ERRTOL) ;
    'MESSAGE' ;
    'MESSAGE' 'GAMMA32' ;
    'MESSAGE' ;
    'ERREUR' 5 ;
 'FINSI' ;

 GAMMA33 = 'XTY' UNVY3 ('KOPS' MAT 'MULT' UNVY3) LMOT LMOT ;
 GAMMA33 = GAMMA33 * DTAU ;
 'SI' (('ABS' (GAMMA33 '-' 0.696000000000000063)) '>' ERRTOL) ;
    'MESSAGE' ;
    'MESSAGE' 'GAMMA33' ;
    'MESSAGE' ;
    'ERREUR' 5 ;
 'FINSI' ;

 GAMMA34 = 'XTY' UNVY3 ('KOPS' MAT 'MULT' UNRETN3) LMOT LMOT ;
 GAMMA34 = GAMMA34 * DTAU ;
 'SI' (('ABS' (GAMMA34 '-' 0.505453818328283733)) '>' ERRTOL) ;
    'MESSAGE' ;
    'MESSAGE' 'GAMMA34' ;
    'MESSAGE' ;
    'ERREUR' 5 ;
 'FINSI' ;


* 'XTY' UNVY3 donne la deuxieme ligne

 GAMMA41 = 'XTY' UNRETN3 ('KOPS' MAT 'MULT' UNRN3) LMOT LMOT ;
 GAMMA41 = GAMMA41 * DTAU ;
 'SI' (('ABS' (GAMMA41 '-' 4.33849594594603438)) '>' ERRTOL) ;
    'MESSAGE' ;
    'MESSAGE' 'GAMMA41' ;
    'MESSAGE' ;
    'ERREUR' 5 ;
 'FINSI' ;

 GAMMA42 = 'XTY' UNRETN3 ('KOPS' MAT 'MULT' UNVX3) LMOT LMOT ;
 GAMMA42 = GAMMA42 * DTAU ;
 'SI' (('ABS' (GAMMA42 '-' -5.77081405705304729)) '>' ERRTOL) ;
    'MESSAGE' ;
    'MESSAGE' 'GAMMA42' ;
    'MESSAGE' ;
    'ERREUR' 5 ;
 'FINSI' ;

 GAMMA43 = 'XTY' UNRETN3 ('KOPS' MAT 'MULT' UNVY3) LMOT LMOT ;
 GAMMA43 = GAMMA43 * DTAU ;
 'SI' (('ABS' (GAMMA43 '-' -11.5416281141060946)) '>' ERRTOL) ;
    'MESSAGE' ;
    'MESSAGE' 'GAMMA43' ;
    'MESSAGE' ;
    'ERREUR' 5 ;
 'FINSI' ;

 GAMMA44 = 'XTY' UNRETN3 ('KOPS' MAT 'MULT' UNRETN3) LMOT LMOT ;
 GAMMA44 = GAMMA44 * DTAU ;
 'SI' (('ABS' (GAMMA44 '-' 20.1899999999999906)) '>' ERRTOL) ;
    'MESSAGE' ;
    'MESSAGE' 'GAMMA44' ;
    'MESSAGE' ;
    'ERREUR' 5 ;
 'FINSI' ;

 
* Gamma with respect to primitive variables

 LMOT1 = LMOT ;
 LMOT2 = 'MOTS' 'RN' 'UX2' 'UY2' 'PN' ;

 MAT = 'KONV' 'VF' 'PERFMONO' 'GAMMCONS'  LMOT ($DOMTOT . 'CENTRE')
       DIAMIN   RN VN PN GAMN UINF  UINF1 ;

 MATGP =  'KONV' 'VF' 'PERFMONO' 'GAMMPRIM'  ($DOMTOT . 'CENTRE')
           LMOT1 LMOT2
           DIAMIN  RN VN PN GAMN UINF UINF1 ;

 MATJAC =  'PRIM' 'CONSPRIM'  ($DOMTOT . 'CENTRE')
           LMOT1 LMOT2
           RN VN PN GAMN ;

* MATGP = MAT2 * MATJAC
* 'KOPS' 'MULT' ne marche pas avec deux matrik
* Mais MAT2 * MATJAC * UN = MATGP * UN

 UN = ('MANUEL' 'CHPO' ($DOMTOT . 'CENTRE') 4
       'RN' 0.10 'UX2' 11.0 'UY2' 31.0 'PN' 41.0 'NATU' 'DISCRET') ;

 UN1 = 'KOPS' 'MULT' MATJAC UN ;

 UN2 = 'KOPS' 'MULT' MAT UN1 ;

 UN3 = 'KOPS' 'MULT' MATGP UN ;

 ERRO = 'MAXIMUM' (UN2 '-' UN3) 'ABS' ;

'SI' (ERRO > 1.0D-12) ;
      'MESSAGE' 'Probleme en ET or en KOPS ' ;
      'ERREUR' 5 ;
 'FINSI' ;

* What happens if UINF = CSON

 MAT1 = 'KOPS' 'MATIDE'  LMOT ($DOMTOT . 'CENTRE') 'MATRIK' ;

 CN = ((GAMN '*' PN) '/' RN) '**' 0.5 ;
 MACHN = 'PSCAL' VN VN ('MOTS' 'UX' 'UY') ('MOTS' 'UX' 'UY') ;
 MACHN = MACHN '**' 0.5 ;
 MACHN = MACHN '/' CN ;

 UINF = CN ;

 MAT2 = 'KONV' 'VF' 'PERFMONO' 'GAMMCONS'  LMOT ($DOMTOT . 'CENTRE')
         DIAMIN RN VN PN GAMN UINF  UINF1 ;
 
* UREF = UINF = CN   ->   DT = DIAMIN '/' (QQ '+' CN)

 QQ = 'PSCAL' VN VN ('MOTS' 'UX' 'UY') ('MOTS' 'UX' 'UY') ;
 QQ = QQ ** 0.5 ;
 UNSDT =  (QQ '+' CN) '/' DIAMIN ;

 UTEST = ('MANUEL' 'CHPO' ($DOMTOT . 'CENTRE') 4
       'R' 11.0 'UX' 12.0 'UY' 13.0 'RET' 14.0 'NATU' 'DISCRET') ;
 UTEST1 = UTEST '*' UNSDT ;

 ERRO = 'MAXIMUM'
   (('KOPS'  'MULT' MAT1 UTEST1) '-' ('KOPS'  'MULT' MAT2 UTEST))
   'ABS'  ;
 
 'SI' (ERRO > 1.0D-12) ;
      'MESSAGE' 'Probleme en ET or en KOPS ' ;
      'ERREUR' 5 ;
 'FINSI' ;

 
  
* What happens if UINF > CSON

 MAT1 = 'KOPS' 'MATIDE' LMOT ($DOMTOT . 'CENTRE') 'MATRIK' ;

 CN = ((GAMN '*' PN) '/' RN) '**' 0.5 ;
 MACHN = 'PSCAL' VN VN ('MOTS' 'UX' 'UY') ('MOTS' 'UX' 'UY') ;
 MACHN = MACHN '**' 0.5 ;
 MACHN = MACHN '/' CN ;

 UINF = 2. '*' CN ;

 MAT2 = 'KONV' 'VF' 'PERFMONO' 'GAMMCONS'  LMOT ($DOMTOT . 'CENTRE')
         DIAMIN RN VN PN GAMN UINF1 UINF  ;

 QQ = 'PSCAL' VN VN ('MOTS' 'UX' 'UY') ('MOTS' 'UX' 'UY') ;
 QQ = QQ ** 0.5 ;
 UNSDT =  (QQ '+' CN) '/' DIAMIN ;

 UTEST = ('MANUEL' 'CHPO' ($DOMTOT . 'CENTRE') 4
       'R' 11.0 'UX' 12.0 'UY' 13.0 'RET' 14.0 'NATU' 'DISCRET') ;
 UTEST1 = UTEST '*' UNSDT ;

 ERRO = 'MAXIMUM'
   (('KOPS'  'MULT' MAT1 UTEST1) '-' ('KOPS'  'MULT' MAT2 UTEST))
   'ABS'  ;
 
 'SI' (ERRO > 1.0D-12) ;
      'MESSAGE' 'Probleme en ET or en KOPS ' ;
      'ERREUR' 5 ;
 'FINSI' ;


* Gamma for primitive = jacobian for primitive 

 LMOT2 = LMOT ;
 LMOT1 = 'MOTS' 'RN' 'UX1' 'UY1' 'PN' ;

 MATGP =  'KONV' 'VF' 'PERFMONO' 'GAMMPRIM'  ($DOMTOT . 'CENTRE')
           LMOT1 LMOT2 DIAMIN
           RN VN PN GAMN UINF1 UINF ;

 MATJAC =  'PRIM' 'CONSPRIM' ($DOMTOT . 'CENTRE')
           LMOT1 LMOT2 
           RN VN PN GAMN ;

 MATR1 = MATJAC 'ET' ('KOPS' 'MULT' MATGP -1.0) ;
 
 UN = ('MANUEL' 'CHPO' ($DOMTOT . 'CENTRE') 4
       'RN' 1.0 'UX1' 1.0 'UY1' 1.0 'PN' 1.0 'NATU' 'DISCRET') ;
 
 RIS1 = ('KOPS' MATR1 'MULT' UN) ;

 ERRO = 'MAXIMUM' RIS1 'ABS' ;

'SI' (ERRO > 1.0D-12) ;
      'MESSAGE' 'Probleme en ET or en KOPS ' ;
      'ERREUR' 5 ;
 'FINSI' ;
 
 'FIN' ;