Télécharger la plaquette Cast3M


fr en

Télécharger equ_chaleur3D_VF.dgibi

Retour à la liste

Numérotation des lignes : 

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

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