SubdomainSet issues

[rhodrin]

The current implementation of SubdomainSet's has the following issues that need to be fixed:

1. ~Many cases do not work with Devito loop optimization and require operators to be created with dle='noop'. See test_subdomains.py for examples.~

2. Currently implicit equations are only added to the expression list 'once' for each Subdomain set. When complex dependencies exist between equation on different subdomains this will lead to incorrect c-code.

3. ~MPI is not currently supported.~

4. The _add_implicit method located in Operator needs refactoring. This is occurring too early in the compilation process leading to certain properties of equations being 'discovered' twice.

5. This functionality needs a notebook tutorial.

[rhodrin]

As of 8/8/2019 issue 1 is no longer present in tests conducted thus far.

[rhodrin]

3. fixed by #1216 .

[FabioLuporini]

Some thoughts below

Instead of generating

for n
  a_m = ...[n];
  a_M = ...[n];
  for a = a_m to a_M ...
     ...

we could treat SubDomainSet as what it actually is -- a set of blocks of, generally, different size. This is a generalisation of loop blocking. Or, IOW, loop blocking is a special case of SubDomainSet where the whole domain is covered by disjoint blocks (those at the domain boundaries potentially being smaller than the inner ones, leading to the so called "remainder loops"). So instead of creating one single Dimension(n), we could create Dimension(nx), Dimension(ny), ..., that is one "block" Dimension per space dimension. These could be -- just like in the case of loop blocking -- IncrDimension(nx, ...). This would lead to generating

for nx = nx_m to nx_M, 1:
  for ny = ny_m to ny_M, 1:
    x_m = array_x_m[nx]
    x_M = array_x_M[nx]
    for x = x_m to x_M, 1:
       ...

This would preserve the block-wise iteration semantics, which I would be able to exploit for optimization in the later stages of the compilation pipeline. For example, with a "flatten" n Dimension it's much more difficult to reconstruct the information that if I hoist some time-invariant code outside of n and time then I should not infer that I'm avoiding n*time*... operations, but rather just time.

Anyway, food for thought...