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