Seismic example re-structuring

[mlange05]

The current seismic examples (Acoustic and TTI) both use explicit wrappers for Forward, Adjoint, Gradient and Born wrapper methods that create new operators for every kernel invocation. This prohibits re-use in tight loops, such as checkpoint-based RTM, because every application of the wrapper bypasses the operator-level caching and triggers a full DSE/DLE/codegen/compilation cascade. We should re-structure the examples to, at the very least, utilise pre-setup operators for the different invocations, or even allow a seamless argument injection into the pre-defined operator objects.

The seismic utility wrappers also previously instantiated all symbolic data object at invocation time, which causes costly memory-copy operations. Some of this has moved to a slightly higher level, but a more concrete model of data flow for the various types of grid data and the various desired use cases (adjoint test, gradient test, RTM, FWI, etc.) is required, in particular with focus on which data lives at which temporal granularity and who instantiates it (user-level, operator-wrapper).

[mlange05]

As an example let's consider the checkpointing case with all data created by the user (either explicitly or through seismic-specific abstractions/utility classes:

# Data allocation
model = Model(velocity, ...)
src = RickerSource(src_coordinates, time_series)
rec = Receivers(rec_coordinates)
u = TimeData(u, ...)
v = TimeData(v, ...)

# Operator instantiation, wrappers return type Operator/StencilKernel
op_fwd = ForwardOperator(model=model, u=u, source=src, receivers=rec, ...)
op_grad = GradientOperator(model=model, v=v, ...)

# Incremental forward loop
ckp = Checkpointer()  # Something clever that handles checkpointing/scheduling
nsteps = 0
while nsteps > 0:
   nsteps = ckp.how_many_timesteps_should_I_run()
   op_fwd(u=u, m=model.m, time_size=nsteps, ...)
   ckp.do_your_thing(u)

for all checkpoints:
    u = ckp.retrieve_forward_wavefield()
    op_grad(v=v, u=u, m=model.m, ...)

The key questions arising from this are

• Can we make op_fwd/op_grad pure Operator/StencilKernel objects?
• When/where do we instantiate the wavefields? Should these have an additional abstraction to capture ABCs, or can we hide the PML on the Model??
• When/where are sources/receivers instantiated and populated?

Note: The above is a rough draft and likely to change.

[mloubout]

• Can we make op_fwd/op_grad pure Operator/StencilKernel objects?

Not sure about this one, a lot of setup is still required.

• When/where do we instantiate the wavefields? Should these have an additional abstraction to capture ABCs, or can we hide the PML on the Model?

I think ABC should be abstracted out. The current version, as pointed out in #140 , as some problems. Also, for PML and ABC we have symbolic expression for it and should be special loops over the small portion of the domain concerned. I had a try at it, but Iteration seems to be broken with 3.0 so I need to rethink it. For the wavefield, this is the tricky question, specially with the Julia interface. I owever think we can discuss it in person with philipp once we get in London, but this shouldn't be required to be solved for the GMD paper.

• When/where are sources/receivers instantiated and populated?

This is usually done at the beginning, as this will never change during the inversion. However, with thousand of sources, it is likely that different sources will have different receiver number and/or positions.

[mlange05]

• Can we make op_fwd/op_grad pure Operator/StencilKernel objects?

Not sure about this one, a lot of setup is still required.

Can you please be more specific? What setup/which fields are required on a per-invocation basis, and which live beyond a single operator run?

I had a try at it, but Iteration seems to be broken with 3.0 so I need to rethink it.

Iteration objects are now used internally only for 3.0 while loop iterations are simply implied by the dimensions used. Can you provide a concrete example of what you tried, and how it failed (possibly in the appropriate issue for ABCs #140)?

For the wavefield, this is the tricky question, specially with the Julia interface. I owever think we can discuss it in person with philipp once we get in London, but this shouldn't be required to be solved for the GMD paper.

I agree this question does not necessarily need to be solved for GMD. Nevertheless I would like to have some idea of what the requirements are (and what rationale we have for them), since this problem still influences a lot of other design decisions. In a similar vain, being clear about the kind of domain-specific utilities and the type of examples code we want to provide is important for further development, and will allow us to make an informed decision on what level of clean-up we can reasonably do for GMD, given a much changed operator API.

[mloubout]

fixed with #243 ?

[mlange05]

Yes, everything not covered in #243 is veing tracked in other issues.