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

Finish test with analytical answer. O(h^2) error pass

Browse Source 
Add function MagSphereFreeSpace in MagAnalytic. Same results as MagSphereAnaFunA but with z-observation flipped.
This commit is contained in:
D Fournier
2015-11-27 16:08:21 -08:00
parent 291cce9d52
commit e9f02e1d77
3 changed files with 175 additions and 58 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simpegPF/Dev/MAG3D_fwr_Test.py
+124 -55
View File
@@ -13,86 +13,155 @@ import matplotlib
from fwr_MAG_obs import fwr_MAG_obs
plt.close('all')
#%% Create survey
B = np.array(([90.,0.,50000.]))
B = np.array(([-45.,315.,50000.]))
M = np.array(([90.,0.]))
M = np.array(([-45.,315.]))
# Sphere radius
R = 0.25
# # Or create juste a plane grid
xr = np.linspace(-1./2., 1./2., 10)
yr = np.linspace(-1./2., 1./2., 10)
xr = np.linspace(-2., 2., 5)
yr = np.linspace(-2., 2., 5)
X, Y = np.meshgrid(xr, yr)
Z = np.ones((xr.size, yr.size))*.75
Z = np.ones((xr.size, yr.size)) * 2.5
rxLoc = np.c_[Utils.mkvc(X), Utils.mkvc(Y), Utils.mkvc(Z)]
ndata = rxLoc.shape[0]
#%%
d_iter = 4
lrl = np.zeros(d_iter)
#%% Loop through decreasing meshes and measure the residual
# Create mesh using simpeg and write out in GIF format
nc = 30.
hxind = [(1./nc, nc)]
hyind = [(1./nc, nc)]
hzind = [(1./nc, nc)]
for ii in range(d_iter):
nc = 3**(ii+1)
hxind = [(1./nc, nc)]
hyind = [(1./nc, nc)]
hzind = [(1./nc, nc)]
mesh = Mesh.TensorMesh([hxind, hyind, hzind], 'CCC')
xn = mesh.vectorNx
yn = mesh.vectorNy
zn = mesh.vectorNz
mcell = mesh.nC
sph_ind = PF.MagAnalytics.spheremodel(mesh, 0, 0, 0, R)
chibkg = 0.
chiblk = 0.01
model = np.ones(mcell)*chibkg
model[sph_ind] = chiblk
#%% Forward mode ldata
d = fwr_MAG_obs(mesh,B,M,rxLoc,model)
fwr_x = mkvc(d[0,:])
fwr_y = mkvc(d[1,:])
fwr_z = mkvc(d[2,:])
#%% Get the analystical answer and compute the residual
bxa,bya,bza = PF.MagAnalytics.MagSphereAnaFunA(rxLoc[:,0],rxLoc[:,1],rxLoc[:,2],R,0.,0.,0.,chiblk, np.array(([0.,0.,B[2]])),'secondary')
Bd = (450.-float(B[1]))%360.
Bi = B[0]; # Convert dip to horizontal to cartesian
Bx = np.cos(np.deg2rad(Bi)) * np.cos(np.deg2rad(Bd)) * B[2]
By = np.cos(np.deg2rad(Bi)) * np.sin(np.deg2rad(Bd)) * B[2]
Bz = np.sin(np.deg2rad(Bi)) * B[2]
Bo = np.c_[Bx, By, Bz]
bxa,bya,bza = PF.MagAnalytics.MagSphereFreeSpace(rxLoc[:,0],rxLoc[:,1],rxLoc[:,2],R,0.,0.,0.,chiblk, Bo)
#bxa,bya,bza = PF.MagAnalytics.MagSphereAnaFunA(rxLoc[:,0],rxLoc[:,1],rxLoc[:,2],R,0.,0.,0.,chiblk, np.array(([0.,0.,B[2]])),'secondary')
r_Bx = fwr_x - bxa
r_By = fwr_y - bya
r_Bz = fwr_z - bza
lrl[ii] = sum( r_Bx**2 + r_By**2 + r_Bz**2 ) **0.5
mesh = Mesh.TensorMesh([hxind, hyind, hzind], 'CCC')
xn = mesh.vectorNx
yn = mesh.vectorNy
zn = mesh.vectorNz
mcell = mesh.nC
sph_ind = PF.MagAnalytics.spheremodel(mesh, 0, 0, 0, 0.25)
chibkg = 0.
chiblk = 0.01
model = np.ones(mcell)*chibkg
model[sph_ind] = chiblk
#%% Forward mode ldata
d = fwr_MAG_obs(mesh,B,M,rxLoc,model)
#%% Get the analystical answer and compute the residual
bxa,bya,bza = PF.MagAnalytics.MagSphereAnaFunA(rxLoc[:,0],rxLoc[:,1],rxLoc[:,2],.25,0.,0.,0.,chiblk, np.array(([0.,0.,B[2]])),'secondary')
r_Bz = mkvc(d) - bza
lrl = sum( r_Bz**2 ) **0.5
print "Residual between analytical and integral= " + str(lrl)
#%% Plot results
for ii in range(d_iter):
nc = 3**(ii+1)
print "Residual= " + str(lrl[ii]) + "\t dx= " + str(1./nc)
#%% Plot fields
plt.figure(1)
ax = plt.subplot(131)
plt.imshow(np.reshape(bxa,X.shape), interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()])
plt.contour(X,Y, np.reshape(bxa,X.shape),10)
plt.scatter(X,Y, c='k', s=5)
plt.imshow(np.reshape(bxa,X.shape).T, interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()], origin = 'lower')
plt.colorbar()
plt.contour(X,Y, np.reshape(bxa,X.shape).T,10)
plt.scatter(X,Y, c=np.reshape(bxa,X.shape).T, s=20)
ax.set_title('Sphere Ana Bx')
ax = plt.subplot(132)
plt.imshow(np.reshape(bya,X.shape), interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()])
plt.contour(X,Y, np.reshape(bya,X.shape),10)
plt.scatter(X,Y, c='k', s=5)
plt.imshow(np.reshape(bya,X.shape).T, interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()], origin = 'lower')
plt.colorbar()
plt.contour(X,Y, np.reshape(bya,X.shape).T,10)
plt.scatter(X,Y, c=np.reshape(bya,X.shape).T, s=20)
ax.set_title('Sphere Ana By')
ax = plt.subplot(133)
plt.imshow(np.reshape(bza,X.shape), interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()])
plt.contour(X,Y, np.reshape(bza,X.shape),10)
plt.scatter(X,Y, c='k', s=5)
plt.imshow(np.reshape(bza,X.shape).T, interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()], origin = 'lower')
plt.colorbar()
plt.contour(X,Y, np.reshape(bza,X.shape).T,10)
plt.scatter(X,Y, c=np.reshape(bza,X.shape).T, s=20)
ax.set_title('Sphere Ana Bz')
#%% Plot foward data
#%% Plot the forward solution from integral
plt.figure(2)
plt.subplot(121)
d2D = np.reshape(d,X.shape)
plt.imshow(d2D, interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()])
plt.contour(X,Y, d2D,10)
plt.scatter(X,Y, c='k', s=5)
ax = plt.subplot(131)
plt.imshow(np.reshape(fwr_x,X.shape).T, interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max() ], origin = 'lower')
plt.colorbar()
plt.contour(X,Y, np.reshape(fwr_x,X.shape).T,10)
plt.scatter(X,Y, c=np.reshape(fwr_x,X.shape).T, s=20)
ax.set_title('Sphere Ana Bx')
#%% Compare fields
ax = plt.subplot(132)
plt.imshow(np.reshape(fwr_y,X.shape).T, interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()], origin = 'lower')
plt.colorbar()
plt.contour(X,Y, np.reshape(fwr_y,X.shape).T,10)
plt.scatter(X,Y, c=np.reshape(fwr_y,X.shape).T, s=20)
ax.set_title('Sphere Ana By')
plt.subplot(122)
plt.imshow(np.reshape(r_Bz,X.shape), interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()])
plt.contour(X,Y, np.reshape(dBz,X.shape),10)
plt.scatter(X,Y, c='k', s=5)
ax = plt.subplot(133)
plt.imshow(np.reshape(fwr_z,X.shape).T, interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()], origin = 'lower')
plt.colorbar()
plt.contour(X,Y, np.reshape(fwr_z,X.shape).T,10)
plt.scatter(X,Y, c=np.reshape(fwr_z,X.shape).T, s=20)
ax.set_title('Sphere Ana Bz')
#%% Plot foward data
plt.figure(3)
ax = plt.subplot(131)
plt.imshow(np.reshape(r_Bx,X.shape).T, interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()], origin = 'lower')
plt.colorbar()
plt.contour(X,Y, np.reshape(r_Bx,X.shape).T,10)
plt.scatter(X,Y, c=np.reshape(r_Bx,X.shape).T, s=20)
ax.set_title('Sphere Ana Bx')
ax = plt.subplot(132)
plt.imshow(np.reshape(r_By,X.shape).T, interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()], origin = 'lower')
plt.colorbar()
plt.contour(X,Y, np.reshape(r_By,X.shape).T,10)
plt.scatter(X,Y, c=np.reshape(r_By,X.shape).T, s=20)
ax.set_title('Sphere Ana By')
ax = plt.subplot(133)
plt.imshow(np.reshape(r_Bz,X.shape).T, interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()], origin = 'lower')
plt.colorbar()
plt.contour(X,Y, np.reshape(r_Bz,X.shape).T,10)
plt.scatter(X,Y, c=np.reshape(r_Bz,X.shape).T, s=20)
ax.set_title('Sphere Ana Bz')
simpegPF/Dev/fwr_MAG_obs.py
+4 -3
View File
@@ -53,7 +53,7 @@ def fwr_MAG_obs(mesh,B,M,rxLoc,model):
Ptmi = mkvc(np.r_[np.cos(np.deg2rad(B[0]))*np.cos(np.deg2rad(D)),np.cos(np.deg2rad(B[0]))*np.sin(np.deg2rad(D)),np.sin(np.deg2rad(B[0]))],2).T;
d = np.zeros((ndata,1))
d = np.zeros((3,ndata))
# Loop through all observations and create forward operator (ndata-by-mcell)
print "Begin forward modeling " +str(int(ndata)) + " data points..."
@@ -64,10 +64,11 @@ def fwr_MAG_obs(mesh,B,M,rxLoc,model):
tx, ty, tz = get_T_mat(xn,yn,zn,rxLoc[ii,:])
G = Ptmi.dot(np.vstack((tx,ty,tz)))*Mxyz
# G = Ptmi.dot(np.vstack((tx,ty,tz)))*Mxyz
Gxyz = np.vstack((tx,ty,tz))*Mxyz
#%%
# Forward operator
d[ii,0] = G.dot(model)
d[:,ii] = Gxyz.dot(model)
d_iter = np.floor(float(ii)/float(ndata)*10.);