Télécharger conj3d.eso

Retour à la liste

Numérotation des lignes :

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

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