Container-over-container broadcasting, broadcasting rules about arraycontext HOT 5 CLOSED

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.

from arraycontext.

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?

from arraycontext.

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.

from arraycontext.

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.

from arraycontext.

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

from arraycontext.