Python Bindings
###############
Fields
******
The µGrid library handles field quantities, i.e. scalar, vectors or tensors,
that vary in space. It supports only a uniform discretization of space. In the
language of µGrid, we call the discrete coordinates **pixels**. Each pixel
is associated with a physical position in space. µGrid supports fields in
two and three dimensional Cartesian grids. Note that a common name for a pixel
in three dimensions is a **voxel**, but we refer to it as a pixel throughout
this documentation.
Each pixel
can be subdivided into logical elements, for example into a number of
support points for numerical quadrature. These subdivisions are called
**sub-points**. Note that while µGrid understands the physical location of
each *pixel*, the *sub-points* are logical subdivisions and not associated
with any position within the pixel. Each *sub-point* carries the field quantity,
which can be a scalar, vector or any type of tensor. The field quantity is
called the **component**.
Multidimensional arrays
***********************
Each µGrid field has a representation in a multidimensional array. This
representation is used when accessing a field with Python via `numpy `_
or when writing the field to a file. The default representation has the shape
.. code-block:: Python
(components, sub-points, pixels)
As a concrete example, a second-rank tensor (for example the deformation
gradient) living on two quadrature points in two dimensions with a spatial
discretization of 11 x 12 grid points would have the following shape:
.. code-block:: Python
(3, 3, 2, 11, 12)
Note that the components are omitted if there is only a single component (i.e.
a scalar value), but the sub-points are always included even if there is only
a single sub-point.
Throughout the code, we call the field quantities **entries**. For example,
`nb_pixels` refers to the number of pixels while `nb_sub_pts` refers to the
number of sub-points per pixel. `nb_entries` is then the product of the two.
The above field has another multidimensional array representation that folds
the sub-point into the last dimension of the components. For the above example,
this means a multidimensional array of shape:
.. code-block:: Python
(3, 6, 11, 12)
Depending on the numerical use case, it can be useful to work with either
representation.
Field collections
*****************
In µGrid, fields are grouped into field collections. A field collection knows
about the spatial discretization of the fields, but each field can differ in number of
sub-points and components.
There are two kinds of field collections:
* A **global** collection that groups fields that have values at all pixels.
* A **local** collection that groups fields that have values only at a subset of pixels.
An example would be the elastic constants of a material that only exists at
certain locations in the domain.
The following example shows how to initialize a field collection and create scalar fields:
.. literalinclude:: ../../examples/field_collection.py
:language: python
The first argument to the constructor of `FieldCollection` is the spatial dimension.
Fields are then registered with the field accessor methods of the field collection,
e.g. `real_field` in the example above. Fields are *named*, and the name needs to be
unique. Accessing a field of the same name yield the same field object.
Note the `FieldCollection` additionally has register methods, e.g. `register_real_field`.
The different to `real_field` method is that the explicit registration of the field fails if
it already exists, while the accessor method `real_field` registers it if it does not exist
but returns the existing field if it does. You can also query a field with `get_field`,
which will raise an exception if the field does not exist.
Components
**********
In the above example, we registered a scalar field, which has no components. Vector or
tensor-valued field can be defined by specifying either simply a number of components or
the shape of the components explicitly. The following example shows how to create a
tensor-valued field that contains 2 x 2 matrices:
.. literalinclude:: ../../examples/components.py
:language: python
Note that there is a difference between a scalar field and a field with a single component.
When a single component is specified, this is reflected in the shape of the field:
.. literalinclude:: ../../examples/scalar.py
:language: python
Sub-points
**********
A pixel can be subdivided into multiple sub-points, each of which holds a value (scalar,
vector or tensor) of the field. Examples for sub-points are elements or quadrature points
in a the finite element method.
Sub-points are named, e.g. common names would be `element` for subdivision into elements or
`quad` for quadrature points. The name of the subdivision must be specified when the field
is created.
The following example initializes a three-dimension grid with
a sub-division of each pixel into five elements:
.. literalinclude:: ../../examples/sub_points.py
:language: python
The example demonstrates two ways of accessing the field. The convenience accessor `p` (that
we also used in the above examples) and the new accessor `s`. Both yield a numpy array that is
a view on the underlying data, but with different shapes. The `s` accessor has the shape
`(components, sub-points, pixels)` and exposes the sub-points explicitly. The `p` accessor folds
the sub-points into the last dimension of the components.
Convolutions
************
A common operation in numerical analysis is the convolution of a field with a short-ranged stencil. µGrid has utility
functions that make this convolution simple. A convolution turns a field defined on nodal points into a field defined
on quadrature points. This means field of shape
.. code-block:: Python
(components, nodal-points, pixels)
is turned into a field of shape
.. code-block:: Python
(components, quadrature-points, pixels)
where the same convolution is applied to each component. Note that the convolution operation itself may have multiple
components (which we will refer to as operators in the following), meaning that the output field may have the shape
.. code-block:: Python
(components, operators, quadrature-points, pixels)
In mathematical notation, the convolution operation can be written as
.. math::
f_{c,o,q,p} = \sum_{n} \sum_{k} s_{o,q,n,k} g_{c,n,p+k}
where :math:`f` is the output field, :math:`g` is the input field, :math:`s` is the stencil, :math:`o` are the operators, :math:`c` are the components,
:math:`q` are the quadrature points, :math:`n` are the nodal points, :math:`p` are the pixels and :math:`k` runs over the stencil width
(possibly in 2 or 3 dimensions).
The convolution operation also has a transpose operation that turns a field defined on quadrature points into a field defined on nodal points.
In mathematical notation, the transpose convolution is written as
.. math::
g_{c,n,p} = \sum_{o} \sum_{q} \sum_{k} s_{o,q,n,k} f_{c,o,q,p-k}
with the same stencil :math:`s` as above.
As an example, we consider the gradient of a two dimensional field with a single nodal point.
If we use two linear finite elements per pixel, we can represent the derivative by two quadrature points per pixel.
The gradient itself has two components, one for each direction of the gradient, which means we need two operators
to represent it. The stencil then has the shape
.. code-block:: Python
(operators, quadrature-points, nodal-points, pixels)
or in this specific case
.. code-block:: Python
(2, 2, 1, 2, 2)
because we are only considering a two-dimensional problem where there are two operators of the gradient.
Note that the stencil can be of a lower dimensional shape, in which case missing dimensions are assumed to
be unity. It is assumed that first `nodel-points`, then `quadrature-points` and finally `operators` are
missing. This means we can represent a simple forward-differences gradient operator in two dimensions as
.. code-block:: Python
stencil = np.array([
[[-1, 1], [0, 0]],
[[-1, 0], [1, 0]]
])
op = muGrid.GenericLinearOperator(2, stencil)
.. literalinclude:: ../../examples/gradient.py
:language: python
Array accessors
***************
Each field provides four accessors that return views into the underlying field data
as multidimensional arrays. For fields on the host (CPU), these return numpy arrays.
For fields on the GPU, these return CuPy arrays. See the :doc:`GPU` chapter
for details on GPU computing.
The accessors differ in two ways:
1. **Layout**: SubPt (``s``) vs Pixel (``p``)
2. **Ghost regions**: Excluded (``s``, ``p``) vs Included (``sg``, ``pg``)
SubPt vs Pixel layout
---------------------
The **SubPt layout** (``s``, ``sg``) exposes the sub-points as an explicit dimension:
.. code-block:: Python
displacement_field.s # shape: (components, sub_pts, *spatial_dims)
The **Pixel layout** (``p``, ``pg``) folds the sub-points into the last
dimension of the components:
.. code-block:: Python
displacement_field.p # shape: (components * sub_pts, *spatial_dims)
The SubPt layout is useful when you need to operate on sub-points separately,
while the Pixel layout is convenient for operations that treat all sub-point
values uniformly.
With vs without ghost regions
-----------------------------
When using domain decomposition with ghost regions (for stencil operations),
the accessors ending in ``g`` include the ghost cells:
.. code-block:: Python
# Excluding ghost regions (interior domain only)
displacement_field.s # SubPt layout, no ghosts
displacement_field.p # Pixel layout, no ghosts
# Including ghost regions (full local domain with ghosts)
displacement_field.sg # SubPt layout, with ghosts
displacement_field.pg # Pixel layout, with ghosts
For fields without ghost regions, ``s`` and ``sg`` (and ``p`` and ``pg``)
return views of the same shape.
The ghost-excluding accessors (``s``, ``p``) are typically used for computations
on the interior domain, while ghost-including accessors (``sg``, ``pg``) are
used for operations that need to access or modify the full buffer including
ghost cells.
The entries of the field occur as the first indices in the multidimensional
because a numerical code is typically vectorized of the spatial domain, i.e.
we carry out the same operation on each pixel but not on each component.
This means for the above displacement field, we can simply get the components
of the field from
.. code-block:: Python
ux, uy, uz = displacement_field.s
Each of the variable `ux`, `uy`, and `uz` is a three-dimensional array with shape
`(2, 11, 12)`. *numpy*'s
`broadcasting rules `_
make it simple to vectorize over pixels,
for example normalizing the displacement field could look like:
.. code-block:: Python
displacement_field.s /= np.sqrt(ux**2 + uy**2 + uz**2)
Note that the default storage order of the field is column-major, which means the
field components are stored next to each other in memory.
I/O
***
Fields can be written to disk in the `NetCDF `_ format.
µGrid uses `Unidata NetCDF `_ when build
with just serial capabilities and `PnetCDF `_ when build
with MPI enabled.
I/O is handled by the `FileIONetCDF` class. The following example shows how to write all fields
from a field collection to disk:
.. literalinclude:: ../../examples/io.py
:language: python
The file has the following structure (output of `ncdump -h`):
.. code-block::
netcdf example {
dimensions:
frame = UNLIMITED ; // (1 currently)
tensor_dim__strain-0 = 3 ;
tensor_dim__strain-1 = 3 ;
subpt__element-5 = 5 ;
nx = 11 ;
ny = 12 ;
nz = 13 ;
variables:
double strain(frame, tensor_dim__strain-0, tensor_dim__strain-1, subpt__element-5, nx, ny, nz) ;
strain:unit = "no unit provided" ;
// global attributes:
:creation_date = "01-06-2024 (d-m-Y)" ;
:creation_time = "23:02:06 (H:M:S)" ;
:last_modified_date = "01-06-2024 (d-m-Y)" ;
:last_modified_time = "23:02:06 (H:M:S)" ;
:muGrid_version_info = "µGrid version: 0.90.1+40-g5bc73b30\n",
"" ;
:muGrid_git_hash = "5bc73b305881ef837ce6598568d41eb3d0307c41" ;
:muGrid_description = "0.90.1+40-g5bc73b30" ;
:muGrid_git_branch_is_dirty = "false" ;
}