def evolve(self, dt):
        """ 
        Evolve the incompressible equations through one timestep. 
        """
    
        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")

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

        myg = self.cc_data.grid

        dtdx = dt/myg.dx
        dtdy = dt/myg.dy

        #---------------------------------------------------------------------
        # create the limited slopes of u and v (in both directions)
        #---------------------------------------------------------------------
        limiter = self.rp.get_param("incompressible.limiter")
        if (limiter == 0): limitFunc = reconstruction_f.nolimit
        elif (limiter == 1): limitFunc = reconstruction_f.limit2
        else: limitFunc = reconstruction_f.limit4
    
        ldelta_ux = limitFunc(1, u, myg.qx, myg.qy, myg.ng)
        ldelta_vx = limitFunc(1, v, myg.qx, myg.qy, myg.ng)

        ldelta_uy = limitFunc(2, u, myg.qx, myg.qy, myg.ng)
        ldelta_vy = limitFunc(2, v, 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
        print("  making MAC velocities")

        u_MAC, v_MAC = incomp_interface_f.mac_vels(myg.qx, myg.qy, myg.ng, 
                                                   myg.dx, myg.dy, dt,
                                                   u, v,
                                                   ldelta_ux, ldelta_vx,
                                                   ldelta_uy, ldelta_vy,
                                                   gradp_x, gradp_y)


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

        # we will solve L phi = D U^MAC, where phi is cell centered, and
        # U^MAC is the MAC-type staggered grid of the advective
        # velocities.

        print("  MAC projection")

        # create the multigrid object
        mg = MG.CellCenterMG2d(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,
                               verbose=0)

        # first compute divU
        divU = mg.soln_grid.scratch_array()

        # MAC velocities are edge-centered.  divU is cell-centered.
        divU[mg.ilo:mg.ihi+1,mg.jlo:mg.jhi+1] = \
            (u_MAC[myg.ilo+1:myg.ihi+2,myg.jlo:myg.jhi+1] - 
             u_MAC[myg.ilo  :myg.ihi+1,myg.jlo:myg.jhi+1])/myg.dx + \
            (v_MAC[myg.ilo:myg.ihi+1,myg.jlo+1:myg.jhi+2] - 
             v_MAC[myg.ilo:myg.ihi+1,myg.jlo  :myg.jhi+1])/myg.dy
    
        # solve the Poisson problem
        mg.init_zeros()
        mg.init_RHS(divU)
        mg.solve(rtol=1.e-12)

        # update the normal velocities with the pressure gradient -- these
        # constitute our advective velocities
        phi_MAC = self.cc_data.get_var("phi-MAC")
        solution = mg.get_solution()

        phi_MAC[myg.ilo-1:myg.ihi+2,myg.jlo-1:myg.jhi+2] = \
            solution[mg.ilo-1:mg.ihi+2,mg.jlo-1:mg.jhi+2]

        # we need the MAC velocities on all edges of the computational domain
        u_MAC[myg.ilo:myg.ihi+2,myg.jlo:myg.jhi+1] -= \
            (phi_MAC[myg.ilo  :myg.ihi+2,myg.jlo:myg.jhi+1] -
             phi_MAC[myg.ilo-1:myg.ihi+1,myg.jlo:myg.jhi+1])/myg.dx

        v_MAC[myg.ilo:myg.ihi+1,myg.jlo:myg.jhi+2] -= \
            (phi_MAC[myg.ilo:myg.ihi+1,myg.jlo  :myg.jhi+2] -
             phi_MAC[myg.ilo:myg.ihi+1,myg.jlo-1:myg.jhi+1])/myg.dy


        #---------------------------------------------------------------------
        # recompute the interface states, using the advective velocity
        # from above
        #---------------------------------------------------------------------
        print("  making u, v edge states")

        u_xint, v_xint, u_yint, v_yint = \
               incomp_interface_f.states(myg.qx, myg.qy, myg.ng, 
                                         myg.dx, myg.dy, dt,
                                         u, v,
                                         ldelta_ux, ldelta_vx,
                                         ldelta_uy, ldelta_vy,
                                         gradp_x, gradp_y,
                                         u_MAC, v_MAC)


        #---------------------------------------------------------------------
        # update U to get the provisional velocity field
        #---------------------------------------------------------------------

        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[myg.ilo:myg.ihi+1,myg.jlo:myg.jhi+1] = \
            0.5*(u_MAC[myg.ilo  :myg.ihi+1,myg.jlo:myg.jhi+1] + 
                 u_MAC[myg.ilo+1:myg.ihi+2,myg.jlo:myg.jhi+1]) * \
            (u_xint[myg.ilo+1:myg.ihi+2,myg.jlo:myg.jhi+1] - 
             u_xint[myg.ilo  :myg.ihi+1,myg.jlo:myg.jhi+1])/myg.dx + \
            0.5*(v_MAC[myg.ilo:myg.ihi+1,myg.jlo  :myg.jhi+1] + 
                 v_MAC[myg.ilo:myg.ihi+1,myg.jlo+1:myg.jhi%2]) * \
            (u_yint[myg.ilo:myg.ihi+1,myg.jlo+1:myg.jhi+2] - 
             u_yint[myg.ilo:myg.ihi+1,myg.jlo  :myg.jhi+1])/myg.dy 

        advect_y[myg.ilo:myg.ihi+1,myg.jlo:myg.jhi+1] = \
            0.5*(u_MAC[myg.ilo  :myg.ihi+1,myg.jlo:myg.jhi+1] + 
                 u_MAC[myg.ilo+1:myg.ihi+2,myg.jlo:myg.jhi+1]) * \
            (v_xint[myg.ilo+1:myg.ihi+2,myg.jlo:myg.jhi+1] - 
             v_xint[myg.ilo  :myg.ihi+1,myg.jlo:myg.jhi+1])/myg.dx + \
            0.5*(v_MAC[myg.ilo:myg.ihi+1,myg.jlo  :myg.jhi+1] + 
                 v_MAC[myg.ilo:myg.ihi+1,myg.jlo+1:myg.jhi+2]) * \
            (v_yint[myg.ilo:myg.ihi+1,myg.jlo+1:myg.jhi+2] - 
             v_yint[myg.ilo:myg.ihi+1,myg.jlo  :myg.jhi+1])/myg.dy 

             
        proj_type = self.rp.get_param("incompressible.proj_type")

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

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

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


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

        # now we solve L phi = D (U* /dt)
        print("  final projection")
    
        # create the multigrid object
        mg = MG.CellCenterMG2d(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,
                               verbose=0)

        # first compute divU

        # u/v are cell-centered, divU is cell-centered    
        divU[mg.ilo:mg.ihi+1,mg.jlo:mg.jhi+1] = \
            0.5*(u[myg.ilo+1:myg.ihi+2,myg.jlo:myg.jhi+1] - 
                 u[myg.ilo-1:myg.ihi  ,myg.jlo:myg.jhi+1])/myg.dx + \
            0.5*(v[myg.ilo:myg.ihi+1,myg.jlo+1:myg.jhi+2] - 
                 v[myg.ilo:myg.ihi+1,myg.jlo-1:myg.jhi  ])/myg.dy
    
        mg.init_RHS(divU/dt)

        # use the old phi as our initial guess
        phiGuess = mg.soln_grid.scratch_array()
        phiGuess[mg.ilo-1:mg.ihi+2,mg.jlo-1:mg.jhi+2] = \
           phi[myg.ilo-1:myg.ihi+2,myg.jlo-1:myg.jhi+2]
        mg.init_solution(phiGuess)

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

        # store the solution
        solution = mg.get_solution()

        phi[myg.ilo-1:myg.ihi+2,myg.jlo-1:myg.jhi+2] = \
            solution[mg.ilo-1:mg.ihi+2,mg.jlo-1:mg.jhi+2]

        # compute the cell-centered gradient of p and update the velocities
        # this differs depending on what we projected.
        gradphi_x = myg.scratch_array()
        gradphi_y = myg.scratch_array()

        gradphi_x[myg.ilo:myg.ihi+1,myg.jlo:myg.jhi+1] = \
            0.5*(phi[myg.ilo+1:myg.ihi+2,myg.jlo:myg.jhi+1] -
                 phi[myg.ilo-1:myg.ihi  ,myg.jlo:myg.jhi+1])/myg.dx

        gradphi_y[myg.ilo:myg.ihi+1,myg.jlo:myg.jhi+1] = \
            0.5*(phi[myg.ilo:myg.ihi+1,myg.jlo+1:myg.jhi+2] -
                 phi[myg.ilo:myg.ihi+1,myg.jlo-1:myg.jhi  ])/myg.dy

        # u = u - grad_x phi dt
        u[:,:] -= dt*gradphi_x
        v[:,:] -= dt*gradphi_y
        
        # store gradp for the next step
        if (proj_type == 1):
            gradp_x[:,:] += gradphi_x[:,:]
            gradp_y[:,:] += gradphi_y[:,:]

        elif (proj_type == 2):
            gradp_x[:,:] = gradphi_x[:,:]
            gradp_y[:,:] = gradphi_y[:,:]
            
        self.cc_data.fill_BC("x-velocity")
        self.cc_data.fill_BC("y-velocity")
dataset/ETHPy150Open zingale/pyro2/lm_combustion/simulation.py/Simulation.evolve/operator `+` -> `%`