* fichier : invide.dgibi ************************************************************************ * Section : Mathematiques Fonctions ************************************************************************ *********************************************************** **** 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' ;