Interactive online version: .

# Grids: Uniformly-Spaced Cartesian Grids¶

Grids are a datastructure that represent one or more physical quantities that share spatial coordinates. For example, the density or magnetic field in a plasma as specified on a Cartesian grid. In addition to storing data, grids have built-in interpolator functions for estimating the values of quantities in between grid vertices.

## Creating a grid¶

[1]:

%matplotlib inline

import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np

from plasmapy.plasma import grids


A grid can be created either by providing three arrays of spatial coordinates for vertices (eg. x,yz positions) or using a np.linspace-like syntax. For example, the two following methods are equivalent:

[2]:

# Method 1
xaxis, yaxis, zaxis = [np.linspace(-1 * u.cm, 1 * u.cm, num=20)] * 3
x, y, z = np.meshgrid(xaxis, yaxis, zaxis, indexing="ij")
grid = grids.CartesianGrid(x, y, z)

# Method 2
grid = grids.CartesianGrid(
np.array([-1, -1, -1]) * u.cm, np.array([1, 1, 1]) * u.cm, num=(150, 150, 150)
)


[3]:

print(f"Is the grid uniformly spaced? {grid.is_uniform}")
print(f"Grid shape: {grid.shape}")
print(f"Grid units: {grid.units}")
print(f"Grid spacing on xaxis: {grid.dax0:.2f}")

Is the grid uniformly spaced? True
Grid shape: (150, 150, 150)
Grid units: [Unit("cm"), Unit("cm"), Unit("cm")]
Grid spacing on xaxis: 0.01 cm


The grid points themselves can be explicitly accessed in one of two forms

[4]:

x, y, z = grid.grids
x.shape

[4]:

(150, 150, 150)

[5]:

xyz = grid.grid
xyz.shape

[5]:

(150, 150, 150, 3)


And the axes can be accessed similarly.

[6]:

xaxis = grid.ax0
xaxis.shape

[6]:

(150,)


Now that the grid has been initialized, we can add quantities to it that represent physical properties defined on the grid vertices. Each quantity must be a u.Quantity array of the same shape as the grid.

[7]:

Ex = np.random.rand(*grid.shape) * u.V / u.m
Ey = np.random.rand(*grid.shape) * u.V / u.m
Ez = np.random.rand(*grid.shape) * u.V / u.m
Bz = np.random.rand(*grid.shape) * u.T
Bz.shape

[7]:

(150, 150, 150)


When quantities are added to the grid, they are associated with a key string (just like a dictionary). Any key string can be used, but PlasmaPy functions use a shared set of recognized quantities to automatically interperet quantities. Each entry is stored as a namedtuple with fields (“key”, “description”, “unit”). The full list of recognized quantities can be accessed in the module:

[8]:

for key in grid.recognized_quantities.keys():
rk = grid.recognized_quantities[key]
key, description, unit = rk.key, rk.description, rk.unit
print(f"{key} -> {description} ({unit})")

x -> x spatial position (m)
y -> y spatial position (m)
z -> z spatial position (m)
rho -> Mass density (kg / m3)
E_x -> Electric field (x component) (V / m)
E_y -> Electric field (y component) (V / m)
E_z -> Electric field (z component) (V / m)
B_x -> Magnetic field (x component) (T)
B_y -> Magnetic field (y component) (T)
B_z -> Magnetic field (z component) (T)
phi -> Electric Scalar Potential (V)


Quantities can be added to the grid as keyword arguments. The keyword becomes the key string for the quantity in the dataset.

[9]:

grid.add_quantities(B_z=Bz)


Adding an unrecognized quantity will lead to a warning, but the quantity will still be added to the grid and can still be accessed by the user later.

[10]:

custom_quantity = np.random.rand(*grid.shape) * u.T * u.mm

/home/docs/checkouts/readthedocs.org/user_builds/plasmapy/envs/stable/lib/python3.7/site-packages/ipykernel_launcher.py:2: UserWarning: Warning: int_B is not recognized quantity key



A summary of the grid, including the currently-defined quantities, can be produced by printing the grid object

[11]:

print(grid)

*** Grid Summary ***
<class 'plasmapy.plasma.grids.CartesianGrid'>
Dimensions: (ax0: 150, ax1: 150, ax2: 150)
Uniformly Spaced: (dax0, dax1, dax2) = (0.013 cm, 0.013 cm, 0.013 cm)
-----------------------------
Coordinates:
-> ax0 (cm) float64 (150,)
-> ax1 (cm) float64 (150,)
-> ax2 (cm) float64 (150,)
-----------------------------
Recognized Quantities:
-> B_z (T) float64 (150, 150, 150)
-> E_x (V / m) float64 (150, 150, 150)
-> E_y (V / m) float64 (150, 150, 150)
-> E_z (V / m) float64 (150, 150, 150)
-----------------------------
Unrecognized Quantities:
-> int_B (mm T) float64 (150, 150, 150)



A simple list of the defined quantity keys on the grid can also be easily accessed

[12]:

print(grid.quantities)

['B_z', 'E_x', 'E_y', 'E_z', 'int_B']


## Methods¶

A number of methods are built into grid objects, and are illustrated here.

The grid.on_grid method determines which points in an array are within the bounds of the grid. Since our example grid is a cube spanning from -1 to -1 cm on each axis, the first of the following points is on the grid while the second is not.

[13]:

pos = np.array([[0.1, -0.3, 0], [3, 0, 0]]) * u.cm
print(grid.on_grid(pos))

[ True False]


Similarly, the grid.vector_intersects function determines whether the line between two points passes through the grid.

[14]:

pt0 = np.array([3, 0, 0]) * u.cm
pt1 = np.array([-3, 0, 0]) * u.cm
pt2 = np.array([3, 10, 0]) * u.cm

print(f"Line from pt0 to pt1 intersects: {grid.vector_intersects(pt0, pt1)}")
print(f"Line from pt0 to pt2 intersects: {grid.vector_intersects(pt0, pt2)}")

Line from pt0 to pt1 intersects: True
Line from pt0 to pt2 intersects: False


## Interpolating Quantities¶

Grid objects contain several interpolator methods to evaluate quantites at positions between the grid vertices. These interpolators use the fact that all of the quantities are defined on the same grid to perform faster interpolations using a nearest-neighbor scheme. When an interpolation at a position is requested, the grid indices closest to that position are calculated. Then, the quantity arrays are evaluated at the interpolated indices. Using this method, many quantities can be interpolated at the same positions in almost the same amount of time as is required to interpolate one quantity.

Positions are provided to the interpolator as a u.Quantity array of shape [N,3] where N is the number of positions and [i,:] represents the x,y,z coordinates of the ith position:

[15]:

pos = np.array([[0.1, -0.3, 0], [0.5, 0.25, 0.8]]) * u.cm
print(f"Pos shape: {pos.shape}")
print(f"Position 1: {pos[0,:]}")
print(f"Position 2: {pos[1,:]}")

Pos shape: (2, 3)
Position 1: [ 0.1 -0.3  0. ] cm
Position 2: [0.5  0.25 0.8 ] cm


The simplest interpolator directly returns the nearest-neighbor values for each quantity. Positions that are out-of-bounds return an interpolated value of zero.

[16]:

Ex_vals = grid.nearest_neighbor_interpolator(pos, "E_x")
print(f"Ex at position 1: {Ex_vals[0]:.2f}")

Ex at position 1: 0.17 V / m


Multiple values can be interpolated at the same time by including additional keys as arguments. In this case, the interpolator returns a tuple of arrays as a result.

[17]:

Ex_vals, Ey_vals, Ez_vals = grid.nearest_neighbor_interpolator(pos, "E_x", "E_y", "E_z")
print(f"E at position 1: ({Ex_vals[0]:.2f}, {Ey_vals[0]:.2f}, {Ez_vals[0]:.2f})")

E at position 1: (0.17 V / m, 0.74 V / m, 0.89 V / m)


For a higher-order interpolation, some grids (such as the CartesianGrid subclass) also include a volume-weighted interpolator. This interpolator averages the values on the eight grid vertices surrounding the position (weighted by their distance). The syntax for using this interpolator is the same:

[18]:

Ex_vals, Ey_vals, Ez_vals = grid.volume_averaged_interpolator(pos, "E_x", "E_y", "E_z")
print(f"E at position 1: ({Ex_vals[0]:.2f}, {Ey_vals[0]:.2f}, {Ez_vals[0]:.2f})")

E at position 1: (0.28 V / m, 0.66 V / m, 0.69 V / m)


If repeated identical calls are being made to the same interpolator at different positions (for example, in a simulation loop), setting the persistent keyword to True will increase performance by not repeatedly re-loading the same quantity arrays from the grid. Setting this keyword for a single interpolation will not improve performace, and is not recommended (and is only done here for illustration).

[19]:

Ex_vals, Ey_vals, Ez_vals = grid.volume_averaged_interpolator(
pos, "E_x", "E_y", "E_z", persistent=True
)
print(f"E at position 1: ({Ex_vals[0]:.2f}, {Ey_vals[0]:.2f}, {Ez_vals[0]:.2f})")

E at position 1: (0.28 V / m, 0.66 V / m, 0.69 V / m)

[ ]: