Helper functions¶

Estimating stable time-steps¶

Helper functions for estimating stable time steps for RKDG methods.

Characteristic lengthscales¶

grudge.dt_utils.characteristic_lengthscales(actx: ArrayContext, dcoll: DiscretizationCollection, dd: DOFDesc | None = None) → DOFArray[source]¶

Computes the characteristic length scale \(h_{\text{loc}}\) at each node. The characteristic length scale is mainly useful for estimating the stable time step size. E.g. for a hyperbolic system, an estimate of the stable time step can be estimated as \(h_{\text{loc}} / c\), where \(c\) is the characteristic wave speed. The estimate is obtained using the following formula:

\[h_{\text{loc}} = \operatorname{min}\left(\Delta r_i\right) r_D\]

where \(\operatorname{min}\left(\Delta r_i\right)\) is the minimum node distance on the reference cell (see dt_non_geometric_factors()), and \(r_D\) is the inradius of the cell (see dt_geometric_factors()).

Returns:: a DOFArray containing a characteristic lengthscale for each element, at each nodal location.

Note

While a prediction of stability is only meaningful in relation to a given time integrator with a known stability region, the lengthscale provided here is not intended to be specific to any one time integrator, though the stability region of standard four-stage, fourth-order Runge-Kutta methods has been used as a guide. Any concrete time integrator will likely require scaling of the values returned by this routine.

Non-geometric quantities¶

grudge.dt_utils.dt_non_geometric_factors(dcoll: DiscretizationCollection, dd: DOFDesc | None = None) → Sequence[float][source]¶

Computes the non-geometric scale factors following [Hesthaven_2008], section 6.4, for each element group in the dd discretization:

\[c_{ng} = \operatorname{min}\left( \Delta r_i \right),\]

where \(\Delta r_i\) denotes the distance between two distinct nodal points on the reference element.

Parameters:: dd – a DOFDesc, or a value convertible to one. Defaults to the base volume discretization if not provided.
Returns:: a list of float values denoting the minimum node distance on the reference element for each group.

Mesh size utilities¶

grudge.dt_utils.dt_geometric_factors(dcoll: DiscretizationCollection, dd: DOFDesc | None = None) → DOFArray[source]¶

Computes a geometric scaling factor for each cell following [Hesthaven_2008], section 6.4, defined as the inradius (radius of an inscribed circle/sphere).

Specifically, the inradius for each element is computed using the following formula from [Shewchuk_2002], Table 1, for simplicial cells (triangles/tetrahedra):

\[r_D = \frac{d V}{\sum_{i=1}^{N_{faces}} F_i},\]

where \(d\) is the topological dimension, \(V\) is the cell volume, and \(F_i\) are the areas of each face of the cell.

Parameters:: dd – a DOFDesc, or a value convertible to one. Defaults to the base volume discretization if not provided.
Returns:: a frozen DOFArray containing the geometric factors for each cell and at each nodal location.

Returns a (maximum) characteristic length based on the volume of the elements. This length may not be representative if the elements have very high aspect ratios.

Parameters:

dim – an integer denoting topological dimension. If None, the spatial dimension specified by grudge.DiscretizationCollection.dim is used.
dd – a DOFDesc, or a value convertible to one. Defaults to the base volume discretization if not provided.

Returns:

a scalar denoting the maximum characteristic length.

Returns a (minimum) characteristic length based on the volume of the elements. This length may not be representative if the elements have very high aspect ratios.

Parameters:

dim – an integer denoting topological dimension. If None, the spatial dimension specified by grudge.DiscretizationCollection.dim is used.
dd – a DOFDesc, or a value convertible to one. Defaults to the base volume discretization if not provided.

Returns:

a scalar denoting the minimum characteristic length.

Array contexts¶

class grudge.array_context.PyOpenCLArrayContext(queue: pyopencl.CommandQueue, allocator: pyopencl.tools.AllocatorBase | None = None, wait_event_queue_length: int | None = None, force_device_scalars: bool = True)[source]¶: Inherits from meshmode.array_context.PyOpenCLArrayContext. Extends it to understand grudge-specific transform metadata. (Of which there isn’t any, for now.)

class grudge.array_context.PytatoPyOpenCLArrayContext(queue, allocator=None, *, compile_trace_callback: Callable[[Any, str, Any], None] | None = None)[source]¶: Inherits from meshmode.array_context.PytatoPyOpenCLArrayContext. Extends it to understand grudge-specific transform metadata. (Of which there isn’t any, for now.)

class grudge.array_context.MPIBasedArrayContext[source]¶

class grudge.array_context.MPIPyOpenCLArrayContext(mpi_communicator, queue: pyopencl.CommandQueue, *, allocator: pyopencl.tools.AllocatorBase | None = None, wait_event_queue_length: int | None = None, force_device_scalars: bool = True)[source]¶

An array context for using distributed computation with pyopencl eager evaluation.

__init__(*args, **kwargs)[source]¶: Initialize self. See help(type(self)) for accurate signature.

class grudge.array_context.MPINumpyArrayContext(mpi_communicator)[source]¶

An array context for using distributed computation with numpy eager evaluation.

__init__(*args, **kwargs)[source]¶: Initialize self. See help(type(self)) for accurate signature.

class grudge.array_context.MPIPytatoArrayContext[source]¶

grudge.array_context.get_reasonable_array_context_class(lazy: bool = True, distributed: bool = True, fusion: bool | None = None, numpy: bool = False) → type[ArrayContext][source]¶: Returns a reasonable ArrayContext currently supported given the constraints of lazy, distributed, and numpy.

Miscellaneous tools¶

grudge.tools.build_jacobian(actx: ArrayContext, f: Callable[[ArrayOrContainerT], ArrayOrContainerT], base_state: ArrayOrContainerT, stepsize: float) → ndarray[source]¶

Returns a Jacobian matrix of f determined by a one-sided finite difference approximation with stepsize.

Parameters:

f – a callable with a single argument, any array or arraycontext.ArrayContainer supported by actx.
base_state – The point at which the Jacobian is to be calculated. May be any array or arraycontext.ArrayContainer supported by actx.

Returns:

a two-dimensional numpy.ndarray.

grudge.tools.map_subarrays(f: Callable[[Any], Any], in_shape: tuple[int, ...], out_shape: tuple[int, ...], ary: Any, *, return_nested: bool = False) → Any[source]¶

Apply a function f to subarrays of shape in_shape of an numpy.ndarray, typically (but not necessarily) of dtype object. Return an numpy.ndarray with the corresponding subarrays replaced by the return values of f, and with the shape adapted to reflect out_shape.

Similar to numpy.vectorize.

Example 1: given a function f that maps arrays of shape (2, 2) to scalars and an input array ary of shape (3, 2, 2), the call:

map_subarrays(f, (2, 2), (), ary)

will produce an array of shape (3,) containing the results of calling f on the 3 subarrays of shape (2, 2) in ary.

Example 2: given a function f that maps arrays of shape (2,) to arrays of shape (2, 2) and an input array ary of shape (3, 2), the call:

map_subarrays(f, (2,), (2, 2), ary)

will produce an array of shape (3, 2, 2) containing the results of calling f on the 3 subarrays of shape (2,) in ary. The call:

map_subarrays(f, (2,), (2, 2), ary, return_nested=True)

will instead produce an array of shape (3,) with each entry containing an array of shape (2, 2).

Parameters:

f – the function to be called.
in_shape – the shape of any inputs to f.
out_shape – the shape of the result of calling f on an array of shape in_shape.
ary – a numpy.ndarray instance.
return_nested – if out_shape is nontrivial, this flag indicates whether to return a nested array (containing one entry for each application of f), or to return a single, higher-dimensional array.

Returns:

an array with the subarrays of shape in_shape replaced with subarrays of shape out_shape resulting from the application of f.

grudge.tools.rec_map_subarrays(f: Callable[[Any], Any], in_shape: tuple[int, ...], out_shape: tuple[int, ...], ary: ArrayOrContainer, *, scalar_cls: type | tuple[type] | None = None, return_nested: bool = False) → ArrayOrContainer[source]¶

Like map_subarrays(), but with support for arraycontext.ArrayContainers.

Parameters:: scalar_cls – An array container of this type will be considered a scalar and arrays of it will be passed to f without further destructuring.