Télécharger la plaquette Cast3M


fr en

Télécharger equ_chaleur3D_VFSYM.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier : equ_chaleur3D_VFSYM.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. ************************************************************************
  5. * NOM         : equ_chaleur3D_VF.dgibi
  6. * ___
  7. *
  8. * DESCRIPTION : Solution stationnaire de l'équation de la chaleur  (3D)
  9. * ___________  
  10. *
  11. *               GEOMETRIE : Un cube de côté 1 maillé avec 'VOLU'.
  12. *
  13. *               EQUATIONS :
  14. *               ----------
  15. *
  16. *               EQUATIONS :
  17. *               ----------
  18. *
  19. *               - Equations :
  20. *
  21. *                   div K GRAD T  = F
  22. *
  23. *                           |          
  24. *                  avec K = |+1  0   0|
  25. *                           |0    1  0|
  26. *                           |0    0  1|
  27. *
  28. *
  29. *               - Conditions aux limites :
  30. *
  31. *                  conditions de Dirichlet
  32. *
  33. *               - Solution exacte :
  34. *
  35. *                  T(x,y)= sin pix sin piy sin pi z 
  36. *
  37. *
  38. *
  39. *
  40. * DISCRETISATION : une méthode de Volume Finis d'ordre 2 en espace, avec la 
  41. *                  methode VFSYM . 
  42. *        SCHEMA PROPOSE PAR Christophe Le Potier
  43. *        Références : {C. Le Potier}
  44. *        \emph{Schema volumes finis pour des operateurs de diffusion
  45. *        fortement anisotropes sur des maillages non structures},
  46. *        C. R. Acad. Sci. Ser. I \textbf{340}, 2005, pp. 921--926.
  47. *        interfaces
  48. *
  49. *
  50. *
  51. *       Le maillage est construit avec l'opérateur VOLU.
  52. *
  53. * Opérateurs utilisés :  PENT (option VF implicite)
  54. *                        LAPN (option VF implicite)
  55. *
  56. *
  57. * RESOLUTION :     - Solveur BiCGStab
  58. * __________       - Préconditionneur ILU(0)
  59. *
  60. *
  61. * TESTS EFFECTUES : Vérification de l'ordre 2 en espace de la méthode
  62. * _______________   (utilisation d'une norme pseudo-L2) et de la 
  63. *                   précision absolue sur le maillage le plus fin.
  64. *
  65. *
  66. *
  67. *
  68. *
  69. * LANGAGE     : GIBIANE-CASTEM 2000
  70. * AUTEUR      : Christophe Le Potier(CEA/DEN/DM2S/SFME/LSET)
  71. *               mél : clepotier@cea.fr 
  72. ************************************************************************
  73. * VERSION    : v1, 24/03/2010, version initiale
  74. * HISTORIQUE :
  75. ************************************************************************
  76. * Prière de PRENDRE LE TEMPS de compléter les commentaires
  77. * en cas de modification de ce sous-programme afin de faciliter
  78. * la maintenance !
  79. ************************************************************************
  80. interact=  FAUX ;
  81. complet =  FAUX ;
  82. graph   = FAUX ;
  83. *
  84. 'OPTION' 'DIME' 3 'ELEM' CUB8   ;
  85. 'OPTION' 'ISOV' 'SULI' ;
  86. 'OPTION' 'ECHO' 1 ;
  87. nbisov = 15 ;
  88. 'SI' ('NON' interact) ;
  89.   'OPTION' 'TRAC' 'PSC' ;
  90. *  'OPTION' 'ISOV' 'LIGNE' ;
  91. 'FINSI' ;  
  92. *******************************************************************
  93. ** Erreur Linfini entre deux Champoints.
  94. *
  95. *******************************************************************
  96. DEBPROC CALCERR vitp1*'CHPOINT' vit*'CHPOINT' ;
  97.   AUX = (vitp1 - vit) ;
  98.   err = MAXI (AUX) 'ABS' ;
  99.   err = err / (MAXI(ABS(vit)));
  100. FINPROC err ;
  101. *
  102. *
  103. *******************************************************************       
  104. ** Erreur Pseudo L2 entre deux Champoints.
  105. *
  106. *******************************************************************       
  107. DEBPROC CALCERR2 vitp1*'CHPOINT' vit*'CHPOINT' ;
  108.   er0 = 'KOPS' vitp1 '-' vit ;
  109.   er2 = 'KOPS' er0 '*' er0 ;
  110.   er3 = 'KOPS' vit '*' vit ;
  111.   err = 0.D0 ;
  112.   erl = 0.D0;
  113.   suppv  = 'EXTRAIRE' er2 'MAIL' ;
  114.   compv  = 'EXTRAIRE' er2 'COMP' ;
  115.   nbcomp = 'DIME' compv ;
  116.   nptd   = 'NBNO' suppv ;
  117.   'REPETER' numcomp nbcomp ;
  118.   lacomp = 'EXTRAIRE' compv &numcomp ;
  119.     'REPETER' numpo nptd ;
  120.       lepoi = suppv 'POIN' &numpo ;
  121.       err = err '+' ('EXTRAIRE' er2 lacomp lepoi) ;
  122.       erl = erl '+' ('EXTRAIRE' er3 lacomp lepoi) ;
  123.     'FIN' numpo ;
  124.   'FIN' numcomp ;
  125.   err = err '/' erl ;
  126.   err = err ** 0.5D0 ;  
  127. FINPROC err ;
  128. *
  129. *******************************************************************       
  130. * Procédure paramétrée (raffinement)
  131. * renvoyant l'erreur en norme L2 sur la température.
  132. * On calcule une solution exacte de l'équation de Laplace
  133. * (équation de la chaleur) ;
  134. *
  135. 'DEBPROC' CALCUL nraff*'ENTIER' ;
  136. *
  137. * titre global pour les dessins
  138. *
  139. titglob = 'CHAINE' ' ; nraff=' nraff ;
  140. *
  141. *  Géométrie
  142. *
  143.  
  144. pA = 0.   0.  0. ;
  145. pB = 1.0   0.  0. ;
  146. pC = 1.0   1.0  0. ;
  147. pD = 0.    1.0  0. ;
  148. pE = 0.   0.  1. ;
  149. *
  150. * Paramètres de la discrétisation de base
  151. *
  152. 'SI' complet ;
  153.   nAB = 1 ;
  154.   nBC = 1 ;
  155.   nCD = 1 ;
  156.   nDA = 1 ;
  157.   nH  = 1 ;
  158. 'SINON' ;
  159.   nAB = 2 '*' (nraff '+' 1) ;
  160.   nBC = 2 '*' (nraff '+' 1) ;
  161.   nCD = 2 '*' (nraff '+' 1) ;
  162.   nDA = 2 '*' (nraff '+' 1) ;
  163.   nH  = 2 '*' (nraff '+' 1) ;
  164. 'FINSI' ;  
  165. *
  166. *******************************************************************       
  167. * Géométrie discrétisée
  168. *******************************************************************       
  169. *
  170. *'OPTION' 'ELEM' 'TRI3'  ;
  171. bas    = 'DROIT' nAB pA pB ;
  172. droite = 'DROIT' nBC pB pC ;
  173. haut   = 'DROIT' nCD pC pD ;
  174. gauche = 'DROIT' nDA pD pA ;
  175. spourt = bas 'ET' droite 'ET' haut 'ET' gauche ;
  176. sbas   = 'SURFACE' spourt 'PLAN' ;
  177.  
  178.  
  179. *'OPTION' 'DIME' 3 'ELEM' TET4   ;
  180. mt = 'VOLUME' sbas 'TRAN' nH  pE ;
  181. elim mt 0.001;
  182. f1 f2 f3  = 'FACE' mt  ;
  183. f1 = 'CHANGER' f1 'TRI3' ;
  184. f2 = 'CHANGER' f2 'TRI3' ;
  185. f3 = 'CHANGER' f3 'TRI3' ;
  186. ss = (f1 'ET' f2 'ET' f3) ;
  187. *'OPTION' 'ELEM' 'TET4' ;
  188. mt = 'VOLUME' ss ;
  189.  
  190. *'OPTION' 'ELEM' 'TET4' ;
  191.  
  192. * VERRUE
  193. pourtour = 'ENVELOPPE' mt ;
  194. *
  195. * Eventuellement, on trace le résultat
  196. *
  197. *   list mt;
  198. 'SI' graph ;
  199.    titgeo = 'CHAINE' 'mt ' 'NBPO=' ('NBNO' mt)
  200.                       ' NBELEM=' ('NBEL' mt) titglob ;
  201.    'TRACER' mt 'NOEUD' 'TITRE' titgeo ;
  202. 'FINSI' ;
  203. *******************************************************************       
  204. *
  205. * Modèle
  206. *
  207. *******************************************************************       
  208. L = 10.D0;
  209. mt = DEPLA (mt) AFFI (L) (0. 0. 0.) (0. 0. 1.);
  210. Mmt   = 'CHAN' mt   'QUAF' ;
  211. 'ELIM' Mmt 1.E-5 ;
  212. MDNS  = 'DARCY' ;
  213. $mt   = 'MODE' Mmt MDNS 'ANISOTROPE' ;
  214. mt    = 'DOMA' $mt 'MAILLAGE' ;
  215.  
  216. MESS 'NOMBRE DE MAILLES' (NBEL (MT)); 
  217.  
  218. modmt = 'MODELISER' ('DOMA' $mt 'MAILLAGE') 'THERMIQUE' ; 
  219. *   list ('DOMA' $mt 'FACEP') ;
  220. mtm = ('DOMA' $mt 'MAILLAGE') ;
  221. mtc = 'COULEUR' ('DOMA' $mt 'CENTRE') 'JAUN' ;
  222. mtf = 'COULEUR' ('DOMA' $mt 'FACE') 'BLEU' ;
  223. cmt = pourtour;
  224. Mcmt   = 'CHAN' cmt   'QUAF' ;
  225. 'ELIM' (Mmt 'ET' Mcmt) 1.E-5 ;
  226. $cmt   = 'MODE' Mcmt MDNS 'ANISOTROPE' ;
  227.  
  228. *******************************************************************       
  229. * Solution exacte aux centres
  230. *******************************************************************       
  231. re = 180.;
  232. ep = 1.D0;
  233. xxc yyc zzc = 'COORDONNEE' ('DOMA' $mt 'CENTRE') ;
  234. solexc = (sin(re*xxc))*(sin(re*yyc))*(sin(re*zzc/L));
  235. solexc = 'KCHT' $mt 'SCAL' 'CENTRE' solexc ;
  236.  
  237.  
  238. *******************************************************************       
  239. * terme source 
  240. *******************************************************************       
  241. sour = 0.0D0 -(((2.D0 + (ep*(L**-2))) *(pi*pi))*solexc);
  242. sour = 'KCHT' $mt 'SCAL' 'CENTRE' sour ;
  243.  
  244.  
  245.  
  246.  
  247. *******************************************************************
  248. * CONDITION DE DIRICHLET
  249. *******************************************************************
  250. xxlim yylim zzlim = 'COORDONNEE' ('DOMA' $cmt  'CENTRE') ;
  251. solim = (sin(re*xxlim))*(sin(re*yylim))*(sin(re*zzlim/L));
  252. solim = 'KCHT' $cmt 'SCAL' 'CENTRE' solim ;
  253.  
  254. *******************************************************************
  255. *
  256. * Calcul du tenseur
  257. *
  258. *******************************************************************
  259. K11 = 'KOPS' xxc  '+' 1.D0 ;
  260. K22 = 'KOPS' yyc  '+' 1.D0 ;
  261. K33 = 'KOPS' zzc  '+' 1.D0 ;
  262. K11 = 1.D0 + (0.0*xxc);
  263. K22 = 1.D0 + (0.0*yyc);
  264. K33 = 1.D0 + (0.0*zzc);
  265. K21 =  (0.0*zzc);
  266. K31 =  (0.0*zzc);
  267. K32 =  (0.0*zzc);
  268.  
  269. K11 = NOMC(K11) 'K11' 'NATU' 'DISCRET';
  270. K22 = NOMC(K22) 'K22' 'NATU' 'DISCRET';
  271. K33 = NOMC(K33) 'K33' 'NATU' 'DISCRET';
  272. K21 = NOMC(K21) 'K21' 'NATU' 'DISCRET';
  273. K31 = NOMC(K31) 'K31' 'NATU' 'DISCRET';
  274. K32 = NOMC(K32) 'K32' 'NATU' 'DISCRET';
  275.  
  276. PERM = (K11 'ET' K22) ;
  277. PERM = PERM 'ET' K33;
  278. PERM = PERM 'ET' K21;
  279. PERM = PERM 'ET' K31;
  280. PERM = PERM 'ET' K32;
  281.  
  282. * Mise en place du calcul numérique
  283. *
  284. * équation de Laplace
  285. *
  286. *
  287. * on utilise une méthode de Newton pour résoudre :
  288. * - \Delta T = 0     (\Delta opérateur laplacien)
  289. * - avec T_{\partial \Omega} donné (CL de Dirichlet)
  290. * T_0 : estimation initiale de la solution
  291. * On a T_1 = T_0 - {\Delta'}^{-1} (\Delta T_0)
  292. *
  293. * L'opérateur 'LAPN' 'VF' nous donne la matrice \Delta' et le
  294. * résidu \Delta T_0.
  295. * On n'inverse bien évidemment pas \Delta' mais on résoud le système:
  296. * \Delta' IncT = \Delta T_0
  297. * => IncT = {\Delta'}^{-1} (\Delta T_0)
  298. *
  299. * La méthode de Newton doit converger en un pas (on vérifie que le
  300. * résidu (\Delta T_1) est nul après le premier pas.
  301. *
  302. *
  303. t0 = 'KCHT' $mt 'SCAL' 'CENTRE' 0.0 ;
  304.  
  305. gradt0 mchamt = 'PENT'  $mt
  306.                        'FACE' 'VFSYM' t0
  307.                         'DISPDIF' PERM 'TIMP' solim
  308.                          ;
  309. jaco chpres dt = 'LAPN' 'VF' 'CLAUDEIS' 'IMPL'
  310. $mt t0 gradt0 mchamt 'TIMP' solim ;
  311. mamat = 'KOPS' 'MULT' -1.0D0 jaco ;
  312. matot =  mamat ;
  313.  
  314.  
  315. AUX = ('EXCO' chpres 'RETN' 'RETN') - (NOMC(sour) 'RETN');
  316. rv = 'EQEX' ;
  317. rv . 'METHINV' . 'TYPINV'   = 3 ;
  318. rv . 'METHINV' . 'PRECOND'  = 3 ;
  319. rv . 'METHINV' . 'MATASS'   = matot ;
  320. rv . 'METHINV' . 'MAPREC'   = matot ;
  321. deltat = 'KRES' matot 'TYPI' (rv . 'METHINV')
  322.                       'SMBR' AUX
  323.                       'IMPR' 0 ;
  324. t1 = t0 '+' ('EXCO' deltat 'RETN' 'SCAL') ;
  325. *******************************************************************
  326. * Verification de la convergence
  327. *******************************************************************
  328. gradt1    = 'PENT'  $mt  'FACE' 'VFSYM' t1
  329.                     'DISPDIF' PERM 'TIMP' solim
  330.                     'GRADGEO' mchamt;
  331. jacbid chpres1 dt = 'LAPN' 'VF' 'CLAUDEISl' 'IMPL'
  332. $mt t0 gradt1 mchamt 'TIMP' solim ;
  333. AUX = ('EXCO' chpres1 'RETN' 'RETN') - (NOMC(sour) 'RETN');
  334. res = 'KRES' matot 'TYPI' (rv . 'METHINV')
  335.                       'SMBR'AUX
  336.                       'IMPR' 0 ;
  337. mres1 = 'MAXIMUM' res  'ABS' ;
  338. 'MESSAGE' ('CHAINE' 'Maxi. mres1 =' mres1) ;
  339. 'SI' ('>' mres1 1.e-5) ;
  340.     'MESSAGE' 'La méthode de Newton na pas converge en un pas.'
  341.     'ERREUR' 5 ;
  342. 'FINSI' ;
  343.  
  344.  
  345.  
  346. *
  347. * Résultats
  348. *
  349. 'SI' graph ;
  350. *
  351. * solutions exactes
  352. *
  353.   tn = solexc ;
  354.   chm_tn = 'KCHA' $mt 'CHAM' tn ;
  355.   titt = 'CHAINE' 'Température exacte' titglob ;
  356.   'TRACER' chm_tn modmt 'TITRE' titt ;
  357. *
  358. * graphe de convergence de la méthode itérative
  359. *
  360.     conver = (rv . 'METHINV' . 'CONVINV') ;
  361.     dimcon = 'DIME' conver ;
  362.     'SI' ('>' dimcon 1) ;
  363.         labs = 'PROG' 1.D0 'PAS' 1.D0 dimcon ;
  364.         lord = ('LOG' conver) '/' ('LOG' 10.D0) ;
  365.         evtot = 'EVOL' 'MANU' 'ITER' labs 'RESID' lord ;
  366.         titev = 'CHAINE' 'Historique de convergence' titglob ;
  367.         'DESSIN' evtot 'TITR' titev 'LEGE' ;
  368.     'FINSI' ;
  369. *
  370. * solutions calculées
  371. *
  372.   tn = t1 ;
  373.   chm_tn = 'KCHA' $mt 'CHAM' tn ;
  374.   titt = 'CHAINE' 'Température calculée' titglob ;
  375.   'TRACER' chm_tn modmt  'TITRE' titt ;
  376.  
  377. * erreur
  378. *
  379.   erl = ABS(solexc '-' tn)/(maxi solexc);
  380.   erlcham = 'KCHA' $mt 'CHAM' erl ;  
  381.   titt = 'CHAINE' 'erreur relative' ;
  382.   'TRACER' (erlcham) modmt  'TITRE' titt ;
  383. 'FINSI' ;
  384. *
  385. * Calcul des erreurs par rapport à la solution analytique
  386. *
  387. echloc = '**' (('MESURE' mt) '/' ('NBEL' mt))
  388.               (1.D0 '/' ('VALEUR' 'DIME')) ;
  389. tn = t1 ;
  390. errlit = CALCERR  tn solexc ;
  391. errl2t = CALCERR2 tn solexc ;
  392. 'MESSAGE' '-------------------------------------------------' ;
  393. 'MESSAGE' ('CHAINE' 'Erreur en norme LINF : ' errlit) ;
  394. 'MESSAGE' ('CHAINE' 'Erreur en norme L2   : ' errl2t) ;
  395. 'MESSAGE' '-------------------------------------------------' ;
  396. *'OPTION' 'DONN' 5 ;
  397. 'FINPROC' echloc errl2t ;
  398. * Fin de la procédure de calcul                             *
  399. *___________________________________________________________*
  400. *-----------------------------------------------------------
  401. *-----------------------------------------------------------
  402. * Paramètres du calcul
  403. *
  404. *
  405. * lraff : nb raffinement du maillage (à chaque fois, on découpe
  406. *         les éléments en quatre).
  407. 'SI' complet ;
  408.   lraff  = 'LECT' 0 PAS 1 2 ;
  409. 'SINON' ;
  410.   lraff  = 'LECT' 0 PAS 1 5 ;
  411. 'FINSI' ;
  412. *
  413. *-----------------------------------------------------------*
  414. * Boucles sur les différents paramètres du calcul           *
  415. ordok = VRAI ;
  416. precok = VRAI ;
  417. * ordre théorique en espace de la méthode
  418. ordtth = 2 ;
  419. * erreur L2 max pour la solution (raffinement 1, complet=FAUX)
  420. errtmax = 1.6605D-2 ;
  421. *
  422. *  Boucle sur les raffinements
  423. *
  424.     lh = 'PROG' ;
  425.     lerl2t = 'PROG' ;
  426.     'REPETER' iraff ('DIME' lraff)  ;
  427.       raff = 'EXTRAIRE' lraff &iraff ;
  428.       h errt = CALCUL raff ;
  429.       lh     = 'ET' lh     ('PROG' h ) ;
  430.       lerl2t = 'ET' lerl2t ('PROG' errt) ;
  431.     'FIN' iraff ;
  432.     lh     = ('LOG' lh)     '/' ('LOG' 10.D0) ;
  433.     lerl2t = ('LOG' lerl2t) '/' ('LOG' 10.D0) ;
  434.     tl2 = 'EVOL' 'MANU' 'Long. carac.' lh
  435.                         'Err. tempér.' lerl2t ;
  436. *
  437. * Recherche des ordres de convergence
  438. *
  439.     cpt tlipol = @POMI tl2 1 'IDEM' ;
  440.     ordt = cpt . 1 ;
  441.     'MESSAGE' ('CHAINE' 'ordre température=' ordt) ;
  442. *
  443. * Tracés graphiques
  444. *  
  445.     'SI' graph ;
  446.       tableg = 'TABLE' ;
  447.       tableg . 'TITRE' = 'TABLE' ;
  448.       tableg . 1 = 'MARQ CROI NOLI' ;
  449.       tableg . 'TITRE' . 1 = 'Erreur L2 calculée' ;
  450.       tableg . 2 = 'TIRR' ;                   
  451.       tableg . 'TITRE' . 2 = 'Erreur L2 interpolée' ;
  452.       titdest= 'CHAINE' 'Err L2 Tempér., ordre=' ordt ;
  453.       'DESSIN' (tl2 'ET' tlipol) 'TITRE' titdest tableg ;
  454.     'FINSI' ;
  455.  
  456. *
  457. * Tests des ordres de convergence
  458. * Tests des erreurs L2 sur le maillage le plus fin dans le
  459. * cas complet=faux
  460. *
  461.      ordok = 'ET' ordok ('EGA' ('ENTIER' ('+' ordt 0.5)) ordtth) ;
  462.      'SI' ('EGA' complet FAUX) ;
  463.        precok = 'ET' precok ('<' errt errtmax) ;
  464.      'FINSI' ;
  465.      'SI' ('NON' ordok) ;
  466.        'MESSAGE' 'On n_obtient pas un ordre de convergence correct' ;
  467.        'ERREUR' 5 ;
  468.      'FINSI' ;
  469.      'SI' ('NON' precok) ;
  470.        'MESSAGE' 'On n_obtient pas une précision correcte' ;
  471.        'ERREUR' 5 ;
  472.      'FINSI' ;
  473. 'SI' interact ;
  474.    'OPTION' 'ECHO' 1 ;
  475.    'OPTION' 'DONN' 5 ;
  476. 'FINSI' ;
  477. 'MESSAGE' 'Tout s_est bien passé' ;
  478. 'FIN' ;
  479.  
  480.  
  481.  
  482.  
  483.  
  484.  
  485.  
  486.  
  487.  
  488.  

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