Final timestep index not clearly defined when `save=False` in TimeData

[vincepandolfo]

This is the stencil generated for test_diffusion.py:

extern "C" int Operator(float *u_vec)
{
  float (*u)[100][100] = (float (*)[100][100]) u_vec;
  {
    int t0;
    int t1;
    for (int i3 = 0; i3<19; i3+=1)
    {
      {
        t0 = (i3)%(2);
        t1 = (t0 + 1)%(2);
      }
      {
        for (int i1 = 1; i1<99; i1++)
        {
          #pragma GCC ivdep
          for (int i2 = 1; i2<99; i2++)
          {
            u[t1][i1][i2] = 1.11022302462516e-16F*u[t0][i1][i2] + 2.5e-1F*u[t0][i1][i2 - 1] + 2.5e-1F*u[t0][i1][i2 + 1] + 2.5e-1F*u[t0][i1 - 1][i2] + 2.5e-1F*u[t0][i1 + 1][i2];
          }
        }
      }
    }
  }
  return 0;
}

The values of t0 and t1 alternate, so it's not clearly defined what index contains the final timestep. If the number of timesteps was even, the final result would be at index 1, while in this code the final result is at index 0

[navjotk]

Which timestep is last written to is a property and also the decision of the Operator that was run on it. I believe we should keep in mind that these data objects would be used for multiple operations and would exist between multiple Operator apply calls. Therefore the last time step written to should be seen as a piece of information the Operator controls and not the data.

[vincepandolfo]

There's no logic in the operator to do that as of right now though

[mloubout]

@navjotk , this data may be used for different operators. As we originally designed it, the PDE (operator) is a solver that get it's input (nt, src/rec pos, h,...) from the data and the model. The data should however not be aware of any time indexing in the C code, and the last time step is the last index in data. In case of cyclic index, the data need to know which index is which time for reuse.

[navjotk]

Yeah now I understand and appreciate this problem better. I think the way to go forward would really be for the data structure to be aware of these cyclic indices indeed. Overriding the array access and redirecting it to the right location is what I suggest. This is probably what you were originally thinking of too, @vincepandolfo

[mlange05]

OK, I think there is am interesting corner case here that we need to consider in the design of this: Suppose we have a 3-way time-buffered grid, and we do 100 time steps on it; we end up writing t + 1 last. Now if we were to apply the same operator again, it assumes we start from t, not t + 1. So what do we do? We could:

• Copy/shuffle data at the end/beginning of Operator.apply()
• Re-generate the code and adjust time indexing
• Keep meta-data about "last accessed time index", pipe this into the auto-generate C-layer and add an implicit offset to the cyclic time indexing?

I'm not quite sure what the solution to this is, so alternative suggestions are welcome.

[mloubout]

@fablup, didn't you fix this already?