Electromagnetic forward Modelling for Plate Models in Cylindrical Coordinates

Here we use SimPEG and simpegEM to forward-model the electromagnetic response of a plate. We use two plate geometries, one a circle, and one a disk with a hole in it. The plates are defined to have equal volumes. For a source, we use a vertical magnetic dipole and for receivers we measure both the in-phase (or real) and quadrature (or imaginary) components. We take advantage of the cylindrical symmetry of this problem and treat it as a 2D cylindrically symmetric one.

The coordinate system is:

z: along the axis of the well
x: radial distance from the well

Using the iPython Notebook

This is the iPython Notebook, an interactive coding and computation environment

To use the notebook:

"Shift + Enter" runs the code within the cell (so does the forward arrow button near the top of the document)
You can alter variables and re-run cells
If you want to start with a clean slate, restart the Kernel either by going to the top, clicking on Kernel: Restart, or by "esc + 00" (if you do this, you will need to re-run Step 0 before running any other cells in the notebook)

Step 0: Import Necessary Packages



In [1]:

    
# import the necessary packages
from SimPEG import *
import simpegEM as EM
%pylab inline









    



Populating the interactive namespace from numpy and matplotlib

Step 1: Model Parameters

Define the size and conductivity of the plate and background. Here, we set up two models, one of a circular fracture and one for a disk that has a hole in it (outer radius Rd1, inner radius Rd2). The thickness of each plate is equal and given by the variable t. The radius of the circular plate is chosen so that the volume of the disk plate and the circular plate are equal.



In [2]:

    
# Model Parameters
sigback = 1e-2 # 100 Ohm-m background
sigfrac = 2000. # S/m Frac: conductivity 

t = 5e-3 # fracture thickness

# for the disk distribution
Rd1 = 20. # inner radius
Rd2 = 40. # outer radius (m)

# for the circle disribution
Rc  = np.sqrt(Rd2**2 -Rd1**2) # conserve volume between examples

Step 2: Survey Parameters

Here we specify the survey parameters. We are using a single vertical magnetic dipole transmitter positioned at the origin and multipple receivers. The parameters we are setting are:

freq : frequency of the transmitter
nrx : number of receivers
rxx : x-location of the receivers (assuming that all receivers lie along
rxz : z-locations of the receivers

Here, we also calculate the skin-depth (delta), which is used to estimate how quickly the fields decay. This is important to note this value because we need to make sure that the mesh is large enough that the fields will have decayed by the time they get to the boundary (typically we want the mesh to extend >2 times this value).



In [4]:

    
# Survey Parameters
freq = np.logspace(1, 5, 20)  # frequency 
print freq
# z - along the well

# receivers
nrx = 101.   # number of receivers 
rxx = 100.  # x location of receivers
rxz = np.linspace(-150,150,nrx) # z location of receivers along the well

# transmitter
txz = 0. # transmitter z location

# transmitter moment
txmom = 2000 # Am^2

# Noise Floor
flr = 1e-14 # smallest signal we can measure
perc_cutoff = 3 # smallest percent difference we are confident in 

# calculate the skin depth
delta = 500./np.sqrt(freq*sigback)
print delta # skin depth









    



[  1.00000000e+01   1.62377674e+01   2.63665090e+01   4.28133240e+01
   6.95192796e+01   1.12883789e+02   1.83298071e+02   2.97635144e+02
   4.83293024e+02   7.84759970e+02   1.27427499e+03   2.06913808e+03
   3.35981829e+03   5.45559478e+03   8.85866790e+03   1.43844989e+04
   2.33572147e+04   3.79269019e+04   6.15848211e+04   1.00000000e+05]
[ 1581.13883008  1240.81446142   973.74151995   764.15336633   599.6769731
   470.60248363   369.3099911    289.81969767   227.43889735   178.48494234
   140.06783806   109.91963244    86.2605275     67.693809      53.12339154
    41.68911117    32.71594565    25.67416454    20.1480566     15.8113883 ]

Step 3. Set up Mesh

Now to solve the problem, we need to discretize it. To do this we define a mesh on which our model and fields will live and on which we will solve the problem.



In [5]:

    
# set up a mesh (all units are meters)
csx, ncx, npadx = 5 , 41, 20 # cell size, number of core cells, number of padding cells in the x- direction
csz, ncz, npadz = 5 , 44, 20 # cell size, number of core cells, number of padding cells in the z- direction
hx = Utils.meshTensor([(csx,npadx,-0.7), (csx,ncx), (csx,npadx,1.3)]) # vector of cell widths in the x-direction
hz = Utils.meshTensor([(csz,npadz,-1.3), (csz,ncz), (csz,npadz,1.3)]) # vector of cell widths in the z-direction
mesh = Mesh.CylMesh([hx,1,hz], [0.,0,-hz.sum()/2.]) # define the cylindrical mesh using SimPEG
mesh.plotGrid() # Plot the mesh
print mesh.nC

Step 4: Put Model on the Mesh

Now that we have a mesh, we need to define our model on that mesh. Since the fracture is so thin compared to the size of the cells, we upscale the conductivity, that is, make a coarse approximation of the conductivity on the coarse mesh. To do this, we treat the cells with the plate running through them as anisotropic, so the conductivity varies with direction. We treat the vertical component of the conductivity like a series circuit and the component parallel to the plate as a parallel circuit.



In [6]:

    
# find where in the mesh the circle and disc should be
fracCircle = (mesh.gridCC[:,0]>mesh.hx[0]) & (mesh.gridCC[:,0]<=Rc) & (mesh.gridCC[:,2]<np.max(np.array([csz,t/2.]))) & (mesh.gridCC[:,2]>-np.max(np.array([csz,t/2.])))
fracDisc   = (mesh.gridCC[:,0]>Rd1) & (mesh.gridCC[:,0]<=Rd2) & (mesh.gridCC[:,2]<np.max(np.array([csz,t/2.]))) & (mesh.gridCC[:,2]>-np.max(np.array([csz,t/2.]))) 

# inline functions for series and parallel circuit approximations of the anisotropic conductivity
scaleByCT = lambda sig1,t1,t2: sig1*t1/t2 

Hxyz = Utils.ndgrid(mesh.hx,np.r_[0.],mesh.hz)

# Background conductivity
sigBack = sigback*np.ones((mesh.nC)) 

# conductivity model for the circular frac
sigfracCircle = sigback*np.ones((mesh.nC)) # anisotropic conductivity
sigfracCircle[fracCircle] = scaleByCT(sigfrac,t,Hxyz[fracCircle,2])

# conductivity model for the disc frac
sigfracDisc = sigback*np.ones((mesh.nC)) # anisotropic conductivity
sigfracDisc[fracDisc] = scaleByCT(sigfrac,t,Hxyz[fracDisc,2])



In [66]:

    
# plot them
fs = 32 # fontsize
xlim = [0.,150.]     #xlimits
zlim = [-200., 200.] #ylimits

def plotIt(img,label=None,title=None):
    plt.colorbar(img,label=label)
    plt.xlim(xlim)
    plt.ylim(zlim)
    plt.plot(0.,txz,'d',markersize=12,markeredgecolor='white',markerfacecolor='white')
    plt.plot(rxx*np.ones(rxz.shape),rxz,'o',markersize=4,markeredgecolor='yellow',markerfacecolor='yellow')
    plt.title(title)
    plt.fontsize = fs
    plt.grid(which='both')
    return

# Circle 
plotIt(mesh.plotImage(np.log10(sigfracCircle))[0],'$\sigma$','Circle, Parallel Conductivity: $\sigma$')

# Disc
plotIt(mesh.plotImage(np.log10(sigfracDisc))[0],'$\sigma$','Disc, Parallel Conductivity: $\sigma$')

Set Problem and Survey

Pair the problem and the survey so that we can create fields and collect data



In [67]:

    
# receivers
XYZ = Utils.ndgrid(np.r_[rxx],np.r_[0.],np.r_[rxz])   # define receivers on grid
Rx0 = EM.FDEM.RxFDEM(XYZ, 'bzr')               # create real / in-phase receivers
Rx1 = EM.FDEM.RxFDEM(XYZ, 'bzi')               # create imag / quadrature receivers

# transmitters
txlist = []
for f in freq:
    Tx =  EM.FDEM.TxFDEM(np.r_[0.,0.,txz], 'VMD',f,[Rx0,Rx1]) # create vertical magnetic dipole transmitter
    txlist.append(Tx)
    
# pair the problem and survey
survey = EM.FDEM.SurveyFDEM(txlist) # set survey
prb = EM.FDEM.ProblemFDEM_e(mesh) # set problem
prb.pair(survey) # pair problem and survey

Solve for Data



In [68]:

    
# background - no fracture
dPrimary = survey.dpred(sigBack)*txmom # in-phase
print 'Done Calculating Primary'

# circle model
dCircle = survey.dpred(sigfracCircle)*txmom
print 'Done Calculating Circle'

# circle model
dDisc = survey.dpred(sigfracDisc)*txmom
print 'Done Calculating Disc'









    



Done Calculating Primary
Done Calculating Circle
Done Calculating Disc



In [69]:

    
# Sort Data
dPrim = np.reshape(dPrimary,(freq.size,2*nrx))
dCirc = np.reshape(dCircle,(freq.size,2*nrx))
dDisc = np.reshape(dDisc,(freq.size,2*nrx))

# Sort into In-Phase and Quadrature
dPrimI = dPrim[:,:nrx]
dPrimQ = dPrim[:,nrx:]

dCircI = dCirc[:,:nrx]
dCircQ = dCirc[:,nrx:]

dDiscI = dDisc[:,:nrx]
dDiscQ = dDisc[:,nrx:]

Plot Results

Plot the data for a single frequency



In [70]:

    
freq_ind = 6
print freq[freq_ind]









    



183.298071083



In [71]:

    
plt.figure(figsize=(11,6))

# Plot In-Phase
plt.subplot(121)
plt.plot(rxz,dCircI[freq_ind,:],rxz,dDiscI[freq_ind,:],linewidth=2)
plt.plot(rxz,flr*np.ones(rxz.shape),'--k', rxz,-flr*np.ones(rxz.shape),'--k',linewidth=2)
plt.grid(which='both')
plt.title('In-Phase')
plt.ylabel('Magnetic Field Response (T)')
plt.xlabel('x (m)')
print 'Max In-Phase Responses: Circle %e T , Disc %e T' %(np.max(dCircI[freq_ind,:]), np.max(dDiscI[freq_ind,:]))

# Plot Quadrature
plt.subplot(122)
plt.plot(rxz,dCircQ[freq_ind,:],rxz,dDiscQ[freq_ind,:],linewidth=2)
plt.plot(rxz,flr*np.ones(rxz.shape),'--k', rxz,-flr*np.ones(rxz.shape),'--k',linewidth=2)
plt.grid(which='both')
plt.title('Quadrature')
plt.legend(['Circle','Disc'],loc=2, bbox_to_anchor = (1.0, 1.0),title='Model')
plt.ylabel('Magnetic Field Response (T)')
plt.xlabel('x (m)')
print 'Max Quadrature Responses: Circle %e T , Disc %e T' %(np.max(dCircQ[freq_ind,:]), np.max(dDiscQ[freq_ind,:]))


plt.tight_layout()
plt.show()









    



Max In-Phase Responses: Circle 3.317402e-11 T , Disc 3.473048e-11 T
Max Quadrature Responses: Circle 4.014403e-11 T , Disc 2.264680e-11 T



In [72]:

    
# Difference between circle and disc signals
Idiff = dCircI[freq_ind,:]-dDiscI[freq_ind,:]
Qdiff = dCircQ[freq_ind,:]-dDiscQ[freq_ind,:]

plt.figure(figsize=(11,6))

# Plot In-phase
plt.subplot(121)
plt.plot(rxz,Idiff,linewidth=2)
plt.plot(rxz,flr*np.ones(rxz.shape),'--k', rxz,-flr*np.ones(rxz.shape),'--k',linewidth=2)
plt.grid(which='both')
plt.title('In-Phase')
plt.ylabel('Difference Magnetic Field Response (T)')
plt.xlabel('x (m)')
print 'Max In-Phase Difference: %e T' %(np.max(dCircI-dDiscI))

# Plot Quadrature
plt.subplot(122)
plt.plot(rxz,Qdiff,linewidth=2)
plt.plot(rxz,flr*np.ones(rxz.shape),'--k', rxz,-flr*np.ones(rxz.shape),'--k',linewidth=2)
plt.grid(which='both')
plt.title('Quadrature')
plt.legend(['Circle - Disc'],loc=2, bbox_to_anchor = (1.0, 1.0))
plt.ylabel('Difference Magnetic Field Response (T)')
plt.xlabel('x (m)')
print 'Max Quadrature Difference: %e T' %(np.max(dCircQ[freq_ind,:]-dDiscQ[freq_ind,:]))

plt.tight_layout()
plt.show()









    



Max In-Phase Difference: 9.494297e-11 T
Max Quadrature Difference: 1.749723e-11 T



In [73]:

    
# Percent difference between circle and disc models
IPerc = (dCircI[freq_ind,:]-dDiscI[freq_ind,:])/np.max(np.abs([dCircI[freq_ind,:],dDiscI[freq_ind,:]]))*100
QPerc = (dCircQ[freq_ind,:]-dDiscQ[freq_ind,:])/np.max(np.abs([dCircQ[freq_ind,:],dDiscQ[freq_ind,:]]))*100

plt.figure(figsize=(11,6))

# Plot In-Phase
plt.subplot(121)
plt.plot(rxz,IPerc,linewidth=2)
plt.plot(rxz,perc_cutoff*np.ones(rxz.shape),'--k', rxz,-perc_cutoff*np.ones(rxz.shape),'--k',linewidth=2)
plt.grid(which='both')
plt.title('In-Phase')
plt.ylabel('Difference Magnetic Field Response (%)')
plt.xlabel('x (m)')
print 'Max In-Phase Difference: %d %s' %(np.max(IPerc), '%')

# Plot Quadrature 
plt.subplot(122)
plt.plot(rxz,QPerc ,linewidth=2)
plt.plot(rxz,perc_cutoff*np.ones(rxz.shape),'--k', rxz,-perc_cutoff*np.ones(rxz.shape),'--k',linewidth=2)
plt.grid(which='both')
plt.title('Quadrature')
plt.legend(['(Circle - Disc)/max(Circle,Disc) *100%'],loc=2, bbox_to_anchor = (1.0, 1.0))
plt.ylabel('Difference Magnetic Field Response (%)')
plt.xlabel('x (m)')
print 'Max Quadrature Difference: %d %s' %(np.max(QPerc), '%')

plt.tight_layout()
plt.show()









    



Max In-Phase Difference: 1 %
Max Quadrature Difference: 43 %

Plot the results as a function of frequency at a specific receiver location



In [74]:

    
# choose receiver position to plot 
rxind = np.floor(nrx/2.) # plot the middle receiver
print rxz[rxind]

0.0



In [75]:

    
plt.figure(figsize=(11,6))

# Plot In-Phase
plt.subplot(121)
plt.semilogx(freq,dCircI[:,rxind],'-o',freq,dDiscI[:,rxind],'-o',linewidth=2)
plt.plot(freq,flr*np.ones(freq.shape),'--k', freq,-flr*np.ones(freq.shape),'--k',linewidth=2)
plt.grid(which='both')
plt.title('In-Phase')
plt.ylabel('Magnetic Field Response (T)')
plt.xlabel('frequency (Hz)')
print 'Max In-Phase Responses: Circle %e T , Disc %e T' %(np.max(dCircI), np.max(dDiscI))

# Plot Quadrature
plt.subplot(122)
plt.semilogx(freq,dCircQ[:,rxind],'-o',freq,dDiscQ[:,rxind],'-o',linewidth=2)
plt.plot(freq,flr*np.ones(freq.shape),'--k', freq,-flr*np.ones(freq.shape),'--k',linewidth=2)
plt.grid(which='both')
plt.title('Quadrature')
plt.legend(['Circle','Disc'],loc=2, bbox_to_anchor = (1.0, 1.0),title='Model')
plt.ylabel('Magnetic Field Response (T)')
plt.xlabel('fraquency (Hz)')
print 'Max Quadrature Responses: Circle %e T , Disc %e T' %(np.max(dCircQ), np.max(dDiscQ))


plt.tight_layout()
plt.show()









    



Max In-Phase Responses: Circle 4.034892e-11 T , Disc 4.035495e-11 T
Max Quadrature Responses: Circle 1.379466e-10 T , Disc 1.039302e-10 T



In [65]:

    
# Difference between circle and disc signals
Idiff = dCircI[:,rxind]-dDiscI[:,rxind]
Qdiff = dCircQ[:,rxind]-dDiscQ[:,rxind]

plt.figure(figsize=(11,6))

# Plot In-phase
plt.subplot(121)
plt.semilogx(freq,Idiff,'-o',linewidth=2)
plt.plot(freq,flr*np.ones(freq.shape),'--k', freq,-flr*np.ones(freq.shape),'--k',linewidth=2)
plt.grid(which='both')
plt.title('In-Phase')
plt.ylabel('Difference Magnetic Field Response (T)')
plt.xlabel('frequency (Hz)')
print 'Max In-Phase Difference: %e T' %(np.max(dCircI-dDiscI))

# Plot Quadrature
plt.subplot(122)
plt.semilogx(freq,Qdiff,'-o',linewidth=2)
plt.plot(freq,flr*np.ones(freq.shape),'--k', freq,-flr*np.ones(freq.shape),'--k',linewidth=2)
plt.grid(which='both')
plt.title('Quadrature')
plt.legend(['Circle - Disc'],loc=2, bbox_to_anchor = (1.0, 1.0))
plt.ylabel('Difference Magnetic Field Response (T)')
plt.xlabel('x (m)')
print 'Max Quadrature Difference: %e T' %(np.max(dCircQ-dDiscQ))

plt.tight_layout()
plt.show()









    



Max In-Phase Difference: 9.494297e-11 T
Max Quadrature Difference: 4.842180e-11 T



In [64]: