Non-time dependent Functions must be passed as Time Function

[speglich]

Hi! We are working on PR #1406 and we are facing some issues with temporaries variables that are non-time dependent, but must be passed as Time Function. When they are declared as Function a seg fault error occurs.

import numpy as np

from devito import (Eq, Operator, Function, VectorFunction, VectorTimeFunction, TimeFunction, NODE,
                    div, grad, norm)

from examples.seismic import Receiver, PointSource, demo_model, setup_geometry, seismic_args

def MFE():

    model = demo_model(preset='layers-viscoacoustic', space_order=4, shape=(50, 50), nbl=40,
                    dtype=np.float32, spacing=(20.0, 20.0))

    geometry = setup_geometry(model, tn=1000.0)

    s = model.grid.stepping_dim.spacing
    f0 = geometry._f0
    vp = model.vp
    b = model.b
    qp = model.qp
    damp = model.damp

    w0 = 2. * np.pi * f0
    rho = 1. / b
    eta = (vp * vp) / (w0 * qp)
    bm = rho * (vp * vp)

    p = TimeFunction(name="p", grid=model.grid,time_order=2, space_order=4,
                     staggered=NODE)

    l = VectorFunction(name="l", grid=model.grid, space_order=4)

    # Works fine
    #h = TimeFunction(name="h", grid=model.grid, time_order=2, space_order=4,
    #                 staggered=NODE)

    #w = VectorTimeFunction(name="w", grid=model.grid, time_order=2,
    #                       space_order=4)

    # When H and W are declared as Function/VectorFunction a seg fault error occurs
    h = Function(name="h", grid=model.grid, space_order=4, staggered=NODE)

    w = VectorFunction(name="w", grid=model.grid, space_order=4)

    u_h = Eq(h, (p - p.backward) / s)

    u_l = Eq(l, b * grad(p))

    u_w = Eq(w, b * grad(h))

    pde_p = 2. * p - damp * p.backward + s * s * bm * div(l) + \
        s * s * eta * rho * div(w)

    u_p = Eq(p.forward, damp * pde_p)

    eqn = [u_h, u_l, u_w, u_p]

    op = Operator(eqn)

    op.apply(time_M=10)

    return 0

[mloubout]

I don't think you need h. With it, the equation u_w is fully independent of time and will split the operator in two different time loops which is not good (but still shouldn't segault will try to figure out why). Runs fine with u_w = Eq(w, b / s * grad(p - p.backward)) to more generally u_w = Eq(w, b * grad(p.dt(fd_order=1, side=left)))

[nogueirapeterson]

@mloubout Do you have any idea why when we pre-compute (using an auxiliary variable), the halo problem disappears?

[nogueirapeterson]

to more generally u_w = Eq(w, b * grad(p.dt(fd_order=1, side=left)))

Is this more general way recommended when using solve to create the stencil?

[mloubout]

Is this more general way recommended when using solve to create the stencil?

Both are the same, the more general one s more readable for users.

[Leitevmd]

Hey guys,

I think we have two topics related to this issue:

1. Similar formulations are resulting in different time-iteration loop arrangements (including more than one time-iteration loop within the operator).
2. Segfaults and compilation errors occur due to more than one time-iteration loop.

Let's get first to the second topic, which is actually producing these errors:

Errors are occurring because substitutions of SteppingDimensions by ModuloDimensions in _lower_stepping_dims(iet) (devito/ir/iet/scheduler.py) are being done in a misleading manner (I suppose): All expressions in iet have their indices replaced with uindices available in the first time-iteration loop (eg. t+2 -> t0+2). That way, when we get to the second time-iteration loop, indices that were supposed to be replaced (eg. t+2 -> t2) are not found because they have already been replaced (eg. t+2 -> t0+2) and so they remain that way in the code.

This is wrong because the second time-iteration loop may not even have that t0 uindice (which causes compilation error) or if it has t0 it can attempt to access wrong index (segfault), since indexes are miscalculated (eg. t0+2 != t2, because t0 = (time)%(3), t2 = (time+2)%(3), but t0+2 could reach 4). That could even not produce any errors at all, but still mess up the numerical result.

I think we can solve this issue changing _lower_stepping_dims(iet) by making uindices substitutions to occur only at child expressions of each iteration at a time, thus building a mapper which further transforms the entire iet (See code here). This way we have no more problems of declaration missing or segfaults. @FabioLuporini what do you think about it?

[Leitevmd]

Almost forgot.. here is a simpler MFE:

from devito import (Eq, Operator, Grid, Function, TimeFunction)

grid = Grid(shape=(10,10))

f = Function(name="f", grid=grid, space_order=4)
g = Function(name="g", grid=grid, space_order=4)
h = TimeFunction(name="h", grid=grid, space_order=4, time_order=2)

eqn = []
eqn.append( Eq(f, h + 1) )
eqn.append( Eq(g, f) )
eqn.append( Eq(h.forward, h + g + 1) )

op = Operator(eqn)

op.apply(time_M=10)

Comments:

• Removing the intermediary equation above, produces only one time-iteration loop. Eg., replace second and third eq lines by the following: eqn.append( Eq(h.forward, h + f + 1) )

• Also, turning g into a TimeFunction, again produces only one time-iteration loop.

[FabioLuporini]

@Leitevmd

This is the list of equations you've produced:

[Eq(f(x, y), h(t, x, y) + 1),
 Eq(g(x, y), f(x, y)),
 Eq(h(t + dt, x, y), g(x, y) + h(t, x, y) + 1)]

The second one has no concept of time, ie, the time dimension doesn't appear in it. What should the compiler do:

• either assume "time" is implicit in the second equation because at least one of the surrounding coupled equations (through f and g) are over time
• or assume the equation really is time-independent and therebefore be scheduled outside the time loop?

there is no simple answer to this question, though someone could prefer point one above over point two. But this would present other issues in other contexts....

so how do we work around this?

Add:

time = grid.time_dim
x, y = grid.dimensions

and decorate the second equation with implicit dimensions (see https://www.devitoproject.org/devito/equation.html#devito.types.equation.Eq)

Eq(g, f, implicit_dims=(time, x, y))

Now you're explicitly telling devito to schedule this expression in the time-x-y iteration space (with loops in this given order)

This should fix the problem

Said that, we should fix the case without implicit_dims anyway, because devito is generating invalid code, and you seem to have an understanding of why. In particular, we get h[t0 + 1][x + 4][y + 4] in the generated code and this is gonna cause an out-of-bounds access as soon as t0==2, thus a potential segfault. This should be fixed, I don't know if you wanna move this to a separate issue or what, but I'd be extremely glad if we found a work around

we could also a short 06_implicit_dims notebook here: https://github.com/devitocodes/devito/tree/master/examples/userapi explaining what the implicit dimensions are for (@speglich , @Leitevmd , ...?)

one final thing: @mloubout could you open an issue where you clearly exaplain why we need these user-level temporaries for now? eg what needs to be fixed with staggering and the differential operators such that this horrible thing becomes unnecessary? probably to be flagged even as bug-py .

[Leitevmd]

That was really clarifying! Thanks a lot! I could work on this short tutorial.

I'm moving the second topic (case without implicit_dims) to another issue then. I just need to organize some other information I have about it, first.

[FabioLuporini]

Prodding @mloubout -- did you read my previous comment ?

[mloubout]

Yes sorry, slipped through. The not time dependent was mentioned earlier and why need the implicit dim indeed.

WIll open issue.