Télécharger la plaquette Cast3M


fr en

Télécharger equ_chaleur2D.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  equ_chaleur2D.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. ************************************************************************
  5. * NOM         : equ_chaleur2D.dgibi
  6. * ___
  7. *
  8. * DESCRIPTION : Solution stationnaire de l'équation de la chaleur  (2D)
  9. * ___________  
  10. *
  11. *               GEOMETRIE : Un carré ! (il sera translaté et tourné)
  12. *               ----------
  13. *
  14. *               y
  15. *                 ^   y=1
  16. *                 |------------
  17. *                 |           |
  18. *                 |           |
  19. *                 |           |
  20. *                 |x = 0      |x=1 
  21. *                 |           |    
  22. *                 |           |    
  23. *                O-----------------------------> x
  24. *                     y=0                                   
  25. *
  26. *
  27. *               EQUATIONS :
  28. *               ----------
  29. *
  30. *               - Equations :
  31. *
  32. *                  mu Laplacien(T) = 0 avec mu = 1 ici
  33. *                  
  34. *               - Conditions aux limites :
  35. *
  36. *                  conditions de Dirichlet sur tout le bord :
  37. *                  on prend la restriction de la solution exacte au bord
  38. *
  39. *               - Solution exacte :
  40. *
  41. *                  T(x,y)= alpha x + beta y + gamma xy + delta
  42. *
  43. *                  pour une discrétisation avec de éléments linéaires,
  44. *                  on prend gamma = 0 (car les triangles linéaires
  45. *                  disponibles ne discrétisent pas de façon exacte le
  46. *                  terme bilinéaire en xy).
  47. *
  48. * DISCRETISATION : On teste tous les types d'éléments continus disponibles
  49. * ______________   aujourd'hui (99/03/04) (voir la notice NAVI) :
  50. *                  * 'QUAF' (température au sommet) ;
  51. *                  * 'MACR' (température au sommet) ;
  52. *                  * 'LINE' (température au sommet).
  53. *                  On devrait également tester les discrétisations
  54. *                  discontinus (linéaires ou Ctes par élément). Ceci
  55. *                  n'est pas fait.
  56. *
  57. *       Le maillage est construit avec l'opérateur SURF, il est perturbé
  58. *       par un bruit blanc gaussien. Il est ensuite translaté d'un
  59. *       vecteur arbitraire et tourné d'un angle arbitraire.
  60. *
  61. * Opérateurs utilisés : LAPN (option EF implicite)
  62. *
  63. *
  64. * RESOLUTION :     - Tous les types de solveurs (voir la notice KRES):
  65. * ___________        * Directs (Crout)
  66. *                    * Itératifs (Gradient Conjugué, BiCGSTAB, GMRES(50))
  67. *                  - Pour les solveurs itératifs, tous les types de
  68. *                    préconditionneurs :
  69. *                    * sans ;
  70. *                    * Jacobi (diagonal) ;
  71. *                    * D-ILU (pour mat. symétrique) ;
  72. *                    * ILU(0) (Choleski incomplet) ;
  73. *                    * MILU(0) (Choleski incomplet relaxé).
  74. *                    5 : préconditionnement ILUT (dual truncation)
  75. *                    6 : préconditionnement ILUT2 (une variante du
  76. *                       précédent qui remplit mieux la mémoire et
  77. *                       fonctionne mieux quelquefois)
  78. *
  79. *
  80. * TESTS EFFECTUES : La solution du problème discrétisé doit être la 
  81. * _______________   même que la solution exacte aux erreurs de
  82. *                   troncature près et à la précision sur le résidu prés
  83. *                   dans le cas des solveurs itératifs.
  84. *                   On impose donc des valeurs faibles pour l'écart en
  85. *                   norme L2 entre solution exacte et calculée :
  86. *                   * 1.D-14 pour le solveur direct ;
  87. *                   * 1.D-12 pour les solveurs itératifs.
  88. *
  89. * LANGAGE     : GIBIANE-CASTEM 2000
  90. * AUTEUR      : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF)
  91. *               mél : gounand@semt2.smts.cea.fr
  92. ************************************************************************
  93. * VERSION    : v1, 09/02/99, version initiale
  94. * HISTORIQUE : v1, 09/02/99, création
  95. * HISTORIQUE : 22/03/00, gounand
  96. *   Ajout du test des préconditionneurs ILUT
  97. * HISTORIQUE :
  98. ************************************************************************
  99. * Prière de PRENDRE LE TEMPS de compléter les commentaires
  100. * en cas de modification de ce sous-programme afin de faciliter
  101. * la maintenance !
  102. ************************************************************************
  103. interact= FAUX ;
  104. complet = FAUX ;
  105. graph   = FAUX ;
  106. *
  107. 'OPTION' 'DIME' 2 'ELEM' QUA4 'MODE' 'PLAN' ;
  108. 'OPTION' 'ISOV' 'SULI' ;
  109. 'OPTION' 'ECHO' 0 ;
  110. nbisov = 15 ;
  111. 'SI' ('NON' interact) ;
  112.   'OPTION' 'TRAC' 'PS' ;
  113.   'OPTION' 'ISOV' 'LIGNE' ;
  114. 'FINSI' ;  
  115. *
  116. ** Erreur Linfini entre deux Champoints.
  117. *
  118. DEBPROC CALCERR vitp1*'CHPOINT' vit*'CHPOINT' ;
  119.   err = MAXI (vitp1 - vit) 'ABS' ;
  120. FINPROC err ;
  121. *
  122. ** Erreur Pseudo L2 entre deux Champoints.
  123. *
  124. DEBPROC CALCERR2 vitp1*'CHPOINT' vit*'CHPOINT' ;
  125.   er0 = 'KOPS' vitp1 '-' vit ;
  126.   er2 = 'KOPS' er0 '*' er0 ;
  127.   err = 0.D0 ;
  128.   suppv  = 'EXTRAIRE' er2 'MAIL' ;
  129.   compv  = 'EXTRAIRE' er2 'COMP' ;
  130.   nbcomp = 'DIME' compv ;
  131.   nptd   = 'NBNO' suppv ;
  132.   'REPETER' numcomp nbcomp ;
  133.   lacomp = 'EXTRAIRE' compv &numcomp ;
  134.     'REPETER' numpo nptd ;
  135.       lepoi = suppv 'POIN' &numpo ;
  136.       err = err '+' ('EXTRAIRE' er2 lacomp lepoi) ;
  137.     'FIN' numpo ;
  138.   'FIN' numcomp ;
  139.   err = err '/' nptd ;
  140.   err = err ** 0.5D0 ;  
  141. FINPROC err ;
  142. *
  143. * Procédure paramétrée (solveur, type de discrétisation)
  144. * renvoyant l'erreur en norme L2 sur la température.
  145. * On calcule une solution exacte de l'équation de Laplace
  146. * (équations de la chaleur) ;
  147. *
  148. 'DEBPROC' CALCUL ktypi*'ENTIER' kprec*'ENTIER' typdis*'MOT' ;
  149. ttypi = 'TABLE' ;
  150. ttypi . 1 = 'CHAINE' 'Crout' ;
  151. ttypi . 2 = 'CHAINE' 'Gradient Conjugué' ;
  152. ttypi . 3 = 'CHAINE' 'BiCGSTAB' ;
  153. ttypi . 5 = 'CHAINE' 'GMRES(50)' ;
  154. tprec = 'TABLE' ;
  155. tprec . 0 = 'sans' ;
  156. tprec . 1 = 'diag' ;
  157. tprec . 2 = 'D-ILU' ;
  158. tprec . 3 = 'ILU(0)' ;
  159. tprec . 4 = 'MILU(0)' ;
  160. tprec . 5 = 'ILUT(2)' ;
  161. tprec . 6 = 'ILUT2(2)' ;
  162. *
  163. * titre global pour les dessins
  164. *
  165. titglob = 'CHAINE' ' typdis=' typdis ' typinv=' (ttypi . ktypi)
  166. ' prec=' (tprec . kprec) ;
  167. *
  168. * Paramètres physiques
  169. *
  170. mu = 1. ;
  171. * Conditions aux limites
  172. alpha = (** 2. 0.5) ;
  173. beta  = (** 3. 0.5) ;
  174. 'SI' ('EGA' typdis 'QUAF') ;
  175.   gamma = (** 5. 0.5)  ;
  176. 'SINON' ;
  177.   gamma = 0. ;
  178. 'FINSI' ;    
  179. delta = (** 1.5 0.5) ;
  180. *
  181. *  Géométrie
  182. *
  183. pA = 0.   0. ;
  184. pB = 1.   0. ;
  185. pC = 1.   1. ;
  186. pD = 0.   1. ;
  187. *
  188. * Paramètres de la discrétisation de base
  189. *
  190. 'SI' complet ;
  191.   nAB = 10 ;
  192.   nBC = 13 ;
  193.   nCD = 15 ;
  194.   nDA = 17 ;
  195. 'SINON' ;
  196.   nAB = 4 ;
  197.   nBC = 5 ;
  198.   nCD = 6 ;
  199.   nDA = 7 ;
  200. 'FINSI' ;  
  201. *
  202. * Rotation et translation aditionnelle + bruit blanc
  203. * + raffinement
  204. *
  205. vtran = ((** 2 0.5) (* -1 (** 3 0.5))) ;
  206. orig  = (0.D0 0.D0) ;
  207. arot  = (** 360. 0.5) ;
  208. *arot  = 0.D0 ;
  209. *arot = 90.0D0 ;
  210. lesdens = 'PROG' (1. '/' nAB) (1. '/' nCD) (1. '/' nBC) (1. '/' nDA) ;
  211. ampbruit = (0.1 * ('MINIMUM' lesdens)) ;
  212. *
  213. * Géométrie discrétisée
  214. *
  215. bas    = 'DROIT' nAB pA pB ;
  216. droite = 'DROIT' nBC pB pC ;
  217. haut   = 'DROIT' nCD pC pD ;
  218. gauche = 'DROIT' nDA pD pA ;
  219. pourtour = bas 'ET' droite 'ET' haut 'ET' gauche ;
  220. mt = 'SURFACE' pourtour 'PLAN' ;
  221. 'TASSER' mt ;
  222. mt = 'ORIENTER' mt ;
  223. *
  224. * On rajoute le bruit
  225. *
  226. pmt = 'CHANGER' mt 'POI1' ;
  227. pcnt= 'CHANGER' pourtour 'POI1' ;
  228. inmt= 'DIFF' pmt pcnt;
  229. brbl1 = 'BRUIT' 'BLAN' 'GAUS' 0.D0 ampbruit inmt ;
  230. brbl2 = 'BRUIT' 'BLAN' 'GAUS' 0.D0 ampbruit inmt ;
  231. brbl = 'ET' ('NOMC' 'UX' brbl1) ('NOMC' 'UY' brbl2) ;
  232. mt = 'PLUS' mt brbl;
  233. *
  234. * Eventuellement, on trace le résultat
  235. *
  236. 'SI' graph ;
  237.    titgeo = 'CHAINE' 'mt ' 'NBPO=' ('NBNO' mt)
  238.                       ' NBELEM=' ('NBEL' mt) titglob ;
  239.    'TRACER' mt 'NOEUD' 'TITRE' titgeo ;
  240. 'FINSI' ;
  241. *
  242. * Modèle
  243. *
  244. _mt     = 'CHANGER' mt 'QUAF' ;
  245. _bas    = 'CHANGER' bas    'QUAF' ;
  246. _droite = 'CHANGER' droite 'QUAF' ;
  247. _haut   = 'CHANGER' haut   'QUAF' ;
  248. _gauche = 'CHANGER' gauche 'QUAF' ;
  249. 'ELIMINATION' 1.D-6 (_gauche 'ET' _haut 'ET' _droite 'ET' _bas
  250.                      'ET' _mt) ;
  251. $mt = 'MODELISER' _mt 'NAVIER_STOKES' typdis ;
  252. *
  253. * Solution exacte bilinéaire
  254. *
  255. xx yy = 'COORDONNEE' ('DOMA' $mt 'SOMMET') ;
  256. solex = 'KOPS' ('KOPS' xx '*' alpha) '+' ('KOPS' yy '*' beta) ;
  257. solex = 'KOPS' solex '+' ('KOPS' ('KOPS' xx '*' yy) '*' gamma) ;
  258. solex = 'KOPS' solex '+' delta ;
  259. solex = 'KCHT' $mt 'SCAL' 'SOMMET' solex ;
  260. *
  261. * On tourne et on translate
  262. *
  263. _tmt solex = 'TOURNER' _mt solex arot orig ;
  264. _tmt solex = 'PLUS'    _tmt solex vtran ;
  265. $tmt = 'MODELISER' _tmt 'NAVIER_STOKES' typdis ;
  266. 'ELIMINATION' 1.D-6 ('ET' ('DOMA' $tmt 'SOMMET')
  267.                           ('EXTRAIRE' solex 'MAIL')) ;
  268. solex = 'KCHT' $tmt 'SCAL' 'SOMMET' solex ;
  269. mtmt = 'DOMA' $tmt 'MAILLAGE' ;
  270. cmtmt = 'CONTOUR' mtmt ;
  271. *
  272. * Mise en place du calcul numérique
  273. *
  274. * équation de Laplace
  275. *
  276. rv = 'EQEX' $tmt 'ITMA' 1
  277. 'OPTI' 'EF' 'IMPL'
  278. 'ZONE' $tmt 'OPER' 'LAPN' alpha 'INCO' 'TN' 'TN' ;
  279. *
  280. * conditions aux limites
  281. *
  282. csolex = 'REDU' solex cmtmt ;
  283. rv = 'EQEX' rv
  284.  'CLIM' 'TN' 'TIMP' cmtmt csolex ;
  285. *                      ________________
  286. *
  287. *                       INITIALISATION
  288. *                      ________________
  289. *
  290. rv . 'INCO' = 'TABLE' 'INCO' ;
  291. rv . 'INCO' . 'TN' = 'KCHT' $tmt 'SCAL' 'SOMMET' 0. ; 
  292. *
  293. rv . 'NITER' = 1 ;
  294. rv . 'METHINV' . 'TYPINV'   = ktypi ;
  295. rv . 'METHINV' . 'PRECOND'  = kprec ;
  296. rv . 'METHINV' . 'ILUTLFIL' = 2.D0 ;
  297. rv . 'METHINV' . 'ILUTDTOL' = 0.D0 ;
  298. rv . 'METHINV' . 'IMPINV'   = 2 ;
  299. rv . 'METHINV' . 'NITMAX'   = 1000 ;
  300. rv . 'METHINV' . 'GMRESTRT' = 50 ;
  301. rv . 'METHINV' . 'RESID'    = 1.D-13 ;
  302. *                         __________
  303. *
  304. *                           CALCUL
  305. *                         __________
  306. *
  307. EXEC rv ;
  308. *
  309. * Résultats
  310. *
  311. 'SI' graph ;
  312. *
  313. * solutions exactes
  314. *
  315.   tn = solex ;
  316.   titt = 'CHAINE' 'Température exacte' titglob ;
  317.   'TRACER' tn mtmt cmtmt nbisov 'TITRE' titt ;
  318. *
  319. * graphe de convergence de la méthode itérative
  320. *
  321.   'SI' ('NEG' ktypi 1) ;
  322.     conver = (rv . 'METHINV' . 'CONVINV') ;
  323.     dimcon = 'DIME' conver ;
  324.     labs = 'PROG' 1.D0 'PAS' 1.D0 dimcon ;
  325.     lord = ('LOG' conver) '/' ('LOG' 10.D0) ;
  326.     evtot = 'EVOL' 'MANU' 'ITER' labs 'RESID' lord ;
  327.     titev = 'CHAINE' 'Historique de convergence' titglob ;
  328.     'DESSIN' evtot 'TITR' titev 'LEGE' ;
  329.   'FINSI' ;
  330.  
  331. *
  332. * solutions calculées
  333. *
  334.   tn = rv . 'INCO' . 'TN' ;
  335.   titt = 'CHAINE' 'Température calculée' titglob ;
  336.   'TRACER' tn mtmt cmtmt nbisov 'TITRE' titt ;
  337. 'FINSI' ;
  338. *
  339. * Calcul des erreurs par rapport à la solution analytique
  340. *
  341. tn = rv . 'INCO' . 'TN' ;
  342. errlit = CALCERR  tn solex ;
  343. errl2t = CALCERR2 tn solex ;
  344. 'MESSAGE' '-------------------------------------------------' ;
  345. 'MESSAGE' ('CHAINE' 'Erreur en norme LINF' titglob ' : ' errlit) ;
  346. 'MESSAGE' ('CHAINE' 'Erreur en norme L2  ' titglob ' : ' errl2t) ;
  347. 'MESSAGE' '-------------------------------------------------' ;
  348. 'FINPROC' errl2t ;
  349. * Fin de la procédure de calcul                             *
  350. *___________________________________________________________*
  351. *
  352. *-----------------------------------------------------------
  353. * Paramètres du calcul
  354. *
  355. * Types d'inversion dans ltypi :
  356. *   1 = Crout (direct) ; 2 = CG ; 3 = BiCGSTAB ; 5 = GMRES
  357. * Types de préconditionnement  dans lprec :
  358. *   0 = sans ; 1 = diag ; 2 = D-ILU ; 3 = ILU(0) ; 4 = MILU(0)
  359. *   5 = ILUT(2) ; 6 = ILUT2(2)
  360. * ldiscr : type d'éléments à tester.
  361. 'SI' complet ;
  362.   ldiscr = 'MOTS' 'QUAF' ;
  363.   ltypi  = 'LECT' 1 2 3 5 ;
  364.   lprec  = 'LECT' 3 ;
  365. 'SINON' ;  
  366.   ldiscr = 'MOTS' 'QUAF' 'MACR' 'LINE' ;
  367.   ltypi  = 'LECT' 1 2 3 5 ;
  368.   lprec  = 'LECT' 0 'PAS' 1 6 ;
  369. 'FINSI' ;  
  370. *
  371. *-----------------------------------------------------------*
  372. * Boucles sur les différents paramètres du calcul           *
  373. ok = VRAI ;
  374. *
  375. *  Boucle sur les discrétisations
  376. *
  377. 'REPETER' idiscr ('DIME' ldiscr) ;
  378.   discr = 'EXTRAIRE' ldiscr &idiscr ;
  379. *
  380. *  Boucle sur les méthodes de résolution
  381. *
  382.   'REPETER' itypi ('DIME' ltypi) ;
  383.     typi = 'EXTRAIRE' ltypi &itypi ;  
  384.     'SI' ('EGA' typi 1) ;
  385.       errl2t = CALCUL typi 0 discr ;
  386.       ok = ('ET' ok (errl2t '<' 1.D-14)) ;
  387.     'SINON' ;
  388.       'REPETER' iprec ('DIME' lprec) ;
  389.         prec = 'EXTRAIRE' lprec &iprec ;
  390.         errl2t = CALCUL typi prec discr ;
  391.         ok = ('ET' ok (errl2t '<' 1.D-12)) ;
  392.       'FIN' iprec ;
  393.     'FINSI' ;    
  394.     'SI' ('ET' ('NON' complet) ('NON' ok)) ;
  395.       'MESSAGE' 'On a dépassé les marges d_erreur :' ;
  396.       'MESSAGE' ' * 1.D-14  pour une méthode directe' ;
  397.       'MESSAGE' ' * 1.D-12  pour une méthode itérative' ;
  398.       'ERREUR' 5 ;
  399.     'FINSI' ;
  400.   'FIN' itypi ;
  401. 'FIN' idiscr ;
  402.  
  403. 'SI' interact ;
  404.    'OPTION' 'DONN' 5 ;
  405. 'FINSI' ;
  406. 'MESSAGE' 'Tout s_est bien passé' ;
  407. 'FIN' ;
  408.  
  409.  
  410.  
  411.  
  412.  
  413.  
  414.  
  415.  
  416.  
  417.  

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