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

Analytic testing and showing online.

Browse Source
This commit is contained in:
rowanc1
2014-07-12 18:23:58 -05:00
parent 77906e2d57
commit adc2360a11
6 changed files with 80 additions and 49 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simpegDC/Examples/Verification.py
+57 -45
View File
@@ -1,51 +1,63 @@
from SimPEG import *
import simpegDC as DC
import matplotlib.pyplot as plt
cs = 25.
hx = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
hy = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
hz = [(cs,7, -1.3),(cs,20)]
mesh = Mesh.TensorMesh([hx, hy, hz], 'CCN')
sighalf = 1e-2
sigma = np.ones(mesh.nC)*sighalf
xtemp = np.linspace(-150, 150, 21)
ytemp = np.linspace(-150, 150, 21)
xyz_rxP = Utils.ndgrid(xtemp-10., ytemp, np.r_[0.])
xyz_rxN = Utils.ndgrid(xtemp+10., ytemp, np.r_[0.])
xyz_rxM = Utils.ndgrid(xtemp, ytemp, np.r_[0.])
fig, ax = plt.subplots(1,1, figsize = (5,5))
mesh.plotSlice(sigma, grid=True, ax = ax)
ax.plot(xyz_rxP[:,0],xyz_rxP[:,1], 'w.')
ax.plot(xyz_rxN[:,0],xyz_rxN[:,1], 'r.', ms = 3)
rx = DC.DipoleRx(xyz_rxP, xyz_rxN)
tx = DC.DipoleTx([-200, 0, -12.5],[+200, 0, -12.5], [rx])
survey = DC.SurveyDC([tx])
problem = DC.ProblemDC(mesh)
problem.pair(survey)
data = survey.dpred(sigma)
def DChalf(txlocP, txlocN, rxloc, sigma, I=1.):
rp = (txlocP.reshape([1,-1])).repeat(rxloc.shape[0], axis = 0)
rn = (txlocN.reshape([1,-1])).repeat(rxloc.shape[0], axis = 0)
rP = np.sqrt(((rxloc-rp)**2).sum(axis=1))
rN = np.sqrt(((rxloc-rn)**2).sum(axis=1))
return I/(sigma*2.*np.pi)*(1/rP-1/rN)
def run(plotIt=False):
cs = 25.
hx = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
hy = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
hz = [(cs,7, -1.3),(cs,20)]
mesh = Mesh.TensorMesh([hx, hy, hz], 'CCN')
sighalf = 1e-2
sigma = np.ones(mesh.nC)*sighalf
xtemp = np.linspace(-150, 150, 21)
ytemp = np.linspace(-150, 150, 21)
xyz_rxP = Utils.ndgrid(xtemp-10., ytemp, np.r_[0.])
xyz_rxN = Utils.ndgrid(xtemp+10., ytemp, np.r_[0.])
xyz_rxM = Utils.ndgrid(xtemp, ytemp, np.r_[0.])
data_analP = DChalf(np.r_[-200, 0, 0.],np.r_[+200, 0, 0.], xyz_rxP, sighalf)
data_analN = DChalf(np.r_[-200, 0, 0.],np.r_[+200, 0, 0.], xyz_rxN, sighalf)
data_anal = data_analP-data_analN
Data_anal = data_anal.reshape((21, 21), order = 'F')
Data = data.reshape((21, 21), order = 'F')
X = xyz_rxM[:,0].reshape((21, 21), order = 'F')
Y = xyz_rxM[:,1].reshape((21, 21), order = 'F')
# if plotIt:
# fig, ax = plt.subplots(1,1, figsize = (5,5))
# mesh.plotSlice(sigma, grid=True, ax = ax)
# ax.plot(xyz_rxP[:,0],xyz_rxP[:,1], 'w.')
# ax.plot(xyz_rxN[:,0],xyz_rxN[:,1], 'r.', ms = 3)
fig, ax = plt.subplots(1,2, figsize = (12, 5))
vmin = np.r_[data, data_anal].min()
vmax = np.r_[data, data_anal].max()
dat1 = ax[1].contourf(X, Y, Data, 60, vmin = vmin, vmax = vmax)
dat0 = ax[0].contourf(X, Y, Data_anal, 60, vmin = vmin, vmax = vmax)
cb0 = plt.colorbar(dat1, orientation = 'horizontal', ax = ax[0])
cb1 = plt.colorbar(dat1, orientation = 'horizontal', ax = ax[1])
ax[1].set_title('Analytic')
ax[0].set_title('Computed')
plt.show()
rx = DC.DipoleRx(xyz_rxP, xyz_rxN)
tx = DC.DipoleTx([-200, 0, -12.5],[+200, 0, -12.5], [rx])
survey = DC.SurveyDC([tx])
problem = DC.ProblemDC(mesh)
problem.pair(survey)
data = survey.dpred(sigma)
def DChalf(txlocP, txlocN, rxloc, sigma, I=1.):
rp = (txlocP.reshape([1,-1])).repeat(rxloc.shape[0], axis = 0)
rn = (txlocN.reshape([1,-1])).repeat(rxloc.shape[0], axis = 0)
rP = np.sqrt(((rxloc-rp)**2).sum(axis=1))
rN = np.sqrt(((rxloc-rn)**2).sum(axis=1))
return I/(sigma*2.*np.pi)*(1/rP-1/rN)
data_analP = DChalf(np.r_[-200, 0, 0.],np.r_[+200, 0, 0.], xyz_rxP, sighalf)
data_analN = DChalf(np.r_[-200, 0, 0.],np.r_[+200, 0, 0.], xyz_rxN, sighalf)
data_anal = data_analP-data_analN
Data_anal = data_anal.reshape((21, 21), order = 'F')
Data = data.reshape((21, 21), order = 'F')
X = xyz_rxM[:,0].reshape((21, 21), order = 'F')
Y = xyz_rxM[:,1].reshape((21, 21), order = 'F')
if plotIt:
fig, ax = plt.subplots(1,2, figsize = (12, 5))
vmin = np.r_[data, data_anal].min()
vmax = np.r_[data, data_anal].max()
dat1 = ax[1].contourf(X, Y, Data, 60, vmin = vmin, vmax = vmax)
dat0 = ax[0].contourf(X, Y, Data_anal, 60, vmin = vmin, vmax = vmax)
cb0 = plt.colorbar(dat1, orientation = 'horizontal', ax = ax[0])
cb1 = plt.colorbar(dat1, orientation = 'horizontal', ax = ax[1])
ax[1].set_title('Analytic')
ax[0].set_title('Computed')
plt.show()
return np.linalg.norm(data-data_anal)/np.linalg.norm(data_anal)
if __name__ == '__main__':
print run(plotIt=True)
simpegDC/Examples/WennerArray.py
+1
View File
@@ -6,6 +6,7 @@ import matplotlib.pyplot as plt
def getTxList(nElecs, aSpacing, in2D=False, plotIt=False):
elocs = np.arange(0,aSpacing*nElecs,aSpacing)
elocs -= (nElecs*aSpacing - aSpacing)/2
space = 1
WENNER = np.zeros((0,),dtype=int)
for ii in range(nElecs):
simpegDC/Examples/__init__.py
+1
View File
@@ -1 +1,2 @@
import WennerArray
import Verification