Inflated geometric factors storage for affine elements

[inducer]

In #310 (comment), @kaushikcfd says:

I noticed a subtle and sub-optimal behavior that is worth noting. If we are using PyOpenCLActx as a setup actx and PytatoActx as an operational actx, then the frozen geometric factors (such as normals) would be stored at per-DOF even for affine geometry mappings. Such use-cases will bump memory footprint.

If true, this wouldn't just increase memory footprint, but also memory bandwidth requirements.

cc @majosm @matthiasdiener

[majosm]

I'm told that we don't use a separate array context for init anymore in Y3, so it looks like this isn't affecting us at the moment.

[inducer]

@kaushikcfd Could you elaborate on the circumstance under which you've observed this? Do you happen to know why this happens?

[kaushikcfd]

I did not observe this in any production-related drivers. But I was writing some tests for #310 and noticed this behavior and thought sharing would be essential so that others do not repeat the same mistake.

I have attached a minimal reproducer script for posterity.

import grudge
import pyopencl as cl
from meshmode.discretization.connection import FACE_RESTR_INTERIOR

from grudge.array_context import get_reasonable_array_context_class
from meshmode.mesh.generation import generate_regular_rect_mesh

ctx = cl.create_some_context()
cq = cl.CommandQueue(ctx)

ref_actx = get_reasonable_array_context_class(lazy=False, distributed=False)(cq)
actx = get_reasonable_array_context_class(lazy=True, distributed=False)(cq)

mesh = generate_regular_rect_mesh(a=(-0.5, -0.5),
                                  b=(0.5, 0.5),
                                  nelements_per_axis=(16, 16))

dcoll = grudge.make_discretization_collection(ref_actx, mesh, order=4)
nx, ny = actx.thaw(dcoll.normal(FACE_RESTR_INTERIOR))

assert nx[0].shape[-1] == 1  # fails

[inducer]

Thanks for providing the reproducer. AFAICT, this affects only the normals. The main issue is that here:

grudge/grudge/discretization.py

Line 770 in 0d501c0

return self._setup_actx.freeze(normal(self._setup_actx, self, dd))

the setup actx is unconditionally used to populate the normal cache, making this an interface bug more than a "bug bug". Perhaps the most reasonable response is to deprecate DiscretizationCollection.normal. grudge.geometry.normal takes an actx and does not have this issue (and it's more convenient to use because there's no freeze/thaw cycle needed).

[inducer]

#315