def METHOD_NAME(steps, t, dt):
    rho = nodes1.massDensity()
    rhoMin = rho.min()
    rhoMax = rho.max()
    #print "Rho range : [%16.12e, %16.12e]" % (rhoMin, rhoMax)
    if abs(rhoMin/rhoMax - 1.0) > tol:
        if graphics == "gnu":
            from SpheralGnuPlotUtilities import plotState
            state = State(db, integrator.physicsPackages())
            rhoPlot, velPlot, epsPlot, PPlot, HPlot = plotState(state, plotGhosts=True)
        pos = nodes1.positions()
        for i in range(nodes1.numInternalNodes):
            if rho[i] == rhoMin:
                sys.stderr.write("rho min @ %i %s\n" % (i, pos[i]))
            if rho[i] == rhoMax:
                sys.stderr.write("rho max @ %i %s\n" % (i, pos[i]))
        raise ValueError("rho outside bounds : [%16.12e, %16.12e]" % (rhoMin, rhoMax))
    return