@@ -23,6 +23,7 @@ def f(p, U_ij, gamma, idens, ixmom, iymom, iener):
2323
2424 return (gamma - 1.0 ) * (tau + D * (1.0 - W ) + p * (1.0 - W ** 2 )) / W ** 2 - p
2525
26+
2627@njit (cache = True )
2728def brentq (x1 , b , U , gamma , idens , ixmom , iymom , iener ,
2829 TOL = 1.e-6 , ITMAX = 100 ):
@@ -41,13 +42,13 @@ def brentq(x1, b, U, gamma, idens, ixmom, iymom, iener,
4142 # root found
4243 if fa * fb >= 0.0 :
4344 return x1
44-
45+
4546 # switch variables
4647 if abs (fa ) < abs (fb ):
4748 d = a
4849 a = b
4950 b = d
50-
51+
5152 d = fa
5253 fa = fb
5354 fb = d
@@ -70,7 +71,7 @@ def brentq(x1, b, U, gamma, idens, ixmom, iymom, iener,
7071 if 0.25 * (3.0 * a + b ) < b :
7172 if s < 0.25 * (3.0 * a + b ) or s > b :
7273 con1 = True
73- elif s < b or s > 0.25 * (3.0 * a + b ):
74+ elif s < b or s > 0.25 * (3.0 * a + b ):
7475 con1 = True
7576
7677 con2 = mflag and abs (s - b ) >= 0.5 * abs (b - c )
@@ -116,6 +117,7 @@ def brentq(x1, b, U, gamma, idens, ixmom, iymom, iener,
116117
117118 return x1
118119
120+
119121@njit (cache = True )
120122def cons_to_prim (U ,
121123 irho , iu , iv , ip , ix , irhox ,
@@ -126,53 +128,52 @@ def cons_to_prim(U,
126128 """
127129
128130 qx , qy , _ = U .shape
129-
131+
130132 for j in range (qy ):
131133 for i in range (qx ):
132134 pmax = max ((gamma - 1.0 )* U [i , j , iener ]* 1.0000000001 , smallp )
133135
134- pmin = max (min (1.0e-6 * pmax , smallp ), np .sqrt (U [i , j , ixmom ]** 2 + U [i , j , iymom ]** 2 ) - U [i , j , iener ] - U [i , j , idens ])
136+ pmin = max (min (1.0e-6 * pmax , smallp ), np .sqrt (U [i , j , ixmom ] **
137+ 2 + U [i , j , iymom ]** 2 ) - U [i , j , iener ] - U [i , j , idens ])
135138
136- fmin = f (pmin , U [i ,j , :], gamma , idens , ixmom , iymom , iener )
137- fmax = f (pmax , U [i ,j , :], gamma , idens , ixmom , iymom , iener )
139+ fmin = f (pmin , U [i , j , :], gamma , idens , ixmom , iymom , iener )
140+ fmax = f (pmax , U [i , j , :], gamma , idens , ixmom , iymom , iener )
138141
139142 if fmin * fmax > 0.0 :
140143 pmin = pmin * 1.0e-2
141- fmin = f (pmin , U [i ,j , :], gamma , idens , ixmom , iymom , iener )
144+ fmin = f (pmin , U [i , j , :], gamma , idens , ixmom , iymom , iener )
142145
143146 if fmin * fmax > 0.0 :
144147 pmax = min (pmax * 1.0e2 , 1.0 )
145-
148+
146149 if fmin * fmax > 0.0 :
147150 q [i , j , ip ] = max ((gamma - 1.0 )* U [i , j , iener ], smallp )
148151 else :
149- q [i , j , ip ] = brentq (pmin , pmax , U [i ,j , :], gamma , idens , ixmom , iymom , iener )
150-
152+ q [i , j , ip ] = brentq (pmin , pmax , U [i , j , :], gamma , idens , ixmom , iymom , iener )
153+
151154 if (q [i , j , ip ] != q [i , j , ip ]) or \
152155 (q [i , j , ip ]- 1.0 == q [i , j , ip ]) or \
153- (abs (q [i , j , ip ]) > 1.0e10 ): # nan or infty alert
156+ (abs (q [i , j , ip ]) > 1.0e10 ): # nan or infty alert
154157 q [i , j , ip ] = max ((gamma - 1.0 )* U [i , j , iener ], smallp )
155-
158+
156159 q [i , j , ip ] = max (q [i , j , ip ], smallp )
157160 if abs (U [i , j , iener ] + U [i , j , idens ] + q [i , j , ip ]) < 1.0e-5 :
158161 q [i , j , iu ] = U [i , j , ixmom ]
159162 q [i , j , iv ] = U [i , j , iymom ]
160163 else :
161164 q [i , j , iu ] = U [i , j , ixmom ]/ (U [i , j , iener ] + U [i , j , idens ] + q [i , j , ip ])
162165 q [i , j , iv ] = U [i , j , iymom ]/ (U [i , j , iener ] + U [i , j , idens ] + q [i , j , ip ])
163-
166+
164167 # nan check
165168 if (q [i , j , iu ] != q [i , j , iu ]):
166169 q [i , j , iu ] = 0.0
167-
170+
168171 if (q [i , j , iv ] != q [i , j , iv ]):
169172 q [i , j , iv ] = 0.0
170-
173+
171174 W = 1.0 / np .sqrt (1.0 - q [:, :, iu ]** 2 - q [:, :, iv ]** 2 )
172175
173176 q [:, :, irho ] = U [:, :, idens ] / W
174177 if naux > 0 :
175178 for i in range (naux ):
176179 q [:, :, ix + i ] = U [:, :, irhox + i ]/ (q [:, :, irho ] * W )
177-
178-
0 commit comments