Télécharger ausmw.eso

Retour à la liste

Numérotation des lignes :

  1. C AUSMW SOURCE CB215821 16/04/21 21:15:10 8920
  2. c----------------------------------------------------------------------
  3. c INPUT:
  4. c
  5. c rog : left density
  6. C
  7. c uxg : left x-component of speed
  8. c
  9. c uyg : left y-component of speed
  10. c
  11. c pg : left pressure
  12. c
  13. c rod : right density
  14. c
  15. c uxd : right x-component of speed
  16. c
  17. c uyd : right y-component of speed
  18. c
  19. c pd : right pressure
  20. c
  21. c ga : gamma of the gas
  22. c
  23. c cnx : x-component of the surface normal
  24. c
  25. c cny : y-component of the surface normal
  26. c
  27. c OUTPUT:
  28. c
  29. c dfrun : derivative of the numerical estimation of (rho u_n u_n +
  30. c p) (flux-function for conservative variable (rho u_n))
  31. c with respect to the primitive variables
  32. c
  33. c----------------------------------------------------------------------
  34. SUBROUTINE AUSMW(ROG,UXG,UYG,PG,ROD,UXD,UYD,PD,
  35. & ga,cnx,cny,dfrun)
  36. IMPLICIT INTEGER(I-N)
  37. real*8 rog,uxg,uyg,pg
  38. real*8 rod,uxd,uyd,pd
  39. real*8 cnx,cny
  40. real*8 jl(4,4),jr(4,4),dfrun(4)
  41. real*8 alpha,beta
  42. real*8 ga, gm1
  43. real*8 n_x,n_y,t_x,t_y
  44. real*8 un_l, un_r
  45. real*8 ml,mr,Mplus,Mmin
  46. real*8 rold_l,uold_l,vold_l,pold_l,eold_l
  47. real*8 rold_r,uold_r,vold_r,pold_r,eold_r
  48. real*8 Pplus,Pmin
  49. real*8 temp_l,temp_r,brac_l,brac_r
  50. real*8 aleft, arigh, am
  51. real*8 damr_l,damr_r,damu_l,damu_r
  52. real*8 damv_l,damv_r,damp_l,damp_r
  53. real*8 dmlr_l,dmlr_r,dmlu_l,dmlu_r
  54. real*8 dmlv_l,dmlv_r,dmlp_l,dmlp_r
  55. real*8 dmrr_l,dmrr_r,dmru_l,dmru_r
  56. real*8 dmrv_l,dmrv_r,dmrp_l,dmrp_r
  57. real*8 dPpr_l,dPpr_r,dPpu_l,dPpu_r
  58. real*8 dPpv_l,dPpv_r,dPpp_l,dPpp_r
  59. real*8 dPmr_l,dPmr_r,dPmu_l,dPmu_r
  60. real*8 dPmv_l,dPmv_r,dPmp_l,dPmp_r
  61. real*8 dpir_l,dpir_r,dpiu_l,dpiu_r
  62. real*8 dpiv_l,dpiv_r,dpip_l,dpip_r
  63. real*8 tcoef,bcoef
  64. integer i
  65. parameter(alpha = 0.1875D0, beta = 0.125D0)
  66. c-----------------------
  67. gm1=ga-1.0d0
  68. c-----------------------
  69. n_x=cnx
  70. n_y=cny
  71. t_x=-cny
  72. t_y=cnx
  73. c----------------------------
  74. c-----------------------
  75. rold_l=rog
  76. uold_l=uxg
  77. vold_l=uyg
  78. pold_l=pg
  79. c-----------------------
  80. rold_r=rod
  81. uold_r=uxd
  82. vold_r=uyd
  83. pold_r=pd
  84. c-----------------------
  85. eold_l=(uold_l*uold_l+vold_l*vold_l)/2.0d0
  86. eold_l=eold_l+pold_l/(gm1*rold_l)
  87. eold_r=(uold_r*uold_r+vold_r*vold_r)/2.0d0
  88. eold_r=eold_r+pold_r/(gm1*rold_r)
  89. c-------------------------------------------------------------------
  90. c Computation of the speed of sound and its derivatives
  91. c---------------------------------------------------------------------
  92. aleft=SQRT(ga*pold_l/rold_l)
  93. arigh=SQRT(ga*pold_r/rold_r)
  94. am=0.5d0*(aleft+arigh)
  95. c--------------------------------------------------------------------
  96. damr_r=-arigh/(4.0d0*rold_r)
  97. damu_r=0.0d0
  98. damv_r=0.0d0
  99. damp_r=ga/(4.0d0*arigh*rold_r)
  100. c-----------------------
  101. damr_l=-aleft/(4.0d0*rold_l)
  102. damu_l=0.0d0
  103. damv_l=0.0d0
  104. damp_l=ga/(4.0d0*aleft*rold_l)
  105. c----------------------------------------------------------------------
  106. c computing numerical Mach number and its derivatives
  107. c----------------------------------------------------------------------
  108. un_l=uold_l*n_x+vold_l*n_y
  109. un_r=uold_r*n_x+vold_r*n_y
  110. c----------------------------------------
  111. ml=un_l/am
  112. mr=un_r/am
  113. c-------------------------------------
  114. c Mplus and Mmin are calligraphic
  115. c lettes M+ and M- from the paper
  116. c-------------------------------------
  117. if(ABS(ml) .ge. 1.0d0) then
  118. Mplus=(ml+ABS(ml))/2.0d0
  119. else
  120. Mplus=(ml+1.0d0)*(ml+1.0d0)/4.0d0
  121. Mplus=Mplus+beta*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  122. endif
  123. c-----------------
  124. if(ABS(mr) .ge. 1.0d0) then
  125. Mmin=(mr-ABS(mr))/2.0d0
  126. else
  127. Mmin=-(mr-1.0d0)*(mr-1.0d0)/4.0d0
  128. Mmin=Mmin-beta*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  129. endif
  130. c---------------------------------------
  131. c Derivatives of ml and mr with respect
  132. c to both sets of primitive variables:
  133. c left and right
  134. c--------------------------------------
  135. temp_l=-un_l/(am*am)
  136. temp_r=-un_r/(am*am)
  137. c--------
  138. dmlr_l=temp_l*damr_l
  139. dmlr_r=temp_l*damr_r
  140. dmrr_l=temp_r*damr_l
  141. dmrr_r=temp_r*damr_r
  142. c--------
  143. dmlu_l=n_x/am+temp_l*damu_l
  144. dmlu_r=temp_l*damu_r
  145. dmru_l=temp_r*damu_l
  146. dmru_r=n_x/am+temp_r*damu_r
  147. c--------
  148. dmlv_l=n_y/am+temp_l*damv_l
  149. dmlv_r=temp_l*damv_r
  150. dmrv_l=temp_r*damv_l
  151. dmrv_r=n_y/am+temp_r*damv_r
  152. c--------
  153. dmlp_l=temp_l*damp_l
  154. dmlp_r=temp_l*damp_r
  155. dmrp_l=temp_r*damp_l
  156. dmrp_r=temp_r*damp_r
  157. c---------------------------------------------------------------
  158. c Computing the calligraphic P+ and P- with their derivatives
  159. c---------------------------------------------------------------
  160. if(ml .ge. 1.0d0) then
  161. Pplus = 1.0d0
  162. else
  163. if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  164. Pplus=(ml+1.0d0)*(ml+1.0d0)*(2.0d0-ml)/4.0d0
  165. Pplus=Pplus+alpha*ml*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  166. else
  167. Pplus = 0.0d0
  168. endif
  169. endif
  170. c---------------------------------------------------------------
  171. if(mr .ge. 1.0d0) then
  172. Pmin = 0.0d0
  173. else
  174. if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  175. Pmin=(mr-1.0d0)*(mr-1.0d0)*(2.0d0+mr)/4.0d0
  176. Pmin=Pmin-alpha*mr*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  177. else
  178. Pmin = 1.0d0
  179. endif
  180. endif
  181. c---------------------------------------------------------------
  182. brac_l=(ml+1.0d0)*(2.0d0-ml)/2.0d0-(ml+1.0d0)*(ml+1.0d0)/4.0d0
  183. brac_l=brac_l+alpha*(ml*ml-1.0d0)*(ml*ml-1.0d0)
  184. brac_l=brac_l+4.0d0*alpha*ml*ml*(ml*ml-1.0d0)
  185. c--------------
  186. brac_r=(mr-1.0d0)*(2.0d0+mr)/2.0d0+(mr-1.0d0)*(mr-1.0d0)/4.0d0
  187. brac_r=brac_r-alpha*(mr*mr-1.0d0)*(mr*mr-1.0d0)
  188. brac_r=brac_r-4.0d0*alpha*mr*mr*(mr*mr-1.0d0)
  189. c---------------------------------------------------------------
  190. if((ml .gt. -1.0d0) .and. (ml .lt. 1.0d0)) then
  191. dPpr_l=brac_l*dmlr_l
  192. dPpr_r=brac_l*dmlr_r
  193. c------------
  194. dPpu_l=brac_l*dmlu_l
  195. dPpu_r=brac_l*dmlu_r
  196. c------------
  197. dPpv_l=brac_l*dmlv_l
  198. dPpv_r=brac_l*dmlv_r
  199. c------------
  200. dPpp_l=brac_l*dmlp_l
  201. dPpp_r=brac_l*dmlp_r
  202. c------------
  203. else
  204. dPpr_l=0.0d0
  205. dPpr_r=0.0d0
  206. c-----------
  207. dPpu_l=0.0d0
  208. dPpu_r=0.0d0
  209. c-----------
  210. dPpv_l=0.0d0
  211. dPpv_r=0.0d0
  212. c-----------
  213. dPpp_l=0.0d0
  214. dPpp_r=0.0d0
  215. c-----------
  216. endif
  217. c---------------------------------------------------------------
  218. if((mr .gt. -1.0d0) .and. (mr .lt. 1.0d0)) then
  219. dPmr_l=brac_r*dmrr_l
  220. dPmr_r=brac_r*dmrr_r
  221. c------------
  222. dPmu_l=brac_r*dmru_l
  223. dPmu_r=brac_r*dmru_r
  224. c------------
  225. dPmv_l=brac_r*dmrv_l
  226. dPmv_r=brac_r*dmrv_r
  227. c------------
  228. dPmp_l=brac_r*dmrp_l
  229. dPmp_r=brac_r*dmrp_r
  230. c------------
  231. else
  232. dPmr_l=0.0d0
  233. dPmr_r=0.0d0
  234. c-----------
  235. dPmu_l=0.0d0
  236. dPmu_r=0.0d0
  237. c-----------
  238. dPmv_l=0.0d0
  239. dPmv_r=0.0d0
  240. c-----------
  241. dPmp_l=0.0d0
  242. dPmp_r=0.0d0
  243. c-----------
  244. endif
  245. c---------------------------------------------------------------------
  246. c computing pmid -- p_{1/2} and its derivatives
  247. c---------------------------------------------------------------------
  248. dpir_l=dPpr_l*pold_l+dPmr_l*pold_r
  249. dpiu_l=dPpu_l*pold_l+dPmu_l*pold_r
  250. dpiv_l=dPpv_l*pold_l+dPmv_l*pold_r
  251. dpip_l=dPpp_l*pold_l+Pplus+dPmp_l*pold_r
  252. c----------------------------
  253. dpir_r=dPpr_r*pold_l+dPmr_r*pold_r
  254. dpiu_r=dPpu_r*pold_l+dPmu_r*pold_r
  255. dpiv_r=dPpv_r*pold_l+dPmv_r*pold_r
  256. dpip_r=dPpp_r*pold_l+Pmin+dPmp_r*pold_r
  257. c---------------------------------------------------------------------
  258. c !!!!!!!!!!!!!!!!!! JACOBIAN !!!!!!!!!!!!!!!!!!!!!!!!!!!!
  259. c---------------------------------------------------------------------
  260. jl(1,1)=0.0D0
  261. jl(1,2)=0.0D0
  262. jl(1,3)=0.0D0
  263. jl(1,4)=0.0D0
  264. c------------------------------------
  265. jr(1,1)=0.0D0
  266. jr(1,2)=0.0D0
  267. jr(1,3)=0.0D0
  268. jr(1,4)=0.0D0
  269. c------------------------------------
  270. ccccc f222222222222 -------------
  271. c------------------------------------
  272. c---------------------------------------------------------
  273. c--------------------------------------------
  274. jl(2,1)=n_x*dpir_l
  275. c-------------------
  276. jl(2,2)=n_x*dpiu_l
  277. c-------------------
  278. jl(2,3)=n_x*dpiv_l
  279. c-------------------
  280. jl(2,4)=n_x*dpip_l
  281. c-------------------------------------------------------------
  282. jr(2,1)=n_x*dpir_r
  283. c-------------------
  284. jr(2,2)=n_x*dpiu_r
  285. c-------------------
  286. jr(2,3)=n_x*dpiv_r
  287. c-------------------
  288. jr(2,4)=n_x*dpip_r
  289. c-------------------------------------------------------------
  290. c------------ f33333333333333333333 ---------------------
  291. c-------------------------------------------------------------
  292. jl(3,1)=n_y*dpir_l
  293. c-------------------
  294. jl(3,2)=n_y*dpiu_l
  295. c-------------------
  296. jl(3,3)=n_y*dpiv_l
  297. c-------------------
  298. jl(3,4)=n_y*dpip_l
  299. c-------------------------------------------------------------
  300. jr(3,1)=n_y*dpir_r
  301. c-------------------
  302. jr(3,2)=n_y*dpiu_r
  303. c-------------------
  304. jr(3,3)=n_y*dpiv_r
  305. c-------------------
  306. jr(3,4)=n_y*dpip_r
  307. c-------------------------------------------------------------
  308. c ------ f44444444444444444444444444444 ---------
  309. c-------------------------------------------------------------
  310. jl(4,1)=0.0D0
  311. c---------------------
  312. jl(4,2)=0.0D0
  313. c---------------------
  314. jl(4,3)=0.0D0
  315. c---------------------
  316. jl(4,4)=0.0D0
  317. c----------------------------------------------------------
  318. c----------------------------------------------------------
  319. jr(4,1)=0.0D0
  320. c---------------------
  321. jr(4,2)=0.0D0
  322. c---------------------
  323. jr(4,3)=0.0D0
  324. c---------------------
  325. jr(4,4)=0.0D0
  326. c----------------------------------------------------------
  327. tcoef=n_x*t_y+t_x*n_y
  328. bcoef=n_x*t_y-t_x*n_y
  329. c-----------------------------
  330. c-----------------------------------
  331. do 11 i=1,4
  332. jl(i,1)=jl(i,1)+jr(i,1)
  333. jl(i,2)=jl(i,2)+jr(i,2)*(-tcoef/bcoef)+
  334. & jr(i,3)*2.0d0*n_x*t_x/bcoef
  335. jl(i,3)=jl(i,3)+jr(i,2)*(-2.0d0*n_y*t_y/bcoef)+
  336. & jr(i,3)*tcoef/bcoef
  337. jl(i,4)=jl(i,4)+jr(i,4)
  338. 11 continue
  339. c--------------------------------------
  340. do 12 i=1,4
  341. dfrun(i)=jl(2,i)*n_x+jl(3,i)*n_y
  342. 12 continue
  343. c----------------------------------------------------------------------
  344. return
  345. end
  346.  
  347.  
  348.  
  349.  
  350.  
  351.  
  352.  
  353.  
  354.  
  355.  
  356.  
  357.  
  358.  
  359.  
  360.  
  361.  
  362.  
  363.  
  364.  
  365.  
  366.  
  367.  

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