.. 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_3_fem_3d.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_3_fem_3d.py: 3D Forward Simulation on a Tree Mesh ==================================== Here we use the module *simpeg.electromagnetics.frequency_domain* to simulate the FDEM response for an airborne survey using an OcTree mesh and a conductivity/resistivity model. To limit computational demant, we simulate airborne data at a single frequency for a vertical coplanar survey geometry. This tutorial can be easily adapted to simulate data at many frequencies. For this tutorial, we focus on the following: - How to define the transmitters and receivers - How to define the survey - How to define the topography - How to solve the FDEM problem on OcTree 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 27-30 Import modules -------------- .. GENERATED FROM PYTHON SOURCE LINES 30-49 .. code-block:: Python from discretize import TreeMesh from discretize.utils import mkvc, refine_tree_xyz, active_from_xyz from simpeg.utils import plot2Ddata from simpeg import maps import simpeg.electromagnetics.frequency_domain as fdem import os import numpy as np import matplotlib as mpl import matplotlib.pyplot as plt save_file = False # sphinx_gallery_thumbnail_number = 2 .. GENERATED FROM PYTHON SOURCE LINES 50-57 Defining Topography ------------------- Here we define surface topography as an (N, 3) numpy array. Topography could also be loaded from a file. Here we define flat topography, however more complex topographies can be considered. .. GENERATED FROM PYTHON SOURCE LINES 57-62 .. code-block:: Python xx, yy = np.meshgrid(np.linspace(-3000, 3000, 101), np.linspace(-3000, 3000, 101)) zz = np.zeros(np.shape(xx)) topo_xyz = np.c_[mkvc(xx), mkvc(yy), mkvc(zz)] .. GENERATED FROM PYTHON SOURCE LINES 63-70 Create Airborne Survey ---------------------- For this example, the survey consists of a uniform grid of airborne measurements. To save time, we will only compute the response for a single frequency. .. GENERATED FROM PYTHON SOURCE LINES 70-111 .. code-block:: Python # Frequencies being predicted frequencies = [100, 500, 2500] # Defining transmitter locations N = 9 xtx, ytx, ztx = np.meshgrid(np.linspace(-200, 200, N), np.linspace(-200, 200, N), [40]) source_locations = np.c_[mkvc(xtx), mkvc(ytx), mkvc(ztx)] ntx = np.size(xtx) # Define receiver locations xrx, yrx, zrx = np.meshgrid(np.linspace(-200, 200, N), np.linspace(-200, 200, N), [20]) 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(len(frequencies)): for jj in range(ntx): # Define receivers of different type at each location bzr_receiver = fdem.receivers.PointMagneticFluxDensitySecondary( receiver_locations[jj, :], "z", "real" ) bzi_receiver = fdem.receivers.PointMagneticFluxDensitySecondary( receiver_locations[jj, :], "z", "imag" ) receivers_list = [bzr_receiver, bzi_receiver] # Must define the transmitter properties and associated receivers source_list.append( fdem.sources.MagDipole( receivers_list, frequencies[ii], source_locations[jj], orientation="z", moment=100, ) ) survey = fdem.Survey(source_list) .. GENERATED FROM PYTHON SOURCE LINES 112-124 Create OcTree Mesh ------------------ Here we define the OcTree mesh that is used for this 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 124-150 .. code-block:: Python dh = 25.0 # base cell width dom_width = 3000.0 # domain width nbc = 2 ** int(np.round(np.log(dom_width / dh) / np.log(2.0))) # num. base cells # Define the base mesh h = [(dh, nbc)] mesh = TreeMesh([h, h, h], x0="CCC") # Mesh refinement based on topography mesh = refine_tree_xyz( mesh, topo_xyz, octree_levels=[0, 0, 0, 1], method="surface", finalize=False ) # Mesh refinement near transmitters and receivers mesh = refine_tree_xyz( mesh, receiver_locations, octree_levels=[2, 4], method="radial", finalize=False ) # Refine core mesh region xp, yp, zp = np.meshgrid([-250.0, 250.0], [-250.0, 250.0], [-300.0, 0.0]) xyz = np.c_[mkvc(xp), mkvc(yp), mkvc(zp)] mesh = refine_tree_xyz(mesh, xyz, octree_levels=[0, 2, 4], method="box", finalize=False) mesh.finalize() .. rst-class:: sphx-glr-script-out .. code-block:: none /home/vsts/work/1/s/tutorials/07-fdem/plot_fwd_3_fem_3d.py:134: DeprecationWarning: The surface option is deprecated as of `0.9.0` please update your code to use the `TreeMesh.refine_surface` functionality. It will be removed in a future version of discretize. /home/vsts/work/1/s/tutorials/07-fdem/plot_fwd_3_fem_3d.py:139: DeprecationWarning: The radial option is deprecated as of `0.9.0` please update your code to use the `TreeMesh.refine_points` functionality. It will be removed in a future version of discretize. /home/vsts/work/1/s/tutorials/07-fdem/plot_fwd_3_fem_3d.py:146: DeprecationWarning: The box option is deprecated as of `0.9.0` please update your code to use the `TreeMesh.refine_bounding_box` functionality. It will be removed in a future version of discretize. .. GENERATED FROM PYTHON SOURCE LINES 151-159 Defining the 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. Here, the model consists of a conductive block within a more resistive background. .. GENERATED FROM PYTHON SOURCE LINES 159-211 .. code-block:: Python # Conductivity in S/m (or resistivity in Ohm m) air_conductivity = 1e-8 background_conductivity = 1e-2 block_conductivity = 1e1 # Find cells that are active in the forward modeling (cells below surface) ind_active = active_from_xyz(mesh, topo_xyz) # Define mapping from model to active cells model_map = maps.InjectActiveCells(mesh, ind_active, air_conductivity) # Define model. Models in SimPEG are vector arrays model = background_conductivity * np.ones(ind_active.sum()) ind_block = ( (mesh.gridCC[ind_active, 0] < 100.0) & (mesh.gridCC[ind_active, 0] > -100.0) & (mesh.gridCC[ind_active, 1] < 100.0) & (mesh.gridCC[ind_active, 1] > -100.0) & (mesh.gridCC[ind_active, 2] > -275.0) & (mesh.gridCC[ind_active, 2] < -75.0) ) model[ind_block] = block_conductivity # Plot Resistivity Model mpl.rcParams.update({"font.size": 12}) fig = plt.figure(figsize=(7, 6)) plotting_map = maps.InjectActiveCells(mesh, ind_active, np.nan) log_model = np.log10(model) ax1 = fig.add_axes([0.13, 0.1, 0.6, 0.85]) mesh.plot_slice( plotting_map * log_model, normal="Y", ax=ax1, ind=int(mesh.h[0].size / 2), grid=False, clim=(np.log10(background_conductivity), np.log10(block_conductivity)), ) ax1.set_title("Conductivity Model at Y = 0 m") ax2 = fig.add_axes([0.75, 0.1, 0.05, 0.85]) norm = mpl.colors.Normalize( vmin=np.log10(background_conductivity), vmax=np.log10(block_conductivity) ) 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_3_fem_3d_001.png :alt: Conductivity Model at Y = 0 m :srcset: /content/tutorials/07-fdem/images/sphx_glr_plot_fwd_3_fem_3d_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 212-221 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 221-226 .. code-block:: Python simulation = fdem.simulation.Simulation3DMagneticFluxDensity( mesh, survey=survey, sigmaMap=model_map ) .. GENERATED FROM PYTHON SOURCE LINES 227-232 Predict and Plot Data --------------------- Here we show how the simulation is used to predict data. .. GENERATED FROM PYTHON SOURCE LINES 232-292 .. code-block:: Python # Compute predicted data for a your model. dpred = simulation.dpred(model) # Data are organized by frequency, transmitter location, then by receiver. We 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 for plotting. bz_real_plotting = np.reshape(bz_real, (len(frequencies), ntx)) bz_imag_plotting = np.reshape(bz_imag, (len(frequencies), ntx)) fig = plt.figure(figsize=(10, 4)) # Real Component frequencies_index = 0 v_max = np.max(np.abs(bz_real_plotting[frequencies_index, :])) ax1 = fig.add_axes([0.05, 0.05, 0.35, 0.9]) plot2Ddata( receiver_locations[:, 0:2], bz_real_plotting[frequencies_index, :], ax=ax1, ncontour=30, clim=(-v_max, v_max), contourOpts={"cmap": "bwr"}, ) ax1.set_title("Re[$B_z$] at 100 Hz") ax2 = fig.add_axes([0.41, 0.05, 0.02, 0.9]) norm = mpl.colors.Normalize(vmin=-v_max, vmax=v_max) cbar = mpl.colorbar.ColorbarBase( ax2, norm=norm, orientation="vertical", cmap=mpl.cm.bwr ) cbar.set_label("$T$", rotation=270, labelpad=15, size=12) # Imaginary Component v_max = np.max(np.abs(bz_imag_plotting[frequencies_index, :])) ax1 = fig.add_axes([0.55, 0.05, 0.35, 0.9]) plot2Ddata( receiver_locations[:, 0:2], bz_imag_plotting[frequencies_index, :], ax=ax1, ncontour=30, clim=(-v_max, v_max), contourOpts={"cmap": "bwr"}, ) ax1.set_title("Im[$B_z$] at 100 Hz") ax2 = fig.add_axes([0.91, 0.05, 0.02, 0.9]) norm = mpl.colors.Normalize(vmin=-v_max, vmax=v_max) cbar = mpl.colorbar.ColorbarBase( ax2, norm=norm, orientation="vertical", cmap=mpl.cm.bwr ) cbar.set_label("$T$", rotation=270, labelpad=15, size=12) plt.show() .. image-sg:: /content/tutorials/07-fdem/images/sphx_glr_plot_fwd_3_fem_3d_002.png :alt: Re[$B_z$] at 100 Hz, Im[$B_z$] at 100 Hz :srcset: /content/tutorials/07-fdem/images/sphx_glr_plot_fwd_3_fem_3d_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/vsts/work/1/s/simpeg/simulation.py:197: DefaultSolverWarning: Using the default solver: Pardiso. If you would like to suppress this notification, add warnings.filterwarnings('ignore', simpeg.utils.solver_utils.DefaultSolverWarning) to your script. .. GENERATED FROM PYTHON SOURCE LINES 293-298 Optional: Export Data --------------------- Write the true model, data and topography .. GENERATED FROM PYTHON SOURCE LINES 298-324 .. code-block:: Python if save_file: dir_path = os.path.dirname(fdem.__file__).split(os.path.sep)[:-3] dir_path.extend(["tutorials", "assets", "fdem"]) dir_path = os.path.sep.join(dir_path) + os.path.sep # Write topography fname = dir_path + "fdem_topo.txt" np.savetxt(fname, np.c_[topo_xyz], fmt="%.4e") # Write data with 2% noise added fname = dir_path + "fdem_data.obs" bz_real = bz_real + 1e-14 * np.random.randn(len(bz_real)) bz_imag = bz_imag + 1e-14 * np.random.randn(len(bz_imag)) f_vec = np.kron(frequencies, np.ones(ntx)) receiver_locations = np.kron(np.ones((len(frequencies), 1)), receiver_locations) np.savetxt(fname, np.c_[f_vec, receiver_locations, bz_real, bz_imag], fmt="%.4e") # Plot true model output_model = plotting_map * model output_model[np.isnan(output_model)] = 1e-8 fname = dir_path + "true_model.txt" np.savetxt(fname, output_model, fmt="%.4e") .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 16.562 seconds) **Estimated memory usage:** 1042 MB .. _sphx_glr_download_content_tutorials_07-fdem_plot_fwd_3_fem_3d.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_fwd_3_fem_3d.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_fwd_3_fem_3d.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_fwd_3_fem_3d.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_