* 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' ;
© Cast3M 2003 - All rights reserved.
Disclaimer