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

`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:

```python
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?

Start at _apply_elementwise_reduction and compare the branches selected by actx.supports_nonscalar_broadcasting, reproducing the dt_geometric_factors result in lazy and eager modes. Check the supported array operations used by the relevant array types. Done means both modes return the expected volume-dd-shaped DOFArray, with coverage for the discrepancy.

Written by the indexing model from the issue text.