Télécharger la plaquette Cast3M


fr en

Télécharger conjp7.eso

Retour à la liste

Numérotation des lignes : 

conjp7

  1. C CONJP7    SOURCE    CHAT      05/01/12    22:17:57     5004
  2.        SUBROUTINE CONJP7(JL,WVEC_L,WVEC_R,NVECT,TVECT,
  3.      &                   ga,v_inf)
  4. C************************************************************************
  5. C
  6. C PROJET            :  CASTEM 2000
  7. C
  8. C NOM               :  CONJP7
  9. C
  10. C DESCRIPTION       :  Voir KONJA5
  11. C
  12. C LANGAGE           :  FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
  13. C
  14. C AUTEUR            :  S. KUDRIAKOV, DM2S/SFME/LTMF
  15. C
  16. C************************************************************************
  17. C
  18. c----------------------------------------------------------------------
  19. c GENERAL DESCRIPTION:
  20. c   This subroutine provides the jacobian which is the derivatives
  21. c   of the numerical flux function defined at the wall
  22. c   with respect to the PRIMITIVE variables of the left
  23. c   cell (relative to the wall).
  24. c   The low-mach number corrections are made for the flux functions
  25. c
  26. c EQUATIONS: 2D Euler equations of gas dynamics
  27. c
  28. c
  29. c REFERENCE: 1) JCP, 129, 364-382 (1996)
  30. c               " A Sequel to AUSM: AUSM+ ".
  31. c                     M.S.Liou
  32. c            2) AIAA Journal, Sept. 1998
  33. c          " Low-Diffusion Flux-Splitting Methods for Flows at All Speeds"
  34. c                J.R.Edwards and M.S.Liou
  35. c----------------------------------------------------------------------
  36. c INPUT:
  37. c
  38. c      alpha   -- parameter of the AUSM+ scheme in the Pressure function;
  39. c                 ( -3/4 <= alpha <= 3/16 ) (imposed as a parameter)
  40. c
  41. c      beta    -- parameter of the AUSM+ scheme in the Mach function;
  42. c                 ( -1/16 <= beta <= 1/2 )  (imposed as a parameter)
  43. c
  44. c      wvec_l  -- vector of the primitive variables (rho,ux,uy,p) at the
  45. c                 left cell;
  46. c
  47. c      wvec_r  -- vector of the primitive variables (rho,ux,uy,p) at the
  48. c                 right cell;
  49. c
  50. c      nvect   -- normal vector to the interface (2 components in 2D);
  51. c
  52. c      tvect   -- tangential vector to the interface;
  53. c
  54. c      ga      -- ratio of the specific heats (assumed constant)
  55. c
  56. c      v_inf   -- parameter for reference velocity, see Ref. 2).
  57. c----------------------------------------------------------------------
  58. c
  59. c OUTPUT:
  60. c
  61. c      jl     -- jakobian matrix 4 by 4 - derivatives of the numerical
  62. c                 flux function with respect to the PRIMITIVE variables
  63. c                 from the left cell;
  64. c
  65. c----------------------------------------------------------------------
  66.       IMPLICIT INTEGER(I-N)
  67.        real*8 wvec_l(4),wvec_r(4)
  68.        real*8 nvect(2),tvect(2)
  69.        real*8 jl(4,4),jr(4,4)
  70.        real*8 alpha,beta
  71.        real*8 ga,gm1,temph
  72.        real*8 n_x,n_y,t_x,t_y
  73.        real*8 un_l, un_r, ut_l, ut_r
  74.        real*8 ml,mr,Mplus,Mmin,mmid
  75.        real*8 mpl_m, mmin_m,am
  76.        real*8 rold_l,uold_l,vold_l,pold_l,eold_l
  77.        real*8 rold_r,uold_r,vold_r,pold_r,eold_r
  78.        real*8 Pplus,Pmin
  79.        real*8 top,bot,bots
  80.        real*8 br3,br4,temp_l,temp_r,brac_l,brac_r
  81.        real*8 aleft, arigh,tcoef,bcoef
  82.        real*8 damr_l,damr_r,damu_l,damu_r
  83.        real*8 damv_l,damv_r,damp_l,damp_r
  84.        real*8 dmlr_l,dmlr_r,dmlu_l,dmlu_r
  85.        real*8 dmlv_l,dmlv_r,dmlp_l,dmlp_r
  86.        real*8 dmrr_l,dmrr_r,dmru_l,dmru_r
  87.        real*8 dmrv_l,dmrv_r,dmrp_l,dmrp_r
  88.        real*8 dMpr_l,dMpr_r,dMpu_l,dMpu_r
  89.        real*8 dMpv_l,dMpv_r,dMpp_l,dMpp_r
  90.        real*8 dMmr_l,dMmr_r,dMmu_l,dMmu_r
  91.        real*8 dMmv_l,dMmv_r,dMmp_l,dMmp_r
  92.        real*8 dmir_l,dmir_r,dmiu_l,dmiu_r
  93.        real*8 dmiv_l,dmiv_r,dmip_l,dmip_r
  94.        real*8 d3mr_l,d3mr_r,d3mu_l,d3mu_r
  95.        real*8 d3mv_l,d3mv_r,d3mp_l,d3mp_r
  96.        real*8 d2mr_l,d2mr_r,d2mu_l,d2mu_r
  97.        real*8 d2mv_l,d2mv_r,d2mp_l,d2mp_r
  98.        real*8 dPpr_l,dPpr_r,dPpu_l,dPpu_r
  99.        real*8 dPpv_l,dPpv_r,dPpp_l,dPpp_r
  100.        real*8 dPmr_l,dPmr_r,dPmu_l,dPmu_r
  101.        real*8 dPmv_l,dPmv_r,dPmp_l,dPmp_r
  102.        real*8 dpir_l,dpir_r,dpiu_l,dpiu_r
  103.        real*8 dpiv_l,dpiv_r,dpip_l,dpip_r
  104.        real*8 epsil,qq,amw,Mmin1,Mplus1
  105.        real*8 fmid,mlw,mrw,termp
  106.        real*8 ur_r,ur_l,urm,mhalf,mhalfr
  107.        real*8 durr_l,durr_r,duru_l,duru_r
  108.        real*8 durv_l,durv_r,durp_l,durp_r
  109.        real*8 dmhr_l,dmhr_r,dmhu_l,dmhu_r
  110.        real*8 dmhv_l,dmhv_r,dmhp_l,dmhp_r
  111.        real*8 dmfr_l,dmfr_r,dmfu_l,dmfu_r
  112.        real*8 dmfv_l,dmfv_r,dmfp_l,dmfp_r
  113.        real*8 dfm_mf,dfm_mh
  114.        real*8 dfmr_l,dfmr_r,dfmu_l,dfmu_r
  115.        real*8 dfmv_l,dfmv_r,dfmp_l,dfmp_r
  116.        real*8 m1mr_l,m1mr_r,m1mu_l,m1mu_r
  117.        real*8 m1mv_l,m1mv_r,m1mp_l,m1mp_r
  118.        real*8 m1pr_l,m1pr_r,m1pu_l,m1pu_r
  119.        real*8 m1pv_l,m1pv_r,m1pp_l,m1pp_r
  120.        real*8 tmpr_l,tmpr_r,tmpu_l,tmpu_r
  121.        real*8 tmpv_l,tmpv_r,tmpp_l,tmpp_r
  122.        real*8 coef,canc,v_inf
  123.        real*8 sr_l,sr_r,su_l,su_r,sv_l,sv_r,sp_l,sp_r
  124.        real*8 rum,rumr_l,rumu_l,rumv_l,rump_l
  125.        real*8 rumr_r,rumu_r,rumv_r,rump_r
  126.        integer i
  127.        parameter(alpha = 0.1875D0, beta = 0.125D0)
  128.        parameter(epsil = 1.0d0)
  129.        canc=1.0d0
  130. c-------------------------------------------------------------
  131.        n_x=nvect(1)
  132.        n_y=nvect(2)
  133.        t_x=tvect(1)
  134.        t_y=tvect(2)
  135. c----------------------------
  136.        gm1=ga-1.0d0
  137. c----------------------------
  138.        rold_l=wvec_l(1)
  139.        uold_l=wvec_l(2)
  140.        vold_l=wvec_l(3)
  141.        pold_l=wvec_l(4)
  142. c-----------------------
  143.        rold_r=wvec_r(1)
  144.        uold_r=wvec_r(2)
  145.        vold_r=wvec_r(3)
  146.        pold_r=wvec_r(4)
  147. c------------------------------------------------------------------
  148. c Computation of the specific total energy on the left and right.
  149. c------------------------------------------------------------------
  150.        eold_l=(uold_l*uold_l+vold_l*vold_l)/2.0d0
  151.        eold_l=eold_l+pold_l/(gm1*rold_l)
  152.        eold_r=(uold_r*uold_r+vold_r*vold_r)/2.0d0
  153.        eold_r=eold_r+pold_r/(gm1*rold_r)
  154. c------------------------------------------------------------------
  155. c   Computing reference velocity and its derivatives
  156. c------------------------------------------------------------------
  157.        aleft=sqrt(ga*pold_l/rold_l)
  158.        arigh=sqrt(ga*pold_r/rold_r)
  159.        qq=sqrt(uold_l*uold_l+vold_l*vold_l)
  160.        if(qq .lt. (epsil*v_inf)) then
  161.           ur_l = epsil*v_inf
  162.           durr_l=0.0d0
  163.           duru_l=0.0d0
  164.           durv_l=0.0d0
  165.           durp_l=0.0d0
  166.        else
  167.           ur_l=qq
  168.           durr_l=0.0d0
  169.           duru_l=uold_l/qq
  170.           durv_l=vold_l/qq
  171.           durp_l=0.0d0
  172.        endif
  173. c-----------------------
  174.        if(ur_l .ge. aleft) then
  175.           ur_l=aleft
  176.           durr_l=-aleft/(2.0d0*rold_l)
  177.           duru_l=0.0d0
  178.           durv_l=0.0d0
  179.           durp_l=ga/(2.0d0*aleft*rold_l)
  180.        endif
  181. c------------------------------------------------------------------
  182. c-------------------------------------------------------------------
  183.        qq=sqrt(uold_r*uold_r+vold_r*vold_r)
  184.        if(qq .lt. (epsil*v_inf)) then
  185.           ur_r = epsil*v_inf
  186.           durr_r=0.0d0
  187.           duru_r=0.0d0
  188.           durv_r=0.0d0
  189.           durp_r=0.0d0
  190.        else
  191.           ur_r=qq
  192.           durr_r=0.0d0
  193.           duru_r=uold_r/qq
  194.           durv_r=vold_r/qq
  195.           durp_r=0.0d0
  196.        endif
  197. c----------------------------
  198.        if(ur_r .ge. arigh) then
  199.           ur_r=arigh
  200.           durr_r=-arigh/(2.0d0*rold_r)
  201.           duru_r=0.0d0
  202.           durv_r=0.0d0
  203.           durp_r=ga/(2.0d0*arigh*rold_r)
  204.        endif
  205. c------------------------------------------------------------------
  206. c  Reference velocity at the interface is taken as an average
  207. c  of the reference velocities of the neighbouring cells
  208. c-------------------------------------------------------------------
  209.        urm=0.5d0*(ur_l+ur_r)
  210.        durr_l=0.5d0*durr_l
  211.        duru_l=0.5d0*duru_l
  212.        durv_l=0.5d0*durv_l
  213.        durp_l=0.5d0*durp_l
  214. c-------------------------
  215.        durr_r=0.5d0*durr_r
  216.        duru_r=0.5d0*duru_r
  217.        durv_r=0.5d0*durv_r
  218.        durp_r=0.5d0*durp_r
  219. c-------------------------------------------------------------------
  220. c Computation of the speed of sound and its derivatives;
  221. c numerical speed of sound at the interface is taken as an average
  222. c of the speeds of sounds of the neighbouring cells
  223. c-------------------------------------------------------------------
  224.        am=0.5d0*(aleft+arigh)
  225. c-------------------------------------------------------------------
  226.        if(abs(urm/am-1.0d0).le.0.000001d0) then
  227.          coef=0.0d0
  228.        else
  229.          coef=1.0d0
  230.        endif
  231. c-------------------------------------------------------------------
  232.        damr_r=-arigh/(4.0d0*rold_r)
  233.        damu_r=0.0d0
  234.        damv_r=0.0d0
  235.        damp_r=ga/(4.0d0*arigh*rold_r)
  236. c-----------------------
  237.        damr_l=-aleft/(4.0d0*rold_l)
  238.        damu_l=0.0d0
  239.        damv_l=0.0d0
  240.        damp_l=ga/(4.0d0*aleft*rold_l)
  241. c-------------------------------------------------------------------
  242. c Computing numerical Mach number and its derivatives,
  243. c see p.370, under (A1).
  244. c-------------------------------------------------------------------
  245.        un_l=uold_l*n_x+vold_l*n_y
  246.        un_r=uold_r*n_x+vold_r*n_y
  247.        ut_l=uold_l*t_x+vold_l*t_y
  248.        ut_r=uold_r*t_x+vold_r*t_y
  249. c-------------------------------------------------------------------
  250.        ml=un_l/am
  251.        mr=un_r/am
  252.        mhalf=0.5d0*(un_l+un_r)/am
  253. c---------------------
  254.        top=0.5d0*(un_l+un_r)/(am*am)
  255.        dmhr_l=-top*damr_l
  256.        dmhu_l=n_x/2.0d0/am-top*damu_l
  257.        dmhv_l=n_y/2.0d0/am-top*damv_l
  258.        dmhp_l=-top*damp_l
  259. c---------------------
  260.        dmhr_r=-top*damr_r
  261.        dmhu_r=n_x/2.0d0/am-top*damu_r
  262.        dmhv_r=n_y/2.0d0/am-top*damv_r
  263.        dmhp_r=-top*damp_r
  264. c--------------------------------
  265.        mhalfr=urm/am
  266. c---------------------
  267.        top=urm/(am*am)
  268.        dmfr_l=durr_l/am-top*damr_l
  269.        dmfu_l=duru_l/am-top*damu_l
  270.        dmfv_l=durv_l/am-top*damv_l
  271.        dmfp_l=durp_l/am-top*damp_l
  272. c---------------------
  273.        dmfr_r=durr_r/am-top*damr_r
  274.        dmfu_r=duru_r/am-top*damu_r
  275.        dmfv_r=durv_r/am-top*damv_r
  276.        dmfp_r=durp_r/am-top*damp_r
  277. c-------------------------------------------------------------------
  278. c Scaling function for the speed of sound and its derivatives
  279. c-------------------------------------------------------------------
  280.        top=(1.0d0-mhalfr*mhalfr)*(1.0d0-mhalfr*mhalfr)
  281.        top=top*mhalf*mhalf+4.0d0*mhalfr*mhalfr
  282.        bot=1.0d0+mhalfr*mhalfr
  283.        if(abs(canc-0.0d0).lt.0.000001d0) then
  284.          fmid=1.0d0
  285.        else
  286.          fmid=sqrt(top)/bot
  287.        endif
  288. c--------------------------
  289.        temph=-4.0d0*(1.0d0-mhalfr*mhalfr)*mhalfr
  290.        temph=temph*mhalf*mhalf+8.0d0*mhalfr
  291.        dfm_mf=temph/(sqrt(top)*2.0d0*bot)
  292.        dfm_mf=dfm_mf-sqrt(top)*2.0d0*mhalfr/(bot*bot)
  293. c--------------------------
  294.        temph=2.0d0*(1.0d0-mhalfr*mhalfr)*(1.0d0-mhalfr*mhalfr)*mhalf
  295.        dfm_mh=temph/(2.0d0*bot*sqrt(top))
  296. c--------------------------
  297.        dfmr_l=dfm_mf*dmfr_l+dfm_mh*dmhr_l
  298.        dfmu_l=dfm_mf*dmfu_l+dfm_mh*dmhu_l
  299.        dfmv_l=dfm_mf*dmfv_l+dfm_mh*dmhv_l
  300.        dfmp_l=dfm_mf*dmfp_l+dfm_mh*dmhp_l
  301. c--------------------------
  302.        dfmr_r=dfm_mf*dmfr_r+dfm_mh*dmhr_r
  303.        dfmu_r=dfm_mf*dmfu_r+dfm_mh*dmhu_r
  304.        dfmv_r=dfm_mf*dmfv_r+dfm_mh*dmhv_r
  305.        dfmp_r=dfm_mf*dmfp_r+dfm_mh*dmhp_r
  306. c--------------------------
  307.        amw=fmid*am
  308.        mlw=un_l/amw
  309.        mrw=un_r/amw
  310. c--------------------------
  311.        damr_l=canc*coef*dfmr_l*am+fmid*damr_l
  312.        damu_l=canc*coef*dfmu_l*am+fmid*damu_l
  313.        damv_l=canc*coef*dfmv_l*am+fmid*damv_l
  314.        damp_l=canc*coef*dfmp_l*am+fmid*damp_l
  315. c--------------------------
  316.        damr_r=canc*coef*dfmr_r*am+fmid*damr_r
  317.        damu_r=canc*coef*dfmu_r*am+fmid*damu_r
  318.        damv_r=canc*coef*dfmv_r*am+fmid*damv_r
  319.        damp_r=canc*coef*dfmp_r*am+fmid*damp_r
  320. c-/-/-/-/-/-/-/-/-/-/-/-/-/-/-/-/-/
  321.        am=amw
  322. c-------------------------------------------------------------------
  323. c  Redefinition of the numerical mach numbers
  324. c-------------------------------------------------------------------
  325.        if(abs(canc-0.0d0).lt.0.000001d0) then
  326.          top=2.0d0
  327.          bot=0.0d0
  328.        else
  329.          top=1.0d0+mhalfr*mhalfr
  330.          bot=1.0d0-mhalfr*mhalfr
  331.        endif
  332.        ml=0.5d0*(top*mlw+bot*mrw)
  333.        mr=0.5d0*(top*mrw+bot*mlw)
  334. c-------------------------------------------------------------------
  335. c  Mplus and Mmin are calligraphic lettes M+ and M- from the paper,
  336. c  see (19a) and (19b), p.367.
  337. c-------------------------------------------------------------------
  338.        if(abs(ml) .ge. 1.0d0) then
  339.          Mplus=(ml+abs(ml))/2.0d0
  340.        else
  341.          Mplus=(ml+1.0d0)*(ml+1.0d0)/4.0d0
  342.          Mplus=Mplus+beta*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  343.        endif
  344.        Mplus1=(ml+abs(ml))/2.0d0
  345. c-------------------------------------------------------------------
  346.        if(abs(mr) .ge. 1.0d0) then
  347.          Mmin=(mr-abs(mr))/2.0d0
  348.        else
  349.          Mmin=-(mr-1.0d0)*(mr-1.0d0)/4.0d0
  350.          Mmin=Mmin-beta*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  351.        endif
  352.        Mmin1=(mr-abs(mr))/2.0d0
  353. c-------------------------------------------------------------------
  354. c Derivatives of ml and mr with respect to both sets of primitive
  355. c variables: left  and right.
  356. c-------------------------------------------------------------------
  357.        temp_l=-un_l/(am*am)
  358.        temp_r=-un_r/(am*am)
  359. c--------
  360.        dmlr_l=temp_l*damr_l
  361.        dmlr_r=temp_l*damr_r
  362.        dmrr_l=temp_r*damr_l
  363.        dmrr_r=temp_r*damr_r
  364. c--------
  365.        dmlu_l=n_x/am+temp_l*damu_l
  366.        dmlu_r=temp_l*damu_r
  367.        dmru_l=temp_r*damu_l
  368.        dmru_r=n_x/am+temp_r*damu_r
  369. c--------
  370.        dmlv_l=n_y/am+temp_l*damv_l
  371.        dmlv_r=temp_l*damv_r
  372.        dmrv_l=temp_r*damv_l
  373.        dmrv_r=n_y/am+temp_r*damv_r
  374. c--------
  375.        dmlp_l=temp_l*damp_l
  376.        dmlp_r=temp_l*damp_r
  377.        dmrp_l=temp_r*damp_l
  378.        dmrp_r=temp_r*damp_r
  379. c-----------------------------
  380.        sr_l=dmlr_l
  381.        sr_r=dmlr_r
  382.        su_l=dmlu_l
  383.        su_r=dmlu_r
  384.        sv_l=dmlv_l
  385.        sv_r=dmlv_r
  386.        sp_l=dmlp_l
  387.        sp_r=dmlp_r
  388. c-----------------------------------------------------------------
  389. c   Redefinition of the derivatives of the ml & mr
  390. c-----------------------------------------------------------------
  391.        temp_l=(mlw-mrw)*mhalfr
  392.        temp_r=-temp_l
  393. c--------
  394.        dmlr_l=0.5d0*(top*dmlr_l+bot*dmrr_l)+canc*coef*temp_l*dmfr_l
  395.        dmlu_l=0.5d0*(top*dmlu_l+bot*dmru_l)+canc*coef*temp_l*dmfu_l
  396.        dmlv_l=0.5d0*(top*dmlv_l+bot*dmrv_l)+canc*coef*temp_l*dmfv_l
  397.        dmlp_l=0.5d0*(top*dmlp_l+bot*dmrp_l)+canc*coef*temp_l*dmfp_l
  398. c--------
  399.        dmlr_r=0.5d0*(top*dmlr_r+bot*dmrr_r)+canc*coef*temp_l*dmfr_r
  400.        dmlu_r=0.5d0*(top*dmlu_r+bot*dmru_r)+canc*coef*temp_l*dmfu_r
  401.        dmlv_r=0.5d0*(top*dmlv_r+bot*dmrv_r)+canc*coef*temp_l*dmfv_r
  402.        dmlp_r=0.5d0*(top*dmlp_r+bot*dmrp_r)+canc*coef*temp_l*dmfp_r
  403. c--------
  404.        dmrr_l=0.5d0*(top*dmrr_l+bot*sr_l)+canc*coef*temp_r*dmfr_l
  405.        dmru_l=0.5d0*(top*dmru_l+bot*su_l)+canc*coef*temp_r*dmfu_l
  406.        dmrv_l=0.5d0*(top*dmrv_l+bot*sv_l)+canc*coef*temp_r*dmfv_l
  407.        dmrp_l=0.5d0*(top*dmrp_l+bot*sp_l)+canc*coef*temp_r*dmfp_l
  408. c--------
  409.        dmrr_r=0.5d0*(top*dmrr_r+bot*sr_r)+canc*coef*temp_r*dmfr_r
  410.        dmru_r=0.5d0*(top*dmru_r+bot*su_r)+canc*coef*temp_r*dmfu_r
  411.        dmrv_r=0.5d0*(top*dmrv_r+bot*sv_r)+canc*coef*temp_r*dmfv_r
  412.        dmrp_r=0.5d0*(top*dmrp_r+bot*sp_r)+canc*coef*temp_r*dmfp_r
  413. c-----------------------------------------------------------
  414. c mmid is m_{1/2} (notation as in the paper, see (13),p.366)
  415. c-----------------------------------------------------------
  416.        mmid=Mplus+Mmin
  417. c-----------------------------------------------------------
  418. c Computing the derivatives of M+ and M-
  419. c-----------------------------------------------------------
  420.        if(ml .ge. 1.0d0) then
  421.          dMpr_l=dmlr_l
  422.          dMpu_l=dmlu_l
  423.          dMpv_l=dmlv_l
  424.          dMpp_l=dmlp_l
  425. c--------------------
  426.          dMpr_r=dmlr_r
  427.          dMpu_r=dmlu_r
  428.          dMpv_r=dmlv_r
  429.          dMpp_r=dmlp_r
  430.        else
  431.         if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  432.          temph=(ml+1.0d0)/2.0d0
  433.          dMpr_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_l
  434.          dMpu_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_l
  435.          dMpv_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_l
  436.          dMpp_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_l
  437. c--------------------
  438.            dMpr_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_r
  439.            dMpu_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_r
  440.            dMpv_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_r
  441.            dMpp_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_r
  442.          else
  443.            dMpr_l=0.0d0
  444.            dMpu_l=0.0d0
  445.            dMpv_l=0.0d0
  446.            dMpp_l=0.0d0
  447. c---------------------
  448.            dMpr_r=0.0d0
  449.            dMpu_r=0.0d0
  450.            dMpv_r=0.0d0
  451.            dMpp_r=0.0d0
  452.          endif
  453.        endif
  454. c-----------------------------------------------------------
  455. c   addition of low Mach number
  456. c-----------------------------------------------------------
  457.        if(ml .ge. 0.0d0) then
  458.          m1pr_l=dmlr_l
  459.          m1pu_l=dmlu_l
  460.          m1pv_l=dmlv_l
  461.          m1pp_l=dmlp_l
  462. c---------------------
  463.          m1pr_r=dmlr_r
  464.          m1pu_r=dmlu_r
  465.          m1pv_r=dmlv_r
  466.          m1pp_r=dmlp_r
  467.        else
  468.          m1pr_l=0.0d0
  469.          m1pu_l=0.0d0
  470.          m1pv_l=0.0d0
  471.          m1pp_l=0.0d0
  472. c--------------------
  473.          m1pr_r=0.0d0
  474.          m1pu_r=0.0d0
  475.          m1pv_r=0.0d0
  476.          m1pp_r=0.0d0
  477.        endif
  478. c-----------------------------------------------------------
  479.        if(mr .ge. 1.0d0) then
  480.          dMmr_l=0.0d0
  481.          dMmu_l=0.0d0
  482.          dMmv_l=0.0d0
  483.          dMmp_l=0.0d0
  484. c---------------------
  485.          dMmr_r=0.0d0
  486.          dMmu_r=0.0d0
  487.          dMmv_r=0.0d0
  488.          dMmp_r=0.0d0
  489.        else
  490.          if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  491.            temph=(-mr+1.0d0)/2.0d0
  492.            dMmr_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_l
  493.            dMmu_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_l
  494.            dMmv_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_l
  495.            dMmp_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_l
  496. c--------------------
  497.            dMmr_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_r
  498.            dMmu_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_r
  499.            dMmv_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_r
  500.            dMmp_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_r
  501.          else
  502.            dMmr_l=dmrr_l
  503.            dMmu_l=dmru_l
  504.            dMmv_l=dmrv_l
  505.            dMmp_l=dmrp_l
  506. c--------------------
  507.            dMmr_r=dmrr_r
  508.            dMmu_r=dmru_r
  509.            dMmv_r=dmrv_r
  510.            dMmp_r=dmrp_r
  511.          endif
  512.        endif
  513. c-----------------------------------------------------------
  514. c   addition of low Mach number
  515. c-----------------------------------------------------------
  516.        if(mr .le. 0.0d0) then
  517.          m1mr_l=dmrr_l
  518.          m1mu_l=dmru_l
  519.          m1mv_l=dmrv_l
  520.          m1mp_l=dmrp_l
  521. c---------------------
  522.          m1mr_r=dmrr_r
  523.          m1mu_r=dmru_r
  524.          m1mv_r=dmrv_r
  525.          m1mp_r=dmrp_r
  526.        else
  527.          m1mr_l=0.0d0
  528.          m1mu_l=0.0d0
  529.          m1mv_l=0.0d0
  530.          m1mp_l=0.0d0
  531. c--------------------
  532.          m1mr_r=0.0d0
  533.          m1mu_r=0.0d0
  534.          m1mv_r=0.0d0
  535.          m1mp_r=0.0d0
  536.        endif
  537. c-----------------------------------------------------------------
  538. c computing the derivatives of m_{1/2} (notation as in the paper)
  539. c-----------------------------------------------------------------
  540.        dmir_l=dMpr_l+dMmr_l
  541.        dmir_r=dMpr_r+dMmr_r
  542. c-------------
  543.        dmiu_l=dMpu_l+dMmu_l
  544.        dmiu_r=dMpu_r+dMmu_r
  545. c-------------
  546.        dmiv_l=dMpv_l+dMmv_l
  547.        dmiv_r=dMpv_r+dMmv_r
  548. c-------------
  549.        dmip_l=dMpp_l+dMmp_l
  550.        dmip_r=dMpp_r+dMmp_r
  551. c----------------------------------------------------------------
  552. c computing the main convective variables and their derivatives
  553. c mpl_m is m^{+}_{1/2} (paper's notation) and
  554. c mmin_m  is m^{-}_{1/2} (paper's notation), see (A2) on p.370.
  555. c----------------------------------------------------------------
  556.        termp=(Mmin1-Mmin+Mplus-Mplus1)*(1.0d0/(mhalfr*mhalfr)-1.0d0)
  557.        termp=termp*(pold_l-pold_r)/(pold_l/rold_l+pold_r/rold_r)
  558. c-------------------------------------------------------------
  559. c    derivatives of the termp
  560. c-------------------------------------------------------------
  561.        top=(Mmin1-Mmin+Mplus-Mplus1)
  562.        bots=1.0d0/(pold_l/rold_l+pold_r/rold_r)
  563.        bot=(pold_l-pold_r)*bots
  564.        temph=1.0d0/(mhalfr*mhalfr)-1.0d0
  565. c---------------------------
  566.        tmpr_l=(m1mr_l-dMmr_l+dMpr_l-m1pr_l)*bot*temph
  567.        tmpr_l=tmpr_l-2.0d0*bot*top*dmfr_l/(mhalfr*mhalfr*mhalfr)
  568.        tmpr_l=tmpr_l+temph*top*bot*bots*(pold_l/rold_l/rold_l)
  569. c---------------------------
  570.        tmpu_l=(m1mu_l-dMmu_l+dMpu_l-m1pu_l)*bot*temph
  571.        tmpu_l=tmpu_l-2.0d0*bot*top*dmfu_l/(mhalfr*mhalfr*mhalfr)
  572. c---------------------------
  573.        tmpv_l=(m1mv_l-dMmv_l+dMpv_l-m1pv_l)*bot*temph
  574.        tmpv_l=tmpv_l-2.0d0*bot*top*dmfv_l/(mhalfr*mhalfr*mhalfr)
  575. c---------------------------
  576.        tmpp_l=(m1mp_l-dMmp_l+dMpp_l-m1pp_l)*bot*temph
  577.        tmpp_l=tmpp_l-2.0d0*bot*top*dmfp_l/(mhalfr*mhalfr*mhalfr)
  578.        tmpp_l=tmpp_l+temph*top*bots*(1.0d0-bot/rold_l)
  579. c------------rrrrrrrr-------
  580. c------------rrrrrrrr-------
  581.        tmpr_r=(m1mr_r-dMmr_r+dMpr_r-m1pr_r)*bot*temph
  582.        tmpr_r=tmpr_r-2.0d0*bot*top*dmfr_r/(mhalfr*mhalfr*mhalfr)
  583.        tmpr_r=tmpr_r+temph*top*bot*bots*(pold_r/rold_r/rold_r)
  584. c---------------------------
  585.        tmpu_r=(m1mu_r-dMmu_r+dMpu_r-m1pu_r)*bot*temph
  586.        tmpu_r=tmpu_r-2.0d0*bot*top*dmfu_r/(mhalfr*mhalfr*mhalfr)
  587. c---------------------------
  588.        tmpv_r=(m1mv_r-dMmv_r+dMpv_r-m1pv_r)*bot*temph
  589.        tmpv_r=tmpv_r-2.0d0*bot*top*dmfv_r/(mhalfr*mhalfr*mhalfr)
  590. c---------------------------
  591.        tmpp_r=(m1mp_r-dMmp_r+dMpp_r-m1pp_r)*bot*temph
  592.        tmpp_r=tmpp_r-2.0d0*bot*top*dmfp_r/(mhalfr*mhalfr*mhalfr)
  593.        tmpp_r=tmpp_r-temph*top*bots*(1.0d0+bot/rold_r)
  594. c-------------------------------------------------------------
  595.        if(mmid .ge. 0.0d0) then
  596.          mpl_m = mmid
  597.          d2mr_l=dmir_l
  598.          d2mu_l=dmiu_l
  599.          d2mv_l=dmiv_l
  600.          d2mp_l=dmip_l
  601. c------------
  602.          d2mr_r=dmir_r
  603.          d2mu_r=dmiu_r
  604.          d2mv_r=dmiv_r
  605.          d2mp_r=dmip_r
  606. c------------
  607.        else
  608.          mpl_m = 0.0d0
  609.          d2mr_l=0.0d0
  610.          d2mu_l=0.0d0
  611.          d2mv_l=0.0d0
  612.          d2mp_l=0.0d0
  613. c------------
  614.          d2mr_r=0.0d0
  615.          d2mu_r=0.0d0
  616.          d2mv_r=0.0d0
  617.          d2mp_r=0.0d0
  618.        endif
  619. c---------------------------------------------
  620. cc       derivatives for the term termm
  621. cc------------------------------------------------------------------
  622.        if(mmid .le. 0.0d0) then
  623.          mmin_m = mmid
  624.          d3mr_l=dmir_l
  625.          d3mu_l=dmiu_l
  626.          d3mv_l=dmiv_l
  627.          d3mp_l=dmip_l
  628. c------------
  629.          d3mr_r=dmir_r
  630.          d3mu_r=dmiu_r
  631.          d3mv_r=dmiv_r
  632.          d3mp_r=dmip_r
  633. c------------
  634.        else
  635.          mmin_m = 0.0d0
  636.          d3mr_l=0.0d0
  637.          d3mu_l=0.0d0
  638.          d3mv_l=0.0d0
  639.          d3mp_l=0.0d0
  640. c------------
  641.          d3mr_r=0.0d0
  642.          d3mu_r=0.0d0
  643.          d3mv_r=0.0d0
  644.          d3mp_r=0.0d0
  645.        endif
  646. c---------------------------------------------------------------
  647. c Computing the calligraphic P+ and P- with their derivatives,
  648. c see (21a) & (21b) on p.368.
  649. c---------------------------------------------------------------
  650.        if(ml .ge. 1.0d0) then
  651.          Pplus = 1.0d0
  652.        else
  653.          if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  654.            Pplus=(ml+1.0d0)*(ml+1.0d0)*(2.0d0-ml)/4.0d0
  655.            Pplus=Pplus+alpha*ml*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  656.          else
  657.            Pplus = 0.0d0
  658.          endif
  659.        endif
  660. c---------------------------------------------------------------
  661.        if(mr .ge. 1.0d0) then
  662.          Pmin = 0.0d0
  663.        else
  664.          if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  665.            Pmin=(mr-1.0d0)*(mr-1.0d0)*(2.0d0+mr)/4.0d0
  666.            Pmin=Pmin-alpha*mr*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  667.          else
  668.            Pmin = 1.0d0
  669.          endif
  670.        endif
  671. c---------------------------------------------------------------
  672.        brac_l=(ml+1.0d0)*(2.0d0-ml)/2.0d0-(ml+1.0d0)*(ml+1.0d0)/4.0d0
  673.        brac_l=brac_l+alpha*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  674.        brac_l=brac_l+4.0d0*alpha*ml*ml*(ml*ml-1.0d0)
  675. c--------------
  676.        brac_r=(mr-1.0d0)*(2.0d0+mr)/2.0d0+(mr-1.0d0)*(mr-1.0d0)/4.0d0
  677.        brac_r=brac_r-alpha*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  678.        brac_r=brac_r-4.0d0*alpha*mr*mr*(mr*mr-1.0d0)
  679. c---------------------------------------------------------------
  680.        if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  681.          dPpr_l=brac_l*dmlr_l
  682.          dPpr_r=brac_l*dmlr_r
  683. c------------
  684.          dPpu_l=brac_l*dmlu_l
  685.          dPpu_r=brac_l*dmlu_r
  686. c------------
  687.          dPpv_l=brac_l*dmlv_l
  688.          dPpv_r=brac_l*dmlv_r
  689. c------------
  690.          dPpp_l=brac_l*dmlp_l
  691.          dPpp_r=brac_l*dmlp_r
  692. c------------
  693.        else
  694.          dPpr_l=0.0d0
  695.          dPpr_r=0.0d0
  696. c-----------
  697.          dPpu_l=0.0d0
  698.          dPpu_r=0.0d0
  699. c-----------
  700.          dPpv_l=0.0d0
  701.          dPpv_r=0.0d0
  702. c-----------
  703.          dPpp_l=0.0d0
  704.          dPpp_r=0.0d0
  705. c-----------
  706.        endif
  707. c---------------------------------------------------------------
  708.        if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  709.          dPmr_l=brac_r*dmrr_l
  710.          dPmr_r=brac_r*dmrr_r
  711. c------------
  712.          dPmu_l=brac_r*dmru_l
  713.          dPmu_r=brac_r*dmru_r
  714. c------------
  715.          dPmv_l=brac_r*dmrv_l
  716.          dPmv_r=brac_r*dmrv_r
  717. c------------
  718.          dPmp_l=brac_r*dmrp_l
  719.          dPmp_r=brac_r*dmrp_r
  720. c------------
  721.        else
  722.          dPmr_l=0.0d0
  723.          dPmr_r=0.0d0
  724. c-----------
  725.          dPmu_l=0.0d0
  726.          dPmu_r=0.0d0
  727. c-----------
  728.          dPmv_l=0.0d0
  729.          dPmv_r=0.0d0
  730. c-----------
  731.          dPmp_l=0.0d0
  732.          dPmp_r=0.0d0
  733. c-----------
  734.        endif
  735. c-------------------------------------------------------------------
  736. c computing pmid -- p_{1/2} and its derivatives, see (20b), p.367.
  737. c-------------------------------------------------------------------
  738.        dpir_l=dPpr_l*pold_l+dPmr_l*pold_r
  739.        dpiu_l=dPpu_l*pold_l+dPmu_l*pold_r
  740.        dpiv_l=dPpv_l*pold_l+dPmv_l*pold_r
  741.        dpip_l=dPpp_l*pold_l+Pplus+dPmp_l*pold_r
  742. c----------------------------
  743.        dpir_r=dPpr_r*pold_l+dPmr_r*pold_r
  744.        dpiu_r=dPpu_r*pold_l+dPmu_r*pold_r
  745.        dpiv_r=dPpv_r*pold_l+dPmv_r*pold_r
  746.        dpip_r=dPpp_r*pold_l+Pmin+dPmp_r*pold_r
  747. c---------------------------------------------------------------------
  748. c  Computation of the mass flux (rho * u)_1/2
  749. c---------------------------------------------------------------
  750.        rum=am*(mpl_m*rold_l+mmin_m*rold_r)+canc*am*termp
  751. c-------------------------------------------------------
  752.        rumr_l=damr_l*(mpl_m*rold_l+mmin_m*rold_r)+
  753.      &        am*(d2mr_l*rold_l+mpl_m+d3mr_l*rold_r)
  754.        rumr_l=rumr_l+canc*coef*(damr_l*termp+am*tmpr_l)
  755.        rumu_l=damu_l*(mpl_m*rold_l+mmin_m*rold_r)+
  756.      &        am*(d2mu_l*rold_l+d3mu_l*rold_r)
  757.        rumu_l=rumu_l+canc*coef*(damu_l*termp+am*tmpu_l)
  758.        rumv_l=damv_l*(mpl_m*rold_l+mmin_m*rold_r)+
  759.      &        am*(d2mv_l*rold_l+d3mv_l*rold_r)
  760.        rumv_l=rumv_l+canc*coef*(damv_l*termp+am*tmpv_l)
  761.        rump_l=damp_l*(mpl_m*rold_l+mmin_m*rold_r)+
  762.      &        am*(d2mp_l*rold_l+d3mp_l*rold_r)
  763.        rump_l=rump_l+canc*coef*(damp_l*termp+am*tmpp_l)
  764. c-------------------------------------------------
  765.        rumr_r=damr_r*(mpl_m*rold_l+mmin_m*rold_r)+
  766.      &        am*(d2mr_r*rold_l+mmin_m+d3mr_r*rold_r)
  767.        rumr_r=rumr_r+canc*coef*(damr_r*termp+am*tmpr_r)
  768.        rumu_r=damu_r*(mpl_m*rold_l+mmin_m*rold_r)+
  769.      &        am*(d2mu_r*rold_l+d3mu_r*rold_r)
  770.        rumu_r=rumu_r+canc*coef*(damu_r*termp+am*tmpu_r)
  771.        rumv_r=damv_r*(mpl_m*rold_l+mmin_m*rold_r)+
  772.      &        am*(d2mv_r*rold_l+d3mv_r*rold_r)
  773.        rumv_r=rumv_r+canc*coef*(damv_r*termp+am*tmpv_r)
  774.        rump_r=damp_r*(mpl_m*rold_l+mmin_m*rold_r)+
  775.      &        am*(d2mp_r*rold_l+d3mp_r*rold_r)
  776.        rump_r=rump_r+canc*coef*(damp_r*termp+am*tmpp_r)
  777. c---------------------------------------------------------------------
  778. c computing JACOBIAN as a derivative of the numerical flux function
  779. c with respect to the primitive variables.
  780. c Notation: jl(i,j) --- is the derivative of the i-component of the
  781. c           flux function with respect to the j-component of the
  782. c           vector of primitive variables of the left state.
  783. c           jr(i,j) --- is the derivative of the i-component of the
  784. c           flux function with respect to the j-component of the
  785. c           vector of primitive variables of the right state.
  786. c---------------------------------------------------------------------
  787.        jl(1,1)=0.0d0
  788.        jl(1,2)=0.0d0
  789.        jl(1,3)=0.0d0
  790.        jl(1,4)=0.0d0
  791. c------------------------------------
  792.        jr(1,1)=0.0d0
  793.        jr(1,2)=0.0d0
  794.        jr(1,3)=0.0d0
  795.        jr(1,4)=0.0d0
  796. c------------------------------------
  797. c---------------------------------------------------------
  798.        if(rum .ge. 0.0d0) then
  799.          br3=rumr_l*un_l
  800.          br4=rumr_l*ut_l
  801.        else
  802.          br3=rumr_l*un_r
  803.          br4=rumr_l*ut_r
  804.        endif
  805.        jl(2,1)=n_x*(br3+dpir_l)+br4*t_x
  806.        jl(3,1)=n_y*(br3+dpir_l)+br4*t_y
  807. c-------------------
  808.        if(rum .ge. 0.0d0) then
  809.          br3=rumu_l*un_l+rum*n_x
  810.          br4=rumu_l*ut_l+rum*t_x
  811.        else
  812.          br3=rumu_l*un_r
  813.          br4=rumu_l*ut_r
  814.        endif
  815.        jl(2,2)=n_x*(br3+dpiu_l)+br4*t_x
  816.        jl(3,2)=n_y*(br3+dpiu_l)+br4*t_y
  817. c-------------------
  818.        if(rum .ge. 0.0d0) then
  819.          br3=rumv_l*un_l+rum*n_y
  820.          br4=rumv_l*ut_l+rum*t_y
  821.        else
  822.          br3=rumv_l*un_r
  823.          br4=rumv_l*ut_r
  824.        endif
  825.        jl(2,3)=n_x*(br3+dpiv_l)+br4*t_x
  826.        jl(3,3)=n_y*(br3+dpiv_l)+br4*t_y
  827. c-------------------
  828.        if(rum .ge. 0.0d0) then
  829.          br3=rump_l*un_l
  830.          br4=rump_l*ut_l
  831.        else
  832.          br3=rump_l*un_r
  833.          br4=rump_l*ut_r
  834.        endif
  835.        jl(2,4)=n_x*(br3+dpip_l)+br4*t_x
  836.        jl(3,4)=n_y*(br3+dpip_l)+br4*t_y
  837. c-------------------------------------------------------------
  838. c-------------------------------------------------------------
  839.        if(rum .ge. 0.0d0) then
  840.          br3=rumr_r*un_l
  841.          br4=rumr_r*ut_l
  842.        else
  843.          br3=rumr_r*un_r
  844.          br4=rumr_r*ut_r
  845.        endif
  846.        jr(2,1)=n_x*(br3+dpir_r)+br4*t_x
  847.        jr(3,1)=n_y*(br3+dpir_r)+br4*t_y
  848. c-------------------
  849.        if(rum .ge. 0.0d0) then
  850.          br3=rumu_r*un_l
  851.          br4=rumu_r*ut_l
  852.        else
  853.          br3=rumu_r*un_r+rum*n_x
  854.          br4=rumu_r*ut_r+rum*t_x
  855.        endif
  856.        jr(2,2)=n_x*(br3+dpiu_r)+br4*t_x
  857.        jr(3,2)=n_y*(br3+dpiu_r)+br4*t_y
  858. c-------------------
  859.        if(rum .ge. 0.0d0) then
  860.          br3=rumv_r*un_l
  861.          br4=rumv_r*ut_l
  862.        else
  863.          br3=rumv_r*un_r+rum*n_y
  864.          br4=rumv_r*ut_r+rum*t_y
  865.        endif
  866.        jr(2,3)=n_x*(br3+dpiv_r)+br4*t_x
  867.        jr(3,3)=n_y*(br3+dpiv_r)+br4*t_y
  868. c-------------------
  869.        if(rum .ge. 0.0d0) then
  870.          br3=rump_r*un_l
  871.          br4=rump_r*ut_l
  872.        else
  873.          br3=rump_r*un_r
  874.          br4=rump_r*ut_r
  875.        endif
  876.        jr(2,4)=n_x*(br3+dpip_r)+br4*t_x
  877.        jr(3,4)=n_y*(br3+dpip_r)+br4*t_y
  878. c-------------------------------------------------------------
  879. c  ------  f44444444444444444444444444444 ---------
  880. c-------------------------------------------------------------
  881.        jl(4,1)=0.0d0
  882.        jl(4,2)=0.0d0
  883.        jl(4,3)=0.0d0
  884.        jl(4,4)=0.0d0
  885. c----------------------------------------------------------
  886.        jr(4,1)=0.0d0
  887.        jr(4,2)=0.0d0
  888.        jr(4,3)=0.0d0
  889.        jr(4,4)=0.0d0
  890. c-------------------------------------------------------------
  891. c---------------------------------------------------------------------
  892.        tcoef=nvect(1)*tvect(2)+tvect(1)*nvect(2)
  893.        bcoef=nvect(1)*tvect(2)-tvect(1)*nvect(2)
  894. c----------------------------------------------------------------------
  895.        do 11 i=1,4
  896.          jl(i,1)=jl(i,1)+jr(i,1)
  897.          jl(i,2)=jl(i,2)+jr(i,2)*(-tcoef/bcoef)+
  898.      &            jr(i,3)*2.0d0*nvect(1)*tvect(1)/bcoef
  899.          jl(i,3)=jl(i,3)+jr(i,2)*(-2.0d0*nvect(2)*tvect(2)/bcoef)+
  900.      &            jr(i,3)*tcoef/bcoef
  901.          jl(i,4)=jl(i,4)+jr(i,4)
  902.  11    continue
  903. c----------------------------------------------------------------------
  904.        return
  905.        end
  906.  
  907.  
  908.  
  909.  
  910.  
  911.  
  912.  
  913.  
  914.  
  915.  
  916.  
  917.  
  918.  
  919.  
  920.  
  921.  
  922.  
  923.  
  924.  
  925.  
  926.  
  927.  
  928.  
  929.  
  930.  
  931.  
  932.  
  933.  
  934.  
  935.  
  936.  
  937.  
  938.  
  939.  

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