Télécharger cwms3d.eso

Retour à la liste

Numérotation des lignes :

  1. C CWMS3D SOURCE CHAT 05/01/12 22:33:32 5004
  2. SUBROUTINE CWMS3D(NSP,JLL,WVEC_L,WVEC_R,NVECT,TVEC1,TVEC2,
  3. & mpyn,lrecp,lrecv,nlcg)
  4. C************************************************************************
  5. C
  6. C PROJET : CASTEM 2000
  7. C
  8. C NOM : CWMS3D ('convection at wall for multispicies en 3D')
  9. C
  10. C DESCRIPTION : Voir KONMS3 (appele par KONMS3.ESO)
  11. C
  12. C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
  13. C
  14. C AUTEUR : S. KUDRIAKOV, DM2S/SFME/LTMF
  15. C
  16. C************************************************************************
  17. C
  18. c----------------------------------------------------------------------
  19. c GENERAL DESCRIPTION:
  20. c This subroutine provides the jacobian which is the derivatives
  21. c of the numerical flux function defined at the wall cell interface
  22. c with respect to the conservative variables of the left (pre-wall)
  23. c cell.
  24. c
  25. c EQUATIONS: 3D Euler equations of gas dynamics - MULTICPECIES GAS
  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 nsp -- number of species (total);
  40. c
  41. c wvec_l -- vector of the primitive variables
  42. c (rho,u_x,u_y,u_z,p) at the left cell;
  43. c
  44. c wvec_r -- vector of the primitive variables
  45. c (rho,u_x,u_y,u_z,p) at the right cell;
  46. c
  47. c nvect -- normal vector to the interface (2 components in 2D);
  48. c
  49. c tvec1 -- first tangential vector to the interface;
  50. c
  51. c tvec2 -- second tangential vector to the interface;
  52. c
  53. c mpyn -- pointer to the vectors of the primitive variables
  54. c (Y_1,Y_2,...Y_(nsp-1)) at the left and the right cells;
  55. c
  56. c lrecp -- pointer to the vector of specific heats at constant pressure
  57. c (size of the vector is equal to number of species (nsp));
  58. c
  59. c lrecv -- pointer to the vector of specific heats at constant volume
  60. c (size of the vector is equal to number of species (nsp));
  61. c
  62. c nlcg -- "local" number corresponding to the cell at wall.
  63. c----------------------------------------------------------------------
  64. c
  65. c OUTPUT:
  66. c
  67. c jll -- jakobian matrix (4+nsp) by (4+nsp) -
  68. c derivatives of the numerical
  69. c flux function with respect to the conservative variables
  70. c from the left cell;
  71. c
  72. c----------------------------------------------------------------------
  73. IMPLICIT INTEGER(I-N)
  74. integer nsp,jll,lrecp,lrecv,nlcg
  75. real*8 wvec_l(5),wvec_r(5)
  76. real*8 nvect(3),tvec1(3),tvec2(3)
  77. real*8 alpha,beta
  78. real*8 gal,gar,gm1l,gm1r
  79. real*8 n_x,n_y,n_z,t1_x,t1_y,t1_z,t2_x,t2_y,t2_z
  80. real*8 un_l, un_r
  81. real*8 ml,mr,Mplus,Mmin,am
  82. real*8 rold_l,uold_l,vold_l,wold_l,pold_l,eold_l
  83. real*8 rold_r,uold_r,vold_r,wold_r,pold_r,eold_r
  84. real*8 Pplus,Pmin
  85. real*8 top,bot
  86. real*8 det,c11,c12,c13,c21,c22,c23,c31,c32,c33
  87. real*8 br1,temp_l,temp_r,brac_l,brac_r
  88. real*8 aleft, arigh
  89. real*8 damr_l,damr_r,damu_l,damu_r
  90. real*8 damv_l,damv_r,damp_l,damp_r
  91. real*8 damg_l,damg_r,damw_l,damw_r
  92. real*8 dmlr_l,dmlr_r,dmlu_l,dmlu_r
  93. real*8 dmlv_l,dmlv_r,dmlp_l,dmlp_r
  94. real*8 dmlw_l,dmlw_r
  95. real*8 dmrr_l,dmrr_r,dmru_l,dmru_r
  96. real*8 dmrv_l,dmrv_r,dmrp_l,dmrp_r
  97. real*8 dmrw_l,dmrw_r
  98. real*8 dPpr_l,dPpr_r,dPpu_l,dPpu_r
  99. real*8 dPpv_l,dPpv_r,dPpp_l,dPpp_r
  100. real*8 dPpw_l,dPpw_r
  101. real*8 dPmr_l,dPmr_r,dPmu_l,dPmu_r
  102. real*8 dPmv_l,dPmv_r,dPmp_l,dPmp_r
  103. real*8 dPmw_l,dPmw_r
  104. real*8 dpir_l,dpir_r,dpiu_l,dpiu_r
  105. real*8 dpiv_l,dpiv_r,dpip_l,dpip_r
  106. real*8 dpiw_l,dpiw_r
  107. real*8 zc11,zc12,zc13,zc21,zc22,zc23,zc31,zc32,zc33
  108. integer i,j,k
  109. parameter(alpha = 0.1875D0, beta = 0.125D0)
  110. C------------------------------------------------------------
  111. -INC SMCHPOI
  112. POINTEUR MPYN.MPOVAL
  113. C-------------------------------------------------------------
  114. -INC SMLREEL
  115. POINTEUR MLRECP.MLREEL, MLRECV.MLREEL
  116. C-------------------------------------------------------------
  117. C******* Les fractionines messiques **************************
  118. C-------------------------------------------------------------
  119. SEGMENT FRAMAS
  120. REAL*8 YET(NSP)
  121. ENDSEGMENT
  122. POINTEUR YL.FRAMAS, YR.FRAMAS
  123. C-------------------------------------------------------
  124. C********** Les CP's and CV's ***********************
  125. C-------------------------------------------------------
  126. SEGMENT GCONST
  127. REAL*8 GC(NSP)
  128. ENDSEGMENT
  129. POINTEUR CP.GCONST, CV.GCONST
  130. C-------------------------------------------------------------
  131. C******** Segments for the elementary matrixes *************
  132. C-------------------------------------------------------------
  133. SEGMENT JACEL
  134. REAL*8 JAC(4+NSP,4+NSP)
  135. ENDSEGMENT
  136. POINTEUR JTL.JACEL, JTR.JACEL, JL.JACEL, JR.JACEL,
  137. & WL.JACEL, WR.JACEL
  138. C-------------------------------------------------------------
  139. C********** Segments for the vectors ***********************
  140. C-------------------------------------------------------------
  141. SEGMENT VECEL
  142. REAL*8 VV(NSP)
  143. ENDSEGMENT
  144. POINTEUR dmly_l.vecel, dmly_r.vecel,
  145. & dmry_l.vecel, dmry_r.vecel,
  146. & dPpy_l.vecel, dPpy_r.vecel,
  147. & dPmy_l.vecel, dPmy_r.vecel,
  148. & dpiy_l.vecel, dpiy_r.vecel,
  149. & dgdyl.vecel, dgdyr.vecel
  150. C-------------------------------------------
  151. SEGINI dmly_l, dmly_r,
  152. & dmry_l, dmry_r,
  153. & dPpy_l, dPpy_r,
  154. & dPmy_l, dPmy_r,
  155. & dpiy_l, dpiy_r,
  156. & dgdyl, dgdyr
  157. C------------------------------------------------------------
  158. SEGINI YL, YR
  159. SEGACT MPYN
  160. DO 4 I=1,(NSP-1)
  161. YL.YET(I)=MPYN.VPOCHA(NLCG,I)
  162. YR.YET(I)=MPYN.VPOCHA(NLCG,I)
  163. 4 CONTINUE
  164. C----------------------------------------
  165. SEGINI CP, CV
  166. MLRECP = LRECP
  167. MLRECV = LRECV
  168. SEGACT MLRECP, MLRECV
  169. DO 5 I=1,(NSP-1)
  170. CP.GC(I)=MLRECP.PROG(I)
  171. CV.GC(I)=MLRECV.PROG(I)
  172. 5 CONTINUE
  173. CP.GC(NSP)=MLRECP.PROG(NSP)
  174. CV.GC(NSP)=MLRECV.PROG(NSP)
  175. c-------------------------------------------------------------
  176. c Computing GAMMA at the left cell and its derivatives
  177. c with respect to the primitive variables Y_i
  178. c-------------------------------------------------------------
  179. top=0.0D0
  180. bot=0.0D0
  181. do 40 i=1,(nsp-1)
  182. top=top+yl.yet(i)*(cp.gc(i)-cp.gc(nsp))
  183. bot=bot+yl.yet(i)*(cv.gc(i)-cv.gc(nsp))
  184. 40 continue
  185. top=cp.gc(nsp)+top
  186. bot=cv.gc(nsp)+bot
  187. gal=top/bot
  188. gm1l=gal-1.0d0
  189. c-------------------------------------------------------------
  190. do 41 i=1,(nsp-1)
  191. dgdyl.vv(i)=(cp.gc(i)-cp.gc(nsp)-
  192. & gal*(cv.gc(i)-cv.gc(nsp)))/bot
  193. 41 continue
  194. c-------------------------------------------------------------
  195. c Computing GAMMA at the right cell and its derivatives
  196. c with respect to the primitive variables Y_i
  197. c-------------------------------------------------------------
  198. top=0.0D0
  199. bot=0.0D0
  200. do 42 i=1,(nsp-1)
  201. top=top+yr.yet(i)*(cp.gc(i)-cp.gc(nsp))
  202. bot=bot+yr.yet(i)*(cv.gc(i)-cv.gc(nsp))
  203. 42 continue
  204. top=cp.gc(nsp)+top
  205. bot=cv.gc(nsp)+bot
  206. gar=top/bot
  207. gm1r=gar-1.0d0
  208. c-------------------------------------------------------------
  209. do 43 i=1,(nsp-1)
  210. dgdyr.vv(i)=(cp.gc(i)-cp.gc(nsp)-
  211. & gar*(cv.gc(i)-cv.gc(nsp)))/bot
  212. 43 continue
  213. c-------------------------------------------------------------
  214. n_x=nvect(1)
  215. n_y=nvect(2)
  216. n_z=nvect(3)
  217. c-------------------
  218. t1_x=tvec1(1)
  219. t1_y=tvec1(2)
  220. t1_z=tvec1(3)
  221. c-------------------
  222. t2_x=tvec2(1)
  223. t2_y=tvec2(2)
  224. t2_z=tvec2(3)
  225. c-------------------------------------
  226. c11=t1_y*t2_z - t1_z*t2_y
  227. c12=t1_z*t2_x - t1_x*t2_z
  228. c13=t1_x*t2_y - t1_y*t2_x
  229. det=n_x*c11 + n_y*c12 + n_z*c13
  230. c--------------------------------------
  231. c21=n_z*t2_y - n_y*t2_z
  232. c22=n_x*t2_z - n_z*t2_x
  233. c23=n_y*t2_x - n_x*t2_y
  234. c--------------------------------------
  235. c31=n_y*t1_z - n_z*t1_y
  236. c32=n_z*t1_x - n_x*t1_z
  237. c33=n_x*t1_y - n_y*t1_x
  238. c----------------------------
  239. c----------------------------
  240. rold_l=wvec_l(1)
  241. uold_l=wvec_l(2)
  242. vold_l=wvec_l(3)
  243. wold_l=wvec_l(4)
  244. pold_l=wvec_l(5)
  245. c-----------------------
  246. rold_r=wvec_r(1)
  247. uold_r=wvec_r(2)
  248. vold_r=wvec_r(3)
  249. wold_r=wvec_r(4)
  250. pold_r=wvec_r(5)
  251. c------------------------------------------------------------------
  252. c Computation of the specific total energy on the left and right.
  253. c------------------------------------------------------------------
  254. eold_l=(uold_l*uold_l+vold_l*vold_l+wold_l*wold_l)/2.0d0
  255. eold_l=eold_l+pold_l/(gm1l*rold_l)
  256. eold_r=(uold_r*uold_r+vold_r*vold_r+wold_r*wold_r)/2.0d0
  257. eold_r=eold_r+pold_r/(gm1r*rold_r)
  258. c-------------------------------------------------------------------
  259. c Computation of the speed of sound and its derivatives;
  260. c numerical speed of sound at the interface is taken as an average
  261. c of the speeds of sounds of the neighbouring cells
  262. c-------------------------------------------------------------------
  263. aleft=sqrt(gal*pold_l/rold_l)
  264. arigh=sqrt(gar*pold_r/rold_r)
  265. am=0.5d0*(aleft+arigh)
  266. c-------------------------------------------------------------------
  267. damr_r=-arigh/(4.0d0*rold_r)
  268. damu_r=0.0d0
  269. damv_r=0.0d0
  270. damw_r=0.0d0
  271. damp_r=gar/(4.0d0*arigh*rold_r)
  272. damg_r=arigh/(4.0d0*gar)
  273. c-----------------------
  274. damr_l=-aleft/(4.0d0*rold_l)
  275. damu_l=0.0d0
  276. damv_l=0.0d0
  277. damw_l=0.0d0
  278. damp_l=gal/(4.0d0*aleft*rold_l)
  279. damg_l=aleft/(4.0d0*gal)
  280. c-------------------------------------------------------------------
  281. c Computing numerical Mach number and its derivatives,
  282. c see p.370, under (A1).
  283. c-------------------------------------------------------------------
  284. un_l=uold_l*n_x+vold_l*n_y+wold_l*n_z
  285. un_r=uold_r*n_x+vold_r*n_y+wold_r*n_z
  286. c-------------------------------------------------------------------
  287. ml=un_l/am
  288. mr=un_r/am
  289. c-------------------------------------------------------------------
  290. c Mplus and Mmin are calligraphic lettes M+ and M- from the paper,
  291. c see (19a) and (19b), p.367.
  292. c-------------------------------------------------------------------
  293. if(abs(ml) .ge. 1.0d0) then
  294. Mplus=(ml+abs(ml))/2.0d0
  295. else
  296. Mplus=(ml+1.0d0)*(ml+1.0d0)/4.0d0
  297. Mplus=Mplus+beta*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  298. endif
  299. c-------------------------------------------------------------------
  300. if(abs(mr) .ge. 1.0d0) then
  301. Mmin=(mr-abs(mr))/2.0d0
  302. else
  303. Mmin=-(mr-1.0d0)*(mr-1.0d0)/4.0d0
  304. Mmin=Mmin-beta*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  305. endif
  306. c-------------------------------------------------------------------
  307. c Derivatives of ml and mr with respect to both sets of primitive
  308. c variables: left and right.
  309. c-------------------------------------------------------------------
  310. temp_l=-un_l/(am*am)
  311. temp_r=-un_r/(am*am)
  312. c--------
  313. dmlr_l=temp_l*damr_l
  314. dmlr_r=temp_l*damr_r
  315. dmrr_l=temp_r*damr_l
  316. dmrr_r=temp_r*damr_r
  317. c--------
  318. dmlu_l=n_x/am+temp_l*damu_l
  319. dmlu_r=temp_l*damu_r
  320. dmru_l=temp_r*damu_l
  321. dmru_r=n_x/am+temp_r*damu_r
  322. c--------
  323. dmlv_l=n_y/am+temp_l*damv_l
  324. dmlv_r=temp_l*damv_r
  325. dmrv_l=temp_r*damv_l
  326. dmrv_r=n_y/am+temp_r*damv_r
  327. c--------
  328. dmlw_l=n_z/am+temp_l*damw_l
  329. dmlw_r=temp_l*damw_r
  330. dmrw_l=temp_r*damw_l
  331. dmrw_r=n_z/am+temp_r*damw_r
  332. c--------
  333. dmlp_l=temp_l*damp_l
  334. dmlp_r=temp_l*damp_r
  335. dmrp_l=temp_r*damp_l
  336. dmrp_r=temp_r*damp_r
  337. c---------------------------------
  338. do 44 i=1,(nsp-1)
  339. dmly_l.vv(i)=temp_l*damg_l*dgdyl.vv(i)
  340. dmly_r.vv(i)=temp_l*damg_r*dgdyr.vv(i)
  341. dmry_l.vv(i)=temp_r*damg_l*dgdyl.vv(i)
  342. dmry_r.vv(i)=temp_r*damg_r*dgdyr.vv(i)
  343. 44 continue
  344. c---------------------------------------------------------------
  345. c Computing the calligraphic P+ and P- with their derivatives,
  346. c see (21a) & (21b) on p.368.
  347. c---------------------------------------------------------------
  348. if(ml .ge. 1.0d0) then
  349. Pplus = 1.0d0
  350. else
  351. if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  352. Pplus=(ml+1.0d0)*(ml+1.0d0)*(2.0d0-ml)/4.0d0
  353. Pplus=Pplus+alpha*ml*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  354. else
  355. Pplus = 0.0d0
  356. endif
  357. endif
  358. c---------------------------------------------------------------
  359. if(mr .ge. 1.0d0) then
  360. Pmin = 0.0d0
  361. else
  362. if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  363. Pmin=(mr-1.0d0)*(mr-1.0d0)*(2.0d0+mr)/4.0d0
  364. Pmin=Pmin-alpha*mr*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  365. else
  366. Pmin = 1.0d0
  367. endif
  368. endif
  369. c---------------------------------------------------------------
  370. brac_l=(ml+1.0d0)*(2.0d0-ml)/2.0d0-(ml+1.0d0)*(ml+1.0d0)/4.0d0
  371. brac_l=brac_l+alpha*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  372. brac_l=brac_l+4.0d0*alpha*ml*ml*(ml*ml-1.0d0)
  373. c--------------
  374. brac_r=(mr-1.0d0)*(2.0d0+mr)/2.0d0+(mr-1.0d0)*(mr-1.0d0)/4.0d0
  375. brac_r=brac_r-alpha*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  376. brac_r=brac_r-4.0d0*alpha*mr*mr*(mr*mr-1.0d0)
  377. c---------------------------------------------------------------
  378. if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  379. dPpr_l=brac_l*dmlr_l
  380. dPpr_r=brac_l*dmlr_r
  381. c------------
  382. dPpu_l=brac_l*dmlu_l
  383. dPpu_r=brac_l*dmlu_r
  384. c------------
  385. dPpv_l=brac_l*dmlv_l
  386. dPpv_r=brac_l*dmlv_r
  387. c------------
  388. dPpw_l=brac_l*dmlw_l
  389. dPpw_r=brac_l*dmlw_r
  390. c------------
  391. dPpp_l=brac_l*dmlp_l
  392. dPpp_r=brac_l*dmlp_r
  393. c----------------------------------------------
  394. do 66 i=1,(nsp-1)
  395. dPpy_l.vv(i)=brac_l*dmly_l.vv(i)
  396. dPpy_r.vv(i)=brac_l*dmly_r.vv(i)
  397. 66 continue
  398. c------------
  399. else
  400. dPpr_l=0.0d0
  401. dPpr_r=0.0d0
  402. c-----------
  403. dPpu_l=0.0d0
  404. dPpu_r=0.0d0
  405. c-----------
  406. dPpv_l=0.0d0
  407. dPpv_r=0.0d0
  408. c-----------
  409. dPpw_l=0.0d0
  410. dPpw_r=0.0d0
  411. c-----------
  412. dPpp_l=0.0d0
  413. dPpp_r=0.0d0
  414. c-----------
  415. do 67 i=1,(nsp-1)
  416. dPpy_l.vv(i)=0.0d0
  417. dPpy_r.vv(i)=0.0d0
  418. 67 continue
  419. c-----------
  420. endif
  421. c---------------------------------------------------------------
  422. if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  423. dPmr_l=brac_r*dmrr_l
  424. dPmr_r=brac_r*dmrr_r
  425. c------------
  426. dPmu_l=brac_r*dmru_l
  427. dPmu_r=brac_r*dmru_r
  428. c------------
  429. dPmv_l=brac_r*dmrv_l
  430. dPmv_r=brac_r*dmrv_r
  431. c------------
  432. dPmw_l=brac_r*dmrw_l
  433. dPmw_r=brac_r*dmrw_r
  434. c------------
  435. dPmp_l=brac_r*dmrp_l
  436. dPmp_r=brac_r*dmrp_r
  437. c------------
  438. do 68 i=1,(nsp-1)
  439. dPmy_l.vv(i)=brac_r*dmry_l.vv(i)
  440. dPmy_r.vv(i)=brac_r*dmry_r.vv(i)
  441. 68 continue
  442. c------------
  443. else
  444. dPmr_l=0.0d0
  445. dPmr_r=0.0d0
  446. c-----------
  447. dPmu_l=0.0d0
  448. dPmu_r=0.0d0
  449. c-----------
  450. dPmv_l=0.0d0
  451. dPmv_r=0.0d0
  452. c-----------
  453. dPmw_l=0.0d0
  454. dPmw_r=0.0d0
  455. c-----------
  456. dPmp_l=0.0d0
  457. dPmp_r=0.0d0
  458. c-----------
  459. do 69 i=1,(nsp-1)
  460. dPmy_l.vv(i)=0.0d0
  461. dPmy_r.vv(i)=0.0d0
  462. 69 continue
  463. c-----------
  464. endif
  465. c-------------------------------------------------------------------
  466. c computing pmid -- p_{1/2} and its derivatives, see (20b), p.367.
  467. c-------------------------------------------------------------------
  468. dpir_l=dPpr_l*pold_l+dPmr_l*pold_r
  469. dpiu_l=dPpu_l*pold_l+dPmu_l*pold_r
  470. dpiv_l=dPpv_l*pold_l+dPmv_l*pold_r
  471. dpiw_l=dPpw_l*pold_l+dPmw_l*pold_r
  472. dpip_l=dPpp_l*pold_l+Pplus+dPmp_l*pold_r
  473. do 70 i=1,(nsp-1)
  474. dpiy_l.vv(i)=dPpy_l.vv(i)*pold_l+dPmy_l.vv(i)*pold_r
  475. 70 continue
  476. c----------------------------
  477. dpir_r=dPpr_r*pold_l+dPmr_r*pold_r
  478. dpiu_r=dPpu_r*pold_l+dPmu_r*pold_r
  479. dpiv_r=dPpv_r*pold_l+dPmv_r*pold_r
  480. dpiw_r=dPpw_r*pold_l+dPmw_r*pold_r
  481. dpip_r=dPpp_r*pold_l+Pmin+dPmp_r*pold_r
  482. do 71 i=1,(nsp-1)
  483. dpiy_r.vv(i)=dPpy_r.vv(i)*pold_l+dPmy_r.vv(i)*pold_r
  484. 71 continue
  485. c---------------------------------------------------------------------
  486. c computing JACOBIAN as a derivative of the numerical flux function
  487. c with respect to the primitive variables.
  488. c Notation: jl(i,j) --- is the derivative of the i-component of the
  489. c flux function with respect to the j-component of the
  490. c vector of primitive variables of the left state.
  491. c jr(i,j) --- is the derivative of the i-component of the
  492. c flux function with respect to the j-component of the
  493. c vector of primitive variables of the right state.
  494. c---------------------------------------------------------------------
  495. SEGINI JL, JR
  496. c---------------------------------------------------------------------
  497. jl.jac(1,1)=0.0d0
  498. jl.jac(1,2)=0.0d0
  499. jl.jac(1,3)=0.0d0
  500. jl.jac(1,4)=0.0d0
  501. jl.jac(1,5)=0.0d0
  502. do 72 i=1,(nsp-1)
  503. jl.jac(1,5+i)=0.0d0
  504. 72 continue
  505. c------------------------------------
  506. jr.jac(1,1)=0.0d0
  507. jr.jac(1,2)=0.0d0
  508. jr.jac(1,3)=0.0d0
  509. jr.jac(1,4)=0.0d0
  510. jr.jac(1,5)=0.0d0
  511. do 73 i=1,(nsp-1)
  512. jr.jac(1,5+i)=0.0d0
  513. 73 continue
  514. c------------------------------------
  515. c------------------------------------
  516. c---------------------------------------------------------
  517. jl.jac(2,1)=dpir_l*n_x
  518. jl.jac(2,2)=dpiu_l*n_x
  519. jl.jac(2,3)=dpiv_l*n_x
  520. jl.jac(2,4)=dpiw_l*n_x
  521. jl.jac(2,5)=dpip_l*n_x
  522. do 77 i=1,(nsp-1)
  523. jl.jac(2,5+i)=dpiy_l.vv(i)*n_x
  524. 77 continue
  525. c-------------------------------------------------
  526. jl.jac(3,1)=dpir_l*n_y
  527. jl.jac(3,2)=dpiu_l*n_y
  528. jl.jac(3,3)=dpiv_l*n_y
  529. jl.jac(3,4)=dpiw_l*n_y
  530. jl.jac(3,5)=dpip_l*n_y
  531. do 78 i=1,(nsp-1)
  532. jl.jac(3,5+i)=dpiy_l.vv(i)*n_y
  533. 78 continue
  534. c-------------------------------------------------
  535. jl.jac(4,1)=dpir_l*n_z
  536. jl.jac(4,2)=dpiu_l*n_z
  537. jl.jac(4,3)=dpiv_l*n_z
  538. jl.jac(4,4)=dpiw_l*n_z
  539. jl.jac(4,5)=dpip_l*n_z
  540. do 79 i=1,(nsp-1)
  541. jl.jac(4,5+i)=dpiy_l.vv(i)*n_z
  542. 79 continue
  543. c-------------------------------------------------------
  544. c derivatives with respect to the right
  545. c set of the primitive variables
  546. c-------------------------------------------------------
  547. jr.jac(2,1)=dpir_r*n_x
  548. jr.jac(2,2)=dpiu_r*n_x
  549. jr.jac(2,3)=dpiv_r*n_x
  550. jr.jac(2,4)=dpiw_r*n_x
  551. jr.jac(2,5)=dpip_r*n_x
  552. do 84 i=1,(nsp-1)
  553. jr.jac(2,5+i)=dpiy_r.vv(i)*n_x
  554. 84 continue
  555. c-------------------------------------------------------
  556. jr.jac(3,1)=dpir_r*n_y
  557. jr.jac(3,2)=dpiu_r*n_y
  558. jr.jac(3,3)=dpiv_r*n_y
  559. jr.jac(3,4)=dpiw_r*n_y
  560. jr.jac(3,5)=dpip_r*n_y
  561. do 85 i=1,(nsp-1)
  562. jr.jac(3,5+i)=dpiy_r.vv(i)*n_y
  563. 85 continue
  564. c--------------------------------------------------------
  565. jr.jac(4,1)=dpir_r*n_z
  566. jr.jac(4,2)=dpiu_r*n_z
  567. jr.jac(4,3)=dpiv_r*n_z
  568. jr.jac(4,4)=dpiw_r*n_z
  569. jr.jac(4,5)=dpip_r*n_z
  570. do 86 i=1,(nsp-1)
  571. jr.jac(4,5+i)=dpiy_r.vv(i)*n_z
  572. 86 continue
  573. c-------------------------------------------------------
  574. jl.jac(5,1)=0.0d0
  575. jl.jac(5,2)=0.0d0
  576. jl.jac(5,3)=0.0d0
  577. jl.jac(5,4)=0.0d0
  578. jl.jac(5,5)=0.0d0
  579. do 87 i=1,(nsp-1)
  580. jl.jac(5,5+i)=0.0d0
  581. 87 continue
  582. c-------------------------------------------------------------
  583. do 180 i=1,(nsp-1)
  584. c---------------------
  585. jl.jac(5+i,1)=0.0d0
  586. jl.jac(5+i,2)=0.0d0
  587. jl.jac(5+i,3)=0.0d0
  588. jl.jac(5+i,4)=0.0d0
  589. jl.jac(5+i,5)=0.0d0
  590. do 181 j=6,(4+nsp)
  591. jl.jac(5+i,j)=0.0d0
  592. 181 continue
  593. 180 continue
  594. c-------------------------------------------------
  595. jr.jac(5,1)=0.0d0
  596. jr.jac(5,2)=0.0d0
  597. jr.jac(5,3)=0.0d0
  598. jr.jac(5,4)=0.0d0
  599. jr.jac(5,5)=0.0d0
  600. do 88 i=1,(nsp-1)
  601. jr.jac(5,5+i)=0.0d0
  602. 88 continue
  603. c----------------------------------------------------------------
  604. do 182 i=1,(nsp-1)
  605. jr.jac(5+i,1)=0.0d0
  606. jr.jac(5+i,2)=0.0d0
  607. jr.jac(5+i,3)=0.0d0
  608. jr.jac(5+i,4)=0.0d0
  609. jr.jac(5+i,5)=0.0d0
  610. do 183 j=1,(nsp-1)
  611. jr.jac(5+i,5+j)=0.0d0
  612. 183 continue
  613. 182 continue
  614. c-------------------------------------------------------------
  615. c matrix wl(i,j) represents the derivative of the i-component
  616. c of the vector of primitive variables of the left state with
  617. c respect to the j-component of the vector of the conservative
  618. c variables of the left state.
  619. c
  620. c Here: (rho, ux, uy, p, Y_1,...,Y_(nsp-1)) -
  621. c vector of primitive variables;
  622. c (rho, rho ux, rho uy, rho e, rho Y_1,..., rho Y_(nsp-1)) -
  623. c vector of conservative variables.
  624. c-------------------------------------------------------------
  625. SEGINI WL, WR
  626. c-------------------------------------------------------------
  627. wl.jac(1,1)=1.0d0
  628. wl.jac(1,2)=0.0d0
  629. wl.jac(1,3)=0.0d0
  630. wl.jac(1,4)=0.0d0
  631. wl.jac(1,5)=0.0d0
  632. do 83 i=1,(nsp-1)
  633. wl.jac(1,5+i)=0.0d0
  634. 83 continue
  635. c------------------------------
  636. wl.jac(2,1)=-uold_l/rold_l
  637. wl.jac(2,2)=1.0d0/rold_l
  638. wl.jac(2,3)=0.0d0
  639. wl.jac(2,4)=0.0d0
  640. wl.jac(2,5)=0.0d0
  641. do 840 i=1,(nsp-1)
  642. wl.jac(2,5+i)=0.0d0
  643. 840 continue
  644. c------------------------------
  645. wl.jac(3,1)=-vold_l/rold_l
  646. wl.jac(3,2)=0.0d0
  647. wl.jac(3,3)=1.0d0/rold_l
  648. wl.jac(3,4)=0.0d0
  649. wl.jac(3,5)=0.0d0
  650. do 850 i=1,(nsp-1)
  651. wl.jac(3,5+i)=0.0d0
  652. 850 continue
  653. c------------------------------
  654. wl.jac(4,1)=-wold_l/rold_l
  655. wl.jac(4,2)=0.0d0
  656. wl.jac(4,3)=0.0d0
  657. wl.jac(4,4)=1.0d0/rold_l
  658. wl.jac(4,5)=0.0d0
  659. do 92 i=1,(nsp-1)
  660. wl.jac(4,5+i)=0.0d0
  661. 92 continue
  662. c------------------------------
  663. br1=0.0d0
  664. do 93 i=1,(nsp-1)
  665. br1=br1+dgdyl.vv(i)*yl.yet(i)
  666. 93 continue
  667. br1=br1*pold_l/(rold_l*gm1l)
  668. wl.jac(5,1)=gm1l*(uold_l*uold_l+vold_l*vold_l+
  669. & wold_l*wold_l)/2.0d0
  670. wl.jac(5,1)=wl.jac(5,1)-br1
  671. wl.jac(5,2)=-uold_l*gm1l
  672. wl.jac(5,3)=-vold_l*gm1l
  673. wl.jac(5,4)=-wold_l*gm1l
  674. wl.jac(5,5)=gm1l
  675. do 94 i=1,(nsp-1)
  676. wl.jac(5,5+i)=dgdyl.vv(i)*pold_l/(rold_l*gm1l)
  677. 94 continue
  678. c------------------------------
  679. do 95 i=1,(nsp-1)
  680. do 96 j=1,5
  681. wl.jac(5+i,j)=0.0d0
  682. if(j.eq.1) wl.jac(5+i,j)=-yl.yet(i)/rold_l
  683. 96 continue
  684. c---------------------
  685. do 960 j=6,(4+nsp)
  686. wl.jac(5+i,j)=0.0d0
  687. if(5+i.eq.j) then
  688. wl.jac(5+i,j)=1.0d0/rold_l
  689. endif
  690. 960 continue
  691. 95 continue
  692. c------------------------------
  693. wr.jac(1,1)=1.0d0
  694. wr.jac(1,2)=0.0d0
  695. wr.jac(1,3)=0.0d0
  696. wr.jac(1,4)=0.0d0
  697. wr.jac(1,5)=0.0d0
  698. do 97 i=1,(nsp-1)
  699. wr.jac(1,5+i)=0.0d0
  700. 97 continue
  701. c------------------------------
  702. wr.jac(2,1)=-uold_r/rold_r
  703. wr.jac(2,2)=1.0d0/rold_r
  704. wr.jac(2,3)=0.0d0
  705. wr.jac(2,4)=0.0d0
  706. wr.jac(2,5)=0.0d0
  707. do 98 i=1,(nsp-1)
  708. wr.jac(2,5+i)=0.0d0
  709. 98 continue
  710. c------------------------------
  711. wr.jac(3,1)=-vold_r/rold_r
  712. wr.jac(3,2)=0.0d0
  713. wr.jac(3,3)=1.0d0/rold_r
  714. wr.jac(3,4)=0.0d0
  715. wr.jac(3,5)=0.0d0
  716. do 99 i=1,(nsp-1)
  717. wr.jac(3,5+i)=0.0d0
  718. 99 continue
  719. c------------------------------
  720. wr.jac(4,1)=-wold_r/rold_r
  721. wr.jac(4,2)=0.0d0
  722. wr.jac(4,3)=0.0d0
  723. wr.jac(4,4)=1.0d0/rold_r
  724. wr.jac(4,5)=0.0d0
  725. do 100 i=1,(nsp-1)
  726. wr.jac(4,5+i)=0.0d0
  727. 100 continue
  728. c------------------------------
  729. br1=0.0d0
  730. do 101 i=1,(nsp-1)
  731. br1=br1+dgdyr.vv(i)*yr.yet(i)
  732. 101 continue
  733. br1=br1*pold_r/(rold_r*gm1r)
  734. wr.jac(5,1)=gm1r*(uold_r*uold_r+vold_r*vold_r+
  735. & wold_r*wold_r)/2.0d0
  736. wr.jac(5,1)=wr.jac(5,1)-br1
  737. wr.jac(5,2)=-uold_r*gm1r
  738. wr.jac(5,3)=-vold_r*gm1r
  739. wr.jac(5,4)=-wold_r*gm1r
  740. wr.jac(5,5)=gm1r
  741. do 102 i=1,(nsp-1)
  742. wr.jac(5,5+i)=dgdyr.vv(i)*pold_r/(rold_r*gm1r)
  743. 102 continue
  744. c----------------------------------------------
  745. do 103 i=1,(nsp-1)
  746. do 104 j=1,5
  747. wr.jac(5+i,j)=0.0d0
  748. if(j.eq.1) wr.jac(5+i,j)=-yr.yet(i)/rold_r
  749. 104 continue
  750. c---------------------
  751. do 1040 j=6,(4+nsp)
  752. wr.jac(5+i,j)=0.0d0
  753. if(5+i.eq.j) wr.jac(5+i,j)=1.0d0/rold_r
  754. 1040 continue
  755. 103 continue
  756. c----------------------------------------------
  757. SEGINI JTL, JTR
  758. c----------------------------------------------
  759. do 1 i=1,(4+nsp)
  760. do 2 j=1,(4+nsp)
  761. jtl.jac(i,j)=0.0d0
  762. jtr.jac(i,j)=0.0d0
  763. do 3 k=1,(4+nsp)
  764. jtl.jac(i,j)=jtl.jac(i,j)+jl.jac(i,k)*wl.jac(k,j)
  765. jtr.jac(i,j)=jtr.jac(i,j)+jr.jac(i,k)*wr.jac(k,j)
  766. 3 continue
  767. 2 continue
  768. 1 continue
  769. c----------------------------------------------------------------------
  770. c11=t1_y*t2_z - t1_z*t2_y
  771. c12=n_y*t2_z - n_z*t2_y
  772. c13=n_y*t1_z - n_z*t1_y
  773. c-------------------------------------
  774. c21=t1_x*t2_z - t1_z*t2_x
  775. c22=n_x*t2_z - n_z*t2_x
  776. c23=n_x*t1_z - n_z*t1_x
  777. c-------------------------------------
  778. c31=t1_x*t2_y - t1_y*t2_x
  779. c32=n_x*t2_y - n_y*t2_x
  780. c33=n_x*t1_y - n_y*t1_x
  781. det=n_x*c11 - n_y*c21 + n_z*c31
  782. c----------------------------------------------------------------------
  783. ZC11=-NVECT(1)*C11-TVEC1(1)*C12+TVEC2(1)*C13
  784. ZC12=-NVECT(2)*C11-TVEC1(2)*C12+TVEC2(2)*C13
  785. ZC13=-NVECT(3)*C11-TVEC1(3)*C12+TVEC2(3)*C13
  786. C---------------------------------
  787. ZC21=NVECT(1)*C21+TVEC1(1)*C22-TVEC2(1)*C23
  788. ZC22=NVECT(2)*C21+TVEC1(2)*C22-TVEC2(2)*C23
  789. ZC23=NVECT(3)*C21+TVEC1(3)*C22-TVEC2(3)*C23
  790. C---------------------------------
  791. ZC31=-NVECT(1)*C31-TVEC1(1)*C32+TVEC2(1)*C33
  792. ZC32=-NVECT(2)*C31-TVEC1(2)*C32+TVEC2(2)*C33
  793. ZC33=-NVECT(3)*C31-TVEC1(3)*C32+TVEC2(3)*C33
  794. c----------------------------------------------------------------------
  795. do 11 i=1,(4+nsp)
  796. jtl.jac(i,1)=jtl.jac(i,1)+jtr.jac(i,1)
  797. jtl.jac(i,2)=jtl.jac(i,2)+jtr.jac(i,2)*zc11/det+
  798. & jtr.jac(i,3)*zc21/det+jtr.jac(i,4)*zc31/det
  799. jtl.jac(i,3)=jtl.jac(i,3)+jtr.jac(i,2)*zc12/det+
  800. & jtr.jac(i,3)*zc22/det+jtr.jac(i,4)*zc32/det
  801. jtl.jac(i,4)=jtl.jac(i,4)+jtr.jac(i,2)*zc13/det+
  802. & jtr.jac(i,3)*zc23/det+jtr.jac(i,4)*zc33/det
  803. jtl.jac(i,5)=jtl.jac(i,5)+jtr.jac(i,5)
  804. do 777 j=6,(4+nsp)
  805. jtl.jac(i,j)=jtl.jac(i,j)+jtr.jac(i,j)
  806. 777 continue
  807. 11 continue
  808. c----------------------------------------------------------------------
  809. SEGDES JL
  810. SEGDES JR
  811. SEGDES WL
  812. SEGDES WR
  813. SEGDES JTR
  814. jll=jtl
  815. SEGDES JTL
  816. SEGDES YL
  817. SEGDES YR
  818. SEGDES CP
  819. SEGDES CV
  820. SEGDES MLRECP, MLRECV
  821. C----------------------------------------------------------------------
  822. SEGDES dmly_l, dmly_r,
  823. & dmry_l, dmry_r,
  824. & dPpy_l, dPpy_r,
  825. & dPmy_l, dPmy_r,
  826. & dpiy_l, dpiy_r,
  827. & dgdyl, dgdyr
  828. C----------------------------------------------------------------------
  829. return
  830. end
  831.  
  832.  
  833.  
  834.  
  835.  
  836.  
  837.  
  838.  
  839.  
  840.  
  841.  
  842.  
  843.  
  844.  
  845.  
  846.  
  847.  
  848.  
  849.  
  850.  
  851.  
  852.  
  853.  
  854.  
  855.  
  856.  
  857.  
  858.  
  859.  
  860.  
  861.  
  862.  
  863.  
  864.  

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