11import mesh .reconstruction as reconstruction
22import numpy as np
33
4+
45def unsplit_fluxes (my_data , rp , dt , scalar_name ):
56 """
67 Construct the fluxes through the interfaces for the linear advection
@@ -58,10 +59,10 @@ def unsplit_fluxes(my_data, rp, dt, scalar_name):
5859
5960 u = my_data .get_var ("x-velocity" )
6061 v = my_data .get_var ("y-velocity" )
61-
62+
6263 cx = myg .scratch_array ()
6364 cy = myg .scratch_array ()
64-
65+
6566 cx .v (buf = 1 )[:, :] = u .v (buf = 1 )* dt / myg .dx
6667 cy .v (buf = 1 )[:, :] = v .v (buf = 1 )* dt / myg .dy
6768
@@ -73,37 +74,37 @@ def unsplit_fluxes(my_data, rp, dt, scalar_name):
7374
7475 ldelta_ax = reconstruction .limit (a , myg , 1 , limiter )
7576 ldelta_ay = reconstruction .limit (a , myg , 2 , limiter )
76-
77+
7778 # x-direction
7879 a_x = myg .scratch_array ()
7980 shift_x = my_data .get_var ("x-shift" ).astype (int )
80-
81- for index ,vel in np .ndenumerate (u .v (buf = 1 )):
81+
82+ for index , vel in np .ndenumerate (u .v (buf = 1 )):
8283 # upwind
8384 if vel < 0 :
8485 a_x .v (buf = 1 )[index ] = a .ip (shift_x .v (buf = 1 )[index ], buf = 1 )[index ] - \
85- 0.5 * (1.0 + cx .v (buf = 1 )[index ])* \
86+ 0.5 * (1.0 + cx .v (buf = 1 )[index ]) * \
8687 ldelta_ax .ip (shift_x .v (buf = 1 )[index ], buf = 1 )[index ]
8788 else :
8889 a_x .v (buf = 1 )[index ] = a .ip (shift_x .v (buf = 1 )[index ], buf = 1 )[index ] + \
89- 0.5 * (1.0 - cx .v (buf = 1 )[index ])* \
90+ 0.5 * (1.0 - cx .v (buf = 1 )[index ]) * \
9091 ldelta_ax .ip (shift_x .v (buf = 1 )[index ], buf = 1 )[index ]
9192
9293 # y-direction
9394 a_y = myg .scratch_array ()
9495 shift_y = my_data .get_var ("y-shift" ).astype (int )
95-
96- for index ,vel in np .ndenumerate (v .v (buf = 1 )):
96+
97+ for index , vel in np .ndenumerate (v .v (buf = 1 )):
9798 # upwind
9899 if vel < 0 :
99100 a_y .v (buf = 1 )[index ] = a .jp (shift_y .v (buf = 1 )[index ], buf = 1 )[index ] - \
100- 0.5 * (1.0 + cy .v (buf = 1 )[index ])* \
101+ 0.5 * (1.0 + cy .v (buf = 1 )[index ]) * \
101102 ldelta_ay .jp (shift_y .v (buf = 1 )[index ], buf = 1 )[index ]
102103 else :
103104 a_y .v (buf = 1 )[index ] = a .jp (shift_y .v (buf = 1 )[index ], buf = 1 )[index ] + \
104- 0.5 * (1.0 - cy .v (buf = 1 )[index ])* \
105+ 0.5 * (1.0 - cy .v (buf = 1 )[index ]) * \
105106 ldelta_ay .jp (shift_y .v (buf = 1 )[index ], buf = 1 )[index ]
106-
107+
107108 # compute the transverse flux differences. The flux is just (u a)
108109 # HOTF
109110 F_xt = u * a_x
@@ -116,15 +117,15 @@ def unsplit_fluxes(my_data, rp, dt, scalar_name):
116117 # depends on the sign of the advective velocity
117118 dtdx2 = 0.5 * dt / myg .dx
118119 dtdy2 = 0.5 * dt / myg .dy
119-
120+
120121 for index , vel in np .ndenumerate (u .v (buf = 1 )):
121122 F_x .v (buf = 1 )[index ] = vel * (a_x .v (buf = 1 )[index ] -
122123 dtdy2 * (F_yt .ip_jp (shift_x .v (buf = 1 )[index ], 1 , buf = 1 )[index ] -
123124 F_yt .ip (shift_x .v (buf = 1 )[index ], buf = 1 )[index ]))
124-
125+
125126 for index , vel in np .ndenumerate (v .v (buf = 1 )):
126127 F_y .v (buf = 1 )[index ] = vel * (a_y .v (buf = 1 )[index ] -
127128 dtdx2 * (F_xt .ip_jp (1 , shift_y .v (buf = 1 )[index ], buf = 1 )[index ] -
128129 F_xt .jp (shift_y .v (buf = 1 )[index ], buf = 1 )[index ]))
130+
129131 return F_x , F_y
130-
0 commit comments