In this example, we set up a solenoid source and position a target, which is a conductive, permeable ellipsoid in the center. We simulate in the frequency domain.
# Basic python packages
import numpy as np
import matplotlib.pyplot as plt
from scipy.constants import mu_0
import ipywidgets
from matplotlib.colors import LogNorm
# Modules of SimPEG we will use for forward modelling
from SimPEG import Mesh, Utils, Maps
from SimPEG.EM import FDEM
from pymatsolver import Pardiso as Solver
%matplotlib inline
Here, we define a model consisting of a uniform background, and a target
# conductivities (S/m) of the wholespace and target
sigma_wholespace = 1e-3
sigma_target = 1
# relative permeability of target and wholespace
mur_wholespace = 1
mur_target = 1
# semi-axes of the target
radius_a = 40 # horizontal axis
radius_c = 20 # vertical axis
# frequencies at which to simulate
freqs = np.logspace(0, 4, 8)
print(freqs)
[1.00000000e+00 3.72759372e+00 1.38949549e+01 5.17947468e+01 1.93069773e+02 7.19685673e+02 2.68269580e+03 1.00000000e+04]
# print the min and max skin depths to make sure mesh is fine enough and
# extends far enough
def skin_depth(sigma, f):
return 500./np.sqrt(sigma * f)
print(
'Minimum skin depth (in target): {:.2e} m '.format(skin_depth(sigma_target, freqs.max()))
)
print(
'Maximum skin depth (in background): {:.2e} m '.format(skin_depth(sigma_wholespace, freqs.min()))
)
Minimum skin depth (in target): 5.00e+00 m Maximum skin depth (in background): 1.58e+04 m
cs = 2.5 # core cell size
ncx = np.ceil(500/cs) # number of core cells in the x direction
ncz = np.ceil(1000/cs) # number of core cells in the z direction
pf = 1.3 # padding factor (how much we expand the cells)
npadx = 29 # number of padding cells in the x direction (need to be > max skin depth)
npadz = 30 # number of padding cells in z-direction
# vectors of cell sizes in each direction
hx = Utils.meshTensor([(cs, ncx), (cs, npadx, pf)])
hz = Utils.meshTensor([(cs, npadz, -pf), (cs, ncz), (cs, npadz, pf)])
# create a cylindrical mesh (has geometry, and differential operators)
mesh = Mesh.CylMesh(
[hx, 1, hz], # cylindrically symmetric mesh
x0='00C' # origin at the axis of symmetry in 0, centered in z
)
# plot the mesh
fig, ax = plt.subplots(1, 1)
mesh.plotGrid(ax=ax)
# ax.set_xlim([-1., 1500.]) # uncomment to zoom in
# ax.set_ylim([-1000., 1000.])
<matplotlib.axes._subplots.AxesSubplot at 0x12b3b24e0>
print("total number of cells: {}".format(mesh.nC))
total number of cells: 105340
# put the model on the mesh
sigma_background = sigma_wholespace*np.ones(mesh.nC)
mur_background = mur_wholespace*np.ones(mesh.nC)
# ellipsoid
sigma = sigma_background.copy()
mur = mur_background.copy()
# indicies in the mesh where the ellipsoid is
inds_ellipsoid = (mesh.gridCC[:,0]**2/radius_a**2 + mesh.gridCC[:,2]**2/radius_c**2) <= 1.
# models with target
sigma[inds_ellipsoid] = sigma_target
mur[inds_ellipsoid] = mur_target
xlim = np.r_[-100., 100.]
ylim = np.r_[-100., 100]
# Plot a cross section of the conductivity model
fig, ax = plt.subplots(1,2, figsize=(10, 4))
# conductivity model
cb = plt.colorbar(mesh.plotImage(np.log10(sigma), ax=ax[0], mirror=True)[0], ax=ax[0])
# plot formatting and titles
cb.set_label('$\log_{10}\sigma$', fontsize=13)
ax[0].axis('equal')
ax[0].set_title('Conductivity Model')
# permeability model
cb = plt.colorbar(mesh.plotImage(mur, ax=ax[1], mirror=True)[0], ax=ax[1])
# plot formatting and titles
cb.set_label('$\mu_r$', fontsize=13)
ax[1].axis('equal')
ax[1].set_title('Permeability Model')
[a.set_xlim(xlim) for a in ax]
[a.set_ylim(ylim) for a in ax]
plt.tight_layout()
# find the indices where we want to put the solenoid - it is on edges
# (electric field is defined on cell edges for EB formulation)
solenoid_x = 500 # radius of the solenoid
solenoid_zlims = np.r_[-2000., 2000.] # 1km long solenoid
solenoid_z = mesh.vectorNz[(mesh.vectorNz>=solenoid_zlims[0]) & (mesh.vectorNz<=solenoid_zlims[1])]
solenoid_inds = Utils.closestPoints(mesh, Utils.ndgrid([np.r_[solenoid_x], np.r_[0], solenoid_z]), 'Ey')
src_vec = np.zeros(mesh.nEy)
src_vec[solenoid_inds] = 1.
# plot the location of the source
fig, ax = plt.subplots(1,1, figsize=(5, 6))
mesh.plotGrid(ax=ax)
ax.plot(mesh.gridEy[solenoid_inds,0], mesh.gridEy[solenoid_inds, 2], 'ro')
ax.set_xlim([-1, 3e3])
ax.set_ylim([-3000, 3000])
(-3000, 3000)
Including the survey and problem
rxList = [] # see FDEM.Rx for receiver types
srcList = [FDEM.Src.RawVec_e(rxList, f, src_vec) for f in freqs]
# define a problem - the statement of which discrete pde system we want to solve
wires = Maps.Wires(('sigma', mesh.nC), ('mu', mesh.nC)) # mapping to work with a vector m = [sigma, mu] as the model
prob = FDEM.Problem3D_e(mesh, sigmaMap=wires.sigma, muMap=wires.mu)
prob.solver = Solver
survey = FDEM.Survey(srcList)
# tell the problem and survey about each other - so the RHS can be constructed
# for the problem and the resulting fields and fluxes can be sampled by the receiver.
prob.pair(survey)
%%time
fields_background = prob.fields(np.r_[sigma_background, mu_0*mur_background])
fields = prob.fields(np.r_[sigma, mu_0*mur])
CPU times: user 52.3 s, sys: 9.88 s, total: 1min 2s Wall time: 51 s
def streamplot_b(field_to_plot, ax):
# Plot the magnetic flux
max_field = np.abs(field_to_plot).max() #use to set colorbar limits
cb_range = 5e2 # dynamic range of colorbar
cb = plt.colorbar(mesh.plotImage(
field_to_plot,
vType='F', view='vec',
range_x=xlim*1.1, range_y=ylim*1.1,
pcolorOpts={
'norm': LogNorm(),
'cmap': plt.get_cmap('viridis')
},
streamOpts={'color': 'w'}, mirror=True, ax=ax,
clim=[max_field/cb_range, max_field]
)[0], ax=ax)
ax.set_xlim(xlim)
ax.set_ylim(ylim)
return ax
def plot_b(
field,
freq_ind=0, # which frequency would you like to look at?
reim='real', # real or imaginary part
xlim=np.r_[-1000, 1000],
ylim=np.r_[-1000, 1000],
ax=None
):
if ax is None:
fig, ax = plt.subplots(1,1, figsize=(6,5))
field_to_plot = getattr(field[srcList[freq_ind], 'b'], reim)
ax = streamplot_b(field_to_plot, ax)
cb.set_label('|B {reim}|'.format(reim=reim))
# give it a title
ax.set_title(
'B {reim}, {freq:10.2f} Hz'.format(
reim=reim,
freq=freqs[freq_ind]
)
)
plt.show()
return ax
def plot_difference(
field,
field_background,
freq_ind=0, # which frequency would you like to look at?
reim='real', # real or imaginary part
ax=None,
xlim=np.r_[-1000, 1000],
ylim=np.r_[-1000, 1000]
):
if ax is None:
fig, ax = plt.subplots(1,1, figsize=(6,5))
field_to_plot = getattr(field[srcList[freq_ind], 'b'], reim)
background = getattr(field_background[srcList[freq_ind], 'b'], reim)
ax = streamplot_b(field_to_plot - background, ax)
cb.set_label('|B {reim}|'.format(reim=reim))
# give it a title
ax.set_title(
'B {reim}, {freq:10.2f} Hz'.format(
reim=reim,
freq=freqs[freq_ind]
)
)
plt.show()
return ax
xlim=np.r_[-200., 200.]
ylim=np.r_[-200., 200.]
def plot_fields(freq_ind, reim, model):
if model == "difference":
return plot_difference(fields, fields_background, freq_ind, reim=reim, xlim=xlim, ylim=ylim)
else:
return plot_b(fields if model == "with target" else fields_background, freq_ind, reim, xlim, ylim)
ipywidgets.interact(
plot_fields,
freq_ind=ipywidgets.IntSlider(min=0, max=len(freqs)-1, value=0),
reim=ipywidgets.ToggleButtons(options=['real', 'imag']),
model=ipywidgets.ToggleButtons(
options=[
'background', # wholespace
'with target', # wholespace with target
'difference' # response due to target (target - background)
]
),
)
interactive(children=(IntSlider(value=0, description='freq_ind', max=7), ToggleButtons(description='reim', opt…
<function __main__.plot_fields(freq_ind, reim, model)>