mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-28 23:58:09 +08:00
Clean up examples
This commit is contained in:
D Fournier
@@ -1,268 +0,0 @@
|
||||
#%%
|
||||
from SimPEG import *
|
||||
import matplotlib.pyplot as plt
|
||||
import simpegPF as PF
|
||||
import scipy.interpolate as interpolation
|
||||
import time
|
||||
#from interpFFT import interpFFT
|
||||
#from fwr_MAG_data import fwr_MAG_data
|
||||
|
||||
import os
|
||||
|
||||
home_dir = 'C:\\LC\\Private\\dominiquef\\Projects\\4414_Minsim\\Modeling\\MAG'
|
||||
|
||||
os.chdir(home_dir)
|
||||
|
||||
plt.close('all')
|
||||
|
||||
topofile = 'Gaussian.topo'
|
||||
|
||||
zoffset = 5
|
||||
#%% Create survey
|
||||
# Load in topofile or create flat surface
|
||||
if not topofile:
|
||||
|
||||
actv = np.ones(mesh.nC)
|
||||
|
||||
else:
|
||||
topo = np.genfromtxt(topofile,skip_header=1)
|
||||
|
||||
|
||||
B = np.array(([90.,0.,50000.]))
|
||||
|
||||
M = np.array(([90.,0.,315.]))
|
||||
|
||||
# Sphere radius
|
||||
R = 25.
|
||||
|
||||
sclx = 100.
|
||||
dx = 5.
|
||||
#%% Loop through decreasing meshes and measure the residual
|
||||
# Create mesh using simpeg and write out in GIF format
|
||||
|
||||
# # Or create juste a plane grid
|
||||
xr = np.linspace(-102.5, 97.5, 41)
|
||||
yr = np.linspace(-52.5, 47.5, 21)
|
||||
X, Y = np.meshgrid(xr, yr)
|
||||
|
||||
nc = int(sclx/dx)
|
||||
|
||||
hxind = [(dx, 2*nc)]
|
||||
hyind = [(dx, nc)]
|
||||
hzind = [(dx, nc)]
|
||||
|
||||
mesh = Mesh.TensorMesh([hxind, hyind, hzind], 'CCN')
|
||||
|
||||
actv = PF.Magnetics.getActiveTopo(mesh,topo,'N')
|
||||
|
||||
# Drape observations on topo + offset
|
||||
if not topofile:
|
||||
Z = np.ones((xr.size, yr.size)) * 2.5
|
||||
|
||||
else:
|
||||
F = interpolation.interp2d(topo[:,0],topo[:,1],topo[:,2])
|
||||
#F = interpolation.NearestNDInterpolator(topo[:,0:2],topo[:,2])
|
||||
Z = F(xr,yr) + zoffset
|
||||
|
||||
rxLoc = np.c_[Utils.mkvc(X.T), Utils.mkvc(Y.T), Utils.mkvc(Z.T)]
|
||||
|
||||
ndata = rxLoc.shape[0]
|
||||
|
||||
xn = mesh.vectorNx
|
||||
yn = mesh.vectorNy
|
||||
zn = mesh.vectorNz
|
||||
|
||||
mcell = mesh.nC
|
||||
|
||||
print 'Mesh size: ' + str(mcell)
|
||||
|
||||
#%% Create model
|
||||
chibkg = 0.0001
|
||||
chiblk = 0.01
|
||||
model = np.ones(mcell)*chibkg
|
||||
|
||||
# Do a three sphere problem for more frequencies
|
||||
sph_ind = PF.MagAnalytics.spheremodel(mesh, 0., 0., -sclx/3, R)
|
||||
model[sph_ind] = 0.5*chiblk
|
||||
|
||||
sph_ind = PF.MagAnalytics.spheremodel(mesh, -sclx/2., 0., -sclx/3., R/3.)
|
||||
model[sph_ind] = 4.*chiblk
|
||||
|
||||
sph_ind = PF.MagAnalytics.spheremodel(mesh, sclx/2., 0., -sclx/2.5, R/2.5)
|
||||
model[sph_ind] = 2.5*chiblk
|
||||
|
||||
# Zero out
|
||||
model[actv==0] = -100
|
||||
Utils.writeUBCTensorMesh('Mesh.msh',mesh)
|
||||
Utils.writeUBCTensorModel('Model.sus',mesh,model)
|
||||
Utils.writeUBCTensorModel('nullcell.dat',mesh,actv)
|
||||
#actv = np.ones(mesh.nC)
|
||||
#%% Forward mode ldata
|
||||
d = PF.Magnetics.Intgrl_Fwr_Data(mesh,B,M,rxLoc,model,actv,'tmi')
|
||||
#fwr_tmi = d[0:ndata]
|
||||
#fwr_y = d[ndata:2*ndata]
|
||||
#fwr_z = d[2*ndata:]
|
||||
|
||||
|
||||
#%% Compute data on a line
|
||||
xx = np.linspace(xr.min(), xr.max(), 200)
|
||||
yy = np.zeros(len(xx))
|
||||
zz = F(xx,0.) + zoffset
|
||||
|
||||
rxLoc = np.c_[xx,yy,zz]
|
||||
|
||||
d_line = PF.Magnetics.Intgrl_Fwr_Data(mesh,B,M,rxLoc,model,actv,'tmi')
|
||||
|
||||
d_iter = 4
|
||||
|
||||
l1_r = np.zeros((d_iter,5))
|
||||
l2_r = np.zeros((d_iter,5))
|
||||
linf_r = np.zeros((d_iter,5))
|
||||
timer = np.zeros((d_iter,5))
|
||||
|
||||
d2d = np.reshape(d, (len(yr),len(xr)))
|
||||
#%% Try different interpolation schemes
|
||||
for ii in range(d_iter):
|
||||
|
||||
indx = ii+1
|
||||
|
||||
|
||||
dsub = d2d[::indx,::indx]
|
||||
|
||||
xsub = xr[::indx]
|
||||
ysub = yr[::indx]
|
||||
|
||||
# Nearest Neighbourg
|
||||
start_time = time.time()
|
||||
F = interpolation.NearestNDInterpolator(np.c_[mkvc(Y[::indx,::indx].T),mkvc(X[::indx,::indx].T)],mkvc(dsub.T))
|
||||
d_i2d_nnb = mkvc( F(0.,xx) )
|
||||
l1_r[ii,0] = np.sum( np.abs(d_line - d_i2d_nnb) )**0.5
|
||||
l2_r[ii,0] = np.sum( (d_line - d_i2d_nnb)**2. )
|
||||
linf_r[ii,0] = np.max( np.abs(d_line - d_i2d_nnb) )
|
||||
timer[ii,0] = (time.time() - start_time)
|
||||
|
||||
# Linear interpolation
|
||||
start_time = time.time()
|
||||
F = interpolation.interp2d(ysub,xsub,mkvc(dsub.T))
|
||||
d_i2d_lin = mkvc( F(0.,xx) )
|
||||
l1_r[ii,1] = np.sum( np.abs(d_line - d_i2d_lin) )**0.5
|
||||
l2_r[ii,1] = np.sum( (d_line - d_i2d_lin)**2. )
|
||||
linf_r[ii,1] = np.max( np.abs(d_line - d_i2d_lin) )
|
||||
timer[ii,1] = (time.time() - start_time)
|
||||
|
||||
|
||||
# Cubic interpolation
|
||||
start_time = time.time()
|
||||
F = interpolation.interp2d(ysub,xsub,mkvc(dsub.T),kind='cubic')
|
||||
d_i2d_cub = mkvc( F(0.,xx) )
|
||||
l1_r[ii,2] = np.sum( np.abs(d_line - d_i2d_cub) )**0.5
|
||||
l2_r[ii,2] = np.sum( (d_line - d_i2d_cub)**2. )
|
||||
linf_r[ii,2] = np.max( np.abs(d_line - d_i2d_cub) )
|
||||
timer[ii,2] = (time.time() - start_time)
|
||||
|
||||
|
||||
# Quintic interpolation
|
||||
start_time = time.time()
|
||||
F = interpolation.interp2d(ysub,xsub,mkvc(dsub.T),kind='quintic')
|
||||
d_i2d_qui = mkvc( F(0.,xx) )
|
||||
l1_r[ii,3] = np.sum( np.abs(d_line - d_i2d_qui) )**0.5
|
||||
l2_r[ii,3] = np.sum( (d_line - d_i2d_qui)**2. )
|
||||
linf_r[ii,3] = np.max( np.abs(d_line - d_i2d_qui) )
|
||||
timer[ii,3] = (time.time() - start_time)
|
||||
|
||||
# CloughTocher interpolation
|
||||
start_time = time.time()
|
||||
F = interpolation.CloughTocher2DInterpolator(np.c_[mkvc(Y[::indx,::indx].T),mkvc(X[::indx,::indx].T)],mkvc(dsub.T))
|
||||
d_i2d_CTI = mkvc( F(0.,xx) )
|
||||
l1_r[ii,4] = np.sum( np.abs(d_line - d_i2d_CTI) )**0.5
|
||||
l2_r[ii,4] = np.sum( (d_line - d_i2d_CTI)**2. )
|
||||
linf_r[ii,4] = np.max( np.abs(d_line - d_i2d_CTI) )
|
||||
timer[ii,4] = (time.time() - start_time)
|
||||
|
||||
# Minimum curvature
|
||||
|
||||
|
||||
#==============================================================================
|
||||
# #%% FFT interpolation
|
||||
# d2d_out = interpFFT(xsub,ysub,dsub)
|
||||
#
|
||||
# # Create new distance vector
|
||||
# XX = np.linspace(np.min(xsub),np.max(xsub),d2d_out.shape[1])
|
||||
# YY = np.linspace(np.min(ysub),np.max(ysub),d2d_out.shape[0])
|
||||
#
|
||||
# start_time = time.time()
|
||||
# F = interpolation.interp2d(XX,YY,d2d_out)
|
||||
# d_i2d_fft = mkvc( F(xx,0.) )
|
||||
# l1_r[ii,4] = np.sum( np.abs(d_line - d_i2d_nnb) )**0.5
|
||||
# l2_r[ii,4] = np.sum( (d_line - d_i2d_nnb)**2. )
|
||||
# linf_r[ii,4] = np.max( np.abs(d_line - d_i2d_nnb) )
|
||||
# timer[ii,4] = (time.time() - start_time)
|
||||
#
|
||||
# print("--- FFT completed in %s seconds ---" % (time.time() - start_time))
|
||||
#
|
||||
#
|
||||
# plt.figure()
|
||||
# ax = plt.subplot()
|
||||
# plt.imshow(d2d_out, interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max() ], origin = 'lower')
|
||||
# plt.colorbar(fraction=0.04)
|
||||
#==============================================================================
|
||||
|
||||
|
||||
#%% Write predicted to file
|
||||
|
||||
#PF.Magnetics.writeUBCobs('Obsloc.loc',B,M,rxLoc,d,np.ones(len(d)))
|
||||
|
||||
#%% Plot results
|
||||
#==============================================================================
|
||||
# print 'Residual between analytical sphere and integral forward'
|
||||
# for ii in range(d_iter):
|
||||
# nc = 3**(ii+1)
|
||||
#
|
||||
# print "||r||= " + str(lrl[ii]) + "\t dx= " + str(1./nc)
|
||||
#==============================================================================
|
||||
|
||||
|
||||
#%% Plot the forward solution from integral
|
||||
plt.figure(figsize=[8,4])
|
||||
ax = plt.subplot()
|
||||
plt.imshow(d2d, interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max() ], origin = 'lower')
|
||||
plt.colorbar(fraction=0.02)
|
||||
plt.contour(X,Y, d2d,10)
|
||||
plt.scatter(X,Y, s=5)
|
||||
|
||||
|
||||
plt.figure(figsize=[8,4])
|
||||
plt.contour(X,Y, np.reshape(d,X.shape),10)
|
||||
plt.scatter(X,Y, c=np.reshape(d,X.shape), s=10)
|
||||
plt.scatter(xx,yy, c='k', s=20, marker='o')
|
||||
ax.set_title('Numerical')
|
||||
|
||||
|
||||
#%%
|
||||
|
||||
plt.figure(figsize=[7,9])
|
||||
ax = plt.subplot()
|
||||
plt.plot(xx,d_line,c='r', linewidth=3)
|
||||
plt.plot(xx,d_i2d_lin)
|
||||
plt.plot(xx,d_i2d_cub)
|
||||
plt.plot(xx,d_i2d_qui)
|
||||
plt.plot(xx,d_i2d_nnb)
|
||||
plt.plot(xx,d_i2d_CTI)
|
||||
plt.xlim(xx.min(),xx.max())
|
||||
|
||||
plt.legend(['True','linear','Cubic','Quintic','NearestN','CloughTorcher'],bbox_to_anchor=(0.75, 0.25))
|
||||
# Plot interpolation from true value on line
|
||||
F = interpolation.interp1d(xx,d_line)
|
||||
dtrue = F(xr[::indx])
|
||||
plt.plot(xr[::indx],dtrue,c='r',linewidth=0.,marker='o')
|
||||
ax.set_title('Interpolated data profile')
|
||||
|
||||
#%% Write result to file
|
||||
with file('Interp_residual.dat','w') as fid:
|
||||
fid.write('NearestN\tLinear\tCubic\tQuintic\tCloughTocher\n')
|
||||
fid.write('\nL2-norm\n')
|
||||
np.savetxt(fid, l2_r, fmt='%5.3e',delimiter='\t',newline='\n')
|
||||
fid.write('\nL1-norm\n')
|
||||
np.savetxt(fid, l1_r, fmt='%5.3e',delimiter='\t',newline='\n')
|
||||
fid.write('\nLinf-norm\n')
|
||||
np.savetxt(fid, linf_r, fmt='%5.3e',delimiter='\t',newline='\n')
|
||||
@@ -1,282 +0,0 @@
|
||||
"""
|
||||
Intgrl_MAG_Fwr_Prism_Z_test.py : Test the discretization error for a semi
|
||||
infinite prism extending from the surface.
|
||||
|
||||
The test runs a sequence of forward calculations of magnetic over a
|
||||
grid of points at varying height. The goal is to find a relation between
|
||||
observation height and mesh cell size.
|
||||
|
||||
Since the integral equation is exact for a rectangular prism, the response
|
||||
is first calculated for a semi-infinite prism aligned with the mesh with
|
||||
inducing field rotated 45d. The prism and data location is then rotate by
|
||||
-45d so prism is rotated with respect to base mesh.
|
||||
|
||||
"""
|
||||
|
||||
#%%
|
||||
from SimPEG import *
|
||||
import matplotlib.pyplot as plt
|
||||
import simpegPF as PF
|
||||
import scipy.interpolate as interpolation
|
||||
import time
|
||||
import os
|
||||
|
||||
home_dir = 'C:\\LC\\Private\\dominiquef\\Projects\\4414_Minsim\\Modeling\\MAG'
|
||||
|
||||
os.chdir(home_dir)
|
||||
|
||||
#from fwr_MAG_data import fwr_MAG_data
|
||||
|
||||
plt.close('all')
|
||||
|
||||
topofile = []
|
||||
|
||||
zoffset = np.array([6., 5. ,3. , 2., 1.,0.5])
|
||||
|
||||
#%% Create survey
|
||||
|
||||
|
||||
B = np.array(([90.,0.,50000.]))
|
||||
|
||||
M = np.array(([90.,0.,315.]))
|
||||
|
||||
# Block side lenght
|
||||
R = 20.
|
||||
|
||||
# # Or create juste a plane grid
|
||||
xr = np.linspace(-19.5, 19.5, 40)
|
||||
yr = np.linspace(-19., 19., 40)
|
||||
X, Y = np.meshgrid(xr, yr)
|
||||
|
||||
|
||||
|
||||
sclx = 90.
|
||||
dx = np.asarray([5., 4. ,3. , 2., 1.])
|
||||
|
||||
d_iter = len(dx)
|
||||
l1_r = np.zeros(d_iter)
|
||||
l2_r = np.zeros(d_iter)
|
||||
linf_r = np.zeros(d_iter)
|
||||
timer = np.zeros(d_iter)
|
||||
mcell = np.zeros(d_iter)
|
||||
#%% First Loop through the observation heights and compute the true response
|
||||
nc = int(sclx/30.)
|
||||
|
||||
hxind = [(30., nc)]
|
||||
hyind = [(30., nc)]
|
||||
hzind = [(1., 1)]
|
||||
|
||||
mesh = Mesh.TensorMesh([hxind, hyind, hzind], 'CCN')
|
||||
|
||||
actv = np.ones(mesh.nC)
|
||||
|
||||
|
||||
# Pre-allocate true data
|
||||
obs_x = np.zeros((X.size,zoffset.size))
|
||||
obs_y = np.zeros((X.size,zoffset.size))
|
||||
obs_z = np.zeros((X.size,zoffset.size))
|
||||
|
||||
plt.figure(1)
|
||||
for ii in range(zoffset.size):
|
||||
|
||||
# Drape observations on topo + offset
|
||||
if not topofile:
|
||||
Z = np.ones((xr.size, yr.size)) * zoffset[ii]
|
||||
|
||||
else:
|
||||
F = interpolation.NearestNDInterpolator(topo[:,0:2],topo[:,2])
|
||||
Z = F(X,Y) + zoffset[ii]
|
||||
|
||||
rxLoc = np.c_[Utils.mkvc(X.T), Utils.mkvc(Y.T), Utils.mkvc(Z.T)]
|
||||
|
||||
ndata = rxLoc.shape[0]
|
||||
|
||||
print 'Mesh size: ' + str(mesh.nC)
|
||||
|
||||
#%% Create model
|
||||
chibkg = 0.
|
||||
chiblk = 0.01
|
||||
model = np.ones((mesh.nCx,mesh.nCy,mesh.nCz))*chibkg
|
||||
|
||||
# Do a three sphere problem for more frequencies
|
||||
model[1:2,1:2,:] = chiblk
|
||||
model = mkvc(model)
|
||||
|
||||
Utils.writeUBCTensorMesh('Mesh.msh',mesh)
|
||||
Utils.writeUBCTensorModel('Model.sus',mesh,model)
|
||||
#actv = np.ones(mesh.nC)
|
||||
#%% Forward mode ldata
|
||||
|
||||
temp = PF.Magnetics.Intgrl_Fwr_Data(mesh,B,M,rxLoc,model,actv,'xyz')
|
||||
|
||||
obs_x[:,ii] = temp[0:ndata]
|
||||
obs_y[:,ii] = temp[ndata:2*ndata]
|
||||
obs_z[:,ii] = temp[2*ndata:]
|
||||
|
||||
#%% Plot fields
|
||||
|
||||
ax = plt.subplot(3,zoffset.size,ii+1)
|
||||
plt.imshow(np.reshape(obs_x[:,ii],X.shape), extent=[xr.min(), xr.max(), yr.min(), yr.max()], origin = 'lower')
|
||||
plt.colorbar(fraction=0.02)
|
||||
plt.contour(X,Y, np.reshape(obs_x[:,ii],X.shape),10)
|
||||
#plt.scatter(X,Y, c='k', s=10)
|
||||
|
||||
ax = plt.subplot(3,zoffset.size,zoffset.size+ii+1)
|
||||
plt.imshow(np.reshape(obs_y[:,ii],X.shape), extent=[xr.min(), xr.max(), yr.min(), yr.max()], origin = 'lower')
|
||||
plt.colorbar(fraction=0.02)
|
||||
plt.contour(X,Y, np.reshape(obs_y[:,ii],X.shape),10)
|
||||
#plt.scatter(X,Y, c='k', s=10)
|
||||
|
||||
ax = plt.subplot(3,zoffset.size,2*zoffset.size + ii+1)
|
||||
plt.imshow(np.reshape(obs_z[:,ii],X.shape), extent=[xr.min(), xr.max(), yr.min(), yr.max()], origin = 'lower')
|
||||
plt.colorbar(fraction=0.02)
|
||||
plt.contour(X,Y, np.reshape(obs_z[:,ii],X.shape),10)
|
||||
#plt.scatter(X,Y, c='k', s=10)
|
||||
|
||||
#%% Rotate the data and discretize the rotated block on various cell size
|
||||
|
||||
Rz = np.array([[np.cos(np.deg2rad(45)), -np.sin(np.deg2rad(45))],
|
||||
[np.sin(np.deg2rad(45)), np.cos(np.deg2rad(45))]])
|
||||
|
||||
temp = np.vstack([rxLoc[:,0].T,rxLoc[:,1].T])
|
||||
|
||||
ROTxy = np.dot(Rz,temp)
|
||||
|
||||
|
||||
|
||||
|
||||
#%%
|
||||
|
||||
blocksurf = ['Block_0.ts',True,True]
|
||||
|
||||
B = np.array(([90.,0.,50000.]))
|
||||
|
||||
M = np.array(([90.,0.,315.]))
|
||||
|
||||
l2 = np.zeros((dx.size,zoffset.size))
|
||||
linf = np.zeros((dx.size,zoffset.size))
|
||||
|
||||
plt.figure()
|
||||
for ii in range(dx.size):
|
||||
|
||||
nc = int(sclx/dx[ii])
|
||||
|
||||
hxind = [(dx[ii], nc)]
|
||||
hyind = [(dx[ii], nc)]
|
||||
hzind = [(1., 1)]
|
||||
|
||||
mesh = Mesh.TensorMesh([hxind, hyind, hzind], 'CCN')
|
||||
|
||||
actv = np.ones(mesh.nC)
|
||||
|
||||
model = np.zeros(mesh.nC)
|
||||
|
||||
indx = PF.BaseMag.gocad2vtk(blocksurf[0],mesh, bcflag = blocksurf[1], inflag = blocksurf[2])
|
||||
|
||||
# Compensate for the volume difference
|
||||
scale = 1.#(30.**2.) / (len(indx)*dx[ii]**2.)
|
||||
model[indx] = chiblk * scale
|
||||
|
||||
Utils.writeUBCTensorMesh('Mesh_' + str(int(dx[ii])) + 'm.msh',mesh)
|
||||
Utils.writeUBCTensorModel('Model_' + str(int(dx[ii])) + 'm.sus',mesh,model)
|
||||
|
||||
|
||||
for jj in range(zoffset.size):
|
||||
|
||||
# Drape observations on topo + offset
|
||||
if not topofile:
|
||||
Z = np.ones((xr.size, yr.size)) * zoffset[jj]
|
||||
|
||||
else:
|
||||
F = interpolation.NearestNDInterpolator(topo[:,0:2],topo[:,2])
|
||||
Z = F(X,Y) + zoffset[jj]
|
||||
|
||||
|
||||
rxLoc_ROT = np.c_[ROTxy.T,mkvc(Z)]
|
||||
|
||||
temp = PF.Magnetics.Intgrl_Fwr_Data(mesh,B,M,rxLoc_ROT,model,actv,'xyz')
|
||||
|
||||
pre_x = temp[0:ndata]
|
||||
pre_y = temp[ndata:2*ndata]
|
||||
pre_z = temp[2*ndata:]
|
||||
|
||||
res_z = (obs_z[:,jj] - pre_z) / np.max(np.abs(pre_z))
|
||||
l2[ii,jj] = np.sum( (res_z) ** 2. )**0.5
|
||||
linf[ii,jj] = np.max( np.abs(res_z) )
|
||||
|
||||
ax = plt.subplot(dx.size,zoffset.size,(ii*(zoffset.size))+jj+1)
|
||||
plt.imshow(np.reshape(res_z,X.shape), extent=[xr.min(), xr.max(), yr.min(), yr.max()], origin = 'lower')
|
||||
plt.colorbar(fraction=0.02)
|
||||
plt.contour(X,Y, np.reshape(res_z,X.shape),10)
|
||||
plt.title('Z:' + str(zoffset[jj]) + ' dx:' + str(dx[ii]))
|
||||
ax.get_xaxis().set_ticks([])
|
||||
ax.get_yaxis().set_ticks([])
|
||||
#plt.scatter(rxLoc[:,0],rxLoc[:,1], c='k', s=10)
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
#%%
|
||||
plt.figure()
|
||||
plt.subplot(2,1,1)
|
||||
for ii in range(l2.shape[0]):
|
||||
plt.plot(zoffset,l2[ii,:])
|
||||
plt.xlabel('Z Heigth')
|
||||
plt.ylabel('L2-residual')
|
||||
plt.legend(dx.astype('str'))
|
||||
|
||||
plt.subplot(2,1,2)
|
||||
for ii in range(l2.shape[1]):
|
||||
plt.plot(dx,l2[:,ii])
|
||||
|
||||
plt.legend(zoffset.astype('str'))
|
||||
plt.xlabel('Cell Size')
|
||||
plt.ylabel('L2-residual')
|
||||
|
||||
|
||||
plt.figure()
|
||||
plt.subplot(2,1,1)
|
||||
for ii in range(linf.shape[0]):
|
||||
plt.plot(zoffset,linf[ii,:])
|
||||
plt.xlabel('Z Heigth')
|
||||
plt.ylabel('L_inf-residual')
|
||||
plt.legend(dx.astype('str'))
|
||||
|
||||
plt.subplot(2,1,2)
|
||||
for ii in range(linf.shape[1]):
|
||||
plt.plot(dx,linf[:,ii])
|
||||
|
||||
plt.legend(zoffset.astype('str'))
|
||||
plt.xlabel('Cell Size')
|
||||
plt.ylabel('Linf-residual')
|
||||
#%% Plot results
|
||||
#==============================================================================
|
||||
# print 'Residual between analytical sphere and integral forward'
|
||||
# print "dx \t nc \t l1 \t l2 \t linf \t Runtime"
|
||||
# for ii in range(d_iter):
|
||||
#
|
||||
# print str(dx[ii]) + "\t" + str(mcell[ii]) + "\t" + str(l1_r[ii]) + "\t" + str(l2_r[ii]) + "\t" + str(linf_r[ii]) + "\t" + str(timer[ii])
|
||||
#
|
||||
#==============================================================================
|
||||
|
||||
|
||||
#%% Plot the forward solution from integral
|
||||
#==============================================================================
|
||||
# plt.figure(2)
|
||||
# ax = plt.subplot()
|
||||
# plt.imshow(np.reshape(d,X.shape), interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max() ], origin = 'lower')
|
||||
# plt.colorbar(fraction=0.02)
|
||||
# plt.contour(X,Y, np.reshape(d,X.shape),10)
|
||||
# plt.scatter(X,Y, c=np.reshape(d,X.shape), s=20)
|
||||
# ax.set_title('Numerical')
|
||||
#
|
||||
# #%% Plot residual data
|
||||
# plt.figure(3)
|
||||
# ax = plt.subplot()
|
||||
# plt.imshow(np.reshape(r_tmi,X.shape), interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()], origin = 'lower')
|
||||
# plt.colorbar(fraction=0.02)
|
||||
# plt.contour(X,Y, np.reshape(r_tmi,X.shape),10)
|
||||
# plt.scatter(X,Y, c=np.reshape(r_tmi,X.shape), s=20)
|
||||
# ax.set_title('Sphere Ana Bx')
|
||||
#==============================================================================
|
||||
@@ -1,112 +0,0 @@
|
||||
import os
|
||||
|
||||
home_dir = 'C:\\LC\\Private\\dominiquef\\Projects\\4414_Minsim\\Modeling\\MAG\\Lalor'
|
||||
|
||||
inpfile = 'PYMAG3C_fwr.inp'
|
||||
|
||||
dsep = '\\'
|
||||
|
||||
os.chdir(home_dir)
|
||||
|
||||
#%%
|
||||
from SimPEG import np, sp, Utils, Mesh
|
||||
import scipy.interpolate as interpolation
|
||||
import simpegPF as PF
|
||||
import pylab as plt
|
||||
import time as tm
|
||||
|
||||
## New scripts to be added to basecode
|
||||
#from fwr_MAG_data import fwr_MAG_data
|
||||
#from read_MAGfwr_inp import read_MAGfwr_inp
|
||||
|
||||
#%%
|
||||
[mshfile, obsfile, modfile, magfile, topofile] = PF.BaseMag.read_MAGfwr_inp(inpfile)
|
||||
|
||||
# Load mesh file
|
||||
mesh = Utils.meshutils.readUBCTensorMesh(mshfile)
|
||||
|
||||
# Load in observation file
|
||||
#[B,M,dobs] = PF.BaseMag.readUBCmagObs(obsfile)
|
||||
|
||||
# Read in topo surface
|
||||
topsurf = "Topography_Local.ts"
|
||||
geosurf = [['_Top_Bollock_Rhyodacite_SUB.ts',False,True],
|
||||
['_Top_Bollock_Rhyodacite_SUB_B.ts',False,True],
|
||||
['_Top_Mafic_Intrusion.ts',True,True],
|
||||
['_Top_Mafic_Volcaniclastic.ts',False,True],
|
||||
['UnitA.ts',True,True],
|
||||
['UnitB.ts',True,True],
|
||||
['UnitC.ts',True,True],
|
||||
['UnitD.ts',True,True],
|
||||
['UnitE.ts',True,True],
|
||||
['UnitF.ts',True,True],
|
||||
]
|
||||
|
||||
vals = np.asarray([0.005,0.005,0.01,0.0025])
|
||||
|
||||
# Background density
|
||||
bkgr = 0.0001
|
||||
|
||||
# Offset data above topo
|
||||
zoffset = 2
|
||||
|
||||
|
||||
#%% Script starts here
|
||||
# # Create a grid of observations and offset the z from topo
|
||||
#xr = np.linspace(mesh.vectorCCx[0], mesh.vectorCCx[-1], 99)
|
||||
#yr = np.linspace(mesh.vectorCCy[0], mesh.vectorCCy[-1], 74)
|
||||
|
||||
xr = mesh.vectorCCx[::3]
|
||||
yr = mesh.vectorCCy[::3]
|
||||
X, Y = np.meshgrid(xr, yr)
|
||||
|
||||
topo = np.genfromtxt(topofile,skip_header=1)
|
||||
F = interpolation.NearestNDInterpolator(topo[:,0:2],topo[:,2])
|
||||
Z = F(X,Y) + zoffset
|
||||
rxLoc = np.c_[Utils.mkvc(X.T), Utils.mkvc(Y.T), Utils.mkvc(Z.T)]
|
||||
|
||||
ndata = rxLoc.shape[0]
|
||||
|
||||
B = np.array(([90.,0.,50000.]))
|
||||
|
||||
M = np.array(([90.,0.,315.]))
|
||||
|
||||
model= np.ones(mesh.nC) * bkgr
|
||||
# Load GOCAD surf
|
||||
#[vrtx, trgl] = PF.BaseMag.read_GOCAD_ts(tsfile)
|
||||
# Find active cells from surface
|
||||
|
||||
for ii in range(len(geosurf)):
|
||||
tin = tm.time()
|
||||
print "Computing indices with VTK: " + geosurf[ii][0]
|
||||
indx = PF.BaseMag.gocad2vtk(geosurf[ii][0],mesh, bcflag = geosurf[ii][1], inflag = geosurf[ii][2])
|
||||
print "VTK operation completed in " + str(tm.time() - tin)
|
||||
|
||||
model[indx] = vals[ii]
|
||||
|
||||
indx = PF.BaseMag.gocad2vtk(topsurf,mesh, bcflag = False, inflag = True)
|
||||
actv = np.zeros(mesh.nC)
|
||||
actv[indx] = 1
|
||||
|
||||
model[actv==0] = -100
|
||||
|
||||
Utils.meshutils.writeUBCTensorModel('VTKout.dat',mesh,model)
|
||||
|
||||
Utils.meshutils.writeUBCTensorMesh('Mesh_temp.msh',mesh)
|
||||
|
||||
start_time = tm.time()
|
||||
|
||||
d = PF.Magnetics.Intgrl_Fwr_Data(mesh,B,M,rxLoc,model,actv,'tmi')
|
||||
|
||||
timer = (tm.time() - start_time)
|
||||
|
||||
#%% Plot data
|
||||
plt.figure()
|
||||
ax = plt.subplot()
|
||||
plt.imshow(np.reshape(d,X.shape), interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()], origin = 'lower')
|
||||
plt.clim(0,25)
|
||||
plt.colorbar(fraction=0.02)
|
||||
plt.contour(X,Y, np.reshape(d,X.shape),10)
|
||||
plt.scatter(X,Y, c=np.reshape(d,X.shape), s=20)
|
||||
|
||||
ax.set_title('Forward data')
|
||||
@@ -1,179 +0,0 @@
|
||||
import os
|
||||
|
||||
home_dir = '.\Test_tria_2_grid'
|
||||
|
||||
inpfile = 'PYMAG3C_fwr.inp'
|
||||
|
||||
dsep = '\\'
|
||||
|
||||
os.chdir(home_dir)
|
||||
|
||||
#%%
|
||||
from SimPEG import np, sp, Utils, Mesh
|
||||
import scipy.interpolate as interpolation
|
||||
import simpegPF as PF
|
||||
import pylab as plt
|
||||
import time as tm
|
||||
|
||||
## New scripts to be added to basecode
|
||||
#from fwr_MAG_data import fwr_MAG_data
|
||||
#from read_MAGfwr_inp import read_MAGfwr_inp
|
||||
|
||||
#%%
|
||||
chibkg = 0.001
|
||||
|
||||
# Read in topo surface
|
||||
tsfile = 'Topo_Gaussian.ts'
|
||||
|
||||
# For now just read both
|
||||
topofile = 'Gaussian.topo'
|
||||
|
||||
topo = np.genfromtxt(topofile,skip_header=1)
|
||||
|
||||
# Offset data above topo
|
||||
zoffset = 2
|
||||
|
||||
# Load mesh file
|
||||
B = np.array(([90.,0.,50000.]))
|
||||
|
||||
M = np.array(([90.,0.,315.]))
|
||||
|
||||
# Sphere radius
|
||||
R = 25.
|
||||
|
||||
|
||||
#%% Script starts here
|
||||
# # Create a grid of observations and offset the z from topo
|
||||
xr = np.linspace(-99., 99., 40)
|
||||
yr = np.linspace(-49., 49., 20)
|
||||
X, Y = np.meshgrid(xr, yr)
|
||||
|
||||
F = interpolation.NearestNDInterpolator(topo[:,0:2],topo[:,2])
|
||||
Z = F(X,Y) + zoffset
|
||||
rxLoc = np.c_[Utils.mkvc(X.T), Utils.mkvc(Y.T), Utils.mkvc(Z.T)]
|
||||
|
||||
ndata = rxLoc.shape[0]
|
||||
|
||||
sclx = 100.
|
||||
|
||||
dx = 10
|
||||
|
||||
nc = int(sclx/dx)
|
||||
|
||||
|
||||
hxind = np.ones(2*nc)*dx
|
||||
hyind = np.ones(nc)*dx
|
||||
hzind = np.ones(nc)*dx
|
||||
|
||||
mesh = Mesh.TensorMesh([hxind, hyind, hzind], 'CCN')
|
||||
|
||||
# Load GOCAD surf
|
||||
#[vrtx, trgl] = PF.BaseMag.read_GOCAD_ts(tsfile)
|
||||
# Find active cells from surface
|
||||
tin = tm.time()
|
||||
print "Computing indices with VTK: "
|
||||
[indx, bc] = PF.BaseMag.gocad2vtk(tsfile,mesh, bcflag = False, inflag = True)
|
||||
print "VTK operation completed in " + str(tm.time() - tin)
|
||||
|
||||
actv = np.zeros(mesh.nC)
|
||||
actv[indx] = 1
|
||||
|
||||
model= np.zeros(mesh.nC)
|
||||
model[indx]= chibkg
|
||||
|
||||
Utils.meshutils.writeUBCTensorModel('VTKout.dat',mesh,model)
|
||||
|
||||
Utils.meshutils.writeUBCTensorMesh('Mesh_temp.msh',mesh)
|
||||
|
||||
start_time = tm.time()
|
||||
|
||||
d = PF.Magnetics.Intgrl_Fwr_Data(mesh,B,M,rxLoc,model,actv,'tmi')
|
||||
|
||||
timer = (tm.time() - start_time)
|
||||
|
||||
#%% Plot data
|
||||
plt.figure()
|
||||
ax = plt.subplot()
|
||||
plt.imshow(np.reshape(d,X.shape), interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()], origin = 'lower')
|
||||
plt.clim(0,25)
|
||||
plt.colorbar(fraction=0.02)
|
||||
plt.contour(X,Y, np.reshape(d,X.shape),10)
|
||||
plt.scatter(X,Y, c=np.reshape(d,X.shape), s=20)
|
||||
|
||||
ax.set_title('Forward data')
|
||||
|
||||
#%% First test, just brake the intersecting cells
|
||||
|
||||
ii = 1
|
||||
ddata = 9999.
|
||||
|
||||
while ddata > 2.:
|
||||
|
||||
ii += 1
|
||||
|
||||
indx = np.unravel_index(bc,(mesh.nCx,mesh.nCy,mesh.nCz), order = 'F')
|
||||
|
||||
xbreak = np.unique(indx[0])
|
||||
ybreak = np.unique(indx[1])
|
||||
zbreak = np.unique(indx[2])
|
||||
|
||||
|
||||
# Compute the new distance
|
||||
dl = np.sum(hxind[xbreak])
|
||||
dx = float(hxind[xbreak].min()/2)
|
||||
nx = dl/dx
|
||||
|
||||
hxind = np.r_[hxind[0:xbreak.min()],np.ones(nx)*dx,hxind[xbreak.max():]]
|
||||
|
||||
dl = np.sum(hyind[ybreak])
|
||||
dy = float(hyind[ybreak].min()/2)
|
||||
ny = dl/dy
|
||||
|
||||
hyind = np.r_[hyind[0:ybreak.min()],np.ones(ny)*dy,hyind[ybreak.max():]]
|
||||
|
||||
dl = np.sum(hzind[zbreak])
|
||||
dz = float(hzind[zbreak].min()/2)
|
||||
nz = dl/dz
|
||||
|
||||
hzind = np.r_[hzind[0:zbreak.min()],np.ones(nz)*dz,hzind[zbreak.max():]]
|
||||
|
||||
mesh = Mesh.TensorMesh([hxind, hyind, hzind], 'CCN')
|
||||
|
||||
# Load GOCAD surf
|
||||
#[vrtx, trgl] = PF.BaseMag.read_GOCAD_ts(tsfile)
|
||||
# Find active cells from surface
|
||||
tin = tm.time()
|
||||
print "Computing indices with VTK: "
|
||||
[indx, bc] = PF.BaseMag.gocad2vtk(tsfile,mesh, bcflag = False, inflag = True)
|
||||
print "VTK operation completed in " + str(tm.time() - tin)
|
||||
|
||||
actv = np.zeros(mesh.nC)
|
||||
actv[indx] = 1
|
||||
|
||||
model= np.zeros(mesh.nC)
|
||||
model[indx]= chibkg
|
||||
|
||||
Utils.meshutils.writeUBCTensorModel('VTKout' + str(ii) + '.dat',mesh,model)
|
||||
|
||||
Utils.meshutils.writeUBCTensorMesh('Mesh_temp' + str(ii) + '.msh',mesh)
|
||||
|
||||
start_time = tm.time()
|
||||
|
||||
d_temp = PF.Magnetics.Intgrl_Fwr_Data(mesh,B,M,rxLoc,model,actv,'tmi')
|
||||
|
||||
timer = (tm.time() - start_time)
|
||||
|
||||
ddata = np.max(abs(d-d_temp))
|
||||
|
||||
d = d_temp
|
||||
|
||||
#%% Plot data
|
||||
plt.figure()
|
||||
ax = plt.subplot()
|
||||
plt.imshow(np.reshape(d,X.shape), interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()], origin = 'lower')
|
||||
plt.clim(0,25)
|
||||
plt.colorbar(fraction=0.02)
|
||||
plt.contour(X,Y, np.reshape(d,X.shape),10)
|
||||
plt.scatter(X,Y, c=np.reshape(d,X.shape), s=20)
|
||||
|
||||
ax.set_title('Forward data')
|
||||
@@ -182,10 +182,10 @@ plt.gca().set_aspect('equal', adjustable='box')
|
||||
#%% Run one more round for sparsity
|
||||
phim = invProb.phi_m_last
|
||||
|
||||
reg = Regularization.SparseRegularization(mesh, mapping=wrMap)
|
||||
reg = Regularization.SparseRegularization(mesh, mapping=wrMap, eps=1e-4)
|
||||
reg.m = mrec
|
||||
reg.mref = mref
|
||||
reg.eps = 1e-4
|
||||
|
||||
|
||||
|
||||
diagA = np.sum(F**2.,axis=0) + beta_in*(reg.W.T*reg.W).diagonal()*(wr**2.0)
|
||||
@@ -195,7 +195,7 @@ PC = Utils.sdiag(diagA**-1.)
|
||||
|
||||
dmis = DataMisfit.l2_DataMisfit(survey)
|
||||
dmis.Wd = wd
|
||||
opt = Optimization.ProjectedGNCG(maxIter=50 ,lower=0.,upper=1.)
|
||||
opt = Optimization.ProjectedGNCG(maxIter=8 ,maxIterLS=10, maxIterCG = 20,tolCG = 1e-4,lower=0.,upper=1.)
|
||||
opt.approxHinv = PC
|
||||
#opt.phim_last = reg.eval(mrec)
|
||||
|
||||
@@ -204,7 +204,7 @@ invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = invProb.beta)
|
||||
beta = Directives.BetaSchedule(coolingFactor=1, coolingRate=1)
|
||||
#betaest = Directives.BetaEstimate_ByEig()
|
||||
target = Directives.TargetMisfit()
|
||||
IRLS =Directives.update_IRLS(factor = 2, eps_min=1e-4, phi_m_last = phim )
|
||||
IRLS =Directives.update_IRLS( eps_min=1e-3, phi_m_last = phim )
|
||||
|
||||
inv = Inversion.BaseInversion(invProb, directiveList=[beta,IRLS])
|
||||
|
||||
|
||||
+7056
-7056
File diff suppressed because it is too large
Load Diff
@@ -1,56 +0,0 @@
|
||||
GOCAD TSurf 1
|
||||
HEADER {
|
||||
name:Crown
|
||||
*solid*color:1 0.447059 0.337255 1
|
||||
ivolmap:false
|
||||
imap:false
|
||||
painted:on
|
||||
*painted*variable:normal
|
||||
last_selected_folder:3D property
|
||||
vectors3d:true
|
||||
*vectors3d*arrow:true
|
||||
*vectors3d*color:0.0980392 0.0980392 0.439216 1
|
||||
*vectors3d*arrow_size:5
|
||||
}
|
||||
GOCAD_ORIGINAL_COORDINATE_SYSTEM
|
||||
NAME Default
|
||||
AXIS_NAME "X" "Y" "Z"
|
||||
AXIS_UNIT "m" "m" "m"
|
||||
ZPOSITIVE Elevation
|
||||
END_ORIGINAL_COORDINATE_SYSTEM
|
||||
PROPERTY_CLASS_HEADER Z {
|
||||
is_z:on
|
||||
}
|
||||
PROPERTY_CLASS_HEADER vector3d {
|
||||
pclip:99
|
||||
low_clip:-0.9956403
|
||||
high_clip:0.9905306
|
||||
}
|
||||
TFACE
|
||||
VRTX 1 17.898199081420898 -24.1842041015625 -70.548362731933594
|
||||
VRTX 2 -24.411655426025391 -24.1842041015625 -69.894088745117188
|
||||
VRTX 3 17.898199081420898 25.8157958984375 -70.548362731933594
|
||||
VRTX 4 -24.411655426025391 25.8157958984375 -69.894088745117188
|
||||
VRTX 5 16.807741165161133 25.8157958984375 -9.4825878143310547
|
||||
VRTX 6 -4.565277099609375 -24.1842041015625 -28.674692153930664
|
||||
VRTX 7 16.807741165161133 -24.1842041015625 -9.4825878143310547
|
||||
VRTX 8 -24.629745483398438 -24.1842041015625 -12.972057342529297
|
||||
VRTX 9 -24.629745483398438 25.8157958984375 -12.972057342529297
|
||||
VRTX 10 -4.565277099609375 25.8157958984375 -28.674692153930664
|
||||
TRGL 8 6 9
|
||||
TRGL 7 1 5
|
||||
TRGL 1 2 3
|
||||
TRGL 9 6 10
|
||||
TRGL 10 7 5
|
||||
TRGL 2 8 4
|
||||
TRGL 5 1 3
|
||||
TRGL 3 2 4
|
||||
TRGL 10 5 3
|
||||
TRGL 4 9 10
|
||||
TRGL 7 6 1
|
||||
TRGL 6 2 1
|
||||
TRGL 4 10 3
|
||||
TRGL 4 8 9
|
||||
TRGL 8 2 6
|
||||
TRGL 6 7 10
|
||||
END
|
||||
File diff suppressed because it is too large
Load Diff
@@ -1,5 +0,0 @@
|
||||
20 20 20
|
||||
0 0 0
|
||||
0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50
|
||||
0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50
|
||||
0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50
|
||||
@@ -1,5 +0,0 @@
|
||||
80 40 40
|
||||
-100.00 -50.00 0.00
|
||||
2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50
|
||||
2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50
|
||||
2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50
|
||||
@@ -1,5 +0,0 @@
|
||||
20 10 10
|
||||
-100.00 -50.00 0.00
|
||||
10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00
|
||||
10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00
|
||||
10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00
|
||||
@@ -1,5 +0,0 @@
|
||||
41 21 15
|
||||
-105.00 -55.00 0.00
|
||||
5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 10.00
|
||||
5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 10.00
|
||||
10.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 5.00 10.00 10.00 10.00 10.00 10.00 10.00
|
||||
@@ -1,5 +0,0 @@
|
||||
85 45 26
|
||||
-110.00 -60.00 0.00
|
||||
2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 10.00
|
||||
2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 10.00
|
||||
10.00 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 2.50 5.00 10.00 10.00 10.00 10.00 10.00 10.00
|
||||
@@ -1,5 +0,0 @@
|
||||
20 10 45
|
||||
-100.00 -50.00 0.00
|
||||
10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00
|
||||
10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00
|
||||
10.00 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 1.25 2.50 2.50 2.50 2.50 5.00 5.00 10.00 10.00 10.00 10.00 10.00 10.00
|
||||
@@ -1,157 +0,0 @@
|
||||
GOCAD TSurf 1
|
||||
HEADER {
|
||||
name:SphereA_test
|
||||
*solid*color:#00ffff
|
||||
ivolmap:false
|
||||
imap:false
|
||||
painted:on
|
||||
*painted*variable:normal
|
||||
mesh:true
|
||||
}
|
||||
GOCAD_ORIGINAL_COORDINATE_SYSTEM
|
||||
NAME Default
|
||||
AXIS_NAME "X" "Y" "Z"
|
||||
AXIS_UNIT "m" "m" "m"
|
||||
ZPOSITIVE Elevation
|
||||
END_ORIGINAL_COORDINATE_SYSTEM
|
||||
PROPERTY_CLASS_HEADER Z {
|
||||
is_z:on
|
||||
}
|
||||
PROPERTY_CLASS_HEADER vector3d {
|
||||
low_clip:-1
|
||||
high_clip:1
|
||||
pclip:99
|
||||
}
|
||||
TFACE
|
||||
VRTX 1 -2.2919261455535889 -19.314384460449219 -17.626590728759766
|
||||
VRTX 2 6.6298894882202148 -23.695545196533203 -28.910013198852539
|
||||
VRTX 3 -15.073513984680176 -19.314384460449219 -38.307632446289063
|
||||
VRTX 4 -12.506924629211426 -11.605335235595703 -51.606067657470703
|
||||
VRTX 5 -19.863996505737305 -8.0631170272827148 -46.194423675537109
|
||||
VRTX 6 -22.580167770385742 3.788722038269043 -23.29432487487793
|
||||
VRTX 7 -11.330929756164551 -18.19896125793457 -20.47222900390625
|
||||
VRTX 8 6.5996670722961426 -16.974418640136719 -50.459728240966797
|
||||
VRTX 9 -17.014554977416992 -8.2651615142822266 -16.98731803894043
|
||||
VRTX 10 21.977100372314453 -7.0079841613769531 -23.695110321044922
|
||||
VRTX 11 17.244131088256836 0 -15.232500076293945
|
||||
VRTX 12 13.95079517364502 10.135844230651855 -51.434154510498047
|
||||
VRTX 13 22.229457855224609 -8.2651615142822266 -41.241447448730469
|
||||
VRTX 14 -5.8128781318664551 14.548799514770508 -13.851467132568359
|
||||
VRTX 15 9.7651596069335937 -2.9769551753997803 -56.153919219970703
|
||||
VRTX 16 -13.742574691772461 5.0091938972473145 -53.607688903808594
|
||||
VRTX 17 -13.711102485656738 20.304235458374023 -28.35902214050293
|
||||
VRTX 18 -7.821871280670166 -18.734676361083984 -47.922031402587891
|
||||
VRTX 19 13.062463760375977 21.316003799438477 -33.333328247070313
|
||||
VRTX 20 -14.062309265136719 4.0251598358154297 -13.05897045135498
|
||||
VRTX 21 24.778470993041992 2.9769551753997803 -31.861873626708984
|
||||
VRTX 22 20.487052917480469 13.627725601196289 -37.756645202636719
|
||||
VRTX 23 -6.8588962554931641 -4.9832801818847656 -9.8147735595703125
|
||||
VRTX 24 9.4985036849975586 16.97282600402832 -49.040065765380859
|
||||
VRTX 25 19.077316284179688 3.788722038269043 -49.040065765380859
|
||||
VRTX 26 -24.44097900390625 -4.9556803703308105 -31.578392028808594
|
||||
VRTX 27 15.63094425201416 11.356545448303223 -17.468194961547852
|
||||
VRTX 28 2.7494394779205322 22.420015335083008 -22.619619369506836
|
||||
VRTX 29 5.0962982177734375 -4.9596676826477051 -9.3660707473754883
|
||||
VRTX 30 -2.6028509140014648 18.735876083374023 -49.679340362548828
|
||||
VRTX 31 7.3593044281005859 6.6037797927856445 -10.371823310852051
|
||||
VRTX 32 -21.796012878417969 -12.156011581420898 -31.861873626708984
|
||||
VRTX 33 -5.0962982177734375 -4.9596676826477051 -57.300586700439453
|
||||
VRTX 34 18.269777297973633 -16.974418640136719 -31.577091217041016
|
||||
VRTX 35 -19.291585922241211 15.273146629333496 -28.910013198852539
|
||||
VRTX 36 -23.036630630493164 3.6853373050689697 -42.318424224853516
|
||||
VRTX 37 8.4685440063476562 9.8534049987792969 -11.974625587463379
|
||||
VRTX 38 16.887338638305664 -14.301387786865234 -44.964668273925781
|
||||
VRTX 39 10.637019157409668 -16.68756103515625 -18.056575775146484
|
||||
VRTX 40 -4.8257155418395996 24.485654830932617 -34.804782867431641
|
||||
VRTX 41 -21.789302825927734 12.130206108093262 -35.089565277099609
|
||||
VRTX 42 -9.5504827499389648 18.734676361083984 -46.853691101074219
|
||||
VRTX 43 -2.6198656558990479 8.0631160736083984 -56.851882934570312
|
||||
VRTX 44 -0.90874701738357544 -24.94017219543457 -34.8037109375
|
||||
VRTX 45 -15.073513984680176 19.314384460449219 -38.307632446289063
|
||||
TRGL 31 29 11
|
||||
TRGL 7 3 44
|
||||
TRGL 36 41 45
|
||||
TRGL 35 41 6
|
||||
TRGL 22 24 19
|
||||
TRGL 45 41 35
|
||||
TRGL 14 23 31
|
||||
TRGL 30 43 42
|
||||
TRGL 39 2 34
|
||||
TRGL 17 35 14
|
||||
TRGL 20 9 23
|
||||
TRGL 2 38 34
|
||||
TRGL 26 5 32
|
||||
TRGL 43 30 24
|
||||
TRGL 11 10 21
|
||||
TRGL 21 13 25
|
||||
TRGL 19 28 27
|
||||
TRGL 36 16 5
|
||||
TRGL 19 40 28
|
||||
TRGL 12 43 24
|
||||
TRGL 1 44 2
|
||||
TRGL 2 44 8
|
||||
TRGL 9 32 7
|
||||
TRGL 22 12 24
|
||||
TRGL 6 26 9
|
||||
TRGL 11 29 39
|
||||
TRGL 33 18 4
|
||||
TRGL 6 9 20
|
||||
TRGL 10 34 13
|
||||
TRGL 33 16 43
|
||||
TRGL 34 38 13
|
||||
TRGL 10 11 39
|
||||
TRGL 35 6 20
|
||||
TRGL 38 15 25
|
||||
TRGL 42 36 45
|
||||
TRGL 27 28 37
|
||||
TRGL 5 4 3
|
||||
TRGL 21 25 22
|
||||
TRGL 16 36 42
|
||||
TRGL 30 42 40
|
||||
TRGL 25 15 12
|
||||
TRGL 33 43 15
|
||||
TRGL 4 18 3
|
||||
TRGL 26 32 9
|
||||
TRGL 13 38 25
|
||||
TRGL 40 42 45
|
||||
TRGL 19 27 22
|
||||
TRGL 14 31 37
|
||||
TRGL 40 17 28
|
||||
TRGL 41 26 6
|
||||
TRGL 43 12 15
|
||||
TRGL 22 25 12
|
||||
TRGL 14 20 23
|
||||
TRGL 8 18 33
|
||||
TRGL 10 13 21
|
||||
TRGL 11 21 27
|
||||
TRGL 40 45 17
|
||||
TRGL 40 19 24
|
||||
TRGL 30 40 24
|
||||
TRGL 1 7 44
|
||||
TRGL 38 8 15
|
||||
TRGL 5 3 32
|
||||
TRGL 2 8 38
|
||||
TRGL 16 4 5
|
||||
TRGL 14 35 20
|
||||
TRGL 17 14 28
|
||||
TRGL 27 37 11
|
||||
TRGL 37 31 11
|
||||
TRGL 33 4 16
|
||||
TRGL 41 36 26
|
||||
TRGL 29 1 39
|
||||
TRGL 8 44 18
|
||||
TRGL 3 18 44
|
||||
TRGL 15 8 33
|
||||
TRGL 22 27 21
|
||||
TRGL 1 2 39
|
||||
TRGL 9 7 23
|
||||
TRGL 45 35 17
|
||||
TRGL 34 10 39
|
||||
TRGL 7 1 23
|
||||
TRGL 36 5 26
|
||||
TRGL 28 14 37
|
||||
TRGL 23 1 29
|
||||
TRGL 32 3 7
|
||||
TRGL 43 16 42
|
||||
TRGL 31 23 29
|
||||
END
|
||||
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
@@ -1,48 +0,0 @@
|
||||
-24.629745 -24.184204 -12.972057
|
||||
-4.565277 -24.184204 -28.674692
|
||||
-24.629745 25.815796 -12.972057
|
||||
16.807741 -24.184204 -9.482588
|
||||
17.898199 -24.184204 -70.548363
|
||||
16.807741 25.815796 -9.482588
|
||||
17.898199 -24.184204 -70.548363
|
||||
-24.411655 -24.184204 -69.894089
|
||||
17.898199 25.815796 -70.548363
|
||||
-24.629745 25.815796 -12.972057
|
||||
-4.565277 -24.184204 -28.674692
|
||||
-4.565277 25.815796 -28.674692
|
||||
-4.565277 25.815796 -28.674692
|
||||
16.807741 -24.184204 -9.482588
|
||||
16.807741 25.815796 -9.482588
|
||||
-24.411655 -24.184204 -69.894089
|
||||
-24.629745 -24.184204 -12.972057
|
||||
-24.411655 25.815796 -69.894089
|
||||
16.807741 25.815796 -9.482588
|
||||
17.898199 -24.184204 -70.548363
|
||||
17.898199 25.815796 -70.548363
|
||||
17.898199 25.815796 -70.548363
|
||||
-24.411655 -24.184204 -69.894089
|
||||
-24.411655 25.815796 -69.894089
|
||||
-4.565277 25.815796 -28.674692
|
||||
16.807741 25.815796 -9.482588
|
||||
17.898199 25.815796 -70.548363
|
||||
-24.411655 25.815796 -69.894089
|
||||
-24.629745 25.815796 -12.972057
|
||||
-4.565277 25.815796 -28.674692
|
||||
16.807741 -24.184204 -9.482588
|
||||
-4.565277 -24.184204 -28.674692
|
||||
17.898199 -24.184204 -70.548363
|
||||
-4.565277 -24.184204 -28.674692
|
||||
-24.411655 -24.184204 -69.894089
|
||||
17.898199 -24.184204 -70.548363
|
||||
-24.411655 25.815796 -69.894089
|
||||
-4.565277 25.815796 -28.674692
|
||||
17.898199 25.815796 -70.548363
|
||||
-24.411655 25.815796 -69.894089
|
||||
-24.629745 -24.184204 -12.972057
|
||||
-24.629745 25.815796 -12.972057
|
||||
-24.629745 -24.184204 -12.972057
|
||||
-24.411655 -24.184204 -69.894089
|
||||
-4.565277 -24.184204 -28.674692
|
||||
-4.565277 -24.184204 -28.674692
|
||||
16.807741 -24.184204 -9.482588
|
||||
-4.565277 25.815796 -28.674692
|
||||
File diff suppressed because it is too large
Load Diff
@@ -1,85 +0,0 @@
|
||||
def interpFFT(x,y,m):
|
||||
"""
|
||||
Load in a 2D grid and resample
|
||||
|
||||
OUTPUT:
|
||||
m_out
|
||||
|
||||
|
||||
"""
|
||||
|
||||
from SimPEG import np, sp
|
||||
import scipy.signal as sn
|
||||
|
||||
|
||||
# Add padding values by reflection (2**n)
|
||||
lenx = np.round( np.log2( 2*len(x) ) )
|
||||
npadx = int(np.floor( ( 2**lenx - len(x) ) /2. ))
|
||||
|
||||
#Create hemming taper
|
||||
if np.mod(npadx*2+len(x),2) != 0:
|
||||
oddx = 1
|
||||
|
||||
else:
|
||||
oddx = 0
|
||||
|
||||
tap0 = sn.hamming(npadx*2)
|
||||
tapl = sp.spdiags(tap0[0:npadx],0,npadx,npadx)
|
||||
tapr = sp.spdiags(tap0[npadx:],0,npadx,npadx+oddx)
|
||||
|
||||
# Mirror the 2d data over the half lenght and apply 0-taper
|
||||
mpad = np.hstack([np.fliplr(m[:,0:npadx]) * tapl, m, np.fliplr(m[:,-npadx:]) * tapr])
|
||||
|
||||
# Repeat along the second dimension
|
||||
leny = np.round( np.log2( 2*len(y) ) )
|
||||
npady = int(np.floor( ( 2**leny - len(y) ) /2. ))
|
||||
|
||||
#Create hemming taper
|
||||
if np.mod(npady*2+len(y),2) != 0:
|
||||
oddy = 1
|
||||
|
||||
else:
|
||||
oddy = 0
|
||||
|
||||
tap0 = sn.hamming(npady*2)
|
||||
tapu = sp.spdiags(tap0[0:npady],0,npady,npady)
|
||||
tapd = sp.spdiags(tap0[npady:],0,npady+oddy,npady)
|
||||
|
||||
mpad = np.vstack([tapu*np.flipud(mpad[0:npady,:]), mpad, tapd*np.flipud(mpad[-npady:,:])])
|
||||
|
||||
# Compute FFT
|
||||
FFTm = np.fft.fft2(mpad)
|
||||
|
||||
# Do an FFT shift
|
||||
FFTshift = np.fft.fftshift(FFTm)
|
||||
|
||||
# Pad high frequencies with zeros to increase the sampling rate
|
||||
py = int(FFTm.shape[0]/2)
|
||||
px = int(FFTm.shape[1]/2)
|
||||
|
||||
FFTshift = np.hstack([np.zeros((FFTshift.shape[0],px)),FFTshift,np.zeros((FFTshift.shape[0],px))])
|
||||
FFTshift = np.vstack([np.zeros((py,FFTshift.shape[1])),FFTshift,np.zeros((py,FFTshift.shape[1]))])
|
||||
|
||||
# Inverse shift
|
||||
FFTm = np.fft.ifftshift(FFTshift)
|
||||
|
||||
|
||||
# Compute inverse FFT
|
||||
IFFTm = np.fft.ifft2(FFTm)*FFTm.size/mpad.size
|
||||
|
||||
|
||||
m_out = np.real(IFFTm)
|
||||
# Extract core
|
||||
#m_out = np.real(IFFTm[npady*2:-(npady*2+oddy+1),npadx*2:-(npadx*2+oddx+1)])
|
||||
|
||||
|
||||
|
||||
m_out = m_out[npady*2:-(npady+oddy)*2,npadx*2:-(npadx+oddx)*2]
|
||||
|
||||
if np.mod(m.shape[0],2) != 0:
|
||||
m_out = m_out[:-1,:]
|
||||
|
||||
if np.mod(m.shape[1],2) != 0:
|
||||
m_out = m_out[:,:-1]
|
||||
|
||||
return m_out
|
||||
+10
-7
@@ -6,19 +6,22 @@ from MagAnalytics import spheremodel, CongruousMagBC
|
||||
class MagneticIntegral(Problem.BaseProblem):
|
||||
|
||||
surveyPair = Survey.LinearSurvey
|
||||
|
||||
def __init__(self, mesh, G, **kwargs):
|
||||
Problem.BaseProblem.__init__(self, mesh, **kwargs)
|
||||
|
||||
def __init__(self, mesh, G, mapping=None, **kwargs):
|
||||
Problem.BaseProblem.__init__(self, mesh, mapping=mapping, **kwargs)
|
||||
self.G = G
|
||||
|
||||
def fields(self, m):
|
||||
return self.G.dot(m)
|
||||
|
||||
|
||||
return self.G.dot(self.mapping*(m))
|
||||
|
||||
def Jvec(self, m, v, u=None):
|
||||
return self.G.dot(v)
|
||||
dmudm = self.mapping.deriv(m)
|
||||
return self.G.dot(dmudm*v)
|
||||
|
||||
def Jtvec(self, m, v, u=None):
|
||||
return self.G.T.dot(v)
|
||||
dmudm = self.mapping.deriv(m)
|
||||
return dmudm.T * (self.G.T.dot(v))
|
||||
|
||||
|
||||
|
||||
|
||||
Reference in New Issue
Block a user