def evolve(self):
        """
        Evolve the low Mach system through one timestep.
        """

        rho = self.cc_data.get_var("density")
        u = self.cc_data.get_var("x-velocity")
        v = self.cc_data.get_var("y-velocity")

        gradp_x = self.cc_data.get_var("gradp_x")
        gradp_y = self.cc_data.get_var("gradp_y")
        
        # note: the base state quantities do not have valid ghost cells
        beta0 = self.base["beta0"]
        beta0_edges = self.base["beta0-edges"]

        rho0 = self.base["rho0"]

        phi = self.cc_data.get_var("phi")

        myg = self.cc_data.grid


        #---------------------------------------------------------------------
        # create the limited slopes of rho, u and v (in both directions)
        #---------------------------------------------------------------------
        limiter = self.rp.get_param("lm-atmosphere.limiter")
        if limiter == 0: limitFunc = reconstruction_f.nolimit
        elif limiter == 1: limitFunc = reconstruction_f.limit2
        else: limitFunc = reconstruction_f.limit4


        ldelta_rx = limitFunc(1, rho.d, myg.qx, myg.qy, myg.ng)
        ldelta_ux = limitFunc(1, u.d, myg.qx, myg.qy, myg.ng)
        ldelta_vx = limitFunc(1, v.d, myg.qx, myg.qy, myg.ng)

        ldelta_ry = limitFunc(2, rho.d, myg.qx, myg.qy, myg.ng)
        ldelta_uy = limitFunc(2, u.d, myg.qx, myg.qy, myg.ng)
        ldelta_vy = limitFunc(2, v.d, myg.qx, myg.qy, myg.ng)

        
        #---------------------------------------------------------------------
        # get the advective velocities
        #---------------------------------------------------------------------

        """
        the advective velocities are the normal velocity through each cell
        interface, and are defined on the cell edges, in a MAC type
        staggered form

                         n+1/2
                        v
                         i,j+1/2
                    +------+------+
                    |             |
            n+1/2   |             |   n+1/2
           u        +     U       +  u
            i-1/2,j |      i,j    |   i+1/2,j
                    |             |
                    +------+------+
                         n+1/2
                        v
                         i,j-1/2

        """

        # this returns u on x-interfaces and v on y-interfaces.  These
        # constitute the MAC grid
        if self.verbose > 0: print("  making MAC velocities")

        # create the coefficient to the grad (pi/beta) term
        coeff = self.aux_data.get_var("coeff")
        coeff.v()[:,:] = 1.0/rho.v()
        coeff.v()[:,:] = coeff.v()*beta0.v2d()
        self.aux_data.fill_BC("coeff")

        # create the source term
        source = self.aux_data.get_var("source_y")

        g = self.rp.get_param("lm-atmosphere.grav")
        rhoprime = self.make_prime(rho, rho0)
        source.v()[:,:] = rhoprime.v()*g/rho.v()
        self.aux_data.fill_BC("source_y")

        _um, _vm = lm_interface_f.mac_vels(myg.qx, myg.qy, myg.ng,
                                           myg.dx, myg.dy, self.dt,
                                           u.d, v.d,
                                           ldelta_ux, ldelta_vx,
                                           ldelta_uy, ldelta_vy,
                                           coeff.d*gradp_x.d, coeff.d*gradp_y.d,
                                           source.d)


        u_MAC = patch.ArrayIndexer(d=_um, grid=myg)
        v_MAC = patch.ArrayIndexer(d=_vm, grid=myg)        


        #---------------------------------------------------------------------
        # do a MAC projection to make the advective velocities divergence
        # free
        #---------------------------------------------------------------------

        # we will solve D (beta_0^2/rho) G phi = D (beta_0 U^MAC), where
        # phi is cell centered, and U^MAC is the MAC-type staggered
        # grid of the advective velocities.

        if self.verbose > 0: print("  MAC projection")

        # create the coefficient array: beta0**2/rho
        # MZ!!!! probably don't need the buf here
        coeff.v(buf=1)[:,:] = 1.0/rho.v(buf=1)
        coeff.v(buf=1)[:,:] = coeff.v(buf=1)*beta0.v2d(buf=1)**2

        # create the multigrid object
        mg = vcMG.VarCoeffCCMG2d(myg.nx, myg.ny,
                                 xl_BC_type=self.cc_data.BCs["phi-MAC"].xlb,
                                 xr_BC_type=self.cc_data.BCs["phi-MAC"].xrb,
                                 yl_BC_type=self.cc_data.BCs["phi-MAC"].ylb,
                                 yr_BC_type=self.cc_data.BCs["phi-MAC"].yrb,
                                 xmin=myg.xmin, xmax=myg.xmax,
                                 ymin=myg.ymin, ymax=myg.ymax,
                                 coeffs=coeff,
                                 coeffs_bc=self.cc_data.BCs["density"],
                                 verbose=0)

        # first compute div{beta_0 U}
        div_beta_U = mg.soln_grid.scratch_array()

        # MAC velocities are edge-centered.  div{beta_0 U} is cell-centered.
        div_beta_U.v()[:,:] = \
            beta0.v2d()*(u_MAC.ip(1) - u_MAC.v())/myg.dx + \
            (beta0_edges.v2dp(1)*v_MAC.jp(1) -
             beta0_edges.v2d()*v_MAC.v())/myg.dy

        # solve the Poisson problem
        mg.init_RHS(div_beta_U.d)
        mg.solve(rtol=1.e-12)


        # update the normal velocities with the pressure gradient -- these
        # constitute our advective velocities.  Note that what we actually
        # solved for here is phi/beta_0
        phi_MAC = self.cc_data.get_var("phi-MAC")
        phi_MAC.d[:,:] = mg.get_solution(grid=myg).d

        coeff = self.aux_data.get_var("coeff")
        coeff.v()[:,:] = 1.0/rho.v()
        coeff.v()[:,:] = coeff.v()*beta0.v2d()
        self.aux_data.fill_BC("coeff")

        coeff_x = myg.scratch_array()
        b = (3, 1, 0, 0)  # this seems more than we need
        coeff_x.v(buf=b)[:,:] = 0.5*(coeff.ip(-1, buf=b) + coeff.v(buf=b))

        coeff_y = myg.scratch_array()
        b = (0, 0, 3, 1)
        coeff_y.v(buf=b)[:,:] = 0.5*(coeff.jp(-1, buf=b) + coeff.v(buf=b))
        
        # we need the MAC velocities on all edges of the computational domain
        # here we do U = U - (beta_0/rho) grad (phi/beta_0)
        b = (0, 1, 0, 0)
        u_MAC.v(buf=b)[:,:] -= \
                coeff_x.v(buf=b)*(phi_MAC.v(buf=b) - phi_MAC.ip(-1, buf=b))/myg.dx

        b = (0, 0, 0, 1)
        v_MAC.v(buf=b)[:,:] -= \
                coeff_y.v(buf=b)*(phi_MAC.v(buf=b) - phi_MAC.jp(-1, buf=b))/myg.dy


        #---------------------------------------------------------------------
        # predict rho to the edges and do its conservative update
        #---------------------------------------------------------------------
        _rx, _ry = lm_interface_f.rho_states(myg.qx, myg.qy, myg.ng,
                                             myg.dx, myg.dy, self.dt,
                                             rho.d, u_MAC.d, v_MAC.d,
                                             ldelta_rx, ldelta_ry)

        rho_xint = patch.ArrayIndexer(d=_rx, grid=myg)
        rho_yint = patch.ArrayIndexer(d=_ry, grid=myg)

        rho_old = rho.copy()

        rho.v()[:,:] -= self.dt*(
            #  (rho u)_x
            (rho_xint.ip(1)*u_MAC.ip(1) - rho_xint.v()*u_MAC.v())/myg.dx +
            #  (rho v)_y
            (rho_yint.jp(1)*v_MAC.jp(1) - rho_yint.v()*v_MAC.v())/myg.dy )

        self.cc_data.fill_BC("density")

        # update eint as a diagnostic
        eint = self.cc_data.get_var("eint")
        gamma = self.rp.get_param("eos.gamma")        
        eint.v()[:,:] = self.base["p0"].v2d()/(gamma - 1.0)/rho.v()
        
        
        #---------------------------------------------------------------------
        # recompute the interface states, using the advective velocity
        # from above
        #---------------------------------------------------------------------
        if self.verbose > 0: print("  making u, v edge states")

        coeff = self.aux_data.get_var("coeff")
        coeff.v()[:,:] = 2.0/(rho.v() + rho_old.v())
        coeff.v()[:,:] = coeff.v()*beta0.v2d()
        self.aux_data.fill_BC("coeff")

        _ux, _vx, _uy, _vy = \
               lm_interface_f.states(myg.qx, myg.qy, myg.ng,
                                     myg.dx, myg.dy, self.dt,
                                     u.d, v.d,
                                     ldelta_ux, ldelta_vx,
                                     ldelta_uy, ldelta_vy,
                                     coeff.d*gradp_x.d, coeff.d*gradp_y.d,
                                     source.d,
                                     u_MAC.d, v_MAC.d)

        u_xint = patch.ArrayIndexer(d=_ux, grid=myg)
        v_xint = patch.ArrayIndexer(d=_vx, grid=myg)
        u_yint = patch.ArrayIndexer(d=_uy, grid=myg)
        v_yint = patch.ArrayIndexer(d=_vy, grid=myg)


        #---------------------------------------------------------------------
        # update U to get the provisional velocity field
        #---------------------------------------------------------------------
        if self.verbose > 0: print("  doing provisional update of u, v")

        # compute (U.grad)U

        # we want u_MAC U_x + v_MAC U_y
        advect_x = myg.scratch_array()
        advect_y = myg.scratch_array()

        advect_x.v()[:,:] = \
            0.5*(u_MAC.v() + u_MAC.ip(1))*(u_xint.ip(1) - u_xint.v())/myg.dx +\
            0.5*(v_MAC.v() + v_MAC.jp(1))*(u_yint.jp(1) - u_yint.v())/myg.dy

        advect_y.v()[:,:] = \
            0.5*(u_MAC.v() + u_MAC.ip(1))*(v_xint.ip(1) - v_xint.v())/myg.dx +\
            0.5*(v_MAC.v() + v_MAC.jp(1))*(v_yint.jp(1) - v_yint.v())/myg.dy


        proj_type = self.rp.get_param("lm-atmosphere.proj_type")

        if proj_type == 1:
            u.v()[:,:] -= (self.dt*advect_x.v() + self.dt*gradp_x.v())
            v.v()[:,:] -= (self.dt*advect_y.v() + self.dt*gradp_y.v())

        elif proj_type == 2:
            u.v()[:,:] -= self.dt*advect_x.v()
            v.v()[:,:] -= self.dt*advect_y.v()


        # add the gravitational source
        rho_half = 0.5*(rho + rho_old)
        rhoprime = self.make_prime(rho_half, rho0)
        source.d[:,:] = (rhoprime*g/rho_half).d
        self.aux_data.fill_BC("source_y")

        v.d[:,:] += self.dt*source.d

        self.cc_data.fill_BC("x-velocity")
        self.cc_data.fill_BC("y-velocity")

        if self.verbose > 0:
            print("min/max rho = {}, {}".format(self.cc_data.min("density"), self.cc_data.max("density")))
            print("min/max u   = {}, {}".format(self.cc_data.min("x-velocity"), self.cc_data.max("x-velocity")))
            print("min/max v   = {}, {}".format(self.cc_data.min("y-velocity"), self.cc_data.max("y-velocity")))


        #---------------------------------------------------------------------
        # project the final velocity
        #---------------------------------------------------------------------

        # now we solve L phi = D (U* /dt)
        if self.verbose > 0: print("  final projection")

        # create the coefficient array: beta0**2/rho
        coeff = 1.0/rho
        coeff.v()[:,:] = coeff.v()*beta0.v2d()**2

        # create the multigrid object
        mg = vcMG.VarCoeffCCMG2d(myg.nx, myg.ny,
                                 xl_BC_type=self.cc_data.BCs["phi"].xlb,
                                 xr_BC_type=self.cc_data.BCs["phi"].xrb,
                                 yl_BC_type=self.cc_data.BCs["phi"].ylb,
                                 yr_BC_type=self.cc_data.BCs["phi"].yrb,
                                 xmin=myg.xmin, xmax=myg.xmax,
                                 ymin=myg.ymin, ymax=myg.ymax,
                                 coeffs=coeff,
                                 coeffs_bc=self.cc_data.BCs["density"],
                                 verbose=0)

        # first compute div{beta_0 U}

        # u/v are cell-centered, divU is cell-centered
        div_beta_U.v()[:,:] = \
            0.5*beta0.v2d()*(u.ip(1) - u.ip(-1))/myg.dx + \
            0.5*(beta0.v2dp(-1)*v.jp(-1) - beta0.v2dp(1)*v.jp(1))/myg.dy

        mg.init_RHS(div_beta_U.d/self.dt)

        # use the old phi as our initial guess
        phiGuess = mg.soln_grid.scratch_array()
        phiGuess.v(buf=1)[:,:] = phi.v(buf=1)
        mg.init_solution(phiGuess.d)

        # solve
        mg.solve(rtol=1.e-12)


        # store the solution in our self.cc_data object -- include a single
        # ghostcell
        phi.d[:,:] = mg.get_solution(grid=myg).d

        # get the cell-centered gradient of p and update the velocities
        # this differs depending on what we projected.
        gradphi_x, gradphi_y = mg.get_solution_gradient(grid=myg)


        # U = U - (beta_0/rho) grad (phi/beta_0)
        coeff = 1.0/rho
        coeff.v()[:,:] = coeff.v()*beta0.v2d()

        u.v()[:,:] -= self.dt*coeff.v()*gradphi_x.v()
        v.v()[:,:] -= self.dt*coeff.v()*gradphi_y.v()

        # store gradp for the next step

        if proj_type == 1:
            gradp_x.v()[:,:] += gradphi_x.v()
            gradp_y.v()[:,:] += gradphi_y.v()

        elif proj_type == 2:
            gradp_x.v()[:,:] = gradphi_x.v()
            gradp_y.v()[:,:] = gradphi_y.v()

        self.cc_data.fill_BC("x-velocity")
        self.cc_data.fill_BC("y-velocity")

        self.cc_data.fill_BC("gradp_x")
        self.cc_data.fill_BC("gradp_y")

        # increment the time
        if not self.in_preevolve:
            self.cc_data.t += self.dt
            self.n += 1
dataset/ETHPy150Open zingale/pyro2/lm_atm/simulation.py/Simulation.evolve/swapped operands of `-`