.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "content/tutorials/07-fdem/plot_fwd_2_fem_cyl.py"
.. LINE NUMBERS ARE GIVEN BELOW.
.. only:: html
.. note::
:class: sphx-glr-download-link-note
:ref:`Go to the end `
to download the full example code
.. rst-class:: sphx-glr-example-title
.. _sphx_glr_content_tutorials_07-fdem_plot_fwd_2_fem_cyl.py:
3D Forward Simulation on a Cylindrical Mesh
===========================================
Here we use the module *SimPEG.electromagnetics.frequency_domain* to simulate the
FDEM response for a borehole survey using a cylindrical mesh and radially symmetric
conductivity model. For this tutorial, we focus on the following:
- How to define the transmitters and receivers
- How to define the survey
- How to solve the FDEM problem on cylindrical meshes
- The units of the conductivity/resistivity model and resulting data
Please note that we have used a coarse mesh to shorten the time of the simulation.
Proper discretization is required to simulate the fields at each frequency with
sufficient accuracy.
.. GENERATED FROM PYTHON SOURCE LINES 23-26
Import modules
--------------
.. GENERATED FROM PYTHON SOURCE LINES 26-47
.. code-block:: default
from discretize import CylindricalMesh
from discretize.utils import mkvc
from SimPEG import maps
import SimPEG.electromagnetics.frequency_domain as fdem
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
try:
from pymatsolver import Pardiso as Solver
except ImportError:
from SimPEG import SolverLU as Solver
write_file = False
# sphinx_gallery_thumbnail_number = 2
.. GENERATED FROM PYTHON SOURCE LINES 48-55
Create Airborne Survey
----------------------
Here we define a x-offset borehole survey that consists of a single vertical line
of source-receiver pairs which measred the secondary magnetic flux density
over a range of frequencies.
.. GENERATED FROM PYTHON SOURCE LINES 55-94
.. code-block:: default
# Frequencies being predicted (10 Hz to 10000 Hz)
frequencies = np.logspace(1, 4, 16)
# Defining transmitter locations
xtx, ytx, ztx = np.meshgrid([0], [0], np.linspace(0, -500, 21))
source_locations = np.c_[mkvc(xtx), mkvc(ytx), mkvc(ztx)]
ntx = np.size(xtx)
# Define receiver locations
xrx, yrx, zrx = np.meshgrid([100], [0], np.linspace(0, -500, 21))
receiver_locations = np.c_[mkvc(xrx), mkvc(yrx), mkvc(zrx)]
source_list = [] # Create empty list to store sources
# Each unique location and frequency defines a new transmitter
for ii in range(ntx):
# Define receivers of different types at each location. Real and imaginary
# measurements require separate receivers. You can define the orientation of
# the transmitters and receivers for different survey geometries.
bzr_receiver = fdem.receivers.PointMagneticFluxDensitySecondary(
receiver_locations[ii, :], "z", "real"
)
bzi_receiver = fdem.receivers.PointMagneticFluxDensitySecondary(
receiver_locations[ii, :], "z", "imag"
)
receivers_list = [bzr_receiver, bzi_receiver] # must be a list
for jj in range(len(frequencies)):
# Must define the transmitter properties and associated receivers
source_list.append(
fdem.sources.MagDipole(
receivers_list, frequencies[jj], source_locations[ii], orientation="z"
)
)
# Define the survey
survey = fdem.Survey(source_list)
.. GENERATED FROM PYTHON SOURCE LINES 95-107
Create Cylindrical Mesh
-----------------------
Here we create the cylindrical mesh that will be used for this tutorial
example. We chose to design a coarser mesh to decrease the run time.
When designing a mesh to solve practical frequency domain problems:
- Your smallest cell size should be 10%-20% the size of your smallest skin depth
- The thickness of your padding needs to be 2-3 times biggest than your largest skin depth
- The skin depth is ~500*np.sqrt(rho/f)
.. GENERATED FROM PYTHON SOURCE LINES 107-117
.. code-block:: default
hr = [(10.0, 30), (10.0, 10, 1.5)] # discretization in the radial direction
hz = [
(10.0, 10, -1.5),
(10.0, 200),
(10.0, 10, 1.5),
] # discretization in vertical direction
mesh = CylindricalMesh([hr, 1, hz], x0="00C")
.. GENERATED FROM PYTHON SOURCE LINES 118-125
Create Conductivity/Resistivity Model and Mapping
-------------------------------------------------
Here, we create the model that will be used to predict frequency domain
data and the mapping from the model to the mesh. The model
consists of several layers. For this example, we will have only flat topography.
.. GENERATED FROM PYTHON SOURCE LINES 125-175
.. code-block:: default
# Conductivity in S/m (or resistivity in Ohm m)
air_conductivity = 1e-8
background_conductivity = 1e-1
layer_conductivity_1 = 1e0
layer_conductivity_2 = 1e-2
# Find cells that are active in the forward modeling (cells below surface)
ind_active = mesh.cell_centers[:, 2] < 0
# Define mapping from model to active cells
model_map = maps.InjectActiveCells(mesh, ind_active, air_conductivity)
# Define the model
model = background_conductivity * np.ones(ind_active.sum())
ind = (mesh.cell_centers[ind_active, 2] > -200.0) & (
mesh.cell_centers[ind_active, 2] < -0
)
model[ind] = layer_conductivity_1
ind = (mesh.cell_centers[ind_active, 2] > -400.0) & (
mesh.cell_centers[ind_active, 2] < -200
)
model[ind] = layer_conductivity_2
# Plot Conductivity Model
mpl.rcParams.update({"font.size": 14})
fig = plt.figure(figsize=(5, 6))
plotting_map = maps.InjectActiveCells(mesh, ind_active, np.nan)
log_model = np.log10(model)
ax1 = fig.add_axes([0.14, 0.1, 0.6, 0.85])
mesh.plot_image(
plotting_map * log_model,
ax=ax1,
grid=False,
clim=(np.log10(layer_conductivity_2), np.log10(layer_conductivity_1)),
)
ax1.set_title("Conductivity Model")
ax2 = fig.add_axes([0.76, 0.1, 0.05, 0.85])
norm = mpl.colors.Normalize(
vmin=np.log10(layer_conductivity_2), vmax=np.log10(layer_conductivity_1)
)
cbar = mpl.colorbar.ColorbarBase(
ax2, norm=norm, orientation="vertical", format="$10^{%.1f}$"
)
cbar.set_label("Conductivity [S/m]", rotation=270, labelpad=15, size=12)
.. image-sg:: /content/tutorials/07-fdem/images/sphx_glr_plot_fwd_2_fem_cyl_001.png
:alt: Conductivity Model
:srcset: /content/tutorials/07-fdem/images/sphx_glr_plot_fwd_2_fem_cyl_001.png
:class: sphx-glr-single-img
.. GENERATED FROM PYTHON SOURCE LINES 176-185
Simulation: Predicting FDEM Data
--------------------------------
Here we define the formulation for solving Maxwell's equations. Since we are
measuring the magnetic flux density and working with a conductivity model,
the EB formulation is the most natural. We must also remember to define
the mapping for the conductivity model. If you defined a resistivity model,
use the kwarg *rhoMap* instead of *sigmaMap*
.. GENERATED FROM PYTHON SOURCE LINES 185-190
.. code-block:: default
simulation = fdem.simulation.Simulation3DMagneticFluxDensity(
mesh, survey=survey, sigmaMap=model_map, solver=Solver
)
.. GENERATED FROM PYTHON SOURCE LINES 191-196
Predict and Plot Data
---------------------
Here we show how the simulation is used to predict data.
.. GENERATED FROM PYTHON SOURCE LINES 196-232
.. code-block:: default
# Compute predicted data for the given model.
dpred = simulation.dpred(model)
# Data are organized by transmitter location, then component, then frequency. We had nFreq
# transmitters and each transmitter had 2 receivers (real and imaginary component). So
# first we will pick out the real and imaginary data
bz_real = dpred[0 : len(dpred) : 2]
bz_imag = dpred[1 : len(dpred) : 2]
# Then we will will reshape the data.
bz_real = np.reshape(bz_real, (ntx, len(frequencies)))
bz_imag = np.reshape(bz_imag, (ntx, len(frequencies)))
# Plot secondary field along the profile at f = 10000 Hz
fig = plt.figure(figsize=(5, 5))
ax1 = fig.add_axes([0.2, 0.15, 0.75, 0.75])
frequencies_index = 0
ax1.plot(bz_real[:, frequencies_index], receiver_locations[:, -1], "b-", lw=3)
ax1.plot(bz_imag[:, frequencies_index], receiver_locations[:, -1], "b--", lw=3)
ax1.set_xlabel("Bz secondary [T]")
ax1.set_ylabel("Elevation [m]")
ax1.set_title("Response at 10000 Hz")
ax1.legend(["Real", "Imaginary"], loc="upper right")
# Plot FEM response for all frequencies
fig = plt.figure(figsize=(5, 5))
ax1 = fig.add_axes([0.2, 0.15, 0.75, 0.75])
location_index = 0
ax1.semilogx(frequencies, bz_real[location_index, :], "b-", lw=3)
ax1.semilogx(frequencies, bz_imag[location_index, :], "b--", lw=3)
ax1.set_xlim((np.min(frequencies), np.max(frequencies)))
ax1.set_xlabel("Frequency [Hz]")
ax1.set_ylabel("Bz secondary [T]")
ax1.set_title("Response at Smallest Offset")
ax1.legend(["Real", "Imaginary"], loc="upper left")
.. rst-class:: sphx-glr-horizontal
*
.. image-sg:: /content/tutorials/07-fdem/images/sphx_glr_plot_fwd_2_fem_cyl_002.png
:alt: Response at 10000 Hz
:srcset: /content/tutorials/07-fdem/images/sphx_glr_plot_fwd_2_fem_cyl_002.png
:class: sphx-glr-multi-img
*
.. image-sg:: /content/tutorials/07-fdem/images/sphx_glr_plot_fwd_2_fem_cyl_003.png
:alt: Response at Smallest Offset
:srcset: /content/tutorials/07-fdem/images/sphx_glr_plot_fwd_2_fem_cyl_003.png
:class: sphx-glr-multi-img
.. rst-class:: sphx-glr-script-out
.. code-block:: none
.. rst-class:: sphx-glr-timing
**Total running time of the script:** (0 minutes 8.336 seconds)
**Estimated memory usage:** 8 MB
.. _sphx_glr_download_content_tutorials_07-fdem_plot_fwd_2_fem_cyl.py:
.. only:: html
.. container:: sphx-glr-footer sphx-glr-footer-example
.. container:: sphx-glr-download sphx-glr-download-python
:download:`Download Python source code: plot_fwd_2_fem_cyl.py `
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download Jupyter notebook: plot_fwd_2_fem_cyl.ipynb `
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery `_