.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/tutorials/timedomain.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_gallery_tutorials_timedomain.py: 09. Transient CSEM ================== The computation of ``emg3d`` happens in the frequency domain (or Laplace domain), each frequency requires a new computation. Using (inverse) Fourier transforms, we can also compute time-domain (transient) CSEM data with ``emg3d``. This is not (yet) implemented in easy, user-friendly functions. It does require quite some input and knowledge from the user, particularly with regards to the gridding. A good starting point to model time-domain data with ``emg3d`` is [WeMS21]_. You can find the paper and all the notebooks in the repo https://github.com/emsig/article-TDEM . The following is a simple example from the above article, the time-domain modelling of a fullspace. It based on the first example (Figures 3-4) of [MuWS08]_. Interactive frequency selection ------------------------------- The most important factor in fast time-domain computation is the frequency selection. You can find an interactive GUI for this in the repo https://github.com/emsig/frequency-selection . A screenshot of the GUI for the interactive frequency selection is shown in the following figure: .. figure:: ../../_static/images/GUI-freqselect.png :scale: 66 % :align: center :alt: Frequency-selection App :name: freqselect The GUI uses the 1D modeller ``empymod`` for a layered model, and internally the ``Fourier`` class of the 3D modeller ``emg3d``. The following parameters can be specified interactively: - points per decade - frequency range (min/max) - offset - Fourier transform (FFTLog or DLF with different filters) - signal (impulse or switch-on/-off) Other parameters have to be specified fix when initiating the widget. .. GENERATED FROM PYTHON SOURCE LINES 51-58 .. code-block:: Python import emg3d import empymod import numpy as np import matplotlib.pyplot as plt plt.style.use('bmh') .. GENERATED FROM PYTHON SOURCE LINES 59-72 Model and Survey ---------------- Model ````` - Homogeneous fullspace of 1 Ohm.m. Survey `````` - Source at origin. - Receiver at an inline-offset of 900 m. - Both source and receiver are x-directed electric dipoles. .. GENERATED FROM PYTHON SOURCE LINES 72-80 .. code-block:: Python src_coo = [0, 0, 0, 0, 0] source = emg3d.TxElectricDipole(src_coo) rec_coo = [900, 0, 0, 0, 0] resistivity = 1 # Fullspace resistivity depth = [] .. GENERATED FROM PYTHON SOURCE LINES 81-89 Fourier Transforms parameters ----------------------------- We only compute frequencies :math:`0.05 < f < 21` Hz, which yields enough precision for our purpose. This means, instead of 30 frequencies from 0.0002 - 126.4 Hz, we only need 14 frequencies from 0.05 - 20.0 Hz. .. GENERATED FROM PYTHON SOURCE LINES 89-110 .. code-block:: Python # Define desired times. time = np.logspace(-2, 1, 201) # Initiate a Fourier instance Fourier = emg3d.Fourier( time=time, fmin=0.05, fmax=21, ft='fftlog', # Fourier transform to use ftarg={'pts_per_dec': 5, 'add_dec': [-2, 1], 'q': 0}, ) # Dense frequencies for comparison reasons freq_dense = np.logspace( np.log10(Fourier.freq_required.min()), np.log10(Fourier.freq_required.max()), 301, ) .. rst-class:: sphx-glr-script-out .. code-block:: none time [s] : 0.01 - 10 : 201 [min-max; #] Fourier : FFTLog > pts_per_dec : 5 > add_dec : [-2. 1.] > q : 0.0 Req. freq [Hz] : 0.000200364 - 126.421 : 30 [min-max; #] Calc. freq [Hz] : 0.0503292 - 20.0364 : 14 [min-max; #] .. GENERATED FROM PYTHON SOURCE LINES 111-113 Frequency-domain computation ---------------------------- .. GENERATED FROM PYTHON SOURCE LINES 113-173 .. code-block:: Python # Automatic gridding settings. grid_opts = { 'center': src_coo[:3], # Source location 'domain': [[-200, 1100], [-50, 50], [-50, 50]], 'properties': resistivity, # Fullspace resistivity. 'min_width_limits': [20., 40.], # Restrict cell width within survey domain 'min_width_pps': 12, # Many points to have small min cell width 'stretching': [1, 1.3], # 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.) data[-fi-1] = efield.get_receiver(rec_coo) 2/14 :: 12.642126 Hz 3/14 :: 7.976643 Hz 4/14 :: 5.032921 Hz 5/14 :: 3.175559 Hz 6/14 :: 2.003642 Hz 7/14 :: 1.264213 Hz 8/14 :: 0.797664 Hz 9/14 :: 0.503292 Hz 10/14 :: 0.317556 Hz 11/14 :: 0.200364 Hz 12/14 :: 0.126421 Hz 13/14 :: 0.079766 Hz 14/14 :: 0.050329 Hz .. GENERATED FROM PYTHON SOURCE LINES 174-190 .. code-block:: Python runtime = 0 for freq in Fourier.freq_compute[::-1]: value = log[str(int(freq*1e6))] print(f" {value['freq']:7.3f} Hz: {value['info']['it_mg']:2g}/" f"{value['info']['it_ssl']:g} it; " f"{value['info']['time']:4.0f} s; " f"max_a: {value['stretching']:.2f}; " f"nC: {value['nC']:8,.0f}; " f"h: {value['dminmax'][0]:5.0f} / {value['dminmax'][1]:7.0f}") runtime += value['info']['time'] print(f"\n **** TOTAL RUNTIME :: {runtime // 60:.0f} min " f"{runtime % 60:.1f} s ****\n") .. rst-class:: sphx-glr-script-out .. code-block:: none 20.036 Hz: 6/1 it; 2 s; max_a: 1.26; nC: 46,080; h: 20 / 208 12.642 Hz: 6/1 it; 4 s; max_a: 1.16; nC: 98,304; h: 20 / 167 7.977 Hz: 6/1 it; 4 s; max_a: 1.19; nC: 98,304; h: 20 / 239 5.033 Hz: 6/1 it; 4 s; max_a: 1.22; nC: 98,304; h: 20 / 340 3.176 Hz: 6/1 it; 3 s; max_a: 1.25; nC: 81,920; h: 24 / 443 2.004 Hz: 6/1 it; 3 s; max_a: 1.23; nC: 81,920; h: 30 / 558 1.264 Hz: 6/1 it; 2 s; max_a: 1.21; nC: 65,536; h: 37 / 620 0.798 Hz: 6/1 it; 2 s; max_a: 1.23; nC: 65,536; h: 40 / 863 0.503 Hz: 6/1 it; 2 s; max_a: 1.25; nC: 65,536; h: 40 / 1200 0.318 Hz: 6/1 it; 2 s; max_a: 1.28; nC: 65,536; h: 40 / 1658 0.200 Hz: 6/1 it; 4 s; max_a: 1.27; nC: 102,400; h: 40 / 1825 0.126 Hz: 6/1 it; 4 s; max_a: 1.30; nC: 102,400; h: 40 / 2564 0.080 Hz: 6/1 it; 5 s; max_a: 1.25; nC: 128,000; h: 40 / 2974 0.050 Hz: 6/1 it; 5 s; max_a: 1.27; nC: 128,000; h: 40 / 3909 **** TOTAL RUNTIME :: 0 min 44.8 s **** .. GENERATED FROM PYTHON SOURCE LINES 191-195 Frequency domain ---------------- Compute analytical result and interpolate missing responses .. GENERATED FROM PYTHON SOURCE LINES 195-208 .. code-block:: Python data_int = Fourier.interpolate(data) # Compute analytical result using empymod epm_req = empymod.bipole(src_coo, rec_coo, depth, resistivity, Fourier.freq_required, verb=1) epm_dense = empymod.bipole(src_coo, rec_coo, depth, resistivity, freq_dense, verb=1) # Compute error err = np.clip(100*abs((data_int.imag-epm_req.imag)/epm_req.imag), 0.1, 100) .. GENERATED FROM PYTHON SOURCE LINES 209-213 Time domain ----------- Do the transform and compute analytical result. .. GENERATED FROM PYTHON SOURCE LINES 213-225 .. code-block:: Python # Compute corresponding time-domain signal. data_time = Fourier.freq2time(data, rec_coo[0]) # Analytical result epm_time = empymod.analytical(src_coo[:3], rec_coo[:3], resistivity, time, solution='dfs', signal=0, verb=1) # Relative error and peak error err_egd = 100*abs((data_time-epm_time)/epm_time) .. GENERATED FROM PYTHON SOURCE LINES 226-228 Plot it ``````` .. GENERATED FROM PYTHON SOURCE LINES 228-293 .. code-block:: Python plt.figure(figsize=(9, 5), constrained_layout=True) # Frequency-domain, imaginary, log-log ax1 = plt.subplot2grid((4, 2), (0, 0), rowspan=3) plt.title('(a) frequency domain') plt.plot(freq_dense, 1e9*abs(epm_dense.imag), 'C3', label='analytical') plt.plot(Fourier.freq_compute, 1e9*abs(data.imag), 'C0o', label='computed') plt.plot(Fourier.freq_required[~Fourier.ifreq_compute], 1e9*abs(data_int[~Fourier.ifreq_compute].imag), 'k.', label='interpolated / 0') plt.ylabel(r'$|\Im\{E_x\}|$ (nV/m)') plt.xscale('log') plt.yscale('symlog', linthresh=5e-9) plt.ylim([-1e-9, 5e-1]) ax1.set_xticklabels([]) plt.legend() plt.grid(axis='y', c='0.9') # Frequency-domain, imaginary, error ax2 = plt.subplot2grid((4, 2), (3, 0)) plt.plot(Fourier.freq_required, err, '.4') plt.plot(Fourier.freq_required[~Fourier.ifreq_compute], err[~Fourier.ifreq_compute], 'k.') plt.plot(Fourier.freq_compute, err[Fourier.ifreq_compute], 'C0o') plt.axhline(1, color='0.4', zorder=1) plt.xscale('log') plt.yscale('log') plt.xlabel('Frequency (Hz)') plt.ylabel('Rel. error %') plt.ylim([8e-3, 120]) plt.yticks([0.01, 0.1, 1, 10, 100], ('0.01', '0.1', '1', '10', '100')) plt.grid(axis='y', c='0.9') # Time-domain ax3 = plt.subplot2grid((4, 2), (0, 1), rowspan=3) plt.title('(b) time domain') plt.plot(time, epm_time*1e9, 'C3', lw=2, label='analytical') plt.plot(time, data_time*1e9, 'k--', label='transformed') plt.xlim([0, 2]) plt.ylabel('$E_x$ (nV/m)') ax3.set_xticklabels([]) plt.legend() ax3.yaxis.tick_right() ax3.yaxis.set_label_position("right") plt.grid(axis='y', c='0.9') # Time-domain, error ax4 = plt.subplot2grid((4, 2), (3, 1)) plt.plot(time, err_egd, 'k') plt.axhline(1, color='0.4', zorder=1) plt.yscale('log') plt.xlabel('Time (s)') plt.ylabel('Rel. error %') plt.xlim([0, 2]) plt.ylim([8e-3, 120]) plt.yticks([0.01, 0.1, 1, 10, 100], ('0.01', '0.1', '1', '10', '100')) ax4.yaxis.tick_right() ax4.yaxis.set_label_position("right") plt.grid(axis='y', c='0.9') .. image-sg:: /gallery/tutorials/images/sphx_glr_timedomain_001.png :alt: (a) frequency domain, (b) time domain :srcset: /gallery/tutorials/images/sphx_glr_timedomain_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 294-296 .. code-block:: Python emg3d.Report() .. raw:: html
Thu Jul 04 15:38:00 2024 CEST
OS Linux (Ubuntu 22.04) CPU(s) 16 Machine x86_64
Architecture 64bit RAM 31.0 GiB Environment Python
File system ext4
Python 3.12.4 | packaged by conda-forge | (main, Jun 17 2024, 10:23:07) [GCC 12.3.0]
numpy 1.26.4 scipy 1.14.0 numba 0.60.0
emg3d 1.8.3 empymod 2.3.1 xarray 2024.6.1.dev31+g6c2d8c33
discretize 0.10.0 h5py 3.11.0 matplotlib 3.8.4
tqdm 4.66.4
Intel(R) oneAPI Math Kernel Library Version 2023.2-Product Build 20230613 for Intel(R) 64 architecture applications


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