Container-over-container broadcasting, broadcasting rules

[inducer]

As pointed out by @mtcam in https://github.com/illinois-ceesd/mirgecom/pull/355#discussion_r637478702:

E TypeError: unsupported operand type(s) for *: 'DOFArray' and 'ConservedVars'

(Here, DOFArrays clearly want to be inside ConservedVars, so there must be a rule to specify who is the "outer" structure.)

The biggest question here is which type should broadcast over which other type, and how the rules are to be communicated.

Some proposals:

• Just give the outer container's with_container_arithmetic a bcast_container_types=(type1, type2) argument.
• Define a _container_broadcast_priority: the higher the "outermost". All containers broadcast over all other containers (or at least ones that have this enabled in with_container_arithmetic. Kind of like numpy's __array_priority__.
• Have containers have a method may_broadcast_over_container_type(other_type).

Some thoughts:

• This is defining a partial order of who broadcasts over who. How do we ensure that this stays consistent? (An easy way to get an inconsistent situation is to tell two types that they should each be the "outer" type with respect to each other. That's easy to catch, but scenarios with three or more "edges" in the cycle are also possible and harder to catch.) Also, maybe it's not necessary to catch this statically...? Joys of multiple dispatch. What does Julia do?
• Technically, object arrays are part of this hierarchy. bcast_obj_array here, unfortunately, casually defines them as "outer". I think we'll have to deprecate that argument, barely a week after it was introduced...
• Of the three proposals above, I think I like bcast_container_types the best. The priority thing feels restrictive and ill-defined. may_broadcast_over_container_type is probably slow.
• Maybe we'll need ways to specify both "inner" and "outer" types. DOFArray has no reason to know about ConservedVars (so it shouldn't have to know who its outer types are). Meanwhile, if we let object arrays particpate in the hierarchy, then many types will want them as their outermost type, but we can't ask object arrays who their inner types are.

Thoughts?

cc @alexfikl @thomasgibson

[thomasgibson]

I will spend some time reading this more carefully, but an initial thought comes to mind. If we want to utilize multiple dispatch to help with this, it might be worth checking out: https://pypi.org/project/multimethod/

The only concerns I have are related to performance. I'm not sure how much heavy-duty compilation (python interpreter) needs to occur. We should carefully examine the performance impacts if we decide to go that route.

[alexfikl]

Just to make sure I understand what's meant to be broadcasted there. So in

flux_avg @ normal - 0.5*lam*(cv_tpair.ext - cv_tpair.int)

the normal is an object array of DOFArray, cv_tpair is a (TracePair of) ConservedVars, where the mass and energy are DOFArray and momentum is an object array, and flux_avg is a ConservedVars, where the mass and energy fluxes are object arrays of DOFArray and the momentum flux is a 2D object array of DOFArray. Is that right?

Then, we'd like to broadcast in the following cases

• DOFArray * object array (already in there), where DOFArray is the inner
• ConservedVars * DOFArray (for variables), where DOFArray is the inner
• ConservedVars * object array (for fluxes), where the object array is the inner Anything else?

Assuming most of that is true, I also prefer the bcast_container_types idea. I don't really see how this can work in a more general setting with something like a priority between container subtrees (?), so having each container say what it can broadcast over and how to do it seems to be the sanest choice.

The inner/outer thing should just work for arithmetic operations that define an __rmul__ (and friends), right? So that leaves the comparisons and bit shifts, which already don't broadcast very well. Do we want to support those?

[alexfikl]

Also, Julia has quite a bit of machinery to customize broadcasting for various types https://docs.julialang.org/en/v1/manual/interfaces/#man-interfaces-broadcasting Can't say I've understood all the edge cases they're trying to cover with it, but one nice thing is that they have that .op notation to explicitly request broadcasting.

[MTCam]

Just to make sure I understand what's meant to be broadcasted there. So in (...)

If what I say here is confusing or is just reiterating what you've already said, then please ignore, but this is also for me, so if I am misguided here about any of these pieces somebody please set me straight!

We are using the ConservedVars object to store quantities for each conservation equation. Most of the time, it is just the fluid (or simulation) "state" - specifically for us, a scalar quantity (1 DOFArray) for each of mass, energy, but an ndim-vector (object array of DOFArrays) for momentum, and an nspecies-vector for species_mass. In all there are neqn = ndim + 2 + nspecies conservation equations.

Sometimes, such as in this particular case, we also store ndim-vector quantities for each conservation equation, e.g. the flux vectors. We actually have both situations in the current expression we are considering. When we store vector quantities, the mass component (of ConservedVars) is the "mass flux" an ndim-vector, and momentum is then the ndim X ndim array with each row corresponding to the momentum component flux vector for each of the ndim components of momentum, and so on.

normal is an ndim-vector and we need flux @ normal (where flux is a ConservedVars storing ndim-vector quantities for each of the neqn conservation equations) to look at flux as like an neqn x ndim array.

lam is a DOFArray and we need lam * [q], where [q] is a ConservedVars object storing the simulation state, to have a "distributive" behavior wrt the components (i.e. conservation eqn) of the ConservedVars.

[inducer]

15 contains a proposal along these lines. It'd be great if you all could take a look.