Reproducing TDEM response of pyGIMLi in empymod.

[ahsan435]

I am trying to convert Pygimli TDEM code into Empymod TDEM but not getting the same results can someone help me to resolve this issue.

Pygimli code

import numpy as np
import matplotlib.pyplot as plt
import pygimli as pg
from pygimli.physics.em import VMDTimeDomainModelling
from pygimli.viewer.mpl import drawModel1D
from pygimli.frameworks import MarquardtInversion
from scipy.interpolate import interp1d

# Define the synthetic model parameters
synThk = [300, 100]
synRes = [500, 10, 100]
nlay = len(synRes)  # number of layers
errorabsTEM_data = 1.  # absolute (ppm of primary signal)

fig, (ax1) = plt.subplots(figsize=(5, 5), ncols=1)

drawModel1D(ax1, synThk, synRes, plot='semilogx', color='C0', label="synth")

# MaxMin-10/Promys instrument, 1m above the ground
times = np.logspace(np.log10(1e-4), np.log10(1e-1), 64)

# Initialize VMDTimeDomainModelling
TEM_modeling = VMDTimeDomainModelling(times=times, nLayers=nlay, txArea=62500, verbose=True)

# Extract data
TEM_data = TEM_modeling(synThk + synRes)

#Plotting subplots
plt.figure(figsize=(7, 5))
plt.loglog(times, TEM_data, label='Observed Data', marker='o')
plt.show()

Empymod Code

import empymod
import numpy as np
import matplotlib.pyplot as plt
from pygimli.viewer.mpl import drawModel1D
plt.style.use('ggplot')

# Model
depth = [300, 100]       # Adjusted depth to match the number of layers
res = [500, 10, 100]      # Adjusted res to match the number of layers
#epermH = [1, 1, 1]        # Adjusted epermH to match the number of layers
model = {
    'src': [0, 0, 62500],           # Src at origin, slightly below interface
    'rec': [0, 0, 62500],          # 6 km off., in-line, slightly bel. interf.
    'depth': depth,                 # Target of 100 m at 2 km depth
    'res': res,                     # Res: [air, overb., target, half-space]
    'epermH': epermH,               # El. permittivity: default values
    'freqtime': t,                  # Times
    'signal': 0,                    # 0: Impulse response
    'ftarg': {'dlf': 'key_81_CosSin_2009'},  # Shorter filter then the default
    'verb': 1,                      # Verbosity; set to 3 to see all parameters
}


# Times (s)
t = np.logspace(np.log10(1e-4), np.log10(1e-1), 64)

# Compute with default eperm_air = 1
res_1 = empymod.dipole(**model)

# Set horizontal and vertical electric permittivity of air to 0
model['epermH'][0] = 0

# Compute with eperm_air = 0
res_0 = empymod.dipole(**model)

plt.figure()

plt.title('Time domain impulse response')
plt.semilogx(t, res_1, label=r'$\varepsilon_{air}=1$')
plt.semilogx(t, res_0, label=r'$\varepsilon_{air}=0$')
plt.xlabel(r'Time $[s]$')
plt.ylabel(r'Amplitude $[V/(m\,s)]$')
plt.legend()

plt.show()

[prisae]

Hi @ahsan435 ! I haven't run your code, just briefly looked at it. You want to run a VMD, right? Then you need to set ab=66 I assume, currently you do not set it what means you model x-directed electric sources and receivers. See https://empymod.emsig.xyz/en/stable/api/empymod.model.dipole.html where the parameter ab is described.

I close it for now, feel free to re-open if you cannot get it to work.

[prisae]

Also, you set your source and receiver at 62500 m depth, that does not make sense. There are various examples in the gallery, have a look how to define it properly (e.g., https://empymod.emsig.xyz/en/stable/gallery/tdomain/, https://empymod.emsig.xyz/en/stable/gallery/reproducing/ward1988.html).

[ahsan435]

Thank you very much for your response

[ahsan435]

@prisae Sir, I am still not getting the same results. I am new to coding, so that's why I'm having problems with it. Can you tell me if, for this code, a dipole will be more suitable than a loop?

[prisae]

I am afraid I do not have the time now to look into your code. But the documentation and the gallery is quite extensive, so you should hopefully find the necessary information!

For such a large area (62500) it might be best to model the loop as a sequence/sum of dipoles.

[ahsan435]

Thank you very much Sir I will try that one

[ahsan435]

@prisae Dear Sir I applied your comments on the code but still not getting the same results, is there anything still missing in the code please guide me. I know that you are very busy but I am really stuck in it and need help. Thank you very much.

import numpy as np
import matplotlib.pyplot as plt
import empymod
from scipy.interpolate import interp1d

# Define the synthetic model parameters
synThk = [300, 100]
synRes = [500, 10, 100]
nlay = len(synRes)  # number of layers
errorabsTEM_data = 1.  # absolute (ppm of primary signal)

fig, (ax1) = plt.subplots(figsize=(5, 5), ncols=1)

# Calculate cumulative thicknesses for plotting
cumulative_thickness = np.cumsum(synThk)
depths = np.concatenate([[0], cumulative_thickness])

# Plot resistivity-depth profile
ax1.semilogx(synRes, depths, color='C0', marker='o', label="synth")
ax1.set_xlabel('Resistivity (Ohm-m)')
ax1.set_ylabel('Depth (m)')
ax1.legend()

# times
times = np.logspace(np.log10(1e-4), np.log10(1e-1), 64)

# Loop modeling as a sequence/sum of dipoles
loop_radius = np.sqrt(62500 / np.pi)  # Calculate radius from area
num_dipoles = 200  # Number of dipoles in the loop

TEM_data = np.zeros_like(times, dtype=complex)  # Initialize as complex array

# Loop over each dipole in the loop
for i in range(num_dipoles):
    # Calculate position of current dipole in loop
    angle = 2 * np.pi * i / num_dipoles
    x_dipole = loop_radius * np.cos(angle)
    y_dipole = loop_radius * np.sin(angle)

    # Calculate response of current dipole
    response_dipole = empymod.dipole(
        src=[x_dipole, y_dipole, 90],  # Location of the dipole
        rec=[0, 0, 0],  # Observation point above the ground
        depth=cumulative_thickness,  # Depths of layers
        res=synRes,  # Resistivities of layers
        freqtime=times,  # Times for response
        ab = 11,
        verb=1,  # Verbosity level
    )

    # Add response of current dipole to total response
    TEM_data += response_dipole

# Plot the synthetic data
plt.figure(figsize=(7, 5))
plt.loglog(times, np.abs(TEM_data), label='Observed Data', marker='o')  # Plot absolute value of complex data
plt.xlabel('Time (s)')
plt.ylabel('Response')
plt.legend()
plt.show()

[prisae]

• What is your survey configuration?
• What do you want to model?

Then you can make the corresponding computation. What does your result with pyGIMLi show? Time versus what?

Looking briefly at the code of VMDTimeDomainModelling it looks like that gives you apparent resistivity (rhoapp) - empymod will give you, however, either the electric or magnetic field, depending on your inputs. So you will have to transfer them into apparent resistivity, if you want to compare it.

[halbmy]

the model response (y axis) for the pyGIMLi case is an apparent resistivity. I like this as it best compresses the data. Try to convert the empymod into apparent resistivity and plot both in the same figure. I am looking forward to see the comparison.

[ahsan435]

Dear Sirs @halbmy, @prisae Thank you very much for your time and suggestions. I only need one thing: what is the formula PYGIMLI and EMPYMOD are using to convert the apparent resistivity values into the magnetic field values and vice versa? Thanks you very much once again for your support.

[prisae]

empymod does not compute apparent resistivities. Apparent resistivity depends on your survey configuration.

In pyGIMLi, you can see here: https://github.com/gimli-org/gimli/blob/master/pygimli/physics/em/vmd.py#L259-L274 how it is computed.

For empymod, to reproduce it, you probably need the empymod.loop function, and setting mrec='loop', to get the response for two loops. To have them horizontally co-planar (which I think corresponds to VMD) you will need to set both to [x, y, z, 0, 90], where x, y, z are the coordinates.

However, you cannot have zero horizontal offset in empymod. One solution is to have the receiver as a point in the middle, and construct a loop around it as source. I suspect that the best starting point would be the WalkTEM example: https://empymod.emsig.xyz/en/stable/gallery/tdomain/tem_walktem.html

Otherwise, please explain in more detail what you exactly want to achieve, and how your survey layout looks like.

[ahsan435]

Dear Sir, @prisae @halbmy I hope both are doing well. I tried a lot to convert the PYGIMLI code to the EMPYMOD code but am still not having success, and I also converted the magnetic field into apparent resistivity, but still the results are not matching each other. According to my understanding, there is still an issue in the configuration, but I am not getting it. Please help me regarding this. This is the current progress of the code

Empymod Code:

import numpy as np
import matplotlib.pyplot as plt
import empymod
from scipy.interpolate import interp1d
from pygimli.viewer.mpl import drawModel1D
import math
from scipy.constants import mu_0
from math import pi

# Define the synthetic model parameters
synThk = [300, 100]
synRes = [500, 10, 100]
nlay = len(synRes)  # number of layers
errorabsTEMdata = 1.  # absolute (ppm of primary signal)

fig, (ax1) = plt.subplots(figsize=(5, 5), ncols=1)

drawModel1D(ax1, synThk, synRes, plot='semilogx', color='C0', label="synth")

# Calculate cumulative thicknesses for plotting
cumulative_thickness = np.cumsum(synThk)
depths = np.concatenate([[0], cumulative_thickness])

# Times
time_s = np.logspace(np.log10(1e-4), np.log10(1e-1), 64)

Loop modeling as a sequence/sum of dipoles
loop_radius = np.sqrt(62500 / np.pi)  # Calculate radius from area

num_dipoles = 100  # Number of dipoles in the loop

TEMdata = np.zeros_like(time_s, dtype=complex)  # Initialize as complex array

# Loop over each dipole in the loop
for i in range (num_dipoles):
    # Calculate position of current dipole in loop
    angle = 2 * np.pi * i / num_dipoles
    x_dipole = loop_radius * np.cos(angle)
    y_dipole = loop_radius * np.sin(angle)

    # Calculate response of current dipole
    response_dipole = empymod.dipole(
        src=[x_dipole, y_dipole, 90],  # Location of the dipole
        rec=[0, 0, 90],  # Observation point above the ground
        depth=cumulative_thickness,  # Depths of layers
        res=synRes,  # Resistivities of Layers
        freqtime=time_s,  # Times for response
        verb=1,  # Verbosity level
        ab=66,
    )

    # Add response of current dipole to total response
    TEMdata += response_dipole
    
    #conversion of B into rhoa
B = np.array(TEMdata)

t = np.array(time_s)

current = 1.0

Tx = loop_radius  # Value of A_Tx

mu0 = 4e-7 * pi
rhoa = (current * Tx * mu0 / 30. / B)**(2. / 3.) * mu0 / pi / t

# Calculate the apparent resistivity
apparent_resistivity = np.abs(rhoa) #** 2 / (np.pi * 10000)  # Compute apparent resistivity

# Plot the synthetic data
plt.figure(figsize=(7, 5))
plt.loglog(time_s, apparent_resistivity, label='Apparent Resistivity', marker='o')  Plot apparent resistivity
plt.xlabel('Time (s)')
plt.ylabel('Apparent Resistivity (Ohm.m)')
plt.legend()
plt.show()" 

Pygimli Code:

import numpy as np
import matplotlib.pyplot as plt
import pygimli as pg
from pygimli.physics.em import VMDTimeDomainModelling
from pygimli.viewer.mpl import drawModel1D
from pygimli.frameworks import MarquardtInversion
from scipy.interpolate import interp1d

# Define the synthetic model parameters
synThk = [300, 100]
synRes = [500, 10, 100]
nlay = len(synRes)  # number of layers
errorabsTEM_data = 1.  # absolute (ppm of primary signal)

fig, (ax1) = plt.subplots(figsize=(5, 5), ncols=1)

drawModel1D(ax1, synThk, synRes, plot='semilogx', color='C0', label="synth")

# MaxMin-10/Promys instrument, 1 meter above the ground
times = np.logspace(np.log10(1e-4), np.log10(1e-1), 64)

# Initialize VMDTimeDomainModelling
TEM_modeling = VMDTimeDomainModelling(times=times, nLayers=nlay, txArea=62500, verbose=True)

# Extract data
TEM_data = TEM_modeling(synThk + synRes)

#Plotting subplots
plt.figure(figsize=(7, 5))
plt.loglog(times, TEM_data, label='Observed Data', marker='o')
plt.show()

[halbmy]

can you please enclose your code lines in two lines with triple backcommas?
```
code
```

[ahsan435]

Dear Sir, Please check what you mean like this?

[ahsan435]

Dear Sir @halbmy @prisae this is another Empymod configuration result; please have a look

import numpy as np
import matplotlib.pyplot as plt
import empymod as epm
import math
from scipy.constants import mu_0
from math import pi

# Define the synthetic model parameters
synThk = [300, 100]
synRes = [500, 10, 100]
nlay = len(synRes)  # number of layers
errorabsTEM_data = 1.  # absolute (ppm of primary signal)

# Times
time_s = np.logspace(np.log10(1e-4), np.log10(1e-1), 64)

# Set up dipole coordinates for loop
src = [62500, 0, 0]  # Source coordinates
rec = [0, 0, 62500]  # Receiver coordinates
ab = 66  # Number of AB combinations as per comment 1

# Initialize response array with complex numbers
TEMdata = np.zeros_like(time_s, dtype=complex)

# Values provided
# B = np.array(TEM_data)

# t = np.array(times)

# current = 1.0
# A_Tx = 62500  # Value of A_Tx

# rho_a=(((A_Tx*current*mu_0)/(30*B))**2)/(3*4e-7/t)

# Loop over each dipole
for i in range(ab):
    # Compute response for each dipole
    response = epm.dipole(src, rec, synThk, synRes, times, verb=1)

    # Add response to total
    TEMdata += response

B = np.array(TEMdata)

t = np.array(time_s)

current = 1.0

Tx = 62500  # Value of A_Tx

mu0 = 4e-7 * pi
rhoa = (current * Tx * mu0 / 30. / B)**(2. / 3.) * mu0 / pi / t

#Plotting
plt.figure(figsize=(7, 5))
plt.loglog(times, np.abs(rhoa), label='Observed Data', marker='o')  # Plot absolute value of complex numbers
plt.show()

[prisae]

Dear Sir, Please check what you mean like this?

I did it for you, but in the future, you need to use three backticks ```, not """

[ahsan435]

Dear Sir Thank you very much, and sorry for my mistake. Next time, I will be more careful.

[prisae]

Hi @ahsan435,

I gave it a quick go. However, later they diverge, and I do not know the reason for it now. To understand this difference we would have to look at the underlying assumption of the formulation used in pyGIMLi (maybe @halbmy has an idea).

import empymod
import numpy as np
import matplotlib.pyplot as plt
from pygimli.physics.em import VMDTimeDomainModelling

# Model and survey parameters
times = np.logspace(-4, -1, 101)
depth = [0, 300, 400]
res = [2e14, 500, 10, 100]
side = 250
area = side*side
radius = np.sqrt(area / np.pi)

# === pyGIMLi ===

tem = VMDTimeDomainModelling(times=times, nLayers=len(depth)-1, txArea=area, verbose=True)(list(np.diff(depth)) + res[1:])

# === empymod ===

# The 'loop'-function multiplies the frequency-domain result
# with iwu, which is equivalent to uH => B
em = empymod.loop(
    src=[radius, 0, 0, 0, -90], # Vertical Magn. Dipole VMD at radius dist.
    rec=[0, 0, 0, 90, 0],       # El. y-directed dipole at center
    depth=depth,
    res=res,
    freqtime=times,
    mrec=False,                 # Electric receiver
    signal=-1,                  # for dB/dt
    ftarg={'dlf': 'key_101_CosSin_2012'},  # Works better in this case
)

# Rho-App calculation following 
# https://github.com/gimli-org/gimli/blob/master/pygimli/physics/em/vmd.py#L259-L274
Arec = 1
mu0 = 4e-7 * np.pi
rhoa = radius**(4/3) * Arec**(2/3) * mu0**(5/3) / (20**(2/3) * np.pi**(1/3))
rhoa /= times**(5/3) * (em * 2 * np.pi * radius)**(2/3)

# === Figure ===
plt.figure()
plt.title('Comparison')
plt.loglog(times, rhoa, '-', label='empymod')
plt.loglog(times, tem, '--', label='pyGIMLi')
plt.xlabel('Time (s)')
plt.ylabel('Apparent Resistivity')
plt.legend()

[prisae]

If I put the last layer resistivity to a low value, the two responses agree. However, if I put it to a high value, they do not. I think it might be related to this uncommented code @halbmy, could that be?
https://github.com/gimli-org/gimli/blob/master/pygimli/physics/em/vmd.py#L132-L134

[prisae]

Here I varied the resistivity of the last layer from 1 Ohm.m to 10'000 Ohm.m. For low resistivities, the results agree, for high resistivities, they don't. pyGIMLi seems to always want to go down at the end, no matter how resistive the lower halfspace is.

[prisae]

OK, I think I know what it is. It is just the point where the used Hankel or Fourier transform breaks down. If I model it for larger times up to 1e5 s, empymod also goes down, but it does get closer to the true resistivity of the lowest layer first. So I think in this example the transform of empymod simply fails a bit later then the one of pyGIMLi, that is all.

So now you should have everything you need to know @ahsan435 ;-)

[ahsan435]

Dear Sir, I am incredibly grateful for your help. It means more to me than I can express.

[prisae]

Out of interest, and only if you are willing to share: What are you working on @ahsan435? What is the reason to want to reproduce the pyGIMLi result with empymod? (It is always good to compare solution, so nothing wrong with it.)

[ahsan435]

Dear Sir, @prisae I am presently engaged in the geologically and structurally constrained inversion of TDEM data to enhance subsurface characterization. Given my novice proficiency in coding, I intend to commence with the basic inversion of TDEM using PYGIMLI. Subsequently, I aim to validate my code outcomes using the EMPYMOD library as an initial step. Your assistance has facilitated this initial phase, enabling me to progress to the subsequent stage, the LCI of TDEM. Following this, I plan to incorporate geological and structural data into the code.

[prisae]

Feel free to join the Mattermost Chat of Software Underground for an easy way to communicate. Many EM folks and code developers are there (from pyGIMLi, SimPEG, etc).

Very soon a paper will come out by Lukas from TU Vienna, and he uses empymod for inversions of FastTem data - his code will be available within days to weeks I believe. Also here with us at TUD María uses currently empymod for LCI inversions together with pyGIMLi - it might be interesting for you to follow her work.

[ahsan435]

Dear Sir, @prisae Thank you sincerely for your assistance and for sharing this valuable information with me.