Télécharger gresp.procedur

Retour à la liste

* GRESP     PROCEDUR  GOUNAND   12/02/01    21:15:14     7266           
************************************************************************
* NOM         : GRESP
* DESCRIPTION : Résout un système par une méthode de projection
*               algébrique incrémentale
*
* LANGAGE     : GIBIANE-CAST3M
* AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
*               mél : gounand@semt2.smts.cea.fr
**********************************************************************
* VERSION    : v1, 22/11/2007, version initiale
* HISTORIQUE : v1, 22/11/2007, création
* HISTORIQUE : 21/12/2011 : correction d'un bug dans le
*              préconditionnement, nettoyage KOPS
* HISTORIQUE :
************************************************************************
*
'DEBPROC' GRESP ;
'ARGUMENT'  mat*'MATRIK' ;
'ARGUMENT' ccli*'CHPOINT' ;
'ARGUMENT' smbr*'CHPOINT' ;
'ARGUMENT' rv*'TABLE' ;
*
* Nom des inconnues : vitesses et pressions
* Dans rvg . 'METHINV', on stocke l'éventuelle table d'inversion
* pour le laplacien de pression
*
debug = FAUX ;
rvg = rv . 'GPROJ' ;
rvm = rv . 'METHINV' ;
'SI' ('EXISTE' rvg 'dblproj') ;
   dblproj = rvg . 'dblproj' ;
'SINON' ;   
   dblproj = VRAI ;
'FINSI' ;   
*
'SI' ('EXISTE' rvg 'NOPREC') ;
   noprec = rvg . 'NOPREC' ;
'SINON' ;
   noprec = faux ;
'FINSI' ;
'SI' debug ;
   'MESSAGE' ('CHAINE' 'noprec=' noprec) ;
'FINSI' ;   
*
'SI' ('NON' noprec) ;
   'SI' ('NON' ('EXISTE' rvg 'preco')) ;
      rvg . 'preco' = 'TABLE' ;
   'FINSI' ;
   rvgp = rvg . 'preco' ;
* Précision relative utilisée pour le préconditionnement du calcul de la
* matrice de pression
   precrel = 1.D-12 ;
*   precrel = 1.D-2 ;
'FINSI' ;   
* Projection incrémentale donc pas besoin de XINIT
'SI' ('EXIS' rvm 'XINIT' ) ;
   'OUBLIER' rvm 'XINIT' ;
'FINSI' ;
*fd 'OUBLIER' rvm 'XINIT' ;
* On a mis chaine sinon plantage qd la pression s'appelle 'PRES' ;
ngvit = 'CHAINE' rvg . 'NOMVIT' ;
ngpre = 'CHAINE' rvg . 'NOMPRES' ;
'SI' ('NON' ('EXISTE' rvg 'METHINV')) ;
   rvg . 'METHINV' = rvm ;
'FINSI' ;
rvgm = rvg . 'METHINV' ;
*
dim = 'VALEUR' 'DIME' ;
*fd
* Si l'utilisateur a transmis les noms de composantes on ne fait rien.
* Sinon (il utilise EXEC) on construit la liste nivit = 'MOTS' '1UN' '2UN'
'SI' ('EXIS' rvg 'COMPVIT') ;
   nivit  = rvg . 'COMPVIT' ;
   nivit2 = nivit           ;
   nipre  = rvg . 'COMPPRES' ;
   nipre2 = nipre ;
'SINON' ;
nivit = 'MOTS' ;
'REPETER' idim dim ;
   niviti = 'CHAINE' &idim ngvit ;
   nivit = 'ET' nivit ('MOTS' niviti) ;
'FIN' idim ;   
nipre = 'MOTS' ngpre ;
nivit2 = 'EXTRAIRE' ('MOTS' 'UX' 'UY' 'UZ') ('LECT' 1 'PAS' 1 dim) ; 
nipre2 = 'MOTS' 'SCAL' ;
'FINSI' ;
* Construit la liste des noms d'inconnues ni vitesse ni pression
nitout = 'EXTR' mat 'COMP' ;
niaut = 'MOTS' ;
'REPETER' itout ('DIME' nitout) ;
   ninc = 'EXTRAIRE' nitout &itout ;
   lexis = FAUX ;
   'REPETER' ivit ('DIME' nivit) ;
      lexis = 'OU' lexis ('EGA' ninc ('EXTRAIRE' nivit &ivit)) ;
   'FIN' ivit ;      
   'REPETER' ipre ('DIME' nipre) ;
      lexis = 'OU' lexis ('EGA' ninc ('EXTRAIRE' nipre &ipre)) ;
   'FIN' ipre ;
   'SI' ('NON' lexis) ;
      niaut = 'ET' niaut ('MOTS' ninc) ;
   'FINSI' ;
'FIN' itout ;
*
* Initialisation de la solution du système total
*
sol = 0. ;
*         
* On résout d'abord la partie ni vitesse ni pression 
*
'SI' ('>' ('DIME' niaut) 0) ;
   mtot  = 'EXTRAIRE' mat niaut niaut ;
   ftot  = 'EXCO' niaut smbr niaut 'NOID' ;
   cltot = 'EXCO' niaut ccli niaut 'NOID' ;
   rvm . 'MATASS' =mtot ;
   rvm . 'MAPREC' =mtot ;
   ct    = 'KRES' mtot ftot 'CLIM' cltot 'TYPI' rvm ;
   sol   = '+' ct sol ;
'FINSI' ;
********************************************************************
*                                                                  *
* On résout la partie vitesse-pression par projection incrémentale *
*                                                                  *
********************************************************************
nivp = 'ET' nivit nipre ;
* Réduction des arguments de la proc à la vitesse-pression
*mat  = 'KOPS' 'EXTRINCO' mat nivp nivp ;
mat  = 'EXTRAIRE' mat nivp nivp ;
ccli = 'EXCO' nivp ccli nivp 'NOID' ;
smbr = 'EXCO' nivp smbr nivp 'NOID' ;
*
un = rv . 'INCO' . ngvit ;
pn = rv . 'INCO' . ngpre ;
un2 = 'NOMC' nivit2 nivit un ;
pn2 = 'NOMC' nipre2 nipre pn ;
incovp = un2 '+' pn2 ;
vpetit = -1.D200 ;
* Comme on va faire une méthode incrémentale, il faut imposer
* les conditions aux limites de Dirichlet sur l'inconnue
ccli1 = 'MASQUE' ccli 'SUPERIEUR' vpetit ;
ccli0 = '-' ccli1 1.D0 ;
niccli = 'EXTRAIRE' ccli 'COMP' ;
dcli  = '-' ccli ('*' incovp ccli1 niccli niccli niccli) ;
incovp  = '+' incovp dcli ;
*
* Matrice : partie diagonale en vitesse et contrainte
*           on est obligé de reconstruire la transposée car KBBT
*           ne stocke que B et pas Bt
*
mkvit    = 'EXTRAIRE' mat nivit nivit ;
mkcontr  = 'EXTRAIRE' mat nivit nipre ;
mkcontrt = 'KOPS' 'TRANSPOS' mkcontr ;
Kt = mkvit 'ET' mkcontr 'ET' mkcontrt ;
*
* Construction du résidu * -1
*
residu = '-' smbr ('*' Kt incovp) ;
*
* Calcul de la matrice de pression et de la matrice masse diagonalisée
*
*  On lumpe la partie diagonale et on l'inverse
chpo1   = 'MASQUE' un2 'SUPERIEUR' vpetit ;
mkvitd  = '*' mkvit chpo1 ;
*
* Choix 1 : si la matrice masse lumpée
*           a tous ses termes positifs, on en prend l'inverse
*           Avantage : simple, autorise le préconditionnement
*                      car ne change pas à chaque itération
*           Inconvénient : ne marche pas en axi quadratique sur l'axe
* Choix 2 : on prend la diagonale de la matrice totale
*           on la met à l'échelle pour avoir la même masse que mkvitd
*           Avantage : marche en axi, sur maillage déformés
*           Inconvénient : change à chaque itération quand le problème
*           est non linéaire (=> recalcul de la matrice de pression),
*           pas forcément positif.
* Choix 3 : valeur absolue du précédent + petite tolérance
*           Avantage : marche peut-être
*           Inconvénient : change à chaque itération.
*          
*
mmi = 'MINIMUM' mkvitd ;
'SI' ('>' mmi 1.D-200) ;
   imkvitd1 = 'INVERSE' mkvitd ;
   'SI' debug ; 'MESSAGE' 'Choix 1 masse lumpée' ; 'FINSI' ;
'SINON' ;
   mkdia = 'EXTR' mkvit 'DIAG' ;
   mmi = 'MINIMUM' mkdia ;
   'SI' ('>' mmi 1.D-200) ;
      'SI' debug ; 'MESSAGE' 'Choix 2 masse lumpée' ; 'FINSI' ;
   'SINON' ;
      mkdia = 'ABS' mkdia ;
      mma = 'MAXIMUM' mkdia ;
      mkdia = '+' mkdia ('*' mma 1.D-8) ;
      'SI' debug ; 'MESSAGE' 'Choix 3 masse lumpée' ; 'FINSI' ;
   'FINSI' ;      
   masvitd = 'MAXIMUM' ('RESULT' mkvitd) ;
*   'MESSAGE' ('CHAINE' 'totmas=' masvitd) ;
   masdia = 'MAXIMUM' ('RESULT' mkdia) ;
   imkvitd1 = '*' ('INVERSE' mkdia) ('/' masdia masvitd) ;
'FINSI' ;   
*  On surcharge par 0. l'inverse de la diagonale là où il y a des clims
*  de Dirichlet
imkvitd  = '-' imkvitd1 ('*' imkvitd1 ccli1 niccli niccli niccli) ;
*
* Attention ! On préconditionne éventuellement le calcul de la matrice
* de pression 
*
precok = FAUX ;
*'MESSAGE' ('CHAINE' 'precok = ' precok) ;
'SI' ('NON' noprec) ;
   exis1 = 'EXISTE' rvgp 'mklapphi' ;
   exis2 = 'EXISTE' rvgp 'imkvitd'   ;
   'SI' (exis1 'ET' exis2) ;
      imkvd = rvgp . 'imkvitd' ;
      mkl   = rvgp . 'mklapphi' ;
* On vérifie l'égalité des mkvitd à un facteur constant près
      ech  = 'MAXIMUM' imkvitd 'ABS' ;
      echp = 'MAXIMUM' imkvd 'ABS' ;
      alfa = '/' ech echp ;
*   'MESSAGE' ('CHAINE' 'alfa = ' alfa) ;
      dimkv = '-' imkvitd ('*' imkvd alfa) ;
      ndimkv = 'EXTRAIRE' dimkv 'COMP' ;
      inimkv = 'INVERSE' ('+' imkvitd ccli1) ;
      dimkvs = '*' dimkv inimkv ndimkv ndimkv ndimkv ;
      mdimkvs = 'MAXIMUM' dimkvs 'ABS' ;
      egaimkv = 'EGA' mdimkvs 0. precrel ;
*
      'SI' debug ;
         'MESSAGE' ('CHAINE' 'egaimkv = ' egaimkv
                       '  a ' ('*' mdimkvs 100.) ' %' ) ;
      'FINSI' ;         
      precok = egaimkv ;
   'FINSI' ;
'FINSI' ;   
*
chdbg = 'CHAINE' 'Matrice de pression MP : ' ;
'SI' precok ;
   'SI' ('NON' ('EGA' alfa 1.D0 precrel)) ;
      chdbg = 'CHAINE' chdbg 'MP(i) = MP(i-1) * ' alfa ;
*      mklapphi = '*' (rvgp . 'mklapphi') alfa ;
*      imkvitd  = '*' (rvgp . 'imkvitd') alfa ;
      mklapphi = rvgp . 'mklapphi' ;
      imkvitd  = rvgp . 'imkvitd' ;
   'SINON' ;
      chdbg = 'CHAINE' chdbg 'MP(i) = MP(i-1)'  ;
      mklapphi = rvgp . 'mklapphi' ;
      imkvitd  = rvgp . 'imkvitd' ;
   'FINSI' ;
'SINON' ;   
* Dans le cas général où mkcontr et mkcontrt sont différents,
* il faudrait mettre :
*mklapphi = 'KOPS' 'CMCT' mkcontr imchd (kops mkcontrt 'TRANSPOS') ;
   chdbg = 'CHAINE' chdbg 'calcul...' ;
   mklapphi = 'KOPS' 'CMCT' mkcontr imkvitd mkcontr ;
   alfa = 1.D0 ;
'FINSI' ;
'SI' debug ; 'MESSAGE' chdbg ; 'FINSI' ;
'SI' ('NON' noprec) ;
   rvgp . 'mklapphi'  = mklapphi ;
   rvgp . 'imkvitd'   = imkvitd  ;
'FINSI' ;
*'MESSAGE' 'Fin du calcul de la matrice de pression' ;
*
* Calcul éventuel d'une estimation de la pression (double projection)
*
*'MESSAGE' ('CHAINE' 'dblproj=' dblproj) ;
'SI' dblproj ;
   desqdm = 'EXCO' nivit residu nivit 'NOID' ;
   idesq  = '*' desqdm imkvitd nivit nivit nivit ;
* 1) Celui-ci se compense avec le 2)
*    donc pas nécessaire
*   'SI' ('NEG' alfa 1.D0 1.D-13) ;
*      idesq = '*' idesq alfa ;
*   'FINSI' ;
   didesq = '*' mkcontr idesq ;
   mtot  = mklapphi ;
*dbg   dmtot = 'EXTRAIRE' mklapphi 'DIAG' ;
*dbg   HCHPO dmtot 'dmtot' ;
   ftot  = didesq ;
   cltot = 'EXCO' nipre ccli0 nipre 'NOID' ;
   rvgm .'MATASS' =mtot ;
   rvgm . 'MAPREC' =mtot ;
*dbg   dphi1  tt = 'KRES' mtot ftot 'CLIM' cltot 'TYPI' rvgm
*dbg      'IMPR' 3 'LTIME' VRAI ;
   dphi1 = 'KRES' mtot ftot 'CLIM' cltot 'TYPI' rvgm ;
* 2) Se compense avec le 1)
*    donc pas nécessaire
*   'SI' ('NEG' alfa 1.D0 1.D-13) ;
*      dphi1 = '/' dphi1 alfa ;
*   'FINSI' ;
*   HCHPO dphi1 'dphi1' ;
   gdphi1 = '*' mkcontrt dphi1 ;
   residu = residu '-' gdphi1 ;
'FINSI' ;
*
* Calcul de l'incrément de vitesse dv* à partir de la qdm
* On splitte les équations par direction
*
dvitstar = 0. ;
ncomp = 'DIME' nivit ;
'REPETER' iicomp ncomp ;
   icomp = &iicomp ;
   licomp = 'LECT' icomp ;
   mcomp  = 'EXTRAIRE' nivit licomp ;
   mtot  = 'EXTRAIRE' mkvit mcomp mcomp ;
   ftot  = 'EXCO' mcomp residu mcomp 'NOID' ;
   cltot = 'EXCO' mcomp ccli0  mcomp 'NOID' ;
   rvm . 'MATASS' =mtot ;
   rvm . 'MAPREC' =mtot ;
*dbg   ct tt   = 'KRES' mtot ftot 'CLIM' cltot 'TYPI' rvm
*dbg      'IMPR' 3 'LTIME' VRAI ;
   ct    = 'KRES' mtot ftot 'CLIM' cltot 'TYPI' rvm ;
   dvitstar = '+' dvitstar ct ;
'FIN' iicomp ;   
*   HCHPO dvitstar 'dvitstar' ;
*
* Calcul de l'incrément de pression dphi pour vérifier la conservation
* de la masse
*
* Pas besoin de mettre dvitstar à 0 sur les clims Dirichlet car c'est
* déjà le cas !
smb1 = '*' mkcontr dvitstar ;
mtot  = mklapphi ;
ftot  = '-' smb1 ('EXCO' nipre residu nipre 'NOID') ;
cltot = 'EXCO' nipre ccli0 nipre 'NOID' ;
rvgm .'MATASS' =mtot ;
rvgm . 'MAPREC' =mtot ;
dphi  = 'KRES' mtot ftot 'CLIM' cltot 'TYPI' rvgm ;
'SI' ('NEG' alfa 1.D0 1.D-13) ;
   dphi = '/' dphi alfa ;
'FINSI' ;
*HCHPO dphi 'dphi' ;
*
* Calcul de l'incrément final de vitesse dvit
*
gdphi = '*' mkcontrt dphi ;
*  On surcharge par 0. là où il y a des clims de Dirichlet
*  Pas besoin !!! vu que imkvitd est là !
*!! gdphi = '-' gdphi ('*' gdphi ccli1 niccli niccli niccli) ;
gdphi = '*' gdphi imkvitd nivit nivit nivit ;
'SI' ('NEG' alfa 1.D0 1.D-13) ;
   gdphi = '*' gdphi alfa ;
'FINSI' ;
dvit  = dvitstar '-' gdphi ;
*HCHPO dvit 'dvit' ;
*
'SI' dblproj ;
   dphi = dphi '+' dphi1 ;
'FINSI' ;   
incrfin = dvit '+' dphi ;
solvp = incovp '+' incrfin ;
*
sol = '+' solvp sol ;
*
'RESPRO' sol ;
*
* End of procedure file GRESP
*
'FINPROC' ;
 

© Cast3M 2003 - Tous droits réservés.
Mentions légales