Accueil | Cast3M

* 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' ;