Télécharger equ_chaleur3D_VFconv.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : equ_chaleur3D_VFconv.dgibi
  2. ************************************************************************
  3. ************************************************************************
  4. ************************************************************************
  5. * NOM : equ_chaleur3D_VFconv.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 - UT = 0
  22. *
  23. * |
  24. * avec K = |x+1 z y|
  25. * |z y+1 x|
  26. * |y x z+1|
  27. *
  28. *
  29. * avec U = |1 + x + y + z|
  30. * |1 + x + y + z|
  31. * |1 + x + y + z|
  32. *
  33. * - Conditions aux limites :
  34. *
  35. * conditions de Dirichlet, conditions de flux,
  36. * conditions mixtes calculées à partir de la
  37. * restriction de la solution exacte sur le bord
  38. * exacte
  39. *
  40. * - Solution exacte :
  41. *
  42. * T(x,y)= e(x+y+z)
  43. *
  44. *
  45. *
  46. *
  47. * DISCRETISATION : une méthode de Volume Finis d'ordre 2 en espace.
  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, 23/01/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. *nraff=1 ;
  136. 'DEBPROC' CALCUL nraff*'ENTIER' ;
  137. *
  138. * titre global pour les dessins
  139. *
  140. titglob = 'CHAINE' ' ; nraff=' nraff ;
  141. *
  142. *
  143. * Géométrie
  144. *
  145. pA = 0. 0. 0. ;
  146. pB = 1.0 0. 0. ;
  147. pC = 1.0 1.0 0. ;
  148. pD = 0. 1.0 0. ;
  149.  
  150. pE = 0. 0. 1. ;
  151. pF = 1. 0. 1. ;
  152. pG = 1. 1. 1. ;
  153. pH = 0. 1. 1. ;
  154. *
  155. * Paramètres de la discrétisation de base
  156. *
  157. 'SI' complet ;
  158. nAB = 6 ;
  159. nBC = 6 ;
  160. nCD = 6 ;
  161. nDA = 6 ;
  162. nH = 6 ;
  163. 'SINON' ;
  164. nAB = 2 '*' (nraff '+' 1) ;
  165. nBC = 2 '*' (nraff '+' 1) ;
  166. nCD = 2 '*' (nraff '+' 1) ;
  167. nDA = 2 '*' (nraff '+' 1) ;
  168. nH = 2 '*' (nraff '+' 1) ;
  169. 'FINSI' ;
  170. *
  171. *******************************************************************
  172. * Géométrie discrétisée
  173. *******************************************************************
  174. *
  175. 'OPTION' 'ELEM' 'TRI3' ;
  176. bas = 'DROIT' nAB pA pB ;
  177. droite = 'DROIT' nBC pB pC ;
  178. haut = 'DROIT' nCD pC pD ;
  179. gauche = 'DROIT' nDA pD pA ;
  180. spourt = bas 'ET' droite 'ET' haut 'ET' gauche ;
  181. sbas = 'SURFACE' spourt 'PLAN' ;
  182.  
  183. 'OPTION' 'DIME' 3 'ELEM' TET4 ;
  184. mt = 'VOLUME' sbas 'TRAN' nH pE ;
  185. elim mt 0.001;
  186. f1 f2 f3 = 'FACE' mt ;
  187. ss = (f1 'ET' f2 'ET' f3) ;
  188. 'OPTION' 'ELEM' 'TET4' ;
  189. mt = 'VOLUME' ss ;
  190. pourtour = 'ENVELOPPE' mt ;
  191. *
  192. * Eventuellement, on trace le résultat
  193. *
  194. * list mt;
  195. 'SI' graph ;
  196. titgeo = 'CHAINE' 'mt ' 'NBPO=' ('NBNO' mt)
  197. ' NBELEM=' ('NBEL' mt) titglob ;
  198. 'TRACER' mt 'NOEUD' 'TITRE' titgeo ;
  199. 'FINSI' ;
  200. *******************************************************************
  201. *
  202. * Modèle
  203. *
  204. *******************************************************************
  205. Mmt = 'CHAN' mt 'QUAF' ;
  206. 'ELIM' Mmt 1.E-5 ;
  207. MDNS = 'NAVIER_STOKES' ;
  208. $mt = 'MODE' Mmt MDNS 'LINE' ;
  209. mt = 'DOMA' $mt 'MAILLAGE' ;
  210. MESS 'NOMBRE DE MAILLES' (NBEL (MT));
  211.  
  212. * list ('DOMA' $mt 'FACEP') ;
  213. mtm = ('DOMA' $mt 'MAILLAGE') ;
  214. mtc = 'COULEUR' ('DOMA' $mt 'CENTRE') 'JAUN' ;
  215. mtf = 'COULEUR' ('DOMA' $mt 'FACE') 'BLEU' ;
  216. *cmt = 'DOMA' $mt 'ENVELOPP' ;
  217. cmt = pourtour;
  218. Mcmt = 'CHAN' cmt 'QUAF' ;
  219. 'ELIM' (Mmt 'ET' Mcmt) 1.E-5 ;
  220. $cmt = 'MODE' Mcmt MDNS 'QUAF';
  221.  
  222. * DEFINITION DES DIFFERENTS BORD DE LA FRONTIERE
  223. bas1 = cmt 'ELEM' 'APPUYE' 'STRICTEMENT' f3 ;
  224. gau1 = cmt 'ELEM' 'APPUYE' 'STRICTEMENT' (f1) ;
  225. titgeo = 'CHAINE' 'mt ' 'NBPO=' ('NBNO' mt);
  226. 'SI' graph ;
  227. trac bas1 'TITRE' 'bas1';
  228. trac gau1 'TITRE' 'gau1';
  229. 'FINSI' ;
  230. basgau = ( bas1 'ET' gau1);
  231. 'SI' graph ;
  232. trac basgau 'TITRE' 'basgau';
  233. 'FINSI' ;
  234. cvn = 'DIFF' cmt basgau ;
  235. 'SI' graph ;
  236. trac cvn 'TITRE' 'cvn';
  237. 'FINSI' ;
  238. basgau = ( bas1 'ET' gau1);
  239. cvn = cmt 'ELEM' 'APPUYE' 'STRICTEMENT' f1 ;
  240. bas1 = cmt 'ELEM' 'APPUYE' 'STRICTEMENT' f2 ;
  241. basgau = ( bas1 'ET' cvn);
  242. gau1 = diff cmt basgau ;
  243. Mbas = 'CHAN' bas1 'QUAF' ;
  244. Mgau = 'CHAN' gau1 'QUAF' ;
  245. Mcvn = 'CHAN' cvn 'QUAF' ;
  246. Mf1 = 'CHAN' f1 'QUAF' ;
  247. 'ELIM' (Mmt 'ET' Mbas 'ET' Mgau 'ET' Mcvn 'ET' Mf1) 1.E-5 ;
  248. $bas = 'MODE' Mbas MDNS 'QUAF' ;
  249. $gau = 'MODE' Mgau MDNS 'QUAF' ;
  250. $cvn = 'MODE' Mcvn MDNS 'QUAF' ;
  251. $f1 = 'MODE' Mf1 MDNS 'QUAF';
  252. cvn = 'DOMA' $cvn 'CENTRE' ;
  253. cvn1 = 'DOMA' $cvn 'MAILLAGE' ;
  254. bas = 'DOMA' $bas 'MAILLAGE' ;
  255. gau = 'DOMA' $gau 'CENTRE' ;
  256. *
  257.  
  258. 'SI' graph ;
  259. pdir1ma = ('DOMA' $bas 'MAILLAGE');
  260. pdir2ma = ('DOMA' $gau 'MAILLAGE');
  261. cvnma = ('DOMA' $cvn 'MAILLAGE');
  262. mai = (cvnma 'COULEUR' 'VERT') et (pdir1ma 'COULEUR' ROUGE)
  263. et (pdir2ma 'COULEUR' BLEU);
  264.  
  265. 'TRACER' (mai) 'TITRE' 'CONDITIONS AUX LIMITES' ;
  266. 'FINSI' ;
  267.  
  268. *******************************************************************
  269. * Solution exacte aux centres
  270. *******************************************************************
  271. xxc yyc zzc = 'COORDONNEE' ('DOMA' $mt 'CENTRE') ;
  272. solexc = (exp ((+1.D0)*(xxc + yyc + zzc))) ;
  273. solexc = 'KCHT' $mt 'SCAL' 'CENTRE' solexc ;
  274.  
  275.  
  276. *******************************************************************
  277. * terme source
  278. *******************************************************************
  279. sour = 0.D0 + (0.0*xxc);
  280. sour = 'KCHT' $mt 'SCAL' 'CENTRE' sour ;
  281.  
  282. *******************************************************************
  283. * calcul de la convection
  284. *******************************************************************
  285. CHYB2 = 'DOMA' $mt 'NORMALE' ;
  286. xxg yyg zzg = 'COORDONNEE' ('DOMA' $mt 'FACE') ;
  287. ux = (1.D0 + xxg + yyg + zzg);
  288. uy = (1.D0 + xxg + yyg + zzg) ;
  289. uz = (1.D0 + xxg + yyg + zzg) ;
  290. ux = 'NOMC' 'UX' ux 'NATU' 'DISCRET' ;
  291. uy = 'NOMC' 'UY' uy 'NATU' 'DISCRET' ;
  292. uz = 'NOMC' 'UZ' uZ 'NATU' 'DISCRET' ;
  293. qchal = (ux + uy + uz);
  294. MOT1 = 'MOTS' 'UX' 'UY' 'UZ' ;
  295. MOT2 = 'MOTS' 'UX' 'UY' 'UZ' ;
  296. FLU1 = 'PSCA' CHYB2 qchal MOT1 MOT2 ;
  297. FLUC = NOMC(FLU1) 'FLUX';
  298. *******************************************************************
  299. * CONDITIONS DE NEUMANN
  300. *******************************************************************
  301. xxf yyf zzf = 'COORDONNEE' ('DOMA' $gau 'CENTRE') ;
  302. gradtx = (+1.D0)*(1.D0 + (xxf + yyf + zzf)) *
  303. (exp ((+1.D0)*(xxf + yyf + zzf))) ;
  304. gradty = (+1.D0)*(1.D0 + (xxf+ yyf + zzf )) *
  305. (exp ((+1.D0)*(xxf + yyf + zzf))) ;
  306. gradtz = (+1.D0)*(1.D0 + (xxf+ yyf + zzf )) *
  307. (exp ((+1.D0)*(xxf + yyf + zzf))) ;
  308. gradtx = 'NOMC' 'UX' gradtx 'NATU' 'DISCRET' ;
  309. gradty = 'NOMC' 'UY' gradty 'NATU' 'DISCRET' ;
  310. gradtz = 'NOMC' 'UZ' gradtz 'NATU' 'DISCRET' ;
  311. CHYB2 = 'DOMA' $mt 'NORMALE' ;
  312. qchal = -1.D0 * (gradtx + gradty + gradtz) ;
  313. qlim2 = qchal;
  314. CHYB2R = 'REDU' CHYB2 gau ;
  315. MOT1 = 'MOTS' 'UX' 'UY' 'UZ' ;
  316. MOT2 = 'MOTS' 'UX' 'UY' 'UZ' ;
  317. FLU1 = 'PSCA' CHYB2R qlim2 MOT1 MOT2 ;
  318. qlim2 = NOMC(FLU1) 'FLUX';
  319. *******************************************************************
  320. * CONDITION DE DIRICHLET
  321. *******************************************************************
  322. xxlim yylim zzlim = 'COORDONNEE' ('DOMA' $bas 'CENTRE') ;
  323. solim = (exp ((+1.D0)*(xxlim + yylim + zzlim))) ;
  324. solim = 'KCHT' $bas 'SCAL' 'CENTRE' solim ;
  325.  
  326. *******************************************************************
  327. * CONDITIONS MIXTES
  328. *******************************************************************
  329. * DEFINITION DES PARAMETRES POUR LES CONDITIONS MIXTES
  330. xxm yym zzm = 'COORDONNEE' ('DOMA' $cvn 'CENTRE') ;
  331. gradmx = (+1.D0)*(1.D0 + (xxm + yym + zzm)) *
  332. (exp ((+1.D0)*(xxm + yym + zzm))) ;
  333. gradmx = 'NOMC' 'UX' gradmx 'NATU' 'DISCRET' ;
  334. gradmy = (+1.D0)*(1.D0 + (xxm + yym + zzm)) *
  335. (exp ((+1.D0)*(xxm + yym + zzm))) ;
  336. gradmy = 'NOMC' 'UY' gradmy 'NATU' 'DISCRET' ;
  337. gradmz = (+1.D0)*(1.D0 + (xxm + yym + zzm)) *
  338. (exp ((+1.D0)*(xxm + yym + zzm))) ;
  339. gradmz = 'NOMC' 'UZ' gradmz 'NATU' 'DISCRET' ;
  340. AUX = 1.D0 + (0.0D0*xxm);
  341. XPAR1 = 'NOMC' 'PAR1' AUX 'NATU' 'DISCRET' ;
  342. XPAR1 = 'REDU' XPAR1 cvn ;
  343. qchal = (-1.0) * (gradmx + gradmy + gradmz) ;
  344. qaux = 'REDU' qchal cvn ;
  345. qlim = (1.D0)* (qaux);
  346. CHYB2R = 'REDU' CHYB2 cvn ;
  347. MOT1 = 'MOTS' 'UX' 'UY' 'UZ' ;
  348. MOT2 = 'MOTS' 'UX' 'UY' 'UZ' ;
  349. FLU1 = 'PSCA' CHYB2R qlim MOT1 MOT2 ;
  350. qlim = NOMC(FLU1) 'FLUX';
  351. MOT1 = 'MOTS' 'UX' 'UY' 'UZ' ;
  352. MOT2 = 'MOTS' 'UX' 'UY' 'UZ' ;
  353. FLU1 = 'PSCA' CHYB2R qaux MOT1 MOT2 ;
  354. XPAR2 = 'NOMC' FLU1 'PAR2';
  355. XPAR2= (0.0D0*XPAR2);
  356. * ON SIMULE UNE CONDITON DE NEUMANN VIA UNE CONDITION MIXTE
  357. *XPAR2 doit ETRE POSITIF
  358. VAL = XPAR1 + XPAR2 + QLIM;
  359.  
  360. *******************************************************************
  361. *
  362. * Calcul du tenseur
  363. *
  364. *******************************************************************
  365. K11 = 'KOPS' xxc '+' 1.D0 ;
  366. K22 = 'KOPS' yyc '+' 1.D0 ;
  367. K33 = 'KOPS' zzc '+' 1.D0 ;
  368. K21 = zzc ;
  369. K31 = (yyc) ;
  370. K32 = (xxc) ;
  371.  
  372. K11 = NOMC(K11) 'K11' 'NATU' 'DISCRET';
  373. K22 = NOMC(K22) 'K22' 'NATU' 'DISCRET';
  374. K33 = NOMC(K33) 'K33' 'NATU' 'DISCRET';
  375. K21 = NOMC(K21) 'K21' 'NATU' 'DISCRET';
  376. K31 = NOMC(K31) 'K31' 'NATU' 'DISCRET';
  377. K32 = NOMC(K32) 'K32' 'NATU' 'DISCRET';
  378.  
  379. PERM = (K11 'ET' K22) ;
  380. PERM = PERM 'ET' K33;
  381. PERM = PERM 'ET' K21;
  382. PERM = PERM 'ET' K31;
  383. PERM = PERM 'ET' K32;
  384.  
  385. * Mise en place du calcul numérique
  386. *
  387. * équation de Laplace
  388. *
  389. *
  390. * on utilise une méthode de Newton pour résoudre :
  391. * - \Delta T = 0 (\Delta opérateur laplacien)
  392. * - avec T_{\partial \Omega} donné (CL de Dirichlet)
  393. * T_0 : estimation initiale de la solution
  394. * On a T_1 = T_0 - {\Delta'}^{-1} (\Delta T_0)
  395. *
  396. * L'opérateur 'LAPN' 'VF' nous donne la matrice \Delta' et le
  397. * résidu \Delta T_0.
  398. * On n'inverse bien évidemment pas \Delta' mais on résoud le système:
  399. * \Delta' IncT = \Delta T_0
  400. * => IncT = {\Delta'}^{-1} (\Delta T_0)
  401. *
  402. * La méthode de Newton doit converger en un pas (on vérifie que le
  403. * résidu (\Delta T_1) est nul après le premier pas.
  404. *
  405. *
  406. t0 = solexc;
  407. t0 = 'KCHT' $mt 'SCAL' 'CENTRE' 0.0 ;
  408. menage;
  409. gradt0 mchamt = 'PENT' $mt
  410. 'FACE' 'MPFA' t0
  411. 'DISPDIF' PERM 'TIMP'
  412. solim 'QIMP' qlim2 'MIXT' val
  413. 'CENTERED' FLUC;
  414. qaux = 'REDU' gradt0 cvn ;
  415. *list qaux;
  416. jaco chpres dt = 'LAPN' 'VF' 'CLAUDEIS' 'IMPL'
  417. $mt t0 gradt0 mchamt 'QIMP' qlim2 'MIXT' val 'TIMP' solim ;
  418. mamat = 'KOPS' 'MULT' -1.0D0 jaco ;
  419. matot = mamat ;
  420.  
  421.  
  422. AUX = ('EXCO' chpres 'RETN' 'RETN') - (NOMC(sour) 'RETN');
  423. rv = 'EQEX' ;
  424. rv . 'METHINV' . 'TYPINV' = 3 ;
  425. rv . 'METHINV' . 'PRECOND' = 3 ;
  426. rv . 'METHINV' . 'MATASS' = matot ;
  427. rv . 'METHINV' . 'MAPREC' = matot ;
  428. deltat = 'KRES' matot 'TYPI' (rv . 'METHINV')
  429. 'SMBR' AUX
  430. 'IMPR' 0 ;
  431. t1 = t0 '+' ('EXCO' deltat 'RETN' 'SCAL') ;
  432. *******************************************************************
  433. * Verification de la convergence
  434. *******************************************************************
  435. gradt1 mchamt = 'PENT' $mt 'FACE' 'MPFA' t1
  436. 'DISPDIF' PERM 'TIMP' solim
  437. 'QIMP' qlim2 'MIXT' val
  438. 'CENTERED' FLUC ;
  439. jacbid chpres1 dt = 'LAPN' 'VF' 'CLAUDEISl' 'IMPL'
  440. $mt t0 gradt1 mchamt 'QIMP' qlim2 'MIXT' val 'TIMP' solim ;
  441. AUX = ('EXCO' chpres1 'RETN' 'RETN') - (NOMC(sour) 'RETN');
  442. res = 'KRES' matot 'TYPI' (rv . 'METHINV')
  443. 'SMBR'AUX
  444. 'IMPR' 0 ;
  445. mres1 = 'MAXIMUM' res 'ABS' ;
  446. 'MESSAGE' ('CHAINE' 'Maxi. mres1 =' mres1) ;
  447. 'SI' ('>' mres1 1.e-5) ;
  448. 'MESSAGE' 'La méthode de Newton na pas converge en un pas.'
  449. 'ERREUR' 5 ;
  450. 'FINSI' ;
  451. *
  452. * Résultats
  453. *
  454. 'SI' graph ;
  455. *
  456. * solutions exactes
  457. *
  458. tn = solexc ;
  459. chm_tn = 'KCHA' $mt 'CHAM' tn ;
  460. titt = 'CHAINE' 'Température exacte' titglob ;
  461. 'TRACER' chm_tn $mt 'TITRE' titt ;
  462. *
  463. * graphe de convergence de la méthode itérative
  464. *
  465. conver = (rv . 'METHINV' . 'CONVINV') ;
  466. dimcon = 'DIME' conver ;
  467. 'SI' ('>' dimcon 1) ;
  468. labs = 'PROG' 1.D0 'PAS' 1.D0 dimcon ;
  469. lord = ('LOG' conver) '/' ('LOG' 10.D0) ;
  470. evtot = 'EVOL' 'MANU' 'ITER' labs 'RESID' lord ;
  471. titev = 'CHAINE' 'Historique de convergence' titglob ;
  472. 'DESSIN' evtot 'TITR' titev 'LEGE' ;
  473. 'FINSI' ;
  474. *
  475. * solutions calculées
  476. *
  477. tn = t1 ;
  478. chm_tn = 'KCHA' $mt 'CHAM' tn ;
  479. titt = 'CHAINE' 'Température calculée' titglob ;
  480. 'TRACER' chm_tn $mt 'TITRE' titt ;
  481. *
  482. * Calcul des erreurs par rapport à la solution analytique
  483. *
  484. * erreur
  485. *
  486. erl = ABS(solexc '-' tn)/(maxi solexc);
  487. erlcham = 'KCHA' $mt 'CHAM' erl ;
  488. titt = 'CHAINE' 'erreur relative' ;
  489. 'TRACER' (erlcham) $mt 'TITRE' titt ;
  490. 'FINSI' ;
  491. echloc = '**' (('MESURE' mt) '/' ('NBEL' mt))
  492. (1.D0 '/' ('VALEUR' 'DIME')) ;
  493. tn = t1 ;
  494. errlit = CALCERR tn solexc ;
  495. errl2t = CALCERR2 tn solexc ;
  496. 'MESSAGE' '-------------------------------------------------' ;
  497. 'MESSAGE' ('CHAINE' 'Erreur en norme LINF : ' errlit) ;
  498. 'MESSAGE' ('CHAINE' 'Erreur en norme L2 : ' errl2t) ;
  499. 'MESSAGE' '-------------------------------------------------' ;
  500. *'OPTION' 'DONN' 5 ;
  501. 'FINPROC' echloc errl2t ;
  502. * Fin de la procédure de calcul *
  503. *___________________________________________________________*
  504. *-----------------------------------------------------------
  505. *-----------------------------------------------------------
  506. * Paramètres du calcul
  507. *
  508. *
  509. * lraff : nb raffinement du maillage (à chaque fois, on découpe
  510. * les éléments en quatre).
  511. 'SI' complet ;
  512. lraff = 'LECT' 0 PAS 1 2 ;
  513. 'SINON' ;
  514. lraff = 'LECT' 1 PAS 1 4 ;
  515. 'FINSI' ;
  516. *
  517. *-----------------------------------------------------------*
  518. * Boucles sur les différents paramètres du calcul *
  519. ordok = VRAI ;
  520. precok = VRAI ;
  521. * ordre théorique en espace de la méthode
  522. ordtth = 2 ;
  523. * erreur L2 max pour la solution (raffinement 1, complet=FAUX)
  524. errtmax = 1.D-2 ;
  525. *
  526. * Boucle sur les raffinements
  527. *
  528. lh = 'PROG' ;
  529. lerl2t = 'PROG' ;
  530. 'REPETER' iraff ('DIME' lraff) ;
  531. raff = 'EXTRAIRE' lraff &iraff ;
  532. h errt = CALCUL raff ;
  533. lh = 'ET' lh ('PROG' h ) ;
  534. lerl2t = 'ET' lerl2t ('PROG' errt) ;
  535. 'FIN' iraff ;
  536. lh = ('LOG' lh) '/' ('LOG' 10.D0) ;
  537. lerl2t = ('LOG' lerl2t) '/' ('LOG' 10.D0) ;
  538. tl2 = 'EVOL' 'MANU' 'Long. carac.' lh
  539. 'Err. tempér.' lerl2t ;
  540. *
  541. * Recherche des ordres de convergence
  542. *
  543. cpt tlipol = @POMI tl2 1 'IDEM' ;
  544. ordt = cpt . 1 ;
  545. 'MESSAGE' ('CHAINE' 'ordre température=' ordt) ;
  546. *
  547. * Tracés graphiques
  548. *
  549. 'SI' graph ;
  550. tableg = 'TABLE' ;
  551. tableg . 'TITRE' = 'TABLE' ;
  552. tableg . 1 = 'MARQ CROI NOLI' ;
  553. tableg . 'TITRE' . 1 = 'Erreur L2 calculée' ;
  554. tableg . 2 = 'TIRR' ;
  555. tableg . 'TITRE' . 2 = 'Erreur L2 interpolée' ;
  556. titdest= 'CHAINE' 'Err L2 Tempér., ordre=' ordt ;
  557. 'DESSIN' (tl2 'ET' tlipol) 'TITRE' titdest tableg ;
  558. 'FINSI' ;
  559. *
  560. * Tests des ordres de convergence
  561. * Tests des erreurs L2 sur le maillage le plus fin dans le
  562. * cas complet=faux
  563. *
  564. ordok = 'ET' ordok ('EGA' ('ENTIER' ('+' ordt 0.5)) ordtth) ;
  565. 'SI' ('EGA' complet FAUX) ;
  566. precok = 'ET' precok ('<' errt errtmax) ;
  567. 'FINSI' ;
  568. 'SI' ('NON' ordok) ;
  569. 'MESSAGE' 'On n_obtient pas un ordre de convergence correct' ;
  570. 'ERREUR' 5 ;
  571. 'FINSI' ;
  572. 'SI' ('NON' precok) ;
  573. 'MESSAGE' 'On n_obtient pas une précision correcte' ;
  574. 'ERREUR' 5 ;
  575. 'FINSI' ;
  576. 'SI' interact ;
  577. 'OPTION' 'ECHO' 1 ;
  578. 'OPTION' 'DONN' 5 ;
  579. 'FINSI' ;
  580. 'MESSAGE' 'Tout s_est bien passé' ;
  581. 'FIN' ;
  582.  
  583.  
  584.  
  585.  
  586.  
  587.  
  588.  
  589.  
  590.  
  591.  

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