Tutorial 1 – Grids and fields
This first part of the tutorial introduces grids and fields. Grids define the space on which fields are discretized. These fields then describe the state of the partial differential equations.
# import packages
import sys
import numpy as np
sys.path.append("../..") # add the pde package to the python path
import pde
We start by defining the space and its discretization. The simplest variant is a rectangular geometry with unit-sized cells, which we call a UnitGrid. Cells with unequal sizes are supported by the more complex CartesianGrid. Finally, we support a small set of curvilinear (yet still orthogonal) coordinate systems: PolarSymGrid on a disk, SphericalSymGrid in a ball, and CylindricalSymGrid in an axisymmetric cylinder. All these grids have in common that they assume that fields do not depend on the angular variables, i.e., they enforce symmetry.
grid = pde.UnitGrid([32, 32])
grid.plot(title=f"Area={grid.volume}")
grid = pde.CartesianGrid([[0, 2], [0, 1]], [16, 32])
grid.plot(title=f"Area={grid.volume}")
grid = pde.PolarSymGrid(5, 16)
grid.plot(title=f"Area={grid.volume:.3g}")
PolarSymGrid and SphericalSymGrid also support annuli. The two radii can be specified by a supplying a tuple to the radius argument.
grid = pde.PolarSymGrid((3, 5), 8)
grid.plot(title=f"Area={grid.volume:.3g}")
Scalar fields represent scalar quantities that depend on position. In this package, scalar fields are represented by their values at the support points of the discretized grids discussed above. Consequently, one first has to construct a grid and pass it to the ScalarField class to construct a scalar field.
grid = pde.CartesianGrid([[0, 4 * np.pi], [0, 2 * np.pi]], [128, 32])
field = pde.ScalarField(grid, data=1)
field.data
Scalar fields support a range of mathematical operations and can thus roughly used like bare numpy arrays. The actual underlying data of the field is accessed by its .data attribute.
field += 4
field.data
Instead of looking at the raw data, fields will be most often visualized in plots. A convenient plotting method is the .plot() method.
grid = pde.CartesianGrid([[0, 4 * np.pi], [0, 2 * np.pi]], [128, 32])
data = np.arange(128 * 32).reshape(grid.shape)
field = pde.ScalarField(grid, data=data)
field.plot(colorbar=True);
Values at individual points can be determined by interpolation
field.interpolate([2.1, 0.3])
There are a range of methods to initialize scalar fields. For instance, various random fields can be initiated using the ScalarField.random_*() methods.
field = pde.ScalarField.random_normal(grid)
field.plot(title=f"Average={field.average:.3f}", colorbar=True);
field = pde.ScalarField.random_colored(grid, exponent=-4, label="Random field")
field.plot();
Fields can be written to files using the optional h5py package to write the data in the Hierarchical Data Format (HDF). Note that the underlying grid is also stored in the file and recreated transparently when the file is read using the ScalarField.from_file() method.
field.to_file("random_field.hdf")
field_loaded = pde.ScalarField.from_file("random_field.hdf")
field_loaded.plot();
Fields can be further analyzed by slicing or projeting them along given axes. Slicing uses interpolation to calculate field values at positions that might not lie on the original grid. In contrast, projecting integrates over the given axes, which are thereby removed.
slice_x = field.slice({"y": np.pi})
slice_x.plot(title="Slice field along line $y=π$");
project_x = field.project("y")
project_x.plot(title="Projection by integration along $y$");
Finally, fields can also be created from mathematical expressions that are parsed using sympy. Note that this general method is unsafe to process user-supplied data, since it uses the exec function.
field = pde.ScalarField.from_expression(grid, "sin(x) * cos(y)")
field.plot(title=f"Average={field.average:f}", colorbar=True);
Various differential operators can be applied to fields. These differential operators are just-in-time compiled using numba. Consequently, their first evaluation can take quite long, since code needs to be analyzed and compiled. However, each subsequent evaluation will be very fast.
Note that differential operators require boundary conditions to be well-defined. Boundary conditions can be specified in a variety of formats, as shown by some examples below. More information on the various formats can be found in the documentation.
grid = pde.CartesianGrid([[0, 4 * np.pi], [0, 2 * np.pi]], [128, 32])
field = pde.ScalarField.from_expression(grid, "sin(x) * cos(y)")
laplace_dir = field.laplace({"value": 0})
laplace_dir.plot(
title="Laplacian of field with Dirichlet boundary conditions", colorbar=True
);
laplace_neu = field.laplace({"derivative": 0})
laplace_neu.plot(
title="Laplacian of field with Neumann boundary conditions", colorbar=True
);
Note that the Laplace operator results in large values at the boundary where the supplied boundary condition is incompatible with the actual field data. In this case, this can be fixed by imposing Dirichlet conditions along the $x$-axis and Neumann conditions along the $y$-axis:
laplace_mix = field.laplace([{"value": 0}, {"derivative": 0}])
laplace_mix.plot(
title="Laplacian of field with mixed boundary conditions", colorbar=True
);
Alternatively, the Laplace operator can be evaluate assuming periodic boundary conditions. Note that periodic boundary conditions are a special case of boundary conditions, since they also affect how distances are measured in the defined space. Consequently, periodic boundary conditions need to be already declared on the grid instance.
grid_per = pde.CartesianGrid([[0, 4 * np.pi], [0, 2 * np.pi]], [128, 32], periodic=True)
field_per = pde.ScalarField.from_expression(grid_per, "sin(x) * cos(y)")
laplace_per = field_per.laplace("periodic")
laplace_per.plot(title="Laplacian of periodic field", colorbar=True);
Periodic and non-periodic axes can also be mixed as in the example below. In this case, the $x$-axis is periodic and thus requires periodic boundary conditions. Conversely, the $y$-axis is non-periodic and any other boundary condition can be specified. The most generic one is a Neumann condition of vanishing derivative. For convenience, we also define natural boundary conditions, which indicate periodic conditions for periodic axes and Neumann conditions otherwise.
grid_mixed = pde.CartesianGrid(
[[0, 4 * np.pi], [0, 2 * np.pi]], [128, 32], periodic=[True, False]
)
field_mixed = pde.ScalarField.from_expression(grid_mixed, "sin(x) * cos(y)")
laplace_mixed = field_mixed.laplace(["periodic", {"derivative": 0}])
# laplace_mixed = field_mixed.laplace('natural')
laplace_mixed.plot(title="Laplacian of mixed field", colorbar=True);
Vector and tensor fields appear often in physical equations, e.g., due to applied differential operators. They offer almost all the methods already introduced for scalar fields above.
grid_per = pde.CartesianGrid([[0, 4 * np.pi], [0, 2 * np.pi]], [128, 32], periodic=True)
vector_field = pde.VectorField.from_expression(grid_per, ["1", "cos(x)"])
vector_field.plot();
Vector fields can also originate from differential operators
field_per = pde.ScalarField.from_expression(grid_per, "sin(x) * cos(y)")
field_grad = field_per.gradient("natural")
field_grad
field_grad.average
res = field_grad.plot(method="quiver", title="Quiver plot of the gradient field");
field_grad.plot(method="streamplot", title="Stream plot of the gradient field");
gradient_norm = field_grad.to_scalar("norm")
gradient_norm.plot(title="Norm of gradient of field", colorbar=True);
Individual components can be extracted by subscripting.
comp_x = field_grad[0]
comp_y = field_grad[1]
comp_x.plot(title="$x$-component of gradient", colorbar=True)
comp_y.plot(title="$y$-component of gradient", colorbar=True);
This can be helpful to calculate quantities component-wise.
gradient_expl_norm = (comp_x**2 + comp_y**2) ** 0.5
np.allclose(gradient_expl_norm.data, gradient_norm.data)
Scalar and vector fields can be interpreted as fields of rank 0 and 1, respectively. The package also supports rank-2 tensors, which are represented by Tensor2Field. Such a tensorial field can also originate from differential operators, e.g., as a Hessian.
field_hess = field_grad.gradient("natural", label="Hessian of field")
field_hess.attributes
field_hess.plot();
Vector and tensor fields can also be used in dot products, either via the .dot() method or using the @ notation denoting matrix products in in python.
scalar_field = field_grad @ field_hess @ field_grad
scalar_field.plot(title="Gradient . Hessian . Gradient", colorbar=True);
scalar_field.slice({"y": np.pi}).plot();
Many partial differential equations combine several fields, which can be represented as a collection.
grid_pol = pde.PolarSymGrid([2, 7], 32)
scalar_field1 = pde.ScalarField.from_expression(grid_pol, "r**2", label="Increasing")
scalar_field2 = pde.ScalarField.from_expression(grid_pol, "1/r", label="Decreasing")
collection = pde.FieldCollection([scalar_field1, scalar_field2])
collection.attributes
collection.plot();
collection.plot("image");
(collection[0] + collection[1]).plot("image");
grid = pde.UnitGrid([96, 64])
sf = pde.ScalarField.random_colored(grid, exponent=-3)
vf = pde.VectorField.random_harmonic(grid, modes=1)
fc = pde.FieldCollection([sf, vf])
fc.plot();