devito creates new time blocks when using conditional dimensions as implicit

[deckerla]

import devito

grid = devito.Grid(shape=(5,5))
time = grid.time_dim
ct = devito.ConditionalDimension(name="ct", parent=time, factor=2)

u = devito.TimeFunction(name="u", grid=grid, time_order=1)
f = devito.Function(name="f", grid=grid, space_order=0)
g = devito.Function(name="g", grid=grid, space_order=0)
h = devito.Function(name="h", grid=grid, space_order=0)

# this simulates a type of imaging condition.  
# f*g should execute every snapshotted time step "ct"


op = devito.Operator([  devito.Eq(u.forward, u + 1.),
                        devito.Eq(f, u.forward, implicit_dims=(ct,)), 
                        devito.Eq(g, u.forward, implicit_dims=(ct,)), 
                        devito.Eq(h, h+f*g,     implicit_dims=(ct,)), 
                      ])


op.apply(time_M=5)

print(h.data)

This generates the following ccode:

#define _POSIX_C_SOURCE 200809L
#define START_TIMER(S) struct timeval start_ ## S , end_ ## S ; gettimeofday(&start_ ## S , NULL);
#define STOP_TIMER(S,T) gettimeofday(&end_ ## S, NULL); T->S += (double)(end_ ## S .tv_sec-start_ ## S.tv_sec)+(double)(end_ ## S .tv_usec-start_ ## S .tv_usec)/1000000;

#include "stdlib.h"
#include "math.h"
#include "sys/time.h"
#include "xmmintrin.h"
#include "pmmintrin.h"
#include "omp.h"

struct dataobj
{
  void *restrict data;
  unsigned long * size;
  unsigned long * npsize;
  unsigned long * dsize;
  int * hsize;
  int * hofs;
  int * oofs;
  void * dmap;
} ;

struct profiler
{
  double section0;
  double section1;
  double section2;
} ;

extern "C" int Kernel(struct dataobj *restrict f_vec, struct dataobj *restrict g_vec, struct dataobj *restrict h_vec, struct dataobj *restrict u_vec, const int time_M, const int time_m, const int x_M, const int x_m, const int y_M, const int y_m, const int nthreads, struct profiler * timers);


int Kernel(struct dataobj *restrict f_vec, struct dataobj *restrict g_vec, struct dataobj *restrict h_vec, struct dataobj *restrict u_vec, const int time_M, const int time_m, const int x_M, const int x_m, const int y_M, const int y_m, const int nthreads, struct profiler * timers)
{
  float (*restrict f)[f_vec->size[1]] __attribute__ ((aligned (64))) = (float (*)[f_vec->size[1]]) f_vec->data;
  float (*restrict g)[g_vec->size[1]] __attribute__ ((aligned (64))) = (float (*)[g_vec->size[1]]) g_vec->data;
  float (*restrict h)[h_vec->size[1]] __attribute__ ((aligned (64))) = (float (*)[h_vec->size[1]]) h_vec->data;
  float (*restrict u)[u_vec->size[1]][u_vec->size[2]] __attribute__ ((aligned (64))) = (float (*)[u_vec->size[1]][u_vec->size[2]]) u_vec->data;

  /* Flush denormal numbers to zero in hardware */
  _MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON);
  _MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);

  for (int time = time_m, t0 = (time)%(2), t1 = (time + 1)%(2); time <= time_M; time += 1, t0 = (time)%(2), t1 = (time + 1)%(2))
  {
    /* Begin section0 */
    START_TIMER(section0)
    #pragma omp parallel num_threads(nthreads)
    {
      #pragma omp for schedule(static,1)
      for (int x = x_m; x <= x_M; x += 1)
      {
        #pragma omp simd aligned(u:64)
        for (int y = y_m; y <= y_M; y += 1)
        {
          u[t1][x + 1][y + 1] = u[t0][x + 1][y + 1] + 1.0F;
        }
      }
    }
    STOP_TIMER(section0,timers)
    /* End section0 */
    if ((time)%(2) == 0)
    {
      /* Begin section1 */
      START_TIMER(section1)
      #pragma omp parallel num_threads(nthreads)
      {
        #pragma omp for schedule(static,1)
        for (int x = x_m; x <= x_M; x += 1)
        {
          #pragma omp simd aligned(f,g,u:64)
          for (int y = y_m; y <= y_M; y += 1)
          {
            f[x][y] = u[t1][x + 1][y + 1];
            g[x][y] = u[t1][x + 1][y + 1];
          }
        }
      }
      STOP_TIMER(section1,timers)
      /* End section1 */
    }
  }

  /* Begin section2 */
  START_TIMER(section2)
  #pragma omp parallel num_threads(nthreads)
  {
    #pragma omp for schedule(static,1)
    for (int x = x_m; x <= x_M; x += 1)
    {
      for (int time = time_m; time <= time_M; time += 1)
      {
        if ((time)%(2) == 0)
        {
          #pragma omp simd aligned(f,g,h:64)
          for (int y = y_m; y <= y_M; y += 1)
          {
            h[x][y] = f[x][y]*g[x][y] + h[x][y];
          }
        }
      }
    }
  }
  STOP_TIMER(section2,timers)
  /* End section2 */

  return 0;
}

As you can see, a new time block has been added after "propagation" of the time function has completed.

[FabioLuporini]

@mloubout fixed by your PR?

[mloubout]

Fixed in #2128 , closing