Télécharger la plaquette Cast3M


fr en

Télécharger conjp4.eso

Retour à la liste

Numérotation des lignes : 

conjp4

  1. C CONJP4    SOURCE    CB215821  16/04/21    21:15:59     8920
  2. c Subroutine appellé par konjp4.eso
  3. c                        cli153.eso
  4. c----------------------------------------------------------------------
  5. c GENERAL DESCRIPTION:
  6. c   This subroutine provides the jacobians which are the derivatives
  7. c   of the numerical flux function defined at the cell interface
  8. c   with respect to the primitive variables of the left and right
  9. c   cells (relative to the cell interface).
  10. c
  11. c EQUATIONS: 3D Euler equations of gas dynamics
  12. c
  13. c
  14. c REFERENCE: JCP, 129, 364-382 (1996)
  15. c            " A Sequel to AUSM: AUSM+ ".
  16. c----------------------------------------------------------------------
  17. c INPUT:
  18. c
  19. c      alpha   -- parameter of the AUSM+ scheme in the Pressure function;
  20. c                 ( -3/4 <= alpha <= 3/16 ) (imposed as a parameter)
  21. c
  22. c      beta    -- parameter of the AUSM+ scheme in the Mach function;
  23. c                 ( -1/16 <= beta <= 1/2 )  (imposed as a parameter)
  24. c
  25. c      wvec_l  -- vector of the primitive variables (rho,ux,uy,uz,p) at the
  26. c                 left cell;
  27. c
  28. c      wvec_r  -- vector of the primitive variables (rho,ux,uy,uz,p) at the
  29. c                 right cell;
  30. c
  31. c      nvect   -- normal vector to the interface (3 components in 3D);
  32. c
  33. c      tvec1   -- tangential vector to the interface;
  34. c
  35. c      tvec2   -- tangential vector to the interface;
  36. c
  37. c      ga      -- ratio of the specific heats (assumed constant)
  38. c----------------------------------------------------------------------
  39. c
  40. c OUTPUT:
  41. c
  42. c      jl     -- jakobian matrix 5 by 5 - derivatives of the numerical
  43. c                 flux function with respect to the primitive variables
  44. c                 from the left cell;
  45. c
  46. c      jr     -- jakobian matrix 5 by 5 - derivatives of the numerical
  47. c                 flux function with respect to the primitive variables
  48. c                 from the right cell.
  49. c
  50. c N.B. The numerical flux is expressed in the global frame
  51. c      (x,y,z); i.e.
  52. c
  53. c      \rho un
  54. c      \rho un ux + p (n.x)
  55. c      \rho un uy + p (n.y)
  56. c      \rho un uz + p (n.z)
  57. c      \rho un ht
  58. c
  59. c----------------------------------------------------------------------
  60.        SUBROUTINE CONJP4(JL,JR,WVEC_L,WVEC_R,NVECT,TVEC1,TVEC2,
  61.      &                   ga)
  62.       IMPLICIT INTEGER(I-N)
  63.        real*8 wvec_l(5),wvec_r(5)
  64.        real*8 nvect(3),tvec1(3),tvec2(3)
  65.        real*8 jl(5,5),jr(5,5),f(5)
  66.        real*8 alpha,beta
  67.        real*8 ga, gm1,temph
  68.        real*8 n_x,n_y,n_z,t1_x,t1_y,t1_z,t2_x,t2_y,t2_z
  69.        real*8 un_l,un_r,ut1_l,ut1_r,ut2_l,ut2_r
  70.        real*8 ml,mr,Mplus,Mmin,mmid
  71.        real*8 mpl_m, mmin_m,am
  72.        real*8 rold_l,uold_l,vold_l,wold_l,pold_l,eold_l
  73.        real*8 rold_r,uold_r,vold_r,wold_r,pold_r,eold_r
  74.        real*8 Pplus,Pmin,pmid
  75.        real*8 hr_l,hr_r,det,c11,c12,c13,c21,c22,c23,c31,c32,c33
  76.        real*8 br1,temp_l,temp_r,brac_l,brac_r
  77.        real*8 aleft, arigh
  78.        real*8 damr_l,damr_r,damu_l,damu_r
  79.        real*8 damv_l,damv_r,damp_l,damp_r
  80.        real*8 damw_l,damw_r
  81.        real*8 dmlr_l,dmlr_r,dmlu_l,dmlu_r
  82.        real*8 dmlv_l,dmlv_r,dmlp_l,dmlp_r
  83.        real*8 dmlw_l,dmlw_r
  84.        real*8 dmrr_l,dmrr_r,dmru_l,dmru_r
  85.        real*8 dmrv_l,dmrv_r,dmrp_l,dmrp_r
  86.        real*8 dmrw_l,dmrw_r
  87.        real*8 dMpr_l,dMpr_r,dMpu_l,dMpu_r
  88.        real*8 dMpv_l,dMpv_r,dMpp_l,dMpp_r
  89.        real*8 dMpw_l,dMpw_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 dMmw_l,dMmw_r
  93.        real*8 dmir_l,dmir_r,dmiu_l,dmiu_r
  94.        real*8 dmiv_l,dmiv_r,dmip_l,dmip_r
  95.        real*8 dmiw_l,dmiw_r
  96.        real*8 d3mr_l,d3mr_r,d3mu_l,d3mu_r
  97.        real*8 d3mv_l,d3mv_r,d3mp_l,d3mp_r
  98.        real*8 d3mw_l,d3mw_r
  99.        real*8 d2mr_l,d2mr_r,d2mu_l,d2mu_r
  100.        real*8 d2mv_l,d2mv_r,d2mp_l,d2mp_r
  101.        real*8 d2mw_l,d2mw_r
  102.        real*8 dPpr_l,dPpr_r,dPpu_l,dPpu_r
  103.        real*8 dPpv_l,dPpv_r,dPpp_l,dPpp_r
  104.        real*8 dPpw_l,dPpw_r
  105.        real*8 dPmr_l,dPmr_r,dPmu_l,dPmu_r
  106.        real*8 dPmv_l,dPmv_r,dPmp_l,dPmp_r
  107.        real*8 dPmw_l,dPmw_r
  108.        real*8 dpir_l,dpir_r,dpiu_l,dpiu_r
  109.        real*8 dpiv_l,dpiv_r,dpip_l,dpip_r
  110.        real*8 dpiw_l,dpiw_r
  111.        real*8 ff,fs,ft,dffr,dffu,dffv,dffw,dffp
  112.        real*8 dfsr,dfsu,dfsv,dfsw,dfsp,dftr,dftu,dftv,dftw,dftp
  113.        parameter(alpha = 0.1875D0, beta = 0.125D0)
  114. c-----------------------
  115.        gm1=ga-1.0d0
  116. c-----------------------
  117.        n_x=nvect(1)
  118.        n_y=nvect(2)
  119.        n_z=nvect(3)
  120. c-------------------
  121.        t1_x=tvec1(1)
  122.        t1_y=tvec1(2)
  123.        t1_z=tvec1(3)
  124. c-------------------
  125.        t2_x=tvec2(1)
  126.        t2_y=tvec2(2)
  127.        t2_z=tvec2(3)
  128. c-------------------------------------------
  129. c These are coefficients for the expressions
  130. c
  131. c  un  = (c11*ux + c12*uy + c13*uz)/det
  132. c  ut1 = (c21*ux + c22*uy + c23*uz)/det
  133. c  ut2 = (c31*ux + c32*uy + c33*uz)/det
  134. c-------------------------------------------
  135.        c11=t1_y*t2_z - t1_z*t2_y
  136.        c12=t1_z*t2_x - t1_x*t2_z
  137.        c13=t1_x*t2_y - t1_y*t2_x
  138.        det=n_x*c11 + n_y*c12 + n_z*c13
  139. c--------------------------------------
  140.        c21=n_z*t2_y - n_y*t2_z
  141.        c22=n_x*t2_z - n_z*t2_x
  142.        c23=n_y*t2_x - n_x*t2_y
  143. c--------------------------------------
  144.        c31=n_y*t1_z - n_z*t1_y
  145.        c32=n_z*t1_x - n_x*t1_z
  146.        c33=n_x*t1_y - n_y*t1_x
  147. c--------------------------------------
  148.        rold_l=wvec_l(1)
  149.        uold_l=wvec_l(2)
  150.        vold_l=wvec_l(3)
  151.        wold_l=wvec_l(4)
  152.        pold_l=wvec_l(5)
  153. c-----------------------
  154.        rold_r=wvec_r(1)
  155.        uold_r=wvec_r(2)
  156.        vold_r=wvec_r(3)
  157.        wold_r=wvec_r(4)
  158.        pold_r=wvec_r(5)
  159. c------------------------------------------------------------------
  160. c Computation of the specific total energy on the left and right.
  161. c------------------------------------------------------------------
  162.        eold_l=(uold_l*uold_l+vold_l*vold_l+wold_l*wold_l)/2.0d0
  163.        eold_l=eold_l+pold_l/(gm1*rold_l)
  164.        eold_r=(uold_r*uold_r+vold_r*vold_r+wold_r*wold_r)/2.0d0
  165.        eold_r=eold_r+pold_r/(gm1*rold_r)
  166. c-------------------------------------------------------------------
  167. c Computation of the speed of sound and its derivatives;
  168. c numerical speed of sound at the interface is taken as an average
  169. c of the speeds of sounds of the neighbouring cells
  170. c---------------------------------------------------------------------
  171.        aleft=SQRT(ga*pold_l/rold_l)
  172.        arigh=SQRT(ga*pold_r/rold_r)
  173.        am=0.5d0*(aleft+arigh)
  174. c--------------------------------------------------------------------
  175.        damr_r=-arigh/(4.0d0*rold_r)
  176.        damu_r=0.0d0
  177.        damv_r=0.0d0
  178.        damw_r=0.0d0
  179.        damp_r=ga/(4.0d0*arigh*rold_r)
  180. c-----------------------
  181.        damr_l=-aleft/(4.0d0*rold_l)
  182.        damu_l=0.0d0
  183.        damv_l=0.0d0
  184.        damw_l=0.0d0
  185.        damp_l=ga/(4.0d0*aleft*rold_l)
  186. c-------------------------------------------------------------------
  187. c Computing numerical Mach number and its derivatives,
  188. c see p.370, under (A1).
  189. c-------------------------------------------------------------------
  190.        un_l=uold_l*n_x+vold_l*n_y+wold_l*n_z
  191.        un_r=uold_r*n_x+vold_r*n_y+wold_r*n_z
  192. c----------
  193.        ut1_l=uold_l*t1_x+vold_l*t1_y+wold_l*t1_z
  194.        ut1_r=uold_r*t1_x+vold_r*t1_y+wold_r*t1_z
  195. c----------
  196.        ut2_l=uold_l*t2_x+vold_l*t2_y+wold_l*t2_z
  197.        ut2_r=uold_r*t2_x+vold_r*t2_y+wold_r*t2_z
  198. c----------------------------------------
  199.        ml=un_l/am
  200.        mr=un_r/am
  201. c-------------------------------------------------------------------
  202. c  Mplus and Mmin are calligraphic lettes M+ and M- from the paper,
  203. c  see (19a) and (19b), p.367.
  204. c-------------------------------------------------------------------
  205.        if(ABS(ml) .ge. 1.0d0) then
  206.          Mplus=(ml+ABS(ml))/2.0d0
  207.        else
  208.          Mplus=(ml+1.0d0)*(ml+1.0d0)/4.0d0
  209.          Mplus=Mplus+beta*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  210.        endif
  211. c-----------------
  212.        if(ABS(mr) .ge. 1.0d0) then
  213.          Mmin=(mr-ABS(mr))/2.0d0
  214.        else
  215.          Mmin=-(mr-1.0d0)*(mr-1.0d0)/4.0d0
  216.          Mmin=Mmin-beta*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  217.        endif
  218. c-------------------------------------------------------------------
  219. c Derivatives of ml and mr with respect to both sets of primitive
  220. c variables: left  and right.
  221. c-------------------------------------------------------------------
  222.        temp_l=-un_l/(am*am)
  223.        temp_r=-un_r/(am*am)
  224. c--------
  225.        dmlr_l=temp_l*damr_l
  226.        dmlr_r=temp_l*damr_r
  227.        dmrr_l=temp_r*damr_l
  228.        dmrr_r=temp_r*damr_r
  229. c--------
  230.        dmlu_l=n_x/am+temp_l*damu_l
  231.        dmlu_r=temp_l*damu_r
  232.        dmru_l=temp_r*damu_l
  233.        dmru_r=n_x/am+temp_r*damu_r
  234. c--------
  235.        dmlv_l=n_y/am+temp_l*damv_l
  236.        dmlv_r=temp_l*damv_r
  237.        dmrv_l=temp_r*damv_l
  238.        dmrv_r=n_y/am+temp_r*damv_r
  239. c--------
  240.        dmlw_l=n_z/am+temp_l*damw_l
  241.        dmlw_r=temp_l*damw_r
  242.        dmrw_l=temp_r*damw_l
  243.        dmrw_r=n_z/am+temp_r*damw_r
  244. c--------
  245.        dmlp_l=temp_l*damp_l
  246.        dmlp_r=temp_l*damp_r
  247.        dmrp_l=temp_r*damp_l
  248.        dmrp_r=temp_r*damp_r
  249. c-----------------------------------------------------------
  250. c mmid is m_{1/2} (notation as in the paper)
  251. c----------------------------------------------------------
  252.        mmid=Mplus+Mmin
  253. c-----------------------------------------------------------
  254. c Computing the derivatives of M+ and M-
  255. c-----------------------------------------------------------
  256.        if(ml .ge. 1.0d0) then
  257.          dMpr_l=dmlr_l
  258.          dMpu_l=dmlu_l
  259.          dMpv_l=dmlv_l
  260.          dMpw_l=dmlw_l
  261.          dMpp_l=dmlp_l
  262. c--------------------
  263.          dMpr_r=dmlr_r
  264.          dMpu_r=dmlu_r
  265.          dMpv_r=dmlv_r
  266.          dMpw_r=dmlw_r
  267.          dMpp_r=dmlp_r
  268.        else
  269.         if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  270.          temph=(ml+1.0d0)/2.0d0
  271.          dMpr_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_l
  272.          dMpu_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_l
  273.          dMpv_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_l
  274.          dMpw_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlw_l
  275.          dMpp_l=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_l
  276. c--------------------
  277.            dMpr_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlr_r
  278.            dMpu_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlu_r
  279.            dMpv_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlv_r
  280.            dMpw_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlw_r
  281.            dMpp_r=(temph+4.0d0*beta*ml*(ml*ml-1.0d0))*dmlp_r
  282.          else
  283.            dMpr_l=0.0d0
  284.            dMpu_l=0.0d0
  285.            dMpv_l=0.0d0
  286.            dMpw_l=0.0d0
  287.            dMpp_l=0.0d0
  288. c---------------------
  289.            dMpr_r=0.0d0
  290.            dMpu_r=0.0d0
  291.            dMpv_r=0.0d0
  292.            dMpw_r=0.0d0
  293.            dMpp_r=0.0d0
  294.          endif
  295.        endif
  296. c-----------------------------------------------------------
  297.        if(mr .ge. 1.0d0) then
  298.          dMmr_l=0.0d0
  299.          dMmu_l=0.0d0
  300.          dMmv_l=0.0d0
  301.          dMmw_l=0.0d0
  302.          dMmp_l=0.0d0
  303. c---------------------
  304.          dMmr_r=0.0d0
  305.          dMmu_r=0.0d0
  306.          dMmv_r=0.0d0
  307.          dMmw_r=0.0d0
  308.          dMmp_r=0.0d0
  309.        else
  310.          if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  311.            temph=(-mr+1.0d0)/2.0d0
  312.            dMmr_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_l
  313.            dMmu_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_l
  314.            dMmv_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_l
  315.            dMmw_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrw_l
  316.            dMmp_l=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_l
  317. c--------------------
  318.            dMmr_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrr_r
  319.            dMmu_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmru_r
  320.            dMmv_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrv_r
  321.            dMmw_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrw_r
  322.            dMmp_r=(temph-4.0d0*beta*mr*(mr*mr-1.0d0))*dmrp_r
  323.          else
  324.            dMmr_l=dmrr_l
  325.            dMmu_l=dmru_l
  326.            dMmv_l=dmrv_l
  327.            dMmw_l=dmrw_l
  328.            dMmp_l=dmrp_l
  329. c--------------------
  330.            dMmr_r=dmrr_r
  331.            dMmu_r=dmru_r
  332.            dMmv_r=dmrv_r
  333.            dMmw_r=dmrw_r
  334.            dMmp_r=dmrp_r
  335.          endif
  336.        endif
  337. c-----------------------------------------------------------------
  338. c computing the derivatives of m_{1/2} (notation as in the paper)
  339. c-----------------------------------------------------------------
  340.        dmir_l=dMpr_l+dMmr_l
  341.        dmir_r=dMpr_r+dMmr_r
  342. c-------------
  343.        dmiu_l=dMpu_l+dMmu_l
  344.        dmiu_r=dMpu_r+dMmu_r
  345. c-------------
  346.        dmiv_l=dMpv_l+dMmv_l
  347.        dmiv_r=dMpv_r+dMmv_r
  348. c-------------
  349.        dmiw_l=dMpw_l+dMmw_l
  350.        dmiw_r=dMpw_r+dMmw_r
  351. c-------------
  352.        dmip_l=dMpp_l+dMmp_l
  353.        dmip_r=dMpp_r+dMmp_r
  354. c-------------------------------------------------------------
  355. c computing the main convective variables and their derivatives
  356. c mpl_m is m^{+}_{1/2} (paper's notation) and
  357. c mmin_m  is m^{-}_{1/2} (paper's notation)
  358. c-------------------------------------------------------------
  359.        if(mmid .gt. 0.0d0) then
  360.          mpl_m = mmid
  361.          d2mr_l=dmir_l
  362.          d2mu_l=dmiu_l
  363.          d2mv_l=dmiv_l
  364.          d2mw_l=dmiw_l
  365.          d2mp_l=dmip_l
  366. c------------
  367.          d2mr_r=dmir_r
  368.          d2mu_r=dmiu_r
  369.          d2mv_r=dmiv_r
  370.          d2mw_r=dmiw_r
  371.          d2mp_r=dmip_r
  372. c------------
  373.        else
  374.          mpl_m = 0.0d0
  375.          d2mr_l=0.0d0
  376.          d2mu_l=0.0d0
  377.          d2mv_l=0.0d0
  378.          d2mw_l=0.0d0
  379.          d2mp_l=0.0d0
  380. c------------
  381.          d2mr_r=0.0d0
  382.          d2mu_r=0.0d0
  383.          d2mv_r=0.0d0
  384.          d2mw_r=0.0d0
  385.          d2mp_r=0.0d0
  386.        endif
  387. c---------------------------------------------
  388.        if(mmid .lt. 0.0d0) then
  389.          mmin_m = mmid
  390.          d3mr_l=dmir_l
  391.          d3mu_l=dmiu_l
  392.          d3mv_l=dmiv_l
  393.          d3mw_l=dmiw_l
  394.          d3mp_l=dmip_l
  395. c------------
  396.          d3mr_r=dmir_r
  397.          d3mu_r=dmiu_r
  398.          d3mv_r=dmiv_r
  399.          d3mw_r=dmiw_r
  400.          d3mp_r=dmip_r
  401. c------------
  402.        else
  403.          mmin_m = 0.0d0
  404.          d3mr_l=0.0d0
  405.          d3mu_l=0.0d0
  406.          d3mv_l=0.0d0
  407.          d3mw_l=0.0d0
  408.          d3mp_l=0.0d0
  409. c------------
  410.          d3mr_r=0.0d0
  411.          d3mu_r=0.0d0
  412.          d3mv_r=0.0d0
  413.          d3mw_r=0.0d0
  414.          d3mp_r=0.0d0
  415.        endif
  416. c---------------------------------------------------------------
  417. c Computing the calligraphic P+ and P- with their derivatives
  418. c see (21a) & (21b) on p.368.
  419. c---------------------------------------------------------------
  420.        if(ml .ge. 1.0d0) then
  421.          Pplus = 1.0d0
  422.        else
  423.          if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  424.            Pplus=(ml+1.0d0)*(ml+1.0d0)*(2.0d0-ml)/4.0d0
  425.            Pplus=Pplus+alpha*ml*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  426.          else
  427.            Pplus = 0.0d0
  428.          endif
  429.        endif
  430. c---------------------------------------------------------------
  431.        if(mr .ge. 1.0d0) then
  432.          Pmin = 0.0d0
  433.        else
  434.          if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  435.            Pmin=(mr-1.0d0)*(mr-1.0d0)*(2.0d0+mr)/4.0d0
  436.            Pmin=Pmin-alpha*mr*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  437.          else
  438.            Pmin = 1.0d0
  439.          endif
  440.        endif
  441. c---------------------------------------------------------------
  442.        brac_l=(ml+1.0d0)*(2.0d0-ml)/2.0d0-(ml+1.0d0)*(ml+1.0d0)/4.0d0
  443.        brac_l=brac_l+alpha*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  444.        brac_l=brac_l+4.0d0*alpha*ml*ml*(ml*ml-1.0d0)
  445. c--------------
  446.        brac_r=(mr-1.0d0)*(2.0d0+mr)/2.0d0+(mr-1.0d0)*(mr-1.0d0)/4.0d0
  447.        brac_r=brac_r-alpha*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  448.        brac_r=brac_r-4.0d0*alpha*mr*mr*(mr*mr-1.0d0)
  449. c---------------------------------------------------------------
  450.        if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  451.          dPpr_l=brac_l*dmlr_l
  452.          dPpr_r=brac_l*dmlr_r
  453. c------------
  454.          dPpu_l=brac_l*dmlu_l
  455.          dPpu_r=brac_l*dmlu_r
  456. c------------
  457.          dPpv_l=brac_l*dmlv_l
  458.          dPpv_r=brac_l*dmlv_r
  459. c------------
  460.          dPpw_l=brac_l*dmlw_l
  461.          dPpw_r=brac_l*dmlw_r
  462. c------------
  463.          dPpp_l=brac_l*dmlp_l
  464.          dPpp_r=brac_l*dmlp_r
  465. c------------
  466.        else
  467.          dPpr_l=0.0d0
  468.          dPpr_r=0.0d0
  469. c-----------
  470.          dPpu_l=0.0d0
  471.          dPpu_r=0.0d0
  472. c-----------
  473.          dPpv_l=0.0d0
  474.          dPpv_r=0.0d0
  475. c-----------
  476.          dPpw_l=0.0d0
  477.          dPpw_r=0.0d0
  478. c-----------
  479.          dPpp_l=0.0d0
  480.          dPpp_r=0.0d0
  481. c-----------
  482.        endif
  483. c---------------------------------------------------------------
  484.        if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  485.          dPmr_l=brac_r*dmrr_l
  486.          dPmr_r=brac_r*dmrr_r
  487. c------------
  488.          dPmu_l=brac_r*dmru_l
  489.          dPmu_r=brac_r*dmru_r
  490. c------------
  491.          dPmv_l=brac_r*dmrv_l
  492.          dPmv_r=brac_r*dmrv_r
  493. c------------
  494.          dPmw_l=brac_r*dmrw_l
  495.          dPmw_r=brac_r*dmrw_r
  496. c------------
  497.          dPmp_l=brac_r*dmrp_l
  498.          dPmp_r=brac_r*dmrp_r
  499. c------------
  500.        else
  501.          dPmr_l=0.0d0
  502.          dPmr_r=0.0d0
  503. c-----------
  504.          dPmu_l=0.0d0
  505.          dPmu_r=0.0d0
  506. c-----------
  507.          dPmv_l=0.0d0
  508.          dPmv_r=0.0d0
  509. c-----------
  510.          dPmw_l=0.0d0
  511.          dPmw_r=0.0d0
  512. c-----------
  513.          dPmp_l=0.0d0
  514.          dPmp_r=0.0d0
  515. c-----------
  516.        endif
  517. c---------------------------------------------------------------------
  518. c computing pmid -- p_{1/2} and its derivatives, see (20b), p.367.
  519. c---------------------------------------------------------------------
  520.        pmid=Pplus*pold_l + Pmin*pold_r
  521.        dpir_l=dPpr_l*pold_l+dPmr_l*pold_r
  522.        dpiu_l=dPpu_l*pold_l+dPmu_l*pold_r
  523.        dpiv_l=dPpv_l*pold_l+dPmv_l*pold_r
  524.        dpiw_l=dPpw_l*pold_l+dPmw_l*pold_r
  525.        dpip_l=dPpp_l*pold_l+Pplus+dPmp_l*pold_r
  526. c----------------------------
  527.        dpir_r=dPpr_r*pold_l+dPmr_r*pold_r
  528.        dpiu_r=dPpu_r*pold_l+dPmu_r*pold_r
  529.        dpiv_r=dPpv_r*pold_l+dPmv_r*pold_r
  530.        dpiw_r=dPpw_r*pold_l+dPmw_r*pold_r
  531.        dpip_r=dPpp_r*pold_l+Pmin+dPmp_r*pold_r
  532. c---------------------------------------------------------------------
  533. c computing JACOBIAN as a derivative of the numerical flux function
  534. c with respect to the primitive variables.
  535. c Notation: jl(i,j) --- is the derivative of the i-component of the
  536. c           flux function with respect to the j-component of the
  537. c           vector of primitive variables of the left state.
  538. c           jr(i,j) --- is the derivative of the i-component of the
  539. c           flux function with respect to the j-component of the
  540. c           vector of primitive variables of the right state.
  541. c---------------------------------------------------------------------
  542.        f(1)=am*(mpl_m*rold_l+mmin_m*rold_r)
  543. c---------------------------------------------------------------------
  544.        jl(1,1)=damr_l*f(1)/am+am*(d2mr_l*rold_l+mpl_m)
  545.        jl(1,1)=jl(1,1)+am*d3mr_l*rold_r
  546.        jl(1,2)=damu_l*f(1)/am+am*(d2mu_l*rold_l+d3mu_l*rold_r)
  547.        jl(1,3)=damv_l*f(1)/am+am*(d2mv_l*rold_l+d3mv_l*rold_r)
  548.        jl(1,4)=damw_l*f(1)/am+am*(d2mw_l*rold_l+d3mw_l*rold_r)
  549.        jl(1,5)=damp_l*f(1)/am+am*(d2mp_l*rold_l+d3mp_l*rold_r)
  550. c------------------------------------
  551.        jr(1,1)=damr_r*f(1)/am+am*(d2mr_r*rold_l+mmin_m)
  552.        jr(1,1)=jr(1,1)+am*d3mr_r*rold_r
  553.        jr(1,2)=damu_r*f(1)/am+am*(d2mu_r*rold_l+d3mu_r*rold_r)
  554.        jr(1,3)=damv_r*f(1)/am+am*(d2mv_r*rold_l+d3mv_r*rold_r)
  555.        jr(1,4)=damw_r*f(1)/am+am*(d2mw_r*rold_l+d3mw_r*rold_r)
  556.        jr(1,5)=damp_r*f(1)/am+am*(d2mp_r*rold_l+d3mp_r*rold_r)
  557. c------------------------------------
  558. c------------------------------------
  559.        ff=mpl_m*rold_l*un_l+mmin_m*rold_r*un_r
  560.        dffr=d2mr_l*rold_l*un_l+mpl_m*un_l+d3mr_l*rold_r*un_r
  561.        dffu=d2mu_l*rold_l*un_l+mpl_m*rold_l*c11/det
  562.        dffu=dffu+d3mu_l*rold_r*un_r
  563.        dffv=d2mv_l*rold_l*un_l+mpl_m*rold_l*c12/det
  564.        dffv=dffv+d3mv_l*rold_r*un_r
  565.        dffw=d2mw_l*rold_l*un_l+mpl_m*rold_l*c13/det
  566.        dffw=dffw+d3mw_l*rold_r*un_r
  567.        dffp=d2mp_l*rold_l*un_l+d3mp_l*rold_r*un_r
  568. c------------------------------------------------
  569.        fs=mpl_m*rold_l*ut1_l+mmin_m*rold_r*ut1_r
  570.        dfsr=d2mr_l*rold_l*ut1_l+mpl_m*ut1_l
  571.        dfsr=dfsr+d3mr_l*rold_r*ut1_r
  572.        dfsu=d2mu_l*rold_l*ut1_l+mpl_m*rold_l*c21/det
  573.        dfsu=dfsu+d3mu_l*rold_r*ut1_r
  574.        dfsv=d2mv_l*rold_l*ut1_l+mpl_m*rold_l*c22/det
  575.        dfsv=dfsv+d3mv_l*rold_r*ut1_r
  576.        dfsw=d2mw_l*rold_l*ut1_l+mpl_m*rold_l*c23/det
  577.        dfsw=dfsw+d3mw_l*rold_r*ut1_r
  578.        dfsp=d2mp_l*rold_l*ut1_l+d3mp_l*rold_r*ut1_r
  579. c-------------------------------------------------
  580.        ft=mpl_m*rold_l*ut2_l+mmin_m*rold_r*ut2_r
  581.        dftr=d2mr_l*rold_l*ut2_l+mpl_m*ut2_l
  582.        dftr=dftr+d3mr_l*rold_r*ut2_r
  583.        dftu=d2mu_l*rold_l*ut2_l+mpl_m*rold_l*c31/det
  584.        dftu=dftu+d3mu_l*rold_r*ut2_r
  585.        dftv=d2mv_l*rold_l*ut2_l+mpl_m*rold_l*c32/det
  586.        dftv=dftv+d3mv_l*rold_r*ut2_r
  587.        dftw=d2mw_l*rold_l*ut2_l+mpl_m*rold_l*c33/det
  588.        dftw=dftw+d3mw_l*rold_r*ut2_r
  589.        dftp=d2mp_l*rold_l*ut2_l+d3mp_l*rold_r*ut2_r
  590. c-------------------------------------------------
  591.        f(2)=am*(ff*n_x+fs*t1_x+ft*t2_x)+pmid*n_x
  592. c-------------------------------------------------
  593.        jl(2,1)=damr_l*(f(2)-pmid*n_x)/am+dpir_l*n_x
  594.        jl(2,1)=jl(2,1)+am*(dffr*n_x+dfsr*t1_x+dftr*t2_x)
  595.        jl(2,2)=damu_l*(f(2)-pmid*n_x)/am+dpiu_l*n_x
  596.        jl(2,2)=jl(2,2)+am*(dffu*n_x+dfsu*t1_x+dftu*t2_x)
  597.        jl(2,3)=damv_l*(f(2)-pmid*n_x)/am+dpiv_l*n_x
  598.        jl(2,3)=jl(2,3)+am*(dffv*n_x+dfsv*t1_x+dftv*t2_x)
  599.        jl(2,4)=damw_l*(f(2)-pmid*n_x)/am+dpiw_l*n_x
  600.        jl(2,4)=jl(2,4)+am*(dffw*n_x+dfsw*t1_x+dftw*t2_x)
  601.        jl(2,5)=damp_l*(f(2)-pmid*n_x)/am+dpip_l*n_x
  602.        jl(2,5)=jl(2,5)+am*(dffp*n_x+dfsp*t1_x+dftp*t2_x)
  603. c-------------------------------------------------
  604.        f(3)=am*(ff*n_y+fs*t1_y+ft*t2_y)+pmid*n_y
  605. c-------------------------------------------------
  606.        jl(3,1)=damr_l*(f(3)-pmid*n_y)/am+dpir_l*n_y
  607.        jl(3,1)=jl(3,1)+am*(dffr*n_y+dfsr*t1_y+dftr*t2_y)
  608.        jl(3,2)=damu_l*(f(3)-pmid*n_y)/am+dpiu_l*n_y
  609.        jl(3,2)=jl(3,2)+am*(dffu*n_y+dfsu*t1_y+dftu*t2_y)
  610.        jl(3,3)=damv_l*(f(3)-pmid*n_y)/am+dpiv_l*n_y
  611.        jl(3,3)=jl(3,3)+am*(dffv*n_y+dfsv*t1_y+dftv*t2_y)
  612.        jl(3,4)=damw_l*(f(3)-pmid*n_y)/am+dpiw_l*n_y
  613.        jl(3,4)=jl(3,4)+am*(dffw*n_y+dfsw*t1_y+dftw*t2_y)
  614.        jl(3,5)=damp_l*(f(3)-pmid*n_y)/am+dpip_l*n_y
  615.        jl(3,5)=jl(3,5)+am*(dffp*n_y+dfsp*t1_y+dftp*t2_y)
  616. c-------------------------------------------------
  617.        f(4)=am*(ff*n_z+fs*t1_z+ft*t2_z)+pmid*n_z
  618. c-------------------------------------------------
  619.        jl(4,1)=damr_l*(f(4)-pmid*n_z)/am+dpir_l*n_z
  620.        jl(4,1)=jl(4,1)+am*(dffr*n_z+dfsr*t1_z+dftr*t2_z)
  621.        jl(4,2)=damu_l*(f(4)-pmid*n_z)/am+dpiu_l*n_z
  622.        jl(4,2)=jl(4,2)+am*(dffu*n_z+dfsu*t1_z+dftu*t2_z)
  623.        jl(4,3)=damv_l*(f(4)-pmid*n_z)/am+dpiv_l*n_z
  624.        jl(4,3)=jl(4,3)+am*(dffv*n_z+dfsv*t1_z+dftv*t2_z)
  625.        jl(4,4)=damw_l*(f(4)-pmid*n_z)/am+dpiw_l*n_z
  626.        jl(4,4)=jl(4,4)+am*(dffw*n_z+dfsw*t1_z+dftw*t2_z)
  627.        jl(4,5)=damp_l*(f(4)-pmid*n_z)/am+dpip_l*n_z
  628.        jl(4,5)=jl(4,5)+am*(dffp*n_z+dfsp*t1_z+dftp*t2_z)
  629. c-------------------------------------------------------
  630. c derivatives with respect to the right
  631. c set of primitive variables
  632. c-------------------------------------------------------
  633.        dffr=d2mr_r*rold_l*un_l+mmin_m*un_r+d3mr_r*rold_r*un_r
  634.        dffu=d2mu_r*rold_l*un_l+mmin_m*rold_r*c11/det
  635.        dffu=dffu+d3mu_r*rold_r*un_r
  636.        dffv=d2mv_r*rold_l*un_l+mmin_m*rold_r*c12/det
  637.        dffv=dffv+d3mv_r*rold_r*un_r
  638.        dffw=d2mw_r*rold_l*un_l+mmin_m*rold_r*c13/det
  639.        dffw=dffw+d3mw_r*rold_r*un_r
  640.        dffp=d2mp_r*rold_l*un_l+d3mp_r*rold_r*un_r
  641. c------------------------------------------------------
  642.        dfsr=d2mr_r*rold_l*ut1_l+mmin_m*ut1_r
  643.        dfsr=dfsr+d3mr_r*rold_r*ut1_r
  644.        dfsu=d2mu_r*rold_l*ut1_l+mmin_m*rold_r*c21/det
  645.        dfsu=dfsu+d3mu_r*rold_r*ut1_r
  646.        dfsv=d2mv_r*rold_l*ut1_l+mmin_m*rold_r*c22/det
  647.        dfsv=dfsv+d3mv_r*rold_r*ut1_r
  648.        dfsw=d2mw_r*rold_l*ut1_l+mmin_m*rold_r*c23/det
  649.        dfsw=dfsw+d3mw_r*rold_r*ut1_r
  650.        dfsp=d2mp_r*rold_l*ut1_l+d3mp_r*rold_r*ut1_r
  651. c------------------------------------------------------
  652.        dftr=d2mr_r*rold_l*ut2_l+mmin_m*ut2_r
  653.        dftr=dftr+d3mr_r*rold_r*ut2_r
  654.        dftu=d2mu_r*rold_l*ut2_l+mmin_m*rold_r*c31/det
  655.        dftu=dftu+d3mu_r*rold_r*ut2_r
  656.        dftv=d2mv_r*rold_l*ut2_l+mmin_m*rold_r*c32/det
  657.        dftv=dftv+d3mv_r*rold_r*ut2_r
  658.        dftw=d2mw_r*rold_l*ut2_l+mmin_m*rold_r*c33/det
  659.        dftw=dftw+d3mw_r*rold_r*ut2_r
  660.        dftp=d2mp_r*rold_l*ut2_l+d3mp_r*rold_r*ut2_r
  661. c-------------------------------------------------------
  662.        jr(2,1)=damr_r*(f(2)-pmid*n_x)/am+dpir_r*n_x
  663.        jr(2,1)=jr(2,1)+am*(dffr*n_x+dfsr*t1_x+dftr*t2_x)
  664.        jr(2,2)=damu_r*(f(2)-pmid*n_x)/am+dpiu_r*n_x
  665.        jr(2,2)=jr(2,2)+am*(dffu*n_x+dfsu*t1_x+dftu*t2_x)
  666.        jr(2,3)=damv_r*(f(2)-pmid*n_x)/am+dpiv_r*n_x
  667.        jr(2,3)=jr(2,3)+am*(dffv*n_x+dfsv*t1_x+dftv*t2_x)
  668.        jr(2,4)=damw_r*(f(2)-pmid*n_x)/am+dpiw_r*n_x
  669.        jr(2,4)=jr(2,4)+am*(dffw*n_x+dfsw*t1_x+dftw*t2_x)
  670.        jr(2,5)=damp_r*(f(2)-pmid*n_x)/am+dpip_r*n_x
  671.        jr(2,5)=jr(2,5)+am*(dffp*n_x+dfsp*t1_x+dftp*t2_x)
  672. c-------------------------------------------------------
  673.        jr(3,1)=damr_r*(f(3)-pmid*n_y)/am+dpir_r*n_y
  674.        jr(3,1)=jr(3,1)+am*(dffr*n_y+dfsr*t1_y+dftr*t2_y)
  675.        jr(3,2)=damu_r*(f(3)-pmid*n_y)/am+dpiu_r*n_y
  676.        jr(3,2)=jr(3,2)+am*(dffu*n_y+dfsu*t1_y+dftu*t2_y)
  677.        jr(3,3)=damv_r*(f(3)-pmid*n_y)/am+dpiv_r*n_y
  678.        jr(3,3)=jr(3,3)+am*(dffv*n_y+dfsv*t1_y+dftv*t2_y)
  679.        jr(3,4)=damw_r*(f(3)-pmid*n_y)/am+dpiw_r*n_y
  680.        jr(3,4)=jr(3,4)+am*(dffw*n_y+dfsw*t1_y+dftw*t2_y)
  681.        jr(3,5)=damp_r*(f(3)-pmid*n_y)/am+dpip_r*n_y
  682.        jr(3,5)=jr(3,5)+am*(dffp*n_y+dfsp*t1_y+dftp*t2_y)
  683. c--------------------------------------------------------
  684.        jr(4,1)=damr_r*(f(4)-pmid*n_z)/am+dpir_r*n_z
  685.        jr(4,1)=jr(4,1)+am*(dffr*n_z+dfsr*t1_z+dftr*t2_z)
  686.        jr(4,2)=damu_r*(f(4)-pmid*n_z)/am+dpiu_r*n_z
  687.        jr(4,2)=jr(4,2)+am*(dffu*n_z+dfsu*t1_z+dftu*t2_z)
  688.        jr(4,3)=damv_r*(f(4)-pmid*n_z)/am+dpiv_r*n_z
  689.        jr(4,3)=jr(4,3)+am*(dffv*n_z+dfsv*t1_z+dftv*t2_z)
  690.        jr(4,4)=damw_r*(f(4)-pmid*n_z)/am+dpiw_r*n_z
  691.        jr(4,4)=jr(4,4)+am*(dffw*n_z+dfsw*t1_z+dftw*t2_z)
  692.        jr(4,5)=damp_r*(f(4)-pmid*n_z)/am+dpip_r*n_z
  693.        jr(4,5)=jr(4,5)+am*(dffp*n_z+dfsp*t1_z+dftp*t2_z)
  694. c-------------------------------------------------------
  695. c-------------------------------------------------------
  696. c-------------------------------------------------------
  697.        hr_l=(eold_l+pold_l/rold_l)*rold_l
  698.        hr_r=(eold_r+pold_r/rold_r)*rold_r
  699.        f(5)=am*(mpl_m*hr_l+mmin_m*hr_r)
  700. c-------------------------------------------------
  701.        br1=d2mr_l*hr_l+d3mr_l*hr_r
  702.        br1=br1+mpl_m*(uold_l*uold_l+vold_l*vold_l+wold_l*wold_l)/2.0d0
  703.        jl(5,1)=damr_l*f(5)/am+am*br1
  704. c-------------------------------------------------
  705.        br1=d2mu_l*hr_l+mpl_m*uold_l*rold_l
  706.        br1=br1+d3mu_l*hr_r
  707.        jl(5,2)=damu_l*f(5)/am+am*br1
  708. c-------------------------------------------------
  709.        br1=d2mv_l*hr_l+mpl_m*vold_l*rold_l
  710.        br1=br1+d3mv_l*hr_r
  711.        jl(5,3)=damv_l*f(5)/am+am*br1
  712. c-------------------------------------------------
  713.        br1=d2mw_l*hr_l+mpl_m*wold_l*rold_l
  714.        br1=br1+d3mw_l*hr_r
  715.        jl(5,4)=damw_l*f(5)/am+am*br1
  716. c-------------------------------------------------
  717.        br1=d2mp_l*hr_l+mpl_m*ga/gm1
  718.        br1=br1+d3mp_l*hr_r
  719.        jl(5,5)=damp_l*f(5)/am+am*br1
  720. c-------------------------------------------------
  721. c-------------------------------------------------
  722.        br1=d2mr_r*hr_l+d3mr_r*hr_r
  723.        br1=br1+mmin_m*(uold_r*uold_r+vold_r*vold_r+wold_r*wold_r)/2.0d0
  724.        jr(5,1)=damr_r*f(5)/am+am*br1
  725. c---------------------
  726.        br1=d2mu_r*hr_l+mmin_m*uold_r*rold_r
  727.        br1=br1+d3mu_r*hr_r
  728.        jr(5,2)=damu_r*f(5)/am+am*br1
  729. c---------------------
  730.        br1=d2mv_r*hr_l+mmin_m*vold_r*rold_r
  731.        br1=br1+d3mv_r*hr_r
  732.        jr(5,3)=damv_r*f(5)/am+am*br1
  733. c---------------------
  734.        br1=d2mw_r*hr_l+mmin_m*wold_r*rold_r
  735.        br1=br1+d3mw_r*hr_r
  736.        jr(5,4)=damw_r*f(5)/am+am*br1
  737. c---------------------
  738.        br1=d2mp_r*hr_l+mmin_m*ga/gm1
  739.        br1=br1+d3mp_r*hr_r
  740.        jr(5,5)=damp_r*f(5)/am+am*br1
  741. c-------------------------------------------------------------
  742. c-------------------------------------------------------------
  743.        return
  744.        end
  745.  
  746.  
  747.  
  748.  
  749.  
  750.  
  751.  
  752.  
  753.  
  754.  
  755.  
  756.  
  757.  
  758.  
  759.  
  760.  
  761.  
  762.  
  763.  
  764.  
  765.  
  766.  
  767.  
  768.  
  769.  

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