"""

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 numpy as np

import pickle

from util import msg

import mesh.array_indexer as ai

# keep track of whether the BCs are solid walls (passed into the

# Riemann solver).

bc_props = {}

bc_props["outflow"] = False

bc_props["periodic"] = False

bc_props["reflect"] = True

bc_props["reflect-even"] = True

bc_props["reflect-odd"] = True

bc_props["dirichlet"] = True

bc_props["neumann"] = False

ext_bcs = {}

def define_bc(bc_type, function, is_solid=False):

"""

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. is_solid indicates whether it should be interpreted as

a solid wall (no flux through the BC)"

"""

bc_props[bc_type] = is_solid

ext_bcs[bc_type] = function

def _set_reflect(odd_reflect_dir, dir_string):

if odd_reflect_dir == dir_string:

return "reflect-odd"

else:

return "reflect-even"

class BCObject(object):

"""Boundary condition container -- hold the BCs on each boundary

for a single variable.

For Neumann and Dirichlet BCs, a function callback can be stored

for inhomogeous BCs. This function should provide the value on

the physical boundary (not cell center). This is evaluated on the

relevant edge when the __init__ routine is called. For this

reason, you need to pass in a grid object. Note: this only

ensures that the first ghost cells is consistent with the BC

value.

"""

def __init__(self,

xlb="outflow", xrb="outflow",

ylb="outflow", yrb="outflow",

xl_func=None, xr_func=None,

yl_func=None, yr_func=None,

grid=None,

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'

xl_func : function, optional

A function, f(y), that provides the value of the

Dirichlet or Neumann BC on the -x physical boundary.

xr_func : function, optional

A function, f(y), that provides the value of the

Dirichlet or Neumann BC on the +x physical boundary.

yl_func : function, optional

A function, f(x), that provides the value of the

Dirichlet or Neumann BC on the -y physical boundary.

yr_func : function, optional

A function, f(x), that provides the value of the

Dirichlet or Neumann BC on the +y physical boundary.

grid : a Grid2d object, optional

The grid object is used for evaluating the function

to define the boundary values for inhomogeneous

Dirichlet and Neumann BCs. It is required if

any functions are passed in.

"""

# 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")

valid = list(bc_props.keys())

# -x boundary

if xlb in valid:

self.xlb = xlb

if self.xlb == "reflect":

self.xlb = _set_reflect(odd_reflect_dir, "x")

else:

msg.fail("ERROR: xlb = %s invalid BC" % (xlb))

# +x boundary

if xrb in valid:

self.xrb = xrb

if self.xrb == "reflect":

self.xrb = _set_reflect(odd_reflect_dir, "x")

else:

msg.fail("ERROR: xrb = %s invalid BC" % (xrb))

# -y boundary

if ylb in valid:

self.ylb = ylb

if self.ylb == "reflect":

self.ylb = _set_reflect(odd_reflect_dir, "y")

else:

msg.fail("ERROR: ylb = %s invalid BC" % (ylb))

# +y boundary

if yrb in valid:

self.yrb = yrb

if self.yrb == "reflect":

self.yrb = _set_reflect(odd_reflect_dir, "y")

else:

msg.fail("ERROR: yrb = %s invalid BC" % (yrb))

# periodic checks

if ((xlb == "periodic" and xrb != "periodic") or

(xrb == "periodic" and xlb != "periodic")):

msg.fail("ERROR: both xlb and xrb must be periodic")

if ((ylb == "periodic" and yrb != "periodic") or

(yrb == "periodic" and ylb != "periodic")):

msg.fail("ERROR: both ylb and yrb must be periodic")

# inhomogeneous functions for Dirichlet or Neumann

self.xl_value = self.xr_value = self.yl_value = self.yr_value = None

if xl_func is not None:

self.xl_value = xl_func(grid.y)

if xr_func is not None:

self.xr_value = xr_func(grid.y)

if yl_func is not None:

self.yl_value = yl_func(grid.x)

if yr_func is not None:

self.yr_value = yr_func(grid.x)

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

# center of the grid (for convenience)

self.ic = self.ilo + nx/2 - 1

self.jc = self.jlo + ny/2 - 1

# define the coordinate information at the left, center, and right

# zone coordinates

self.dx = (xmax - xmin)/nx

self.xl = (np.arange(self.qx) - ng)*self.dx + xmin

self.xr = (np.arange(self.qx) + 1.0 - ng)*self.dx + xmin

self.x = 0.5*(self.xl + self.xr)

self.dy = (ymax - ymin)/ny

self.yl = (np.arange(self.qy) - ng)*self.dy + ymin

self.yr = (np.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 = np.repeat(self.x, self.qy)

x2d.shape = (self.qx, self.qy)

self.x2d = x2d

y2d = np.repeat(self.y, self.qx)

y2d.shape = (self.qy, self.qx)

y2d = np.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:

_tmp = np.zeros((self.qx, self.qy), dtype=np.float64)

else:

_tmp = np.zeros((self.qx, self.qy, nvar), dtype=np.float64)

return ai.ArrayIndexer(d=_tmp, grid=self)

def norm(self, d):

"""

find the norm of the quantity d defined on the same grid, in the

domain's valid region

"""

return np.sqrt(self.dx*self.dy*

np.sum((d[self.ilo:self.ihi+1,self.jlo:self.jhi+1]**2).flat))

def coarse_like(self, N):

"""

return a new grid object coarsened by a factor n, but with

all the other properties the same

"""

return Grid2d(self.nx/N, self.ny/N, ng=self.ng,

xmin=self.xmin, xmax=self.xmax,

ymin=self.ymin, ymax=self.ymax)

def fine_like(self, N):

"""

return a new grid object finer by a factor n, but with

all the other properties the same

"""

return Grid2d(self.nx*N, self.ny*N, ng=self.ng,

xmin=self.xmin, xmax=self.xmax,

ymin=self.ymin, ymax=self.ymax)

def __str__(self):

""" print out some basic information about the grid object """

return "2-d grid: nx = {}, ny = {}, ng = {}".format(

self.nx, self.ny, self.ng)

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=np.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

np.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 = np.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:

my_str = "CellCenterData2d object not yet initialized"

return my_str

my_str = "cc data: nx = {}, ny = {}, ng = {}\n".format(

self.grid.nx, self.grid.ny, self.grid.ng)

my_str += " nvars = {}\n".format(self.nvar)

my_str += " variables:\n"

ilo = self.grid.ilo

ihi = self.grid.ihi

jlo = self.grid.jlo

jhi = self.grid.jhi

for n in range(self.nvar):

my_str += "%16s: min: %15.10f max: %15.10f\n" % \

(self.vars[n],

np.min(self.data[n,ilo:ihi+1,jlo:jhi+1]),

np.max(self.data[n,ilo:ihi+1,jlo:jhi+1]) )

my_str += "%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)

return my_str

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 ai.ArrayIndexer(d=self.data[n,:,:], grid=self.grid)

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 ai.ArrayIndexer(d=self.data[n,:,:], grid=self.grid)

def get_vars(self):

"""

Return the entire data array. Any changes made to this

are automatically reflected in the CellCenterData2d object.

Returns

-------

out : ndarray

The array of data

"""

return ai.ArrayIndexer(d=self.data, grid=self.grid)

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 in ["outflow", "neumann"]:

if self.BCs[name].xl_value is None:

for i in range(self.grid.ilo):

self.data[n,i,:] = self.data[n,self.grid.ilo,:]

else:

self.data[n,self.grid.ilo-1,:] = \

self.data[n,self.grid.ilo,:] - self.grid.dx*self.BCs[name].xl_value[:]

elif self.BCs[name].xlb == "reflect-even":

for i in range(self.grid.ilo):

self.data[n,i,:] = self.data[n,2*self.grid.ng-i-1,:]

elif self.BCs[name].xlb in ["reflect-odd", "dirichlet"]:

if self.BCs[name].xl_value is None:

for i in range(self.grid.ilo):

self.data[n,i,:] = -self.data[n,2*self.grid.ng-i-1,:]

else:

self.data[n,self.grid.ilo-1,:] = \

2*self.BCs[name].xl_value[:] - self.data[n,self.grid.ilo,:]

elif self.BCs[name].xlb == "periodic":

for i in range(self.grid.ilo):

self.data[n,i,:] = self.data[n,self.grid.ihi-self.grid.ng+i+1,:]

# +x boundary

if self.BCs[name].xrb in ["outflow", "neumann"]:

if self.BCs[name].xr_value is None:

for i in range(self.grid.ihi+1, self.grid.nx+2*self.grid.ng):

self.data[n,i,:] = self.data[n,self.grid.ihi,:]

else:

self.data[n,self.grid.ihi+1,:] = \

self.data[n,self.grid.ihi,:] + self.grid.dx*self.BCs[name].xr_value[:]

elif self.BCs[name].xrb == "reflect-even":

for i in range(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,:]

elif self.BCs[name].xrb in ["reflect-odd", "dirichlet"]:

if self.BCs[name].xr_value is None:

for i in range(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,:]

else:

self.data[n,self.grid.ihi+1,:] = \

2*self.BCs[name].xr_value[:] - self.data[n,self.grid.ihi,:]

elif self.BCs[name].xrb == "periodic":

for i in range(self.grid.ihi+1, 2*self.grid.ng + self.grid.nx):

self.data[n,i,:] = self.data[n,i-self.grid.ihi-1+self.grid.ng,:]

# -y boundary

if self.BCs[name].ylb in ["outflow", "neumann"]:

if self.BCs[name].yl_value is None:

for j in range(self.grid.jlo):

self.data[n,:,j] = self.data[n,:,self.grid.jlo]

else:

self.data[n,:,self.grid.jlo-1] = \

self.data[n,:,self.grid.jlo] - self.grid.dy*self.BCs[name].yl_value[:]

elif self.BCs[name].ylb == "reflect-even":

for j in range(self.grid.jlo):

self.data[n,:,j] = self.data[n,:,2*self.grid.ng-j-1]

elif self.BCs[name].ylb in ["reflect-odd", "dirichlet"]:

if self.BCs[name].yl_value is None:

for j in range(self.grid.jlo):

self.data[n,:,j] = -self.data[n,:,2*self.grid.ng-j-1]

else:

self.data[n,:,self.grid.jlo-1] = \

2*self.BCs[name].yl_value[:] - self.data[n,:,self.grid.jlo]

elif self.BCs[name].ylb == "periodic":

for j in range(self.grid.jlo):

self.data[n,:,j] = self.data[n,:,self.grid.jhi-self.grid.ng+j+1]

else:

if self.BCs[name].ylb in ext_bcs.keys():

ext_bcs[self.BCs[name].ylb](self.BCs[name].ylb, "ylb", name, self)

# +y boundary

if self.BCs[name].yrb in ["outflow", "neumann"]:

if self.BCs[name].yr_value is None:

for j in range(self.grid.jhi+1, self.grid.ny+2*self.grid.ng):

self.data[n,:,j] = self.data[n,:,self.grid.jhi]

else:

self.data[n,:,self.grid.jhi+1] = \

self.data[n,:,self.grid.jhi] + self.grid.dy*self.BCs[name].yr_value[:]

elif self.BCs[name].yrb == "reflect-even":

for j in range(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]

elif self.BCs[name].yrb in ["reflect-odd", "dirichlet"]:

if self.BCs[name].yr_value is None:

for j in range(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]

else:

self.data[n,:,self.grid.jhi+1] = \

2*self.BCs[name].yr_value[:] - self.data[n,:,self.grid.jhi]

elif self.BCs[name].yrb == "periodic":

for j in range(self.grid.jhi+1, 2*self.grid.ng + self.grid.ny):

self.data[n,:,j] = self.data[n,:,j-self.grid.jhi-1+self.grid.ng]

else:

if self.BCs[name].yrb in ext_bcs.keys():

ext_bcs[self.BCs[name].yrb](self.BCs[name].yrb, "yrb", name, self)

def min(self, name, ng=0):

"""

return the minimum of the variable name in the domain's valid region

"""

n = self.vars.index(name)

g = self.grid

return np.min(self.data[n,g.ilo-ng:g.ihi+1+ng,g.jlo-ng:g.jhi+1+ng])

def max(self, name, ng=0):

"""

return the maximum of the variable name in the domain's valid region

"""

n = self.vars.index(name)

g = self.grid

return np.max(self.data[n,g.ilo-ng:g.ihi+1+ng,g.jlo-ng:g.jhi+1+ng])

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)

"""

fine_grid = self.grid

fdata = self.get_var(varname)

# allocate an array for the coarsely gridded data

coarse_grid = fine_grid.coarse_like(2)

cdata = coarse_grid.scratch_array()

# fill the coarse array with the restricted data -- just

# average the 4 fine cells into the corresponding coarse cell

# that encompasses them.

cdata.v()[:,:] = \

0.25*(fdata.v(s=2) + fdata.ip(1, s=2) +

fdata.jp(1, s=2) + fdata.ip_jp(1, 1, s=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 + <f>

where the m's are the limited differences in each direction.

When averaged over the parent cell, this reproduces <f>.

Each zone's reconstrution will be averaged over 4 children.

+-----------+ +-----+-----+

| | | | |

| | | 3 | 4 |

| <f> | --> +-----+-----+

| | | | |

| | | 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.

"""

coarse_grid = self.grid

cdata = self.get_var(varname)

# allocate an array for the finely gridded data

fine_grid = coarse_grid.fine_like(2)

fdata = fine_grid.scratch_array()

# slopes for the coarse data

m_x = coarse_grid.scratch_array()

m_x.v()[:,:] = 0.5*(cdata.ip(1) - cdata.ip(-1))

m_y = coarse_grid.scratch_array()

m_y.v()[:,:] = 0.5*(cdata.jp(1) - cdata.jp(-1))

# fill the children

fdata.v(s=2)[:,:] = cdata.v() - 0.25*m_x.v() - 0.25*m_y.v() # 1 child

fdata.ip(1, s=2)[:,:] = cdata.v() + 0.25*m_x.v() - 0.25*m_y.v() # 2

fdata.jp(1, s=2)[:,:] = cdata.v() - 0.25*m_x.v() + 0.25*m_y.v() # 3

fdata.ip_jp(1, 1, s=2)[:,:] = cdata.v() + 0.25*m_x.v() + 0.25*m_y.v() # 4

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, var, fmt=None):

a = self.get_var(var)

a.pretty_print(fmt=fmt)

# 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)

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 do_demo():

# illustrate basic mesh operations

myg = Grid2d(8, 16, xmax=1.0, ymax=2.0)

mydata = CellCenterData2d(myg)

bc = BCObject()

mydata.register_var("a", bc)

mydata.create()

a = mydata.get_var("a")

a.d[:,:] = np.exp(-(myg.x2d - 0.5)**2 - (myg.y2d - 1.0)**2)

print(mydata)

# output