Different code generation when using subdomains

[ZoeLeibowitz]

Sub expression elimination behaves strangely on subdomains if you replace a dimension with a Symbol.

This code:

from devito import Grid, TimeFunction, Eq, solve, Operator
from devito.types import Symbol
grid = Grid(shape=(5, 5), extent=(1., 1.))

f = TimeFunction(name='f', grid=grid, space_order=2)
g = TimeFunction(name='g', grid=grid, space_order=2)

mapper = {grid.stepping_dim: Symbol('tau')}

eq = Eq(f.dt, g.laplace)
stencil = solve(eq, f.forward)
update = Eq(f.forward, stencil, subdomain=grid.interior).subs(mapper)
op = Operator([update])

produces the following loop:

  START(section0)
  for (int x = x_ltkn0 + x_m; x <= x_M - x_rtkn0; x += 1)
  {
    for (int y = y_ltkn0 + y_m; y <= y_M - y_rtkn0; y += 1)
    {
      f[tau + 1][x + 2][y + 2] = dt*(g[tau][x + 1][y + 2]/((h_x*h_x)) - 2.0F*g[tau][x + 2][y + 2]/(h_x*h_x) + g[tau][x + 3][y + 2]/((h_x*h_x)) + g[tau][x + 2][y + 1]/((h_y*h_y)) - 2.0F*g[tau][x + 2][y + 2]/(h_y*h_y) + g[tau][x + 2][y + 3]/((h_y*h_y)) + f[tau][x + 2][y + 2]/dt);
    }
  }
  STOP(section0,timers)

The issues are:

• None of the Temp expressions are generated
• Double brackets are generated
• SIMD alignment doesn't happen

If you remove subdomain=grid.interior then the generated code is as expected:

  float r0 = 1.0F/dt;
  float r1 = 1.0F/(h_x*h_x);
  float r2 = 1.0F/(h_y*h_y);

  START(section0)
  for (int x = x_m; x <= x_M; x += 1)
  {
    #pragma omp simd aligned(f,g:32)
    for (int y = y_m; y <= y_M; y += 1)
    {
      f[tau + 1][x + 2][y + 2] = dt*(r0*f[tau][x + 2][y + 2] + r1*g[tau][x + 1][y + 2] + r1*g[tau][x + 3][y + 2] + r2*g[tau][x + 2][y + 1] + r2*g[tau][x + 2][y + 3] + (-2.0F)*(r1*g[tau][x + 2][y + 2] + r2*g[tau][x + 2][y + 2]));
    }
  }
  STOP(section0,timers)