Télécharger poiseuille2D.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : poiseuille2D.dgibi
  2. ************************************************************************
  3. * NOM : poiseuille2D.dgibi
  4. * ___
  5. *
  6. * DESCRIPTION : Ecoulement de poiseuille 2D classique
  7. * ___________ (écoulement parallèle en conduit bidimensionnel infini)
  8. *
  9. * GEOMETRIE (problème dimensionnel) :
  10. * -----------------------------------
  11. *
  12. * y
  13. * ^ vitesse=0
  14. * --------------|---------------------- y = b ------------
  15. * Pression=P0| |
  16. * | | Pression = 0
  17. * | |
  18. * -----------------------------> x
  19. * O|x = 0 |x = a
  20. * | |
  21. * | |
  22. * ------------------------------------- y = -b -----------
  23. * vitesse=0
  24. *
  25. *
  26. * EQUATIONS ADIMENSIONNEES :
  27. * --------------------------
  28. *
  29. * - Echelles :
  30. *
  31. * longueur (suivant x) : a
  32. * longueur (suivant y) : b
  33. * pression : p0
  34. * vitesse : ((b*b) / (2*mu)) * (p0 / a)
  35. * (mu : viscosité dynamique)
  36. *
  37. * - Equations : Navier-Stokes se réduit à
  38. *
  39. *
  40. * (1/2) * (d²u/dy²) = dp/dx avec : u=u(y) et p=p(x)
  41. *
  42. * - Conditions aux limites :
  43. *
  44. *
  45. * u(b) = u(-b) = 0 (non-glissement)
  46. *
  47. * - Solution exacte :
  48. *
  49. * On déduit des équations : dp/dx = (1/2) * (d²u/dy²) = constante
  50. * On impose donc le gradient de pression et on trouve (avec
  51. * l'adimensionnement ci-dessus, constante=-1) :
  52. *
  53. * u = 1 - y²
  54. *
  55. *
  56. * DISCRETISATION : On aimerait les tester toutes mais on rencontre les
  57. * ______________ problèmes suivants (au 01/03/1999):
  58. *
  59. * - TOIMP fonctionne mal : on ne peut pas tourner le
  60. * maillage et on doit imposer une contrainte pipo
  61. * sur le bord gauche : (0 p0) au lieu de (0 -p0).
  62. *
  63. * - Les éléments MACRO triangulaires ne fonctionnent
  64. * pas.
  65. *
  66. * Pour l'instant, on se cantonne à la discrétisation
  67. * suivante :
  68. * - éléments quadratiques (triangle Crouzeix-Raviart
  69. * et quadrangle Bercovier-Pironneau)
  70. * - discrétisation centrée des termes de convection
  71. * (qui sont analytiquement nulles à convergence
  72. * -------------
  73. * sur les itérations pour résoudre la non-linéarité)
  74. *
  75. * Le maillage est construit avec l'opérateur SURF, il est perturbé
  76. * par un bruit blanc gaussien. Il est ensuite translaté d'un
  77. * vecteur arbitraire et tourné
  78. * d'un angle arbitraire.
  79. *
  80. *
  81. * Opérateurs utilisés : KBBT, NS, TOIMP
  82. * options : EF, IMPL, CENTREE
  83. *
  84. * RESOLUTION : - Solveur direct Crout
  85. * __________ - Solveurs itératifs pour matrice non-symétrique :
  86. * * Bi-CGSTAB
  87. * * GMRES(50)
  88. * NB : - Pour les méthodes itératives, on garde le premier
  89. * préconditionneur (Choleski incomplet) calculé pour la boucle
  90. * sur les non-linéarités.
  91. *
  92. *
  93. * TESTS EFFECTUES : - On vérifie en gros si l'ordre de convergence est
  94. * _______________ respecté (car le maillage reste grossier donc
  95. * l'ordre est approché)
  96. * - On compare l'écart en norme L2 à la solution
  97. * analytique pour le maillage le plus fin.
  98. *
  99. *
  100. *
  101. * LANGAGE : GIBIANE-CASTEM 2000
  102. * AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF)
  103. * mél : gounand@semt2.smts.cea.fr
  104. ************************************************************************
  105. * VERSION : v1, 09/02/99, version initiale
  106. * HISTORIQUE : v1, 09/02/99, création
  107. * HISTORIQUE : 14/06/99, modif : TOIMP fonctionne comme
  108. * le dit la notice (sauf pour le signe).
  109. * HISTORIQUE :
  110. ************************************************************************
  111. * Prière de PRENDRE LE TEMPS de compléter les commentaires
  112. * en cas de modification de ce sous-programme afin de faciliter
  113. * la maintenance !
  114. ************************************************************************
  115. *
  116. interact= FAUX ;
  117. complet = FAUX ;
  118. graph = FAUX ;
  119. *
  120. 'OPTION' 'DIME' 2 'ELEM' QUA4 'MODE' 'PLAN' ;
  121. 'OPTION' 'ECHO' 0 ;
  122. 'OPTION' 'ISOV' 'SULI' ;
  123. nbisov = 15 ;
  124. 'SI' ('NON' interact) ;
  125. 'OPTION' 'TRAC' 'PS' ;
  126. 'OPTION' 'ISOV' 'LIGNE' ;
  127. 'FINSI' ;
  128. *
  129. ** Erreur Linfini entre deux Champoints.
  130. *
  131. DEBPROC CALCERR vitp1*'CHPOINT' vit*'CHPOINT' ;
  132. err = MAXI (vitp1 - vit) 'ABS' ;
  133. FINPROC err ;
  134. *
  135. ** Erreur Pseudo L2 entre deux Champoints.
  136. *
  137. DEBPROC CALCERR2 vitp1*'CHPOINT' vit*'CHPOINT' ;
  138. er0 = 'KOPS' vitp1 '-' vit ;
  139. er2 = 'KOPS' er0 '*' er0 ;
  140. err = 0.D0 ;
  141. suppv = 'EXTRAIRE' er2 'MAIL' ;
  142. compv = 'EXTRAIRE' er2 'COMP' ;
  143. nbcomp = 'DIME' compv ;
  144. nptd = 'NBNO' suppv ;
  145. 'REPETER' numcomp nbcomp ;
  146. lacomp = 'EXTRAIRE' compv &numcomp ;
  147. 'REPETER' numpo nptd ;
  148. lepoi = suppv 'POIN' &numpo ;
  149. err = err '+' ('EXTRAIRE' er2 lacomp lepoi) ;
  150. 'FIN' numpo ;
  151. 'FIN' numcomp ;
  152. err = err '/' nptd ;
  153. err = err ** 0.5D0 ;
  154. FINPROC err ;
  155. *
  156. * Procédure paramétrée (raffinement, type de discrétisation)
  157. * renvoyant l'échelle de longueur locale, les erreurs en
  158. * norme L2 sur la vitesse et la pression.
  159. * L'écoulement calculé est un écoulement de Poiseuille 2D.
  160. *
  161. 'DEBPROC' CALCUL ktypi*'ENTIER' nraff*'ENTIER' typdis*'MOT ' ;
  162. *
  163. ttypi = 'TABLE' ;
  164. ttypi . 1 = 'CHAINE' 'Crout' ;
  165. ttypi . 3 = 'CHAINE' 'BiCGSTAB' ;
  166. ttypi . 5 = 'CHAINE' 'GMRES(50)' ;
  167. *
  168. * titre global pour les dessins
  169. *
  170. titglob = 'CHAINE' ' ; nraff=' nraff ' typdis=' typdis ' typinv='
  171. (ttypi . ktypi);
  172. 'MESSAGE' '----------------------------------------------' ;
  173. 'MESSAGE' ('CHAINE' 'Cas' titglob) ;
  174. *
  175. * Paramètres physiques
  176. *
  177. mu = 1 ;
  178. *
  179. * Les échelles du problème adimensionné
  180. *
  181. * Echelle de longueur suivant x et y
  182. a = 1.D0 ;
  183. b = 1.D0 ;
  184. * Echelle de pression + contrainte en entrée et en sortie
  185. * du tube
  186. p0 = 1.D0 ;
  187. *!
  188. *! Si TOIMP était conforme à la notice, on devrait avoir
  189. *! qqch comme cgau = (0. -1.*p0)
  190. *!
  191. cgau = (0 p0) ;
  192. cdro = (0.D0 0.D0) ;
  193. * Echelle de vitesse
  194. u0 = ((b * b) '/' (2 * mu)) * (p0 '/' a) ;
  195. *
  196. * Géométrie
  197. *
  198. pA = 0. -1. ;
  199. pB = 1. -1. ;
  200. pC = 1. 1. ;
  201. pD = 0. 1. ;
  202. *
  203. * Paramètres de la discrétisation de base
  204. *
  205. 'SI' complet ;
  206. nAB = 4 ;
  207. nBC = 5 ;
  208. nCD = 6 ;
  209. nDA = 7 ;
  210. 'SINON' ;
  211. nAB = 1 ;
  212. nBC = 2 ;
  213. nCD = 3 ;
  214. nDA = 4 ;
  215. 'FINSI' ;
  216. *
  217. * Rotation et translation aditionnelle + bruit blanc
  218. * + raffinement
  219. *
  220. vtran = ((** 2 0.5) (* -1 (** 3 0.5))) ;
  221. orig = (0.D0 0.D0) ;
  222. arot = 11.3D0 ;
  223. lesdens = 'PROG' (1. '/' nAB) (1. '/' nCD) (1. '/' nBC) (1. '/' nDA) ;
  224. ampbruit = (0.1 * ('MINIMUM' lesdens)) ;
  225. *
  226. * Géométrie discrétisée
  227. *
  228. bas = 'DROIT' nAB pA pB ;
  229. droite = 'DROIT' nBC pB pC ;
  230. haut = 'DROIT' nCD pC pD ;
  231. gauche = 'DROIT' nDA pD pA ;
  232. pourtour = bas 'ET' droite 'ET' haut 'ET' gauche ;
  233. mt = 'SURFACE' pourtour 'PLAN' ;
  234. 'TASSER' mt ;
  235. mt = 'ORIENTER' mt ;
  236. *
  237. * On rajoute le bruit
  238. *
  239. pmt = 'CHANGER' mt 'POI1' ;
  240. pcnt= 'CHANGER' pourtour 'POI1' ;
  241. inmt= 'DIFF' pmt pcnt;
  242. brbl1 = 'BRUIT' 'BLAN' 'GAUS' 0.D0 ampbruit inmt ;
  243. brbl2 = 'BRUIT' 'BLAN' 'GAUS' 0.D0 ampbruit inmt ;
  244. brbl = 'ET' ('NOMC' 'UX' brbl1) ('NOMC' 'UY' brbl2) ;
  245. mt = 'PLUS' mt brbl;
  246. *
  247. * On raffine nraff fois
  248. *
  249. 'SI' (nraff > 0) ;
  250. 'REPETER' bcl nraff ;
  251. mt = 'CHANGER' mt 'QUADRATIQUE' ;
  252. $mt = 'DOMA' mt 'MACRO' ;
  253. mt = ($mt . 'MAILLAGE') ;
  254. bas = 'CHANGER' bas 'QUADRATIQUE' ;
  255. $bas = 'DOMA' bas 'MACRO' ;
  256. bas = ($bas . 'MAILLAGE') ;
  257. droite = 'CHANGER' droite 'QUADRATIQUE' ;
  258. $droite = 'DOMA' droite 'MACRO' ;
  259. droite = ($droite . 'MAILLAGE') ;
  260. haut = 'CHANGER' haut 'QUADRATIQUE' ;
  261. $haut = 'DOMA' haut 'MACRO' ;
  262. haut = ($haut . 'MAILLAGE') ;
  263. gauche = 'CHANGER' gauche 'QUADRATIQUE' ;
  264. $gauche = 'DOMA' gauche 'MACRO' ;
  265. gauche = ($gauche . 'MAILLAGE') ;
  266. 'MENAGE' ;
  267. 'FIN' bcl ;
  268. 'FINSI' ;
  269. *
  270. * Eventuellement, on trace le résultat
  271. *
  272. 'SI' graph ;
  273. titgeo = 'CHAINE' 'mt ' 'NBPO=' ('NBNO' mt)
  274. ' NBELEM=' ('NBEL' mt) titglob ;
  275. 'TRACER' mt 'NOEUD' 'TITRE' titgeo ;
  276. 'FINSI' ;
  277. *
  278. * Modèle
  279. *
  280. _mt = 'CHANGER' mt 'QUAF' ;
  281. _bas = 'CHANGER' bas 'QUAF' ;
  282. _droite = 'CHANGER' droite 'QUAF' ;
  283. _haut = 'CHANGER' haut 'QUAF' ;
  284. _gauche = 'CHANGER' gauche 'QUAF' ;
  285. 'ELIMINATION' 1.D-6 (_gauche 'ET' _haut 'ET' _droite 'ET' _bas
  286. 'ET' _mt) ;
  287. $mt = 'MODELISER' _mt 'NAVIER_STOKES' typdis ;
  288. $gauche = 'MODELISER' _gauche 'NAVIER_STOKES' typdis ;
  289. $droite = 'MODELISER' _droite 'NAVIER_STOKES' typdis ;
  290. *
  291. * Solution exacte
  292. *
  293. chy = 'COORDONNEE' 2 ('DOMA' $mt 'SOMMET') ;
  294. uyexa = 'NOMC' 'UY' ('KCHT' $mt 'SCAL' 'SOMMET' 0.D0) 'NATU' 'DISCRET' ;
  295. uxexa = 'NOMC' 'UX' ( '-' ('KCHT' $mt 'SCAL' 'SOMMET' 1.D0)
  296. (chy * chy))
  297. 'NATU' 'DISCRET' ;
  298. uexact = 'KCHT' $mt 'VECT' 'SOMMET' (uxexa 'ET' uyexa) ;
  299. chxcp1= 'COORDONNEE' 1 ('DOMA' $mt 'CENTREP1') ;
  300. pexa = 'NOMC' 'SCAL' ( '-' ('KCHT' $mt 'SCAL' 'CENTREP1' 1.D0)
  301. (chxcp1)) ;
  302. pexact = 'KCHT' $mt 'SCAL' 'CENTREP1' pexa ;
  303. *
  304. * On tourne et on translate
  305. *
  306. _tmt _tbas _tdroite _thaut _tgauche uexact pexact = 'TOURNER'
  307. _mt _bas _droite _haut _gauche uexact pexact arot orig ;
  308. * On n'a pas besoin de tourner les contraintes imposées avec TOIM
  309. * car elles sont exprimées dans le repère local.
  310. *cgau cdro = 'TOURNER' cgau cdro arot orig ;
  311. _tmt _tbas _tdroite _thaut _tgauche uexact pexact = 'PLUS'
  312. _tmt _tbas _tdroite _thaut _tgauche uexact pexact vtran ;
  313. $tmt = 'MODELISER' _tmt 'NAVIER_STOKES' typdis ;
  314. $tgauche = 'MODELISER' _tgauche 'NAVIER_STOKES' typdis ;
  315. $tdroite = 'MODELISER' _tdroite 'NAVIER_STOKES' typdis ;
  316. 'ELIMINATION' 1.D-6 ('ET' ('DOMA' $tmt 'SOMMET')
  317. ('EXTRAIRE' uexact 'MAIL')) ;
  318. uexact = 'KCHT' $tmt 'VECT' 'SOMMET' uexact ;
  319. 'ELIMINATION' 1.D-6 ('ET' ('DOMA' $tmt 'CENTREP1')
  320. ('EXTRAIRE' pexact 'MAIL')) ;
  321. pexact = 'KCHT' $tmt 'SCAL' 'CENTREP1' pexact ;
  322. mtmt = 'DOMA' $tmt 'MAILLAGE' ;
  323. cmtmt = 'CONTOUR' mtmt ;
  324. *
  325. * Mise en place du calcul numérique
  326. *
  327. * équations de quantité de mouvement
  328. *
  329. rv = 'EQEX' $tmt 'ITMA' 1
  330. 'OPTI' 'EF' 'IMPL' 'CENTREE' 'CENTREP1'
  331. 'ZONE' $tmt 'OPER' 'KBBT' (+1.D0) 'INCO' 'UN' 'PN'
  332. 'ZONE' $tmt 'OPER' 'NS' 1. 'UN' (+0.5D0) 'INCO' 'UN' ;
  333. *
  334. * conditions aux limites
  335. *
  336. * vitesse
  337. rv = 'EQEX' rv
  338. 'CLIM' 'UN' 'UIMP' (_thaut 'ET' _tbas) 0.
  339. 'UN' 'VIMP' (_thaut 'ET' _tbas) 0.
  340. ;
  341. * pression
  342. rv = 'EQEX' rv
  343. 'OPTI' 'EF' 'IMPL' 'CENTREE' 'CENTREP1'
  344. 'ZONE' $tgauche 'OPER' 'TOIM' cgau 'INCO' 'UN'
  345. 'ZONE' $tdroite 'OPER' 'TOIM' cdro 'INCO' 'UN'
  346. ;
  347. * ________________
  348. *
  349. * INITIALISATION
  350. * ________________
  351. *
  352. rv . 'INCO' = 'TABLE' 'INCO' ;
  353. rv . 'INCO' . 'UN' = 'KCHT' $tmt 'VECT' 'SOMMET' (0. 0.) ;
  354. rv . 'INCO' . 'PN' = 'KCHT' $tmt 'SCAL' 'CENTREP1' 0. ;
  355.  
  356. rv . 'NITER' = 5 ;
  357. rv . 'METHINV' . 'TYPINV' = ktypi ;
  358. rv . 'METHINV' . 'PRECOND'= 3 ;
  359. * Les fréquences de recalcul du préconditionneur, respectivement
  360. * en fonction du pas de temps et des itérations dans la boucle pour
  361. * résoudre les non-linéarités.
  362. rv . 'METHINV' . 'FCPRECT' = 1 ;
  363. rv . 'METHINV' . 'FCPRECI' = 1000 ;
  364. rv . 'METHINV' . 'IMPINV' = 0 ;
  365. rv . 'METHINV' . 'NITMAX' = 1000 ;
  366. rv . 'METHINV' . 'GMRESTRT' = 50 ;
  367. rv . 'METHINV' . 'RESID' = 1.D-12 ;
  368. * __________
  369. *
  370. * CALCUL
  371. * __________
  372. *
  373. EXEC rv ;
  374. *
  375. * Résultats
  376. *
  377. 'SI' graph ;
  378. *
  379. * solutions exactes
  380. *
  381. un = uexact ;
  382. uref = (1.D0 * ('MINIMUM' lesdens));
  383. ung1 = 'VECTEUR' un uref 'UX' 'UY' 'JAUNE' ;
  384. titv = 'CHAINE' 'Vitesse UX UY exactes' titglob ;
  385. 'TRACER' ung1 cmtmt 'TITRE' titv ;
  386. pnn = 'ELNO' $tmt pexact 'CENTREP1' ;
  387. titp = 'CHAINE' 'Pression exacte ELNOisée' titglob ;
  388. 'TRACER' pnn mtmt cmtmt nbisov 'TITRE' titp ;
  389. *
  390. * graphe de convergence de la méthode itérative (dernière
  391. * itération)
  392. *
  393. 'SI' ('NEG' ktypi 1) ;
  394. conver = (rv . 'METHINV' . 'CONVINV') ;
  395. dimcon = 'DIME' conver ;
  396. labs = 'PROG' 1.D0 'PAS' 1.D0 dimcon ;
  397. lord = ('LOG' conver) '/' ('LOG' 10.D0) ;
  398. evtot = 'EVOL' 'MANU' 'ITER' labs 'RESID' lord ;
  399. titev = 'CHAINE' 'Historique de convergence' titglob ;
  400. 'DESSIN' evtot 'TITR' titev 'LEGE' ;
  401. 'FINSI' ;
  402. *
  403. * solutions calculées
  404. *
  405. un = rv . 'INCO' . 'UN' ;
  406. ung1 = 'VECTEUR' un uref 'UX' 'UY' 'JAUNE' ;
  407. titv = 'CHAINE' 'Vitesse UX UY calculées' titglob ;
  408. 'TRACER' ung1 cmtmt 'TITRE' titv ;
  409. pnn = 'ELNO' $tmt (rv . 'INCO' . 'PN') 'CENTREP1' ;
  410. pnn = 'KCHT' $tmt 'SCAL' 'SOMMET' pnn ;
  411. titp = 'CHAINE' 'Pression ELNOisée' titglob ;
  412. 'TRACER' pnn mtmt cmtmt nbisov 'TITRE' titp ;
  413. 'FINSI' ;
  414. *
  415. * Calcul des erreurs par rapport à la solution analytique
  416. *
  417. echloc = (('MESURE' mt) '/' ('NBEL' mt)) ** 0.5 ;
  418. un = (rv . 'INCO' . 'UN') ;
  419. pn = (rv . 'INCO' . 'PN') ;
  420. errutot = CALCERR2 un uexact ;
  421. errl2p = CALCERR2 pn pexact ;
  422. 'MESSAGE' ('CHAINE' 'errutot=' errutot) ;
  423. 'MESSAGE' ('CHAINE' 'errl2p=' errl2p) ;
  424. 'FINPROC' echloc errutot errl2p ;
  425. * Fin de la procédure de calcul *
  426. *___________________________________________________________*
  427. *
  428. *-----------------------------------------------------------
  429. * Paramètres du calcul
  430. *
  431. * Types d'inversion dans ltypi :
  432. * 1 = Crout (direct) ; 3 = BiCGSTAB ; 5 = GMRES
  433. * lraff : nb raffinement du maillage (à chaque fois, on découpe
  434. * les éléments en quatre).
  435. * ldiscr : type d'éléments à tester.
  436. 'SI' complet ;
  437. lraff = 'LECT' 0 PAS 1 2 ;
  438. ldiscr = 'MOTS' 'QUAF' 'MACR' ;
  439. ltypi = 'LECT' 3 ;
  440. 'SINON' ;
  441. lraff = 'LECT' 0 PAS 1 2 ;
  442. *!
  443. *! Ici, on devra mettre aussi 'MACR' quand ca marchera
  444. *!
  445. ldiscr = 'MOTS' 'QUAF' ;
  446. ltypi = 'LECT' 1 3 5 ;
  447. 'FINSI' ;
  448. *
  449. *-----------------------------------------------------------*
  450. * Boucles sur les différents paramètres du calcul *
  451. ordok = VRAI ;
  452. precok = VRAI ;
  453. *
  454. * Boucle sur les discrétisations
  455. *
  456. 'REPETER' idiscr ('DIME' ldiscr) ;
  457. discr = 'EXTRAIRE' ldiscr &idiscr ;
  458. * Ordres de convergence des éléments
  459. 'SI' ('EGA' discr 'QUAF') ;
  460. orduth = 3 ;
  461. ordpth = 2 ;
  462. 'SI' ('EGA' complet FAUX) ;
  463. errumax = 5.D-5 ;
  464. errpmax = 6.D-4 ;
  465. 'FINSI' ;
  466. 'SINON' ;
  467. 'SI' ('EGA' discr 'MACR') ;
  468. orduth = 2 ;
  469. ordpth = 1 ;
  470. 'SI' ('EGA' complet FAUX) ;
  471. 'MESSAGE' 'Définissez les erreurs admissibles pour MACR' ;
  472. 'ERREUR' 5 ;
  473. 'FINSI' ;
  474. 'SINON' ;
  475. 'MESSAGE' ('CHAINE' 'Type d_éléments non reconnus : ' discr ) ;
  476. 'ERREUR' 5 ;
  477. 'FINSI' ;
  478. 'FINSI' ;
  479. *
  480. * Boucle sur les méthodes de résolution
  481. *
  482. 'REPETER' itypi ('DIME' ltypi) ;
  483. typi = 'EXTRAIRE' ltypi &itypi ;
  484. *
  485. * Boucle sur les raffinements
  486. *
  487. lh = 'PROG' ;
  488. lerl2u = 'PROG' ;
  489. lerl2p = 'PROG' ;
  490. 'REPETER' iraff ('DIME' lraff) ;
  491. raff = 'EXTRAIRE' lraff &iraff ;
  492. h erru errp = CALCUL typi raff discr ;
  493. lh = 'ET' lh ('PROG' h ) ;
  494. lerl2u = 'ET' lerl2u ('PROG' erru) ;
  495. lerl2p = 'ET' lerl2p ('PROG' errp) ;
  496. 'FIN' iraff ;
  497. lh = ('LOG' lh) '/' ('LOG' 10.D0) ;
  498. lerl2u = ('LOG' lerl2u) '/' ('LOG' 10.D0) ;
  499. lerl2p = ('LOG' lerl2p) '/' ('LOG' 10.D0) ;
  500. ul2 = 'EVOL' 'MANU' 'Long. carac.' lh
  501. 'Err Vitesse' lerl2u ;
  502. pl2 = 'EVOL' 'MANU' 'Long. carac' lh
  503. 'Err Pression' lerl2p ;
  504. *
  505. * Recherche des ordres de convergence
  506. *
  507. cpu ulipol = @POMI ul2 1 'IDEM' ;
  508. ordu = cpu . 1 ;
  509. cpp plipol = @POMI pl2 1 'IDEM' ;
  510. ordp = cpp . 1 ;
  511. 'MESSAGE' ('CHAINE' 'ordre vitesse=' ordu
  512. ' ; ordre pression=' ordp) ;
  513. *
  514. * Tests des ordres de convergence
  515. * Tests des erreurs L2 sur le maillage le plus fin dans le
  516. * cas complet=faux
  517. *
  518. ordok = 'ET' ordok ('EGA' ('ENTIER' ('+' ordu 0.5)) orduth) ;
  519. ordok = 'ET' ordok ('EGA' ('ENTIER' ('+' ordp 0.5)) ordpth) ;
  520. 'SI' ('EGA' complet FAUX) ;
  521. precok = 'ET' precok ('<' erru errumax) ;
  522. precok = 'ET' precok ('<' errp errpmax) ;
  523. 'FINSI' ;
  524. 'SI' ('NON' ordok) ;
  525. 'MESSAGE' 'On n_obtient pas un ordre de convergence correct' ;
  526. 'ERREUR' 5 ;
  527. 'FINSI' ;
  528. 'SI' ('NON' precok) ;
  529. 'MESSAGE' 'On n_obtient pas une précision correcte' ;
  530. 'ERREUR' 5 ;
  531. 'FINSI' ;
  532. *
  533. * Tracés graphiques
  534. *
  535. 'SI' graph ;
  536. tableg = 'TABLE' ;
  537. tableg . 'TITRE' = 'TABLE' ;
  538. tableg . 1 = 'MARQ CROI NOLI' ;
  539. tableg . 'TITRE' . 1 = 'Erreur L2 calculée' ;
  540. tableg . 2 = 'TIRR' ;
  541. tableg . 'TITRE' . 2 = 'Erreur L2 interpolée' ;
  542. titdesu= 'CHAINE' 'Err L2 Vitesse, ordre=' ordu ' ; '
  543. 'discr=' discr ;
  544. 'DESSIN' (ul2 'ET' ulipol) 'TITRE' titdesu tableg ;
  545. titdesp= 'CHAINE' 'Err L2 Pression, ordre=' ordp ' ; '
  546. 'discr=' discr ;
  547. 'DESSIN' (pl2 'ET' plipol) 'TITRE' titdesp tableg ;
  548. 'FINSI' ;
  549. 'FIN' itypi ;
  550. 'FIN' idiscr ;
  551. 'SI' interact ;
  552. 'OPTION' 'ECHO' 1 ;
  553. 'OPTION' 'DONN' 5 ;
  554. 'FINSI' ;
  555. 'MESSAGE' 'Tout s_est bien passé' ;
  556. 'FIN' ;
  557.  
  558.  
  559.  
  560.  
  561.  
  562.  
  563.  
  564.  
  565.  
  566.  
  567.  
  568.  
  569.  
  570.  
  571.  
  572.  
  573.  
  574.  
  575.  
  576.  

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