Télécharger la plaquette Cast3M


fr en

Télécharger invdiag.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  invdiag.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. ***********************************************************
  5. **** OPERATEURS 'KOPS' et 'KRES'                       ****
  6. **** Creation et inversion de la matrice identité      ****
  7. ****                                                   ****
  8. **** A. BECCANTINI, TTMF    AOUT    2000               ****
  9. ***********************************************************
  10.  
  11.  
  12.  'OPTION'  'ECHO' 1 ; 
  13.  'OPTION'  'DIME' 2 ;
  14.  'OPTION'  'ELEM' 'TRI3' ;
  15.  'OPTION'  'TRAC' 'X' ;
  16.  
  17.  GRAPH = FAUX ;
  18.  
  19. * GRAPH = VRAI ;
  20.  
  21. *
  22. *** MAILLAGE
  23. *
  24.  
  25.  
  26.  P0 = 0.0D0 0.0D0 ;
  27.  P1 = 3.0D0 0.0D0 ;
  28.  P2 = 3.0D0 3.0D0 ;
  29.  P3 = 0.0D0 3.0D0 ;
  30.  P4 = 6.0D0 0.0D0 ;
  31.  P5 = 6.0D0 3.0D0 ;
  32.  
  33.  RAF = 5 ;
  34.  
  35.  N1 = 1 '*' RAF ;
  36.  N2 = 1 '*' RAF ;
  37.  N3 = 1 '*' RAF ;
  38.  N4 = 1 '*' RAF ;
  39.  N5 = 1 '*' RAF ;
  40.  N6 = 2 '*' RAF ;
  41.  N7 = 2 '*' RAF ;
  42.  
  43.  
  44.  LINEXT1 = ((P0 'DROIT'  N1 P1) 'ET' (P1 'DROIT'  N2 P2) 'ET' 
  45.            (P2 'DROIT'  N3 P3) 'ET' (P3 'DROIT'  N4 P0)) ;
  46.  
  47.  LINEXT2 = ((P1 'DROIT'  N2 P2) 'ET' (P2 'DROIT'  N5 P5) 'ET' 
  48.            (P5 'DROIT'  N6 P4) 'ET' (P4 'DROIT'  N7 P1)) ;          
  49.  
  50.  'OPTION' 'ELEM' QUA4 ;
  51.  DOM1 = 'SURFACE'  LINEXT1 'PLAN' ;
  52.  
  53.  'OPTION' 'ELEM' TRI3 ;
  54.  DOM2 = 'SURFACE'  LINEXT2 'PLAN' ;
  55.  
  56.  DOMTOT = DOM1 'ET' DOM2;
  57.  'ELIMINATION' 1D-6 DOMTOT;
  58.  
  59.  'SI' GRAPH ;
  60.     'TRACER' DOMTOT ;
  61.  'FINSI' ;
  62.  
  63.  $DOMTOT = 'DOMA' DOMTOT ;
  64.  $DOM1 = 'DOMA' DOM1 'INCL' $DOMTOT 0.001 ;
  65.  $DOM2 = 'DOMA' DOM2 'INCL' $DOMTOT 0.001 ;
  66.  
  67.  
  68. * LMOT = 'MOTS' 'UX' 'UY' ;
  69. * MAT1 = 'KOPS' 'MATDIAGO' LMOT (kcht $DOMTOT SCAL SOMMET 1.D0) ;
  70.  MAT1 = 'KOPS' 'MATDIAGO' (kcht $DOMTOT 'VECT' SOMMET (1.D0 1.D0))
  71.    'MATRIK' ;
  72. * 'OPTION' DONN 5 ;
  73. * OBJ1 = MANU 'OBJE' 'MAILLAGE'  75661 ;
  74. * MATOLD =  'KOPS' 'MATIDE' LMOT DOMTOT ;
  75.  
  76. * On test 'KRES'
  77. * On cree la table de sous-type 'METHINV'
  78.  
  79. * TAB1 = 'TABLE' 'METHINV' ;
  80.  TAB1 = ('EQEX') . 'METHINV' ;
  81.  
  82. * Methode d'inversion
  83.  
  84.  TAB1 .   'TYPINV' = 1 ;
  85.  
  86. * Indice d'impression
  87.  
  88.  TAB1 .  'IMPINV' = 1 ; 
  89.  
  90. * Matrice pour assurer que la matrice à inverser est correctement
  91. * assemblé
  92.  
  93. *  TAB1  . 'MATASS' = MAT1 ;
  94.  
  95. * Methode de numerotation de DDL
  96.  
  97.   TAB1  . 'TYRENU' = 'SLOA' ;
  98.  
  99. * Prise en compte des multiplicateurs de Lagrange 
  100. *
  101.   TAB1 . 'PCMLAG' = 'APR2' ;
  102.  
  103. ***********************************************************************
  104. **** Indices spécifiques aux méthodes itératives (TYPINV=2..5) ********       
  105. ***********************************************************************
  106.  
  107. * Champoint d'initialisation de la méthode                     
  108. * TAB1 .  'XINIT' = 
  109.  
  110. * Matrice de même structure que MA1 (éventuellement égale)     
  111. *            servant au calcul du préconditionnement.
  112.  
  113. *  TAB1  . 'MAPREC' = MAT1 ;
  114.  
  115. * Nombre maximum d'itérations à effectuer.
  116.  
  117.   TAB1  . 'NITMAX' = 1000 ;
  118.  
  119. * Norme maximum (L2 normé par le second membre) du             
  120. * résidu b-Ax calculé                                          
  121.  
  122.   TAB1  . 'RESID' = 1.D-10 ;
  123.  
  124. * Type de préconditionnement :                                 
  125. *                      0 : pas de préconditionnement                      
  126. *                      1 : préconditionnement par la diagonale            
  127. *                      2 : préconditionnement D-ILU (ca ne sers a rien!)
  128. *            défaut -> 3 : préconditionnement ILU(0) (Crout)
  129. *                      4 : préconditionnement MILU(0) (Crout modifié) (ca ne sers a rien!)
  130. *                      5 : préconditionnement ILUT (dual truncation)
  131. *                          (plus cher que ILU, mais plus stable si on
  132. *                          incombre bien la memeoire !)
  133. *                      6 : préconditionnement ILUT2 (une variante du      
  134. *                          précédent qui remplit mieux la mémoire et      
  135. *                          fonctionne mieux quelquefois)                  
  136. *  
  137.  
  138.   TAB1  . 'PRECOND' = 3 ;
  139.  
  140.  
  141. * ILUTLFIL   (type REEL) :                                        
  142. *            Pour un préconditionnement ILUT ou ILUT2 :                   
  143. *            encombrement maximal du préconditionneur par rapport à la    
  144. *            matrice.                                                     
  145. *            Par défaut : 2.D0                                            
  146. *  
  147.  
  148.   TAB1  . 'ILUTLFIL' = 2.0 ;
  149.  
  150. * ILUTDTOL (type REEL) :                                        
  151. *            Pour un préconditionnement ILUT ou ILUT2 :                   
  152. *            "drop tolerance" pour le préconditionneur, i.e. en-dessous de
  153. *            cette valeur relative, les termes de la factorisation        
  154. *            incomplète seront oubliés.                                   
  155. *            Par défaut : 0.D0 (on garde tous les termes).                
  156.  
  157.   TAB1  . 'ILUTDTOL' = 0.0D0 ;
  158.  
  159. * BCGSBTOL (type REEL) :                                        
  160. *            'Breakdown tolerance' pour les méthodes de type              
  161. *            BiCGSTAB. Si un certain produit scalaire de vecteurs         
  162. *            "direction" est inférieur à cette tolérance, la              
  163. *            méthode s'arrete.                                            
  164. *            Par défaut : 1.D-40
  165.  
  166.  TAB1 . 'BCGSBTOL' =  1.D-40 ;
  167.  
  168. * MILURELX (type REEL) :                                        
  169. *           Paramètre de relaxation pour le préconditionnement           
  170. *            MILU(0) compris entre 0. et 1.                               
  171. *            S'il est égal à 0, on se ramène à ILU(0)                     
  172. *            S'il est égal à 1, MILU(0) est dit non relaxé                
  173. *            Par défaut : 1.D0                                            
  174.  
  175. *GMRESTRT (type ENTIER) :
  176. *
  177. *            Paramètre m de redémarrage pour la méthode GMRES(m).         
  178. *            Par défaut : 50                                              
  179.  
  180. * CONVINV (type LISTREEL) :                                     
  181. *            Une fois la résolution achevée, contient                     
  182. *            l'historique de la valeur de ||b-Ax||/||b|| en               
  183. *            fonction du nombre d'itérations effectuées.                  
  184. *            (|| || est la norme L2)                                      
  185.  
  186.  
  187.  CHPB1 =  'BRUI' 'BLAN' 'UNIF' 17.1 15.0 ($DOMTOT .'SOMMET') ;
  188.  CHPB2 =  'BRUI' 'BLAN' 'UNIF' 1.1 115.0 ($DOMTOT . 'SOMMET') ;
  189.  
  190.  CHPB = ('NOMC' 'UX' CHPB1 'NATURE' 'DISCRET') 'ET'
  191.         ('NOMC' 'UY' CHPB2 'NATURE' 'DISCRET') ;
  192.  
  193. *
  194. * CHPOINT VIDE
  195. *
  196.  
  197.  CHPVID = '*' MAT1 ('MANUEL' 'CHPO' ($DOMTOT . SOMMET) 1
  198.     'SCAL' 1.0) ;
  199.  
  200. *
  201. **** Resolution directe (Crout)
  202. *
  203.  
  204.  
  205.  CHPR = 'KRES' MAT1 'TYPI' TAB1
  206.                     'CLIM' CHPVID
  207.                     'SMBR' CHPB
  208.                     'IMPR' 1 ;
  209.  
  210.  
  211.  ERRO = 'MAXIMUM' (CHPR '-' CHPB) 'ABS' ;
  212.  'MESSAGE' 'ERRO = ' ERRO  ;
  213.  
  214.  'SI' (ERRO > 1.0D-14) ;
  215.      'MESSAGE' ;
  216.      'MESSAGE' 'Probleme dedans KRES. 1' ;
  217.      'MESSAGE' ;
  218.      'ERREUR' 5 ;
  219.  'FINSI' ;
  220.  
  221.  
  222.  
  223. *
  224. **** CG
  225. *
  226.  
  227.  TAB1 . 'TYPINV' = 2 ;
  228.  
  229.  TAB1 . 'XINIT' = 0.0 '*' CHPB ;
  230.  
  231.  CHPR = 'KRES' MAT1 'TYPI' TAB1
  232.                     'CLIM' CHPVID
  233.                     'SMBR' CHPB
  234.                     'IMPR' 1 ;
  235.  
  236.  ERRO = 'MAXIMUM' (CHPR '-' CHPB) 'ABS' ; 
  237.  'MESSAGE' 'ERRO = ' ERRO  ;
  238.  'SI' (ERRO > 1.0D-14) ;
  239.      'MESSAGE' ;
  240.      'MESSAGE' 'Probleme dedans KRES. 2' ;
  241.      'MESSAGE' ;
  242.      'ERREUR' 5 ;
  243.  'FINSI' ;
  244.  
  245. *
  246. **** BiCGSTAB
  247. *
  248.  
  249.  TAB1 . 'TYPINV' = 3;
  250.  
  251.  TAB1 . 'XINIT' = 0.0 '*' CHPB ;
  252.  
  253.  CHPR = 'KRES' MAT1 'TYPI' TAB1
  254.                     'CLIM' CHPVID
  255.                     'SMBR' CHPB
  256.                     'IMPR' 1 ;
  257.  
  258.  ERRO = 'MAXIMUM' (CHPR '-' CHPB) 'ABS' ;                     
  259.  'MESSAGE' 'ERRO = ' ERRO  ;
  260.  'SI' (ERRO > 1.0D-14) ;
  261.      'MESSAGE' ;
  262.      'MESSAGE' 'Probleme dedans KRES. 3' ;
  263.      'MESSAGE' ;
  264.      'ERREUR' 5 ;
  265.  'FINSI' ;
  266.  
  267.  
  268.  'SI' FAUX ;
  269.  
  270. *
  271. **** BiCGSTAB(2)
  272. *
  273. *    Il ne faut pas les preconditioner!
  274.  
  275.  
  276.     TAB1  . 'PRECOND' = 0 ;
  277.  
  278.     TAB1 . 'TYPINV' = 4;
  279.  
  280.     TAB1 . 'XINIT' = 0.0 '*' CHPB ;
  281.  
  282.     CHPR = 'KRES' MAT1 'TYPI' TAB1
  283.                        'CLIM' CHPVID
  284.                        'SMBR' CHPB
  285.                        'IMPR' 1 ;
  286.  
  287.     ERRO = 'MAXIMUM' (CHPR '-' CHPB) 'ABS' ;
  288.     'MESSAGE' 'ERRO = ' ERRO  ;
  289.  
  290.     'SI' (ERRO > 1.0D-14) ;
  291.         'MESSAGE' ;
  292.         'MESSAGE' 'Probleme dedans KRES. 4' ;
  293.         'MESSAGE' ;
  294.         'ERREUR' 5 ;
  295.         'MESSAGE' ;
  296.     'FINSI' ;
  297.  
  298. *
  299. **** Il ne marche pas
  300. *
  301.  
  302.  'FINSI' ;
  303.  
  304. *
  305. **** GMRES(m)
  306. *
  307.  
  308.  TAB1 . 'TYPINV' = 5 ;
  309.  
  310.  TAB1 . 'XINIT' = 0.0 '*' CHPB ;
  311.  
  312.  TAB1 . 'GMRESTRT' = 50 ;
  313.  
  314. * 50 = 50 directions dans l'espace de KRILOW 
  315.  
  316.  CHPR = 'KRES' MAT1 'TYPI' TAB1
  317.                     'CLIM' CHPVID
  318.                     'SMBR' CHPB
  319.                     'IMPR' 1 ;
  320.  
  321.  ERRO = 'MAXIMUM' (CHPR '-' CHPB) 'ABS' ;                     
  322.  'MESSAGE' 'ERRO = ' ERRO  ;
  323.  'SI' (ERRO > 1.0D-12) ;
  324.      'MESSAGE' ;
  325.      'MESSAGE' 'Probleme dedans KRES. 5' ;
  326.      'MESSAGE' ;
  327.      'ERREUR' 5 ;
  328.  'FINSI' ;
  329.  
  330.  'FIN' ;
  331.  
  332.  
  333.  
  334.  
  335.  
  336.  
  337.  
  338.  
  339.  
  340.  
  341.  
  342.  
  343.  
  344.  

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