Hi,
It may be a couple of small bugs in the receivers for the H and B-fields in the electromagnetics module:
-
SimPEG.electromagnetics.frequency_domain.receivers.PointMagneticFluxDensitySecondary gives the same real part (in-phase) as .receivers.PointMagneticFluxDensity.
-
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
from scipy.integrate import quad
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
as separate receivers.
receiver_location = np.array([0.0, 0.0, hRx])
receiver_orientation = “z” # “x”, “y” or “z”
data_type = “field” # “secondary”, “total” or “ppm”
Hz_Rx_low = fdem.receivers.PointMagneticFluxDensitySecondary(receiver_location, orientation=receiver_orientation, component=“both”, data_type= “field”)
Hz_Rx_low_tot = fdem.receivers.PointMagneticFluxDensity(receiver_location, orientation=receiver_orientation, component=“both”, data_type= “field” )
receiver_list = [Hz_Rx_low]
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(
fdem.sources.CircularLoop(receiver_list=receiver_list, frequency=freq, location=source_location, orientation=source_orientation, radius = a , current= I))
Define a 1D FDEM survey
survey = fdem.survey.Survey(source_list)
Layer thicknesses
thicknesses = np.array([20, 40])
n_layer = len(thicknesses) + 1
In SimPEG, the Cole-Cole model is used to define a frequency-dependent
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)
Here, we let the infinite conductivity be the model. As a result, we only
need to define the mapping for this parameter. All other parameters used
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)
#######################################################################
Numerical integration and exact solution
-------------------------------------------------
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 = ( mu1e1*(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 = ( mu1e1*(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/( (k12)(a3) ) )( 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()
#plt.savefig(‘exact_num_int_SimPEG_Won_loop_a%.2f_rho%.1f_dr%.3f.pdf’%(a,rho,dr), format=‘pdf’)