devitocodes / devito

DSL and compiler framework for automated finite-differences and stencil computation
http://www.devitoproject.org
MIT License
566 stars 229 forks source link

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

Closed vincepandolfo closed 5 years ago

vincepandolfo commented 8 years ago

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 commented 8 years ago

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 commented 8 years ago

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

mloubout commented 8 years ago

@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 commented 8 years ago

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 commented 8 years ago

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:

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

mloubout commented 6 years ago

@fablup, didn't you fix this already?

navjotk commented 5 years ago

I think this is almost not relevant any more.

FabioLuporini commented 5 years ago

wdym ?

FabioLuporini commented 5 years ago

@navjotk

ggorman commented 5 years ago

This appears to have been resolved.