Télécharger la plaquette Cast3M


fr en

Télécharger equ_chaleur2D_tenseur_VF2vfsym.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  equ_chaleur2D_tenseur_VF2.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. ***********************************************************************
  5. * NOM         : equ_chaleur2D_VF2_tenseur.dgibi
  6. * ___
  7. *
  8. * DESCRIPTION : Solution stationnaire de l'équation de la chaleur  (2D)
  9. * ___________
  10. *
  11. *               GEOMETRIE : Un carré
  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. *                   div K GRAD T = f   
  33. *
  34. *                           |         |
  35. *                  avec K = |x+2  x   |
  36. *                           |x    y+2 |
  37. *                           |         |
  38. *
  39. *                  f = (9x^2) + (12y^2) + (12x) + (12y)
  40. *
  41. *               - Conditions aux limites :
  42. *
  43. *                  conditions de Dirichlet sur tout le bord except en
  44. *                  (0;0);
  45. *                  on prend la restriction de la solution exacte au bord
  46. *
  47. *               - Solution exacte :
  48. *
  49. *                  T(x,y)= x**3    + y**3 
  50. *
  51. *
  52. * DISCRETISATION : une méthode de Volume Finis d'ordre 2 en espace.
  53. * ______________
  54. *
  55. *
  56. *
  57. *
  58. *       Le maillage est construit avec l'opérateur SURF, il est perturbé
  59. *       par un bruit blanc gaussien.
  60. *
  61. * Opérateurs utilisés : PENT (option VF implicite)
  62. *                       LAPN (option VF implicite)
  63. *
  64. * RESOLUTION :     - Solveur BiCGStab
  65. * __________       - Préconditionneur ILU(0)
  66. *
  67. * TESTS EFFECTUES : Vérification de l'ordre 2 en espace de la méthode
  68. * _______________   (utilisation d'une norme pseudo-L2) et de la
  69. *                   précision absolue sur le maillage le plus fin.
  70. *
  71. *
  72. *
  73. * LANGAGE     : GIBIANE-CASTEM 2000
  74. * AUTEUR      : C LE POTIER 
  75. *               mél : clepotier@cea.fr 
  76. ************************************************************************
  77. * VERSION    : v1, 12/02, version initiale
  78. * HISTORIQUE :
  79. ************************************************************************
  80. * Prière de PRENDRE LE TEMPS de compléter les commentaires
  81. * en cas de modification de ce sous-programme afin de faciliter
  82. * la maintenance !
  83. ******************************************************************
  84.  
  85. interact= FAUX ;
  86. complet = FAUX ;
  87. graph   = FAUX ;
  88. logxmgr = FAUX ;
  89. *
  90. 'OPTION' 'DIME' 2 'ELEM' TRI3 'MODE' 'PLAN' ;
  91. 'OPTION' 'ISOV' 'SULI' ;
  92. 'OPTION' 'ECHO' 0 ;
  93. nbisov = 15 ;
  94. 'SI' ('NON' interact) ;
  95.   'OPTION' 'TRAC' 'PS' ;
  96.   'OPTION' 'ISOV' 'LIGNE' ;
  97. 'FINSI' ;
  98. *
  99. ** Erreur Linfini entre deux Champoints.
  100. *
  101. DEBPROC CALCERR vitp1*'CHPOINT' vit*'CHPOINT' ;
  102.   err = MAXI (vitp1 - vit) 'ABS' ;
  103. FINPROC err ;
  104. *
  105. ** Erreur L2 entre deux Champoints.
  106. *
  107. DEBPROC CALCERR2 vitp1*'CHPOINT' vit*'CHPOINT' vol*'CHPOINT' ;
  108.   er2 = vitp1 '-' vit ;
  109.   compv  = 'EXTRAIRE' er2 'COMP' ;
  110.   er2 = 'PSCAL' er2 er2 compv compv ;
  111.   suppv  = 'EXTRAIRE' vol 'MAIL' ;
  112.   chpun = 'MANUEL' 'CHPO' suppv 1 'SCAL' 1 ;
  113.   voltot = 'XTY' chpun ('MOTS' 'SCAL') vol ('MOTS' 'SCAL') ;
  114.   error = 'XTY' er2 ('MOTS' 'SCAL') vol ('MOTS' 'SCAL') ;
  115.   error = error '/' voltot ;
  116.   error = error '**' 0.5 ;
  117. FINPROC error ;
  118. *
  119. * Procédure paramétrée (raffinement)
  120. * renvoyant l'erreur en norme L2 sur la température.
  121. * On calcule une solution de l'équation de Laplace
  122. * (équations de la chaleur) ;
  123. *
  124. 'DEBPROC' CALCUL nraff*'ENTIER' ;
  125. *nraff=2 ;
  126. *
  127. * titre global pour les dessins
  128. *
  129. titglob = 'CHAINE' ' ; nraff=' nraff ;
  130. *
  131. * Paramètres physiques
  132. *
  133. cv=  1.0 ;
  134. mu = 0.0 ;
  135. kappa = 1.0 ;
  136. rho = 1.0 ;
  137. *
  138. *  Géométrie
  139. *
  140. L = 1.;
  141. H = 1;
  142. pA = 0.   0. ;
  143. pB = L   0. ;
  144. pC = L   1. ;
  145. pD = 0.   1. ;
  146. *
  147. * Paramètres de la discrétisation de base
  148. *
  149. 'SI' complet ;
  150.   nAB = 8 ;
  151.   nBC = 8 ;
  152.   nCD = 8 ;
  153.   nDA = 8 ;
  154. 'SINON' ;
  155.   nAB = 4 ;
  156.   nBC = 4 ;
  157.   nCD = 4 ;
  158.   nDA = 4 ;
  159.  
  160. 'FINSI' ;  
  161. *
  162. * Rotation et translation aditionnelle + bruit blanc
  163. * + raffinement
  164. *
  165. vtran = ((** 2 0.5) (* -1 (** 3 0.5))) ;
  166. orig  = (0.D0 0.D0) ;
  167. arot  = 11.3D0 ;
  168. lesdens = 'PROG' (1. '/' nAB) (1. '/' nCD) (1. '/' nBC) (1. '/' nDA) ;
  169. ampbruit = (0.15 * ('MINIMUM' lesdens)) ;
  170. *
  171. * Géométrie discrétisée
  172. *
  173. bas    = 'DROIT' nAB pA pB ;
  174. droite = 'DROIT' nBC pB pC ;
  175. haut   = 'DROIT' nCD pC pD ;
  176. gauche = 'DROIT' nDA pD pA ;
  177. pourtour = bas 'ET' droite 'ET' haut 'ET' gauche ;
  178. mt = 'SURFACE' pourtour 'PLAN' ;
  179. *
  180. * On rajoute le bruit
  181. *
  182. pmt = 'CHANGER' mt 'POI1' ;
  183. pcnt= 'CHANGER' pourtour 'POI1' ;
  184. inmt= 'DIFF' pmt pcnt;
  185. brbl1 = 'BRUIT' 'BLAN' 'GAUS' 0.D0 ampbruit inmt ;
  186. brbl2 = 'BRUIT' 'BLAN' 'GAUS' 0.D0 ampbruit inmt ;
  187. brbl = 'ET' ('NOMC' 'UX' brbl1) ('NOMC' 'UY' brbl2) ;
  188. 'DEPLACER' mt 'PLUS' brbl;
  189. *
  190. * On raffine nraff fois
  191. *
  192. 'SI' (nraff > 0) ;
  193.   'REPETER' bcl nraff ;
  194.     mt  = 'CHANGER' mt 'QUADRATIQUE' ;
  195.     $mt = 'DOMA' mt 'MACRO' ;
  196.     mt  = ('DOMA '$mt 'MAILLAGE') ;
  197.     'MENAGE' ;
  198.   'FIN' bcl ;
  199. 'FINSI' ;
  200. *
  201. * Eventuellement, on trace le résultat
  202. *
  203. 'SI' graph ;
  204.    titgeo = 'CHAINE' 'mt ' 'NBPO=' ('NBNO' mt)
  205.                       ' NBELEM=' ('NBEL' mt) titglob ;
  206.    'TRACER' mt 'NOEUD' 'TITRE' titgeo ;
  207. 'FINSI' ;
  208. *
  209. * Modèles
  210. *
  211. Mmt   = 'CHAN' mt   'QUAF' ;
  212. 'ELIM' Mmt 1.E-5 ;
  213. MDNS  = 'NAVIER_STOKES' ;
  214. $mt   = 'MODE' Mmt MDNS 'LINE' ;
  215. mt    = 'DOMA' $mt 'MAILLAGE' ;
  216.  
  217. modmt = 'MODELISER' ('DOMA' $mt 'MAILLAGE') 'THERMIQUE' ;
  218. mtm = ('DOMA' $mt 'MAILLAGE') ;
  219. mtc = 'COULEUR' ('DOMA' $mt 'CENTRE') 'JAUN' ;
  220. mtf = 'COULEUR' ('DOMA' $mt 'FACE') 'BLEU' ;
  221. cmt = 'DOMA' $mt 'ENVELOPP' ;
  222. Mcmt   = 'CHAN' cmt   'QUAF' ;
  223. 'ELIM' (Mmt 'ET' Mcmt) 1.E-5 ;
  224. $cmt   = 'MODE' Mcmt MDNS 'QUAF' ;
  225.  
  226. *
  227. *
  228. *
  229. * Solution exacte aux centres
  230. *
  231. xxc yyc = 'COORDONNEE' ('DOMA' $mt 'CENTRE') ;
  232. solexc = 'KOPS' ('KOPS' xxc '**' 3) '+' ('KOPS' yyc '**' 3) ;
  233. solexc = 'KCHT' $mt 'SCAL' 'CENTRE' solexc ;
  234. *
  235. * Gradient exact (aux faces)
  236. *
  237. xxf yyf = 'COORDONNEE' ('DOMA' $mt 'FACE') ;
  238. gradtx = 'NOMC' 'UX' (3 '*' (xxf '**' 2)) 'NATU' 'DISCRET' ;
  239. gradty = 'NOMC' 'UY' (3 '*' (yyf '**' 2)) 'NATU' 'DISCRET' ;
  240. *
  241. *
  242. * Second membre de l'équation 
  243. *
  244. sour = 'KOPS' ('KOPS' xxc '**' 2) '*' 9. ;
  245. sour = 'KOPS' sour '+' ('KOPS' ('KOPS' yyc '**'2) '*' 12);
  246. sour = 'KOPS' sour '+'  ('KOPS' yyc '*' 12.) ;
  247. sour = 'KOPS' sour '+' ('KOPS' xxc '*' 12.)  ;
  248. sour = 'KCHT' $mt 'SCAL' 'CENTRE' sour ;
  249.  
  250.  
  251. xxlim yylim = 'COORDONNEE' ('DOMA' $cmt  'SOMMET') ;
  252. solim = 'KOPS' ('KOPS' xxlim '**' 3) '+' ('KOPS' yylim '**' 3) ;
  253. solim = 'KCHT' $cmt 'SCAL' 'SOMMET' solim ;
  254.  
  255. *
  256. * Calcul du tenseur
  257. *
  258. K11 = 'KOPS' xxc  '+' 2. ;
  259. K22 = 'KOPS' yyc  '+' 2. ;
  260. K21 =  xxc  ;
  261.  
  262. K11 = NOMC(K11) 'K11' 'NATU' 'DISCRET';
  263. K22 = NOMC(K22) 'K22' 'NATU' 'DISCRET';
  264. K21 = NOMC(K21) 'K21' 'NATU' 'DISCRET';
  265. PERM = (K11 'ET' K22) ;
  266. PERM = PERM 'ET' K21;
  267.  
  268. *
  269. * Mise en place du calcul numérique
  270. *
  271. * Mise en place du calcul numérique
  272. *
  273. * équation de Laplace
  274. *
  275. * on utilise une méthode de Newton pour résoudre :
  276. * - \Delta T = 0     (\Delta opérateur laplacien)
  277. * - avec T  donné sur gauche (CL de Dirichlet)
  278. *        grad(T) \cdot next donné ailleurs
  279. *
  280. * T_0 : estimation initiale de la solution
  281. * On a T_1 = T_0 - {\Delta'}^{-1} (\Delta T_0)
  282. *
  283. * L'opérateur 'LAPN' 'VF' nous donne la matrice \Delta' et le
  284. * résidu \Delta T_0.
  285. * On n'inverse bien évidemment pas \Delta' mais on résoud le système:
  286. * \Delta' IncT = \Delta T_0
  287. * => IncT = {\Delta'}^{-1} (\Delta T_0)
  288. *
  289. * La méthode de Newton doit converger en un pas (on vérifie que le
  290. * résidu (\Delta T_1) est nul après le premier pas.
  291. *
  292. *
  293. tmp    = 0.0 ;
  294. *deltat = 5.D-3 ;
  295. t0 = 'KCHT' $mt 'SCAL' 'CENTRE' 0.D0 ;
  296. t0 = solexc; 
  297.  
  298. gradt0 mchamt = 'PENT'  $mt  'FACE' 'VFSYM' t0 
  299.                     'DISPDIF' PERM 'TIMP' solim ;
  300.  
  301.  
  302.  
  303. listinco = 'MOTS' 'RN' 'RUXN' 'RUYN' 'RETN' ;
  304. jaco chpres dt = 'LAPN' 'VF' 'CLAUDEIS'  'IMPL' 
  305. $mt t0 gradt0 mchamt 'TIMP' solim ;
  306.  
  307.  
  308. mamat = 'KOPS' 'MULT' -1.000D0 jaco ;
  309. matot =  mamat  ;
  310.  
  311. AUX = ('EXCO' chpres 'RETN' 'RETN') - (NOMC(sour) 'RETN');
  312. rv = 'EQEX' ;
  313. rv . 'METHINV' . 'TYPINV'   = 3 ;
  314. rv . 'METHINV' . 'PRECOND'  = 3 ;
  315. rv . 'METHINV' . 'MATASS'   = matot ;
  316. rv . 'METHINV' . 'MAPREC'   = matot ;
  317. rv . 'METHINV' . 'IMPINV'   = 3 ;
  318. res = 'KRES' matot 'TYPI' (rv . 'METHINV')
  319.                       'SMBR'AUX 
  320.                       'IMPR' 0 ;
  321. t1 = t0 '+' ('EXCO' res 'RETN' 'SCAL') ;
  322. gradt1 mchamt = 'PENT'  $mt  'FACE' 'VFSYM' t1 
  323.                     'DISPDIF' PERM 'TIMP' solim ;
  324. *listinco = 'MOTS' 'RN' 'RUXN' 'RUYN' 'RETN' ;
  325. jacbid chpres1 dt = 'LAPN' 'VF' 'CLAUDEIS' 'IMPL' 
  326. $mt  t0 gradt1 mchamt 'TIMP' solim ;
  327. AUX = ('EXCO' chpres1 'RETN' 'RETN') - (NOMC(sour) 'RETN');
  328. res = 'KRES' matot 'TYPI' (rv . 'METHINV')
  329.                       'SMBR'AUX 
  330.                       'IMPR' 0 ;
  331. mres1 = 'MAXIMUM' res  'ABS' ;
  332. 'MESSAGE' ('CHAINE' 'Maxi. mres1 =' mres1) ;
  333. 'SI' ('>' mres1 1.e-5) ;
  334.     'MESSAGE' 'La méthode de Newton na pas converge en un pas.'
  335.     'ERREUR' 5 ;
  336. 'FINSI' ;
  337. * 
  338. *
  339. * Résultats
  340. *
  341. 'SI' graph ;
  342. *
  343. * solutions exactes
  344. *
  345.   tn = solexc ;
  346.   chm_tnex = 'KCHA' $mt 'CHAM' tn ;
  347.   titt = 'CHAINE' 'Température exacte' titglob ;
  348.   'TRACER' chm_tnex modmt nbisov 'TITRE' titt ;
  349. *
  350. * graphe de convergence de la méthode itérative
  351. *
  352.     conver = (rv . 'METHINV' . 'CONVINV') ;
  353.     dimcon = 'DIME' conver ;
  354.     labs = 'PROG' 1.D0 'PAS' 1.D0 dimcon ;
  355.     lord = ('LOG' conver) '/' ('LOG' 10.D0) ;
  356.     evtot = 'EVOL' 'MANU' 'ITER' labs 'RESID' lord ;
  357.     titev = 'CHAINE' 'Historique de convergence' titglob ;
  358.     'DESSIN' evtot 'TITR' titev 'LEGE' ;
  359. *
  360. * solutions calculées
  361. *
  362.   tn = t1 ;
  363.   chm_tn = 'KCHA' $mt 'CHAM' tn ;
  364.   titt = 'CHAINE' 'Température calculée' titglob ;
  365.   'TRACER' chm_tn modmt nbisov 'TITRE' titt ;
  366. *
  367. * erreur
  368. *
  369.   titt = 'CHAINE' 'Abs (Température calculée - Température exacte)'
  370.          titglob ;
  371.   'TRACER' ('ABS' (chm_tnex '-' chm_tn)) modmt nbisov 'TITRE' titt ;
  372. 'FINSI' ;
  373. *
  374. * Calcul des erreurs par rapport à la solution analytique
  375. *
  376. vol = 'DOMA' $mt 'VOLUME' ;
  377. echloc = (('MESURE' mt) '/' ('NBEL' mt)) ** 0.5 ;
  378. tn = t1 ;
  379. errlit = CALCERR  tn solexc ;
  380. errl2t = CALCERR2 tn solexc vol ;
  381. 'MESSAGE' '-------------------------------------------------' ;
  382. 'MESSAGE' ('CHAINE' 'Erreur en norme LINF : ' errlit) ;
  383. 'MESSAGE' ('CHAINE' 'Erreur en norme L2   : ' errl2t) ;
  384. 'MESSAGE' '-------------------------------------------------' ;
  385. *'OPTION' 'DONN' 5 ;
  386. 'FINPROC' echloc errl2t ;
  387. * Fin de la procédure de calcul                             *
  388. *___________________________________________________________*
  389. *-----------------------------------------------------------
  390. * Paramètres du calcul
  391. *
  392. * lraff : nb raffinement du maillage (à chaque fois, on découpe
  393. *         les éléments en quatre).
  394. 'SI' complet ;
  395.   lraff  = 'LECT' 0 PAS 1 3 ;
  396. 'SINON' ;
  397.   lraff  = 'LECT' 0 PAS 1 3 ;
  398. 'FINSI' ;
  399. *
  400. *-----------------------------------------------------------*
  401. * Boucles sur les différents paramètres du calcul           *
  402. ordok = VRAI ;
  403. precok = VRAI ;
  404. * ordre théorique en espace de la méthode
  405. ordtth = 2 ;
  406. * erreur L2 max pour la solution (raffinement 2, complet=FAUX)
  407. errtmax = 4.D-3 ;
  408. *
  409. *  Boucle sur les raffinements
  410. *
  411.     lh = 'PROG' ;
  412.     lerl2t = 'PROG' ;
  413.     'REPETER' iraff ('DIME' lraff)  ;
  414.       raff = 'EXTRAIRE' lraff &iraff ;
  415.       h errt = CALCUL raff ;
  416.       lh     = 'ET' lh     ('PROG' h ) ;
  417.       lerl2t = 'ET' lerl2t ('PROG' errt) ;
  418.     'FIN' iraff ;
  419.     lh = lh '/' ('EXTRAIRE' lh 1) ;
  420.     lh     = ('LOG' lh)     '/' ('LOG' 10.D0) ;
  421.     lerl2t = ('LOG' lerl2t) '/' ('LOG' 10.D0) ;
  422.     tl2 = 'EVOL' 'MANU' 'Long. carac.' lh
  423.                         'Err. tempér.' lerl2t ;
  424.  
  425. *
  426. * Recherche des ordres de convergence
  427. *
  428.     cpt tlipol = @POMI tl2 1 'IDEM' ;
  429.     ordt = cpt . 1 ;
  430.     'MESSAGE' ('CHAINE' 'ordre température=' ordt) ;
  431. *
  432. * Tracés graphiques
  433. *
  434.     'SI' graph ;
  435.       tableg = 'TABLE' ;
  436.       tableg . 'TITRE' = 'TABLE' ;
  437.       tableg . 1 = 'MARQ CROI NOLI' ;
  438.       tableg . 'TITRE' . 1 = 'Erreur L2 calculée' ;
  439.       tableg . 2 = 'TIRR' ;
  440.       tableg . 'TITRE' . 2 = 'Erreur L2 interpolée' ;
  441.       titdest= 'CHAINE' 'Err L2 Tempér., ordre=' ordt ;
  442.       'DESSIN' (tl2 'ET' tlipol) 'TITRE' titdest tableg ;
  443.     'FINSI' ;
  444. *
  445. * Tests des ordres de convergence
  446. * Tests des erreurs L2 sur le maillage le plus fin dans le
  447. * cas complet=faux
  448. *
  449.      ordok = 'ET' ordok ('EGA' ('ENTIER' ('+' ordt 0.5)) ordtth) ;
  450.      'SI' ('EGA' complet FAUX) ;
  451.        precok = 'ET' precok ('<' errt errtmax) ;
  452.      'FINSI' ;
  453.      'FINSI' ;
  454.      'SI' ('NON' ordok) ;
  455.        'MESSAGE' 'On n_obtient pas un ordre de convergence correct' ;
  456.        'ERREUR' 5 ;
  457.      'FINSI' ;
  458.      'SI' ('NON' precok) ;
  459.        'MESSAGE' 'On n_obtient pas une précision correcte' ;
  460.        'ERREUR' 5 ;
  461.      'FINSI' ;
  462. 'SI' interact ;
  463.    'OPTION' 'ECHO' 1 ;
  464.    'OPTION' 'DONN' 5 ;
  465.    'SI' logxmgr ;
  466. *     Sortie pour xmgr
  467.       'OPTION' ECHO 0 ;
  468.       'OPTION' 'IMPR' 'file.txt' ;
  469.       ndim = 'DIME' lh ;
  470.       'REPETER' BL1 ndim ;
  471.          'MESSAGE' ('EXTRAIRE' lh &BL1) '  '  ('EXTRAIRE' lerl2t &BL1) ;
  472.       'FIN' BL1 ;
  473.       lh1 = 'EXTRAIRE' tlipol 'ABSC' ;
  474.       lerr = 'EXTRAIRE' tlipol 'ORDO' ;
  475.       ndim = 'DIME' lh1 ;
  476.       'OPTION' 'IMPR' 'file2.txt' ;
  477.       'REPETER' BL1 ndim ;
  478.          'MESSAGE' ('EXTRAIRE' lh1 &BL1) '  '  ('EXTRAIRE' lerr &BL1) ;
  479.       'FIN' BL1 ;
  480.       'OPTION' 'IMPR' 'poubelle.txt' ;
  481.       'FIN' ;
  482.    'FINSI' ;
  483. 'FINSI' ;
  484.  
  485. 'MESSAGE' 'Tout s_est bien passé' ;
  486. 'FIN' ;
  487.  
  488.  
  489.  
  490.  
  491.  
  492.  
  493.  
  494.  
  495.  
  496.  

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