Wavefield indexing causing compilation error when saving the wavefield

[kgullikson88]

I am trying to use hybrid boundary conditions while saving the forward wavefield. I found the reference notebook for these boundary conditions here, and updated in two ways:

1. I used the simpler A1 conditions to rule out some of the more complicated equations
2. The wavefield definition includes the keyword argument save=time_range.num
3. Running the operation is via op(time=nt-1,dt=dt0)

However, this now produces compilation errors as shown below:

/tmp/devito-jitcache-uid885437/3be37c299ad0452c96d1f5507e8bab91d2ff5ec7.c: In function ‘Kernel’:
/tmp/devito-jitcache-uid885437/3be37c299ad0452c96d1f5507e8bab91d2ff5ec7.c:216:23: error: ‘posx’ undeclared (first use in this function)
  216 |           if (rrecx + posx >= x_m - 1 && rrecy + posy >= y_m - 1 && rrecx + posx <= x_M + 1 && rrecy + posy <= y_M + 1)
      |                       ^~~~
/tmp/devito-jitcache-uid885437/3be37c299ad0452c96d1f5507e8bab91d2ff5ec7.c:216:23: note: each undeclared identifier is reported only once for each function it appears in
/tmp/devito-jitcache-uid885437/3be37c299ad0452c96d1f5507e8bab91d2ff5ec7.c:216:50: error: ‘posy’ undeclared (first use in this function)
  216 |           if (rrecx + posx >= x_m - 1 && rrecy + posy >= y_m - 1 && rrecx + posx <= x_M + 1 && rrecy + posy <= y_M + 1)
      |                                                  ^~~~
/tmp/devito-jitcache-uid885437/3be37c299ad0452c96d1f5507e8bab91d2ff5ec7.c:218:13: error: ‘sum’ undeclared (first use in this function)
  218 |             sum += (rrecx*px + (1 - rrecx)*(1 - px))*(rrecy*py + (1 - rrecy)*(1 - py))*u[time][rrecx + posx + 2][rrecy + posy + 2];
      |             ^~~
/tmp/devito-jitcache-uid885437/3be37c299ad0452c96d1f5507e8bab91d2ff5ec7.c:218:27: error: ‘px’ undeclared (first use in this function)
  218 |             sum += (rrecx*px + (1 - rrecx)*(1 - px))*(rrecy*py + (1 - rrecy)*(1 - py))*u[time][rrecx + posx + 2][rrecy + posy + 2];
      |                           ^~
/tmp/devito-jitcache-uid885437/3be37c299ad0452c96d1f5507e8bab91d2ff5ec7.c:218:61: error: ‘py’ undeclared (first use in this function)
  218 |             sum += (rrecx*px + (1 - rrecx)*(1 - px))*(rrecy*py + (1 - rrecy)*(1 - py))*u[time][rrecx + posx + 2][rrecy + posy + 2];
      |                                                             ^~
FAILED compiler invocation: gcc -march=native -O3 -g -fPIC -Wall -std=c99 -Wno-unused-result -Wno-unused-variable -Wno-unused-but-set-variable -ffast-math -shared -fopenmp /tmp/devito-jitcache-uid885437/3be37c299ad0452c96d1f5507e8bab91d2ff5ec7.c -lm -o /tmp/devito-jitcache-uid885437/3be37c299ad0452c96d1f5507e8bab91d2ff5ec7.so

I tracked it down to the free surface boundary conditions, defined below:

bc  = [Eq(u[t+1,x,0], u[t+1,x,1])]

If I exclude those from the operator definition, it works. Is indexing like this not allowed when saving the wavefield? Is there a different way to implement the corners?

I've also attached the full notebook I used, for reference/reproducibility.

habc_with_save (1).ipynb.zip

[mloubout]

You are indexing u with the wrong time dimension. You create u as save=nt so its time dimension is time (grid.time_dim) not t (grid.stepping_dim). By using the wrong time dimension it is breaking the time loop due to mismatch.

If you replace your bc by bc = [Eq(u[time+1,x,0], u[time+1,x,1])] then it works fine

[kgullikson88]

Huh, that does work but now I get a similar error if I do not save. What is the conceptual difference between grid.time_dim and grid.stepping_dim? Do I just need to know if I am saving the wavefield to know which thing to use?

[mloubout]

You just need to always use the dimension of u which you can easily do by using some of the API tools such as using u.forward._subs(y, 0) rather than u[time+1,x, 0] as it will always be the correct dimension, or getting the dimension from u i.e ut = u.dimensions[0] then do u[ut+1, ...]

[kgullikson88]

Ah, that works. Thanks for the help and the prompt reply!