this notebook illustrates the basic ways of interacting with the pyro2 mesh module. We create some data that lives on a grid and show how to fill the ghost cells. The pretty_print() function shows us that they work as expected.
In [1]:
from __future__ import print_function
import numpy as np
import mesh.boundary as bnd
import mesh.patch as patch
import matplotlib.pyplot as plt
%matplotlib inline
# for unit testing, we want to ensure the same random numbers
np.random.seed(100)
There are a few core classes that we deal with when creating a grid with associated variables:
Grid2d : this holds the size of the grid (in zones) and the physical coordinate information, including coordinates of cell edges and centers
BC : this is a container class that simply holds the type of boundary condition on each domain edge.
ArrayIndexer : this is an array of data along with methods that know how to access it with different offsets into the data that usually arise in stencils (like {i+1, j})
CellCenterData2d : this holds the data that lives on a grid. Each variable that is part of this class has its own boundary condition type.
We start by creating a Grid2d object with 4 x 6 cells and 2 ghost cells
In [2]:
g = patch.Grid2d(4, 6, ng=2)
print(g)
In [3]:
help(g)
Then create a dataset that lives on this grid and add a variable name. For each variable that lives on the grid, we need to define the boundary conditions -- this is done through the BC object.
In [4]:
bc = bnd.BC(xlb="periodic", xrb="periodic", ylb="reflect", yrb="outflow")
print(bc)
In [5]:
d = patch.CellCenterData2d(g)
d.register_var("a", bc)
d.create()
print(d)
Now we fill the grid with random data. get_var() returns an ArrayIndexer object that has methods for accessing views into the data. Here we use a.v() to get the "valid" region, i.e. excluding ghost cells.
In [6]:
a = d.get_var("a")
a.v()[:,:] = np.random.rand(g.nx, g.ny)
when we pretty_print() the variable, we see the ghost cells colored red. Note that we just filled the interior above.
In [7]:
a.pretty_print()
pretty_print() can also take an argumet, specifying the format string to be used for the output.
In [8]:
a.pretty_print(fmt="%7.3g")
now fill the ghost cells -- notice that the left and right are periodic, the upper is outflow, and the lower is reflect, as specified when we registered the data above.
In [9]:
d.fill_BC("a")
a.pretty_print()
We can find the L2 norm of the data easily
In [10]:
a.norm()
Out[10]:
and the min and max
In [11]:
print(a.min(), a.max())
We we access the data, an ArrayIndexer object is returned. The ArrayIndexer sub-classes the NumPy ndarray, so it can do all of the methods that a NumPy array can, but in addition, we can use the ip(), jp(), or ipjp() methods to the ArrayIndexer object shift our view in the x, y, or x & y directions.
To make this clearer, we'll change our data set to be nicely ordered numbers. We index the ArrayIndex the same way we would a NumPy array. The index space includes ghost cells, so the ilo and ihi attributes from the grid object are useful to index just the valid region. The .v() method is a shortcut that also gives a view into just the valid data.
Note: when we use one of the ip(), jp(), ipjp(), or v() methods, the result is a regular NumPy ndarray, not an ArrayIndexer object. This is because it only spans part of the domain (e.g., no ghost cells), and therefore cannot be associated with the Grid2d object that the ArrayIndexer is built from.
In [12]:
type(a)
Out[12]:
In [13]:
type(a.v())
Out[13]:
In [14]:
a[:,:] = np.arange(g.qx*g.qy).reshape(g.qx, g.qy)
In [15]:
a.pretty_print()
We index our arrays as {i,j}, so x (indexed by i) is the row and y (indexed by j) is the column in the NumPy array. Note that python arrays are stored in row-major order, which means that all of the entries in the same row are adjacent in memory. This means that when we simply print out the ndarray, we see constant-x horizontally, which is the transpose of what we are used to.
In [16]:
a.v()
Out[16]:
We can offset our view into the array by one in x -- this would be like {i+1, j} when we loop over data. The ip() method is used here, and takes an argument which is the (positive) shift in the x (i) direction. So here's a shift by 1
In [17]:
a.ip(-1, buf=1)
Out[17]:
A shifted view is necessarily smaller than the original array, and relies on ghost cells to bring new data into view. Because of this, the underlying data is no longer the same size as the original data, so we return it as an ndarray (which is actually just a view into the data in the ArrayIndexer object, so no copy is made.
To see that it is simply a view, lets shift and edit the data
In [18]:
d = a.ip(1)
d[1,1] = 0.0
a.pretty_print()
Here, since d was really a view into $a_{i+1,j}$, and we accessed element (1,1) into that view (with 0,0 as the origin), we were really accessing the element (2,1) in the valid region
ArrayIndexer objects are easy to use to construct differences, like those that appear in a stencil for a finite-difference, without having to explicitly loop over the elements of the array.
Here's we'll create a new dataset that is initialized with a sine function
In [19]:
g = patch.Grid2d(8, 8, ng=2)
d = patch.CellCenterData2d(g)
bc = bnd.BC(xlb="periodic", xrb="periodic", ylb="periodic", yrb="periodic")
d.register_var("a", bc)
d.create()
a = d.get_var("a")
a[:,:] = np.sin(2.0*np.pi*a.g.x2d)
d.fill_BC("a")
Our grid object can provide us with a scratch array (an ArrayIndexer object) define on the same grid
In [20]:
b = g.scratch_array()
type(b)
Out[20]:
We can then fill the data in this array with differenced data from our original array -- since b has a separate data region in memory, its elements are independent of a. We do need to make sure that we have the same number of elements on the left and right of the =. Since by default, ip() will return a view with the same size as the valid region, we can use .v() on the left to accept the differences.
Here we compute a centered-difference approximation to the first derivative
In [21]:
b.v()[:,:] = (a.ip(1) - a.ip(-1))/(2.0*a.g.dx)
# normalization was 2.0*pi
b[:,:] /= 2.0*np.pi
In [22]:
plt.plot(g.x[g.ilo:g.ihi+1], a[g.ilo:g.ihi+1,a.g.jc])
plt.plot(g.x[g.ilo:g.ihi+1], b[g.ilo:g.ihi+1,b.g.jc])
print (a.g.dx)
we can get a new ArrayIndexer object on a coarser grid for one of our variables
In [23]:
c = d.restrict("a")
In [24]:
c.pretty_print()
or a finer grid
In [25]:
f = d.prolong("a")
In [26]:
f.pretty_print(fmt="%6.2g")