DISPATCH
hlld.py
1 from numpy import sqrt, zeros
2 
3 #-------------------------------------------------------------------------------
4 def find_speed_fast (qvar, gamma):
5  half=0.5
6  zero=0.0
7  d=qvar[0]; P=qvar[1]; A=qvar[3]; B=qvar[5]; C=qvar[7]
8  B2 = A*A+B*B+C*C
9  c2 = gamma*P/d
10  d2 = half*(B2/d+c2)
11  cf = sqrt(d2 + sqrt(max(d2*d2-c2*A*A/d,zero)))
12  vel_info = cf
13  return vel_info
14 
15 #-------------------------------------------------------------------------------
16 def hlld (qleft, qright, gamma=1.4, dbg=0):
17  one=1.0
18  entho=one/(gamma-one)
19  half=0.5
20  zero=0.0
21 
22  # Enforce continuity of normal component
23  A=half*(qleft[3]+qright[3])
24  sgnm=A/(abs(A)+1e-60)
25  qleft[3]=A; qright[3]=A
26 
27  # Left variables
28  rl=qleft[0]; Pl=qleft[1]; ul=qleft[2]
29  vl=qleft[4]; Bl=qleft[5]; wl=qleft[6]
30  Cl=qleft[7]
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
34  Ptotl = Pl + emagl
35  vdotBl= ul*A+vl*Bl+wl*Cl
36 
37  # Right variables
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
43  Ptotr = Pr + emagr
44  vdotBr= ur*A+vr*Br+wr*Cr
45 
46  # Find the largest eigenvalues in the normal direction to the interface
47  cfastl = find_speed_fast (qleft, gamma)
48  cfastr = find_speed_fast (qright, gamma)
49 
50  # Compute HLL wave speed
51  SL=min([ul,ur])-max([cfastl,cfastr])
52  SR=max([ul,ur])+max([cfastl,cfastr])
53 
54  lim = 1.0
55 
56  # Compute lagrangian sound speed
57  rcl=rl*(ul-SL)
58  rcr=rr*(SR-ur)
59 
60  # Compute acoustic star state
61  ustar =(rcr*ur +rcl*ul + (Ptotl-Ptotr))/(rcr+rcl)
62  Ptotstar=(rcr*Ptotl+rcl*Ptotr+rcl*rcr*(ul-ur))/(rcr+rcl)
63 
64  # Left star region variables
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
68 
69  if (estar == 0):
70  vstarl=vl
71  Bstarl=Bl
72  wstarl=wl
73  Cstarl=Cl
74  else:
75  vstarl=vl-A*Bl*(ustar-ul)/estar
76  Bstarl=Bl*el/estar
77  wstarl=wl-A*Cl*(ustar-ul)/estar
78  Cstarl=Cl*el/estar
79 
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
84  SAL=ustar-calfvenl
85 
86  # Right star region variables
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
90  if (estar <= zero):
91  vstarr=vr
92  Bstarr=Br
93  wstarr=wr
94  Cstarr=Cr
95  else:
96  vstarr=vr-A*Br*(ustar-ur)/estar
97  Bstarr=Br*er/estar
98  wstarr=wr-A*Cr*(ustar-ur)/estar
99  Cstarr=Cr*er/estar
100 
101 
102  vdotBstarr=ustar*A+vstarr*Bstarr+wstarr*Cstarr
103 
104  etotstarr=((SR-ur)*etotr-Ptotr*ur+Ptotstar*ustar+A*(vdotBr-vdotBstarr))/(SR-ustar)
105 
106  sqrrstarr=sqrt(rstarr)
107  calfvenr=abs(A)/sqrrstarr
108  SAR=ustar+calfvenr
109 
110  # Double star region variables
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)
118 
119  # Sample the solution at x/t=0
120  if (SL > 0.0):
121  ro=rl
122  uo=ul
123  vo=vl
124  wo=wl
125  Bo=Bl
126  Co=Cl
127  Ptoto=Ptotl
128  etoto=etotl
129  vdotBo=vdotBl
130  emago=emagl
131  elif (SAL > 0.0):
132  ro=rstarl
133  uo=ustar
134  vo=vstarl
135  wo=wstarl
136  Bo=Bstarl
137  Co=Cstarl
138  Ptoto=Ptotstar
139  etoto=etotstarl
140  vdotBo=vdotBstarl
141  emago=half*(A*A+Bo*Bo+Co*Co)
142  elif (ustar > 0.0):
143  ro=rstarl
144  uo=ustar
145  vo=vstarstar
146  wo=wstarstar
147  Bo=Bstarstar
148  Co=Cstarstar
149  Ptoto=Ptotstar
150  etoto=etotstarstarl
151  vdotBo=vdotBstarstar
152  emago=half*(A*A+Bo*Bo+Co*Co)
153  elif (SAR > 0.0):
154  ro=rstarr
155  uo=ustar
156  vo=vstarstar
157  wo=wstarstar
158  Bo=Bstarstar
159  Co=Cstarstar
160  Ptoto=Ptotstar
161  etoto=etotstarstarr
162  help,etotstarstarr
163  vdotBo=vdotBstarstar
164  emago=half*(A*A+Bo*Bo+Co*Co)
165  elif (SR > 0.0):
166  ro=rstarr
167  uo=ustar
168  vo=vstarr
169  wo=wstarr
170  Bo=Bstarr
171  Co=Cstarr
172  Ptoto=Ptotstar
173  etoto=etotstarr
174  vdotBo=vdotBstarr
175  emago=half*(A*A+Bo*Bo+Co*Co)
176  else:
177  ro=rr
178  uo=ur
179  vo=vr
180  wo=wr
181  Bo=Br
182  Co=Cr
183  Ptoto=Ptotr
184  etoto=etotr
185  vdotBo=vdotBr
186  emago=emagr
187 
188 
189  # Compute the Godunov flux
190  fgdnv = zeros(len(qleft))
191  fgdnv[0] = ro*uo
192  fgdnv[1] = (etoto+Ptoto)*uo-A*vdotBo
193  fgdnv[2] = ro*uo*uo+Ptoto+emago*(lim-1.0)-A*A*lim
194  fgdnv[3] = zero
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
199  nvar = len(qleft)
200  if nvar > 8:
201  for ivar in (9,nvar):
202  if fgdnv[0] > 0:
203  fgdnv[ivar] = fgdnv[0]*qleft [ivar]
204  else:
205  fgdnv[ivar] = fgdnv[0]*qright[ivar]
206 
207  return fgdnv
Definition: hlld.py:1