# Possible bug in Electromagnetics frequency domain module

Hi,

It may be a couple of small bugs in the receivers for the H and B-fields in the electromagnetics module:

2. Also, it may be that there is missing a scaling since the B-field is the same as the H-field, this appears when …PointMagneticFieldSecondary and …PointMagneticField are used instead of PointMagneticFluxDensitySecondary and PointMagneticFluxDensity. It seems that it is just mu that is missing in B = mu H.

These artifacts appear when using the example for 1D

Thanks in advance for help, or a verification of the bug.

Best regards
Erlend

The code that reproduce these effects are pasted below. The example is a 50 m radius loop lying on a unform halfspace of 100 Ohm m, which has an exact solution that also are included in this code and compared to the SimPEG results. The particular example is from Ward and Hohmann, 1988.

import numpy as np
from matplotlib import pyplot as plt
import scipy.special as special
from SimPEG import maps
import SimPEG.electromagnetics.frequency_domain as fdem
from SimPEG.electromagnetics.utils.em1d_utils import ColeCole
mu0 = 1.25663706212e-06 # 4np.pi1e-07
e0 = 8.85e-12

def pos(data):
“”“Return positive data; set negative data to NaN.”“”
return np.array([x if x > 0 else np.nan for x in data])
def neg(data):
“”“Return -negative data; set positive data to NaN.”“”
return np.array([-x if x < 0 else np.nan for x in data])

a = 50 # Loop source radius
I = 1 # Loop current

hTx = 0.02 # Tx height
hRx = 0.01 # Rx …

Kappa = 0

z = hRx

# Frequencies being observed in Hz

f = np.logspace(0, 6, 41)

# Define a list of receivers. The real and imaginary components are defined

receiver_orientation = “z” # “x”, “y” or “z”
data_type = “field” # “secondary”, “total” or “ppm”

# Define a source list. A source must defined for each frequency.

source_location = np.array([0.0, 0.0, hTx])
source_orientation = “z” # “x”, “y” or “z”
#moment = 1.0 # dipole moment

source_list =
for freq in f:
source_list.append(

# Define a 1D FDEM survey

survey = fdem.survey.Survey(source_list)

# Layer thicknesses

thicknesses = np.array([20, 40])
n_layer = len(thicknesses) + 1

# electrical conductivity when the Earth is chargeable.

rho = 100 # infinite resisivity

e1 = e0*1
sigma = 1/rho
sigma1 = sigma

mu = mu0*(1+Kappa)
mu1 = mu

sigma_model = sigma * np.ones(n_layer)

# to define physical properties will be fixed when creating the simulation.

model_mapping = maps.IdentityMap(nP=n_layer)

# Response for conductive Earth

simulation = fdem.Simulation1DLayered(
survey=survey, thicknesses=thicknesses, sigmaMap=model_mapping
)

dpred = simulation.dpred(sigma_model)

#######################################################################

W = 2np.pif

# Numerical integration fuctions:

def integrand_Rye_real(lam, w, z, Kappa, rho, a, mu0=mu0, e0=e0):
mu1 = mu0*(1+Kappa)
e1 = e01
sigma1 = 1/rho
k1 = ( mu1
e1*(w**2) - 1jmu1sigma1w)**0.5 # sigma >> e_1w

# k1 = (- 1jwmu_1sigma_1)**2 # sigma >> e_1w

``````k0 = (mu0*e0*(w**2) )**0.5

u1 = (lam**2 - k1**2)**0.5
u0 = (lam**2 - k0**2)**0.5

Z0 = (-1j*w*mu0)/u0
Z1 = (-1j*w*mu1)/u1

int_real = np.real((   ( np.exp(-u0*z)/u0 )*( (Z1*lam) / (Z0 + Z1) ) - np.exp(-lam*z)/2  )*special.jv(1,a*lam)*lam)

return int_real
``````

def loopint_Rye_real(w, z, Kappa):
return quad(integrand_Rye_real, 0, np.inf, args = (w, z, Kappa, rho, a, mu0, e0) )

def integrand_Rye_imag(lam, w, z, Kappa, rho, a, mu0=mu0, e0=e0):
mu1 = mu0*(1+Kappa)
e1 = e01
sigma1 = 1/rho
k1 = ( mu1
e1*(w**2) - 1jmu1sigma1w)**0.5 # sigma >> e_1w

# k1 = (- 1jwmu_1sigma_1)**2 # sigma >> e_1w

``````k0 = (mu0*e0*(w**2) )**0.5

u1 = (lam**2 - k1**2)**0.5
u0 = (lam**2 - k0**2)**0.5

Z0 = (-1j*w*mu0)/u0
Z1 = (-1j*w*mu1)/u1

int_imag = np.imag((   ( np.exp(-u0*z)/u0 )*( (Z1*lam) / (Z0 + Z1) ) - np.exp(-lam*z)/2  )*special.jv(1,a*lam)*lam)

return int_imag
``````

def loopint_Rye_imag(w, z, Kappa):
return quad(integrand_Rye_imag, 0, np.inf, args = (w, z, Kappa, rho, a, mu0, e0) )

Hz_surface_exact_real = np.full((1, len(f)), np.nan)
Hz_surface_exact_imag = np.full((1, len(f)), np.nan)
Hz_numeric_int_real = np.full((1, len(f)), np.nan)
Hz_numeric_int_imag = np.full((1, len(f)), np.nan)

wdx = 0
for w in W:
# Analytisk
k1 = ( mu1e1(w2) - 1jmu1sigma1w)0.5
Hz_surface_exact_real[0,wdx] = np.real( (-I/( (k1
2)
(a
3) ) )( 3 - (3 + 31jk1a - (k12)*(a2))np.exp(-1jk1a) ) )
Hz_surface_exact_imag[0,wdx] = np.imag( (-I/( (k1**2)
(a3) ) )( 3 - (3 + 31jk1a - (k12)(a**2))np.exp(-1jk1a) ) )

# Numerisk integrasjon

``````results_Rye_real = loopint_Rye_real(w, z, Kappa)
results_Rye_imag = loopint_Rye_imag(w, z, Kappa)

Hz_numeric_int_real[0,wdx] = I*a*results_Rye_real[0] + (((a**2)*I)/( 2*((a**2 + z**2 )**(3/2)) ))
Hz_numeric_int_imag[0,wdx] = I*a*results_Rye_imag[0]

wdx = wdx + 1
``````

#######################################################################

# -------------------------------------------------

dr=1

fig= plt.figure(1)
plt.clf()
ax1 = fig.add_axes([0.15, 0.1, 0.8, 0.85])

plt.title(“SimPEG. res=%.1f \$\Omega\$m. a=%.3f m” %(rho,a) )

plt.plot(f, pos( Hz_surface_exact_real[0,:]), ‘-k’, label=‘Real exact > 0’, linewidth = 3)
plt.plot(f, pos( Hz_numeric_int_real[0,:]), ‘+k’, label=‘Real numeric integration > 0’ , linewidth = 3)
plt.plot(f, pos((dpred[0::2])), “r–”, label=‘Real SimPEG > 0’, lw=3)

plt.plot(f, neg( Hz_surface_exact_real[0,:]), ‘-m’, label=‘Real exact < 0’, linewidth = 3)
plt.plot(f, neg( Hz_numeric_int_real[0,:]), ‘+’, color = ‘tab:orange’, label=‘Real numeric integration < 0’, linewidth = 3)
plt.plot(f, neg((dpred[0::2])), “g–”, label=‘Real SimPEG < 0’, lw=3)

plt.plot(f, neg( Hz_surface_exact_imag[0,:]), ‘-b’, label=‘Imag exact < 0’, linewidth = 3)
plt.plot(f, neg( Hz_numeric_int_imag[0,:]), ‘+c’, label=‘Imag numeric integration < 0’, linewidth = 3)
plt.plot(f, neg((dpred[1::2])), “y–”, label=‘Imag SimPEG < 0’, lw=3)

plt.plot(f, pos( Hz_surface_exact_imag[0,:]), ‘-g’, label=‘Imag exact > 0’ , linewidth = 3)
plt.plot(f, pos( Hz_numeric_int_imag[0,:]), ‘+y’, label=‘Imag numeric integration > 0’, linewidth = 3)
plt.plot(f, pos((dpred[1::2])), ‘–’, color = ‘tab:orange’, label=‘Imag SimPEG > 0’, lw=3)

plt.yscale(‘log’)
plt.xscale(‘log’)
plt.xlabel(“f (Hz)”)
plt.ylabel(“Hz (A/m)”)
plt.legend(frameon = False)
plt.show()
#plt.savefig(‘exact_num_int_SimPEG_Won_loop_a%.2f_rho%.1f_dr%.3f.pdf’%(a,rho,dr), format=‘pdf’)

Hi again,

Sorry, the pasted code did not appear very nice in the last post, and as a new user I am not allowed to upload files. I paste the working code without so many comments below.

To summarize, the receivers for B-secondary, B-primary, and H-secondary gives exactly the same results. These are:

To reproduce the bugs, just uncomment the different receivers in one of the four receiver_lists:

import numpy as np
from matplotlib import pyplot as plt
import scipy.special as special
from SimPEG import maps
import SimPEG.electromagnetics.frequency_domain as fdem
from SimPEG.electromagnetics.utils.em1d_utils import ColeCole
mu0 = 1.25663706212e-06 # 4np.pi1e-07
e0 = 8.85e-12

def pos(data):
“”“Return positive data; set negative data to NaN.”“”
return np.array([x if x > 0 else np.nan for x in data])
def neg(data):
“”“Return -negative data; set positive data to NaN.”“”
return np.array([-x if x < 0 else np.nan for x in data])

a = 50 # Loop source radius
I = 1 # Loop current

hTx = 0.02 # Tx height
hRx = 0.01 # Rx …

Kappa = 0

z = hRx

f = np.logspace(0, 6, 41)

receiver_orientation = “z” # “x”, “y” or “z”
data_type = “field” # “secondary”, “total” or “ppm”

source_location = np.array([0.0, 0.0, hTx])
source_orientation = “z” # “x”, “y” or “z”

source_list =
for freq in f:
source_list.append(

survey = fdem.survey.Survey(source_list)

thicknesses = np.array([20, 40])
n_layer = len(thicknesses) + 1
rho = 100 # infinite resisivity

e1 = e0*1
sigma = 1/rho
sigma1 = sigma

mu = mu0*(1+Kappa)
mu1 = mu

sigma_model = sigma * np.ones(n_layer)

model_mapping = maps.IdentityMap(nP=n_layer)

simulation = fdem.Simulation1DLayered(
survey=survey, thicknesses=thicknesses, sigmaMap=model_mapping
)

dpred = simulation.dpred(sigma_model)

#######################################################################

W = 2np.pif

# Numerical integration fuctions:

def integrand_Rye_real(lam, w, z, Kappa, rho, a, mu0=mu0, e0=e0):
mu1 = mu0*(1+Kappa)
e1 = e01
sigma1 = 1/rho
k1 = ( mu1
e1*(w2) - 1jmu1sigma1w)**0.5 # sigma >> e_1w
k0 = (mu0e0(w
2) )**0.5

``````u1 = (lam**2 - k1**2)**0.5
u0 = (lam**2 - k0**2)**0.5

Z0 = (-1j*w*mu0)/u0
Z1 = (-1j*w*mu1)/u1

int_real = np.real((   ( np.exp(-u0*z)/u0 )*( (Z1*lam) / (Z0 + Z1) ) - np.exp(-lam*z)/2  )*special.jv(1,a*lam)*lam)

return int_real
``````

def loopint_Rye_real(w, z, Kappa):
return quad(integrand_Rye_real, 0, np.inf, args = (w, z, Kappa, rho, a, mu0, e0) )

def integrand_Rye_imag(lam, w, z, Kappa, rho, a, mu0=mu0, e0=e0):
mu1 = mu0*(1+Kappa)
e1 = e01
sigma1 = 1/rho
k1 = ( mu1
e1*(w2) - 1jmu1sigma1w)**0.5 # sigma >> e_1w
k0 = (mu0e0(w
2) )**0.5

``````u1 = (lam**2 - k1**2)**0.5
u0 = (lam**2 - k0**2)**0.5

Z0 = (-1j*w*mu0)/u0
Z1 = (-1j*w*mu1)/u1

int_imag = np.imag((   ( np.exp(-u0*z)/u0 )*( (Z1*lam) / (Z0 + Z1) ) - np.exp(-lam*z)/2  )*special.jv(1,a*lam)*lam)

return int_imag
``````

def loopint_Rye_imag(w, z, Kappa):
return quad(integrand_Rye_imag, 0, np.inf, args = (w, z, Kappa, rho, a, mu0, e0) )

Hz_surface_exact_real = np.full((1, len(f)), np.nan)
Hz_surface_exact_imag = np.full((1, len(f)), np.nan)
Hz_numeric_int_real = np.full((1, len(f)), np.nan)
Hz_numeric_int_imag = np.full((1, len(f)), np.nan)

wdx = 0
for w in W:
# Analytisk
k1 = ( mu1e1(w2) - 1jmu1sigma1w)0.5
Hz_surface_exact_real[0,wdx] = np.real( (-I/( (k1
2)
(a
3) ) )( 3 - (3 + 31jk1a - (k12)*(a2))np.exp(-1jk1a) ) )
Hz_surface_exact_imag[0,wdx] = np.imag( (-I/( (k1**2)
(a3) ) )( 3 - (3 + 31jk1a - (k12)(a**2))np.exp(-1jk1a) ) )

# Numerisk integrasjon

``````results_Rye_real = loopint_Rye_real(w, z, Kappa)
results_Rye_imag = loopint_Rye_imag(w, z, Kappa)

Hz_numeric_int_real[0,wdx] = I*a*results_Rye_real[0] + (((a**2)*I)/( 2*((a**2 + z**2 )**(3/2)) ))
Hz_numeric_int_imag[0,wdx] = I*a*results_Rye_imag[0]

wdx = wdx + 1
``````

#######################################################################

## Plotting Results

dr=1

fig= plt.figure(1)
plt.clf()
ax1 = fig.add_axes([0.15, 0.1, 0.8, 0.85])

plt.title(“SimPEG. res=%.1f \$\Omega\$m. a=%.3f m” %(rho,a) )

plt.plot(f, pos( Hz_surface_exact_real[0,:]), ‘-k’, label=‘Real exact > 0’, linewidth = 3)
plt.plot(f, pos( Hz_numeric_int_real[0,:]), ‘+k’, label=‘Real numeric integration > 0’ , linewidth = 3)
plt.plot(f, pos((dpred[0::2])), “r–”, label=‘Real SimPEG > 0’, lw=3)

plt.plot(f, neg( Hz_surface_exact_real[0,:]), ‘-m’, label=‘Real exact < 0’, linewidth = 3)
plt.plot(f, neg( Hz_numeric_int_real[0,:]), ‘+’, color = ‘tab:orange’, label=‘Real numeric integration < 0’, linewidth = 3)
plt.plot(f, neg((dpred[0::2])), “g–”, label=‘Real SimPEG < 0’, lw=3)

plt.plot(f, neg( Hz_surface_exact_imag[0,:]), ‘-b’, label=‘Imag exact < 0’, linewidth = 3)
plt.plot(f, neg( Hz_numeric_int_imag[0,:]), ‘+c’, label=‘Imag numeric integration < 0’, linewidth = 3)
plt.plot(f, neg((dpred[1::2])), “y–”, label=‘Imag SimPEG < 0’, lw=3)

plt.plot(f, pos( Hz_surface_exact_imag[0,:]), ‘-g’, label=‘Imag exact > 0’ , linewidth = 3)
plt.plot(f, pos( Hz_numeric_int_imag[0,:]), ‘+y’, label=‘Imag numeric integration > 0’, linewidth = 3)
plt.plot(f, pos((dpred[1::2])), ‘–’, color = ‘tab:orange’, label=‘Imag SimPEG > 0’, lw=3)

plt.yscale(‘log’)
plt.xscale(‘log’)
plt.xlabel(“f (Hz)”)
plt.ylabel(“Hz (A/m)”)
plt.legend(frameon = False)
plt.show()

Best regards,
Erlend