Skip to content

Commit

Permalink
revert advection files
Browse files Browse the repository at this point in the history
  • Loading branch information
zhichen3 committed Jul 12, 2024
1 parent 829d164 commit 4f2d04e
Show file tree
Hide file tree
Showing 2 changed files with 17 additions and 34 deletions.
28 changes: 11 additions & 17 deletions pyro/advection/interface.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,9 @@ def linear_interface(a, myg, rp, dt):
u = rp.get_param("advection.u")
v = rp.get_param("advection.v")

cx = u*dt/myg.dx
cy = v*dt/myg.dy

# --------------------------------------------------------------------------
# monotonized central differences
# --------------------------------------------------------------------------
Expand All @@ -16,34 +19,25 @@ def linear_interface(a, myg, rp, dt):
ldelta_ax = reconstruction.limit(a, myg, 1, limiter)
ldelta_ay = reconstruction.limit(a, myg, 2, limiter)

ldelta_Fx = reconstruction.limit(u*a, myg, 1, limiter)
ldelta_Fy = reconstruction.limit(v*a, myg, 2, limiter)

# x-direction
a_x = myg.scratch_array()

# upwind
if u < 0:
# a_x[i,j] = a[i,j] - 0.5*ldelta_ax[i,j] - 0.5*dt/dx*ldelta_Fx[i,j]
a_x.v(buf=1)[:, :] = a.v(buf=1) - 0.5*ldelta_ax.v(buf=1) \
- 0.5*dt/myg.dx*ldelta_Fx.v(buf=1)

# a_x[i,j] = a[i,j] - 0.5*(1.0 + cx)*ldelta_a[i,j]
a_x.v(buf=1)[:, :] = a.v(buf=1) - 0.5*(1.0 + cx)*ldelta_ax.v(buf=1)
else:
# a_x[i,j] = a[i-1,j] + 0.5*ldelta_ax[i-1,j] - 0.5*dt/dx*ldelta_Fx[i-1,j]
a_x.v(buf=1)[:, :] = a.ip(-1, buf=1) + 0.5*ldelta_ax.ip(-1, buf=1) \
- 0.5*dt/myg.dx*ldelta_Fx.ip(-1, buf=1)
# a_x[i,j] = a[i-1,j] + 0.5*(1.0 - cx)*ldelta_a[i-1,j]
a_x.v(buf=1)[:, :] = a.ip(-1, buf=1) + 0.5*(1.0 - cx)*ldelta_ax.ip(-1, buf=1)

# y-direction
a_y = myg.scratch_array()

# upwind
if v < 0:
# a_y[i,j] = a[i,j] - 0.5*ldelta_ay[i,j] - 0.5*dt/dy*ldelta_Fy[i,j]
a_y.v(buf=1)[:, :] = a.v(buf=1) - 0.5*ldelta_ay.v(buf=1) \
- 0.5*dt/myg.dy*ldelta_Fy.v(buf=1)
# a_y[i,j] = a[i,j] - 0.5*(1.0 + cy)*ldelta_a[i,j]
a_y.v(buf=1)[:, :] = a.v(buf=1) - 0.5*(1.0 + cy)*ldelta_ay.v(buf=1)
else:
# a_y[i,j] = a[i,j-1] + 0.5*ldelta_ay[i,j-1] - 0.5*dt/dy*ldelta_Fy[i,j-1]
a_y.v(buf=1)[:, :] = a.jp(-1, buf=1) + 0.5*ldelta_ay.jp(-1, buf=1) \
- 0.5*dt/myg.dy*ldelta_Fy.jp(-1, buf=1)
# a_y[i,j] = a[i,j-1] + 0.5*(1.0 - cy)*ldelta_a[i,j-1]
a_y.v(buf=1)[:, :] = a.jp(-1, buf=1) + 0.5*(1.0 - cy)*ldelta_ay.jp(-1, buf=1)

return u, v, a_x, a_y
23 changes: 6 additions & 17 deletions pyro/advection/simulation.py
Original file line number Diff line number Diff line change
Expand Up @@ -52,14 +52,8 @@ def method_compute_timestep(self):
v = self.rp.get_param("advection.v")

# the timestep is min(dx/|u|, dy/|v|)
dx_min = np.min(self.cc_data.grid.V() / \
self.cc_data.grid.A_x())

dy_min = np.min(self.cc_data.grid.V() / \
self.cc_data.grid.A_y())

xtmp = dx_min/max(abs(u), self.SMALL)
ytmp = dy_min/max(abs(v), self.SMALL)
xtmp = self.cc_data.grid.dx/max(abs(u), self.SMALL)
ytmp = self.cc_data.grid.dy/max(abs(v), self.SMALL)

self.dt = cfl*min(xtmp, ytmp)

Expand All @@ -70,7 +64,8 @@ def evolve(self):
is part of the Simulation.
"""

myg = self.cc_data.grid
dtdx = self.dt/self.cc_data.grid.dx
dtdy = self.dt/self.cc_data.grid.dy

flux_x, flux_y = flx.unsplit_fluxes(self.cc_data, self.rp, self.dt, "density", linear_interface)

Expand All @@ -85,14 +80,8 @@ def evolve(self):

dens = self.cc_data.get_var("density")

dens.v()[:, :] = dens.v() + self.dt / myg.V() * \
(myg.A_x()*flux_x.v() - myg.A_x(1)*flux_x.ip(1) + \
myg.A_y()*flux_y.v() - myg.A_y(1)*flux_y.jp(1))


# dens.v()[:, :] = dens.v() \
# + dtdx*(flux_x.v() - flux_x.ip(1)) + \
# dtdy*(flux_y.v() - flux_y.jp(1))
dens.v()[:, :] = dens.v() + dtdx*(flux_x.v() - flux_x.ip(1)) + \
dtdy*(flux_y.v() - flux_y.jp(1))

if self.particles is not None:
myg = self.cc_data.grid
Expand Down

0 comments on commit 4f2d04e

Please sign in to comment.