Télécharger cwmsp.eso

Retour à la liste

Numérotation des lignes :

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

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