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(nsp8sPan style="color: #009900;">)+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)=äpiv_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
  • wr.jac(1,2)=0.0d0
  • wr.jac(1,3)=0.0d0
  • wr.jac(1,4)=0.0d0
  • wr.jac(1,5)=0.0d0
  • do 97 i=1,(nsp-1)
  • wr.jac(1,5+i)=0.0d0
  • 97 continue
  • c------------------------------
  • wr.jac(2,1)=-uold_r/rold_r
  • wr.jac(2,2)=1.0d0/rold_r
  • wr.jac(2,3)=0.0d0
  • wr.jac(2,4)=0.0d0
  • wr.jac(2,5)=0.0d0
  • do 98 i=1,(nsp-1)
  • wr.jac(2,5+i)=0.0d0
  • 98 continue
  • c------------------------------
  • wr.jac(3,1)=-vold_r/rold_r
  • wr.jac(3,2)=0.0d0
  • wr.jac(3,3)=1.0d0/rold_r
  • wr.jac(3,4)=0.0d0
  • wr.jac(3,5)=0.0d0
  • do 99 i=1,(nsp-1)
  • wr.jac(3,5+i)=0.0d0
  • 99 continue
  • c------------------------------
  • wr.jac(4,1)=-wold_r/rold_r
  • wr.jac(4,2)=0.0d0
  • wr.jac(4,3)=0.0d0
  • wr.jac(4,4)=1.0d0/rold_r
  • wr.jac(4,5)=0.0d0
  • do 100 i=1,(nsp-1)
  • wr.jac(4,5+i)=0.0d0
  • 100 continue
  • c------------------------------
  • br1=0.0d0
  • do 101 i=1,(nsp-1)
  • br1=br1+dgdyr.vv(i)*yr.yet(i)
  • 101 continue
  • br1=br1*pold_r/(rold_r*gm1r)
  • wr.jac(5,1)=gm1r*(uold_r*uold_r+vold_r*vold_r+
  • & wold_r*wold_r)/2.0d0
  • wr.jac(5,1)=wr.jac(5,1)-br1
  • wr.jac(5,2)=-uold_r*gm1r
  • wr.jac(5,3)=-vold_r*gm1r
  • wr.jac(5,4)=-wold_r*gm1r
  • wr.jac(5,5)=gm1r
  • do 102 i=1,(nsp-1)
  • wr.jac(5,5+i)=dgdyr.vv(i)*pold_r/(rold_r*gm1r)
  • 102 continue
  • c----------------------------------------------
  • do 103 i=1,(nsp-1)
  • do 104 j=1,5
  • wr.jac(5+i,j)=0.0d0
  • if(j.eq.1) wr.jac(5+i,j)=-yr.yet(i)/rold_r
  • 104 continue
  • c---------------------
  • do 1040 j=6,(4+nsp)
  • wr.jac(5+i,j)=0.0d0
  • if(5+i.eq.j) wr.jac(5+i,j)=1.0d0/rold_r
  • 1040 continue
  • 103 continue
  • c----------------------------------------------
  • SEGINI JTL, JTR
  • c----------------------------------------------
  • do 1 i=1,(4+nsp)
  • do 2 j=1,(4+nsp)
  • jtl.jac(i,j)=0.0d0
  • jtr.jac(i,j)=0.0d0
  • do 3 k=1,(4+nsp)
  • jtl.jac(i,j)=jtl.jac(i,j)+jl.jac(i,k)*wl.jac(k,j)
  • jtr.jac(i,j)=jtr.jac(i,j)+jr.jac(i,k)*wr.jac(k,j)
  • 3 continue
  • 2 continue
  • 1 continue
  • c----------------------------------------------------------------------
  • c11=t1_y*t2_z - t1_z*t2_y
  • c12=n_y*t2_z - n_z*t2_y
  • c13=n_y*t1_z - n_z*t1_y
  • c-------------------------------------
  • c21=t1_x*t2_z - t1_z*t2_x
  • c22=n_x*t2_z - n_z*t2_x
  • c23=n_x*t1_z - n_z*t1_x
  • c-------------------------------------
  • c31=t1_x*t2_y - t1_y*t2_x
  • c32=n_x*t2_y - n_y*t2_x
  • c33=n_x*t1_y - n_y*t1_x
  • det=n_x*c11 - n_y*c21 + n_z*c31
  • c----------------------------------------------------------------------
  • ZC11=-NVECT(1)*C11-TVEC1(1)*C12+TVEC2(1)*C13
  • ZC12=-NVECT(2)*C11-TVEC1(2)*C12+TVEC2(2)*C13
  • ZC13=-NVECT(3)*C11-TVEC1(3)*C12+TVEC2(3)*C13
  • C---------------------------------
  • ZC21=NVECT(1)*C21+TVEC1(1)*C22-TVEC2(1)*C23
  • ZC22=NVECT(2)*C21+TVEC1(2)*C22-TVEC2(2)*C23
  • ZC23=NVECT(3)*C21+TVEC1(3)*C22-TVEC2(3)*C23
  • C---------------------------------
  • ZC31=-NVECT(1)*C31-TVEC1(1)*C32+TVEC2(1)*C33
  • ZC32=-NVECT(2)*C31-TVEC1(2)*C32+TVEC2(2)*C33
  • ZC33=-NVECT(3)*C31-TVEC1(3)*C32+TVEC2(3)*C33
  • c----------------------------------------------------------------------
  • do 11 i=1,(4+nsp)
  • jtl.jac(i,1)=jtl.jac(i,1)+jtr.jac(i,1)
  • jtl.jac(i,2)=jtl.jac(i,2)+jtr.jac(i,2)*zc11/det+
  • & jtr.jac(i,3)*zc21/det+jtr.jac(i,4)*zc31/det
  • jtl.jac(i,3)=jtl.jac(i,3)+jtr.jac(i,2)*zc12/det+
  • & jtr.jac(i,3)*zc22/det+jtr.jac(i,4)*zc32/det
  • jtl.jac(i,4)=jtl.jac(i,4)+jtr.jac(i,2)*zc13/det+
  • & jtr.jac(i,3)*zc23/det+jtr.jac(i,4)*zc33/det
  • jtl.jac(i,5)=jtl.jac(i,5)+jtr.jac(i,5)
  • do 777 j=6,(4+nsp)
  • jtl.jac(i,j)=jtl.jac(i,j)+jtr.jac(i,j)
  • 777 continue
  • 11 continue
  • c----------------------------------------------------------------------
  • SEGDES JL
  • SEGDES JR
  • SEGDES WL
  • SEGDES WR
  • SEGDES JTR
  • jll=jtl
  • SEGDES JTL
  • SEGDES YL
  • SEGDES YR
  • SEGDES CP
  • SEGDES CV
  • SEGDES MLRECP, MLRECV
  • C----------------------------------------------------------------------
  • SEGDES dmly_l, dmly_r,
  • & dmry_l, dmry_r,
  • & dPpy_l, dPpy_r,
  • & dPmy_l, dPmy_r,
  • & dpiy_l, dpiy_r,
  • & dgdyl, dgdyr
  • C----------------------------------------------------------------------
  • return
  • end
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  •  
  • © Cast3M 2003 - Tous droits réservés.
    Mentions légales