wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-26 13:37:21 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge remote-tracking branch 'origin/Examples' into ex/mt1d

Browse Source 
# Conflicts:
#	SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py
#	SimPEG/Examples/__init__.py
#	SimPEG/Examples/sphereElectrostatic_example.py
#	SimPEG/Optimization.py
This commit is contained in:
Thibaut Astic
2016-04-07 11:19:54 -07:00
parent 45c4fa0d95 f3c9626133
commit c86b9bdd6a
25 changed files with 1089 additions and 765 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/Examples/EM_Schenkel_Morrison_Casing.py
+275
View File
@@ -0,0 +1,275 @@
from SimPEG import *
from SimPEG.EM import FDEM, Analytics, mu_0
import time
try:
from pymatsolver import MumpsSolver
solver = MumpsSolver
except Exception:
solver = SolverLU
pass
def run(plotIt=True):
"""
EM: Schenkel and Morrison Casing Model
======================================
Here we create and run a FDEM forward simulation to calculate the vertical
current inside a steel-cased. The model is based on the Schenkel and
Morrison Casing Model, and the results are used in a 2016 SEG abstract by
Yang et al.
- Schenkel, C.J., and H.F. Morrison, 1990, Effects of well casing on potential field measurements using downhole current sources: Geophysical prospecting, 38, 663-686.
The model consists of:
- Air: Conductivity 1e-8 S/m, above z = 0
- Background: conductivity 1e-2 S/m, below z = 0
- Casing: conductivity 1e6 S/m
- 300m long
- radius of 0.1m
- thickness of 6e-3m
Inside the casing, we take the same conductivity as the background.
We are using an EM code to simulate DC, so we use frequency low enough
that the skin depth inside the casing is longer than the casing length (f
= 1e-6 Hz). The plot produced is of the current inside the casing.
These results are shown in the SEG abstract by Yang et al., 2016: 3D DC
resistivity modeling of steel casing for reservoir monitoring using
equivalent resistor network. The solver used to produce these results and
achieve the CPU time of ~30s is Mumps, which was installed using pymatsolver_
.. _pymatsolver: https://github.com/rowanc1/pymatsolver
This example is on figshare: https://dx.doi.org/10.6084/m9.figshare.3126961.v1
If you would use this example for a code comparison, or build upon it, a
citation would be much appreciated!
"""
if plotIt:
import matplotlib.pylab as plt
# ------------------ MODEL ------------------
sigmaair = 1e-8 # air
sigmaback = 1e-2 # background
sigmacasing = 1e6 # casing
sigmainside = sigmaback # inside the casing
casing_t = 0.006 # 1cm thickness
casing_l = 300 # length of the casing
casing_r = 0.1
casing_a = casing_r - casing_t/2. # inner radius
casing_b = casing_r + casing_t/2. # outer radius
casing_z = np.r_[-casing_l,0.]
# ------------------ SURVEY PARAMETERS ------------------
freqs = np.r_[1e-6] #[1e-1, 1, 5] # frequencies
dsz = -300 # down-hole z source location
src_loc = np.r_[0.,0.,dsz]
inf_loc = np.r_[0.,0.,1e4]
print 'Skin Depth: ', [(500./np.sqrt(sigmaback*_)) for _ in freqs]
# ------------------ MESH ------------------
# fine cells near well bore
csx1, csx2 = 2e-3, 60.
pfx1, pfx2 = 1.3, 1.3
ncx1 = np.ceil(casing_b/csx1+2)
# pad nicely to second cell size
npadx1 = np.floor(np.log(csx2/csx1) / np.log(pfx1))
hx1a,hx1b = Utils.meshTensor([(csx1,ncx1)]),Utils.meshTensor([(csx1,npadx1,pfx1)])
dx1 = sum(hx1a)+sum(hx1b)
dx1 = np.floor(dx1/csx2)
hx1b *= (dx1*csx2 - sum(hx1a))/sum(hx1b)
# second chunk of mesh
dx2 = 300. # uniform mesh out to here
ncx2 = np.ceil((dx2 - dx1)/csx2)
npadx2 = 45
hx2a, hx2b = Utils.meshTensor([(csx2,ncx2)]), Utils.meshTensor([(csx2,npadx2,pfx2)])
hx = np.hstack([hx1a,hx1b,hx2a,hx2b])
# z-direction
csz = 0.05
nza = 10
ncz, npadzu, npadzd = np.int(np.ceil(np.diff(casing_z)[0]/csz))+10, 68, 68 # cell size, number of core cells, number of padding cells in the x- direction
hz = Utils.meshTensor([(csz,npadzd,-1.3), (csz,ncz), (csz,npadzu,1.3)]) # vector of cell widths in the z-direction
# Mesh
mesh = Mesh.CylMesh([hx,1.,hz], [0.,0.,-np.sum(hz[:npadzu+ncz-nza])])
print 'Mesh Extent xmax: %f,: zmin: %f, zmax: %f'%(mesh.vectorCCx.max(), mesh.vectorCCz.min(), mesh.vectorCCz.max())
print 'Number of cells', mesh.nC
if plotIt is True:
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
ax.set_title('Simulation Mesh')
mesh.plotGrid(ax=ax)
plt.show()
# Put the model on the mesh
sigWholespace = sigmaback*np.ones((mesh.nC))
sigBack = sigWholespace.copy()
sigBack[mesh.gridCC[:,2] > 0.] = sigmaair
sigCasing = sigBack.copy()
iCasingZ = (mesh.gridCC[:,2] <= casing_z[1]) & (mesh.gridCC[:,2] >= casing_z[0])
iCasingX = (mesh.gridCC[:,0] >= casing_a) & (mesh.gridCC[:,0] <= casing_b)
iCasing = iCasingX & iCasingZ
sigCasing[iCasing] = sigmacasing
if plotIt is True:
# plotting parameters
xlim = np.r_[0., 0.2]
zlim = np.r_[-350., 10.]
clim_sig = np.r_[-8,6]
# plot models
fig, ax = plt.subplots(1,1,figsize=(4,4))
f = plt.colorbar(mesh.plotImage(np.log10(sigCasing),ax=ax)[0], ax=ax)
ax.grid(which='both')
ax.set_title('Log_10 (Sigma)')
ax.set_xlim(xlim)
ax.set_ylim(zlim)
f.set_clim(clim_sig)
plt.show()
# -------------- Sources --------------------
# Define Custom Current Sources
# surface source
sg_x = np.zeros(mesh.vnF[0],dtype=complex)
sg_y = np.zeros(mesh.vnF[1],dtype=complex)
sg_z = np.zeros(mesh.vnF[2],dtype=complex)
nza = 2 # put the wire two cells above the surface
ncin = 2
# vertically directed wire
sgv_indx = (mesh.gridFz[:,0] > casing_a) & (mesh.gridFz[:,0] < casing_a + csx1) # hook it up to casing at the surface
sgv_indz = (mesh.gridFz[:,2] <= +csz*nza) & (mesh.gridFz[:,2] >= -csz*2)
sgv_ind = sgv_indx & sgv_indz
sg_z[sgv_ind] = -1.
# horizontally directed wire
sgh_indx = (mesh.gridFx[:,0] > casing_a) & (mesh.gridFx[:,0] <= inf_loc[2])
sgh_indz = (mesh.gridFx[:,2] > csz*(nza-0.5)) & (mesh.gridFx[:,2] < csz*(nza+0.5))
sgh_ind = sgh_indx & sgh_indz
sg_x[sgh_ind] = -1.
sgv2_indx = (mesh.gridFz[:,0] >= mesh.gridFx[sgh_ind,0].max()) & (mesh.gridFz[:,0] <= inf_loc[2]*1.2) # hook it up to casing at the surface
sgv2_indz = (mesh.gridFz[:,2] <= +csz*nza) & (mesh.gridFz[:,2] >= -csz*2)
sgv2_ind = sgv2_indx & sgv2_indz
sg_z[sgv2_ind] = 1.
# assemble the source
sg = np.hstack([sg_x,sg_y,sg_z])
sg_p = [FDEM.Src.RawVec_e([],_,sg/mesh.area) for _ in freqs]
# downhole source
dg_x = np.zeros(mesh.vnF[0],dtype=complex)
dg_y = np.zeros(mesh.vnF[1],dtype=complex)
dg_z = np.zeros(mesh.vnF[2],dtype=complex)
# vertically directed wire
dgv_indx = (mesh.gridFz[:,0] < csx1) # go through the center of the well
dgv_indz = (mesh.gridFz[:,2] <= +csz*nza) & (mesh.gridFz[:,2] > dsz + csz/2.)
dgv_ind = dgv_indx & dgv_indz
dg_z[dgv_ind] = -1.
# couple to the casing downhole
dgh_indx = mesh.gridFx[:,0] < casing_a + csx1
dgh_indz = (mesh.gridFx[:,2] < dsz + csz) & (mesh.gridFx[:,2] >= dsz)
dgh_ind = dgh_indx & dgh_indz
dg_x[dgh_ind] = 1.
# horizontal part at surface
dgh2_indx = mesh.gridFx[:,0] <= inf_loc[2]*1.2
dgh2_indz = sgh_indz.copy()
dgh2_ind = dgh2_indx & dgh2_indz
dg_x[dgh2_ind] = -1.
# vertical part at surface
dgv2_ind = sgv2_ind.copy()
dg_z[dgv2_ind] = 1.
# assemble the source
dg = np.hstack([dg_x,dg_y,dg_z])
dg_p = [FDEM.Src.RawVec_e([],_,dg/mesh.area) for _ in freqs]
# ------------ Problem and Survey ---------------
survey = FDEM.Survey(sg_p + dg_p)
mapping = [('sigma', Maps.IdentityMap(mesh))]
problem = FDEM.Problem_h(mesh, mapping=mapping)
problem.pair(survey)
# ------------- Solve ---------------------------
t0 = time.time()
fieldsCasing = problem.fields(sigCasing)
print 'Time to solve 2 sources', time.time() - t0
# Plot current
# current density
jn0 = fieldsCasing[dg_p,'j']
jn1 = fieldsCasing[sg_p,'j']
# current
in0 = [mesh.area*fieldsCasing[dg_p,'j'][:,i] for i in range(len(freqs))]
in1 = [mesh.area*fieldsCasing[sg_p,'j'][:,i] for i in range(len(freqs))]
in0 = np.vstack(in0).T
in1 = np.vstack(in1).T
# integrate to get z-current inside casing
inds_inx = (mesh.gridFz[:,0] >= casing_a) & (mesh.gridFz[:,0] <= casing_b)
inds_inz = (mesh.gridFz[:,2] >= dsz ) & (mesh.gridFz[:,2] <= 0)
inds_fz = inds_inx & inds_inz
indsx = [False]*mesh.nFx
inds = list(indsx) + list(inds_fz)
in0_in = in0[np.r_[inds]]
in1_in = in1[np.r_[inds]]
z_in = mesh.gridFz[inds_fz,2]
in0_in = in0_in.reshape([in0_in.shape[0]/3,3])
in1_in = in1_in.reshape([in1_in.shape[0]/3,3])
z_in = z_in.reshape([z_in.shape[0]/3,3])
I0 = in0_in.sum(1).real
I1 = in1_in.sum(1).real
z_in = z_in[:,0]
if plotIt is True:
fig, ax = plt.subplots(1,2,figsize=(12,4))
ax[0].plot(z_in,np.absolute(I0), z_in,np.absolute(I1))
ax[0].legend(['top casing', 'bottom casing'],loc='best')
ax[0].set_title('Magnitude of Vertical Current in Casing')
ax[1].semilogy(z_in,np.absolute(I0), z_in,np.absolute(I1))
ax[1].legend(['top casing', 'bottom casing'],loc='best')
ax[1].set_title('Magnitude of Vertical Current in Casing')
ax[1].set_ylim([1e-2, 1.])
plt.show()
if __name__ == '__main__':
run()
SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py
+87 -101
View File
@@ -1,7 +1,8 @@
from scipy.constants import epsilon_0, mu_0
import matplotlib.pyplot as plt
import numpy as np
#from SimPEG.EM.Utils import k, omega
from ipywidgets import *
from SimPEG.EM.Utils import k, omega
"""
MT1D: n layered earth problem
@@ -15,51 +16,45 @@ This code compute the analytic response of a n-layered Earth to a plane wave (Ma
We start by looking at Maxwell's equations in the electric
field \\\(\\\mathbf{E}\\) and the magnetic flux
\\\(\\\mathbf{H}\\) to write the wave equations
\\\(\\\mathbf{H}\\) to write the wave equations
\\(\\ \nabla ^2 \mathbf{E_x} + k^2 \mathbf{E_x} = 0 \\) &
\\(\\ \nabla ^2 \mathbf{H_y} + k^2 \mathbf{H_y} = 0 \\)
Then solving the equations in each layer "j" between z_{j-1} and z_j in the form of
Then solving the equations in each layer "j" between z_{j-1} and z_j in the form of
\\(\\ E_{x,j} (z) = U_j e^{i k (z-z_{j-1})} + D_j e^{-i k (z-z_{j-1})} \\)
\\(\\ H_{y,j} (z) = \frac{1}{Z_j} (D_j e^{-i k (z-z_{j-1})} - U_j e^{i k (z-z_{j-1})}) \\)
With U and D the Up and Down components of the E-field.
The iteration from one layer to another is ensure by:
The iteration from one layer to another is ensure by:
\\(\\ \left(\begin{matrix} E_{x,j} \\ H_{y,j} \end{matrix} \right) =
\\(\\ \left(\begin{matrix} E_{x,j} \\ H_{y,j} \end{matrix} \right) =
P_j T_j P^{-1}_J \left(\begin{matrix} E_{x,j+1} \\ H_{y,j+1} \end{matrix} \right) \\)
And the Boundary Condition is set for the E-field in the last layer, with no Up component (=0)
And the Boundary Condition is set for the E-field in the last layer, with no Up component (=0)
and only a down component (=1 then normalized by the highest amplitude to ensure numeric stability)
The layer 0 is assumed to be the air layer.
"""
#Frequency conversion
omega = lambda f: 2.*np.pi*f
#Evaluate k wavenumber
k = lambda mu,sig,eps,f: np.sqrt(mu*mu_0*eps*epsilon_0*(2.*np.pi*f)**2.-1.j*mu*mu_0*sig*omega(f))
#Define a frquency range for a survey
frange = lambda minfreq, maxfreq, step: np.logspace(minfreq,maxfreq,num = step, base = 10.)
#Functions to create random physical Perties for a n-layered earth
thick = lambda minthick, maxthick, nlayer: np.append(np.array([1.2*10.**5]),
thick = lambda minthick, maxthick, nlayer: np.append(np.array([1.2*10.**5]),
np.ndarray.round(minthick + (maxthick-minthick)* np.random.rand(nlayer-1,1)
,decimals =1))
sig = lambda minsig, maxsig, nlayer: np.append(np.array([0.]),
sig = lambda minsig, maxsig, nlayer: np.append(np.array([0.]),
np.ndarray.round(10.**minsig + (10.**maxsig-10.**minsig)* np.random.rand(nlayer,1)
,decimals=3))
mu = lambda minmu, maxmu, nlayer: np.append(np.array([1.]),
mu = lambda minmu, maxmu, nlayer: np.append(np.array([1.]),
np.ndarray.round(minmu + (maxmu-minmu)* np.random.rand(nlayer,1)
,decimals=1))
eps = lambda mineps, maxeps, nlayer: np.append(np.array([1.]),
eps = lambda mineps, maxeps, nlayer: np.append(np.array([1.]),
np.ndarray.round(mineps + (maxeps-mineps)* np.random.rand(nlayer,1)
,decimals=1))
@@ -69,17 +64,8 @@ ImpZ = lambda f, mu, k: omega(f)*mu*mu_0/k
#Complex Cole-Cole Conductivity - EM utils
PCC= lambda siginf,m,t,c,f: siginf*(1.-(m/(1.+(1j*omega(f)*t)**c)))
#Converted thickness array into top of layer array
def top(thick):
topv= np.zeros(len(thick)+1)
topv[0]=-thick[0]
for i in range(1,len(topv),1):
topv[i] = topv[i-1] + thick[i-1]
return topv
top = lambda thick: np.cumsum(thick)
#Propagation Matrix and theirs inverses
@@ -104,36 +90,36 @@ H_ZT = lambda U,D,Z,f,t : (1./Z)*np.exp(1j*omega(f)*t)*(D-U)
#Plot the configuration of the problem
def PlotConfiguration(thick,sig,eps,mu,ax,widthg,z):
topn = top(thick)
widthn = np.arange(-widthg,widthg+widthg/10.,widthg/10.)
ax.set_ylim([z.min(),z.max()])
ax.set_xlim([-widthg,widthg])
ax.set_ylabel("Depth (m)", fontsize=16.)
ax.yaxis.tick_right()
ax.yaxis.set_label_position("right")
#define filling for the different layers
hatches=['/' , '+', 'x', '|' , '\\', '-' , 'o' , 'O' , '.' , '*' ]
hatches=['/' , '+', 'x', '|' , '\\', '-' , 'o' , 'O' , '.' , '*' ]
#Write the physical properties of air
ax.annotate(("Air, $\sigma$ =%1.0f mS/m")%(sig[0]*10**(3)),
xy=(-widthg/2., -np.abs(z.max())/2.), xycoords='data',
xytext=(-widthg/2., -np.abs(z.max())/2.), textcoords='data',
fontsize=14.)
ax.annotate(("$\epsilon_r$= %1i")%(eps[0]),
xy=(-widthg/2., -np.abs(z.max())/3.), xycoords='data',
xytext=(-widthg/2., -np.abs(z.max())/3.), textcoords='data',
fontsize=14.)
ax.annotate(("$\mu_r$= %1i")%(mu[0]),
xy=(-widthg/2., -np.abs(z.max())/3.), xycoords='data',
xytext=(0, -np.abs(z.max())/3.), textcoords='data',
fontsize=14.)
#Write the physical properties of the differents layers up to the (n-1)-th and fill it with pattern
for i in range(1,len(topn)-1,1):
if topn[i] == topn[i+1]:
@@ -143,39 +129,39 @@ def PlotConfiguration(thick,sig,eps,mu,ax,widthg,z):
xy=(0., (2.*topn[i]+topn[i+1])/3), xycoords='data',
xytext=(0., (2.*topn[i]+topn[i+1])/3), textcoords='data',
fontsize=14.)
ax.annotate(("$\epsilon_r$= %1i")%(eps[i]),
xy=(-widthg/1.1, (2.*topn[i]+topn[i+1])/3), xycoords='data',
xytext=(-widthg/1.1, (2.*topn[i]+topn[i+1])/3), textcoords='data',
fontsize=14.)
ax.annotate(("$\mu_r$= %1.2f")%(mu[i]),
xy=(-widthg/2., (2.*topn[i]+topn[i+1])/3), xycoords='data',
xytext=(-widthg/2., (2.*topn[i]+topn[i+1])/3), textcoords='data',
fontsize=14.)
ax.plot(widthn,topn[i]*np.ones_like(widthn),color='black')
ax.plot(widthn,topn[i]*np.ones_like(widthn),color='black')
ax.fill_between(widthn,topn[i],topn[i+1],alpha=0.3,color="none",edgecolor='black', hatch=hatches[(i-1)%10])
#Write the physical properties of the n-th layer and fill it with pattern
ax.plot(widthn,topn[-1]*np.ones_like(widthn),color='black')
ax.plot(widthn,topn[-1]*np.ones_like(widthn),color='black')
ax.fill_between(widthn,topn[-1],z.max(),alpha=0.3,color="none",edgecolor='black', hatch=hatches[(len(topn)-2)%10])
ax.annotate(("$\sigma$ =%3.3f mS/m")%(sig[-1]*10**(3)),
xy=(0., (2.*topn[-1]+z.max())/3), xycoords='data',
xytext=(0., (2.*topn[-1]+z.max())/3), textcoords='data',
fontsize=14.)
ax.annotate(("$\epsilon_r$= %1i")%(eps[-1]),
xy=(-widthg/1.1, (2.*topn[-1]+z.max())/3), xycoords='data',
xytext=(-widthg/1.1, (2.*topn[-1]+z.max())/3), textcoords='data',
fontsize=14.)
ax.annotate(("$\mu_r$= %1.2f")%(mu[-1]),
xy=(-widthg/2., (2.*topn[-1]+z.max())/3), xycoords='data',
xytext=(-widthg/2., (2.*topn[-1]+z.max())/3), textcoords='data',
fontsize=14.)
#plot Trees!
ax.annotate("",
xy=(widthg/2., -1.*z.max()/5.), xycoords='data',
@@ -194,7 +180,7 @@ def PlotConfiguration(thick,sig,eps,mu,ax,widthg,z):
xytext=(widthg/2., 0.), textcoords='data',
arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.)
)
ax.annotate("",
xy=(1.2*widthg/2., -1.*z.max()/5.), xycoords='data',
xytext=(1.2*widthg/2., 0.), textcoords='data',
@@ -231,7 +217,7 @@ def PlotConfiguration(thick,sig,eps,mu,ax,widthg,z):
arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.)
)
ax.invert_yaxis()
return ax
@@ -239,83 +225,83 @@ def PlotConfiguration(thick,sig,eps,mu,ax,widthg,z):
#Propagate Up and Down component for a certain frequency & evaluate E and H field
def Propagate(f,H,sig,chg,taux,c,mu,eps,n):
sigcm = np.zeros_like(sig,dtype='complex_')
for j in range(1,len(sig)):
sigcm[j]=PCC(sig[j],chg[j],taux[j],c[j],f)
K = k(mu,sigcm,eps,f)
K = k(f, sigcm, mu, eps)
Z = ImpZ(f,mu,K)
EH = np.matrix(np.zeros((2,n+1),dtype = 'complex_'),dtype = 'complex_')
UD = np.matrix(np.zeros((2,n+1),dtype = 'complex_'),dtype = 'complex_')
UD[1,-1] = 1.
for i in range(-2,-(n+2),-1):
UD[:,i] = Tinv(H[i+1],K[i])*Pinv(Z[i])*P(Z[i+1])*UD[:,i+1]
UD = UD/((np.abs(UD[0,:]+UD[1,:])).max())
for j in range(0,n+1):
for j in range(0,n+1):
EH[:,j] = np.matrix([[1.,1,],[-1./Z[j],1./Z[j]]])*UD[:,j]
return UD, EH, Z ,K
#Evaluate the apparent resistivity and phase for a frequency range
def appres(F,H,sig,chg,taux,c,mu,eps,n):
Res = np.zeros_like(F)
Phase = np.zeros_like(F)
App_ImpZ= np.zeros_like(F,dtype='complex_')
for i in range(0,len(F)):
UD,EH,Z ,K = Propagate(F[i],H,sig,chg,taux,c,mu,eps,n)
App_ImpZ[i] = EH[0,1]/EH[1,1]
Res[i] = np.abs(App_ImpZ[i])**2./(mu_0*omega(F[i]))
Phase[i] = np.angle(App_ImpZ[i], deg = True)
return Res,Phase
#Evaluate Up, Down components, E and H field, for a frequency range,
#a discretized depth range and a time range (use to calculate envelope)
#Evaluate Up, Down components, E and H field, for a frequency range,
#a discretized depth range and a time range (use to calculate envelope)
def calculateEHzt(F,H,sig,chg,taux,c,mu,eps,n,zsample,tsample):
topc = top(H)
layer = np.zeros(len(zsample),dtype=np.int)-1
Exzt = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
Hyzt = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
Uz = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
Dz = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
UDaux = np.matrix(np.zeros((2,len(zsample)),dtype = 'complex_'),dtype = 'complex_')
for i in range(0,n+1,1):
layer = layer+(zsample>=topc[i])*1
for j in range(0,len(F)):
UD,EH,Z ,K = Propagate(F[j],H,sig,chg,taux,c,mu,eps,n)
for p in range(0,len(zsample)):
UDaux[:,p] = UD_Z(UD[:,layer[p]],zsample[p],topc[layer[p]],K[layer[p]])
for q in range(0,len(tsample)):
Exzt[p,q] = Exzt[p,q] + E_ZT(UDaux[0,p],UDaux[1,p],F[j],tsample[q])/len(F)
Hyzt[p,q] = Hyzt[p,q] + H_ZT(UDaux[0,p],UDaux[1,p],Z[layer[p]],F[j],tsample[q])/len(F)
Uz[p,q] = Uz[p,q] + UDaux[0,p]*np.exp(1j*omega(F[j])*tsample[q])/len(F)
Dz[p,q] = Dz[p,q] + UDaux[1,p]*np.exp(1j*omega(F[j])*tsample[q])/len(F)
return Exzt,Hyzt,Uz,Dz,UDaux,layer
#Function to Plot Apparent Resistivity and Phase
def PlotAppRes(F,H,sig,chg,taux,c,mu,eps,n,fenvelope,PlotEnvelope):
@@ -334,44 +320,44 @@ def PlotAppRes(F,H,sig,chg,taux,c,mu,eps,n,fenvelope,PlotEnvelope):
ax[0].grid(which='major')
ax0 = ax[0].twiny()
ax0.set_xlim([0.,90.])
ax0.set_ylim([F.min(),F.max()])
ax0.scatter(Phase,F,color='purple')
ax0.set_xlabel('Phase (Degrees)',fontsize=16.,color="purple")
zc=np.arange(-(H[1:].max()+10)*n,(H[1:].max()+10)*n,10.)
ax[0].tick_params(labelsize=16)
ax[1].tick_params(labelsize=16)
ax0.tick_params(labelsize=16)
if PlotEnvelope:
widthn=np.logspace(np.log10(Res.min())-1., np.log10(Res.max())+1., num=100, endpoint=True, base=10.0)
fenvelope1n=np.ones(100)*fenvelope
ax[0].plot(widthn,fenvelope1n,linestyle='dashed',color='black')
tc=np.arange(0.,1./fenvelope,0.01/(fenvelope))
Exzt,Hyzt,Uz,Dz,UDaux,layer = calculateEHzt(np.array([fenvelope]),H,sig,chg,taux,c,mu,eps,n,zc,tc)
ax1=ax[1].twiny()
ax[1].tick_params(labelsize=16)
ax1.tick_params(labelsize=16)
ax[1].set_xlabel('Amplitude Electric Field E (V/m)',color='blue',fontsize=16)
ax1.set_xlabel('Amplitude Magnetic Field H (A/m)',color='red',fontsize=16)
ax[1].fill_betweenx(zc,np.squeeze(np.asarray(np.real(Exzt.min(axis=1)))),
np.squeeze(np.asarray(np.real(Exzt.max(axis=1)))),
np.squeeze(np.asarray(np.real(Exzt.max(axis=1)))),
color='blue', alpha=0.1)
ax1.fill_betweenx(zc,np.squeeze(np.asarray(np.real(Hyzt.min(axis=1)))),
np.squeeze(np.asarray(np.real(Hyzt.max(axis=1)))),
np.squeeze(np.asarray(np.real(Hyzt.max(axis=1)))),
color='red', alpha=0.1)
ax[1] = PlotConfiguration(H,sig,eps,mu,ax[1],(1.5*np.abs(Exzt).max()),zc)
ax1.set_xlim([-1.5*np.abs(Hyzt).max(),1.5*np.abs(Hyzt).max()])
ax1.set_xlim([-1.5*np.abs(Hyzt).max(),1.5*np.abs(Hyzt).max()])
@@ -379,12 +365,12 @@ def PlotAppRes(F,H,sig,chg,taux,c,mu,eps,n,fenvelope,PlotEnvelope):
print 'No envelop (if True, might be slow)'
ax[1] = PlotConfiguration(H,sig,eps,mu,ax[1],1.,zc)
ax[1].get_xaxis().set_ticks([])
plt.show()
#Interactive MT for Notebook
def PlotAppRes3LayersInteract(h1,h2,sigl1,sigl2,sigl3,mul1,mul2,mul3,epsl1,epsl2,epsl3,PlotEnvelope,F_Envelope):
frangn=frange(-5,5,100.)
sig3= np.array([0.,0.001,0.1, 0.001])
thick3 = np.array([120000.,50.,50.])
@@ -394,7 +380,7 @@ def PlotAppRes3LayersInteract(h1,h2,sigl1,sigl2,sigl3,mul1,mul2,mul3,epsl1,epsl2
chg3_0=np.array([0.,0.1,0.,0.])
taux3=np.array([0.,0.1,0.,0.1])
c3=np.array([1.,1.,1.,1.])
sig3[1]=sigl1
sig3[1]=10.**sig3[1]
sig3[2]=sigl2
@@ -409,11 +395,11 @@ def PlotAppRes3LayersInteract(h1,h2,sigl1,sigl2,sigl3,mul1,mul2,mul3,epsl1,epsl2
eps3[3]=epsl3
thick3[1]=h1
thick3[2]=h2
PlotAppRes(frangn,thick3,sig3,chg3_0,taux3,c3,mu3,eps3,3,F_Envelope,PlotEnvelope)
def run(n=3,plotIt=True):
PlotAppRes(frangn,thick3,sig3,chg3_0,taux3,c3,mu3,eps3,3,F_Envelope,PlotEnvelope)
def run(n,plotIt=True):
# something to make a plot
F = frange(-5.,5.,20)
@@ -429,14 +415,14 @@ def run(n=3,plotIt=True):
if plotIt:
PlotAppRes(F, H, sign, chg, taux, c, mun, epsn, n, fenvelope=1000., PlotEnvelope=True)
PlotAppRes(F, H, sign, chg, taux, c, mun, epsn, n, fenvelope=1000., PlotEnvelope=True)
return Res, Phase
if __name__ == '__main__':
run()
run(3)
SimPEG/Examples/__init__.py
+5 -7
View File
@@ -3,8 +3,11 @@
##### AUTOIMPORTS #####
import DC_Analytic_Dipole
import DC_Forward_PseudoSection
import DC_PseudoSection_Simulation
import EM_FDEM_1D_Inversion
import EM_FDEM_Analytic_MagDipoleWholespace
import EM_FDEM_SusEffects
import EM_Schenkel_Morrison_Casing
import EM_TDEM_1D_Inversion
import FLOW_Richards_1D_Celia1990
import Forward_BasicDirectCurrent
@@ -16,17 +19,12 @@ import Mesh_QuadTree_Creation
import Mesh_QuadTree_FaceDiv
import Mesh_QuadTree_HangingNodes
import Mesh_Tensor_Creation
<<<<<<< HEAD
import MT_1D_analytic_nlayer_Earth
import sphereElectrostatic_example
__examples__ = ["EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_analytic_nlayer_Earth", "sphereElectrostatic_example"]
=======
import MT_1D_ForwardAndInversion
import MT_3D_Foward
import sphereElectrostatic_example
__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward"]
>>>>>>> master
__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "DC_PseudoSection_Simulation", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_FDEM_SusEffects", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_analytic_nlayer_Earth", "MT_1D_ForwardAndInversion", "MT_3D_Foward", "sphereElectrostatic_example"]
##### AUTOIMPORTS #####
SimPEG/Examples/sphereElectrostatic_example.py
+181 -197
View File
@@ -16,8 +16,6 @@ finally the charges accumulation.
Several plotting functions are defined for data visualisation.
Please visit http://em.geosci.xyz/en/latest/content/maxwell2_steady_state/electrostatic_sphere.html
for more examples using this code.
'''
@@ -34,7 +32,7 @@ sigf = lambda sig0,sig1: (sig1-sig0)/(sig1+2.*sig0)
def conductivity_log_wrapper(log_sig0,log_sig1):
sig0 = 10.**log_sig0
sig1 = 10.**log_sig1
return sig0,sig1
# Examples
@@ -56,7 +54,7 @@ def get_Setup(XYZ,sig0,sig1,R,E0,ax,label,colorsphere):
dx = xr[1]-xr[0]
top = np.sqrt(R**2-xplt**2)
bot = -np.sqrt(R**2-xplt**2)
if R != 0:
ax.plot(xplt, top, xplt, bot, color=colorsphere,linewidth=1.5)
ax.fill_between(xplt,bot,top,color=colorsphere,alpha=0.5 )
@@ -87,7 +85,7 @@ def get_Setup(XYZ,sig0,sig1,R,E0,ax,label,colorsphere):
ax.set_xticklabels([])
ax.set_yticklabels([])
ax.text(-1.,-np.sqrt(R)/2.-10.,'$\sigma_1$',fontsize=14)
ax.text(-0.05,-R-10,'$\sigma_0$',fontsize=14)
ax.text(-0.05,-R-10,'$\sigma_0$',fontsize=14)
ax.annotate(('$\mathbf{E_0} = E_0 \mathbf{\hat{x}}$ V/m'),
xy=(xr.min()+np.abs(xr.max()-xr.min())/20.,0), xycoords='data',
xytext=(xr.min()+np.abs(xr.max()-xr.min())/20.,0), textcoords='data',
@@ -98,7 +96,7 @@ def get_Setup(XYZ,sig0,sig1,R,E0,ax,label,colorsphere):
fontsize=14.)
ax.set_xlabel('x',fontsize=12)
ax.set_ylabel('y',fontsize=12)
else:
if label:
ax.annotate(("$\sigma_0$= %3.3f mS/m")%(sig0*10.**(3.)),
@@ -116,7 +114,7 @@ def get_Setup(XYZ,sig0,sig1,R,E0,ax,label,colorsphere):
else:
ax.set_xticklabels([])
ax.set_yticklabels([])
ax.text(-0.05,-10,'$\sigma_0$',fontsize=14)
ax.text(-0.05,-10,'$\sigma_0$',fontsize=14)
ax.text(xr.min()+np.abs(xr.max()-xr.min())/20., 0, '$\mathbf{E_0} = E_0 \mathbf{\hat{x}}$ V/m', fontsize=14)
ax.set_xlabel('x',fontsize=12)
ax.set_ylabel('y',fontsize=12)
@@ -130,8 +128,8 @@ def get_Setup(XYZ,sig0,sig1,R,E0,ax,label,colorsphere):
ax.set_aspect('equal')
return ax
def get_Conductivity(XYZ,sig0,sig1,R):
@@ -140,61 +138,61 @@ def get_Conductivity(XYZ,sig0,sig1,R):
'''
x,y,z = XYZ[:,0],XYZ[:,1],XYZ[:,2]
r_view=r(x,y,z)
ind0= (r_view>R)
ind1= (r_view<=R)
assert (ind0 + ind1).all(), 'Some indicies not included'
Sigma = np.zeros_like(x)
Sigma[ind0] = sig0
Sigma[ind1] = sig1
return Sigma
def get_Potential(XYZ,sig0,sig1,R,E0):
def get_Potential(XYZ,sig0,sig1,R,E0):
'''
Function that returns the total, the primary and the secondary potentials, assumes an x-oriented inducing field and that the sphere is at the origin
:input: grid, outer sigma, inner sigma, radius of the sphere, strength of the electric field
'''
x,y,z = XYZ[:,0],XYZ[:,1],XYZ[:,2]
sig_cur = sigf(sig0,sig1)
r_cur = r(x,y,z) # current radius
ind0 = (r_cur > R)
ind1 = (r_cur <= R)
assert (ind0 + ind1).all(), 'Some indicies not included'
Vt = np.zeros_like(x)
Vp = np.zeros_like(x)
Vs = np.zeros_like(x)
Vt[ind0] = -E0*x[ind0]*(1.-sig_cur*R**3./r_cur[ind0]**3.) # total potential outside the sphere
Vt[ind1] = -E0*x[ind1]*3.*sig0/(sig1+2.*sig0) # inside the sphere
Vp = - E0*x # primary potential
Vs = Vt - Vp # secondary potential
return Vt,Vp,Vs
#plot the primary potential on ax
def Plot_Primary_Potential(XYZ,sig0,sig1,R,E0,ax):
Vt,Vp,Vs = get_Potential(XYZ,sig0,sig1,R,E0)
xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
xcirc = xr[np.abs(xr) <= R]
Pplot = ax.pcolor(xr,yr,Vp.reshape(xr.size,yr.size))
ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k')
ax.set_title('Primary Potential',fontsize=ftsize_title)
@@ -207,19 +205,19 @@ def Plot_Primary_Potential(XYZ,sig0,sig1,R,E0,ax):
ax.set_xlabel('X coordinate ($m$)',fontsize = ftsize_label)
ax.set_aspect('equal')
ax.tick_params(labelsize=ftsize_axis)
return ax
#plot the total potential on ax
def Plot_Total_Potential(XYZ,sig0,sig1,R,E0,ax):
Vt,Vp,Vs = get_Potential(XYZ,sig0,sig1,R,E0)
xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
xcirc = xr[np.abs(xr) <= R]
Pplot = ax.pcolor(xr,yr,Vt.reshape(xr.size,yr.size))
ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k')
ax.set_title('Total Potential',fontsize=ftsize_title)
@@ -232,16 +230,16 @@ def Plot_Total_Potential(XYZ,sig0,sig1,R,E0,ax):
ax.set_xlabel('X coordinate ($m$)',fontsize = ftsize_label)
ax.set_aspect('equal')
ax.tick_params(labelsize=ftsize_axis)
return ax
#plot the secondary potential on ax
def Plot_Secondary_Potential(XYZ,sig0,sig1,R,E0,ax):
Vt,Vp,Vs = get_Potential(XYZ,sig0,sig1,R,E0)
xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
xcirc = xr[np.abs(xr) <= R]
Pplot = ax.pcolor(xr,yr,Vs.reshape(xr.size,yr.size))
@@ -256,30 +254,30 @@ def Plot_Secondary_Potential(XYZ,sig0,sig1,R,E0,ax):
ax.set_xlabel('X coordinate ($m$)',fontsize = ftsize_label)
ax.set_aspect('equal')
ax.tick_params(labelsize=ftsize_axis)
return ax
def get_ElectricField(XYZ,sig0,sig1,R,E0):
'''
Function that returns the total, the primary and the secondary electric fields,
Function that returns the total, the primary and the secondary electric fields,
input: grid, outer sigma, inner sigma, radius of the sphere, strength of the electric field
'''
x,y,z= XYZ[:,0], XYZ[:,1], XYZ[:,2]
r_cur=r(x,y,z) # current radius
ind0= (r_cur>R)
ind1= (r_cur<=R)
assert (ind0 + ind1).all(), 'Some indicies not included'
Ep = np.zeros(shape=(len(x),3))
Ep[:,0] = E0
Et = np.zeros(shape=(len(x),3))
Et[ind0,0] = E0 + E0*R**3./(r_cur[ind0]**5.)*sigf(sig0,sig1)*(2.*x[ind0]**2.-y[ind0]**2.-z[ind0]**2.);
Et[ind0,1] = E0*R**3./(r_cur[ind0]**5.)*3.*x[ind0]*y[ind0]*sigf(sig0,sig1);
Et[ind0,2] = E0*R**3./(r_cur[ind0]**5.)*3.*x[ind0]*z[ind0]*sigf(sig0,sig1);
@@ -287,16 +285,16 @@ def get_ElectricField(XYZ,sig0,sig1,R,E0):
Et[ind1,0] = 3.*sig0/(sig1+2.*sig0)*E0;
Et[ind1,1] = 0.;
Et[ind1,2] = 0.;
Es = Et - Ep
return Et, Ep, Es
#plot the total electric field on ax
def Plot_Total_ElectricField(XYZ,sig0,sig1,R,E0,ax):
Et, Ep, Es = get_ElectricField(XYZ,sig0,sig1,R,E0)
xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
xcirc = xr[np.abs(xr) <= R]
@@ -304,7 +302,7 @@ def Plot_Total_ElectricField(XYZ,sig0,sig1,R,E0,ax):
EtXr = Et[:,0].reshape(xr.size, yr.size)
EtYr = Et[:,1].reshape(xr.size, yr.size)
EtAmp = np.sqrt(Et[:,0]**2+Et[:,1]**2 + Et[:,2]**2).reshape(xr.size, yr.size)
ax.set_xlim([xr.min(),xr.max()])
ax.set_ylim([yr.min(),yr.max()])
ax.set_ylabel('Y coordinate ($m$)',fontsize = ftsize_label)
@@ -312,22 +310,22 @@ def Plot_Total_ElectricField(XYZ,sig0,sig1,R,E0,ax):
ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k')
ax.tick_params(labelsize=ftsize_axis)
ax.set_aspect('equal')
Eplot = ax.pcolor(xr,yr,EtAmp)
cb = plt.colorbar(Eplot,ax=ax)
cb.set_label(label= 'Amplitude ($V/m$)',size=ftsize_label) #weight='bold')
cb.ax.tick_params(labelsize=ftsize_axis)
ax.streamplot(xr,yr,EtXr,EtYr,color='gray',linewidth=2.,density=0.75)#angles='xy',scale_units='xy',scale=0.05)
ax.set_title('Total Field',fontsize=ftsize_title)
return ax
#plot the secondary electric field on ax
#plot the secondary electric field on ax
def Plot_Secondary_ElectricField(XYZ,sig0,sig1,R,E0,ax):
Et, Ep, Es = get_ElectricField(XYZ,sig0,sig1,R,E0)
xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
xcirc = xr[np.abs(xr) <= R]
@@ -335,7 +333,7 @@ def Plot_Secondary_ElectricField(XYZ,sig0,sig1,R,E0,ax):
EsXr = Es[:,0].reshape(xr.size, yr.size)
EsYr = Es[:,1].reshape(xr.size, yr.size)
EsAmp = np.sqrt(Es[:,0]**2+Es[:,1]**2+Es[:,2]**2).reshape(xr.size, yr.size)
ax.set_xlim([xr.min(),xr.max()])
ax.set_ylim([yr.min(),yr.max()])
ax.set_ylabel('Y coordinate ($m$)',fontsize = ftsize_label)
@@ -343,7 +341,7 @@ def Plot_Secondary_ElectricField(XYZ,sig0,sig1,R,E0,ax):
ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k')
ax.tick_params(labelsize=ftsize_axis)
ax.set_aspect('equal')
Eplot = ax.pcolor(xr,yr,EsAmp)
cb = plt.colorbar(Eplot,ax=ax)
cb.set_label(label= 'Amplitude ($V/m$)',size=ftsize_label) #weight='bold')
@@ -351,53 +349,53 @@ def Plot_Secondary_ElectricField(XYZ,sig0,sig1,R,E0,ax):
ax.streamplot(xr,yr,EsXr,EsYr,color='gray',linewidth=2.,density=0.75)#,angles='xy',scale_units='xy',scale=0.05)
ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k')
ax.set_title('Secondary Field',fontsize=ftsize_title)
return ax
def get_Current(XYZ,sig0,sig1,R,Et,Ep,Es):
'''
Function that returns the total, the primary and the secondary current densities,
Function that returns the total, the primary and the secondary current densities,
:input: grid, outer sigma, inner sigma, radius of the sphere, total, the primary and the seconadry electric fields,
'''
x,y,z= XYZ[:,0], XYZ[:,1], XYZ[:,2]
r_cur=r(x,y,z)
ind0= (r_cur>R)
ind1= (r_cur<=R)
assert (ind0 + ind1).all(), 'Some indicies not included'
Jt = np.zeros(shape=(len(x),3))
J0 = np.zeros(shape=(len(x),3))
Js = np.zeros(shape=(len(x),3))
Jp = sig0*Ep
Jt[ind0,:] = sig0*Et[ind0,:]
Jt[ind0,:] = sig0*Et[ind0,:]
Jt[ind1,:] = sig1*Et[ind1,:]
Js[ind0,:] = sig0*(Et[ind0,:]-Ep[ind0,:])
Js[ind1,:] = sig1*Et[ind1,:]-sig0*Ep[ind1,:]
return Jt,Jp,Js
#plot the total currents density on ax
def Plot_Total_Currents(XYZ,sig0,sig1,R,E0,ax):
Et,Ep,Es = get_ElectricField(XYZ,sig0,sig1,R,E0)
Jt,Jp,Js = get_Current(XYZ,sig0,sig1,R,Et,Ep,Es)
xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
xcirc = xr[np.abs(xr) <= R]
JtXr = Jt[:,0].reshape(xr.size, yr.size)
JtYr = Jt[:,1].reshape(xr.size, yr.size)
JtAmp = np.sqrt(Jt[:,0]**2+Jt[:,1]**2+Jt[:,2]**2).reshape(xr.size, yr.size)
ax.set_xlim([xr.min(),xr.max()])
ax.set_ylim([yr.min(),yr.max()])
ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k')
@@ -405,30 +403,30 @@ def Plot_Total_Currents(XYZ,sig0,sig1,R,E0,ax):
ax.set_xlabel('X coordinate ($m$)',fontsize=ftsize_label)
ax.tick_params(labelsize=ftsize_axis)
ax.set_aspect('equal')
Jplot = ax.pcolor(xr,yr,JtAmp.reshape(xr.size,yr.size))
cb = plt.colorbar(Jplot,ax=ax)
cb.set_label(label= 'Current Density ($A/m^2$)',size=ftsize_label) #weight='bold')
cb.ax.tick_params(labelsize=ftsize_axis)
ax.streamplot(xr,yr,JtXr,JtYr,color='gray',linewidth=2.,density=0.75)#,angles='xy',scale_units='xy',scale=1)
ax.set_title('Total Current Density',fontsize=ftsize_title)
return ax
#plot the secondary currents density on ax
def Plot_Secondary_Currents(XYZ,sig0,sig1,R,E0,ax):
Et,Ep,Es = get_ElectricField(XYZ,sig0,sig1,R,E0)
Jt,Jp,Js = get_Current(XYZ,sig0,sig1,R,Et,Ep,Es)
xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
xcirc = xr[np.abs(xr) <= R]
JsXr = Js[:,0].reshape(xr.size, yr.size)
JsYr = Js[:,1].reshape(xr.size, yr.size)
JsAmp = np.sqrt(Js[:,1]**2+Js[:,0]**2+Jt[:,2]**2).reshape(xr.size,yr.size)
ax.set_xlim([xr.min(),xr.max()])
ax.set_ylim([yr.min(),yr.max()])
ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k')
@@ -436,52 +434,52 @@ def Plot_Secondary_Currents(XYZ,sig0,sig1,R,E0,ax):
ax.set_xlabel('X coordinate ($m$)',fontsize=ftsize_label)
ax.tick_params(labelsize=ftsize_axis)
ax.set_aspect('equal')
Jplot = ax.pcolor(xr,yr,JsAmp.reshape(xr.size,yr.size))
cb = plt.colorbar(Jplot,ax=ax)
cb.set_label(label= 'Current Density ($A/m^2$)',size=ftsize_label) #weight='bold')
cb.ax.tick_params(labelsize=ftsize_axis)
ax.streamplot(xr,yr,JsXr,JsYr,color='gray',linewidth=2.,density=0.75)#,angles='xy',scale_units='xy',scale=1)
ax.set_title('Secondary Current Density',fontsize=ftsize_title)
return ax
def get_ChargesDensity(XYZ,sig0,sig1,R,Et,Ep):
'''
Function that returns the charges accumulation at the background/sphere interface,
Function that returns the charges accumulation at the background/sphere interface,
:input: grid, outer sigma, inner sigma, radius of the sphere, total and the primary electric fields,
'''
x,y,z= XYZ[:,0], XYZ[:,1], XYZ[:,2]
dx = x[1]-x[0]
r_cur=r(x,y,z)
ind0 = (r_cur > R)
ind1 = (r_cur < R)
ind2 = ((r_cur < (R+dx/2)) & (r_cur > (R-dx/2)) )
assert (ind0 + ind1 + ind2).all(), 'Some indicies not included'
rho = np.zeros_like(x)
rho[ind0] = 0
rho[ind1] = 0
rho[ind2] = epsilon_0*3.*Ep[ind2,0]*sigf(sig0,sig1)*x[ind2]/(np.sqrt(x[ind2]**2.+y[ind2]**2.))
return rho
#Plot charges density on ax
def Plot_ChargesDensity(XYZ,sig0,sig1,R,E0,ax):
xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
xcirc = xr[np.abs(xr) <= R]
Et, Ep, Es = get_ElectricField(XYZ,sig0,sig1,R,E0)
rho = get_ChargesDensity(XYZ,sig0,sig1,R,Et,Ep)
ax.set_xlim([xr.min(),xr.max()])
ax.set_ylim([yr.min(),yr.max()])
ax.set_aspect('equal')
@@ -494,11 +492,11 @@ def Plot_ChargesDensity(XYZ,sig0,sig1,R,E0,ax):
ax.set_xlabel('X coordinate ($m$)',fontsize=ftsize_label)
ax.tick_params(labelsize=ftsize_axis)
ax.set_title('Charges Density', fontsize=ftsize_title)
return ax
def MN_Potential_total(sig0,sig1,R,E0,start,end,nbmp,mn):
'''
Function that return array of midpoints electrodes, electrodes positions,
potentials differences for total and secondary potentials fields, unormalized and
@@ -515,20 +513,20 @@ def MN_Potential_total(sig0,sig1,R,E0,start,end,nbmp,mn):
#D: total distance from start to end
D = np.sqrt((start[0]-end[0])**2.+(start[1]-end[1])**2.)
#MP: dipoles'midpoint positions (x,y)
MP = np.zeros(shape=(nbmp,2))
MP = np.zeros(shape=(nbmp,2))
MP[:,0] = np.linspace(start[0],end[0],nbmp)
MP[:,1] = np.linspace(start[1],end[1],nbmp)
#Dipoles'Electrodes positions around each midpoints
EL = np.zeros(shape=(2*nbmp,2))
EL = np.zeros(shape=(2*nbmp,2))
for n in range(0,len(EL),2):
EL[n,0] = MP[n/2,0] - ((end[0]-start[0])/D)*mn/2.
EL[n+1,0] = MP[n/2,0] + ((end[0]-start[0])/D)*mn/2.
EL[n,1] = MP[n/2,1] - ((end[1]-start[1])/D)*mn/2.
EL[n+1,1] = MP[n/2,1] + ((end[1]-start[1])/D)*mn/2.
VtEL = np.zeros(2*nbmp) #Total Potential (Vt-) at each electrode (-EL)
VsEL = np.zeros(2*nbmp) #Secondary Potential (Vt-) at each electrode (-EL)
dVtMP = np.zeros(nbmp) #Diffence (d-) of Total Potential (Vt-) at each dipole (-MP)
@@ -537,109 +535,109 @@ def MN_Potential_total(sig0,sig1,R,E0,start,end,nbmp,mn):
dVsMPn = np.zeros(nbmp) #Diffence (d-) of Secondary Potential (Vt-) at each dipole (-MP) normalized for the mn spacing (n)
dVpMP = np.zeros(nbmp) #Diffence (d-) of Primary Potential (Vt-) at each dipole (-MP)
dVpMPn = np.zeros(nbmp) #Diffence (d-) of Primary Potential (Vt-) at each dipole (-MP) normalized for the mn spacing (n)
#Computing VtEL
#Computing VtEL
for m in range(0,2*nbmp):
if (r(EL[m,0],EL[m,1],0) > R):
VtEL[m] = -E0*EL[m,0]*(1.-sigf(sig0,sig1)*R**3./r(EL[m,0],EL[m,1],0)**3.)
else:
VtEL[m] = -E0*EL[m,0]*3.*sig0/(sig1+2.*sig0)
#Computing VsEL
VsEL = VtEL + E0*EL[:,0]
#Computing dVtMP, dVsMP
for p in range(0,nbmp):
dVtMP[p] = VtEL[2*p]-VtEL[2*p+1]
dVtMPn[p] = dVtMP[p]/mn
dVsMP[p] = VsEL[2*p]-VsEL[2*p+1]
dVsMPn[p] = dVsMP[p]/mn
return MP,EL,dVtMP,dVtMPn,dVsMP,dVsMPn
#Compare the DC response of two configurations
def two_configurations_comparison(XYZ,sig0,sig1,sig2,R0,R1,E0,xstart,ystart,xend,yend,nb_dipole,electrode_spacing,PlotOpt,ax):
def two_configurations_comparison(XYZ,sig0,sig1,sig2,R0,R1,E0,xstart,ystart,xend,yend,nb_dipole,electrode_spacing,PlotOpt):#,linearcolor):
#Define the mesh
xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
#Defining the Profile
start = np.array([xstart,ystart])
end = np.array([xend,yend])
#Calculating the data from the defined survey line for Configuration 0 and 1
MP0,EL0,VtdMP0,VtdMPn0,VsdMP0,VsdMPn0 = MN_Potential_total(sig0,sig1,R0,E0,start,end,nb_dipole,electrode_spacing)
MP1,EL1,VtdMP1,VtdMPn1,VsdMP1,VsdMPn1 = MN_Potential_total(sig0,sig2,R1,E0,start,end,nb_dipole,electrode_spacing)
# Initializing the figure
#fig = plt.figure(figsize=(20,20))
#ax0 = plt.subplot2grid((20,12), (0, 0),colspan=6,rowspan=6)
#ax1 = plt.subplot2grid((20,12), (0, 6),colspan=6,rowspan=6)
#ax2 = plt.subplot2grid((20,12), (16, 2), colspan=9,rowspan=4)
#ax3 = plt.subplot2grid((20,12), (8, 0),colspan=6,rowspan=6)
#ax4 = plt.subplot2grid((20,12), (8, 6),colspan=6,rowspan=6)
fig = plt.figure(figsize=(20,20))
ax0 = plt.subplot2grid((20,12), (0, 0),colspan=6,rowspan=6)
ax1 = plt.subplot2grid((20,12), (0, 6),colspan=6,rowspan=6)
ax2 = plt.subplot2grid((20,12), (16, 2), colspan=9,rowspan=4)
ax3 = plt.subplot2grid((20,12), (8, 0),colspan=6,rowspan=6)
ax4 = plt.subplot2grid((20,12), (8, 6),colspan=6,rowspan=6)
#Plotting the Configuration 0
ax[0] = get_Setup(XYZ,sig0,sig1,R0,E0,ax[0],True,[0.6,0.1,0.1])
ax0 = get_Setup(XYZ,sig0,sig1,R0,E0,ax0,True,[0.6,0.1,0.1])
#Plotting the Configuration 1
ax[1] = get_Setup(XYZ,sig0,sig2,R1,E0,ax[1],True,[0.1,0.1,0.6])
ax1 = get_Setup(XYZ,sig0,sig2,R1,E0,ax1,True,[0.1,0.1,0.6])
#Plotting the Data (Legends)
ax[2].set_title('Potential Differences',fontsize=ftsize_title)
ax[2].set_ylabel('Potential difference ($V$)',fontsize=ftsize_label)
ax[2].set_xlabel('Distance from start point ($m$)',fontsize=ftsize_label)
ax[2].tick_params(labelsize=ftsize_axis)
ax[2].grid()
ax2.set_title('Potential Differences',fontsize=ftsize_title)
ax2.set_ylabel('Potential difference ($V$)',fontsize=ftsize_label)
ax2.set_xlabel('Distance from start point ($m$)',fontsize=ftsize_label)
ax2.tick_params(labelsize=ftsize_axis)
ax2.grid()
if PlotOpt == 'Total':
ax[3]= Plot_Total_Potential(XYZ,sig0,sig1,R0,E0,ax[3])
ax[4]= Plot_Total_Potential(XYZ,sig0,sig2,R1,E0,ax[4])
#Plot the Data (from Configuration 0)
gphy0 = ax[2].plot(np.sqrt((MP0[0,0]-MP0[:,0])**2+(MP0[:,1]-MP0[0,1])**2),VtdMP0
ax3= Plot_Total_Potential(XYZ,sig0,sig1,R0,E0,ax3)
ax4= Plot_Total_Potential(XYZ,sig0,sig2,R1,E0,ax4)
#Plot the Data (from Configuration 0)
gphy0 = ax2.plot(np.sqrt((MP0[0,0]-MP0[:,0])**2+(MP0[:,1]-MP0[0,1])**2),VtdMP0
,marker='o',color='blue',linewidth=3.,label ='Left Model Response' )
#Plot the Data (from Configuration 1)
gphy1 = ax[2].plot(np.sqrt((MP1[0,0]-MP1[:,0])**2+(MP1[:,1]-MP1[0,1])**2),VtdMP1
gphy1 = ax2.plot(np.sqrt((MP1[0,0]-MP1[:,0])**2+(MP1[:,1]-MP1[0,1])**2),VtdMP1
,marker='o',color='red',linewidth=2.,label ='Right Model Response' )
ax[2].legend(('Left Model Response','Right Model Response'),loc=4)
ax2.legend(('Left Model Response','Right Model Response'),loc=4)
elif PlotOpt == 'Secondary':
#plot the secondary potentials
ax[3]= Plot_Secondary_Potential(XYZ,sig0,sig1,R0,E0,ax[3])
ax[4]= Plot_Secondary_Potential(XYZ,sig0,sig2,R1,E0,ax[4])
ax3= Plot_Secondary_Potential(XYZ,sig0,sig1,R0,E0,ax3)
ax4= Plot_Secondary_Potential(XYZ,sig0,sig2,R1,E0,ax4)
#Plot the data(from configuration 0)
gphy0 = ax[2].plot(np.sqrt((MP0[0,0]-MP0[:,0])**2+(MP0[:,1]-MP0[0,1])**2),VsdMP0,color='blue'
gphy0 = ax2.plot(np.sqrt((MP0[0,0]-MP0[:,0])**2+(MP0[:,1]-MP0[0,1])**2),VsdMP0,color='blue'
,marker='o',linewidth=3.,label ='Left Model Response' )
#Plot the Data (from Configuration 1)
gphy1 = ax[2].plot(np.sqrt((MP1[0,0]-MP1[:,0])**2+(MP1[:,1]-MP1[0,1])**2),VsdMP1
gphy1 = ax2.plot(np.sqrt((MP1[0,0]-MP1[:,0])**2+(MP1[:,1]-MP1[0,1])**2),VsdMP1
,marker='o',color='red',linewidth=2.,label ='Right Model Response' )
ax[2].legend(('Left Model Response','Right Model Response'),loc=4 )
ax2.legend(('Left Model Response','Right Model Response'),loc=4 )
else:
print('What dont you get? Total or Secondary?')
#Legends
ax[3].plot(MP0[:,0],MP0[:,1],color='gray')
Dip_Midpoint0 = ax[3].scatter(MP0[:,0],MP0[:,1],color='black')
Electrodes0 = ax[3].scatter(EL0[:,0],EL0[:,1],color='red')
ax[3].legend([Dip_Midpoint0,Electrodes0], ["Dipole Midpoint", "Electrodes"],scatterpoints=1)
ax[4].plot(MP1[:,0],MP1[:,1],color='gray')
Dip_Midpoint1 = ax[4].scatter(MP1[:,0],MP1[:,1],color='black')
Electrodes1 = ax[4].scatter(EL1[:,0],EL1[:,1],color='red')
ax[4].legend([Dip_Midpoint1,Electrodes1], ["Dipole Midpoint", "Electrodes"],scatterpoints=1)
return ax
ax3.plot(MP0[:,0],MP0[:,1],color='gray')
Dip_Midpoint0 = ax3.scatter(MP0[:,0],MP0[:,1],color='black')
Electrodes0 = ax3.scatter(EL0[:,0],EL0[:,1],color='red')
ax3.legend([Dip_Midpoint0,Electrodes0], ["Dipole Midpoint", "Electrodes"],scatterpoints=1)
ax4.plot(MP1[:,0],MP1[:,1],color='gray')
Dip_Midpoint1 = ax4.scatter(MP1[:,0],MP1[:,1],color='black')
Electrodes1 = ax4.scatter(EL1[:,0],EL1[:,1],color='red')
ax4.legend([Dip_Midpoint1,Electrodes1], ["Dipole Midpoint", "Electrodes"],scatterpoints=1)
return fig
#Function to visualise and compare any two meaningful plots for the sphere in a uniform backgound with an unifom Electric Field
def interact_conductiveSphere(R,log_sig0,log_sig1,Figure1a,Figure1b,Figure2a,Figure2b):
sig0,sig1 = conductivity_log_wrapper(log_sig0,log_sig1)
E0 = 1. # inducing field strength in V/m
n = 100 #level of discretisation
@@ -650,45 +648,45 @@ def interact_conductiveSphere(R,log_sig0,log_sig1,Figure1a,Figure1b,Figure2a,Fig
fig, ax = plt.subplots(1,2,figsize=(18,6))
#Setup figure 1 with options Configuration, Total or Secondary,
#Setup figure 1 with options Configuration, Total or Secondary,
#then Potential, ElectricField, Current Density or Charges Density
if Figure1a == 'Configuration':
ax[0] = get_Setup(XYZ,sig0,sig1,R,E0,ax[0],True,[0.1,0.1,0.6])
elif Figure1a == 'Total':
if Figure1b == 'Potential':
ax[0] = Plot_Total_Potential(XYZ,sig0,sig1,R,E0,ax[0])
elif Figure1b == 'ElectricField':
ax[0] = Plot_Total_ElectricField(XYZ,sig0,sig1,R,E0,ax[0])
elif Figure1b == 'CurrentDensity':
ax[0] = Plot_Total_Currents(XYZ,sig0,sig1,R,E0,ax[0])
elif Figure1b == 'ChargesDensity':
ax[0] = Plot_ChargesDensity(XYZ,sig0,sig1,R,E0,ax[0])
elif Figure1a == 'Secondary':
if Figure1b == 'Potential':
ax[0] = Plot_Secondary_Potential(XYZ,sig0,sig1,R,E0,ax[0])
elif Figure1b == 'ElectricField':
ax[0] = Plot_Secondary_ElectricField(XYZ,sig0,sig1,R,E0,ax[0])
elif Figure1b == 'CurrentDensity':
ax[0] = Plot_Secondary_Currents(XYZ,sig0,sig1,R,E0,ax[0])
elif Figure1b == 'ChargesDensity':
ax[0] = Plot_ChargesDensity(XYZ,sig0,sig1,R,E0,ax[0])
if Figure1a== 'Configuration':
ax[1] = Plot_Primary_Potential(XYZ,sig0,sig1,R,E0,ax[1])
print 'While figure1 is plotting Configuration, figure2 plots the primary field'
elif Figure2a == 'Total':
elif Figure2a == 'Total':
if Figure2b == 'Potential':
ax[1] = Plot_Total_Potential(XYZ,sig0,sig1,R,E0,ax[1])
@@ -701,8 +699,8 @@ def interact_conductiveSphere(R,log_sig0,log_sig1,Figure1a,Figure1b,Figure2a,Fig
elif Figure2b == 'ChargesDensity':
ax[1] = Plot_ChargesDensity(XYZ,sig0,sig1,R,E0,ax[1])
elif Figure2a == 'Secondary':
elif Figure2a == 'Secondary':
if Figure2b == 'Potential':
ax[1] = Plot_Secondary_Potential(XYZ,sig0,sig1,R,E0,ax[1])
@@ -717,10 +715,10 @@ def interact_conductiveSphere(R,log_sig0,log_sig1,Figure1a,Figure1b,Figure2a,Fig
plt.tight_layout(True)
plt.show()
#Interactive Visualisation of the responses of two configurations to a (pseudo) DC resistivity survey
def interactive_two_configurations_comparison(log_sig0,log_sig1,log_sig2,R0,R1,xstart,ystart,xend,yend,dipole_number,electrode_spacing,matching_spheres_example):
sig0,sig1 = conductivity_log_wrapper(log_sig0,log_sig1)
sig2 = 10.**log_sig2
E0 = 1. # inducing field strength in V/m
@@ -731,32 +729,21 @@ def interactive_two_configurations_comparison(log_sig0,log_sig1,log_sig2,R0,R1,x
XYZ = ndgrid(xr,yr,zr) # Space Definition
PlotOpt = 'Total'
#Initializing the figure
fig = plt.figure(figsize=(20,20))
ax0 = plt.subplot2grid((20,12), (0, 0),colspan=6,rowspan=6) #Configuration Conductive Sphere
ax1 = plt.subplot2grid((20,12), (0, 6),colspan=6,rowspan=6) #Configuration Resistive Sphere
ax2 = plt.subplot2grid((20,12), (16, 2), colspan=9,rowspan=4) # Data
ax3 = plt.subplot2grid((20,12), (8, 0),colspan=6,rowspan=6) #Potential Conductive Sphere
ax4 = plt.subplot2grid((20,12), (8, 6),colspan=6,rowspan=6) #Potential Resistive Potential
ax = [ax0,ax1,ax2,ax3,ax4]
if matching_spheres_example:
sig0 = 10.**(-3)
sig1 = 10.**(-2)
sig0 = 10.**(-3)
sig1 = 10.**(-2)
sig2 = 1.310344828 * 10**(-3)
R0 = 20.
R0 = 20.
R1 = 40.
two_configurations_comparison(XYZ,sig0,sig1,sig2,R0,R1,E0,xstart,ystart,xend,yend,dipole_number,electrode_spacing,PlotOpt,ax)
two_configurations_comparison(XYZ,sig0,sig1,sig2,R0,R1,E0,xstart,ystart,xend,yend,dipole_number,electrode_spacing,PlotOpt)
else:
two_configurations_comparison(XYZ,sig0,sig1,sig2,R0,R1,E0,xstart,ystart,xend,yend,dipole_number,electrode_spacing,PlotOpt,ax)
two_configurations_comparison(XYZ,sig0,sig1,sig2,R0,R1,E0,xstart,ystart,xend,yend,dipole_number,electrode_spacing,PlotOpt)
plt.tight_layout(True)
plt.show()
def run(plotIt=True):
sig0 = -3. # conductivity of the wholespace
sig1 = -1. # conductivity of the sphere
@@ -769,11 +756,6 @@ def run(plotIt=True):
zr = np.r_[0] # identical to saying `zr = np.array([0])`
XYZ = ndgrid(xr,yr,zr) # Space Definition
Vt,Vp,Vs = get_Potential(XYZ,sig0,sig1,R,E0)
Et, Ep, Es = get_ElectricField(XYZ,sig0,sig1,R,E0)
Jt,Jp,Js = get_Current(XYZ,sig0,sig1,R,Et,Ep,Es)
rho = get_ChargesDensity(XYZ,sig0,sig1,R,Et,Ep)
if plotIt:
fig, ax = plt.subplots(2,5,figsize=(50,10))
ax[0,0] = get_Setup(XYZ,sig0,sig1,R,E0,ax[0,0],True,[0.6,0.1,0.1])
@@ -786,11 +768,13 @@ def run(plotIt=True):
ax[1,3] = Plot_Secondary_Currents(XYZ,sig0,sig1,R,E0,ax[1,3])
ax[0,4] = Plot_Primary_Potential(XYZ,sig0,sig1,R,E0,ax[0,4])
ax[1,4] = Plot_ChargesDensity(XYZ,sig0,sig1,R,E0,ax[1,4])
else:
get_Potential(XYZ,sig0,sig1,R,E0) # This is so travis tests it
plt.show()
return Vt,Vp,Vs,Et,Ep,Es,Jt,Jp,Js,rho
if __name__ == '__main__':
run()