diff --git a/pyro/incompressible/incomp_interface.py b/pyro/incompressible/incomp_interface.py index 3eaad0129..751c72172 100644 --- a/pyro/incompressible/incomp_interface.py +++ b/pyro/incompressible/incomp_interface.py @@ -32,12 +32,29 @@ def mac_vels(grid, dt, # get the full u and v left and right states (including transverse # terms) on both the x- and y-interfaces - # pylint: disable-next=unused-variable - u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr = get_interface_states(grid, dt, - u, v, - ldelta_ux, ldelta_vx, - ldelta_uy, ldelta_vy, - gradp_x, gradp_y) + + u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr = \ + burgers_interface.get_interface_states(grid, dt, + u, v, + ldelta_ux, ldelta_vx, + ldelta_uy, ldelta_vy) + + # apply transverse terms on both x- and y- interfaces + u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr = \ + burgers_interface.apply_transverse_corrections(grid, dt, + u_xl, u_xr, + u_yl, u_yr, + v_xl, v_xr, + v_yl, v_yr) + + # apply pressure gradient correction terms to the interface state + u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr = \ + apply_gradp_corrections(dt, + u_xl, u_xr, + u_yl, u_yr, + v_xl, v_xr, + v_yl, v_yr, + gradp_x, gradp_y) # Riemann problem -- this follows Burger's equation. We don't use # any input velocity for the upwinding. Also, we only care about @@ -82,13 +99,29 @@ def states(grid, dt, x and y velocities predicted to the interfaces """ - # get the full u and v left and right states (including transverse - # terms) on both the x- and y-interfaces - u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr = get_interface_states(grid, dt, - u, v, - ldelta_ux, ldelta_vx, - ldelta_uy, ldelta_vy, - gradp_x, gradp_y) + # get the full u and v left and right states without transverse terms and gradp + u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr = \ + burgers_interface.get_interface_states(grid, dt, + u, v, + ldelta_ux, ldelta_vx, + ldelta_uy, ldelta_vy) + + # apply transverse terms on both x- and y- interfaces + u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr = \ + burgers_interface.apply_transverse_corrections(grid, dt, + u_xl, u_xr, + u_yl, u_yr, + v_xl, v_xr, + v_yl, v_yr) + + # apply pressure gradient correction terms to the interface state + u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr = \ + apply_gradp_corrections(dt, + u_xl, u_xr, + u_yl, u_yr, + v_xl, v_xr, + v_yl, v_yr, + gradp_x, gradp_y) # upwind using the MAC velocity to determine which state exists on # the interface @@ -100,128 +133,51 @@ def states(grid, dt, return u_xint, v_xint, u_yint, v_yint -def get_interface_states(grid, dt, - u, v, - ldelta_ux, ldelta_vx, - ldelta_uy, ldelta_vy, - gradp_x, gradp_y): +def apply_gradp_corrections(dt, + u_xl, u_xr, + u_yl, u_yr, + v_xl, v_xr, + v_yl, v_yr, + gradp_x, gradp_y): r""" - Compute the unsplit predictions of u and v on both the x- and - y-interfaces. This includes the transverse terms. - Parameters ---------- grid : Grid2d The grid object dt : float The timestep - u, v : ndarray - x-velocity and y-velocity - ldelta_ux, ldelta_uy: ndarray - Limited slopes of the x-velocity in the x and y directions - ldelta_vx, ldelta_vy: ndarray - Limited slopes of the y-velocity in the x and y directions - gradp_x, gradp_y : ndarray - Pressure gradients in the x and y directions - + u_xl, u_xr : ndarray ndarray + left and right states of x-velocity in x-interface. + u_yl, u_yr : ndarray ndarray + left and right states of x-velocity in y-interface. + v_xl, v_xr : ndarray ndarray + left and right states of y-velocity in x-interface. + v_yl, u_yr : ndarray ndarray + left and right states of y-velocity in y-interface. Returns ------- out : ndarray, ndarray, ndarray, ndarray, ndarray, ndarray, ndarray, ndarray - unsplit predictions of u and v on both the x- and - y-interfaces + correct the interface states of the left and right states of u and v on + both the x- and y-interfaces interface states with the pressure gradient + terms. """ - u_xl = grid.scratch_array() - u_xr = grid.scratch_array() - u_yl = grid.scratch_array() - u_yr = grid.scratch_array() - - v_xl = grid.scratch_array() - v_xr = grid.scratch_array() - v_yl = grid.scratch_array() - v_yr = grid.scratch_array() - - # first predict u and v to both interfaces, considering only the normal - # part of the predictor. These are the 'hat' states. - - dtdx = dt / grid.dx - dtdy = dt / grid.dy - - # u on x-edges - u_xl.ip(1, buf=2)[:, :] = u.v(buf=2) + \ - 0.5 * (1.0 - dtdx * u.v(buf=2)) * ldelta_ux.v(buf=2) - u_xr.v(buf=2)[:, :] = u.v(buf=2) - \ - 0.5 * (1.0 + dtdx * u.v(buf=2)) * ldelta_ux.v(buf=2) - - # v on x-edges - v_xl.ip(1, buf=2)[:, :] = v.v(buf=2) + \ - 0.5 * (1.0 - dtdx * u.v(buf=2)) * ldelta_vx.v(buf=2) - v_xr.v(buf=2)[:, :] = v.v(buf=2) - \ - 0.5 * (1.0 + dtdx * u.v(buf=2)) * ldelta_vx.v(buf=2) - - # u on y-edges - u_yl.jp(1, buf=2)[:, :] = u.v(buf=2) + \ - 0.5 * (1.0 - dtdy * v.v(buf=2)) * ldelta_uy.v(buf=2) - u_yr.v(buf=2)[:, :] = u.v(buf=2) - \ - 0.5 * (1.0 + dtdy * v.v(buf=2)) * ldelta_uy.v(buf=2) - - # v on y-edges - v_yl.jp(1, buf=2)[:, :] = v.v(buf=2) + \ - 0.5 * (1.0 - dtdy * v.v(buf=2)) * ldelta_vy.v(buf=2) - v_yr.v(buf=2)[:, :] = v.v(buf=2) - \ - 0.5 * (1.0 + dtdy * v.v(buf=2)) * ldelta_vy.v(buf=2) - - # now get the normal advective velocities on the interfaces by solving - # the Riemann problem. - uhat_adv = burgers_interface.riemann(grid, u_xl, u_xr) - vhat_adv = burgers_interface.riemann(grid, v_yl, v_yr) - - # now that we have the advective velocities, upwind the left and right - # states using the appropriate advective velocity. - - # on the x-interfaces, we upwind based on uhat_adv - u_xint = burgers_interface.upwind(grid, u_xl, u_xr, uhat_adv) - v_xint = burgers_interface.upwind(grid, v_xl, v_xr, uhat_adv) - - # on the y-interfaces, we upwind based on vhat_adv - u_yint = burgers_interface.upwind(grid, u_yl, u_yr, vhat_adv) - v_yint = burgers_interface.upwind(grid, v_yl, v_yr, vhat_adv) - - # at this point, these states are the `hat' states -- they only - # considered the normal to the interface portion of the predictor. - - ubar = grid.scratch_array() - vbar = grid.scratch_array() - - ubar.v(buf=2)[:, :] = 0.5 * (uhat_adv.v(buf=2) + uhat_adv.ip(1, buf=2)) - vbar.v(buf=2)[:, :] = 0.5 * (vhat_adv.v(buf=2) + vhat_adv.jp(1, buf=2)) - - # v du/dy is the transerse term for the u states on x-interfaces - vu_y = grid.scratch_array() - vu_y.v(buf=2)[:, :] = vbar.v(buf=2) * (u_yint.jp(1, buf=2) - u_yint.v(buf=2)) - - u_xl.ip(1, buf=2)[:, :] += -0.5 * dtdy * vu_y.v(buf=2) - 0.5 * dt * gradp_x.v(buf=2) - u_xr.v(buf=2)[:, :] += -0.5 * dtdy * vu_y.v(buf=2) - 0.5 * dt * gradp_x.v(buf=2) - - # v dv/dy is the transverse term for the v states on x-interfaces - vv_y = grid.scratch_array() - vv_y.v(buf=2)[:, :] = vbar.v(buf=2) * (v_yint.jp(1, buf=2) - v_yint.v(buf=2)) - - v_xl.ip(1, buf=2)[:, :] += -0.5 * dtdy * vv_y.v(buf=2) - 0.5 * dt * gradp_y.v(buf=2) - v_xr.v(buf=2)[:, :] += -0.5 * dtdy * vv_y.v(buf=2) - 0.5 * dt * gradp_y.v(buf=2) + # Apply pressure gradient correction terms - # u dv/dx is the transverse term for the v states on y-interfaces - uv_x = grid.scratch_array() - uv_x.v(buf=2)[:, :] = ubar.v(buf=2) * (v_xint.ip(1, buf=2) - v_xint.v(buf=2)) + # transverse term for the u states on x-interfaces + u_xl.ip(1, buf=2)[:, :] += - 0.5 * dt * gradp_x.v(buf=2) + u_xr.v(buf=2)[:, :] += - 0.5 * dt * gradp_x.v(buf=2) - v_yl.jp(1, buf=2)[:, :] += -0.5 * dtdx * uv_x.v(buf=2) - 0.5 * dt * gradp_y.v(buf=2) - v_yr.v(buf=2)[:, :] += -0.5 * dtdx * uv_x.v(buf=2) - 0.5 * dt * gradp_y.v(buf=2) + # transverse term for the v states on x-interfaces + v_xl.ip(1, buf=2)[:, :] += - 0.5 * dt * gradp_y.v(buf=2) + v_xr.v(buf=2)[:, :] += - 0.5 * dt * gradp_y.v(buf=2) - # u du/dx is the transverse term for the u states on y-interfaces - uu_x = grid.scratch_array() - uu_x.v(buf=2)[:, :] = ubar.v(buf=2) * (u_xint.ip(1, buf=2) - u_xint.v(buf=2)) + # transverse term for the v states on y-interfaces + v_yl.jp(1, buf=2)[:, :] += - 0.5 * dt * gradp_y.v(buf=2) + v_yr.v(buf=2)[:, :] += - 0.5 * dt * gradp_y.v(buf=2) - u_yl.jp(1, buf=2)[:, :] += -0.5 * dtdx * uu_x.v(buf=2) - 0.5 * dt * gradp_x.v(buf=2) - u_yr.v(buf=2)[:, :] += -0.5 * dtdx * uu_x.v(buf=2) - 0.5 * dt * gradp_x.v(buf=2) + # transverse term for the u states on y-interfaces + u_yl.jp(1, buf=2)[:, :] += - 0.5 * dt * gradp_x.v(buf=2) + u_yr.v(buf=2)[:, :] += - 0.5 * dt * gradp_x.v(buf=2) return u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr diff --git a/pyro/incompressible/simulation.py b/pyro/incompressible/simulation.py index fb3e4cc83..0b9b07d0f 100644 --- a/pyro/incompressible/simulation.py +++ b/pyro/incompressible/simulation.py @@ -4,14 +4,15 @@ import numpy as np import pyro.mesh.array_indexer as ai +from pyro.burgers import Simulation as burgers_simulation from pyro.incompressible import incomp_interface from pyro.mesh import patch, reconstruction from pyro.multigrid import MG from pyro.particles import particles -from pyro.simulation_null import NullSimulation, bc_setup, grid_setup +from pyro.simulation_null import bc_setup, grid_setup -class Simulation(NullSimulation): +class Simulation(burgers_simulation): def initialize(self): """ @@ -51,27 +52,6 @@ def initialize(self): problem = importlib.import_module(f"pyro.incompressible.problems.{self.problem_name}") problem.init_data(self.cc_data, self.rp) - def method_compute_timestep(self): - """ - The timestep() function computes the advective timestep - (CFL) constraint. The CFL constraint says that information - cannot propagate further than one zone per timestep. - - We use the driver.cfl parameter to control what fraction of the CFL - step we actually take. - """ - - cfl = self.rp.get_param("driver.cfl") - - u = self.cc_data.get_var("x-velocity") - v = self.cc_data.get_var("y-velocity") - - # the timestep is min(dx/|u|, dy|v|) - xtmp = self.cc_data.grid.dx/(abs(u)) - ytmp = self.cc_data.grid.dy/(abs(v)) - - self.dt = cfl*float(min(xtmp.min(), ytmp.min())) - def preevolve(self): """ preevolve is called before we being the timestepping loop. For