Télécharger la plaquette Cast3M


fr en

Télécharger sinebum_fmm4.dgibi

Retour à la liste

Numérotation des lignes : 

  1. * fichier :  sinebum_fmm4.dgibi
  2. *************************************************************************
  3. *                                                                       *
  4. * SINE-SHAPED BUMP                                                      *
  5. *                                                                       *
  6. * CALCUL DE L'ECOULEMENT SUBSONIQUE ISENTROPIQUE STATIONNAIRE DANS      *
  7. * UN CANAL                                                              *
  8. *                                                                       *
  9. * Methode implicite sans matrice pour les equations d'Euler             *
  10. * (bas Mach)                                                            *
  11. *                                                                       *
  12. * BECCANTINI A., SFME/LTMF,  DEC 2003                                   *
  13. *                                                                       *
  14. *                                                                       *
  15. *************************************************************************
  16.  
  17. *************************************************************************
  18. *************************************************************************
  19. *********** SOLUTION OF THE EULER EQUATIONS  ****************************
  20. *********** FMM                              ****************************
  21. *************************************************************************
  22. *************************************************************************
  23. *
  24. *
  25. * RVX . 'MODEL'           : model object
  26. * RVX . 'RN0'             : density
  27. * RVX . 'GN0'             : qdm
  28. * RVX . 'RET0'            : total energy
  29. *
  30. * RVX . 'PGAS'            : gas model 
  31. * RVX . 'PGAS' . 'GAMN'   : gamma
  32. *
  33. * RVX . 'GRAVITY'         : gravity
  34. *
  35. * RVX . 'LISTCONS'        : name of the conservative variables
  36. * RVX . 'LISTPRIM'        : name of the primitive variables
  37. *
  38. * RVX . 'LISTERR'         : name of the error variables
  39. *
  40. * RVX . 'METHOD'          : numerical scheme
  41. * RVX . 'CUTOFF'          : cut off speed
  42. *
  43. * RVX . 'SPACEA'          : space accuracy
  44. * RVX . 'LIMITER'         : limiter type
  45. * RVX . 'TIMEA'           : time accuracy
  46. * 
  47. * RVX . 'T0'              : initial time
  48. * RVX . 'TFINAL'          : final time
  49. * RVX . 'DTPS' ('CFL')    : time step or CFL number
  50. *
  51. * RVX . 'DCFL'            : CFL number for dual time
  52. *
  53. *                          Error criteria for dual time loop:
  54. * RVX . 'NDTITER'         : number of iterations
  55. * RVX . 'EPSDT'           : relative error
  56. *
  57. * RVX . 'NJAC'            : Jacobi iterations
  58. *
  59.  'DEBPROC' PEUSM ;
  60.  'ARGUMENT' RVX*TABLE ;
  61. *
  62.  'SI' ('EXISTE' RVX 'RESULTS') ;
  63.     'MESS' 'Table RESULTS already exists' ;
  64.  'SINON' ;
  65.     RVX . 'RESULTS' = 'TABLE' ;
  66.     RVX . 'RESULTS' . 'TPS' = RVX . 'T0' ;
  67.     RVX . 'RESULTS' . 'RN'  = 'COPIER' (RVX . 'RN0') ;
  68.     RVX . 'RESULTS' . 'GN'  = 'COPIER' (RVX . 'GN0') ;
  69.     RVX . 'RESULTS' . 'RET' = 'COPIER' (RVX . 'RET0') ;
  70.     RVX . 'RESULTS' . 'LISTLINF' = 'PROG' ;
  71.     RVX . 'RESULTS' . 'LISTITDT' = 'LECT' ;
  72.     RVX . 'RESULTS' . 'LISTITER' = 'LECT' ;
  73.     RVX . 'RESULTS' . 'NITER' = 0 ;
  74.  'FINSI' ;
  75.  
  76.  MDOMINT = RVX . 'MODEL' ;
  77.  NELT = 'NBEL' ('DOMA' MDOMINT 'CENTRE') ;
  78.  GAMN = RVX . 'PGAS' . 'GAMN' ;
  79. * 
  80.  LISTINCO = RVX . 'LISTCONS' ;
  81.  LISTPRIM = RVX . 'LISTPRIM' ;
  82.  LISTERR  = RVX . 'LISTERR'  ;
  83. *  
  84. * Upwind scheme
  85. *
  86.  METO = RVX . 'METHOD' ;
  87.  
  88. * Space accuracy (1 or 2) and limiter
  89. * Time accuracy (1 or 2)
  90.  
  91.  ORDESP  = RVX . 'SPACEA' ;
  92.  TYPELIM = RVX . 'LIMITER' ;
  93.  ORDTPS  = RVX . 'TIMEA' ;
  94.  
  95. * Initial/final time
  96. * Deltat or CFL
  97.  TPS     = RVX . 'RESULTS' . 'TPS' ;
  98.  TFINAL  = RVX . 'TFINAL' ;
  99.  'SI' ('EXISTE' RVX 'DTPS') ;
  100.     'SI' ('EXISTE' RVX 'CFL') ;
  101.         'MESSAGE' 'DTPS or CFL ???' ;
  102.         'ERREUR' 21 ;
  103.     'FINSI' ;
  104.     DTPS = RVX . 'DTPS' ;
  105.  'SINON' ;
  106.     CFL = RVX . 'CFL' ;
  107.  'FINSI' ;
  108.  
  109. * Dual time iterations
  110.  NDT = RVX . 'NDTITER' ;
  111. * Relative error 
  112.  EPSDT = RVX . 'EPSDT' ;
  113. * Jacobi iterations
  114. * NJAC = RVX . 'NJAC' ;
  115. * Cut off speed
  116.  ICO = RVX . 'CUTOFF' ;
  117.  
  118. *
  119. **** Conservative variables
  120. *
  121. *
  122.  MOT1 = 'EXTRAIRE' LISTINCO 1 ;
  123.  'SI' ('EGA' ('VALE' 'DIME') 2) ;
  124.     NOMMOM = 'EXTRAIRE' LISTINCO
  125.        ('LECT' 2 3 ) ;
  126.     NOMVEL = 'MOTS' 'UX' 'UY' ;
  127.     MOT2 = 'EXTRAIRE' LISTINCO 4 ;
  128.  'SINON' ;
  129.     NOMMOM = 'EXTRAIRE' LISTINCO
  130.       ('LECT' 2 3 4) ;
  131.     NOMVEL = 'MOTS' 'UX' 'UY' 'UZ' ;
  132.     MOT2 = 'EXTRAIRE' LISTINCO 5 ;
  133.  'FINSI' ;
  134.  
  135.  RN0   = 'REDU' (RVX . 'RESULTS' . 'RN')  ('DOMA' MDOMINT 'CENTRE') ;
  136.  GN0   = 'REDU' (RVX . 'RESULTS' . 'GN')  ('DOMA' MDOMINT 'CENTRE') ;
  137.  RET0  = 'REDU' (RVX . 'RESULTS' . 'RET') ('DOMA' MDOMINT 'CENTRE') ;
  138. *
  139. **** Primitive variables
  140. *
  141.  VN0 PN0 = 'PRIM' 'PERFMONO' RN0 GN0 RET0 GAMN ;
  142.  
  143.  MOTRN = 'EXTRAIRE' LISTPRIM 1 ;
  144.  'SI' ('EGA' ('VALE' 'DIME') 2) ;
  145.     MOTVN = 'EXTRAIRE' LISTPRIM
  146.        ('LECT' 2 3 ) ;
  147.     MOTPN = 'EXTRAIRE' LISTPRIM 4 ;
  148.  'SINON' ;
  149.     MOTVN = 'EXTRAIRE' LISTPRIM
  150.        ('LECT' 2 3 4) ;
  151.     MOTPN = 'EXTRAIRE' LISTPRIM 5 ;
  152.  'FINSI' ;
  153. *
  154. **** Coeff to compute gradients (MCHSCA, MCHVEC)
  155. *    Boundary conditions have to be taken into account
  156. *
  157.  RCHLIM RESLIM = PROLIM (RVX . 'PROLIM') MDOMINT LISTINCO LISTPRIM
  158.     RN0 VN0 PN0 GAMN  ;
  159.  
  160.  MAILLIM = 'EXTRAIRE' RCHLIM 'MAILLAGE' ;
  161.  
  162.  'SI' (MAILLIM 'EGA' 0) ;
  163.     GRADRN ALRN MCHSCA = 'PENT' MDOMINT 'CENTRE' 'EULESCAL' 'NOLIMITE'
  164.       ('MOTS' 'SCAL') RN0 ;
  165.     GRADVN ALVN MCHVEC = 'PENT' MDOMINT 'CENTRE' 'EULEVECT' 'NOLIMITE'
  166.       NOMVEL GN0 ;
  167.  'SINON' ;
  168.     SCALBC = 'MANUEL' 'CHPO' MAILLIM 1 'SCAL' 1.0 ;
  169.     'SI' (('VALEUR' 'DIME') 'EGA' 2) ; 
  170.         VECTBC = 'MANUEL' 'CHPO' MAILLIM
  171.             2 'UX' 0.0 'UY' 0.0 ;
  172. *       Same name as NOMVEL            
  173.     'SINON' ;
  174.         VECTBC = 'MANUEL' 'CHPO' MAILLIM
  175.             3 'UX' 0.0 'UY' 0.0 'UZ' ;
  176.     'FINSI' ;
  177.     GRADRN ALRN MCHSCA = 'PENT' MDOMINT 'CENTRE' 'EULESCAL' 'NOLIMITE'
  178.         ('MOTS' 'SCAL') RN0
  179.        'CLIM' SCALBC ;
  180.  
  181.     GRADVN ALVN MCHVEC = 'PENT' MDOMINT 'CENTRE' 'EULEVECT' 'NOLIMITE'
  182.         NOMVEL VN0
  183.        'CLIM' VECTBC ;
  184.  'FINSI' ;
  185. *
  186. *
  187. *** After each dual time loop, we could display
  188. *   the evolution of the error in the dual time loop
  189. *   the evolution of the dual time step (the safety factor)
  190. *
  191.  
  192.  LISTLINF = RVX . 'RESULTS' . 'LISTLINF' ;
  193.  LISTITDT = RVX . 'RESULTS' . 'LISTITDT' ; 
  194.  LISTITER = RVX . 'RESULTS' . 'LISTITER' ;
  195.  
  196. *  
  197.  
  198.  'MESSAGE' ;
  199.  'MESSAGE' ('CHAINE' 'Methode = ' METO) ;
  200.  'MESSAGE' ;
  201.  
  202.  'TEMPS' ;
  203.  
  204. ************************************************************************
  205. ************************************************************************
  206. **** Temporal loop *****************************************************
  207. ************************************************************************
  208. ************************************************************************
  209.  
  210.  RN_N1M1  = 'COPIER' RN0 ;
  211.  GN_N1M1  = 'COPIER' GN0 ;
  212.  RET_N1M1 = 'COPIER' RET0 ;
  213.  
  214.  AA = 'DIME' LISTITER ;
  215.  'SI' (AA > 0) ;
  216.     PTITER = 'EXTRAIRE' LISTITER AA ;
  217.  'SINON' ;
  218.     PTITER = 0 ;
  219.  'FINSI' ;
  220.  
  221.  DUSDT = 0.0D0 ;
  222.  
  223.  'REPETER' BLITER -1 ;
  224.  
  225.     PTITER = PTITER '+' 1 ;
  226. *
  227. **** Personal procedure
  228. *
  229.     PROCPT RVX ;
  230. *
  231. *
  232. ****  _N1M  = (t^n,\tau^m)
  233. *     _N1M1 = (t^n,\tau^{m+1})
  234. *
  235. *
  236. ************************************************************************
  237. ****** Loop on dual time***********************************************
  238. ************************************************************************
  239. *
  240. *
  241. *** DUSDT0 is the increment of DUSDT in the previous (physical) time
  242. *   iteration.
  243. *
  244.     DUSDT0 = DUSDT ;
  245.     DUSDT = 0.0D0 ;
  246.  
  247.     'REPETER' BLDT NDT ;
  248.  
  249.        RN_N1M  = RN_N1M1 ;
  250.        GN_N1M  = GN_N1M1 ;
  251.        RET_N1M = RET_N1M1 ;
  252.  
  253. *
  254. **** Primitive variables
  255. *
  256.  
  257.        VN_N1M PN_N1M = 'PRIM' 'PERFMONO' RN_N1M GN_N1M RET_N1M GAMN ;
  258. *
  259. *** Boundary conditions
  260. *
  261.  
  262.     RCHLIM RESLIM = PROLIM (RVX . 'PROLIM') MDOMINT LISTINCO LISTPRIM
  263.            RN_N1M VN_N1M PN_N1M GAMN  ;
  264.  
  265. *    
  266. ****** First/second order reconstruction
  267. *
  268.        'SI' (ORDESP 'EGA' 2) ;
  269.           GRADRN ALRN0  = 'PENT' MDOMINT 'CENTRE' 'EULESCAL' TYPELIM
  270.               ('MOTS' 'SCAL')  RN_N1M
  271.              'CLIM' ('EXCO' MOTRN RCHLIM 'SCAL') 'GRADGEO' MCHSCA ;
  272.           GRADPN ALPN0  = 'PENT' MDOMINT 'CENTRE' 'EULESCAL' TYPELIM
  273.               ('MOTS' 'SCAL')  PN_N1M
  274.              'CLIM' ('EXCO' MOTPN RCHLIM 'SCAL') 'GRADGEO' MCHSCA ;
  275.           GRADVN ALVN0  = 'PENT' MDOMINT 'CENTRE' 'EULEVECT' TYPELIM
  276.               ('MOTS' 'UX' 'UY') VN_N1M
  277.              'CLIM' ('EXCO' MOTVN RCHLIM MOTVN) 'GRADGEO' MCHVEC ;
  278.  
  279. *          'SI' (&BLDT < NLCB) ;
  280. *              ALRN0  =  'COPIER' ALRN ;
  281. *              ALPN0  =  'COPIER' ALPN ;
  282. *              ALVN0  =  'COPIER' ALVN ;
  283. *          'SINON' ;
  284. *             'SI' (&BLDT 'EGA' NLCB) ;
  285. *                'MESSAGE' ;
  286. *                'MESSAGE' 'On gele les limiteurs!!!' ;
  287. *                'MESSAGE' ;
  288. *             'FINSI' ;
  289. *          'FINSI' ;
  290.  
  291.           ROF VITF PF GAMF = 'PRET' 'PERFMONO'  2 1 
  292.                               MDOMINT
  293.                               RN_N1M  GRADRN ALRN0 
  294.                               VN_N1M  GRADVN ALVN0  
  295.                               PN_N1M   GRADPN ALPN0  
  296.                               GAMN  ;
  297.        'SINON' ;
  298.           ROF VITF PF GAMF = 'PRET' 'PERFMONO'  1 1 
  299.                               MDOMINT
  300.                               RN_N1M VN_N1M PN_N1M GAMN  ;
  301.        'FINSI' ;
  302.  
  303.        RESIDU DELTAT =  'KONV' 'VF' 'PERFMONO' 'RESI' METO
  304.             MDOMINT LISTINCO  ROF VITF PF GAMF MAILLIM 
  305.             ICO ICO ;
  306.        RESIDU = RESIDU '+' RESLIM ;
  307.  
  308. *
  309. **** La gravite
  310. *
  311.  
  312.        RESGRA = 'FIMP' 'VF' 'GRAVMONO' 'RESI' LISTINCO
  313.                  RN_N1M GN_N1M1 (RVX . 'GRAVITY') ;
  314.  
  315.        RESIDU = RESIDU '+' RESGRA ;
  316.  
  317. *
  318. ****** Residuum for dual tims stepping also involved the
  319. *      variation of the conserved variables with respect
  320. *      to time
  321.  
  322.        'SI' ((&BLITER 'EGA' 1) 'OU' (ORDTPS 'EGA' 1)) ;
  323.            RESIDU = RESIDU '-' DUSDT ;
  324.        'SINON' ;
  325.            RESIDU = RESIDU '-' ((1.5 '*' DUSDT) '-' (0.5 '*' DUSDT0)) ;
  326.        'FINSI' ;
  327.  
  328. *
  329. *** Time step at the first iteration/jacobi iteration
  330. *
  331.  
  332.        'SI' (&BLDT 'EGA' 1) ;
  333.           'SI' ('EXISTE' RVX 'CFL') ;
  334.              DTPS = (RVX . 'CFL') '*' 2.0D0 '*' DELTAT ;
  335.           'SINON' ;
  336.              DTPS = RVX . 'DTPS' ;
  337.           'FINSI' ;
  338.           DTPS = 'MINIMUM' ('PROG' DTPS ((TFINAL '-' TPS) '*' 1.001)) ;
  339.           TPS = TPS '+' DTPS ;
  340. *
  341.           NJAC = 'ENTIER' ('MINIMUM' (RVX . 'NJACITER')) ;
  342.           NJAC0 = NJAC ;
  343.        'SINON' ;
  344.           NJAC = 'IPOL' (('LOG' ERRINF) '/' ('LOG' 10))
  345.                  (RVX . 'NJACLERR') (RVX . 'NJACITER') ;
  346.           NJAC = 'MAXIMUM' ('PROG' NJAC0 NJAC) ;
  347.           NJAC = 'ENTIER' NJAC ;
  348.           NJAC0 = NJAC ;
  349.        'FINSI' ;
  350.  
  351. *
  352. *** JACOBI
  353. *   
  354.  
  355.         'SI' ('NON' ('EXISTE' RVX 'TYPEJAC')) ;
  356.            LOGBF = (&BLDT '/' 2 '*' 2) 'EGA' &BLDT ;
  357.            'SI' LOGBF ;
  358.               RVX . 'TYPEJAC' = 'LJACOF' ;
  359.            'SINON' ;
  360.               RVX . 'TYPEJAC' = 'LJACOB' ;
  361.            'FINSI' ;
  362.         'FINSI' ;
  363. *
  364. **** CFL dual
  365. *
  366.          'SI' (&BLDT 'EGA' 1) ;
  367.              SAFFACD = ('MINIMUM' (RVX . 'DCFL')) '*' 2 ;
  368.          'SINON' ;
  369.              SAFFACD = ('IPOL' (('LOG' ERRINF) '/' ('LOG' 10))
  370.                  (RVX . 'DCFLERR') (RVX . 'DCFL')) '*' 2 ;
  371.          'FINSI' ;
  372.  
  373.           'SI' ((&BLITER 'EGA' 1) 'OU' (ORDTPS 'EGA' 1)) ;
  374. *              DUN IPRO = 'DETO'  (RVX . 'TYPEJAC')
  375.               DUN IPRO = 'KONV' 'VF' 'PMON1FMM'  (RVX . 'TYPEJAC')
  376.                  LISTINCO MDOMINT  RESIDU RN_N1M GN_N1M RET_N1M GAMN ICO
  377.                  DTPS SAFFACD
  378.                  NJAC 'CLIM' LISTPRIM RCHLIM (RN_N1M '*' 0) ;
  379.           'SINON' ;
  380. *              DUN IPRO = 'DETO'  (RVX . 'TYPEJAC')
  381.               DUN IPRO = 'KONV' 'VF' 'PMON1FMM'  (RVX . 'TYPEJAC') 
  382.                  LISTINCO MDOMINT  RESIDU RN_N1M GN_N1M RET_N1M GAMN ICO
  383.                  (DTPS '/' 1.5)
  384.                  SAFFACD NJAC 'CLIM' LISTPRIM RCHLIM (RN_N1M '*' 0) ;
  385.           'FINSI' ;
  386. *       'FINSI' ;
  387.  
  388.        'SI' (IPRO 'NEG' 0) ;
  389.           'MESSAGE' ;
  390.           'MESSAGE' 'Probleme dans FMM' ;
  391.           'MESSAGE' ;
  392.           'ERREUR' 21 ;
  393.        'FINSI' ;
  394.  
  395. *
  396. **** We compute DUSDT for the future loop
  397. *
  398.  
  399.        DUSDT = DUSDT '+' (DUN '/' DTPS) ;
  400.  
  401. *
  402. **** We evaluate the conservative variables at t^{n+1}, \tau^{m+1}
  403. *
  404.  
  405.        DRN = 'EXCO' MOT1 DUN 'SCAL' ;
  406.        DGN = 'EXCO' NOMMOM DUN NOMVEL ;
  407.        DRET = 'EXCO' MOT2 DUN 'SCAL' ;
  408.  
  409.        RN_N1M1   = RN_N1M '+' DRN ;
  410.        GN_N1M1   = GN_N1M '+' DGN ;
  411.        RET_N1M1  = RET_N1M '+' DRET ;
  412.  
  413.        ERRINF =  'MAXIMUM' DUN 'ABS' LISTERR ;
  414.        LISTLINF = LISTLINF 'ET' ('PROG' ERRINF) ;
  415.        LISTITDT = LISTITDT 'ET' ('LECT'  &BLDT) ;
  416.        LISTITER = LISTITER 'ET' ('LECT' PTITER) ;
  417.  
  418.        'SI' (RVX . 'RELERR') ;
  419. *         Relative error       
  420.           'SI' (&BLDT 'EGA' 1) ;
  421.              ERRINF0 = ERRINF ;
  422.           'SINON' ;
  423.              'SI' (ERRINF < (EPSDT '*' ERRINF0)) ;
  424.                 'QUITTER' BLDT ;
  425.              'FINSI' ;
  426.           'FINSI' ;
  427.        'SINON' ;
  428.           'SI' (ERRINF < EPSDT) ;
  429.              'QUITTER' BLDT ;
  430.           'FINSI' ;
  431.        'FINSI' ;
  432.  
  433.        'SI' ((&BLDT 'EGA' 1) 'OU'
  434.              (((&BLDT '/' 20) '*' 20) 'EGA' &BLDT)) ;
  435.           'MESSAGE' ;
  436.           'MESSAGE'
  437.           ('CHAINE' 'ITER =' PTITER ' TPS =' TPS
  438.            ' DTITER ='  &BLDT  ' LINF =' ERRINF
  439.            ' DCFL =' SAFFACD ' NJAC =' NJAC) ;
  440.           'MESSAGE' ;
  441.        'FINSI' ;
  442.  
  443.     'FIN' BLDT ;
  444.  
  445. ************************************************************************
  446. ****** End of the loop on dual time*************************************
  447. ************************************************************************
  448. *
  449. *** Update of RVX . 'RESULTS'
  450. *
  451.     RVX . 'RESULTS' . 'RN'   = RN_N1M1 ;
  452.     RVX . 'RESULTS' . 'GN'   = GN_N1M1 ;
  453.     RVX . 'RESULTS' . 'RET' = RET_N1M1 ;
  454.     RVX . 'RESULTS' . 'LISTLINF' = LISTLINF ;
  455.     RVX . 'RESULTS' . 'LISTITDT' = LISTITDT ;
  456.     RVX . 'RESULTS' . 'LISTITER' = LISTITER ;
  457.     RVX . 'RESULTS' . 'TPS' = TPS ;
  458.     RVX . 'RESULTS' . 'NITER' = (RVX . 'RESULTS' . 'NITER') '+' 1 ;
  459. *
  460.     'SI' (TPS '>EG' TFINAL) ;
  461.        'QUITTER' BLITER ;
  462.     'FINSI' ;
  463.  
  464.  'FIN' BLITER ;
  465.  
  466.  TCPU = 'TEMPS' 'NOEC' ;
  467.  RVX . 'RESULTS' . 'TCPU' = TCPU ;
  468.  
  469.  'FINPROC' ;
  470.  
  471. *************************************************************************
  472. *************************************************************************
  473. *************************************************************************
  474. *************************************************************************
  475.  
  476.  TYEL = 'TRI3' ;
  477. * TYEL = 'QUA4' ;
  478.  
  479.  'OPTION'  'DIME' 2  'ELEM' TYEL  'ISOV' 'SULI'
  480.     'ECHO' 0 'TRAC' 'X' ;
  481.  
  482. * GRAPH = VRAI ;
  483.  GRAPH = FAUX ;
  484.  
  485. ******************
  486. ******************
  487. ******************
  488. **** MAILLAGE ****
  489. ******************
  490. ******************
  491. ******************
  492. ******************
  493. *
  494. *
  495.  
  496.  RAF = 2 ;
  497.  
  498.  NY = 5 '*' RAF ;
  499.  NX1 = 4 '*' RAF ;
  500.  NX2 = 2 '*' NX1 ;
  501.  NX3 = NX1  ;
  502.  NX = (NX1 '+' NX2 '+' NX3)  ;
  503.  DX = (4.0 '/' NX)  ;
  504.  
  505.  A0 = -2.0 0.0 ;
  506.  A1 = -1.0 0.0 ;
  507.  A2 = 1.0 0.0 ;
  508.  A3 = 2.0 0.0 ;
  509.  A4 = 2.0 1.0 ;
  510.  A5 = -2.0 1.0 ;
  511.  
  512. *
  513. **** LIGB
  514. *
  515.  
  516.  LIGB1 = A0 'DROIT' NX1 A1 ;
  517.  
  518. * LIGB2 
  519.  
  520.  xcel = ('COORDONNEE' 1 A1) '+' DX ;
  521.  ycel = 0.1 '*' ( 1.0 '+' ('COS' (180 '*' xcel)));
  522.  ACEL = xcel ycel ;
  523.  LIGB2 = A1 'DROIT' 1 ACEL ;
  524.  'REPETER' BL1 (NX2 '-' 2) ;
  525.     ACEL0 = ACEL ;
  526.     xcel = xcel '+' DX ;
  527.     ycel = 0.1 '*' ( 1.0 '+' ('COS' (180 '*' xcel)));
  528.     ACEL = xcel ycel ;
  529.     LIGB2 = LIGB2 'ET' (ACEL0 'DROIT' 1 ACEL) ;
  530.  'FIN' BL1;
  531.  LIGB2 = LIGB2 'ET' (ACEL 'DROIT' 1 A2) ;
  532.  
  533.  
  534.  LIGB3 = A2 'DROIT' NX3 A3 ;
  535.  
  536.  LIGB = LIGB1 'ET' LIGB2 'ET' LIGB3 ;
  537.  
  538. *
  539. **** LIGH
  540. *
  541.  
  542.  LIGH = A4 'DROIT' NX A5 ;
  543.  
  544. *
  545. **** DOMINT
  546. *
  547.  
  548.  DOMINT =  LIGH 'REGLER' NY ('INVERSE' LIGB) ;
  549. * DOMINT =  LIGB 'REGLER' NY ('INVERSE' LIGH) ;
  550.  LIGCON = 'CONTOUR' DOMINT ;
  551.  
  552. *
  553. *** LIGG
  554. *
  555.  
  556.  LIGG = LIGCON 'ELEM' 'COMP' A5  A0 ;
  557.  
  558. *
  559. **** LIGD
  560. *
  561.  
  562.  LIGD = LIGCON 'ELEM' 'COMP' A3  A4 ;
  563.  
  564.  
  565.  'SI' GRAPH ;
  566.     'TRACER' DOMINT  'TITRE'
  567.        ('CHAINE' 'Domaine, nbel=' ('NBEL' DOMINT)) ;
  568.  'FINSI' ;
  569.  
  570. *
  571. **** Creation of MODELS
  572. * 
  573.  
  574.  MDOMINT  = 'MODELISER' DOMINT 'EULER' ;
  575.  MLIGG    = 'MODELISER' LIGG   'EULER' ;
  576.  MLIGD    = 'MODELISER' LIGD   'EULER' ;
  577.  
  578.  TDOMINT  = 'DOMA' MDOMINT 'VF' ;
  579.  TLIGG    = 'DOMA' MLIGG   'VF' ;
  580.  TLIGD    = 'DOMA' MLIGD   'VF' ;
  581.  
  582.  QDOMINT = TDOMINT . 'QUAF' ;
  583.  QLIGD = TLIGD . 'QUAF' ;
  584.  QLIGG = TLIGG . 'QUAF' ;
  585.  
  586.  'ELIMINATION' QDOMINT (1.0D-3 '/' RAF)  QLIGG ;
  587.  'ELIMINATION' QDOMINT (1.0D-3 '/' RAF)  QLIGD ;
  588.  
  589. ***************************************************************
  590. ***************************************************************
  591. ***************************************************************
  592. ************** Initial conditions *****************************
  593. ***************************************************************
  594. ***************************************************************
  595. ***************************************************************
  596.  
  597.  GAMAIR = 1.4 ;
  598.  RO = 1.4 ;
  599.  P  = 1.0 ;
  600.  U  = 0.01 ;
  601.  UCO = 1.0D-3*U ;
  602.  
  603. ***************************************************************
  604. ***************************************************************
  605. ***************************************************************
  606. ************** Boundary Conditions ****************************
  607. ***************************************************************
  608. ***************************************************************
  609. ***************************************************************
  610. *
  611. * Subsonic inlet.
  612. * We impose the total enthalpy and the entropy.
  613. * The velocity direction is supposed to be 90 degrees
  614. *
  615.   HTIN = (GAMAIR '/' (GAMAIR '-' 1.0)) '*' P ;
  616.   ECIN = 0.5 '*' RO '*' U '*' U ;
  617.   HTIN = HTIN '+' ECIN ;
  618.   HTIN = HTIN '/' RO ;
  619.   SINL = P  '/' (RO '**' GAMAIR) ;
  620.  
  621.   CHPINL = 'MANUEL' 'CHPO' ('DOMA' MLIGG 'CENTRE') 2 'HT' HTIN
  622.            'S' sinl ;
  623. *           
  624. * Subsonic outlet
  625. * We impose the pressure
  626. *
  627.   CHPOUT = 'MANUEL' 'CHPO' ('DOMA' MLIGD 'CENTRE') 1 'PN' P ;
  628.  
  629. ***************************************************************
  630. ***** PROCEDURE POUR CALCULER LES VARIABLES CONSERVATIVES *****
  631. ***************************************************************
  632.  
  633.  'DEBPROC' CONS ;
  634.  'ARGUMENT' RN*'CHPOINT' VN*'CHPOINT' PN*'CHPOINT' GAMN*'CHPOINT' ;
  635.  
  636.  RVN  =  RN '*' VN ('MOTS' 'SCAL' 'SCAL') ('MOTS' 'UX' 'UY')
  637.     ('MOTS' 'UX' 'UY') ;
  638.  CELL = 'PSCAL' RVN VN ('MOTS' 'UX' 'UY') ('MOTS' 'UX' 'UY') ;
  639.  RECIN = 0.5 '*' CELL ;
  640.  REIN = PN '/' (GAMN '-' 1.0) ;
  641.  RETN = RECIN '+' REIN ;
  642.  
  643.  DETR CELL ;
  644.  DETR RECIN ;
  645.  DETR REIN ;
  646.  'RESPRO' RVN RETN ;
  647.  
  648.  'FINPROC' ;
  649.  
  650. ***************************************************************
  651. ***************************************************************
  652. ***************************************************************
  653.  
  654.  RN0  = 'MANU' 'CHPO' ('DOMA' MDOMINT 'CENTRE') 1 'SCAL' ro
  655.         'NATURE' 'DISCRET' ;
  656.  
  657.  VN0  = 'MANU' 'CHPO' ('DOMA' MDOMINT 'CENTRE') 2 'UX'  u
  658.         'UY' 0.0  'NATURE' 'DISCRET' ;
  659.  
  660.  PN0  = 'MANU' 'CHPO' ('DOMA' MDOMINT 'CENTRE') 1 'SCAL' p
  661.         'NATURE' 'DISCRET' ;
  662.  
  663.  GAMN =  'MANU' 'CHPO' ('DOMA' MDOMINT 'CENTRE') 1 'SCAL' GAMAIR
  664.         'NATURE' 'DISCRET' ;
  665.  
  666.  GN0 RET0 = CONS RN0 VN0 PN0 GAMN  ;
  667.  
  668.  VN1 PN1 = 'PRIM' 'PERFMONO' RN0 GN0 RET0 GAMN ;
  669.  
  670.  ERRO = 'MAXIMUM' (PN1 '-' PN0) 'ABS' ;
  671.  
  672.  'SI' (ERRO > 1.0D-6) ;
  673.     'MESSAGE' 'Problem in the ic file!!!'
  674.     'ERREUR' 5 ;
  675.  'FINSI' ;
  676.  
  677. *
  678. **** Plot of IC
  679. *
  680.  
  681.  MOD1  =  'MODELISER' ('DOMA' MDOMINT 'MAILLAGE' ) 'THERMIQUE';
  682.  
  683.  'SI' GRAPH ;
  684.  
  685.     CHM_RN   =  'KCHA' MDOMINT 'CHAM' RN0 ;
  686.     CHM_PN   =  'KCHA' MDOMINT 'CHAM' PN0 ;
  687.     CHM_VN  =  'KCHA' MDOMINT 'CHAM'  VN0 ;
  688.  
  689.     'TRACER'  CHM_RN  MOD1 
  690.        'TITR'  ('CHAINE'  'RN at t='  0.0);
  691.     'TRACER'  CHM_PN  MOD1 
  692.        'TITR'  ('CHAINE'  'PN at t='  0.0);
  693.     'TRACER'  CHM_VN MOD1
  694.        'TITR'  ('CHAINE' 'MACHN at t='  0.0);
  695.  
  696.  'FINSI' ;
  697.  
  698. *************************************************************************
  699. *************************************************************************
  700. *************************************************************************
  701. *************************************************************************
  702. ************   COMPUTATION OF THE SOLUTION  *****************************
  703. *************************************************************************
  704. *************************************************************************
  705. *************************************************************************
  706. *************************************************************************
  707.  
  708. *************************************************************************
  709. ******PROCEDURE PROCPT **************************************************
  710. *************************************************************************
  711. *
  712. *
  713.  'DEBPROC' PROCPT ;
  714.  'FINPROC' ;
  715. *
  716.  
  717. *************************************************************************
  718. ******PROCEDURE PROLIM **************************************************
  719. *************************************************************************
  720. *
  721.  'DEBPROC' PROLIM ;
  722.  'ARGUMENT' RVX*'TABLE' MDOMINT*'MMODEL' LISTINCO*'LISTMOTS'
  723.             LISTPRIM*'LISTMOTS'
  724.             RN*'CHPOINT' VN*'CHPOINT' PN*'CHPOINT' GAMN*'CHPOINT' ;
  725. *
  726. *** Boundary conditions
  727. *   Subsonic inlet and subsonic outlet
  728. *    
  729.  
  730.  RCHLIM1 RESLIM1 = 'KONV' 'VF' 'PERFMONO' 'CLIM'  'RESI'
  731.     MDOMINT (RVX . 'MODLI1') LISTINCO LISTPRIM
  732.     RN VN PN GAMN (RVX. 'CHPLI1') 'INSU' ;
  733.  
  734.  RCHLIM2 RESLIM2 = 'KONV' 'VF' 'PERFMONO' 'CLIM'  'RESI'
  735.     MDOMINT (RVX . 'MODLI2') LISTINCO LISTPRIM
  736.     RN VN PN GAMN (RVX. 'CHPLI2') 'OUTP' ;
  737.  
  738.  
  739. *   RCHLIM  = primitive variables on the boundary
  740. *   RESLIM  = contribution of boundary conditions to the residuum
  741.  
  742.  RESLIM  = RESLIM1 '+' RESLIM2 ;
  743.  RCHLIM  = RCHLIM1 '+' RCHLIM2 ;
  744.  
  745.  'RESPRO' RCHLIM RESLIM ;
  746.  'FINPROC' ;
  747.  
  748. *************************************************************************
  749. ****** FIN PROCEDURE PROLIM *********************************************
  750. *************************************************************************
  751.  
  752. **********************************************************
  753. **********************************************************
  754. **********************************************************
  755. ************** COMPUTATION OF THE SOLUTION ***************
  756. **********************************************************
  757. **********************************************************
  758. **********************************************************
  759. **********************************************************
  760. *
  761.  RV = 'TABLE' ;
  762.  RV . 'MODEL' = MDOMINT ;
  763. *
  764. **** Conservative variables / primitive variables
  765. *
  766.  MOTRN   = 'RN' ;
  767.  MOTGNX  = 'RUX' ;
  768.  MOTGNY  = 'RUY' ;
  769.  MOTVNX  = 'UX' ;
  770.  MOTVNY  = 'UY' ;
  771.  MOTVN   = 'MOTS' 'UX' 'UY' ;
  772.  NOMMOM  = 'MOTS' MOTGNX MOTGNY ;
  773.  MOTRET = 'RETN' ;
  774.  MOTPN   = 'PN' ;
  775. *
  776.  RV . 'LISTCONS' = 'MOTS' MOTRN MOTGNX MOTGNY MOTRET ;
  777.  RV . 'LISTPRIM' = 'MOTS' MOTRN MOTVNX MOTVNY MOTPN  ;
  778. *
  779.  RV . 'RN0' = RN0 ;
  780.  RV . 'GN0' = GN0 ;
  781.  RV . 'RET0' = RET0 ;
  782. *
  783. **** Gas property/gravity
  784. *
  785.  RV . 'PGAS' = 'TABLE' ;
  786.  RV . 'PGAS' . 'GAMN' = GAMN ;
  787.  RV . 'GRAVITY' = 'MANUEL' 'CHPO' ('DOMA' MDOMINT 'CENTRE') 2 'UX' 0.0
  788.                   'UY' 0.0 ;
  789. *
  790. * Table for BC
  791. *
  792.  RV . 'PROLIM' = 'TABLE' ;
  793.  RV . 'PROLIM' . 'MODLI1' = MLIGG ;
  794.  RV . 'PROLIM' . 'MODLI2' = MLIGD ;
  795.  RV . 'PROLIM' . 'CHPLI1' = CHPINL ;
  796.  RV . 'PROLIM' . 'CHPLI2' = CHPOUT ;
  797. *
  798. **** Numerical parameters
  799. *
  800. *
  801. * Variable to compute L2 error
  802. *
  803.  RV . 'LISTERR' = 'MOTS' MOTGNX MOTGNY  ;
  804. * RV . 'LISTERR' = 'MOTS' MOTRN ;
  805. * RV . 'LISTERR' = 'MOTS' MOTRET ;
  806. *  
  807. * Upwind scheme
  808. *
  809. * RV . 'METHOD' = 'RUSANOLM' ;
  810. * RV . 'METHOD' = 'HLLCLM' ;
  811. * RV . 'METHOD' = 'ROELM' ;
  812.  RV . 'METHOD' = 'AUSMPUP' ;
  813.  RV . 'CUTOFF' = 'MANU' 'CHPO' ('DOMA' MDOMINT 'CENTRE') 1 'SCAL' uco
  814.         'NATURE' 'DISCRET' ;
  815. *
  816. * Reconstruction/limiter
  817. * Time accuracy (1 or 2)
  818. * Iterations
  819. * Final time
  820.  RV . 'SPACEA' = 2 ;
  821.  RV . 'LIMITER' = 'NOLIMITE' ;
  822. * RV . 'LIMITER' = 'LIMITEUR' ;
  823.  RV . 'TIMEA' = 1 ;
  824. *
  825. **** Phisical time
  826. *
  827.  RV . 'T0' = 0 ;
  828.  RV . 'TFINAL' = 1.0D6 ;
  829. * RV . 'CFL' = 20. ;
  830.  RV . 'DTPS' = 1.0D16 ;  
  831. *
  832. **** Dual time
  833. * 
  834. *    Safety factor for the dual time step
  835. *    Max. Dual time iterations
  836. *    Relative error 
  837. *
  838.  RV . 'DCFLERR' = 'PROG' 16. -32. ;
  839.  RV . 'DCFL'    = 'PROG' 1.0D2 1.0D2 ;
  840.  RV . 'NDTITER' = 2000 ;
  841.  RV . 'RELERR' = FAUX ;
  842.  RV . 'EPSDT' = 1.0D-13 ;
  843. **** Jacobi iterations
  844.  RV . 'TYPEJAC' = 'PJACO' ;
  845. * RV . 'TYPEJAC' = 'LJACOF' ;
  846. * RV . 'TYPEJAC' = 'LJACOB' ;
  847. * RV . 'TYPEJAC' = 'LJACOFB' ;
  848.  RV . 'NJACITER' = 'PROG' 15   15 ;
  849.  RV . 'NJACLERR' = 'PROG' -1  -16 ;
  850. *
  851. **** Parameters for PROCPT
  852. *
  853. *
  854.  'TEMPS' 'ZERO' ;
  855.  PEUSM RV ;
  856.  'TEMPS' ;
  857.  
  858. * 'TEMPS' 'IMPR' ;
  859.  
  860.  RN  = RV . 'RESULTS' . 'RN'  ;
  861.  GN  = RV . 'RESULTS' . 'GN'  ;
  862.  RET = RV . 'RESULTS' . 'RET' ;
  863.  
  864.  LISTITER = RV . 'RESULTS' . 'LISTITER' ;
  865.  LISTITDT = RV . 'RESULTS' . 'LISTITDT' ;
  866.  LISTLINF = RV . 'RESULTS' . 'LISTLINF' ;
  867.  
  868. **********************************************************
  869. **********************************************************
  870. **********************************************************
  871. ************** PLOTS *************************************
  872. **********************************************************
  873. **********************************************************
  874. **********************************************************
  875. **********************************************************
  876.  
  877. *
  878. *** Convergence evolution inside of each iteration
  879. *
  880.  
  881.  'SI' (VRAI) ;
  882.     NITERE = 'DIME' LISTITER ;
  883.     I1 = 1 ;
  884.     AA = 'PROG' ;
  885.     BB = 'PROG' ;
  886.     CC = 'PROG' ;
  887.     'REPETER' BLITER NITERE ;
  888.        I0 = I1 ;
  889.        I1 = 'EXTRAIRE' LISTITER &BLITER ;
  890.        'SI' (I1 'EGA' I0) ;
  891.           AA = AA 'ET' ('PROG' ('EXTRAIRE' LISTITDT &BLITER)) ;
  892.           BB = BB 'ET' ('PROG' ('EXTRAIRE' LISTLINF &BLITER)) ;
  893.        'SINON' ;
  894.           everr = 'EVOL' 'MANU' 'niter' AA 'Log(Linf)'
  895.             (('LOG' (BB '+' ('PROG' ('DIME' BB) '*' 1.0D-20))) '/'
  896.              ('LOG' 10.)) ;
  897.           'SI' GRAPH ;
  898.              'DESSIN' everr 'TITRE' ('CHAINE' 'Convergence at iter '
  899.              (I1 '-' 1)) ;
  900.           'FINSI' ;
  901.           AA = 'PROG' ;
  902.           BB = 'PROG' ;
  903.        'FINSI' ;
  904.     'FIN' BLITER ;
  905.      everr = 'EVOL' 'MANU' 'niter' AA 'Log(Linf)'
  906.             (('LOG' (BB '+' ('PROG' ('DIME' BB) '*' 1.0D-20))) '/'
  907.              ('LOG' 10.)) ;
  908.     'SI' GRAPH ;
  909.         'DESSIN' everr 'TITRE' ('CHAINE' 'Convergence at iter '
  910.           (I1)) ;
  911.     'FINSI' ;
  912.   'FINSI' ;
  913.  
  914. *
  915. **** The 2D graphics
  916. *
  917.  
  918.   VN PN = 'PRIM' 'PERFMONO' RN GN RET GAMN 'TRICHE' ;
  919.   CN2 = GAMN '*' (PN '/' RN) ;
  920.   VN2 = 'PSCAL' VN VN MOTVN MOTVN ;
  921.   MACHN2 = VN2 '/' CN2 ;
  922.   MACHN = MACHN2 '**' 0.5 ;
  923.  
  924.   HTN = (GAMN '/' (GAMN '-' 1.0)) '*' (PN '/' RN) ;
  925.   ECIN = 0.5 '*' ('PSCAL' VN VN MOTVN MOTVN) ;
  926.   HTN = HTN '+' ECIN ;
  927.  
  928.   SN = PN '/' (RN '**' GAMAIR) ;
  929.  
  930.   tps = RV . 'RESULTS' . 'TPS' ;
  931.  
  932.  'SI' GRAPH ;
  933.  
  934.     CHM_RN  =  'KCHA'  MDOMINT  'CHAM'  RN  ;
  935.     CHM_VN  =  'KCHA'  MDOMINT  'CHAM'  VN  ;
  936.     CHM_PN  =  'KCHA'  MDOMINT  'CHAM'  PN  ;
  937.     CHM_MN  =  'KCHA'  MDOMINT  'CHAM'  MACHN  ;
  938.     CHM_HTN =  'KCHA'  MDOMINT  'CHAM'  HTN  ;
  939.     CHM_SN  =  'KCHA'  MDOMINT  'CHAM'  SN  ;
  940.  
  941.     'TRAC' CHM_RN  MDOMINT ('CONTOUR' DOMINT)
  942.        'TITR' ('CHAINE' 'rho at t=' TPS)   ;
  943.     'TRAC' CHM_VN  MDOMINT  ('CONTOUR' DOMINT)
  944.        'TITR' ('CHAINE' 'v at t= '  TPS)   ;
  945.     'TRAC' CHM_PN  MDOMINT ('CONTOUR' DOMINT)
  946.        'TITR' ('CHAINE' 'p at t= '  TPS)   ;
  947.     'TRAC' CHM_MN  MDOMINT ('CONTOUR' DOMINT)
  948.        'TITR' ('CHAINE' 'Mach at t= '  TPS)   ;
  949.     'TRAC' CHM_HTN  MDOMINT ('CONTOUR' DOMINT)
  950.        'TITR' ('CHAINE' 'ht at t= '  TPS ' hl =' HTL) ;
  951.     'TRAC' CHM_SN  MDOMINT ('CONTOUR' DOMINT)
  952.        'TITR' ('CHAINE' 'sn at t= '  TPS ' hl =' SL) ;
  953.  
  954.     'OPTION' 'ISOV' 'LIGN' ;
  955.     RNV = 'ELNO' TDOMINT ('KCHT' TDOMINT 'SCAL' 'CENTRE' RN) ;
  956.     'TRACER' DOMINT RNV ('CONTOUR' DOMINT) 15 'TITRE' 'ro';
  957.     PNV = 'ELNO' TDOMINT ('KCHT' TDOMINT 'SCAL' 'CENTRE' PN) ;
  958.     'TRACER' DOMINT PNV ('CONTOUR' DOMINT) 15 'TITRE' 'p' ;
  959.     MNV = 'ELNO' TDOMINT ('KCHT' TDOMINT 'SCAL' 'CENTRE' MACHN) ;
  960.     'TRACER' DOMINT MNV ('CONTOUR' DOMINT) 15 'TITRE' 'Mach' ;
  961.  
  962.  'FINSI' ;
  963.  
  964. * TEST
  965.  
  966.   LERR = 'EXTRAIRE' EVERR 'ORDO' ;
  967.   AA = 10. '**' ('EXTRAIRE' LERR 1) ;
  968.   BB = 10. '**' ('EXTRAIRE' LERR ('DIME' LERR)) ;
  969.   'SI' ((BB '/' AA) '>' 1.0D-5) ;
  970.       'ERREUR' 'Probleme de convergence' ;
  971.   'FINSI' ;
  972.   'FIN' ;
  973.  
  974.  
  975.  
  976.  
  977.  
  978.  
  979.  

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