5.1. Arrays and Array Containers in mirgecom

mirgecom is quite flexible in terms of what arrays it can use; in fact, no aspect of the code prescribes the use of a specific type of array. Instead, we use a arraycontext.ArrayContext to determine the types of arrays on which the computation should take place. At its simplest, this can be a arraycontext.PyOpenCLArrayContext, which uses pyopencl.array.Arrays that are eagerly evaluated. “Eager evaluation” means that when you type a + b with two arrays a and b, then the result is computed as soon as the expression is evaluated.

5.1.1. Lazy/deferred evaluation

For efficiency, it is sometimes better to defer the evaluation of the expression, to have an opportunity to apply some optimizations. We call this lazy or deferred evaluation; it is a common trick in efficiency-minded software (e.g. Google’s Tensorflow, or expression templates in C++). This is realized in mirgecom via an ArrayContext based on pytato (not merged to main) at the time of this writing.

5.1.2. Frozen and thawed arrays

All array contexts have a notion of thawed and frozen arrays that is important to understand, see Freezing and thawing for more details.

You might be wondering whether a given function or method will return frozen or thawed data. In general, the documentation should state that, however since a fair bit of functionality predates these concepts, we are still catching up in terms of updating the documentation. (Help welcome!)

In the meantime, these rules of thumb should cover most cases:

  • If you did not provide an array context to the function (explicitly or implicitly via an input array), you will receive frozen data.

  • Any data that is cached/long-lived/”at rest” is going to be frozen.

  • If the data is the result of a memoized function, then it will likely be frozen.

To demonstrate the effect of this, first we need some setup:

>>> import pyopencl as cl
>>> from arraycontext import PyOpenCLArrayContext, thaw
>>> ctx = cl.create_some_context()
>>> queue = cl.CommandQueue(ctx)
>>> actx = PyOpenCLArrayContext(queue)
>>> from meshmode.mesh.generation import generate_regular_rect_mesh
>>> mesh = generate_regular_rect_mesh(a=(0, 0), b=(1, 1), nelements_per_axis=(10, 10))
>>> from grudge import DiscretizationCollection
>>> dcoll = DiscretizationCollection(actx, mesh, order=5)

Most quantities that are maintained by the discretization will be frozen. For example, if one wanted to grab the nodes of the mesh or normals of a named surface in the mesh:

>>> from grudge.dof_desc import DOFDesc
>>> nodes = thaw(dcoll.nodes(), actx)
>>> dd = DOFDesc("all_faces")
>>> nhat = thaw(dcoll.normal(dd), actx)

What can go wrong? Attempts to operate on frozen data will yield errors similar to the following:

>>> dcoll.nodes() * 5
Traceback (most recent call last):
 ...
ValueError: array containers with frozen arrays cannot be operated upon

Fortunately, recovering from this is straightforward:

>>> nodes = thaw(dcoll.nodes(), actx)
>>> result = nodes * 5

5.1.3. Array Containers

Arrays in mirgecom live in (somewhat) deeply nested data structures that are ArrayContainers. Array containers typically support arithmetic and can be passed to most methods in ArrayContext that take arrays: they simply get applied to all arrays in the container. The same goes for many of the discretization-focused functions in grudge and meshmode.

For example, think of the solver state for Euler’s equations of gas dynamics: