1 from numpy
import sqrt, zeros
4 def find_speed_fast (qvar, gamma):
7 d=qvar[0]; P=qvar[1]; A=qvar[3]; B=qvar[5]; C=qvar[7]
11 cf = sqrt(d2 + sqrt(max(d2*d2-c2*A*A/d,zero)))
16 def hlld (qleft, qright, gamma=1.4, dbg=0):
23 A=half*(qleft[3]+qright[3])
25 qleft[3]=A; qright[3]=A
28 rl=qleft[0]; Pl=qleft[1]; ul=qleft[2]
29 vl=qleft[4]; Bl=qleft[5]; wl=qleft[6]
31 ecinl = half*(ul*ul+vl*vl+wl*wl)*rl
32 emagl = half*(A*A+Bl*Bl+Cl*Cl)
33 etotl = Pl*entho+ecinl+emagl
35 vdotBl= ul*A+vl*Bl+wl*Cl
38 rr=qright[0]; Pr=qright[1]; ur=qright[2]
39 vr=qright[4]; Br=qright[5]; wr=qright[6]; Cr=qright[7]
40 ecinr = half*(ur*ur+vr*vr+wr*wr)*rr
41 emagr = half*(A*A+Br*Br+Cr*Cr)
42 etotr = Pr*entho+ecinr+emagr
44 vdotBr= ur*A+vr*Br+wr*Cr
47 cfastl = find_speed_fast (qleft, gamma)
48 cfastr = find_speed_fast (qright, gamma)
51 SL=min([ul,ur])-max([cfastl,cfastr])
52 SR=max([ul,ur])+max([cfastl,cfastr])
61 ustar =(rcr*ur +rcl*ul + (Ptotl-Ptotr))/(rcr+rcl)
62 Ptotstar=(rcr*Ptotl+rcl*Ptotr+rcl*rcr*(ul-ur))/(rcr+rcl)
65 rstarl=rl*(SL-ul)/(SL-ustar)
66 estar =rl*(SL-ul)*(SL-ustar)-A**2
67 el =rl*(SL-ul)*(SL-ul )-A**2
75 vstarl=vl-A*Bl*(ustar-ul)/estar
77 wstarl=wl-A*Cl*(ustar-ul)/estar
80 vdotBstarl=ustar*A+vstarl*Bstarl+wstarl*Cstarl
81 etotstarl=((SL-ul)*etotl-Ptotl*ul+Ptotstar*ustar+A*(vdotBl-vdotBstarl))/(SL-ustar)
82 sqrrstarl=sqrt(rstarl)
83 calfvenl=abs(A)/sqrrstarl
87 rstarr=rr*(SR-ur)/(SR-ustar)
88 estar =rr*(SR-ur)*(SR-ustar)-A**2
89 er =rr*(SR-ur)*(SR-ur )-A**2
96 vstarr=vr-A*Br*(ustar-ur)/estar
98 wstarr=wr-A*Cr*(ustar-ur)/estar
102 vdotBstarr=ustar*A+vstarr*Bstarr+wstarr*Cstarr
104 etotstarr=((SR-ur)*etotr-Ptotr*ur+Ptotstar*ustar+A*(vdotBr-vdotBstarr))/(SR-ustar)
106 sqrrstarr=sqrt(rstarr)
107 calfvenr=abs(A)/sqrrstarr
111 vstarstar=(sqrrstarl*vstarl+sqrrstarr*vstarr+sgnm*(Bstarr-Bstarl))/(sqrrstarl+sqrrstarr)
112 wstarstar=(sqrrstarl*wstarl+sqrrstarr*wstarr+sgnm*(Cstarr-Cstarl))/(sqrrstarl+sqrrstarr)
113 Bstarstar=(sqrrstarl*Bstarr+sqrrstarr*Bstarl+sgnm*sqrrstarl*sqrrstarr*(vstarr-vstarl))/(sqrrstarl+sqrrstarr)
114 Cstarstar=(sqrrstarl*Cstarr+sqrrstarr*Cstarl+sgnm*sqrrstarl*sqrrstarr*(wstarr-wstarl))/(sqrrstarl+sqrrstarr)
115 vdotBstarstar=ustar*A+vstarstar*Bstarstar+wstarstar*Cstarstar
116 etotstarstarl=etotstarl-sgnm*sqrrstarl*(vdotBstarl-vdotBstarstar)
117 etotstarstarr=etotstarr+sgnm*sqrrstarr*(vdotBstarr-vdotBstarstar)
141 emago=half*(A*A+Bo*Bo+Co*Co)
152 emago=half*(A*A+Bo*Bo+Co*Co)
164 emago=half*(A*A+Bo*Bo+Co*Co)
175 emago=half*(A*A+Bo*Bo+Co*Co)
190 fgdnv = zeros(len(qleft))
192 fgdnv[1] = (etoto+Ptoto)*uo-A*vdotBo
193 fgdnv[2] = ro*uo*uo+Ptoto+emago*(lim-1.0)-A*A*lim
195 fgdnv[4] = ro*uo*vo-A*Bo*lim
196 fgdnv[5] = Bo*uo-A*vo
197 fgdnv[6] = ro*uo*wo-A*Co*lim
198 fgdnv[7] = Co*uo-A*wo
201 for ivar
in (9,nvar):
203 fgdnv[ivar] = fgdnv[0]*qleft [ivar]
205 fgdnv[ivar] = fgdnv[0]*qright[ivar]