Télécharger sinebump.dgibi

Retour à la liste

Numérotation des lignes :

  1. * fichier : sinebump.dgibi
  2. *********************************************************************
  3. * *
  4. * SINE-SHAPED BUMP *
  5. * *
  6. * CALCUL DE L'ECOULEMENT SUBSONIQUE ISENTROPIQUE STATIONNAIRE DANS *
  7. * UN CANAL *
  8. * *
  9. * BECCANTNI/PAILLERE, *
  10. * UPWIND FLUX SPLITTING SCHEMES FOR THE MULTIDIMENSIONAL EULER *
  11. * EQUATIONS, 1998 *
  12. * *
  13. * BECCANTINI A., DRN/DMT/SEMT/TTMF, NOVEMBRE 98 *
  14. * Modif, 10/07/01, syntaxe de PENT changée *
  15. *********************************************************************
  16.  
  17. 'OPTION' 'ECHO' 0 ;
  18. 'OPTION' 'DIME' 2 ;
  19. 'OPTION' 'ELEM' 'TRI3' ;
  20. 'OPTION' 'TRAC' 'PSC' ;
  21.  
  22. 'MESSAGE' 'A mettre a jours' ;
  23. 'FIN' ;
  24.  
  25. COMPLET = FAUX ;
  26. GRAPH = FAUX ;
  27.  
  28. NY = 4 ;
  29. NX1 = NY ;
  30. NX2 = 2 '*' NX1 ;
  31. NX3 = NX1 ;
  32. NX = NX1 '+' NX2 '+' NX3 ;
  33. DX = 4.0 '/' NX ;
  34.  
  35. CFL = 1.0D0 ;
  36. TFIN = 10.0 ;
  37.  
  38.  
  39. 'SI' COMPLET ;
  40. NRAFF = 3 ;
  41. 'SINON' ;
  42. NRAFF = 0 ;
  43. 'FINSI' ;
  44.  
  45.  
  46. *
  47. *** Methodes possibles :
  48. *
  49. * 'VANLEER'
  50. * 'VLH '
  51. * 'HUSVL '
  52. * 'HUSVLH '
  53. * 'GODUNOV'
  54.  
  55.  
  56.  
  57. METO = 'HUSVL ' ;
  58.  
  59. ***************************************************************
  60. ***** PROCEDURE POUR CALCULER LES VARIABLES CONSERVATIVES *****
  61. ***************************************************************
  62.  
  63. 'DEBPROC' CONS ;
  64. 'ARGUMENT' RN*'CHPOINT' VN* 'CHPOINT' PN*'CHPOINT' GAMN*'CHPOINT'
  65. $MODE*'MMODEL' ;
  66.  
  67. RN = 'KCHT' $MODE 'SCAL' 'CENTRE' RN ;
  68. PN = 'KCHT' $MODE 'SCAL' 'CENTRE' PN ;
  69. VN = 'KCHT' $MODE 'VECT' 'CENTRE' VN ;
  70. GAMN = 'KCHT' $MODE 'SCAL' 'CENTRE' GAMN ;
  71.  
  72. RVN = 'KOPS' RN '*' VN ;
  73. lcg=mots 'UX' 'UY';
  74. RECIN = 'KOPS' 0.5 '*' (RVN lcg 'PSCA' VN lcg) ;
  75. RETER = 'KOPS' PN '/' ('KOPS' GAMN '-' 1) ;
  76. REN = 'KOPS' RETER '+' RECIN ;
  77.  
  78. 'RESPRO' RVN REN ;
  79. 'FINPROC' ;
  80.  
  81. ****************************************************************
  82. ***** PROCEDURE POUR CALCULER LES VARIABLES ENTROPIE, MACH *****
  83. ****************************************************************
  84.  
  85. 'DEBPROC' EMMA ;
  86. 'ARGUMENT' RN*'CHPOINT' VN* 'CHPOINT' PN*'CHPOINT' GAMN*'CHPOINT'
  87. $MODE*'MMODEL' ;
  88.  
  89. RN = 'KCHT' $MODE 'SCAL' 'CENTRE' RN ;
  90. PN = 'KCHT' $MODE 'SCAL' 'CENTRE' PN ;
  91. VN = 'KCHT' $MODE 'VECT' 'CENTRE' VN ;
  92. GAMN = 'KCHT' $MODE 'SCAL' 'CENTRE' GAMN ;
  93.  
  94. lcg=mots 'UX' 'UY';
  95. MV = 'KOPS' (VN lcg 'PSCA' VN lcg) '**' 0.5 ;
  96. ASON = 'KOPS' ('KOPS' ('KOPS' PN '/' RN) '*' GAMN) '**' 0.5 ;
  97. MACHN = 'KOPS' MV '/' ASON ;
  98.  
  99. AN = 'KOPS' PN '/' ('KOPS' RN '**' GAMN) ;
  100.  
  101. 'RESPRO' MACHN AN ;
  102. 'FINPROC' ;
  103.  
  104.  
  105. *******************************************************
  106. ***** PROCEDURE POUR TESTER CONVERGENCE *****
  107. *******************************************************
  108.  
  109. 'DEBPROC' CALCUL ;
  110. 'ARGUMENT' RVX*'TABLE' ;
  111.  
  112. RV = RVX . 'EQEX' ;
  113. KINCO = RV . 'INCO' ;
  114. *KDOMA = RV . 'DOMAINE' ;
  115. KDOMA = RV . 'MODTOT' ;
  116.  
  117. RNi = 'COPIER' (RV . 'INCO' . 'RNI') ;
  118. RN0 = 'COPIER' (RV . 'INCO' . 'RN0') ;
  119.  
  120. *
  121. *** RN0 = solution à t = tn_M1;
  122. *
  123.  
  124. TPS = RV . 'PASDETPS' . 'TPS' ;
  125. DD = RV . 'PASDETPS' . 'NUPASDT' ;
  126. NINTE = 10 ;
  127. NN = DD '/' NINTE ;
  128. LO = ( DD '-' (NINTE '*' NN)) 'EGA' 0 ;
  129.  
  130. 'SI' ( LO ) ;
  131.  
  132. RNi = 'KCHT' KDOMA 'SCAL' 'CENTRE' RNi ;
  133. RN0 = 'KCHT' KDOMA 'SCAL' 'CENTRE' RN0 ;
  134. ERR = 'KOPS' RNi '-' RN0 ;
  135. ERR = 'ABS' ERR ;
  136. ERR = 'KOPS' ERR '/' RN0;
  137. KINCO . 'ERR' = ERR ;
  138. ELI = 'MAXIMUM' ERR 'ABS' ;
  139. ELI = ('LOG' (ELI + 1.0E-20))/('LOG' 10.) ;
  140. 'MESSAGE' ('CHAINE' 'ITER ' DD
  141. ' TPS ' TPS ' ERREUR LINF ' ELI );
  142. IT = 'PROG' DD ;
  143. ER = 'PROG' ELI ;
  144. RV . 'INCO' . 'IT' = (RV . 'INCO' . 'IT') 'ET' IT ;
  145. RV . 'INCO' . 'ER' = (RV . 'INCO' . 'ER') 'ET' ER ;
  146. 'FINSI' ;
  147.  
  148. RV . 'INCO' . 'RN0' = 'COPIER' RNi;
  149.  
  150. 'FINPROC' ;
  151.  
  152. *******************************************************
  153. ***** PROCEDURE POUR CALCULER L2 entropie *****
  154. *******************************************************
  155.  
  156. 'DEBPROC' CALCUL2 ;
  157. 'ARGUMENT' RVX*'TABLE' ;
  158.  
  159. RV = RVX . 'EQEX' ;
  160. KINCO = RV . 'INCO' ;
  161. *KDOMA = RV . 'DOMAINE' ;
  162. KDOMA = RV . 'MODTOT' ;
  163.  
  164. RN = KINCO . 'RNI' ;
  165. GN = KINCO . 'GNI' ;
  166. EN = KINCO . 'ENI' ;
  167. GAMMA = KINCO . 'GAMMA' ;
  168.  
  169. VN PN = 'PRIM' 'PERFMONO' RN GN EN GAMMA ;
  170. MACHN AN = EMMA RN VN PN GAMMA KDOMA ;
  171.  
  172. SNi = AN ;
  173. SN0 = 'COPIER' (RV . 'INCO' . 'SN0') ;
  174.  
  175. *
  176.  
  177. DD = RV . 'PASDETPS' . 'NUPASDT' ;
  178. NINTE = 10 ;
  179. NN = DD '/' NINTE ;
  180. LO = ( DD '-' (NINTE '*' NN)) 'EGA' 0 ;
  181.  
  182. 'SI' ( LO ) ;
  183.  
  184. SNi = 'KCHT' KDOMA 'SCAL' 'CENTRE' SNi ;
  185. SN0 = 'KCHT' KDOMA 'SCAL' 'CENTRE' SN0 ;
  186. ERR = 'KOPS' ('KOPS' SNi '-' SN0) '/' SN0 ;
  187. ERR = 'KOPS' ERR '**' 2 ;
  188. UNO = 'KCHT' KDOMA 'SCAL' 'CENTRE' 1.0 ;
  189. CHVOLU = 'DOMA' KDOMA 'XXVOLUM' ;
  190. TOTVOL = 'XTY' UNO CHVOLU ('MOTS' 'SCAL') ('MOTS' 'SCAL') ;
  191. ERRVOL = 'XTY' ERR CHVOLU ('MOTS' 'SCAL') ('MOTS' 'SCAL') ;
  192. EL2 = (ERRVOL '/' TOTVOL) '**' 0.5 ;
  193. 'MESSAGE' ('CHAINE'
  194. 'ITER ' DD ' ERREUR ENTROPIE L2 ' EL2 );
  195. IT = 'PROG' DD ;
  196. ER = 'PROG' EL2 ;
  197. RV . 'INCO' . 'IT1' = (RV . 'INCO' . 'IT1') 'ET' IT ;
  198. RV . 'INCO' . 'ER1' = (RV . 'INCO' . 'ER1') 'ET' ER ;
  199. 'FINSI' ;
  200.  
  201. RV . 'INCO' . 'SN0' = 'COPIER' SNi;
  202.  
  203. 'FINPROC' ;
  204.  
  205.  
  206. *****************************************************
  207. *****************************************************
  208. ** PROCEDURE EXEX POUR FORMULATION VF COMPRESSIBLE **
  209. ** CAS MONOESPECE **
  210. *****************************************************
  211. *****************************************************
  212.  
  213. 'DEBPROC' EXEX ;
  214. 'ARGUMENT' RV*TABLE ;
  215.  
  216. *******************************************
  217. **** Recherche de RV . *KONV . IDEUL ****
  218. *******************************************
  219.  
  220. *
  221. **** Nom de la table RV . *'KONV' -> NOMT
  222. *
  223.  
  224. NBOP = 'DIME' (RV . 'LISTOPER' ) ;
  225.  
  226. 'REPETER' BL1 NBOP;
  227. MCEL = 'EXTRAIRE' &BL1 RV . 'LISTOPER';
  228. 'SI' ( 'EGA' MCEL 'KONV ');
  229. NOMT = 'MOT' ('TEXTE' ('CHAINE' &BL1 MCEL));
  230. 'QUITTER' BL1;
  231. 'FINSI' ;
  232. 'FIN' BL1;
  233.  
  234. IEUL = RV . NOMT . 'KOPT' . 'IDEUL';
  235. 'SI' ('NON' (IEUL 'EGA' 1));
  236. 'MESSAGE' 'EULER 1E ???';
  237. 'ERREUR' 21;
  238. 'FINSI' ;
  239.  
  240.  
  241. *
  242. **** CL
  243. *
  244.  
  245.  
  246. LOGLIM = RV . 'INCO' . 'CLIM' ;
  247.  
  248.  
  249. ******************************************
  250. **** Ordre en espace, ordre en temps ****
  251. ******************************************
  252.  
  253. ORD_ESP = RV . 'ORDREESP' ;
  254. ORD_TPS = RV . 'ORDRETPS' ;
  255.  
  256. 'MESSAGE' '--------------------------';
  257. 'MESSAGE' 'Ordre en Espace :' ord_esp;
  258. 'MESSAGE' 'Ordre en Temps :' ord_tps;
  259. 'MESSAGE' '--------------------------';
  260.  
  261. 'SI' ((ORD_ESP 'EGA' 1) 'ET' (ORD_TPS 'EGA' 2));
  262. 'MESSAGE' ;
  263. 'MESSAGE' (CHAINE 'Ordre en Espace doit etre 2');
  264. 'MESSAGE' (CHAINE 'On impose ça.');
  265. 'MESSAGE' ;
  266. RV . 'ORDREESP' = 2 ;
  267. 'MESSAGE' ;
  268. 'MESSAGE' '--------------------------';
  269. 'MESSAGE' 'Ordre en Espace :' ord_esp;
  270. 'MESSAGE' 'Ordre en Temps :' ord_tps;
  271. 'MESSAGE' '--------------------------';
  272. 'FINSI' ;
  273.  
  274.  
  275. ******************************
  276. **** La table 'PASDETPS' ****
  277. ******************************
  278.  
  279. TPSI = RV . 'TPSI' ;
  280. TFIN = RV . 'TFINAL';
  281. RV . 'PASDETPS' . 'TPS' = TPSI;
  282.  
  283. *
  284. **** DELTAT-1 est un argument de PRET (prediction)
  285. * Donc on doit l'initialiser.
  286. *
  287.  
  288. RV . 'PASDETPS' . 'DELTAT-1' = 0.0D0;
  289. CFL = rv.'ALFA' ;
  290.  
  291.  
  292. *********************
  293. **** Les TABLES *****
  294. *********************
  295.  
  296. *
  297. **** RV . 'INCO'
  298. * RV . 'DOMAINE'
  299. * RV . 'KIZD'
  300. * RV . 'KIZG'
  301.  
  302.  
  303. *
  304. **** RV . 'INCO' -> KINCO
  305. *
  306.  
  307. KINCO = (RV . 'INCO') ;
  308.  
  309. *
  310. **** RV . 'DOMAINE' -> KDOMA
  311. *
  312.  
  313. *KDOMA = (RV . 'DOMAINE') ;
  314. KDOMA = (RV . 'MODTOT') ;
  315. KDOMA2 = (RV . 'DOMAINE') ;
  316.  
  317. *
  318. **** RV . 'KIZD' contient les volumes des ELTs
  319. *
  320.  
  321. 'SI' ('NON' ('EXISTE' RV 'KIZD')) ;
  322. 'KDIA' RV ;
  323. 'FINSI' ;
  324.  
  325. *
  326. ***** RV . 'KIZG' contient les flux aux interfaces.
  327. *
  328.  
  329. 'SI' ('NON' ('EXISTE' RV 'KIZG')) ;
  330. RV . 'KIZG' = 'TABLE' 'KIZG' ;
  331. 'FINSI' ;
  332.  
  333.  
  334. ******************************************
  335. **** Boucle Sur les Pas de Temps ****
  336. ******************************************
  337.  
  338. * No limiteurs des pentes
  339.  
  340.  
  341. ALSF0 = 'KCHT' (KDOMA2.'FRONT') 'SCAL' 'CENTRE' 'COMP' 'P1' 0.0;
  342. ALVF0 = 'KCHT' (KDOMA2.'FRONT') 'VECT'
  343. 'CENTRE' 'COMP' 'P1' 'P2' (0.0 0.0);
  344. ALR0 = 'KCHT' KDOMA 'SCAL' 'CENTRE' 'COMP' 'P1' 1.0 ;
  345. ALR0 = 'KCHT' KDOMA 'SCAL' 'CENTRE' 'COMP' 'P1'ALR0 ALSF0;
  346. ALV0 = 'KCHT' KDOMA 'VECT' 'CENTRE' 'COMP' 'P1' 'P2' (1.0 1.0) ;
  347. ALV0 = 'KCHT' KDOMA 'VECT' 'CENTRE' 'COMP' 'P1' 'P2' ALV0 ALVF0;
  348. ALP0 = 'KCHT' KDOMA 'SCAL' 'CENTRE' 'COMP' 'P1' 1.0 ;
  349. ALP0 = 'KCHT' KDOMA 'SCAL' 'CENTRE' 'COMP' 'P1' ALP0 ALSF0;
  350.  
  351. *
  352. **** Evaluation de coeff pour le calcule des pentes
  353. *
  354.  
  355. KINCO . 'V' KINCO . 'P' = 'PRIM' 'PERFMONO'
  356. (KINCO . 'RNI') (KINCO . 'GNI') (KINCO . 'ENI')
  357. (KINCO. 'GAMMA');
  358.  
  359.  
  360. GRADR CHPLIM COEFR = 'PENT' KDOMA 'CENTRE' 'EULESCAL' 'NOLIMITE'
  361. (KINCO . 'RNI');
  362.  
  363. GRADP CHPLIM COEFP = 'PENT' KDOMA 'CENTRE' 'EULESCAL' 'NOLIMITE'
  364. (KINCO . 'P');
  365.  
  366. GRADV CHPLIM COEFV = 'PENT' KDOMA 'CENTRE' 'EULEVECT' 'NOLIMITE'
  367. (KINCO . 'V');
  368.  
  369.  
  370. I = 0 ;
  371.  
  372. 'REPETER' BLOC1 (RV . 'ITMA') ;
  373.  
  374. I = I + 1 ;
  375.  
  376. *
  377. ***** Les variables primitives
  378. *
  379.  
  380. KINCO . 'V' KINCO . 'P' = 'PRIM' 'PERFMONO'
  381. (KINCO . 'RNI') (KINCO . 'GNI') (KINCO . 'ENI')
  382. (KINCO. 'GAMMA');
  383.  
  384. 'SI' (ORD_ESP 'EGA' 1) ;
  385.  
  386. ROF VITF PF GAMF =
  387. 'PRET' 'PERFMONO' ORD_ESP ORD_TPS KDOMA
  388. (KINCO . 'RNI') (KINCO . 'V')
  389. (KINCO . 'P') (KINCO . 'GAMMA');
  390.  
  391. 'SINON';
  392.  
  393. *
  394. ***** Ordre 2 en espace => calcul des pentes
  395. *
  396.  
  397. ALR = 'COPIER' ALR0 ;
  398. ALP = 'COPIER' ALP0 ;
  399. ALV = 'COPIER' ALV0 ;
  400.  
  401. GRADR CHPLIM = 'PENT' KDOMA 'CENTRE' 'EULESCAL' 'NOLIMITE'
  402. (KINCO . 'RNI') 'GRADGEO' COEFR ;
  403.  
  404. GRADP CHPLIM = 'PENT' KDOMA 'CENTRE' 'EULESCAL' 'NOLIMITE'
  405. (KINCO . 'P') 'GRADGEO' COEFP ;
  406.  
  407. GRADV CHPLIM = 'PENT' KDOMA 'CENTRE' 'EULEVECT' 'NOLIMITE'
  408. (KINCO . 'V') 'GRADGEO' COEFV ;
  409.  
  410. 'SI' (ORD_TPS 'EGA' 1);
  411.  
  412. ROF VITF PF GAMF = 'PRET' 'PERFMONO' ORD_ESP ORD_TPS
  413. KDOMA
  414. (KINCO . 'RNI') GRADR ALR
  415. (KINCO . 'V') GRADV ALV
  416. (KINCO . 'P') GRADP ALP
  417. (KINCO . 'GAMMA') ;
  418. 'SINON' ;
  419. *
  420. ********* Ordre 2 en temps
  421. *
  422. ROF VITF PF GAMF = 'PRET' 'PERFMONO' ORD_ESP ORD_TPS
  423. KDOMA
  424. (KINCO . 'RNI') GRADR ALR
  425. (KINCO . 'V') GRADV ALV
  426. (KINCO . 'P') GRADP ALP
  427. (KINCO . 'GAMMA')
  428. ((RV . 'PASDETPS' . 'DELTAT-1')/2.0);
  429. 'FINSI' ;
  430.  
  431. 'FINSI' ;
  432.  
  433. *
  434. *********** Creation de MCHAML de type 'FACEL' pour les
  435. * calcul de flux aux interfaces
  436.  
  437.  
  438.  
  439. KINCO . 'RNF' = ROF;
  440. KINCO . 'VITNF' = VITF;
  441. KINCO . 'PNF' = PF;
  442. KINCO . 'GAMMAF' = GAMF;
  443.  
  444.  
  445. *
  446. ********* Boucle sur les operateurs
  447. *
  448.  
  449. 'REPETER' BLOC2 NBOP ;
  450. NOMPER = 'EXTRAIRE' &BLOC2 (RV . 'LISTOPER');
  451. NOTABLE = 'MOT' ('TEXTE' ('CHAINE' &BLOC2 NOMPER) ) ;
  452. ('TEXTE' NOMPER) (RV . NOTABLE) ;
  453. 'FIN' BLOC2 ;
  454.  
  455.  
  456. ************************************************************
  457. ********* Mise a jour de la table RV . 'PASDETPS' ********
  458. ************************************************************
  459.  
  460. *
  461. **** Calcul du DT
  462. *
  463.  
  464. 'SI' ('EXISTE' RV 'DTI');
  465. *
  466. * Il faut comparer DTIMP, (DTCONV '*' CFL)
  467. *
  468. DTI = 'MINIMUM' ('PROG'
  469. (((RV . 'DTI') '/' CFL) (RV . 'PASDETPS' . 'DTCONV')
  470. ));
  471. 'SINON';
  472. DTI = (RV . 'PASDETPS' . 'DTCONV') ;
  473. 'FINSI';
  474.  
  475. RV . 'PASDETPS' . 'DELTAT' = DTI ;
  476.  
  477. *
  478. **** On verifie que (DTI * CFL) < (TFIN '-' TMPS)
  479. *
  480.  
  481. TMPS = RV . 'PASDETPS' . 'TPS';
  482. DTI0 = TFIN '-' TMPS;
  483. DTI0 = DTI0 '/' CFL;
  484.  
  485. 'SI' (DTI0 '&lt;EG' DTI);
  486. DTI = DTI0;
  487. RV . 'PASDETPS' . 'DELTAT' = DTI ;
  488. LOGQUIT = VRAI;
  489. 'SINON' ;
  490. LOGQUIT = FAUX ;
  491. 'FINSI' ;
  492.  
  493.  
  494. ***********************************************************
  495. ******** On avance au pas de temps suivant *********
  496. ***********************************************************
  497. *
  498. * N.B. Tn+1 = Tn + (CFL * RV . 'PASDETPS' . 'DELTAT');
  499. *
  500.  
  501. 'AVCT' RV CFL 'IMPR' (RV.'FIDT') ;
  502.  
  503. ******************************
  504. ******** On impose les CL ****
  505. ******************************
  506. *
  507. *
  508. 'SI' LOGLIM;
  509. PROLIM RV ;
  510. 'FINSI' ;
  511.  
  512.  
  513. 'SI' LOGQUIT;
  514. 'QUITTER' BLOC1;
  515. 'FINSI';
  516.  
  517.  
  518. 'FIN' BLOC1 ;
  519.  
  520.  
  521. ********************************************************
  522. ******** On detrui les choses qui ne servent plus ****
  523. ********************************************************
  524. *
  525. * Les variables primitives
  526. *
  527. 'DETR' ( KINCO . 'V');
  528. 'DETR' ( KINCO . 'P');
  529. 'OUBL' KINCO 'V';
  530. 'OUBL' KINCO 'P';
  531. *
  532. * Les MCHAML faces
  533. *
  534. 'DETR' ROF ;
  535. 'DETR' VITF ;
  536. 'DETR' PF ;
  537. 'DETR' GAMF;
  538. 'OUBL' KINCO 'RNF';
  539. 'OUBL' KINCO 'VITNF';
  540. 'OUBL' KINCO 'PNF';
  541. 'OUBL' KINCO 'GAMMAF';
  542. *
  543. * Les pentes 'ET' les limiteurs
  544. *
  545. 'SI' (ORD_ESP 'EGA' 2);
  546. 'DETR' GRADR;
  547. 'DETR' GRADP;
  548. 'DETR' GRADV;
  549. 'DETR' ALR;
  550. 'DETR' ALP;
  551. 'DETR' ALV;
  552. 'FINSI' ;
  553.  
  554.  
  555. ********************************************
  556. **** FIN de Boucle Sur les Pas de Temps ***
  557. ********************************************
  558.  
  559. 'FINPROC' ;
  560.  
  561. *****************************************************
  562. *****************************************************
  563. ** FIN PROCEDURE EXEX **
  564. *****************************************************
  565. *****************************************************
  566.  
  567. *****************************************************
  568. *****************************************************
  569. ***** PROCEDURE PROLIM *****
  570. *****************************************************
  571. *****************************************************
  572. *
  573. *
  574. **** Cas Euler Mono especes
  575. *
  576.  
  577.  
  578. 'DEBPROC' PROLIM ;
  579. 'ARGUMENT' RV*TABLE ;
  580.  
  581. *
  582. **** RV . 'DOMAINE' -> KDOMA
  583. *
  584.  
  585. KDOMA = RV . 'DOMAINE' ;
  586. KINCO = RV . 'INCO' ;
  587. KFRONT = RV . 'MODFRONT' ;
  588. * Domaine et connectivite entree subsonic
  589. KENTSUB = KDOMA . 'ENTSUB' ;
  590. KCONENT = KDOMA . 'CONENT' ;
  591. * Domaine et connectivite sortie subsonic
  592. KSORSUB = KDOMA . 'SORSUB' ;
  593. KCONSOR = KDOMA . 'CONSOR' ;
  594.  
  595. RNPE = 'KCHT' KFRONT 'SCAL' 'CENTRE'
  596. ('KPRO' (KINCO . 'RNI') KCONENT)
  597. ('KPRO' (KINCO . 'RNI') KCONSOR) ;
  598. GNPE = 'KCHT' KFRONT 'VECT' 'CENTRE'
  599. ('KPRO' (KINCO . 'GNI') KCONSOR)
  600. ('KPRO' (KINCO . 'GNI') KCONENT) ;
  601. ENPE = 'KCHT' KFRONT 'SCAL' 'CENTRE'
  602. ('KPRO' (KINCO . 'ENI') KCONENT)
  603. ('KPRO' (KINCO . 'ENI') KCONSOR) ;
  604. GAPE = 'KCHT' KFRONT 'SCAL' 'CENTRE'
  605. (KINCO . 'GAMMA' );
  606. SCAGAM = KINCO . 'SCALGAM' ;
  607. GM1PE = SCAGAM '-' 1.0 ;
  608. UNSGM1 = 1.0 '/' GM1PE ;
  609. VNPE PNPE = 'PRIM' 'PERFMONO' RNPE GNPE ENPE GAPE ;
  610.  
  611. *
  612. *
  613. * R1 = u '-' (2 '/' (gamma '-' 1)) '*' C (L2 = u '-' 'C)
  614. * R2 = s (L1 = u)
  615. * R3 = u '+' (2 '/' (gamma '-' 1)) '*' C (L2 = u '+' C)
  616. *
  617.  
  618. UXPE = 'EXCO' 'UX' VNPE ;
  619. UYPE = 'EXCO' 'UY' VNPE 'UY' ;
  620. ASPE = (2.0 '/' GM1PE) '*' ((GAPE '*' PNPE '/' RNPE ) '**' 0.5) ;
  621.  
  622. R1PE = UXPE '-' ASPE ;
  623. R2PE = PNPE '/' (RNPE '**' SCAGAM) ;
  624. R3PE = UXPE '+' ASPE ;
  625.  
  626. *
  627. **** Entree subsonic (u > 0)
  628. *
  629. * R2, R3, Uy sortent
  630. * R1 sort
  631. *
  632. **** Sortie subsonic (u > 0)
  633. *
  634. * R1 entre
  635. * R2, R3, Uy sortent
  636. *
  637. *
  638.  
  639. R1PE = 'KCHT' KFRONT 'SCAL' 'CENTRE' R1PE
  640. (KINCO . 'R1PE') ;
  641. R2PE = 'KCHT' KFRONT 'SCAL' 'CENTRE' R2PE
  642. (KINCO . 'R2PE') ;
  643. R3PE = 'KCHT' KFRONT 'SCAL' 'CENTRE' R3PE
  644. (KINCO . 'R3PE') ;
  645. UYPE = 'KCHT' KFRONT 'SCAL' 'CENTRE' 'COMP' 'UY'
  646. UYPE (KINCO . 'UYPE') ;
  647.  
  648. VNE = 'KCHT' KFRONT 'VECT' 'CENTRE'
  649. ('NOMC' 'UX' (0.5 '*' (R1PE '+' R3PE))) UYPE ;
  650. ASE2 = 'KCHT' KFRONT 'SCAL' 'CENTRE'
  651. (((GM1PE '/' 4.0) '*' (R3PE '-' R1PE)) '**' 2 ) ;
  652. RNE = 'KCHT' KFRONT 'SCAL' 'CENTRE'
  653. ((ASE2 '/' (R2PE '*' GAPE)) '**' UNSGM1) ;
  654. PNE = 'KOPS' ('KOPS' ASE2 '*' RNE) '/' GAPE ;
  655.  
  656. GNE ENE = CONS RNE VNE PNE GAPE KFRONT ;
  657.  
  658. KINCO . 'RNI' = 'KCHT' KDOMA 'SCAL' 'CENTRE'
  659. (KINCO . 'RNI') RNE ;
  660. KINCO . 'GNI' = 'KCHT' KDOMA 'VECT' 'CENTRE'
  661. (KINCO . 'GNI') GNE ;
  662. KINCO . 'ENI' = 'KCHT' KDOMA 'SCAL' 'CENTRE'
  663. (KINCO . 'ENI') ENE ;
  664.  
  665. 'FINPROC' ;
  666.  
  667.  
  668. *****************************************************
  669. *****************************************************
  670. ***** FIN PROCEDURE PROLIM *****
  671. *****************************************************
  672. *****************************************************
  673.  
  674.  
  675.  
  676.  
  677. ******************
  678. **** MAILLAGE ****
  679. ******************
  680.  
  681.  
  682. A0 = -2.0 0.0 ;
  683. A1 = -1.0 0.0 ;
  684. A2 = 1.0 0.0 ;
  685. A3 = 2.0 0.0 ;
  686. A4 = 2.0 1.0 ;
  687. A5 = -2.0 1.0 ;
  688.  
  689. *
  690. **** LIGB
  691. *
  692. LIGB1 = A0 'DROIT' NX1 A1 ;
  693.  
  694. * LIGB2 (On utilise un propriete de 'ET' ; 'SI' 'ET' change ?)
  695. xcel = ('COORDONNEE' 1 A1) '+' DX ;
  696. ycel = 0.1 '*' ( 1.0 '+' ('COS' (180 '*' xcel)));
  697. ACEL = xcel ycel ;
  698. LIGB2 = A1 'DROIT' 1 ACEL ;
  699. 'REPETER' BL1 (NX2 '-' 2) ;
  700. ACEL0 = ACEL ;
  701. xcel = xcel '+' DX ;
  702. ycel = 0.1 '*' ( 1.0 '+' ('COS' (180 '*' xcel)));
  703. ACEL = xcel ycel ;
  704. LIGB2 = LIGB2 'ET' (ACEL0 'DROIT' 1 ACEL) ;
  705. 'FIN' BL1;
  706. LIGB2 = LIGB2 'ET' (ACEL 'DROIT' 1 A2) ;
  707.  
  708.  
  709. LIGB3 = A2 'DROIT' NX3 A3 ;
  710.  
  711. LIGB = LIGB1 'ET' LIGB2 'ET' LIGB3 ;
  712.  
  713. *
  714. **** LIGH
  715. *
  716. LIGH = A4 'DROIT' NX A5 ;
  717.  
  718. *
  719. *** LIGG
  720. *
  721. LIGG = A5 'DROIT' NY A0 ;
  722.  
  723. *
  724. **** LIGD
  725. *
  726. LIGD = A3 'DROIT' NY A4 ;
  727.  
  728.  
  729. LIGT = LIGB 'ET' LIGD 'ET' LIGH 'ET' LIGG ;
  730.  
  731. *
  732. **** Le domaine
  733. *
  734. DOMINT = LIGT 'SURFACE' 'PLANE' ;
  735. DOMINT = 'DEDU' DOMINT LIGT LIGT 'REGU' ;
  736.  
  737. *
  738. **** Le rafinement
  739. *
  740. *
  741. **** C'est la Perle qui a fait ça quand il etait encore tres brillant!!!
  742. * (Avant la lettre de l'Arme!!!)
  743. * Mais il va bientot etre un vrai homme!!!
  744. *
  745.  
  746.  
  747. DOMCEL = DOMINT ;
  748.  
  749. 'SI' (nraff > 0) ;
  750. 'REPETER' bcl nraff ;
  751. mttemp = 'CHANGER' DOMCEL 'QUADRATIQUE' ;
  752. $mtt = 'DOMA' mttemp 'MACRO' ;
  753. DOMCEL = ($mtt . 'MAILLAGE') ;
  754. 'OUBLIER' mttemp ;
  755. 'OUBLIER' $mtt ;
  756. DOMINT = DOMCEL ;
  757. LIGCON = 'CONTOUR' DOMINT ;
  758.  
  759. *
  760. **** Je doit transler le points du BUMP (operateur FORME)
  761. *
  762. NLIGB2 = LIGCON 'ELEM' 'COMP' A1 A2 ;
  763. MAILB = 'CHANGER' NLIGB2 'POI1';
  764. MAILB0 = MAILB 'PLUS' (0.0D0 0.0D0);
  765.  
  766. * Listreels des deplacements
  767. NP = 'NBNO' MAILB ;
  768. LISTX = 'PROG' NP '*' 0.0 ;
  769. LISTY = 'PROG' NP '*' 0.0 ;
  770.  
  771. 'REPETER' BL1 ('NBNO' MAILB) ;
  772. xcel ycel0 = 'COORDONNEE' (MAILB 'POIN' &BL1);
  773. ycel = 0.1 '*' ( 1.0 '+' ('COS' (180 '*' xcel)));
  774. 'REMPLACER' LISTY &BL1 (ycel '-' ycel0) ;
  775. 'FIN' BL1 ;
  776.  
  777. DEPCH = 'MANUEL' 'CHPO' MAILB 2 'UX' LISTX 'UY' LISTY ;
  778. 'FORME' DEPCH ;
  779.  
  780. 'FIN' bcl ;
  781. DX = DX '/' (NRAFF '*' 2.0 ) ;
  782. 'SINON' ;
  783. LIGCON = 'CONTOUR' DOMINT ;
  784. 'FINSI' ;
  785. 'MENAGE' ;
  786.  
  787.  
  788.  
  789.  
  790. *
  791. **** Les C.L.
  792. *
  793.  
  794. NLIGD = LIGCON 'ELEM' 'COMP' A3 A4 ;
  795. NLIGG = LIGCON 'ELEM' 'COMP' A5 A0 ;
  796.  
  797. 'OPTION' 'ELEM' 'QUA4' ;
  798. FG = (NLIGG 'TRANSLATION' 1 ((-1.0 '*' DX) 0.0)) 'COULEUR' 'ROUGE' ;
  799. FD = (NLIGD 'TRANSLATION' 1 (DX 0.0)) 'COULEUR' 'ROUGE' ;
  800. FRONT = FG 'ET' FD ;
  801. *
  802. **** Domaine '+' CL
  803. *
  804.  
  805. DOMTOT = DOMINT 'ET' FG 'ET' FD ;
  806. 'ELIMINATION' DOMTOT 0.0001 ;
  807.  
  808. **********************
  809. *** OBJETS MODELES ***
  810. **********************
  811. MDOMTOT = 'CHANGER' DOMTOT 'QUAF' ;
  812. MDOMINT = 'CHANGER' DOMINT 'QUAF' ;
  813. MFG = 'CHANGER' FG 'QUAF' ;
  814. MFD = 'CHANGER' FD 'QUAF' ;
  815. MFRONT = 'CHANGER' FRONT 'QUAF' ;
  816. 'ELIM' (MDOMTOT 'ET' MDOMINT 'ET' MFG 'ET' MFD 'ET' MFRONT) 1.E-5;
  817. MDNS = 'EULER' ;
  818. $DOMTOT = 'MODE' MDOMTOT MDNS ;
  819. $DOMINT = 'MODE' MDOMINT MDNS ;
  820. $FG = 'MODE' MFG MDNS ;
  821. $FD = 'MODE' MFD MDNS ;
  822. $FRONT = 'MODE' MFRONT MDNS ;
  823.  
  824. 'SI' GRAPH ;
  825. MOD1 = 'MODELISER' ('DOMA' $DOMTOT 'MAILLAGE' ) 'THERMIQUE';
  826. 'FINSI' ;
  827. *
  828. **** On doit imposer un'entree subsonic (u>0 i.e. 3 condition a imposer) et une
  829. * sortie subsonic (1 condition)
  830. *
  831.  
  832. *
  833. **** Maillage entree subsonic
  834. *
  835.  
  836. POIN0 = NLIGG 'POIN' 1;
  837. X1 Y1 = 'COORDONNEE' POIN0;
  838. 'REPETER' BLLIM (('NBNO' NLIGG) '-' 1) ;
  839. X0 = X1 ;
  840. Y0 = Y1 ;
  841. POIN0 = NLIGG 'POIN' (&BLLIM '+' 1) ;
  842. X1 Y1 = 'COORDONNEE' POIN0 ;
  843. XFAC = (X0 '+' X1) '/' 2 ;
  844. YFAC = (Y0 '+' Y1) '/' 2 ;
  845. PFAC = ('DOMA' $DOMINT 'FACE') 'POIN' 'PROC' (XFAC YFAC);
  846. GEOFAC1 = ('DOMA' $DOMINT 'FACEL') 'ELEM' 'APPUYE'
  847. 'LARGEMENT' PFAC ;
  848. GEOFAC2 = ('DOMA' $FG 'FACEL') 'ELEM' 'APPUYE'
  849. 'LARGEMENT' PFAC ;
  850. *
  851. **** Tranformation en POI1
  852. *
  853. GEO1POI1 = 'CHANGER' 'POI1' GEOFAC1 ;
  854. PCEL11 = GEO1POI1 'POIN' 1 ;
  855. PCEL12 = GEO1POI1 'POIN' 2 ;
  856. GEO2POI1 = 'CHANGER' 'POI1' GEOFAC2 ;
  857. PCEL21 = GEO2POI1 'POIN' 1 ;
  858. PCEL22 = GEO2POI1 'POIN' 2 ;
  859. *
  860. **** Il faur verifier que PFAC = PCEL12 = PCEL22
  861. * ('NBEL' GEO1POI1) = ('NBEL' GEO2POI1) = 2
  862. *
  863. 'SI' ( ('NBEL' GEO1POI1) 'NEG' 2);
  864. 'MESSAGE' ;
  865. 'MESSAGE'
  866. 'Probleme dans la creation du domaine entree subsonique';
  867. 'MESSAGE' ;
  868. 'ERREUR' 21 ;
  869. 'FINSI' ;
  870. 'SI' ( ('NBEL' GEO2POI1) 'NEG' 2);
  871. 'MESSAGE' ;
  872. 'MESSAGE'
  873. 'Probleme dans la creation du domaine entree subsonique';
  874. 'MESSAGE' ;
  875. 'ERREUR' 21 ;
  876. 'FINSI' ;
  877. 'SI' ( PCEL12 'NEG' PFAC);
  878. 'MESSAGE' ;
  879. 'MESSAGE'
  880. 'Probleme dans la creation du domaine entree subsonique';
  881. 'MESSAGE' ;
  882. 'ERREUR' 21 ;
  883. 'FINSI' ;
  884. 'SI' ( PCEL22 'NEG' PFAC);
  885. 'MESSAGE' ;
  886. 'MESSAGE'
  887. 'Probleme dans la creation du domaine entree subsonique';
  888. 'MESSAGE' ;
  889. 'ERREUR' 21 ;
  890. 'FINSI' ;
  891. *
  892. *** Creation d'un maillage SEG2
  893. *
  894. 'SI' (&BLLIM 'EGA' 1);
  895. ENTSUB = 'MANUEL' 'SEG2' PCEL11 PCEL21 'COULEUR' 'VERT' ;
  896. 'SINON' ;
  897. ENTSUB = ENTSUB 'ET'
  898. ( 'MANUEL' 'SEG2' PCEL11 PCEL21 'COULEUR' 'VERT' );
  899. 'FINSI' ;
  900. 'FIN' BLLIM ;
  901.  
  902.  
  903.  
  904. *
  905. **** Maillage sortie subsonic
  906. *
  907.  
  908. POIN0 = NLIGD 'POIN' 1;
  909. X1 Y1 = 'COORDONNEE' POIN0;
  910. 'REPETER' BLLIM (('NBNO' NLIGD) '-' 1) ;
  911. X0 = X1 ;
  912. Y0 = Y1 ;
  913. POIN0 = NLIGD 'POIN' (&BLLIM '+' 1) ;
  914. X1 Y1 = 'COORDONNEE' POIN0 ;
  915. XFAC = (X0 '+' X1) '/' 2 ;
  916. YFAC = (Y0 '+' Y1) '/' 2 ;
  917. PFAC = ('DOMA' $DOMINT 'FACE') 'POIN' 'PROC' (XFAC YFAC);
  918. GEOFAC1 = ('DOMA' $DOMINT 'FACEL') 'ELEM' 'APPUYE'
  919. 'LARGEMENT' PFAC ;
  920. GEOFAC2 = ('DOMA' $FD 'FACEL') 'ELEM' 'APPUYE'
  921. 'LARGEMENT' PFAC ;
  922. *
  923. **** Tranformation en POI1
  924. *
  925. GEO1POI1 = 'CHANGER' 'POI1' GEOFAC1 ;
  926. PCEL11 = GEO1POI1 'POIN' 1 ;
  927. PCEL12 = GEO1POI1 'POIN' 2 ;
  928. GEO2POI1 = 'CHANGER' 'POI1' GEOFAC2 ;
  929. PCEL21 = GEO2POI1 'POIN' 1 ;
  930. PCEL22 = GEO2POI1 'POIN' 2 ;
  931. *
  932. **** Il faur verifier que PFAC = PCEL12 = PCEL22
  933. * ('NBEL' GEO1POI1) = ('NBEL' GEO2POI1) = 2
  934. *
  935. 'SI' ( ('NBEL' GEO1POI1) 'NEG' 2);
  936. 'MESSAGE' ;
  937. 'MESSAGE'
  938. 'Probleme dans la creation du domaine entree subsonique';
  939. 'MESSAGE' ;
  940. 'ERREUR' 21 ;
  941. 'FINSI' ;
  942. 'SI' ( ('NBEL' GEO2POI1) 'NEG' 2);
  943. 'MESSAGE' ;
  944. 'MESSAGE'
  945. 'Probleme dans la creation du domaine entree subsonique';
  946. 'MESSAGE' ;
  947. 'ERREUR' 21 ;
  948. 'FINSI' ;
  949. 'SI' ( PCEL12 'NEG' PFAC);
  950. 'MESSAGE' ;
  951. 'MESSAGE'
  952. 'Probleme dans la creation du domaine entree subsonique';
  953. 'MESSAGE' ;
  954. 'ERREUR' 21 ;
  955. 'FINSI' ;
  956. 'SI' ( PCEL22 'NEG' PFAC);
  957. 'MESSAGE' ;
  958. 'MESSAGE'
  959. 'Probleme dans la creation du domaine entree subsonique';
  960. 'MESSAGE' ;
  961. 'ERREUR' 21 ;
  962. 'FINSI' ;
  963. *
  964. *** Creation d'un maillage SEG2
  965. *
  966. 'SI' (&BLLIM 'EGA' 1);
  967. SORSUB = 'MANUEL' 'SEG2' PCEL11 PCEL21 'COULEUR' 'VERT' ;
  968. 'SINON' ;
  969. SORSUB = SORSUB 'ET'
  970. ( 'MANUEL' 'SEG2' PCEL11 PCEL21 'COULEUR' 'VERT' );
  971. 'FINSI' ;
  972. 'FIN' BLLIM ;
  973.  
  974. 'SI' GRAPH ;
  975. 'TRACER' (DOMTOT 'ET' SORSUB 'ET' ENTSUB)
  976. 'TITRE' 'Maillage et conditions aux bords' ;
  977. 'FINSI' ;
  978.  
  979. *
  980. *** C.L. 'ET' initiales
  981. *
  982.  
  983. RO_INF = 1.4D0 ;
  984. P_INF = 1.0D0 ;
  985. U_INF = 0.5D0 ;
  986. GAM = 1.4D0;
  987.  
  988.  
  989. RN0 = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' RO_INF ;
  990. VN0 = 'KCHT' $DOMTOT 'VECT' 'CENTRE' (U_INF 0.0) ;
  991. PN0 = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' P_INF ;
  992. GAMMA = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' GAM ;
  993.  
  994. GN0 REN0 = CONS RN0 VN0 PN0 GAMMA $DOMTOT ;
  995. MACHN0 AN0 = EMMA RN0 VN0 PN0 GAMMA $DOMTOT ;
  996.  
  997. ******************
  998. * TABLE EQUATION *
  999. ******************
  1000.  
  1001. NITER = 500000 ;
  1002.  
  1003. RV = 'EQEX' ('DOMA' $DOMTOT 'TABLE') 'ITMA' NITER 'ALFA' CFL
  1004. 'TPSI' 0. 'TFINAL' TFIN
  1005. *
  1006. ***** Options : VF = volumes finis
  1007. * CONS = conservative
  1008. * EULER = euler monoespece
  1009. * VANLEER = methode
  1010.  
  1011. 'OPTI' 'VF' 'CONS' 'EULER' METO
  1012.  
  1013. *
  1014. ***** Procedure 'CALCUL'
  1015. *
  1016. 'ZONE' $DOMTOT 'OPER' 'CALCUL'
  1017. *
  1018. ***** Procedure 'CALCUL2'
  1019. *
  1020. 'ZONE' $DOMTOT 'OPER' 'CALCUL2'
  1021. *
  1022. ***** Operateur 'KONV'
  1023. *
  1024. 'ZONE' $DOMTOT 'OPER' 'KONV'
  1025. *
  1026. ***** Arguments de 'KONV' (maximum 8 chatacters)
  1027. *
  1028. 'RNF' 'VITNF' 'PNF' 'GAMMAF' 'REPF'
  1029. *
  1030. ***** LIST des inconnues
  1031. *
  1032. 'INCO' 'RNI' 'GNI' 'ENI';
  1033.  
  1034. *
  1035. *** La table RV . INCO (de soustype INCO);
  1036. *
  1037.  
  1038. RV . 'INCO' = TABLE 'INCO' ;
  1039.  
  1040. *
  1041. *** Stokage des inconnues des equations d'Euler.
  1042. *
  1043.  
  1044. RV . 'INCO' . 'RNI' = 'COPIER' RN0 ;
  1045. RV . 'INCO' . 'GNI' = 'COPIER' GN0 ;
  1046. RV . 'INCO' . 'ENI' = 'COPIER' REN0 ;
  1047.  
  1048. *
  1049. *** GAMMA
  1050. *
  1051.  
  1052. RV . 'INCO' . 'GAMMA' = GAMMA ;
  1053. RV . 'INCO' . 'SCALGAM' = GAM ;
  1054.  
  1055.  
  1056. *
  1057. *** CONDITIONS LIMITS
  1058. *
  1059.  
  1060. * Invariants de Riemann 1D
  1061. *
  1062.  
  1063. A_INF = (((GAM '*' P_INF) '/' RO_INF) '**' 0.5) '*'
  1064. (2.0 '/' (GAM '-' 1.0)) ;
  1065. R1 = U_INF '-' A_INF ;
  1066. R2 = P_INF '/' (RO_INF '**' GAM) ;
  1067. R3 = U_INF '+' A_INF ;
  1068. RV . 'INCO' . 'CLIM' = VRAI;
  1069. ('DOMA' $DOMTOT 'TABLE') . 'FRONT' = 'DOMA' $FRONT 'TABLE' ;
  1070.  
  1071. * Entree subsonique
  1072.  
  1073. ('DOMA' $DOMTOT 'TABLE') . 'ENTSUB' = 'DOMA' $FG 'TABLE' ;
  1074. ('DOMA' $DOMTOT 'TABLE') . 'CONENT' = ENTSUB ;
  1075. RV . 'INCO' . 'R2PE' = 'KCHT' $FG 'SCAL' 'CENTRE' R2 ;
  1076. RV . 'INCO' . 'R3PE' = 'KCHT' $FG 'SCAL' 'CENTRE' R3 ;
  1077. RV . 'INCO' . 'UYPE' = 'KCHT' $FG 'SCAL' 'CENTRE' 'COMP' 'UY' 0.0 ;
  1078.  
  1079. *Sortie subsonique
  1080.  
  1081. ('DOMA' $DOMTOT 'TABLE') . 'SORSUB' = 'DOMA' $FD 'TABLE' ;
  1082. ('DOMA' $DOMTOT 'TABLE') . 'CONSOR' = SORSUB ;
  1083. RV . 'INCO' . 'R1PE' = 'KCHT' $FD 'SCAL' 'CENTRE' R1 ;
  1084.  
  1085. *
  1086. **** Ordre en espace
  1087. * Ordre en temps
  1088. *
  1089.  
  1090. IT = 1 ;
  1091. IE = 1 ;
  1092. RV . 'ORDRETPS' = IT ;
  1093. RV . 'ORDREESP' = IE ;
  1094.  
  1095. *
  1096. **** Convergence
  1097. *
  1098.  
  1099. RV . 'INCO' . 'IT' = 'PROG' ;
  1100. RV . 'INCO' . 'ER' = 'PROG' ;
  1101. RV . 'INCO' . 'IT1' = 'PROG' ;
  1102. RV . 'INCO' . 'ER1' = 'PROG' ;
  1103. RV . 'INCO' . 'RN0' = 'COPIER' RN0 ;
  1104. RV . 'INCO' . 'SN0' = 'KCHT' $DOMTOT 'SCAL' 'CENTRE'
  1105. (P_INF '/' (RO_INF '**' 1.4));
  1106.  
  1107. *???
  1108. RV . 'MODTOT' = $DOMTOT ;
  1109. RV . 'MODFRONT' = $FRONT ;
  1110.  
  1111. *
  1112. *
  1113. ***********************************
  1114. *** Execution EXEX ***
  1115. ***********************************
  1116. *
  1117.  
  1118.  
  1119. 'MESSAGE' ;
  1120. 'MESSAGE' '**************************';
  1121. 'MESSAGE' ('CHAINE' 'METHODE : ' METO) ;
  1122. 'MESSAGE' '**************************';
  1123. 'MESSAGE' ;
  1124.  
  1125.  
  1126.  
  1127. RV . 'INCO' . 'V' = 'COPIER' VN0 ;
  1128. RV . 'INCO' . 'P' = 'COPIER' PN0 ;
  1129.  
  1130.  
  1131. 'TEMPS' ;
  1132. EXEX RV ;
  1133. 'TEMPS' ;
  1134.  
  1135.  
  1136.  
  1137. TFIN = RV . 'PASDETPS' . 'TPS' ;
  1138. ERR = RV . 'INCO' . 'ERR' ;
  1139. RN = RV . 'INCO' . 'RNI' ;
  1140. GN = RV . 'INCO' . 'GNI' ;
  1141. REN = RV . 'INCO' . 'ENI' ;
  1142.  
  1143. VN PN = 'PRIM' 'PERFMONO' RN GN REN (RV . 'INCO'. 'GAMMA');
  1144.  
  1145. MACHN AN = EMMA RN VN PN GAMMA $DOMTOT ;
  1146.  
  1147. EVOL4 = EVOL 'MANU' 'ITERATIONS' (RV.INCO.'IT') 'LOG|E|inf'
  1148. (RV.INCO.'ER') ;
  1149.  
  1150. AS = ((GAMMA '*' PN '/' RN) '**' 0.5) '*' ( 2.0 '/' (GAM '-' 1.0)) ;
  1151. UX = 'EXCO' 'UX' VN ;
  1152. UY = 'EXCO' 'UY' VN ;
  1153. R1 = UX '-' AS ;
  1154. R2 = AN ;
  1155. R3 = UX '+' AS ;
  1156.  
  1157. 'SI' GRAPH ;
  1158. *
  1159. *** Invariants de Riemann pour les equations 1D
  1160. *
  1161. CHM_R1 = 'KCHA' $DOMTOT 'CHAM' R1 ;
  1162. CHM_R2 = 'KCHA' $DOMTOT 'CHAM' R2 ;
  1163. CHM_R3 = 'KCHA' $DOMTOT 'CHAM' R3 ;
  1164. 'TRACER' MOD1 CHM_R1 'TITRE' 'R1';
  1165. 'TRACER' MOD1 CHM_R2 'TITRE' 'R2' ;
  1166. 'TRACER' MOD1 CHM_R3 'TITRE' 'R3' ;
  1167. *
  1168.  
  1169. CHM_RN = 'KCHA' $DOMTOT 'CHAM' RN ;
  1170. CHM_PN = 'KCHA' $DOMTOT 'CHAM' PN ;
  1171. CHM_AN = 'KCHA' $DOMTOT 'CHAM' AN ;
  1172. CHM_MN = 'KCHA' $DOMTOT 'CHAM' MACHN ;
  1173. CHM_UNY = 'KCHA' $DOMTOT 'CHAM' ('EXCO' 'UY' VN) ;
  1174. 'TRACER' MOD1 CHM_RN 'TITRE' 'Ro';
  1175. 'TRACER' MOD1 CHM_PN 'TITRE' 'P' ;
  1176. 'TRACER' MOD1 CHM_AN 'TITRE' 'Entropie' ;
  1177. 'TRACER' MOD1 CHM_MN 'TITRE' 'Mach' ;
  1178. 'TRACER' MOD1 CHM_UNY 'TITRE' 'UY' ;
  1179. 'DESSIN' EVOL4 'XBOR' 0. (1.5 '*' ('MAXI' (RV . 'INCO' .'IT') ))
  1180. 'TITRE' 'Erreur Linf sur ro' ;
  1181.  
  1182. 'FINSI' ;
  1183.  
  1184. *
  1185. **** Controle sur L_inf, L_2 ENTROPIE
  1186. *
  1187.  
  1188. AN_INF = 1.0D0 '/' (1.4 '**' 1.4) ;
  1189. SIGMA = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'SIGMA'
  1190. ('NOMC' 'SIGMA' ((AN '/' AN_INF) '-' 1.0)) ;
  1191.  
  1192. LINF_SIG = 'MAXIMUM' SIGMA 'ABS' ;
  1193.  
  1194. 'SI' COMPLET ;
  1195. 'SI' (LINF_SIG > 1.5D-2) ;
  1196. 'MESSAGE' ('CHAINE' METO ' ' 'Ordre ' IE ' ' IT);
  1197. 'MESSAGE' ('CHAINE' 'Norme infinite sur entropie ' LINF_SIG) ;
  1198. 'ERREUR' 5;
  1199. 'FINSI' ;
  1200. 'SINON' ;
  1201. 'SI' (LINF_SIG > 5D-2) ;
  1202. 'MESSAGE' ('CHAINE' METO ' ' 'Ordre ' IE ' ' IT);
  1203. 'MESSAGE' ('CHAINE' 'Norme infinite sur entropie ' LINF_SIG) ;
  1204. 'ERREUR' 5;
  1205. 'FINSI' ;
  1206. 'FINSI' ;
  1207.  
  1208. SIGMAI2 = 'KCHT' $DOMINT 'SCAL' 'CENTRE' 'COMP' 'SIGMA'
  1209. (SIGMA '**' 2);
  1210. 'DOMA' $DOMINT 'VOLUME' ;
  1211. CHPUN = 'KCHT' $DOMINT 'SCAL' 'CENTRE' 1.0 ;
  1212. M1 = 'MOTS' 'SCAL' ;
  1213. M2 = 'MOTS' 'SIGMA' ;
  1214. TOTVOL = 'XTY' ('DOMA' $DOMINT 'XXVOLUM') CHPUN M1 M1;
  1215. L2SIGMA = ('XTY' ('DOMA' $DOMINT 'XXVOLUM') SIGMAI2 M1 M2 )
  1216. '/' TOTVOL ;
  1217. L2SIGMA = L2SIGMA '**' 0.5 ;
  1218.  
  1219. 'SI' COMPLET ;
  1220. 'SI' (L2SIGMA > 5D-3) ;
  1221. 'MESSAGE' ('CHAINE' METO ' ' 'Ordre ' IE ' ' IT);
  1222. 'MESSAGE' ('CHAINE' 'Norme 2 sur entropie ' L2SIGMA) ;
  1223. 'ERREUR' 5;
  1224. 'FINSI' ;
  1225. 'SINON' ;
  1226. 'SI' (L2SIGMA > 1.2D-2) ;
  1227. 'MESSAGE' ('CHAINE' METO ' ' 'Ordre ' IE ' ' IT);
  1228. 'MESSAGE' ('CHAINE' 'Norme 2 sur entropie ' L2SIGMA) ;
  1229. 'ERREUR' 5;
  1230. 'FINSI' ;
  1231. 'FINSI' ;
  1232.  
  1233. ERR_CON = 'EXTRAIRE' (RV . 'INCO' .'ER') ('DIME' (RV . 'INCO' .'ER')) ;
  1234.  
  1235. 'SI' (ERR_CON > -3.5) ;
  1236. 'MESSAGE' 'Probleme de convergence' ;
  1237. 'ERREUR' 5 ;
  1238. 'FINSI' ;
  1239.  
  1240. 'SI' COMPLET ;
  1241. *
  1242. **** Calcul aux deuxieme ordre en espace et en temps
  1243. *
  1244.  
  1245. 'OUBLIER' RV;
  1246. MENAGE ;
  1247.  
  1248.  
  1249. ******************
  1250. * TABLE EQUATION *
  1251. ******************
  1252.  
  1253. NITER = 500000 ;
  1254.  
  1255. RV = 'EQEX' ('DOMA' $DOMTOT 'TABLE') 'ITMA' NITER 'ALFA' CFL
  1256. 'TPSI' 0. 'TFINAL' TFIN
  1257. *
  1258. ***** Options : VF = volumes finis
  1259. * CONS = conservative
  1260. * EULER = euler monoespece
  1261. * VANLEER = methode
  1262.  
  1263. 'OPTI' 'VF' 'CONS' 'EULER' METO
  1264.  
  1265. *
  1266. ***** Procedure 'CALCUL'
  1267. *
  1268. 'ZONE' $DOMTOT 'OPER' 'CALCUL'
  1269. *
  1270. ***** Procedure 'CALCUL2'
  1271. *
  1272. 'ZONE' $DOMTOT 'OPER' 'CALCUL2'
  1273. *
  1274. ***** Operateur 'KONV'
  1275. *
  1276. 'ZONE' $DOMTOT 'OPER' 'KONV'
  1277. *
  1278. ***** Arguments de 'KONV' (maximum 8 chatacters)
  1279. *
  1280. 'RNF' 'VITNF' 'PNF' 'GAMMAF' 'REPF'
  1281. *
  1282. ***** LIST des inconnues
  1283. *
  1284. 'INCO' 'RNI' 'GNI' 'ENI';
  1285.  
  1286. *
  1287. *** La table RV . INCO (de soustype INCO);
  1288. *
  1289.  
  1290. RV . 'INCO' = TABLE 'INCO' ;
  1291.  
  1292. *
  1293. *** Stokage des inconnues des equations d'Euler.
  1294. *
  1295.  
  1296. RV . 'INCO' . 'RNI' = 'COPIER' RN ;
  1297. RV . 'INCO' . 'GNI' = 'COPIER' GN ;
  1298. RV . 'INCO' . 'ENI' = 'COPIER' REN ;
  1299.  
  1300. *
  1301. *** GAMMA
  1302. *
  1303.  
  1304. RV . 'INCO' . 'GAMMA' = GAMMA ;
  1305. RV . 'INCO' . 'SCALGAM' = GAM ;
  1306.  
  1307.  
  1308. *
  1309. *** CONDITIONS LIMITS
  1310. *
  1311.  
  1312. * Invariants de Riemann 1D
  1313. *
  1314.  
  1315. A_INF = (((GAM '*' P_INF) '/' RO_INF) '**' 0.5) '*'
  1316. (2.0 '/' (GAM '-' 1.0)) ;
  1317. R1 = U_INF '-' A_INF ;
  1318. R2 = P_INF '/' (RO_INF '**' GAM) ;
  1319. R3 = U_INF '+' A_INF ;
  1320. RV . 'INCO' . 'CLIM' = VRAI;
  1321. $DOMTOT . 'FRONT' = $FRONT ;
  1322.  
  1323. * Entree subsonique
  1324.  
  1325. $DOMTOT . 'ENTSUB' = $FG ;
  1326. $DOMTOT . 'CONENT' = ENTSUB ;
  1327. RV . 'INCO' . 'R2PE' = 'KCHT' $FG 'SCAL' 'CENTRE' R2 ;
  1328. RV . 'INCO' . 'R3PE' = 'KCHT' $FG 'SCAL' 'CENTRE' R3 ;
  1329. RV . 'INCO' . 'UYPE' = 'KCHT' $FG 'SCAL' 'CENTRE' 'COMP' 'UY' 0.0 ;
  1330.  
  1331. *Sortie subsonique
  1332.  
  1333. $DOMTOT . 'SORSUB' = $FD ;
  1334. $DOMTOT . 'CONSOR' = SORSUB ;
  1335. RV . 'INCO' . 'R1PE' = 'KCHT' $FD 'SCAL' 'CENTRE' R1 ;
  1336.  
  1337. *
  1338. **** Ordre en espace
  1339. * Ordre en temps
  1340. *
  1341.  
  1342. IT = 2 ;
  1343. IE = 2 ;
  1344. RV . 'ORDRETPS' = IT ;
  1345. RV . 'ORDREESP' = IE ;
  1346.  
  1347. *
  1348. **** Convergence
  1349. *
  1350.  
  1351. RV . 'INCO' . 'IT' = 'PROG' ;
  1352. RV . 'INCO' . 'ER' = 'PROG' ;
  1353. RV . 'INCO' . 'IT1' = 'PROG' ;
  1354. RV . 'INCO' . 'ER1' = 'PROG' ;
  1355. RV . 'INCO' . 'RN0' = 'COPIER' RN0 ;
  1356. RV . 'INCO' . 'SN0' = 'KCHT' $DOMTOT 'SCAL' 'CENTRE'
  1357. (P_INF '/' (RO_INF '**' 1.4));
  1358.  
  1359.  
  1360. *
  1361. *
  1362. ***********************************
  1363. *** Execution EXEX ***
  1364. ***********************************
  1365. *
  1366.  
  1367.  
  1368. 'MESSAGE' ;
  1369. 'MESSAGE' '**************************';
  1370. 'MESSAGE' ('CHAINE' 'METHODE : ' METO) ;
  1371. 'MESSAGE' '**************************';
  1372. 'MESSAGE' ;
  1373.  
  1374.  
  1375.  
  1376. RV . 'INCO' . 'V' = 'COPIER' VN0 ;
  1377. RV . 'INCO' . 'P' = 'COPIER' PN0 ;
  1378.  
  1379.  
  1380. 'TEMPS' ;
  1381. EXEX RV ;
  1382. 'TEMPS' ;
  1383.  
  1384.  
  1385.  
  1386. TFIN = RV . 'PASDETPS' . 'TPS' ;
  1387. ERR = RV . 'INCO' . 'ERR' ;
  1388. RN = RV . 'INCO' . 'RNI' ;
  1389. GN = RV . 'INCO' . 'GNI' ;
  1390. REN = RV . 'INCO' . 'ENI' ;
  1391.  
  1392. VN PN = 'PRIM' 'PERFMONO' RN GN REN (RV . 'INCO'. 'GAMMA');
  1393.  
  1394. MACHN AN = EMMA RN VN PN GAMMA $DOMTOT ;
  1395.  
  1396. EVOL4 = EVOL 'MANU' 'ITERATIONS' (RV.INCO.'IT') 'LOG|E|inf'
  1397. (RV.INCO.'ER') ;
  1398.  
  1399. AS = ((GAMMA '*' PN '/' RN) '**' 0.5) '*' ( 2.0 '/' (GAM '-' 1.0)) ;
  1400. UX = 'EXCO' 'UX' VN ;
  1401. UY = 'EXCO' 'UY' VN ;
  1402. R1 = UX '-' AS ;
  1403. R2 = AN ;
  1404. R3 = UX '+' AS ;
  1405.  
  1406. 'SI' GRAPH ;
  1407. *
  1408. *** Invariants de Riemann pour les equations 1D
  1409. *
  1410. CHM_R1 = 'KCHA' $DOMTOT 'CHAM' R1 ;
  1411. CHM_R2 = 'KCHA' $DOMTOT 'CHAM' R2 ;
  1412. CHM_R3 = 'KCHA' $DOMTOT 'CHAM' R3 ;
  1413. 'TRACER' MOD1 CHM_R1 'TITRE' 'R1';
  1414. 'TRACER' MOD1 CHM_R2 'TITRE' 'R2' ;
  1415. 'TRACER' MOD1 CHM_R3 'TITRE' 'R3' ;
  1416. *
  1417.  
  1418. CHM_RN = 'KCHA' $DOMTOT 'CHAM' RN ;
  1419. CHM_PN = 'KCHA' $DOMTOT 'CHAM' PN ;
  1420. CHM_AN = 'KCHA' $DOMTOT 'CHAM' AN ;
  1421. CHM_MN = 'KCHA' $DOMTOT 'CHAM' MACHN ;
  1422. CHM_UNY = 'KCHA' $DOMTOT 'CHAM' ('EXCO' 'UY' VN) ;
  1423. 'TRACER' MOD1 CHM_RN 'TITRE' 'Ro';
  1424. 'TRACER' MOD1 CHM_PN 'TITRE' 'P' ;
  1425. 'TRACER' MOD1 CHM_AN 'TITRE' 'Entropie' ;
  1426. 'TRACER' MOD1 CHM_MN 'TITRE' 'Mach' ;
  1427. 'TRACER' MOD1 CHM_UNY 'TITRE' 'UY' ;
  1428. 'DESSIN' EVOL4 'XBOR' 0. (1.5 '*' ('MAXI' (RV . 'INCO' .'IT') ))
  1429. 'TITRE' 'Erreur Linf sur ro' ;
  1430.  
  1431. 'FINSI' ;
  1432.  
  1433. *
  1434. **** Controle sur L_inf, L_2 ENTROPIE
  1435. *
  1436.  
  1437. AN_INF = 1.0D0 '/' (1.4 '**' 1.4) ;
  1438. SIGMA = 'KCHT' $DOMTOT 'SCAL' 'CENTRE' 'COMP' 'SIGMA'
  1439. ('NOMC' 'SIGMA' ((AN '/' AN_INF) '-' 1.0)) ;
  1440.  
  1441. LINF_SIG = 'MAXIMUM' SIGMA 'ABS' ;
  1442. 'SI' COMPLET ;
  1443. 'SI' (LINF_SIG > 6D-4) ;
  1444. 'MESSAGE' ('CHAINE' METO ' ' 'Ordre ' IE ' ' IT);
  1445. 'MESSAGE' ('CHAINE' 'Norme infinite sur entropie ' LINF_SIG) ;
  1446. 'ERREUR' 5;
  1447. 'FINSI' ;
  1448. 'SINON' ;
  1449. 'SI' (LINF_SIG > 8D-3) ;
  1450. 'MESSAGE' ('CHAINE' METO ' ' 'Ordre ' IE ' ' IT);
  1451. 'MESSAGE' ('CHAINE' 'Norme infinite sur entropie ' LINF_SIG) ;
  1452. 'ERREUR' 5;
  1453. 'FINSI' ;
  1454. 'FINSI' ;
  1455.  
  1456. SIGMAI2 = 'KCHT' $DOMINT 'SCAL' 'CENTRE' 'COMP' 'SIGMA'
  1457. (SIGMA '**' 2);
  1458. 'DOMA' $DOMINT 'VOLUME' ;
  1459. CHPUN = 'KCHT' $DOMINT 'SCAL' 'CENTRE' 1.0 ;
  1460. M1 = 'MOTS' 'SCAL' ;
  1461. M2 = 'MOTS' 'SIGMA' ;
  1462. TOTVOL = 'XTY' ($DOMINT . 'XXVOLUM') CHPUN M1 M1;
  1463. L2SIGMA = ('XTY' ($DOMINT . 'XXVOLUM') SIGMAI2 M1 M2 ) '/' TOTVOL ;
  1464. L2SIGMA = L2SIGMA '**' 0.5 ;
  1465.  
  1466. 'SI' COMPLET ;
  1467. 'SI' (L2SIGMA > 6D-5) ;
  1468. 'MESSAGE' ('CHAINE' METO ' ' 'Ordre ' IE ' ' IT);
  1469. 'MESSAGE' ('CHAINE' 'Norme 2 sur entropie ' L2SIGMA) ;
  1470. 'ERREUR' 5;
  1471. 'FINSI' ;
  1472. 'SINON' ;
  1473. 'SI' (L2SIGMA > 5D-3) ;
  1474. 'MESSAGE' ('CHAINE' METO ' ' 'Ordre ' IE ' ' IT);
  1475. 'MESSAGE' ('CHAINE' 'Norme 2 sur entropie ' L2SIGMA) ;
  1476. 'ERREUR' 5;
  1477. 'FINSI' ;
  1478. 'FINSI' ;
  1479.  
  1480. 'FINSI' ;
  1481.  
  1482. 'FIN' ;
  1483.  
  1484.  
  1485.  
  1486.  
  1487.  
  1488.  
  1489.  
  1490.  
  1491.  
  1492.  
  1493.  
  1494.  
  1495.  
  1496.  
  1497.  
  1498.  
  1499.  
  1500.  
  1501.  

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