def cell_center_data_clone(old):
    """
    Create a new CellCenterData2d object that is a copy of an existing
    one

    Parameters
    ----------
    old : CellCenterData2d object
        The CellCenterData2d object we wish to copy

    Note
    ----
    It may be that this whole thing can be replaced with a copy.deepcopy()

    """

    if not isinstance(old, CellCenterData2d):
        msg.fail("Can't clone object")

    new = CellCenterData2d(old.grid, dtype=old.dtype)

    for n in range(old.nvar):
        new.register_var(old.vars[n], old.BCs[old.vars[n]])

    new.create()

    new.aux = old.aux.copy()
    new.data = old.data.copy()

    return new

    def pretty_print(self, varname):
        """
        Print out a small dataset to the screen with the ghost cells
        a different color, to make things stand out
        """

        a = self.get_var(varname)

        if self.dtype == numpy.int:
            fmt = "%4d"
        elif self.dtype == numpy.float64:
            fmt = "%10.5g"
        else:
            msg.fail("ERROR: dtype not supported")
        
        j = 0
        while j < self.grid.qy:
            i = 0
            while i < self.grid.qx:

                if (j < self.grid.jlo or j > self.grid.jhi or
                    i < self.grid.ilo or i > self.grid.ihi):
                    gc = 1
                else:
                    gc = 0

                if gc:
                    print "\033[31m" + fmt % (a[i,j]) + "\033[0m" ,
                else:
                    print fmt % (a[i,j]) ,

                i += 1

            print " "
            j += 1

    def __init__(self, solver_name):
        """
        Constructor

        Parameters
        ----------
        solver_name : str
            Name of solver to use
        """

        msg.bold('pyro ...')

        if solver_name not in valid_solvers:
            msg.fail("ERROR: %s is not a valid solver" % solver_name)

        self.pyro_home = os.path.dirname(os.path.realpath(__file__)) + '/'

        # import desired solver under "solver" namespace
        self.solver = importlib.import_module(solver_name)
        self.solver_name = solver_name

        # -------------------------------------------------------------------------
        # runtime parameters
        # -------------------------------------------------------------------------

        # parameter defaults
        self.rp = runparams.RuntimeParameters()
        self.rp.load_params(self.pyro_home + "_defaults")
        self.rp.load_params(self.pyro_home + solver_name + "/_defaults")

        self.tc = profile.TimerCollection()

        self.is_initialized = False

    def print_paramfile(self):
        """
        Create a file, inputs.auto, that has the structure of a pyro
        inputs file, with all known parameters and values
        """

        all_keys = list(self.params.keys())

        try:
            f = open('inputs.auto', 'w')
        except IOError:
            msg.fail("ERROR: unable to open inputs.auto")

        f.write('# automagically generated parameter file\n')

        # find all the sections
        secs = set([q for (q, _) in [k.split(".") for k in all_keys]])

        for sec in sorted(secs):
            keys = [q for q in all_keys if q.startswith("{}.".format(sec))]

            f.write("\n[{}]\n".format(sec))

            for key in keys:
                _, option = key.split('.')

                value = self.params[key]

                if self.param_comments[key] != '':
                    f.write("{} = {}    ; {}\n".format(option, value, self.param_comments[key]))
                else:
                    f.write("{} = {}\n".format(option, value))

        f.close()

    def load_params(self, pfile, no_new=0):
        """
        Reads line from file and makes dictionary pairs from the data
        to store.

        Parameters
        ----------
        file : str
            The name of the file to parse
        no_new : int, optional
            If no_new = 1, then we don't add any new paramters to the
            dictionary of runtime parameters, but instead just override
            the values of existing ones.

        """

        # check to see whether the file exists
        try: f = open(pfile, 'r')
        except IOError:
            msg.fail("ERROR: parameter file does not exist: {}".format(pfile))

        # we could use the ConfigParser, but we actually want to
        # have our configuration files be self-documenting, of the
        # format key = value ; comment
        sec = re.compile(r'^\[(.*)\]')
        eq = re.compile(r'^([^=#]+)=([^;]+);{0,1}(.*)')

        for line in f.readlines():

            if sec.search(line):
                lbracket, section, rbracket = sec.split(line)
                section = section.strip().lower()

            elif eq.search(line):
                left, item, value, comment, right = eq.split(line)
                item = item.strip().lower()

                # define the key
                key = section + "." + item

                # if we have no_new = 1, then we only want to override existing
                # key/values
                if no_new:
                    if not key in self.params.keys():
                        msg.warning("warning, key: %s not defined" % (key))
                        continue

                self.params[key] = _get_val(value)

                # if the comment already exists (i.e. from reading in
                # _defaults) and we are just resetting the value of
                # the parameter (i.e.  from reading in inputs), then
                # we don't want to destroy the comment
                if comment.strip() == "":
                    try:
                        comment = self.param_comments[key]
                    except KeyError:
                        comment = ""

                self.param_comments[key] = comment.strip()

def LoadParams(file, noNew=0):
    """
    reads lines from file and makes dictionary pairs from the data
    to store in globalParams.
    """
    global globalParams

    # check to see whether the file exists
    try: f = open(file, 'r')
    except IOError:
        msg.fail("ERROR: parameter file does not exist: %s" % (file))


    # we could use the ConfigParser, but we actually want to have
    # our configuration files be self-documenting, of the format
    # key = value     ; comment
    sec = re.compile(r'^\[(.*)\]')
    eq = re.compile(r'^([^=#]+)=([^;]+);{0,1}(.*)')

    for line in f.readlines():

        if sec.search(line): 
            lbracket, section, rbracket = sec.split(line)
            section = string.lower(section.strip())
            
        elif eq.search(line):
            left, item, value, comment, right = eq.split(line) 		
            item = string.lower(item.strip())

            # define the key
            key = section + "." + item
            
            # if we have noNew = 1, then we only want to override existing
            # key/values
            if (noNew):
                if (not key in globalParams.keys()):
                    msg.warning("warning, key: %s not defined" % (key))
                    continue

            # check in turn whether this is an interger, float, or string
            if (isInt(value)):
                globalParams[key] = int(value)
            elif (isFloat(value)):
                globalParams[key] = float(value)
            else:
                globalParams[key] = value.strip()

            # if the comment already exists (i.e. from reading in _defaults)
            # and we are just resetting the value of the parameter (i.e.
            # from reading in inputs), then we don't want to destroy the
            # comment
            if comment.strip() == "":
                try:
                    comment = globalParamComments[key]
                except KeyError:
                    comment = ""
                    
            globalParamComments[key] = comment.strip()

    def get_var(self, v):
        """
        Alias for cc_data's get_var routine, returns the cell-centered data
        given the variable name v.
        """

        if not self.is_initialized:
            msg.fail("ERROR: problem has not been initialized")

        return self.sim.cc_data.get_var(v)

def initData(my_data):
    """ initialize the incompressible shear problem """

    msg.bold("initializing the incompressible shear problem...")

    rp = my_data.rp

    # make sure that we are passed a valid patch object
    if not isinstance(my_data, patch.CellCenterData2d):
        print my_data.__class__
        msg.fail("ERROR: patch invalid in shear.py")


    # get the necessary runtime parameters
    rho_s = rp.get_param("shear.rho_s")
    delta_s = rp.get_param("shear.delta_s")

    
    # get the velocities
    u = my_data.get_var("x-velocity")
    v = my_data.get_var("y-velocity")

    myg = my_data.grid

    if (myg.xmin != 0 or myg.xmax != 1 or
        myg.ymin != 0 or myg.ymax != 1):
        msg.fail("ERROR: domain should be a unit square")
        
    y_half = 0.5*(myg.ymin + myg.ymax)

    print 'y_half = ', y_half
    print 'delta_s = ', delta_s
    print 'rho_s = ', rho_s
    
    # there is probably an easier way to do this without loops, but
    # for now, we will just do an explicit loop.
    i = myg.ilo
    while i <= myg.ihi:

        j = myg.jlo
        while j <= myg.jhi:

            if (myg.y[j] <= y_half):
                u[i,j] = numpy.tanh(rho_s*(myg.y[j] - 0.25))
            else:
                u[i,j] = numpy.tanh(rho_s*(0.75 - myg.y[j]))
            
            v[i,j] = delta_s*numpy.sin(2.0*math.pi*myg.x[i])
            
            j += 1
        i += 1
        
    
    print "extrema: ", numpy.min(u.flat), numpy.max(u.flat)

    def create(self):
        """
        Called after all the variables are registered and allocates
        the storage for the state data.
        """

        if self.initialized == 1:
            msg.fail("ERROR: grid already initialized")

        self.data = np.zeros((self.nvar, self.grid.qx, self.grid.qy),
                                dtype=self.dtype)
        self.initialized = 1

    def initialize(self):
        """ 
        Initialize the grid and variables for diffusion and set the initial
        conditions for the chosen problem.
        """

        # setup the grid
        nx = self.rp.get_param("mesh.nx")
        ny = self.rp.get_param("mesh.ny")
        
        xmin = self.rp.get_param("mesh.xmin")
        xmax = self.rp.get_param("mesh.xmax")
        ymin = self.rp.get_param("mesh.ymin")
        ymax = self.rp.get_param("mesh.ymax")
    
        my_grid = patch.Grid2d(nx, ny, 
                               xmin=xmin, xmax=xmax, 
                               ymin=ymin, ymax=ymax, ng=1)


        # create the variables

        # first figure out the boundary conditions -- we allow periodic,
        # Dirichlet, and Neumann.

        xlb_type = self.rp.get_param("mesh.xlboundary")
        xrb_type = self.rp.get_param("mesh.xrboundary")
        ylb_type = self.rp.get_param("mesh.ylboundary")
        yrb_type = self.rp.get_param("mesh.yrboundary")

        bcparam = []
        for bc in [xlb_type, xrb_type, ylb_type, yrb_type]:
            if bc == "periodic": bcparam.append("periodic")
            elif bc == "neumann":  bcparam.append("neumann")
            elif bc == "dirichlet":  bcparam.append("dirichlet")
            else:
                msg.fail("invalid BC")


        bc = patch.BCObject(xlb=bcparam[0], xrb=bcparam[1], 
                            ylb=bcparam[2], yrb=bcparam[3])    


        my_data = patch.CellCenterData2d(my_grid)

        my_data.register_var("phi", bc)

        my_data.create()

        self.cc_data = my_data

        # now set the initial conditions for the problem           
        exec(self.problem_name + '.init_data(self.cc_data, self.rp)')

    def __init__(self, sim_data, bc, n_particles, particle_generator="grid",
                 pos_array=None, init_array=None):
        """
        Initialize the Particles object.

        Particles are stored as a dictionary, with their keys being tuples
        of their initial position. This was done in order to have a simple way
        to access the initial particle positions when plotting.

        However, this assumes that no two particles are
        initialised with the same initial position, which is fine for the
        massless particle case, however could no longer be a sensible thing
        to do if have particles have other properties (e.g. mass).

        Parameters
        ----------
        sim_data : CellCenterData2d object
            The cell-centered simulation data
        bc : BC object
            Boundary conditions
        n_particles : int
            Number of particles
        particle_generator : string or function
            String with generator name of custom particle generator function
        pos_array : float array
            Array of particle positions to use with particle initialization
        init_array : float array
            Array of initial particle positions required for plotting from file.
        """

        self.sim_data = sim_data
        self.bc = bc
        self.particles = dict()

        if n_particles <= 0:
            msg.fail("ERROR: n_particles = %s <= 0" % (n_particles))

        if callable(particle_generator):  # custom particle generator function
            self.particles = particle_generator(n_particles)
        else:
            if particle_generator == "random":
                self.randomly_generate_particles(n_particles)
            elif particle_generator == "grid":
                self.grid_generate_particles(n_particles)
            elif particle_generator == "array":
                self.array_generate_particles(pos_array, init_array)
            else:
                msg.fail("ERROR: do not recognise particle generator %s"
                         % (particle_generator))

        self.n_particles = len(self.particles)

    def registerVar(self, name, bcObject):
        """ 
        register a variable with ccData2d object.  Here we pass in a
        bcObject that describes the boundary conditions for that
        variable.
        """

        if (self.initialized == 1):
            msg.fail("ERROR: grid already initialized")

        self.vars.append(name)
        self.nvar += 1

        self.BCs[name] = bcObject

    def create(self):
        """
        called after all the variables are registered and allocates
        the storage for the state data
        """

        if (self.initialized) == 1:
            msg.fail("ERROR: grid already initialized")

        self.data = numpy.zeros((self.nvar,
                                 2*self.grid.ng+self.grid.nx, 
                                 2*self.grid.ng+self.grid.ny),
                                dtype=self.dtype)
        self.initialized = 1

    def store_as_benchmark(self):
        """ Are we storing a benchmark? """

        if not os.path.isdir(self.solver_name + "/tests/"):
            try:
                os.mkdir(self.solver_name + "/tests/")
            except (FileNotFoundError, PermissionError):
                msg.fail(
                    "ERROR: unable to create the solver's tests/ directory")

        basename = self.rp.get_param("io.basename")
        bench_file = self.pyro_home + self.solver_name + "/tests/" + \
            basename + "%4.4d" % (self.sim.n)
        msg.warning("storing new benchmark: {}\n".format(bench_file))
        self.sim.write(bench_file)

def init_data(my_data, rp):
    """ initialize the incompressible shear problem """

    msg.bold("initializing the incompressible shear problem...")

    # make sure that we are passed a valid patch object
    if not isinstance(my_data, patch.CellCenterData2d):
        print(my_data.__class__)
        msg.fail("ERROR: patch invalid in shear.py")

    # get the necessary runtime parameters
    eps = rp.get_param("vortex.eps")

    print('eps = ', eps)

    # get the velocities
    u = my_data.get_var("x-velocity")
    v = my_data.get_var("y-velocity")

    myg = my_data.grid

    u.d[:,:] = -np.sin(math.pi*myg.y2d)
    v.d[:,:] = np.sin(math.pi*myg.x2d)
    #u.d[:,:] = -np.sin(2.0*math.pi*myg.x2d)*np.cos(2.0*math.pi*myg.y2d)*ran
    #v.d[:,:] = np.cos(2.0*math.pi*myg.x2d)*np.sin(2.0*math.pi*myg.y2d)*ran

    if eps != 0.0:
    #perturbed velocity1 at (0,0)
      r2 = myg.x2d**2+myg.y2d**2
      dvx1l = -eps**3*myg.y2d/r2*(1-np.exp(-r2/eps**2))
      dvy1l = eps**3*myg.x2d/r2*(1-np.exp(-r2/eps**2))

    #perturbed velocity1 at (2pi,0)
      r2 = (myg.x2d - 2.0)**2+myg.y2d**2
      dvx1r = -eps**3*myg.y2d/r2*(1-np.exp(-r2/eps**2))
      dvy1r = eps**3*(myg.x2d-2.0)/r2*(1-np.exp(-r2/eps**2))


    #perturbed velocity2 at (pi,0)
      r2 = (myg.x2d - 1.0)**2+myg.y2d**2
      dvx2 = eps**3*myg.y2d/r2*(1-np.exp(-r2/eps**2))
      dvy2 = -eps**3*(myg.x2d-1.0)/r2*(1-np.exp(-r2/eps**2))

      u.d[:,:] = u.d[:,:] + dvx1l + dvx1r + dvx2
      v.d[:,:] = v.d[:,:] + dvy1l + dvy1r + dvy2

    print("extrema: ", u.min(), u.max())

def initialize(rp):
    """ 
    initialize the grid and variables for diffusion
    """

    # setup the grid
    nx = rp.get_param("mesh.nx")
    ny = rp.get_param("mesh.ny")

    xmin = rp.get_param("mesh.xmin")
    xmax = rp.get_param("mesh.xmax")
    ymin = rp.get_param("mesh.ymin")
    ymax = rp.get_param("mesh.ymax")
    
    my_grid = patch.Grid2d(nx, ny, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, ng=1)


    # create the variables

    # first figure out the boundary conditions -- we allow periodic,
    # Dirichlet, and Neumann.

    xlb_type = rp.get_param("mesh.xlboundary")
    xrb_type = rp.get_param("mesh.xrboundary")
    ylb_type = rp.get_param("mesh.ylboundary")
    yrb_type = rp.get_param("mesh.yrboundary")

    bcparam = []
    for bc in [xlb_type, xrb_type, ylb_type, yrb_type]:
        if   (bc == "periodic"): bcparam.append("periodic")
        elif (bc == "neumann"):  bcparam.append("neumann")
        elif (bc == "dirichlet"):  bcparam.append("dirichlet")
        else:
            msg.fail("invalid BC")


    bc = patch.BCObject(xlb=bcparam[0], xrb=bcparam[1], 
                        ylb=bcparam[2], yrb=bcparam[3])    


    my_data = patch.CellCenterData2d(my_grid, runtime_parameters=rp)

    my_data.register_var("phi", bc)

    my_data.create()

    return my_grid, my_data

def getParam(key):
    """
    returns the value of the runtime parameter corresponding to the
    input key
    """
    if globalParams == {}:
        msg.warning("WARNING: runtime parameters not yet initialized")
        LoadParams("_defaults")

    # debugging
    if not key in usedParams:
        usedParams.append(key)
        
    if key in globalParams.keys():
        return globalParams[key]
    else:
        msg.fail("ERROR: runtime parameter %s not found" % (key))

    def sdc_integral(self, m_start, m_end, As):
        """Compute the integral over the sources from m to m+1 with a
        Simpson's rule"""

        integral = self.cc_data.grid.scratch_array(nvar=self.ivars.nvar)

        if m_start == 0 and m_end == 1:
            for n in range(self.ivars.nvar):
                integral.v(n=n)[:, :] = self.dt/24.0 * (5.0*As[0].v(n=n) + 8.0*As[1].v(n=n) - As[2].v(n=n))

        elif m_start == 1 and m_end == 2:
            for n in range(self.ivars.nvar):
                integral.v(n=n)[:, :] = self.dt/24.0 * (-As[0].v(n=n) + 8.0*As[1].v(n=n) + 5.0*As[2].v(n=n))

        else:
            msg.fail("invalid quadrature range")

        return integral

def init_data(my_data, rp):
    """ initialize the slotted advection problem """
    msg.bold("initializing the slotted advection problem...")

    # make sure that we are passed a valid patch object
    if not isinstance(my_data, patch.CellCenterData2d):
        print(my_data.__class__)
        msg.fail("ERROR: patch invalid in slotted.py")

    offset = rp.get_param("slotted.offset")
    omega = rp.get_param("slotted.omega")

    myg = my_data.grid

    xctr_dens = 0.5*(myg.xmin + myg.xmax)
    yctr_dens = 0.5*(myg.ymin + myg.ymax) + offset

    # setting initial condition for density
    dens = my_data.get_var("density")
    dens[:, :] = 0.0

    R = 0.15
    slot_width = 0.05

    inside = (myg.x2d - xctr_dens)**2 + (myg.y2d - yctr_dens)**2 < R**2

    slot_x = np.logical_and(myg.x2d > (xctr_dens - slot_width*0.5),
                            myg.x2d < (xctr_dens + slot_width*0.5))
    slot_y = np.logical_and(myg.y2d > (yctr_dens - R),
                            myg.y2d < (yctr_dens))
    slot = np.logical_and(slot_x, slot_y)

    dens[inside] = 1.0
    dens[slot] = 0.0

    # setting initial condition for velocity
    u = my_data.get_var("x-velocity")
    v = my_data.get_var("y-velocity")

    u[:, :] = omega*(myg.y2d - xctr_dens)
    v[:, :] = -omega*(myg.x2d - (yctr_dens-offset))

    print("extrema: ", np.amax(u), np.amin(u))

    def pretty_print(self, varname):
        """
        Print out a small dataset to the screen with the ghost cells
        a different color, to make things stand out
        """

        a = self.get_var(varname)

        if self.dtype == numpy.int:
            fmt = "%4d"
        elif self.dtype == numpy.float64:
            fmt = "%10.5g"
        else:
            msg.fail("ERROR: dtype not supported")
        
        # print j descending, so it looks like a grid (y increasing with height)
        j = self.grid.qy-1
        while j >= 0:
            i = 0
            while i < self.grid.qx:

                if (j < self.grid.jlo or j > self.grid.jhi or
                    i < self.grid.ilo or i > self.grid.ihi):
                    gc = 1
                else:
                    gc = 0

                if gc:
                    print("\033[31m" + fmt % (a[i,j]) + "\033[0m", end="")
                else:
                    print (fmt % (a[i,j]), end="")

                i += 1

            print(" ")
            j -= 1

        leg = """
         ^ y
         |
         +---> x
        """
        print(leg)