Télécharger la plaquette Cast3M


fr en

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.5D-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.  
  577.  
  578.  

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