Typical FWI workflow

To perform full-waveform inversion, the following actions must be performed:

  1. Set some observed data (we did this in the previous cell)

  2. Forward simulate all sources

  3. Calculate the misfit

  4. Calculate the adjoint sources

  5. Empty the kernels (from any previous calculation)

  6. Adjoint simulate all sources

  7. Transform the sensitivity kernels (in Lamé’s paramters) to velocity

In Notebook 2, these actions are described. If you simply want a class that combines everything, use the following code, given you have a fdModel instance stored in solver:

class _fwi:

    last_model = numpy.empty_like(solver.get_model_vector())
    misfit = None
    g = None
    smoothing = None

    def __init__(self):
        pass

    def misfit(self, m):

        if numpy.allclose(m, self.last_model):
            return self.misfit
        else:
            _ = self.grad(m)
            return self.misfit

    def grad(self, m):

        if numpy.allclose(m, self.last_model):
            return self.g
        else:
            solver.set_model_vector(m)

            # Simulate forward
            for i_shot in range(solver.n_shots):
                solver.forward_simulate(i_shot, omp_threads_override=6)

            # Calculate misfit and adjoint sources
            solver.calculate_l2_misfit()
            solver.calculate_l2_adjoint_sources()
            self.misfit = solver.misfit

            # Simulate adjoint
            solver.reset_kernels()
            for i_shot in range(solver.n_shots):
                solver.adjoint_simulate(i_shot, omp_threads_override=6)
            solver.map_kernels_to_velocity()

            g = solver.get_gradient_vector()

            if self.smoothing is not None:
                g = numpy.hstack(
                    [
                        gaussian_filter(
                            igs.reshape((60, 180)), sigma=self.smoothing
                        ).flatten()
                        for igs in numpy.split(g, 3)
                    ]
                )

            self.g = g

            self.last_model = m

            return self.g