mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-10 01:08:14 +08:00
On the fly forward modeling and pseudo section for 3D conductivity model
This commit is contained in:
D Fournier
@@ -10,49 +10,19 @@ from SimPEG import np, Utils, Mesh, mkvc, SolverLU, sp
|
||||
import simpegDCIP as DC
|
||||
import pylab as plt
|
||||
import time
|
||||
from scipy.interpolate import griddata
|
||||
import numpy.matlib as npm
|
||||
|
||||
#from scipy.linalg import solve_banded
|
||||
|
||||
# Load UBC mesh 3D
|
||||
#mesh = Utils.meshutils.readUBCTensorMesh('Mesh_20m.msh')
|
||||
mesh = Utils.meshutils.readUBCTensorMesh('Mesh_40m.msh')
|
||||
|
||||
# Load model
|
||||
#model = Utils.meshutils.readUBCTensorModel('MtIsa_3D.con',mesh)
|
||||
model = Utils.meshutils.readUBCTensorModel('Synthetic.con',mesh)
|
||||
|
||||
#%%
|
||||
# Display top section
|
||||
top = int(mesh.nCz)-1
|
||||
mesh.plotSlice(model, ind=top, normal='Z', grid=True, pcolorOpts={'alpha':0.8})
|
||||
ylim=(546000,546750)
|
||||
xlim=(422900,423675)
|
||||
# Takes two points from ginput and create survey
|
||||
temp = plt.ginput(2)
|
||||
|
||||
# Add z coordinate
|
||||
nz = mesh.vectorNz
|
||||
endp = np.c_[np.asarray(temp),np.ones(2).T*nz[-1]]
|
||||
|
||||
# Create dipole survey receivers and plot
|
||||
ab = 40
|
||||
a = 20
|
||||
|
||||
# Evenly distribute transmitters for now and put on surface
|
||||
dplen = np.sqrt( np.sum((endp[1,:] - endp[0,:])**2) )
|
||||
dp_x = ( endp[1,0] - endp[0,0] ) / dplen
|
||||
dp_y = ( endp[1,1] - endp[0,1] ) / dplen
|
||||
|
||||
nstn = np.floor( dplen / ab )
|
||||
nrx = nstn-1
|
||||
|
||||
stn_x = endp[0,0] + np.cumsum( np.ones(nstn)*dp_x*ab )
|
||||
stn_y = endp[0,1] + np.cumsum( np.ones(nstn)*dp_y*ab )
|
||||
|
||||
plt.scatter(stn_x,stn_y,s=100, c='w')
|
||||
|
||||
M = np.c_[stn_x-a*dp_x, stn_y-a*dp_y, np.ones(nstn).T*nz[-1]]
|
||||
N = np.c_[stn_x+a*dp_x, stn_y+a*dp_y, np.ones(nstn).T*nz[-1]]
|
||||
|
||||
plt.scatter(M[:,0],M[:,1],s=10,c='r')
|
||||
plt.scatter(N[:,0],N[:,1],s=10,c='b')
|
||||
|
||||
#%% Create system
|
||||
#Set boundary conditions
|
||||
mesh.setCellGradBC('neumann')
|
||||
@@ -65,24 +35,68 @@ A = Div*Msig*Grad
|
||||
|
||||
# Change one corner to deal with nullspace
|
||||
A[0,0] = 1.
|
||||
A = sp.csc_matrix(A)
|
||||
|
||||
start_time = time.time()
|
||||
|
||||
# Factor A matrix
|
||||
Ainv = sp.linalg.splu(A)
|
||||
|
||||
print("LU DECOMP--- %s seconds ---" % (time.time() - start_time))
|
||||
|
||||
#%%
|
||||
# Display top section
|
||||
top = int(mesh.nCz)-1
|
||||
mesh.plotSlice(model, ind=top, normal='Z', grid=True, pcolorOpts={'alpha':0.8})
|
||||
|
||||
|
||||
# Takes two points from ginput and create survey
|
||||
temp = plt.ginput(2)
|
||||
|
||||
# Add z coordinate
|
||||
nz = mesh.vectorNz
|
||||
endp = np.c_[np.asarray(temp),np.ones(2).T*nz[-1]]
|
||||
|
||||
# Create dipole survey receivers and plot
|
||||
a = 40
|
||||
n = 8
|
||||
|
||||
# Evenly distribute transmitters for now and put on surface
|
||||
dplen = np.sqrt( np.sum((endp[1,:] - endp[0,:])**2) )
|
||||
dp_x = ( endp[1,0] - endp[0,0] ) / dplen
|
||||
dp_y = ( endp[1,1] - endp[0,1] ) / dplen
|
||||
|
||||
nstn = np.floor( dplen / a )
|
||||
nrx = nstn-1
|
||||
|
||||
stn_x = endp[0,0] + np.cumsum( np.ones(nstn)*dp_x*a )
|
||||
stn_y = endp[0,1] + np.cumsum( np.ones(nstn)*dp_y*a )
|
||||
|
||||
plt.scatter(stn_x,stn_y,s=100, c='w')
|
||||
|
||||
M = np.c_[stn_x-a*dp_x/2, stn_y-a*dp_y/2, np.ones(nstn).T*nz[-1]]
|
||||
N = np.c_[stn_x+a*dp_x/2, stn_y+a*dp_y/2, np.ones(nstn).T*nz[-1]]
|
||||
|
||||
plt.scatter(M[:,0],M[:,1],s=10,c='r')
|
||||
plt.scatter(N[:,0],N[:,1],s=10,c='b')
|
||||
|
||||
|
||||
|
||||
#%% Forward model data
|
||||
data = np.zeros( nstn*nrx )
|
||||
data = []#np.zeros( nstn*nrx )
|
||||
problem = DC.ProblemDC_CC(mesh)
|
||||
fig = plt.figure()
|
||||
|
||||
|
||||
|
||||
|
||||
for ii in range(0, int(nstn)):
|
||||
|
||||
for ii in range(0, int(nstn)-2):
|
||||
start_time = time.time()
|
||||
|
||||
rxloc_M = np.r_[M[0:ii,:],M[ii+1:,:]]
|
||||
rxloc_N = np.r_[N[0:ii,:],N[ii+1:,:]]
|
||||
# Select dipole locations for receiver: n || end of line
|
||||
|
||||
idx = int( np.min([ii+n+1,nstn+1]) )
|
||||
rxloc_M = M[ii+2:ii+n+1,:]#np.r_[M[0:ii,:],M[ii+1:,:]]
|
||||
rxloc_N = N[ii+2:ii+n+1,:]#np.r_[N[0:ii,:],N[ii+1:,:]]
|
||||
|
||||
nrx = rxloc_M.shape[0]
|
||||
|
||||
Rx = DC.RxDipole(rxloc_M,rxloc_N)
|
||||
|
||||
@@ -98,25 +112,70 @@ for ii in range(0, int(nstn)):
|
||||
P1 = mesh.getInterpolationMat(rxloc_M, 'CC')
|
||||
P2 = mesh.getInterpolationMat(rxloc_N, 'CC')
|
||||
|
||||
#Direct Solve
|
||||
phi = Ainv.solve(RHS)
|
||||
d = P1*phi - P2*phi
|
||||
|
||||
data[(nrx*(ii)):nrx+(nrx*(ii))] = d.T#survey.dpred(model)
|
||||
# Iterative Solve
|
||||
#Ainvb = sp.linalg.bicgstab(A,RHS, tol=1e-5)
|
||||
#phi = mkvc(Ainvb[0])
|
||||
|
||||
# Compute potential at each electrode
|
||||
d = P1*phi - P2*phi
|
||||
|
||||
# Convert 3D location to distance along survey line for 2D
|
||||
# Plot pseudo section along line
|
||||
txmidx = endp[0,0] - np.mean(np.c_[M[ii,0],N[ii,0]])
|
||||
rxmidx = endp[0,0] - np.mean( np.c_[rxloc_M[:,0], rxloc_N[:,0]], axis=1 )
|
||||
|
||||
txmidy = endp[0,1] - np.mean(np.c_[M[ii,1],N[ii,1]])
|
||||
rxmidy = endp[0,1] - np.mean( np.c_[rxloc_M[:,1], rxloc_N[:,1]], axis=1 )
|
||||
|
||||
rxmid = np.sqrt(rxmidx**2 + rxmidy**2)
|
||||
txmid = np.sqrt(txmidx**2 + txmidy**2)
|
||||
|
||||
midp = ( rxmid + txmid )/2
|
||||
rP1 = np.sqrt( np.sum( ( endp[0,:] - M[ii,:] )**2 , axis=0))
|
||||
rP2 = np.sqrt( np.sum( ( endp[0,:] - N[ii,:] )**2 , axis=0))
|
||||
rC1 = np.sqrt( np.sum( ( npm.repmat(endp[0,:],nrx, 1) - rxloc_M )**2 , axis=1))
|
||||
rC2 = np.sqrt( np.sum( ( npm.repmat(endp[0,:],nrx, 1) - rxloc_N )**2 , axis=1))
|
||||
|
||||
if ii == 0:
|
||||
data = np.c_[np.ones(nrx)*rP1, np.ones(nrx)*rP2, rC1, rC2, mkvc(d), np.ones(nrx)*1e-2]
|
||||
|
||||
else:
|
||||
temp = np.c_[np.ones(nrx)*rP1, np.ones(nrx)*rP2, rC1, rC2, mkvc(d), np.ones(nrx)*1e-2]
|
||||
data = np.r_[data,temp]#survey.dpred(model)
|
||||
|
||||
|
||||
print("--- %s seconds ---" % (time.time() - start_time))
|
||||
|
||||
|
||||
|
||||
plt.scatter(midp,-np.abs(txmid-midp),s=50,c=data[(nrx*(ii)):nrx+(nrx*(ii))])
|
||||
# Write data to UBC-2D format
|
||||
#temp = np.c_[np.ones(nrx)*txmid-a/2, np.ones(nrx)*txmid+a/2,
|
||||
# rxmid-a/2, rxmid+a/2,
|
||||
# mkvc(d) , np.ones(nrx)*1e-2]
|
||||
|
||||
|
||||
|
||||
fid = open(home_dir + '\FWR_data.dat','w')
|
||||
fid.write('SIMPEG FORWARD\n')
|
||||
np.savetxt(fid, data, fmt='%e',delimiter=' ',newline='\n')
|
||||
fid.close()
|
||||
|
||||
#%% Plot pseudo section
|
||||
|
||||
# Get distances between each poles
|
||||
rC1P1 = data[:,0] - data[:,2]
|
||||
rC2P1 = data[:,0] - data[:,3]
|
||||
rC1P2 = data[:,1] - data[:,2]
|
||||
rC2P2 = data[:,1] - data[:,3]
|
||||
|
||||
# Compute apparent resistivity
|
||||
rho = data[:,4] * 2*np.pi / ( 1/rC1P1 - 1/rC2P1 - 1/rC1P2 + 1/rC2P2 )
|
||||
|
||||
Cmid = (data[:,0] + data[:,1])/2
|
||||
Pmid = (data[:,2] + data[:,3])/2
|
||||
|
||||
midp = ( Cmid + Pmid )/2
|
||||
midz = -np.abs(Cmid-Pmid)
|
||||
|
||||
# Grid points
|
||||
grid_x, grid_z = np.mgrid[np.min(midp):np.max(midp), np.min(midz):np.max(midz)]
|
||||
grid_rho = griddata(np.c_[midp,midz], np.log10(abs(rho.T)), (grid_x, grid_z), method='linear')
|
||||
plt.imshow(grid_rho.T, extent = (np.min(midp),np.max(midp),np.min(midz),np.max(midz)), origin='lower')
|
||||
plt.colorbar()
|
||||
|
||||
# Plot apparent resistivity
|
||||
plt.scatter(midp,midz,s=50,c=np.log10(abs(rho.T)))
|
||||
@@ -0,0 +1,78 @@
|
||||
SIMPEG FORWARD
|
||||
2.000000e+01 6.000000e+01 1.000000e+02 1.400000e+02 -1.829346e-01 1.000000e-02
|
||||
2.000000e+01 6.000000e+01 1.400000e+02 1.800000e+02 -7.168519e-02 1.000000e-02
|
||||
2.000000e+01 6.000000e+01 1.800000e+02 2.200000e+02 -3.011343e-02 1.000000e-02
|
||||
2.000000e+01 6.000000e+01 2.200000e+02 2.600000e+02 -1.433145e-02 1.000000e-02
|
||||
2.000000e+01 6.000000e+01 2.600000e+02 3.000000e+02 -9.242969e-03 1.000000e-02
|
||||
2.000000e+01 6.000000e+01 3.000000e+02 3.400000e+02 -5.889792e-03 1.000000e-02
|
||||
2.000000e+01 6.000000e+01 3.400000e+02 3.800000e+02 -2.128083e-04 1.000000e-02
|
||||
6.000000e+01 1.000000e+02 1.400000e+02 1.800000e+02 -6.222483e-01 1.000000e-02
|
||||
6.000000e+01 1.000000e+02 1.800000e+02 2.200000e+02 -1.592011e-01 1.000000e-02
|
||||
6.000000e+01 1.000000e+02 2.200000e+02 2.600000e+02 -6.560129e-02 1.000000e-02
|
||||
6.000000e+01 1.000000e+02 2.600000e+02 3.000000e+02 -3.832817e-02 1.000000e-02
|
||||
6.000000e+01 1.000000e+02 3.000000e+02 3.400000e+02 -2.357037e-02 1.000000e-02
|
||||
6.000000e+01 1.000000e+02 3.400000e+02 3.800000e+02 -8.335221e-04 1.000000e-02
|
||||
6.000000e+01 1.000000e+02 3.800000e+02 4.200000e+02 -6.171719e-04 1.000000e-02
|
||||
1.000000e+02 1.400000e+02 1.800000e+02 2.200000e+02 -1.767638e-01 1.000000e-02
|
||||
1.000000e+02 1.400000e+02 2.200000e+02 2.600000e+02 -7.728604e-02 1.000000e-02
|
||||
1.000000e+02 1.400000e+02 2.600000e+02 3.000000e+02 -3.834522e-02 1.000000e-02
|
||||
1.000000e+02 1.400000e+02 3.000000e+02 3.400000e+02 -2.333121e-02 1.000000e-02
|
||||
1.000000e+02 1.400000e+02 3.400000e+02 3.800000e+02 -8.301106e-04 1.000000e-02
|
||||
1.000000e+02 1.400000e+02 3.800000e+02 4.200000e+02 -4.325777e-04 1.000000e-02
|
||||
1.000000e+02 1.400000e+02 4.200000e+02 4.600000e+02 -5.007358e-04 1.000000e-02
|
||||
1.400000e+02 1.800000e+02 2.200000e+02 2.600000e+02 -1.604845e-01 1.000000e-02
|
||||
1.400000e+02 1.800000e+02 2.600000e+02 3.000000e+02 -6.987283e-02 1.000000e-02
|
||||
1.400000e+02 1.800000e+02 3.000000e+02 3.400000e+02 -3.263975e-02 1.000000e-02
|
||||
1.400000e+02 1.800000e+02 3.400000e+02 3.800000e+02 -9.927530e-04 1.000000e-02
|
||||
1.400000e+02 1.800000e+02 3.800000e+02 4.200000e+02 -7.313568e-04 1.000000e-02
|
||||
1.400000e+02 1.800000e+02 4.200000e+02 4.600000e+02 -8.461857e-04 1.000000e-02
|
||||
1.400000e+02 1.800000e+02 4.600000e+02 5.000000e+02 -6.191896e-04 1.000000e-02
|
||||
1.800000e+02 2.200000e+02 2.600000e+02 3.000000e+02 -5.233396e-01 1.000000e-02
|
||||
1.800000e+02 2.200000e+02 3.000000e+02 3.400000e+02 -1.694418e-01 1.000000e-02
|
||||
1.800000e+02 2.200000e+02 3.400000e+02 3.800000e+02 -4.695903e-03 1.000000e-02
|
||||
1.800000e+02 2.200000e+02 3.800000e+02 4.200000e+02 -2.408726e-03 1.000000e-02
|
||||
1.800000e+02 2.200000e+02 4.200000e+02 4.600000e+02 -2.525300e-03 1.000000e-02
|
||||
1.800000e+02 2.200000e+02 4.600000e+02 5.000000e+02 -1.797016e-03 1.000000e-02
|
||||
1.800000e+02 2.200000e+02 5.000000e+02 5.400000e+02 -1.307715e-03 1.000000e-02
|
||||
2.200000e+02 2.600000e+02 3.000000e+02 3.400000e+02 -1.029885e+00 1.000000e-02
|
||||
2.200000e+02 2.600000e+02 3.400000e+02 3.800000e+02 -1.760991e-02 1.000000e-02
|
||||
2.200000e+02 2.600000e+02 3.800000e+02 4.200000e+02 -6.266768e-03 1.000000e-02
|
||||
2.200000e+02 2.600000e+02 4.200000e+02 4.600000e+02 -5.460624e-03 1.000000e-02
|
||||
2.200000e+02 2.600000e+02 4.600000e+02 5.000000e+02 -3.411365e-03 1.000000e-02
|
||||
2.200000e+02 2.600000e+02 5.000000e+02 5.400000e+02 -2.254223e-03 1.000000e-02
|
||||
2.200000e+02 2.600000e+02 5.400000e+02 5.800000e+02 -1.571200e-03 1.000000e-02
|
||||
2.600000e+02 3.000000e+02 3.400000e+02 3.800000e+02 -2.644781e-02 1.000000e-02
|
||||
2.600000e+02 3.000000e+02 3.800000e+02 4.200000e+02 -1.198780e-02 1.000000e-02
|
||||
2.600000e+02 3.000000e+02 4.200000e+02 4.600000e+02 -1.224957e-02 1.000000e-02
|
||||
2.600000e+02 3.000000e+02 4.600000e+02 5.000000e+02 -8.806925e-03 1.000000e-02
|
||||
2.600000e+02 3.000000e+02 5.000000e+02 5.400000e+02 -6.551031e-03 1.000000e-02
|
||||
2.600000e+02 3.000000e+02 5.400000e+02 5.800000e+02 -4.896339e-03 1.000000e-02
|
||||
2.600000e+02 3.000000e+02 5.800000e+02 6.200000e+02 -3.707667e-03 1.000000e-02
|
||||
3.000000e+02 3.400000e+02 3.800000e+02 4.200000e+02 -2.027985e-02 1.000000e-02
|
||||
3.000000e+02 3.400000e+02 4.200000e+02 4.600000e+02 -1.305516e-02 1.000000e-02
|
||||
3.000000e+02 3.400000e+02 4.600000e+02 5.000000e+02 -6.743298e-03 1.000000e-02
|
||||
3.000000e+02 3.400000e+02 5.000000e+02 5.400000e+02 -3.821798e-03 1.000000e-02
|
||||
3.000000e+02 3.400000e+02 5.400000e+02 5.800000e+02 -2.494772e-03 1.000000e-02
|
||||
3.000000e+02 3.400000e+02 5.800000e+02 6.200000e+02 -1.676624e-03 1.000000e-02
|
||||
3.000000e+02 3.400000e+02 6.200000e+02 6.600000e+02 -1.132107e-03 1.000000e-02
|
||||
3.400000e+02 3.800000e+02 4.200000e+02 4.600000e+02 -1.993678e-02 1.000000e-02
|
||||
3.400000e+02 3.800000e+02 4.600000e+02 5.000000e+02 -4.958931e-03 1.000000e-02
|
||||
3.400000e+02 3.800000e+02 5.000000e+02 5.400000e+02 -1.775192e-03 1.000000e-02
|
||||
3.400000e+02 3.800000e+02 5.400000e+02 5.800000e+02 -8.165377e-04 1.000000e-02
|
||||
3.400000e+02 3.800000e+02 5.800000e+02 6.200000e+02 -4.006810e-04 1.000000e-02
|
||||
3.400000e+02 3.800000e+02 6.200000e+02 6.600000e+02 -2.235690e-04 1.000000e-02
|
||||
3.800000e+02 4.200000e+02 4.600000e+02 5.000000e+02 -2.050915e-02 1.000000e-02
|
||||
3.800000e+02 4.200000e+02 5.000000e+02 5.400000e+02 -1.214802e-02 1.000000e-02
|
||||
3.800000e+02 4.200000e+02 5.400000e+02 5.800000e+02 -8.207648e-03 1.000000e-02
|
||||
3.800000e+02 4.200000e+02 5.800000e+02 6.200000e+02 -5.784972e-03 1.000000e-02
|
||||
3.800000e+02 4.200000e+02 6.200000e+02 6.600000e+02 -4.210745e-03 1.000000e-02
|
||||
4.200000e+02 4.600000e+02 5.000000e+02 5.400000e+02 -2.770260e-01 1.000000e-02
|
||||
4.200000e+02 4.600000e+02 5.400000e+02 5.800000e+02 -1.050325e-01 1.000000e-02
|
||||
4.200000e+02 4.600000e+02 5.800000e+02 6.200000e+02 -4.296203e-02 1.000000e-02
|
||||
4.200000e+02 4.600000e+02 6.200000e+02 6.600000e+02 -2.089104e-02 1.000000e-02
|
||||
4.600000e+02 5.000000e+02 5.400000e+02 5.800000e+02 -5.583148e-01 1.000000e-02
|
||||
4.600000e+02 5.000000e+02 5.800000e+02 6.200000e+02 -1.542567e-01 1.000000e-02
|
||||
4.600000e+02 5.000000e+02 6.200000e+02 6.600000e+02 -6.312412e-02 1.000000e-02
|
||||
5.000000e+02 5.400000e+02 5.800000e+02 6.200000e+02 -1.649706e-01 1.000000e-02
|
||||
5.000000e+02 5.400000e+02 6.200000e+02 6.600000e+02 -7.170996e-02 1.000000e-02
|
||||
5.400000e+02 5.800000e+02 6.200000e+02 6.600000e+02 -1.452292e-01 1.000000e-02
|
||||
@@ -0,0 +1,5 @@
|
||||
100 21 55
|
||||
-915.00 10725 20.00
|
||||
250 200 175 150 103.00 86.00 72.00 60.00 50.00 40.00 35.00 30.00 24.00 74*20.00 24.00 30.00 35.00 40.00 50.00 60.00 72.00 86.00 103.00 150 175 200 250
|
||||
300 250 200 175 150 103.00 86.00 72.00 5*60 72.00 86.00 103.00 150 175 200 250 300
|
||||
20*10.00 12 16 20*20 24.00 30.00 35.00 40.00 50.00 60.00 72.00 86.00 103.00 150 175 200 250
|
||||
+115500
File diff suppressed because it is too large
Load Diff
Reference in New Issue
Block a user