.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/tutorials/check_calc_domain.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note Click :ref:`here ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_gallery_tutorials_check_calc_domain.py: 08. Ensure computation domain is big enough =========================================== Ensure the boundary in :math:`\pm x`, :math:`\pm y`, and :math:`+ z` is big enough for :math:`\rho_\text{air}`. The air is very resistive, and EM waves propagate at the speed of light as a wave, not diffusive any longer. The whole concept of skin depth does therefore not apply to the air layer. The only attenuation is caused by geometrical spreading. In order to not have any effects from the boundary one has to choose the air layer appropriately. The important bit is that the resistivity of air has to be taken into account also for the horizontal directions, not only for positive :math:`z` (upwards into the sky). This is an example to test boundaries on a simple marine model (air, water, subsurface) and compare them to the 1D result. .. GENERATED FROM PYTHON SOURCE LINES 19-27 .. code-block:: default import emg3d import empymod import numpy as np import matplotlib.pyplot as plt from matplotlib.colors import LogNorm plt.style.use('ggplot') .. GENERATED FROM PYTHON SOURCE LINES 29-31 Model, Survey, and Analytical Solution -------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 31-45 .. code-block:: default water_depth = 500 # 500 m water depth offsets = np.linspace(2000, 7000, 501) # Offsets src_coo = [0, 0, -water_depth+50, 0, 0] # Src at origin, 50 m above seafloor rec_coo = [offsets, offsets*0, -water_depth, 0, 0] # Receivers on the seafloor depth = [-water_depth, 0] # Simple model resistivity = [1, 0.3, 1e8] # Simple model frequency = 0.1 # Frequency source = emg3d.TxElectricDipole(src_coo) # Compute analytical solution epm = empymod.bipole(src_coo, rec_coo, depth, resistivity, frequency) .. rst-class:: sphx-glr-script-out .. code-block:: none :: empymod END; runtime = 0:00:00.132017 :: 1 kernel call(s) .. GENERATED FROM PYTHON SOURCE LINES 46-48 3D Modelling ------------ .. GENERATED FROM PYTHON SOURCE LINES 48-64 .. code-block:: default # Parameter we keep the same for both grids grid_inp = { 'frequency': frequency, 'center': [src_coo[0], src_coo[1], -water_depth-100], 'domain': ([src_coo[0]-500, offsets[-1]+500], [src_coo[1], src_coo[1]], [-600, 100]), 'seasurface': 0, 'min_width_limits': 100, 'stretching': [1, 1.25], 'lambda_from_center': True, 'center_on_edge': True, 'verb': 1, } .. GENERATED FROM PYTHON SOURCE LINES 65-81 1st grid, only considering air resistivity for +z ------------------------------------------------- Here we are in the water, so the signal is attenuated before it enters the air. So we don't use the resistivity of air to compute the required boundary, but 100 Ohm.m instead. (100 is the result of a quick parameter test with :math:`\rho=1e4, 1e3, 1e2, 1e1`, and the result was that after 100 there is not much improvement any longer.) Also note that the function :func:`emg3d.meshes.construct_mesh` internally uses six times the skin depth for the boundary. For :math:`\rho` = 100 Ohm.m and :math:`f` = 0.1 Hz, the skin depth :math:`\delta` is roughly 16 km, which therefore results in a boundary of roughly 96 km. See the documentation of :func:`emg3d.meshes.get_hx_h0` for more information on how the grid is created. .. GENERATED FROM PYTHON SOURCE LINES 81-88 .. code-block:: default grid_1 = emg3d.construct_mesh( properties=[resistivity[1], resistivity[0], resistivity[0], 100], **grid_inp, ) grid_1 .. rst-class:: sphx-glr-script-out .. code-block:: none == GRIDDING IN X == Skin depth [m] : 872 / 1592 / 1592 [corr. to `properties`] Survey dom. DS [m] : -500 - 7500 Comp. dom. DC [m] : -10250 - 13750 Final extent [m] : -10405 - 14205 Cell widths [m] : 100 / 100 / 904 [min(DS) / max(DS) / max(DC)] Number of cells : 128 (80 / 48 / 0) [Total (DS/DC/remain)] Max stretching : 1.000 (1.000) / 1.088 [DS (seasurface) / DC] == GRIDDING IN Y == Skin depth [m] : 872 / 1592 / 1592 [corr. to `properties`] Survey dom. DS [m] : 0 - 0 Comp. dom. DC [m] : -10000 - 10000 Final extent [m] : -10226 - 10226 Cell widths [m] : 100 / 100 / 1907 [min(DS) / max(DS) / max(DC)] Number of cells : 32 (2 / 30 / 0) [Total (DS/DC/remain)] Max stretching : 1.000 (1.000) / 1.217 [DS (seasurface) / DC] == GRIDDING IN Z == Skin depth [m] : 872 / 1592 / 15915 [corr. to `properties`] Survey dom. DS [m] : -600 - 100 Comp. dom. DC [m] : -10600 - 99400 Final extent [m] : -12620 - 102179 Cell widths [m] : 100 / 100 / 19514 [min(DS) / max(DS) / max(DC)] Number of cells : 48 (8 / 40 / 0) [Total (DS/DC/remain)] Max stretching : 1.000 (1.000) / 1.235 [DS (seasurface) / DC] .. raw:: html
TensorMesh 196,608 cells
MESH EXTENT CELL WIDTH FACTOR
dir nC min max min max max
x 128 -10,404.94 14,204.82 100.00 904.35 1.09
y 32 -10,225.81 10,225.81 100.00 1,906.68 1.22
z 48 -12,620.23 102,179.01 100.00 19,513.80 1.23


.. GENERATED FROM PYTHON SOURCE LINES 89-104 .. code-block:: default # Create corresponding model res_1 = resistivity[0]*np.ones(grid_1.n_cells) res_1[grid_1.cell_centers[:, 2] > -water_depth] = resistivity[1] res_1[grid_1.cell_centers[:, 2] > 0] = resistivity[2] model_1 = emg3d.Model(grid_1, property_x=res_1, mapping='Resistivity') # QC grid_1.plot_3d_slicer( np.log10(model_1.property_x), zlim=(-2000, 100), clim=[-1, 2]) # Solve the system efield_1 = emg3d.solve_source(model_1, source, frequency, verb=3) .. image-sg:: /gallery/tutorials/images/sphx_glr_check_calc_domain_001.png :alt: check calc domain :srcset: /gallery/tutorials/images/sphx_glr_check_calc_domain_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none :: emg3d START :: 21:35:33 :: v1.8.0 MG-cycle : 'F' sslsolver : 'bicgstab' semicoarsening : True [1 2 3] tol : 1e-06 linerelaxation : True [4 5 6] maxit : 50 (3) nu_{i,1,c,2} : 0, 2, 1, 2 verb : 3 Original grid : 128 x 32 x 48 => 196,608 cells Coarsest grid : 2 x 2 x 3 => 12 cells Coarsest level : 6 ; 4 ; 4 :: emg3d :: 1.2e-02; 0(1); 0:00:04 :: emg3d :: 1.4e-03; 0(2); 0:00:07 :: emg3d :: 1.8e-04; 0(3); 0:00:10 :: emg3d :: 6.6e-05; 0(4); 0:00:14 :: emg3d :: 7.9e-06; 0(5); 0:00:17 :: emg3d :: 1.2e-06; 0(6); 0:00:20 :: emg3d :: 7.8e-07; 1(6); 0:00:20 > CONVERGED > Solver steps : 1 > MG prec. steps : 6 > Final rel. error : 7.844e-07 :: emg3d END :: 21:35:53 :: runtime = 0:00:20 .. GENERATED FROM PYTHON SOURCE LINES 105-109 2nd grid, considering air resistivity for +/- x, +/- y, and +z -------------------------------------------------------------- See comments below the heading of the 1st grid regarding boundary. .. GENERATED FROM PYTHON SOURCE LINES 109-116 .. code-block:: default grid_2 = emg3d.construct_mesh( properties=[resistivity[1], resistivity[0], 100], **grid_inp, ) grid_2 .. rst-class:: sphx-glr-script-out .. code-block:: none == GRIDDING IN X == Skin depth [m] : 872 / 15915 / 15915 [corr. to `properties`] Survey dom. DS [m] : -500 - 7500 Comp. dom. DC [m] : -100000 - 100000 Final extent [m] : -101682 - 108682 Cell widths [m] : 100 / 100 / 20173 [min(DS) / max(DS) / max(DC)] Number of cells : 128 (80 / 48 / 0) [Total (DS/DC/remain)] Max stretching : 1.000 (1.000) / 1.247 [DS (seasurface) / DC] == GRIDDING IN Y == Skin depth [m] : 872 / 15915 / 15915 [corr. to `properties`] Survey dom. DS [m] : 0 - 0 Comp. dom. DC [m] : -100000 - 100000 Final extent [m] : -102877 - 102877 Cell widths [m] : 100 / 100 / 15539 [min(DS) / max(DS) / max(DC)] Number of cells : 64 (2 / 62 / 0) [Total (DS/DC/remain)] Max stretching : 1.000 (1.000) / 1.177 [DS (seasurface) / DC] == GRIDDING IN Z == Skin depth [m] : 872 / 1592 / 15915 [corr. to `properties`] Survey dom. DS [m] : -600 - 100 Comp. dom. DC [m] : -10600 - 99400 Final extent [m] : -12620 - 102179 Cell widths [m] : 100 / 100 / 19514 [min(DS) / max(DS) / max(DC)] Number of cells : 48 (8 / 40 / 0) [Total (DS/DC/remain)] Max stretching : 1.000 (1.000) / 1.235 [DS (seasurface) / DC] .. raw:: html
TensorMesh 393,216 cells
MESH EXTENT CELL WIDTH FACTOR
dir nC min max min max max
x 128 -101,682.39 108,682.39 100.00 20,172.62 1.25
y 64 -102,877.15 102,877.15 100.00 15,538.63 1.18
z 48 -12,620.23 102,179.01 100.00 19,513.80 1.23


.. GENERATED FROM PYTHON SOURCE LINES 117-132 .. code-block:: default # Create corresponding model res_2 = resistivity[0]*np.ones(grid_2.n_cells) res_2[grid_2.cell_centers[:, 2] > -water_depth] = resistivity[1] res_2[grid_2.cell_centers[:, 2] > 0] = resistivity[2] model_2 = emg3d.Model(grid_2, property_x=res_2, mapping='Resistivity') # QC # grid_2.plot_3d_slicer( # np.log10(model_2.property_x), zlim=(-2000, 100), clim=[-1, 2]) # Define source and solve the system efield_2 = emg3d.solve_source(model_2, source, frequency, verb=3) .. rst-class:: sphx-glr-script-out .. code-block:: none :: emg3d START :: 21:35:54 :: v1.8.0 MG-cycle : 'F' sslsolver : 'bicgstab' semicoarsening : True [1 2 3] tol : 1e-06 linerelaxation : True [4 5 6] maxit : 50 (3) nu_{i,1,c,2} : 0, 2, 1, 2 verb : 3 Original grid : 128 x 64 x 48 => 393,216 cells Coarsest grid : 2 x 2 x 3 => 12 cells Coarsest level : 6 ; 5 ; 4 :: emg3d :: 1.0e-02; 0(1); 0:00:07 :: emg3d :: 1.9e-03; 0(2); 0:00:13 :: emg3d :: 5.0e-04; 0(3); 0:00:19 :: emg3d :: 6.0e-05; 0(4); 0:00:25 :: emg3d :: 1.0e-05; 0(5); 0:00:31 :: emg3d :: 3.2e-06; 0(6); 0:00:38 :: emg3d :: 2.5e-06; 1(6); 0:00:38 :: emg3d :: 2.2e-07; 1(7); 0:00:44 :: emg3d :: 2.3e-07; 2(7); 0:00:44 > CONVERGED > Solver steps : 2 > MG prec. steps : 7 > Final rel. error : 2.292e-07 :: emg3d END :: 21:36:38 :: runtime = 0:00:44 .. GENERATED FROM PYTHON SOURCE LINES 133-135 Plot receiver responses ----------------------- .. GENERATED FROM PYTHON SOURCE LINES 135-141 .. code-block:: default # Interpolate fields at receiver positions emg_1 = efield_1.get_receiver(tuple(rec_coo)) emg_2 = efield_2.get_receiver(tuple(rec_coo)) .. GENERATED FROM PYTHON SOURCE LINES 142-216 .. code-block:: default plt.figure(figsize=(10, 7)) # Real, log-lin ax1 = plt.subplot(321) plt.title('(a) lin-lin Real') plt.plot(offsets/1e3, epm.real, 'k', lw=2, label='analytical') plt.plot(offsets/1e3, emg_1.real, 'C0--', label='grid 1') plt.plot(offsets/1e3, emg_2.real, 'C1:', label='grid 2') plt.ylabel('$E_x$ (V/m)') plt.legend() # Real, log-symlog ax3 = plt.subplot(323, sharex=ax1) plt.title('(c) lin-symlog Real') plt.plot(offsets/1e3, epm.real, 'k') plt.plot(offsets/1e3, emg_1.real, 'C0--') plt.plot(offsets/1e3, emg_2.real, 'C1:') plt.ylabel('$E_x$ (V/m)') plt.yscale('symlog', linthresh=1e-15) # Real, error ax5 = plt.subplot(325, sharex=ax3) plt.title('(e) clipped 0.01-10') # Compute the error err_real_1 = np.clip(100*abs((epm.real-emg_1.real)/epm.real), 0.01, 10) err_real_2 = np.clip(100*abs((epm.real-emg_2.real)/epm.real), 0.01, 10) plt.ylabel('Rel. error %') plt.plot(offsets/1e3, err_real_1, 'C0--') plt.plot(offsets/1e3, err_real_2, 'C1:') plt.axhline(1, color='.4') plt.yscale('log') plt.ylim([0.008, 12]) plt.xlabel('Offset (km)') # Imaginary, log-lin ax2 = plt.subplot(322) plt.title('(b) lin-lin Imag') plt.plot(offsets/1e3, epm.imag, 'k') plt.plot(offsets/1e3, emg_1.imag, 'C0--') plt.plot(offsets/1e3, emg_2.imag, 'C1:') # Imaginary, log-symlog ax4 = plt.subplot(324, sharex=ax2) plt.title('(d) lin-symlog Imag') plt.plot(offsets/1e3, epm.imag, 'k') plt.plot(offsets/1e3, emg_1.imag, 'C0--') plt.plot(offsets/1e3, emg_2.imag, 'C1:') plt.yscale('symlog', linthresh=1e-15) # Imaginary, error ax6 = plt.subplot(326, sharex=ax2) plt.title('(f) clipped 0.01-10') # Compute error err_imag_1 = np.clip(100*abs((epm.imag-emg_1.imag)/epm.imag), 0.01, 10) err_imag_2 = np.clip(100*abs((epm.imag-emg_2.imag)/epm.imag), 0.01, 10) plt.plot(offsets/1e3, err_imag_1, 'C0--') plt.plot(offsets/1e3, err_imag_2, 'C1:') plt.axhline(1, color='.4') plt.yscale('log') plt.ylim([0.008, 12]) plt.xlabel('Offset (km)') plt.tight_layout() plt.show() .. image-sg:: /gallery/tutorials/images/sphx_glr_check_calc_domain_002.png :alt: (a) lin-lin Real, (c) lin-symlog Real, (e) clipped 0.01-10, (b) lin-lin Imag, (d) lin-symlog Imag, (f) clipped 0.01-10 :srcset: /gallery/tutorials/images/sphx_glr_check_calc_domain_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 217-227 Plot entire fields to analyze and compare ----------------------------------------- 1st grid ```````` Upper plot shows the entire grid. One can see that the airwave attenuates to amplitudes in the order of 1e-17 at the boundary, absolutely good enough. However, the amplitudes in the horizontal directions are very high even at the boundaries :math:`\pm x` and :math:`\pm y`. .. GENERATED FROM PYTHON SOURCE LINES 227-239 .. code-block:: default grid_1.plot_3d_slicer( efield_1.fx.ravel('F'), view='abs', v_type='Ex', xslice=src_coo[0], yslice=src_coo[1], zslice=rec_coo[2], pcolor_opts={'norm': LogNorm(vmin=1e-17, vmax=1e-9)}) grid_1.plot_3d_slicer( efield_1.fx.ravel('F'), view='abs', v_type='Ex', zlim=[-5000, 1000], xslice=src_coo[0], yslice=src_coo[1], zslice=rec_coo[2], pcolor_opts={'norm': LogNorm(vmin=1e-17, vmax=1e-9)}) .. rst-class:: sphx-glr-horizontal * .. image-sg:: /gallery/tutorials/images/sphx_glr_check_calc_domain_003.png :alt: check calc domain :srcset: /gallery/tutorials/images/sphx_glr_check_calc_domain_003.png :class: sphx-glr-multi-img * .. image-sg:: /gallery/tutorials/images/sphx_glr_check_calc_domain_004.png :alt: check calc domain :srcset: /gallery/tutorials/images/sphx_glr_check_calc_domain_004.png :class: sphx-glr-multi-img .. GENERATED FROM PYTHON SOURCE LINES 240-246 2nd grid ```````` Again, upper plot shows the entire grid. One can see that the field attenuated sufficiently in all directions. Lower plot shows the same cut-out as the lower plot for the first grid, our zone of interest. .. GENERATED FROM PYTHON SOURCE LINES 246-260 .. code-block:: default grid_2.plot_3d_slicer( efield_2.fx.ravel('F'), view='abs', v_type='Ex', xslice=src_coo[0], yslice=src_coo[1], zslice=rec_coo[2], pcolor_opts={'norm': LogNorm(vmin=1e-17, vmax=1e-9)}) grid_2.plot_3d_slicer( efield_2.fx.ravel('F'), view='abs', v_type='Ex', xlim=[grid_1.nodes_x[0], grid_1.nodes_x[-1]], # Same square as grid_1 ylim=[grid_1.nodes_y[0], grid_1.nodes_y[-1]], # Same square as grid_1 zlim=[-5000, 1000], xslice=src_coo[0], yslice=src_coo[1], zslice=rec_coo[2], pcolor_opts={'norm': LogNorm(vmin=1e-17, vmax=1e-9)}) .. rst-class:: sphx-glr-horizontal * .. image-sg:: /gallery/tutorials/images/sphx_glr_check_calc_domain_005.png :alt: check calc domain :srcset: /gallery/tutorials/images/sphx_glr_check_calc_domain_005.png :class: sphx-glr-multi-img * .. image-sg:: /gallery/tutorials/images/sphx_glr_check_calc_domain_006.png :alt: check calc domain :srcset: /gallery/tutorials/images/sphx_glr_check_calc_domain_006.png :class: sphx-glr-multi-img .. GENERATED FROM PYTHON SOURCE LINES 261-263 .. code-block:: default emg3d.Report() .. raw:: html
Wed Aug 31 21:36:42 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


.. rst-class:: sphx-glr-timing **Total running time of the script:** ( 1 minutes 10.291 seconds) **Estimated memory usage:** 270 MB .. _sphx_glr_download_gallery_tutorials_check_calc_domain.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: check_calc_domain.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: check_calc_domain.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_