Télécharger gravite.dgibi

Retour à la liste

* 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 - Tous droits réservés.
Mentions légales