# 04. Gradient of the misfit function#

A basic example how to use the `emg3d.simulations.Simulation.gradient` routine to compute the adjoint-state gradient of the misfit function. Here we just show its usage.

For this example we use the survey and data as obtained in the example 03. Simulation.

```import os
import pooch
import emg3d
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm, SymLogNorm

```

First we load the survey and accompanying data as obtained in the example 03. Simulation.

```fname = 'GemPy-II-survey-A.h5'
pooch.retrieve(
'https://raw.github.com/emsig/data/2021-05-21/emg3d/surveys/'+fname,
'5f2ed0b959a4f80f5378a071e6f729c6b7446898be7689ddc9bbd100a8f5bce7',
fname=fname,
path=data_path,
)

# Let's have a look
survey
```
```Data loaded from «/home/dtr/3-GitHub/emsig/emg3d-gallery/examples/download/GemPy-II-survey-A.h5»
[emg3d v0.17.1.dev22+ge23c468 (format 1.0) on 2021-04-14T22:24:03.229285].
```

### Survey «GemPy-II Survey A»

```<xarray.Dataset>
Dimensions:    (src: 3, rec: 45, freq: 2)
Coordinates:
* src        (src) <U6 'TxED-1' 'TxED-2' 'TxED-3'
* rec        (rec) <U7 'RxEP-01' 'RxEP-02' 'RxEP-03' ... 'RxEP-44' 'RxEP-45'
* freq       (freq) <U3 'f-1' 'f-2'
Data variables:
observed   (src, rec, freq) complex128 (-5.965512862515066e-14-5.52219071...
synthetic  (src, rec, freq) complex128 (-1.821164983623268e-14-5.42793941...
Attributes:
noise_floor:     1e-15
relative_error:  0.05```

We can see that the survey consists of three sources, 45 receivers, and two frequencies.

## Create an initial model#

To create an initial model we load the true model, but set all subsurface resistivities to 1 Ohm.m. So we are left with a homogeneous three-layer model air-seawater-subsurface, which includes the topography of the seafloor.

```# Load true model
fname = "GemPy-II.h5"
pooch.retrieve(
'https://raw.github.com/emsig/data/2021-05-21/emg3d/models/'+fname,
'ea8c23be80522d3ca8f36742c93758370df89188816f50cb4e1b2a6a3012d659',
fname=fname,
path=data_path,
)
grid = model.grid
```
```Data loaded from «/home/dtr/3-GitHub/emsig/emg3d-gallery/examples/download/GemPy-II.h5»
[emg3d v1.0.0rc3.dev5+g0cd9e09 (format 1.0) on 2021-05-21T18:40:16.721968].
```
```# Overwrite all subsurface resistivity values with 1.0
res = model.property_x
subsurface = (res > 0.5) & (res < 1000)
res[subsurface] = 1.0
model.property_x = res

# QC the initial model and the survey.
grid.plot_3d_slicer(model.property_x, xslice=12000, yslice=7500,
pcolor_opts={'norm': LogNorm(vmin=0.3, vmax=200)})

# Plot survey in figure above
fig = plt.gcf()
fig.suptitle('Initial resistivity model (Ohm.m)')
axs = fig.get_children()
src_coords = survey.source_coordinates()
axs[1].plot(rec_coords[0], rec_coords[1], 'bv')
axs[2].plot(rec_coords[0], rec_coords[2], 'bv')
axs[3].plot(rec_coords[2], rec_coords[1], 'bv')
axs[1].plot(src_coords[0], src_coords[1], 'r*')
axs[2].plot(src_coords[0], src_coords[2], 'r*')
axs[3].plot(src_coords[2], src_coords[1], 'r*')
plt.show()
```

## Options for automatic gridding#

```gridding_opts = {
'center': (src_coords[0][1], src_coords[1][1], -2200),
'properties': [0.3, 10, 1, 0.3],
'domain': (
[rec_coords[0].min()-100, rec_coords[0].max()+100],
[rec_coords[1].min()-100, rec_coords[1].max()+100],
[-5500, -2000]
),
'min_width_limits': (100, 100, 50),
'stretching': (None, None, [1.05, 1.5]),
'center_on_edge': False,
}
```

## Create the Simulation#

```simulation = emg3d.simulations.Simulation(
name="Initial Model",    # A name for this simulation
survey=survey,           # Our survey instance
model=model,             # The model
gridding='both',         # Src- and freq-dependent grids
max_workers=4,           # How many parallel jobs
# solver_opts=...,       # Any parameter to pass to emg3d.solve
gridding_opts=gridding_opts,
)

# Let's QC our Simulation instance
simulation
```

### Simulation «Initial Model»

• Survey «GemPy-II Survey A»: 3 sources; 45 receivers; 2 frequencies
• Model: resistivity; isotropic; 100 x 100 x 102 (1,020,000)
• Gridding: Frequency- and source-dependent grids; 192 x 48 x 48 (442,368)

```grad = simulation.gradient
```
```Compute efields            0/6  [00:00]
Compute efields #6         1/6  [01:21]
Compute efields #####      3/6  [01:21]
Compute efields ########3  5/6  [01:53]
Compute efields ########## 6/6  [01:53]

Back-propagate            0/6  [00:00]
Back-propagate #6         1/6  [01:28]
Back-propagate ########3  5/6  [01:59]
Back-propagate ########## 6/6  [01:59]
```

```# Set the gradient of air and water to NaN.
# This will eventually move directly into emgd3 (active and inactive cells).

grid.plot_3d_slicer(
pcolor_opts={'cmap': 'RdBu_r',
'norm': SymLogNorm(
linthresh=1e-2, base=10, vmin=-1e1, vmax=1e1)}
)

fig = plt.gcf()
axs = fig.get_children()
axs[1].plot(rec_coords[0], rec_coords[1], 'bv')
axs[2].plot(rec_coords[0], rec_coords[2], 'bv')
axs[3].plot(rec_coords[2], rec_coords[1], 'bv')
axs[1].plot(src_coords[0], src_coords[1], 'r*')
axs[2].plot(src_coords[0], src_coords[2], 'r*')
axs[3].plot(src_coords[2], src_coords[1], 'r*')
plt.show()
```
```emg3d.Report()
```
 Wed Aug 31 21:31:36 2022 CEST OS Linux CPU(s) 4 Machine x86_64 Architecture 64bit RAM 15.5 GiB Environment Python File system ext4 Python 3.9.12 | packaged by conda-forge | (main, Mar 24 2022, 23:22:55) [GCC 10.3.0] numpy 1.22.4 scipy 1.9.0 numba 0.55.2 emg3d 1.8.0 empymod 2.2.0 xarray 2022.6.0 discretize 0.8.2 h5py 3.7.0 matplotlib 3.4.3 tqdm 4.64.0 IPython 8.4.0 Intel(R) oneAPI Math Kernel Library Version 2022.0-Product Build 20211112 for Intel(R) 64 architecture applications

Total running time of the script: ( 4 minutes 3.304 seconds)

Estimated memory usage: 610 MB

Gallery generated by Sphinx-Gallery