""" The patch module defines the classes necessary to describe finite-volume data and the grid that it lives on. Typical usage: -- create the grid grid = Grid2d(nx, ny) -- create the data that lives on that grid data = CellCenterData2d(grid) bc = BCObject(xlb="reflect", xrb="reflect", ylb="outflow", yrb="outflow") data.register_var("density", bc) ... data.create() -- initialize some data dens = data.get_var("density") dens[:,:] = ... -- fill the ghost cells data.fill_BC("density") """ from __future__ import print_function import copy import numpy import pickle from util import msg valid = ["outflow", "periodic", "reflect", "reflect-even", "reflect-odd", "dirichlet", "neumann"] extBCs = {} def define_bc(type, function): """ use this to extend the types of boundary conditions supported on a solver-by-solver basis. Here we pass in the reference to a function that can be called with the data that needs to be filled. """ valid.append(type) extBCs[type] = function class BCObject: """ Boundary condition container -- hold the BCs on each boundary for a single variable """ def __init__ (self, xlb="outflow", xrb="outflow", ylb="outflow", yrb="outflow", odd_reflect_dir=""): """ Create the BCObject. Parameters ---------- xlb : {'outflow', 'periodic', 'reflect', 'reflect-even', 'reflect-odd', 'dirichlet', 'neumann', user-defined}, optional The type of boundary condition to enforce on the lower x boundary. user-defined requires one to have defined a new boundary condition type using define_bc() xrb : {'outflow', 'periodic', 'reflect', 'reflect-even', 'reflect-odd', 'dirichlet', 'neumann', user-defined}, optional The type of boundary condition to enforce on the upper x boundary. user-defined requires one to have defined a new boundary condition type using define_bc() ylb : {'outflow', 'periodic', 'reflect', 'reflect-even', 'reflect-odd', 'dirichlet', 'neumann', user-defined}, optional The type of boundary condition to enforce on the lower y boundary. user-defined requires one to have defined a new boundary condition type using define_bc() yrb : {'outflow', 'periodic', 'reflect', 'reflect-even', 'reflect-odd', 'dirichlet', 'neumann', user-defined}, optional The type of boundary condition to enforce on the upper y boundary. user-defined requires one to have defined a new boundary condition type using define_bc() odd_reflect_dir : {'x', 'y'}, optional The direction along which reflection should be odd (sign changes). If not specified, a boundary condition of 'reflect' will always be set to 'reflect-even' """ # note: "reflect" is ambiguous and will be converted into # either reflect-even (the default) or reflect-odd if # odd_reflect_dir specifies the corresponding direction ("x", # "y") # -x boundary if xlb in valid: self.xlb = xlb if self.xlb == "reflect": if odd_reflect_dir == "x": self.xlb = "reflect-odd" else: self.xlb = "reflect-even" else: msg.fail("ERROR: xlb = %s invalid BC" % (xlb)) # +x boundary if xrb in valid: self.xrb = xrb if self.xrb == "reflect": if odd_reflect_dir == "x": self.xrb = "reflect-odd" else: self.xrb = "reflect-even" else: msg.fail("ERROR: xrb = %s invalid BC" % (xrb)) # -y boundary if ylb in valid: self.ylb = ylb if self.ylb == "reflect": if odd_reflect_dir == "y": self.ylb = "reflect-odd" else: self.ylb = "reflect-even" else: msg.fail("ERROR: ylb = %s invalid BC" % (ylb)) # +y boundary if yrb in valid: self.yrb = yrb if self.yrb == "reflect": if odd_reflect_dir == "y": self.yrb = "reflect-odd" else: self.yrb = "reflect-even" else: msg.fail("ERROR: yrb = %s invalid BC" % (yrb)) # periodic checks if ((xlb == "periodic" and not xrb == "periodic") or (xrb == "periodic" and not xlb == "periodic")): msg.fail("ERROR: both xlb and xrb must be periodic") if ((ylb == "periodic" and not yrb == "periodic") or (yrb == "periodic" and not ylb == "periodic")): msg.fail("ERROR: both ylb and yrb must be periodic") def __str__(self): """ print out some basic information about the BC object """ string = "BCs: -x: %s +x: %s -y: %s +y: %s" % \ (self.xlb, self.xrb, self.ylb, self.yrb) return string class Grid2d: """ the 2-d grid class. The grid object will contain the coordinate information (at various centerings). A basic (1-d) representation of the layout is: | | | X | | | | X | | | +--*--+- // -+--*--X--*--+--*--+- // -+--*--+--*--X--*--+- // -+--*--+ 0 ng-1 ng ng+1 ... ng+nx-1 ng+nx 2ng+nx-1 ilo ihi |<- ng guardcells->|<---- nx interior zones ----->|<- ng guardcells->| The '*' marks the data locations. """ def __init__ (self, nx, ny, ng=1, \ xmin=0.0, xmax=1.0, ymin=0.0, ymax=1.0): """ Create a Grid2d object. The only data that we require is the number of points that make up the mesh in each direction. Optionally we take the extrema of the domain (default is [0,1]x[0,1]) and number of ghost cells (default is 1). Note that the Grid2d object only defines the discretization, it does not know about the boundary conditions, as these can vary depending on the variable. Parameters ---------- nx : int Number of zones in the x-direction ny : int Number of zones in the y-direction ng : int, optional Number of ghost cells xmin : float, optional Physical coordinate at the lower x boundary xmax : float, optional Physical coordinate at the upper x boundary ymin : float, optional Physical coordinate at the lower y boundary ymax : float, optional Physical coordinate at the upper y boundary """ # size of grid self.nx = nx self.ny = ny self.ng = ng self.qx = 2*ng+nx self.qy = 2*ng+ny # domain extrema self.xmin = xmin self.xmax = xmax self.ymin = ymin self.ymax = ymax # compute the indices of the block interior (excluding guardcells) self.ilo = ng self.ihi = ng+nx-1 self.jlo = ng self.jhi = ng+ny-1 # define the coordinate information at the left, center, and right # zone coordinates self.dx = (xmax - xmin)/nx self.xl = (numpy.arange(self.qx) - ng)*self.dx + xmin self.xr = (numpy.arange(self.qx) + 1.0 - ng)*self.dx + xmin self.x = 0.5*(self.xl + self.xr) self.dy = (ymax - ymin)/ny self.yl = (numpy.arange(self.qy) - ng)*self.dy + ymin self.yr = (numpy.arange(self.qy) + 1.0 - ng)*self.dy + ymin self.y = 0.5*(self.yl + self.yr) # 2-d versions of the zone coordinates (replace with meshgrid?) x2d = numpy.repeat(self.x, self.qy) x2d.shape = (self.qx, self.qy) self.x2d = x2d y2d = numpy.repeat(self.y, self.qx) y2d.shape = (self.qy, self.qx) y2d = numpy.transpose(y2d) self.y2d = y2d def scratch_array(self, nvar=1): """ return a standard numpy array dimensioned to have the size and number of ghostcells as the parent grid """ if nvar == 1: return numpy.zeros((self.qx, self.qy), dtype=numpy.float64) else: return numpy.zeros((self.qx, self.qy, nvar), dtype=numpy.float64) def __str__(self): """ print out some basic information about the grid object """ string = "2-d grid: nx = " + repr(self.nx) + \ ", ny = " + repr(self.ny) + \ ", ng = " + repr(self.ng) return string def __eq__(self, other): """ are two grids equivalent? """ result = (self.nx == other.nx) and \ (self.ny == other.ny) and \ (self.ng == other.ng) and \ (self.xmin == other.xmin) and \ (self.xmax == other.xmax) and \ (self.ymin == other.ymin) and \ (self.ymax == other.ymax) return result class CellCenterData2d: """ A class to define cell-centered data that lives on a grid. A CellCenterData2d object is built in a multi-step process before it can be used. -- Create the object. We pass in a grid object to describe where the data lives: my_data = patch.CellCenterData2d(myGrid) -- Register any variables that we expect to live on this patch. Here BCObject describes the boundary conditions for that variable. my_data.register_var('density', BCObject) my_data.register_var('x-momentum', BCObject) ... -- Register any auxillary data -- these are any parameters that are needed to interpret the data outside of the simulation (for example, the gamma for the equation of state). my_data.set_aux(keyword, value) -- Finish the initialization of the patch my_data.create() This last step actually allocates the storage for the state variables. Once this is done, the patch is considered to be locked. New variables cannot be added. """ def __init__ (self, grid, dtype=numpy.float64): """ Initialize the CellCenterData2d object. Parameters ---------- grid : Grid2d object The grid upon which the data will live dtype : NumPy data type, optional The datatype of the data we wish to create (defaults to numpy.float64 runtime_parameters : RuntimeParameters object, optional The runtime parameters that go along with this data """ self.grid = grid self.dtype = dtype self.data = None self.vars = [] self.nvar = 0 self.aux = {} self.BCs = {} # time self.t = -1.0 self.initialized = 0 def register_var(self, name, bc_object): """ Register a variable with CellCenterData2d object. Parameters ---------- name : str The variable name bc_object : BCObject object The boundary conditions that describe the actions to take for this variable at the physical domain boundaries. """ if self.initialized == 1: msg.fail("ERROR: grid already initialized") self.vars.append(name) self.nvar += 1 self.BCs[name] = bc_object def set_aux(self, keyword, value): """ Set any auxillary (scalar) data. This data is simply carried along with the CellCenterData2d object Parameters ---------- keyword : str The name of the datum value : any time The value to associate with the keyword """ self.aux[keyword] = value 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, self.grid.qx, self.grid.qy), dtype=self.dtype) self.initialized = 1 def __str__(self): """ print out some basic information about the CellCenterData2d object """ if self.initialized == 0: myStr = "CellCenterData2d object not yet initialized" return myStr myStr = "cc data: nx = " + repr(self.grid.nx) + \ ", ny = " + repr(self.grid.ny) + \ ", ng = " + repr(self.grid.ng) + "\n" + \ " nvars = " + repr(self.nvar) + "\n" + \ " variables: \n" ilo = self.grid.ilo ihi = self.grid.ihi jlo = self.grid.jlo jhi = self.grid.jhi n = 0 while n < self.nvar: myStr += "%16s: min: %15.10f max: %15.10f\n" % \ (self.vars[n], numpy.min(self.data[n,ilo:ihi+1,jlo:jhi+1]), numpy.max(self.data[n,ilo:ihi+1,jlo:jhi+1]) ) myStr += "%16s BCs: -x: %-12s +x: %-12s -y: %-12s +y: %-12s\n" %\ (" " , self.BCs[self.vars[n]].xlb, self.BCs[self.vars[n]].xrb, self.BCs[self.vars[n]].ylb, self.BCs[self.vars[n]].yrb) n += 1 return myStr def get_var(self, name): """ Return a data array for the variable described by name. Any changes made to this are automatically reflected in the CellCenterData2d object. Parameters ---------- name : str The name of the variable to access Returns ------- out : ndarray The array of data corresponding to the variable name """ n = self.vars.index(name) return self.data[n,:,:] def get_var_by_index(self, n): """ Return a data array for the variable with index n in the data array. Any changes made to this are automatically reflected in the CellCenterData2d object. Parameters ---------- n : int The index of the variable to access Returns ------- out : ndarray The array of data corresponding to the index """ return self.data[n,:,:] def get_aux(self, keyword): """ Get the auxillary data associated with keyword Parameters ---------- keyword : str The name of the auxillary data to access Returns ------- out : variable type The value corresponding to the keyword """ if keyword in self.aux.keys(): return self.aux[keyword] else: return None def zero(self, name): """ Zero out the data array associated with variable name. Parameters ---------- name : str The name of the variable to zero """ n = self.vars.index(name) self.data[n,:,:] = 0.0 def fill_BC_all(self): """ Fill boundary conditions on all variables. """ for name in self.vars: self.fill_BC(name) def fill_BC(self, name): """ Fill the boundary conditions. This operates on a single state variable at a time, to allow for maximum flexibility. We do periodic, reflect-even, reflect-odd, and outflow Each variable name has a corresponding BCObject stored in the CellCenterData2d object -- we refer to this to figure out the action to take at each boundary. Parameters ---------- name : str The name of the variable for which to fill the BCs. """ # there is only a single grid, so every boundary is on # a physical boundary (except if we are periodic) # Note: we piggy-back on outflow and reflect-odd for # Neumann and Dirichlet homogeneous BCs respectively, but # this only works for a single ghost cell n = self.vars.index(name) # -x boundary if self.BCs[name].xlb == "outflow" or self.BCs[name].xlb == "neumann": i = 0 while i < self.grid.ilo: self.data[n,i,:] = self.data[n,self.grid.ilo,:] i += 1 elif self.BCs[name].xlb == "reflect-even": i = 0 while i < self.grid.ilo: self.data[n,i,:] = self.data[n,2*self.grid.ng-i-1,:] i += 1 elif (self.BCs[name].xlb == "reflect-odd" or self.BCs[name].xlb == "dirichlet"): i = 0 while i < self.grid.ilo: self.data[n,i,:] = -self.data[n,2*self.grid.ng-i-1,:] i += 1 elif self.BCs[name].xlb == "periodic": i = 0 while i < self.grid.ilo: self.data[n,i,:] = self.data[n,self.grid.ihi-self.grid.ng+i+1,:] i += 1 # +x boundary if self.BCs[name].xrb == "outflow" or self.BCs[name].xrb == "neumann": i = self.grid.ihi+1 while i < self.grid.nx+2*self.grid.ng: self.data[n,i,:] = self.data[n,self.grid.ihi,:] i += 1 elif self.BCs[name].xrb == "reflect-even": i = 0 while i < self.grid.ng: i_bnd = self.grid.ihi+1+i i_src = self.grid.ihi-i self.data[n,i_bnd,:] = self.data[n,i_src,:] i += 1 elif (self.BCs[name].xrb == "reflect-odd" or self.BCs[name].xrb == "dirichlet"): i = 0 while i < self.grid.ng: i_bnd = self.grid.ihi+1+i i_src = self.grid.ihi-i self.data[n,i_bnd,:] = -self.data[n,i_src,:] i += 1 elif self.BCs[name].xrb == "periodic": i = self.grid.ihi+1 while i < 2*self.grid.ng + self.grid.nx: self.data[n,i,:] = self.data[n,i-self.grid.ihi-1+self.grid.ng,:] i += 1 # -y boundary if self.BCs[name].ylb == "outflow" or self.BCs[name].ylb == "neumann": j = 0 while j < self.grid.jlo: self.data[n,:,j] = self.data[n,:,self.grid.jlo] j += 1 elif self.BCs[name].ylb == "reflect-even": j = 0 while j < self.grid.jlo: self.data[n,:,j] = self.data[n,:,2*self.grid.ng-j-1] j += 1 elif (self.BCs[name].ylb == "reflect-odd" or self.BCs[name].ylb == "dirichlet"): j = 0 while j < self.grid.jlo: self.data[n,:,j] = -self.data[n,:,2*self.grid.ng-j-1] j += 1 elif self.BCs[name].ylb == "periodic": j = 0 while j < self.grid.jlo: self.data[n,:,j] = self.data[n,:,self.grid.jhi-self.grid.ng+j+1] j += 1 else: if self.BCs[name].ylb in extBCs.keys(): extBCs[self.BCs[name].ylb](self.BCs[name].ylb, "ylb", name, self) # +y boundary if self.BCs[name].yrb == "outflow" or self.BCs[name].yrb == "neumann": j = self.grid.jhi+1 while j < self.grid.ny+2*self.grid.ng: self.data[n,:,j] = self.data[n,:,self.grid.jhi] j += 1 elif self.BCs[name].yrb == "reflect-even": j = 0 while j < self.grid.ng: j_bnd = self.grid.jhi+1+j j_src = self.grid.jhi-j self.data[n,:,j_bnd] = self.data[n,:,j_src] j += 1 elif (self.BCs[name].yrb == "reflect-odd" or self.BCs[name].yrb == "dirichlet"): j = 0 while j < self.grid.ng: j_bnd = self.grid.jhi+1+j j_src = self.grid.jhi-j self.data[n,:,j_bnd] = -self.data[n,:,j_src] j += 1 elif self.BCs[name].yrb == "periodic": j = self.grid.jhi+1 while j < 2*self.grid.ng + self.grid.ny: self.data[n,:,j] = self.data[n,:,j-self.grid.jhi-1+self.grid.ng] j += 1 else: if self.BCs[name].yrb in extBCs.keys(): extBCs[self.BCs[name].yrb](self.BCs[name].yrb, "yrb", name, self) def restrict(self, varname): """ Restrict the variable varname to a coarser grid (factor of 2 coarser) and return an array with the resulting data (and same number of ghostcells) """ fG = self.grid fData = self.get_var(varname) # allocate an array for the coarsely gridded data ng_c = fG.ng nx_c = fG.nx/2 ny_c = fG.ny/2 cData = numpy.zeros((2*ng_c+nx_c, 2*ng_c+ny_c), dtype=self.dtype) ilo_c = ng_c ihi_c = ng_c+nx_c-1 jlo_c = ng_c jhi_c = ng_c+ny_c-1 # fill the coarse array with the restricted data -- just # average the 4 fine cells into the corresponding coarse cell # that encompasses them. # This is done by shifting our view into the fData array and # using a stride of 2 in the indexing. cData[ilo_c:ihi_c+1,jlo_c:jhi_c+1] = \ 0.25*(fData[fG.ilo :fG.ihi+1:2,fG.jlo :fG.jhi+1:2] + fData[fG.ilo+1:fG.ihi+1:2,fG.jlo :fG.jhi+1:2] + fData[fG.ilo :fG.ihi+1:2,fG.jlo+1:fG.jhi+1:2] + fData[fG.ilo+1:fG.ihi+1:2,fG.jlo+1:fG.jhi+1:2]) return cData def prolong(self, varname): """ Prolong the data in the current (coarse) grid to a finer (factor of 2 finer) grid. Return an array with the resulting data (and same number of ghostcells). Only the data for the variable varname will be operated upon. We will reconstruct the data in the zone from the zone-averaged variables using the same limited slopes as in the advection routine. Getting a good multidimensional reconstruction polynomial is hard -- we want it to be bilinear and monotonic -- we settle for having each slope be independently monotonic: (x) (y) f(x,y) = m x/dx + m y/dy + where the m's are the limited differences in each direction. When averaged over the parent cell, this reproduces . Each zone's reconstrution will be averaged over 4 children. +-----------+ +-----+-----+ | | | | | | | | 3 | 4 | | | --> +-----+-----+ | | | | | | | | 1 | 2 | +-----------+ +-----+-----+ We will fill each of the finer resolution zones by filling all the 1's together, using a stride 2 into the fine array. Then the 2's and ..., this allows us to operate in a vector fashion. All operations will use the same slopes for their respective parents. """ cG = self.grid cData = self.get_var(varname) # allocate an array for the coarsely gridded data ng_f = cG.ng nx_f = cG.nx*2 ny_f = cG.ny*2 fData = numpy.zeros((2*ng_f+nx_f, 2*ng_f+ny_f), dtype=self.dtype) ilo_f = ng_f ihi_f = ng_f+nx_f-1 jlo_f = ng_f jhi_f = ng_f+ny_f-1 # slopes for the coarse data m_x = cG.scratch_array() m_x[cG.ilo:cG.ihi+1,cG.jlo:cG.jhi+1] = \ 0.5*(cData[cG.ilo+1:cG.ihi+2,cG.jlo:cG.jhi+1] - cData[cG.ilo-1:cG.ihi ,cG.jlo:cG.jhi+1]) m_y = cG.scratch_array() m_y[cG.ilo:cG.ihi+1,cG.jlo:cG.jhi+1] = \ 0.5*(cData[cG.ilo:cG.ihi+1,cG.jlo+1:cG.jhi+2] - cData[cG.ilo:cG.ihi+1,cG.jlo-1:cG.jhi ]) # fill the '1' children fData[ilo_f:ihi_f+1:2,jlo_f:jhi_f+1:2] = \ cData[cG.ilo:cG.ihi+1,cG.jlo:cG.jhi+1] \ - 0.25*m_x[cG.ilo:cG.ihi+1,cG.jlo:cG.jhi+1] \ - 0.25*m_y[cG.ilo:cG.ihi+1,cG.jlo:cG.jhi+1] # fill the '2' children fData[ilo_f+1:ihi_f+1:2,jlo_f:jhi_f+1:2] = \ cData[cG.ilo:cG.ihi+1,cG.jlo:cG.jhi+1] \ + 0.25*m_x[cG.ilo:cG.ihi+1,cG.jlo:cG.jhi+1] \ - 0.25*m_y[cG.ilo:cG.ihi+1,cG.jlo:cG.jhi+1] # fill the '3' children fData[ilo_f:ihi_f+1:2,jlo_f+1:jhi_f+1:2] = \ cData[cG.ilo:cG.ihi+1,cG.jlo:cG.jhi+1] \ - 0.25*m_x[cG.ilo:cG.ihi+1,cG.jlo:cG.jhi+1] \ + 0.25*m_y[cG.ilo:cG.ihi+1,cG.jlo:cG.jhi+1] # fill the '4' children fData[ilo_f+1:ihi_f+1:2,jlo_f+1:jhi_f+1:2] = \ cData[cG.ilo:cG.ihi+1,cG.jlo:cG.jhi+1] \ + 0.25*m_x[cG.ilo:cG.ihi+1,cG.jlo:cG.jhi+1] \ + 0.25*m_y[cG.ilo:cG.ihi+1,cG.jlo:cG.jhi+1] return fData def write(self, filename): """ write out the CellCenterData2d object to disk, stored in the file filename. We use a python binary format (via pickle). This stores a representation of the entire object. """ pF = open(filename + ".pyro", "wb") pickle.dump(self, pF, pickle.HIGHEST_PROTOCOL) pF.close() 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) # backwards compatibility ccData2d = CellCenterData2d grid2d = Grid2d bcObject = BCObject def read(filename): """ Read a CellCenterData object from a file and return it and the grid info and data. """ # if we come in with .pyro, we don't need to add it again if filename.find(".pyro") < 0: filename += ".pyro" pF = open(filename, "rb") data = pickle.load(pF) pF.close() return data.grid, data 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) n = 0 while n < old.nvar: new.register_var(old.vars[n], old.BCs[old.vars[n]]) n += 1 new.create() new.aux = old.aux.copy() new.data = old.data.copy() return new if __name__== "__main__": # illustrate basic mesh operations myg = Grid2d(16,32, xmax=1.0, ymax=2.0) mydata = CellCenterData2d(myg) bc = BCObject() mydata.register_var("a", bc) mydata.create() a = mydata.get_var("a") a[:,:] = numpy.exp(-(myg.x2d - 0.5)**2 - (myg.y2d - 1.0)**2) print(mydata) # output print("writing\n") mydata.write("mesh_test") print("reading\n") myg2, myd2 = read("mesh_test") print(myd2) mydata.pretty_print("a")