Télécharger la plaquette Cast3M


fr en

Télécharger conjak.eso

Retour à la liste

Numérotation des lignes : 

conjak

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

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