Target in a solenoid

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.

Install packages

In [1]:
# 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

model parameters

Here, we define a model consisting of a uniform background, and a target

In [2]:
# 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
In [3]:
# semi-axes of the target
radius_a = 40  # horizontal axis
radius_c = 20  # vertical axis 
In [4]:
# frequencies at which to simulate
freqs = np.logspace(0, 4, 8)
[1.00000000e+00 3.72759372e+00 1.38949549e+01 5.17947468e+01
 1.93069773e+02 7.19685673e+02 2.68269580e+03 1.00000000e+04]
In [5]:
# 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)

    'Minimum skin depth (in target): {:.2e} m '.format(skin_depth(sigma_target, freqs.max()))
    '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 

Set up a mesh

In [6]:
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
In [7]:
# plot the mesh
fig, ax = plt.subplots(1, 1)
# ax.set_xlim([-1., 1500.])  # uncomment to zoom in
# ax.set_ylim([-1000., 1000.])
<matplotlib.axes._subplots.AxesSubplot at 0x12b3b24e0>
In [8]:
print("total number of cells: {}".format(mesh.nC))
total number of cells: 105340
In [9]:
# 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

Plot the model

In [10]:
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].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].set_title('Permeability Model')

[a.set_xlim(xlim) for a in ax]
[a.set_ylim(ylim) for a in ax]


set up the source

In [11]:
# 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])]
In [12]:
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.
In [13]:
# plot the location of the source
fig, ax = plt.subplots(1,1, figsize=(5, 6))
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)

set up the forward simulation

Including the survey and problem

In [14]:
rxList = [] # see FDEM.Rx for receiver types
srcList = [FDEM.Src.RawVec_e(rxList, f, src_vec) for f in freqs]
In [15]:
# 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,  
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. 

Run the simulation

In [16]:
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

Plot the magnetic fields

In [17]:
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(
        vType='F', view='vec', 
        range_x=xlim*1.1, range_y=ylim*1.1,
                'norm': LogNorm(), 
                'cmap': plt.get_cmap('viridis')
        streamOpts={'color': 'w'}, mirror=True, ax=ax, 
        clim=[max_field/cb_range, max_field]
    )[0], ax=ax)

    return ax
def plot_b(
    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],
    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
        'B {reim}, {freq:10.2f} Hz'.format(
    return ax

def plot_difference(
    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]
    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
        'B {reim}, {freq:10.2f} Hz'.format(
    return ax
In [18]:
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)
        return plot_b(fields if model == "with target" else fields_background, freq_ind, reim, xlim, ylim)

    freq_ind=ipywidgets.IntSlider(min=0, max=len(freqs)-1, value=0), 
    reim=ipywidgets.ToggleButtons(options=['real', 'imag']), 
            'background', # wholespace
            'with target',  # wholespace with target
            'difference'  # response due to target (target - background)
<function __main__.plot_fields(freq_ind, reim, model)>
In [ ]: