Télécharger la plaquette Cast3M


fr en

Télécharger invide.dgibi

Retour à la liste

Numérotation des lignes : 

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

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