* fichier : gravite.dgibi
'OPTI' 'ECHO' 0 ;
************************************************************************
* NOM : GRAVITE
* DESCRIPTION : Cas-test gravité servant à tester la méthode
* méthode de projection incrémentale.
* On doit trouver vitesse nulle et pression
* hydrostatique
*
* LANGAGE : GIBIANE-CAST3M
* AUTEUR : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
* mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION : v1, 21/11/2007, version initiale
* HISTORIQUE : v1, 21/11/2007, création
* HISTORIQUE :
* HISTORIQUE :
************************************************************************
* Prière de PRENDRE LE TEMPS de compléter les commentaires
* en cas de modification de ce sous-programme afin de faciliter
* la maintenance !
************************************************************************
'SAUTER' 2 'LIGNE' ;
'MESSAGE' ' Execution de gravite.dgibi' ;
'SAUTER' 2 'LIGNE' ;
*
interact = FAUX ;
graph = FAUX ;
*
*
* Procédure pour afficher la résultante d'un chpoint et de sa valeur
* absolue (utile pour debugger sans lister le CHPOINT en question)
*
*'DEBPROC' HCHPO ;
*'ARGUMENT' chp*'CHPOINT' ;
*'ARGUMENT' tit*'MOT' ;
**
*achp = 'ABS' chp ;
*nachp = 'MAXIMUM' achp ;
**'SI' ('<EG' nachp 1.D-60) ; nachp= 1.D0 ; 'FINSI' ;
*nachp = 1.D0 ;
*cc = '/' achp nachp ;
**hh = '+' cc 1. ;
*hh = 'ABS' cc ;
*rcc = 'RESULT' cc ;
*rhh = 'RESULT' hh ;
*mrrcc = 'MAXIMUM' rcc 'ABS' ;
*'MESSAGE' ('CHAINE' 'Resultante de ' tit ' ' mrrcc) ;
**'LISTE' rcc ;
**'LISTE' rhh ;
*'FINPROC' ;
*
'DEBP' TRACVP;
'ARGUMENT' RVX*'TABLE';
rv = rvx . 'EQEX' ;
un = rv . 'INCO' . 'UN' ;
pn = rv . 'INCO' . 'PN' ;
maxun = 'MAXIMUM' un 'ABS' ;
echvit = 'MAXIMUM' ('PROG' maxun 1.D-10) ;
cm = '/' 1. ('*' nmail echvit) ;
vn = 'VECT' un cm 'UX' 'UY' 'JAUN' ;
ppn = 'ELNO' $mt pn kpre ;
tit = 'CHAINE' 'V-P echvit=' echvit ;
'TRACER' ppn mt cmt vn 'TITR' tit 'NCLK' ;
as2 ama1 = 'KOPS' 'MATRIK';
'RESPRO' as2 ama1;
'FINP';
*
'DEBP' MULDT ;
'ARGUMENT' RVX*'TABLE';
rv = rvx . 'EQEX' ;
fdt = rvx . 'ARG1' ;
rv . 'INCO' . 'DT' = '*' (rv . 'INCO' . 'DT') fdt ;
'MESSAGE' ('CHAINE' 'dt = ' (rv . 'INCO' . 'DT')) ;
as2 ama1 = 'KOPS' 'MATRIK';
'RESPRO' as2 ama1;
'FINP';
*
*
'OPTION' 'DIME' 2 ;
'SI' ('NON' interact) ;
'OPTION' 'TRAC' 'PS' ;
'SINON' ;
'OPTION' 'TRAC' 'X' ;
'FINSI' ;
*
* Paramètres du calcul
*
ttit = 'TABLE' ;
tresu = 'TABLE' ;
tresp = 'TABLE' ;
tresd = 'TABLE' ;
tymesh = 3 ;
kvit = 'QUAF' ;
kpre = 'CENTREP1' ;
Re = 1. ;
dtini = 0.01 ;
fdt = 0.98 ;
*fdt = 0.8 ;
'REPETER' imesh 2 ;
'SI' ('EGA' &imesh 1) ;
tymesh = 3 ;
'MESSAGE' '**************' ;
'MESSAGE' 'Maillage carré' ;
'MESSAGE' '**************' ;
'FINSI' ;
'SI' ('EGA' &imesh 2) ;
tymesh = 4 ;
'MESSAGE' '**************' ;
'MESSAGE' 'Maillage rond ' ;
'MESSAGE' '**************' ;
'FINSI' ;
*
* Maillage
*
'SI' ('EGA' tymesh 1) ;
'OPTI' 'ELEM' 'QUA8' ;
nmail = 2 ;
'FINSI' ;
'SI' ('EGA' tymesh 2) ;
'OPTI' 'ELEM' 'QUA8' ;
nmail = 5 ;
'FINSI' ;
'SI' ('EGA' tymesh 3) ;
'OPTI' 'ELEM' 'TRI6' ;
nmail = 10 ;
'FINSI' ;
'SI' ('EGA' tymesh 4) ;
'OPTI' 'ELEM' 'QUA8' 'MODE' 'PLAN' ;
nmail = 10 ;
'FINSI' ;
p0 = 0. 0. ; p1 = 1. 0. ;
lt = 'DROIT' nmail p0 p1 ;
mt = 'TRANSLATION' lt nmail (0. 1.) ;
*mt = 'TRANSLATION' lt nmail (0. iang) ;
cmt = 'CONTOUR' mt ;
bas dro hau gau = 'COTE' mt ;
'SI' ('EGA' tymesh 4) ;
xhau = 'COORDONNEE' 1 hau ;
sxhau = '*' ('SIN' ('*' xhau 180.)) 0.2 ;
cshau = 'NOMC' 'UY' sxhau ;
'FORME' cshau ;
'FINSI' ;
'SI' ('OU' ('EGA' tymesh 3) ('EGA' tymesh 4)) ;
mt = 'SURFACE' cmt ;
'FINSI' ;
*'TRACER' mt ;
hau = 'INVERSE' hau ;
_bas = 'CHANGER' bas 'QUAF' ; _dro = 'CHANGER' dro 'QUAF' ;
_hau = 'CHANGER' hau 'QUAF' ; _gau = 'CHANGER' gau 'QUAF' ;
_mt = 'CHANGER' mt 'QUAF' ;
'ELIMINATION' (_mt 'ET' _bas 'ET' _gau 'ET' _dro 'ET' _hau) 1.D-6 ;
$mt = 'MODE' _mt 'NAVIER_STOKES' kvit ;
$gau = 'MODE' _gau 'NAVIER_STOKES' kvit ;
$dro = 'MODE' _dro 'NAVIER_STOKES' kvit ;
$hau = 'MODE' _hau 'NAVIER_STOKES' kvit ;
$bas = 'MODE' _bas 'NAVIER_STOKES' kvit ;
mt = 'DOMA' $mt 'MAILLAGE' ; gau = 'DOMA' $gau 'MAILLAGE' ;
dro = 'DOMA' $dro 'MAILLAGE' ; hau = 'DOMA' $hau 'MAILLAGE' ;
bas = 'DOMA' $bas 'MAILLAGE' ;
cmt = bas 'ET' dro 'ET' hau 'ET' gau ;
mpres = 'DOMA' $mt kpre ;
mp1 = 'ELEM' mpres 1 ;
pp1 = 'POIN' mp1 1 ;
* Calcul de pexact
pexact = '-' ('COORDONNEE' 2 mpres) ('COORDONNEE' 2 pp1) ;
*
lok = VRAI ;
*
'REPETER' iaxi 2 ;
'SAUTER' 2 'LIGN' ;
'SI' ('EGA' &iaxi 1) ;
'OPTI' 'MODE' 'PLAN' ;
'MESSAGE' ' *********' ;
'MESSAGE' ' Mode plan' ;
'MESSAGE' ' *********' ;
'FINSI' ;
'SI' ('EGA' &iaxi 2) ;
'OPTI' 'MODE' 'AXIS' ;
'MESSAGE' ' ********' ;
'MESSAGE' ' Mode axi' ;
'MESSAGE' ' ********' ;
'FINSI' ;
nmet = 5 ;
'REPETER' imet nmet ;
*GOON'REPETER' imet 3 ;
'SAUTER' 1 'LIGN' ;
'SI' ('EGA' &imet 1) ;
tit = 'CHAINE' 'Methode directe' ;
mproj = FAUX ; lprec = FAUX ;
'FINSI' ;
'SI' ('EGA' &imet 2) ;
tit = 'CHAINE' 'Methode projection' ;
mproj = VRAI ; lprec = FAUX ; doubl = FAUX ;
'FINSI' ;
'SI' ('EGA' &imet 3) ;
tit= 'CHAINE' 'Methode projection preconditionnee' ;
mproj = VRAI ; lprec = VRAI ; doubl = FAUX ;
'FINSI' ;
'SI' ('EGA' &imet 4) ;
tit = 'CHAINE' 'Methode double projection' ;
mproj = VRAI ; lprec = FAUX ; doubl = VRAI ;
'FINSI' ;
'SI' ('EGA' &imet 5) ;
tit = 'CHAINE' 'Methode double projection preconditionnee' ;
mproj = VRAI ; lprec = VRAI ; doubl = VRAI ;
'FINSI' ;
'MESSAGE' tit ;
'SAUTER' 1 'LIGN' ;
ttit . &imet = tit ;
*
* table EQEX
*
'SI' mproj ;
'SI' ('NON' doubl) ;
* Simple proj
*GOON
* niter = 3 ;
niter = 5 ;
'SINON' ;
* Double proj : 1 suffit mais on met 3 car ca teste
* le préconditionnement du calcul de la
* matrice de pression
* niter = 3 ;
niter = 3 ;
'FINSI' ;
'SINON' ;
niter = 1 ;
'FINSI' ;
rv = 'EQEX' 'NITER' niter 'OMEGA' 1. 'ITMA' 1 ;
'SI' graph ;
rv = 'EQEX' rv 'ZONE' $mt 'OPER' TRACVP ;
'FINSI' ;
rv = 'EQEX' rv
'OPTI' 'EF' 'IMPL' 'CENTREE'
'ZONE' $mt 'OPER' 'DFDT' 1. 'UNM1' 'DT' 'INCO' 'UN'
'ZONE' $mt 'OPER' MULDT fdt
'OPTI' 'EF' 'IMPL' 'CENTREE'
'ZONE' $mt 'OPER' 'NS' 1. 'UN' ('/' 1. Re)
(0. -1.) 'TN' 0. 'INCO' 'UN'
'OPTI' 'EF' 'IMPL' 'CENTREE' kpre
'ZONE' $mt 'OPER' 'KBBT' 1. 'INCO' 'UN' 'PN'
;
rv = 'EQEX' rv
'CLIM' 'UN' 'UIMP' cmt 0.
'CLIM' 'UN' 'VIMP' cmt 0.
'CLIM' 'PN' 'TIMP' mp1 0.
;
rv . 'INCO' = 'TABLE' 'INCO' ;
rv . 'INCO' . 'UNM1' = 'KCHT' $mt 'VECT' 'SOMMET' (0. 0.) ;
rv . 'INCO' . 'UN' = 'KCHT' $mt 'VECT' 'SOMMET' (0. 0.) ;
rv . 'INCO' . 'PN' = 'KCHT' $mt 'SCAL' kpre 0. ;
rv . 'INCO' . 'TN' = 'KCHT' $mt 'SCAL' 'SOMMET' 1. ;
rv . 'INCO' . 'DT' = dtini ;
*
'SI' mproj ;
rv . 'GPROJ' = 'TABLE' ;
rv . 'GPROJ' . 'NOMVIT' = 'UN' ;
rv . 'GPROJ' . 'NOMPRES' = 'PN' ;
'SI' ('NON' lprec) ;
rv . 'GPROJ' . 'NOPREC' = VRAI ;
'FINSI' ;
'SI' ('NON' doubl) ;
rv . 'GPROJ' . 'dblproj' = doubl ;
'FINSI' ;
'FINSI' ;
*
EXEC rv ;
*
un = rv . 'INCO' . 'UN' ;
pn = '-' (rv . 'INCO' . 'PN') pexact ;
'SI' graph ;
bidon matkbbt = 'KBBT' (rv . '5KBBT') ;
'SINON' ;
bidon matkbbt = 'KBBT' (rv . '4KBBT') ;
'FINSI' ;
nomavv = 'MOTS' 'UX' 'UY' ;
nomapv = 'MOTS' '1UN' '2UN' ;
nomavp = 'MOTS' 'SCAL' ;
nomapp = 'MOTS' 'PN' ;
vit = 'NOMC' nomavv nomapv un ;
fordiv = 'KOPS' matkbbt '*' vit ;
*'LISTE' ('EXTRAIRE' fordiv 'COMP') ;
*'LISTE' ('MAXIMUM' fordiv 'ABS') ;
*'LISTE' ('MAXIMUM' ('RESULT' ('ABS' fordiv))) ;
tresu . &imet = un ;
tresp . &imet = pn ;
tresd . &imet = fordiv ;
*
maxun = 'MAXIMUM' un 'ABS' ;
echvit = 'MAXIMUM' ('PROG' maxun 1.D-10) ;
cm = '/' 1. ('*' nmail echvit) ;
*vn = 'VECT' un cm 'DEPL' 'JAUN' ;
vn = 'VECT' un cm 'UX' 'UY' 'JAUN' ;
ppn = 'ELNO' $mt pn kpre ;
*
'SI' graph ;
titg = 'CHAINE' 'Erreur V-P echvit=' echvit ' mproj=' mproj
' lprec=' lprec ;
'TRACER' ppn mt mt vn 'TITR' titg ;
* 'TRACER' ('ELNO' $mt (rv . 'INCO' . 'PN') kpre) mt mt 'TITR' tit ;
* 'TRACER' ('ELNO' $mt pexact kpre) mt mt 'TITR' tit ;
'FINSI' ;
'FIN' imet ;
*
* Tests
*
*
* Tolérance pour chaque méthode
*
'SI' ('EGA' &imesh 1) ;
'SI' ('EGA' &iaxi 1) ;
lerru = 'PROG' 1.D-12 1.D-4 1.D-12 1.D-12 1.D-12 ;
lerrp = 'PROG' 2.D-12 4.D-2 1.D-12 2.D-12 1.D-12 ;
lerrd = 'PROG' 1.D-14 1.D-14 1.D-14 1.D-14 1.D-14 ;
'FINSI' ;
'SI' ('EGA' &iaxi 2) ;
lerru = 'PROG' 1.D-12 3.D-4 1.D-12 1.D-12 1.D-12 ;
lerrp = 'PROG' 2.D-12 5.D-2 1.D-12 2.D-12 1.D-10 ;
lerrd = 'PROG' 1.D-14 1.D-14 1.D-14 1.D-14 1.D-14 ;
'FINSI' ;
'FINSI' ;
'SI' ('EGA' &imesh 2) ;
'SI' ('EGA' &iaxi 1) ;
lerru = 'PROG' 1.5D-4 1.5D-4 1.D-12 1.5D-4 1.D-12 ;
lerrp = 'PROG' 1.D-2 4.0D-2 1.D-12 1.D-2 1.D-12 ;
lerrd = 'PROG' 1.D-14 1.D-14 1.D-14 1.D-14 1.D-14 ;
'FINSI' ;
'SI' ('EGA' &iaxi 2) ;
lerru = 'PROG' 1.5D-4 5.D-4 1.D-12 1.5D-4 1.D-12 ;
lerrp = 'PROG' 1.D-2 4.0D-2 1.D-12 1.D-2 1.D-12 ;
lerrd = 'PROG' 3.D-10 1.D-9 1.D-9 3.D-10 3.D-10 ;
'FINSI' ;
'FINSI' ;
*
*
*
'SAUTER' 1 'LIGN' ;
*
'REPETER' imet nmet ;
und = tresu . &imet ;
'SI' ('OU' ('EGA' &imet 3) ('EGA' &imet 5)) ;
undd = tresu . ('-' &imet 1) ;
errud = ('MAXIMUM' ('-' und undd) 'ABS') ;
vv = 'MAXIMUM' undd 'ABS' ;
'SI' ('>' vv 1.D-10) ; errud = '/' errud vv ; 'FINSI' ;
tit = 'CHAINE' (ttit . &imet) '/non prec' ;
'SINON' ;
errud = 'MAXIMUM' und 'ABS' ;
tit = 'CHAINE' (ttit . &imet) ;
'FINSI' ;
chmes = 'CHAINE' tit ' max err. vitesse =' errud ;
'MESSAGE' chmes ;
erruda = 'EXTRAIRE' lerru &imet ;
tst = ('<' errud erruda) ;
'SI' ('NON' tst) ;
chmes = 'CHAINE' '!!! Erreur, pas inferieure a ' erruda ;
'MESSAGE' chmes ;
'FINSI' ;
lok = lok 'ET' tst ;
*
pnd = tresp . &imet ;
'SI' ('OU' ('EGA' &imet 3) ('EGA' &imet 5)) ;
pndd = tresp . ('-' &imet 1) ;
errpd = 'MAXIMUM' ('-' pnd pndd) 'ABS' ;
vv = 'MAXIMUM' pndd 'ABS' ;
'SI' ('>' vv 1.D-10) ; errpd = '/' errpd vv ; 'FINSI' ;
tit = 'CHAINE' (ttit . &imet) '/non prec' ;
'SINON' ;
errpd = 'MAXIMUM' pnd 'ABS' ;
tit = 'CHAINE' (ttit . &imet) ;
'FINSI' ;
chmes = 'CHAINE' tit ' max err. pression =' errpd ;
'MESSAGE' chmes ;
errpda = 'EXTRAIRE' lerrp &imet ;
tst = ('<' errpd errpda) ;
'SI' ('NON' tst) ;
chmes = 'CHAINE' '!!! Erreur, pas inferieure a ' errpda ;
'MESSAGE' chmes ;
'FINSI' ;
lok = lok 'ET' tst ;
*
dnd = tresd . &imet ;
errdd = 'MAXIMUM' dnd 'ABS' ;
chmes = 'CHAINE' (ttit . &imet) ' max err. divergence =' errdd ;
'MESSAGE' chmes ;
errdda = 'EXTRAIRE' lerrd &imet ;
tst = ('<' errdd errdda) ;
'SI' ('NON' tst) ;
chmes = 'CHAINE' '!!! Erreur, pas inferieure a ' errdda ;
'MESSAGE' chmes ;
'FINSI' ;
lok = lok 'ET' tst ;
'FIN' imet ;
*
'FIN' iaxi ;
'FIN' imesh ;
*
'SAUTER' 2 'LIGNE' ;
'SI' lok ;
'MESSAGE' 'Tout sest bien passe' ;
'SINON' ;
'MESSAGE' 'Il y a eu des erreurs' ;
'FINSI' ;
'SAUTER' 2 'LIGNE' ;
*
'SI' interact ;
'OPTION' 'DONN' 5 'ECHO' 1 ;
'FINSI' ;
'SI' ('NON' lok) ;
'ERREUR' 5 ;
'FINSI' ;
*
* End of dgibi file GRAVITE
*
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer