Why does the symbol $dt$ disappear in C/C++ code?

[jiwei08]

I wanna use the symbol dt given by model.grid.stepping_dim.spacing in the index of the matrix. But the symbol dt disappears in the C/C++ code when I check it using print(op). So what should I do to fix it? I am also confused with model.grid.stepping_dim.spacing and model.grid.time_dim.spacing. It seems like they both represent dt.

Any suggestions are appreciated!

[jiwei08]

I have fixed this bug by introducing a devito.Constant index and define index by devito.Eq(). But I would still like to know why this bug arise.

[FabioLuporini]

are you using Operator(..., subs=my_own_subs, ...) where my_own_subs includes a numerical substitution for dt?

[jiwei08]

I change to use index=my_own_subs where index is a devito.Constant of a int16 dtype. This indeed works. But I also want to know why using the subs directly doesn’t work.

[FabioLuporini]

can you write an MFE?

see here: https://github.com/devitocodes/devito/blob/master/FAQ.md#how-do-i-find-the-source-of-my-bug-quickly-and-get-support

[jiwei08]

OK. Here is an MFE.

from devito import Grid, Operator, TimeFunction, Eq, Constant
from examples.seismic import Model
import numpy as np

model = Model(
    origin=(0, 0, 0),
    spacing=(0.1, 0.1, 0.1),
    shape=(11, 11, 11),
    space_order=2,
    vp=1500.0,
)
t = model.grid.time_dim
dt = model.grid.time_dim.spacing
time_func = TimeFunction(name="time_func", shape=(100,), dimensions=(t,))
a = TimeFunction(name="a", shape=(100,), dimensions=(t,))
# dt disappears in the C-code when using t/dt as subs directly, but it works through Constant index.
op = Operator(Eq(a, time_func[t / dt]))
# index = Constant(name="index", dtype=np.int16)
# op = Operator([Eq(index, t / dt)] + [Eq(a, time_func[index])])
print(op)

[mloubout]

Hi

This is actually the expected behavior. When used as an index the spacing is replaced by 1 at lowering so that FD indices are correct, for example f[x + dx] become as indexed f[x +1] so in your case it does lower to what is intended.

Separate question, what is the intended behavior for t/dt ? t is the dimension i.e the index not the physical time, the physical time is t*dt so not sure if it's the confusion here

[jiwei08]

Thank you very much for helping me!!! For the question of t/dt. I just want to show this failure case which does not have any physical meaning. In fact, I want to calculate the delay time of the emitted source wave in my project so t-(d[0]*x*dx+d[1]*y*dy+d[2]*z*dz)/vp/dt are used where the vector d is the direction of emitted wave.

Thanks again for helping me a lot!

[jiwei08]

I also have another question. Does model.grid.time_dim.spacing and model.grid.stepping_dim.spacing represent the same thing? If not, what‘s the diffenences between them?

[mloubout]

They are the same thing as time_dim is the parent of stepping_dim

[jiwei08]

Emmmm, so why only dt is replaced by 1 and dx,dy,dz are replaced by h_x,h_y,h_z. Is it because waveform is a TimeFunction variable?

[mloubout]

They are only replaced by 1 when used as pure indices not as expressions in equations. For example dt is used during propagation as a dt**2 * ... and is replaced by its value in the wave stencil.

In your case they are within a floor function so they are considered an expression not an index but will double check there is no odd corner case we are missing

[jiwei08]

But dt and dx,dy,dz are all within a floor function. Why only dt is replaced by 1?