def METHOD_NAME(vp_in):
    # Create symbols to hold the gradient
    grad = Function(name="grad", grid=model.grid)
    objective = 0.
    for i in range(nshots):
        # Create placeholders for the data residual and data
        residual = Receiver(name='residual', grid=model.grid,
                            time_range=geometry.time_axis,
                            coordinates=geometry.rec_positions)
        d_obs = Receiver(name='d_obs', grid=model.grid,
                         time_range=geometry.time_axis,
                         coordinates=geometry.rec_positions)
        d_syn = Receiver(name='d_syn', grid=model.grid,
                         time_range=geometry.time_axis,
                         coordinates=geometry.rec_positions)
        # Update source location
        solver.geometry.src_positions[0, :] = source_locations[i, :]
        # Generate synthetic data from true model
        solver.forward(vp=model.vp, rec=d_obs)
        # Compute smooth data and full forward wavefield u0
        _, u0, _ = solver.forward(vp=vp_in, save=True, rec=d_syn)
        # Compute gradient from data residual and update objective function
        residual = compute_residual(residual, d_obs, d_syn)
        objective += .5*norm(residual)**2
        solver.jacobian_adjoint(rec=residual, u=u0, vp=vp_in, grad=grad)
    return objective, grad