ORB does not allow arbitrary decomposition axis

[aliemen]

If I understand it correctly, the ORB algorithm only decomposes in all dimensions (or none). In OPALX we sometimes have bunches that are (in their rest frame) highly stretched in z-direction, which is why we would want to be able to decompose e.g. only in z direction.

The idea is to generalize binaryRepartition(...) by passing an additional const std::array<bool, Dim>& allowedAxes, where true indicates which axis to bisect. The old binaryRepartition signature will remain available as an overload. "Legacy" behavior is preserved by implicitly setting allowedAxes = {true, true, true} where necessary.

[aaadelmann]

Highly stretched mean that hr_m[2] is different i.e. larger. I fail to see how that changes the decomposition.

[aaadelmann]

The ORB can decompose in 1D, slab 2D or brick 3D

[srikrrish]

Yeah I am a bit confused as well. If you only set parallel decomposition along the z-direction (i.e. serial along other two which leads to slab as Andreas said) and perform ORB on this I think it should automatically work. But I am not sure if we have tested this though.

[aaadelmann]

But maybe I am confusing this with parallelisation.

[aliemen]

@aaadelmann regarding your first comment: true, just changing hr_m probably doesn't change much. But in these case we could e.g. have something like a (64, 64, 1024) grid. And then it would make sense to decompose only in z, right? If I understand correctly, that's also what's usually done in old OPAL (I think...).

Here is how I read the code: FieldLayout::initialize has this isParallel, and uses it for the initial partitioning. But FieldLayout::updateLayout(domains) then overwrites it with whatever the ORB algorithm computed. The current ORB always cuts the currently longest domain axis at each recursive step, meaning ALL axis are considered (and it ignores whatever isParallel in the field says). My addition keeps this behavior, but now ignores axis that are set to "false".

Old OPAL does more or less the same (reuse the PARFFTX/Y/Z from the solver for initial decomposition, runtime load balancing and solver decomposition).

[aaadelmann]

@aliemen did you check https://amas.pages.psi.ch/ETH/cse/ORBMichael.pdf ?

[aliemen]

@aaadelmann no, I didn't read it before 😅. But don't think the thesis mentions anything about selecting the axis during the ORB. In Algorithm 1 it says

cutAxis <- longest axis of doms[it]

...which I think is exactly what I described is happening in our ORB class: we always split along the longest axis. You can see the actual code implementation here:

return std::distance(dom.begin(), std::max_element(dom.begin(), dom.end(),
                   [&](const Index& a, const Index& b) {
                       return a.length() < b.length();
                   }));

My addition to that part is something like the following (in that reduction loop; a bit simplified):

if (!allowedAxes[d]) {
    continue;
}

Currently there is nothing in IPPL that would prevent ORB from choosing non-desired axes.

[aaadelmann]

Ok I see. My take here is we use first what we have and make sure the implemented domain decomposition is working in OPALX on large problems. I agree a full 3D domain decomposition is desirable but not something we need to consider now. We have much more pressing issues at hand. So in this respect this is a enhancement, as you have labeled it!

I guess this remains me to clearly make a list of things we have to attack next!