* fichier : invide.dgibi
***********************************************************
**** OPERATEURS 'KOPS' et 'KRES' ****
**** Creation et inversion de la matrice identité ****
**** ****
**** A. BECCANTINI, TTMF AOUT 2000 ****
***********************************************************
'OPTION' 'ECHO' 0 ;
'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 = 5 ;
N1 = 1 '*' RAF ;
N2 = 1 '*' RAF ;
N3 = 1 '*' RAF ;
N4 = 1 '*' 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;
'SI' GRAPH ;
'TRACER' DOMTOT ;
'FINSI' ;
$DOMTOT = 'DOMA' DOMTOT ;
$DOM1 = 'DOMA' DOM1 'INCL' $DOMTOT 0.001 ;
$DOM2 = 'DOMA' DOM2 'INCL' $DOMTOT 0.001 ;
LMOT = 'MOTS' 'UX' 'UY' ;
MAT1 = 'KOPS' 'MATIDE' LMOT DOMTOT 'MATRIK' ;
* 'OPTION' DONN 5 ;
* OBJ1 = MANU 'OBJE' 'MAILLAGE' 75661 ;
* MATOLD = 'KOPS' 'MATIDE' LMOT DOMTOT ;
* On test 'KRES'
* On cree la table de sous-type 'METHINV'
* TAB1 = 'TABLE' 'METHINV' ;
TAB1 = ('EQEX') . 'METHINV' ;
* Methode d'inversion
TAB1 . 'TYPINV' = 1 ;
* Indice d'impression
TAB1 . 'IMPINV' = 1 ;
* Matrice pour assurer que la matrice à inverser est correctement
* assemblé
TAB1 . 'MATASS' = MAT1 ;
* Methode de numerotation de DDL
TAB1 . 'TYRENU' = 'SLOA' ;
* Prise en compte des multiplicateurs de Lagrange
*
TAB1 . 'PCMLAG' = 'APR2' ;
***********************************************************************
**** Indices spécifiques aux méthodes itératives (TYPINV=2..5) ********
***********************************************************************
* Champoint d'initialisation de la méthode
* TAB1 . 'XINIT' =
* Matrice de même structure que MA1 (éventuellement égale)
* servant au calcul du préconditionnement.
TAB1 . 'MAPREC' = MAT1 ;
* Nombre maximum d'itérations à effectuer.
TAB1 . 'NITMAX' = 1000 ;
* Norme maximum (L2 normé par le second membre) du
* résidu b-Ax calculé
TAB1 . 'RESID' = 1.D-10 ;
* Type de préconditionnement :
* 0 : pas de préconditionnement
* 1 : préconditionnement par la diagonale
* 2 : préconditionnement D-ILU (ca ne sers a rien!)
* défaut -> 3 : préconditionnement ILU(0) (Crout)
* 4 : préconditionnement MILU(0) (Crout modifié) (ca ne sers a rien!)
* 5 : préconditionnement ILUT (dual truncation)
* (plus cher que ILU, mais plus stable si on
* incombre bien la memeoire !)
* 6 : préconditionnement ILUT2 (une variante du
* précédent qui remplit mieux la mémoire et
* fonctionne mieux quelquefois)
*
TAB1 . 'PRECOND' = 3 ;
* ILUTLFIL (type REEL) :
* Pour un préconditionnement ILUT ou ILUT2 :
* encombrement maximal du préconditionneur par rapport à la
* matrice.
* Par défaut : 2.D0
*
TAB1 . 'ILUTLFIL' = 2.0 ;
* ILUTDTOL (type REEL) :
* Pour un préconditionnement ILUT ou ILUT2 :
* "drop tolerance" pour le préconditionneur, i.e. en-dessous de
* cette valeur relative, les termes de la factorisation
* incomplète seront oubliés.
* Par défaut : 0.D0 (on garde tous les termes).
TAB1 . 'ILUTDTOL' = 0.0D0 ;
* BCGSBTOL (type REEL) :
* 'Breakdown tolerance' pour les méthodes de type
* BiCGSTAB. Si un certain produit scalaire de vecteurs
* "direction" est inférieur à cette tolérance, la
* méthode s'arrete.
* Par défaut : 1.D-40
TAB1 . 'BCGSBTOL' = 1.D-40 ;
* MILURELX (type REEL) :
* Paramètre de relaxation pour le préconditionnement
* MILU(0) compris entre 0. et 1.
* S'il est égal à 0, on se ramène à ILU(0)
* S'il est égal à 1, MILU(0) est dit non relaxé
* Par défaut : 1.D0
*GMRESTRT (type ENTIER) :
*
* Paramètre m de redémarrage pour la méthode GMRES(m).
* Par défaut : 50
* CONVINV (type LISTREEL) :
* Une fois la résolution achevée, contient
* l'historique de la valeur de ||b-Ax||/||b|| en
* fonction du nombre d'itérations effectuées.
* (|| || est la norme L2)
CHPB1 = 'BRUI' 'BLAN' 'UNIF' 17.1 15.0 ($DOMTOT .'SOMMET') ;
CHPB2 = 'BRUI' 'BLAN' 'UNIF' 1.1 115.0 ($DOMTOT . 'SOMMET') ;
CHPB = ('NOMC' 'UX' CHPB1 'NATURE' 'DISCRET') 'ET'
('NOMC' 'UY' CHPB2 'NATURE' 'DISCRET') ;
*
* CHPOINT VIDE
*
CHPVID = 'KOPS' 'MULT' MAT1 ('MANUEL' 'CHPO' ($DOMTOT . SOMMET) 1
'SCAL' 1.0) ;
*
**** Resolution directe (Crout)
*
CHPR = 'KRES' MAT1 'TYPI' TAB1
'CLIM' CHPVID
'SMBR' CHPB
'IMPR' 1 ;
ERRO = 'MAXIMUM' (CHPR '-' CHPB) 'ABS' ;
'SI' (ERRO > 1.0D-14) ;
'MESSAGE' ;
'MESSAGE' 'Probleme dedans KRES. 1' ;
'MESSAGE' 'ERRO = ' ERRO ;
'MESSAGE' ;
'ERREUR' 5 ;
'FINSI' ;
*
**** CG
*
TAB1 . 'TYPINV' = 2 ;
TAB1 . 'XINIT' = 0.0 '*' CHPB ;
CHPR = 'KRES' MAT1 'TYPI' TAB1
'CLIM' CHPVID
'SMBR' CHPB
'IMPR' 1 ;
ERRO = 'MAXIMUM' (CHPR '-' CHPB) 'ABS' ;
'SI' (ERRO > 1.0D-14) ;
'MESSAGE' ;
'MESSAGE' 'Probleme dedans KRES. 2' ;
'MESSAGE' 'ERRO = ' ERRO ;
'MESSAGE' ;
'ERREUR' 5 ;
'FINSI' ;
*
**** BiCGSTAB
*
TAB1 . 'TYPINV' = 3;
TAB1 . 'XINIT' = 0.0 '*' CHPB ;
CHPR = 'KRES' MAT1 'TYPI' TAB1
'CLIM' CHPVID
'SMBR' CHPB
'IMPR' 1 ;
ERRO = 'MAXIMUM' (CHPR '-' CHPB) 'ABS' ;
'SI' (ERRO > 1.0D-14) ;
'MESSAGE' ;
'MESSAGE' 'Probleme dedans KRES. 3' ;
'MESSAGE' 'ERRO = ' ERRO ;
'MESSAGE' ;
'ERREUR' 5 ;
'FINSI' ;
'SI' FAUX ;
*
**** BiCGSTAB(2)
*
* Il ne faut pas les preconditioner!
TAB1 . 'PRECOND' = 0 ;
TAB1 . 'TYPINV' = 4;
TAB1 . 'XINIT' = 0.0 '*' CHPB ;
CHPR = 'KRES' MAT1 'TYPI' TAB1
'CLIM' CHPVID
'SMBR' CHPB
'IMPR' 1 ;
ERRO = 'MAXIMUM' (CHPR '-' CHPB) 'ABS' ;
'SI' (ERRO > 1.0D-14) ;
'MESSAGE' ;
'MESSAGE' 'Probleme dedans KRES. 4' ;
'MESSAGE' 'ERRO = ' ERRO ;
'MESSAGE' ;
'ERREUR' 5 ;
'MESSAGE' ;
'FINSI' ;
*
**** Il ne marche pas
*
'FINSI' ;
*
**** GMRES(m)
*
TAB1 . 'TYPINV' = 5 ;
TAB1 . 'XINIT' = 0.0 '*' CHPB ;
TAB1 . 'GMRESTRT' = 50 ;
* 50 = 50 directions dans l'espace de KRILOW
CHPR = 'KRES' MAT1 'TYPI' TAB1
'CLIM' CHPVID
'SMBR' CHPB
'IMPR' 1 ;
ERRO = 'MAXIMUM' (CHPR '-' CHPB) 'ABS' ;
'SI' (ERRO > 1.0D-12) ;
'MESSAGE' ;
'MESSAGE' 'Probleme dedans KRES. 5' ;
'MESSAGE' 'ERRO = ' ERRO ;
'MESSAGE' ;
'ERREUR' 5 ;
'FINSI' ;
'FIN' ;
© Cast3M 2003 - All rights reserved.
Disclaimer