Introduce DISCR_TAG_CONSTANT (e.g. for outputs of elementwise reductions)

[majosm]

dt_geometric_factors is returning different array shapes in lazy vs. eager. In eager it returns a volume-dd-shaped DOFArray, but in lazy it returns a DOFArray with 1 value per element. This discrepancy appears to originate in _apply_elementwise_reduction, where different code is executed based on the value of actx.supports_scalar_broadcasting. I can get the expected (volume-dd-shaped) result in both lazy and eager if I comment out the True case.

The True case looks like this:

    if actx.supports_nonscalar_broadcasting:
        return DOFArray(
            actx,
            data=tuple(
                getattr(actx.np, op_name)(vec_i, axis=1).reshape(-1, 1)
                for vec_i in vec
            )
        )

I suspect the reshape(-1, 1) should be changed to something else, but I'm not sure what numpy-like constructs are supported for the array types being operated on here. Any suggestions?

[inducer]

Everything you describe, at least for now, is by design. I'll try to describe the reasoning for this, and then we can see whether this is something that should change.

The intention is for the output of the routine to be treated as if it had a normal DOFArray shape. In the lazy case, broadcasting will automatically "make it so" once the first "normal" DOFArray is encountered in arithmetic. Until then, there is no real reason to redundantly represent what is just a single value per element.

Since the PyOpenCL Array types used with eager are not smart enough to broadcast, for them instead a redundantly represented array is returned.

If nothing else, I suspect this behavior should be documented better. It is the same for all elementwise reductions.

[majosm]

Makes sense. I think something extra might need to be done in the case where someone tries to project one of these elementwise DOFArrays though. I currently get an error when I try it:

meshmode.dof_array.InconsistentDOFArray: DOFArray group 0 array has unexpected shape. (observed: (1976, 1), expected: (1976, 10))

(I can just multiply by ones beforehand, but I'm not sure if that's optimal.)

[inducer]

I've spent some time thinking about this (typing different variants of this reply), and I actually don't think the current approach of trying to pass off a single-DOF-per-element array as a "full" DOFArray makes a great deal of sense. There's one situation where this works transparently (broadcasting), everywhere else it's trouble. If we wanted to honestly support this, we would need to litter an unbounded number of places (not metaphorically speaking!) with conditionals, which (IMO) makes this an obvious nonstarter.

I also don't want to go back to returning the same value plastered across an element; in addition to inefficiency (which lazy may or may not be able to remove), it's redundant representation which is not nice.

So my proposal is that we introduce another discretization tag, DISCR_TAG_CONSTANT (or so, name up for debate) that describes these arrays. There exists a discretization (P^0 or Q^0) for which they're valid after all! We could document that arithmetic between DISCR_TAG_CONSTANT and other discretization tags is allowed (even modal I think!), preserving the status quo in that regard. With that change, I think everything should more or less just work unchanged, as long as we pass the right DOFDescs, even the project you mention.

How does that sound?