In [1]:
# Setting up a model and a mesh for the MT forward problem
In [2]:
import SimPEG as simpeg, sys
import numpy as np
from SimPEG import NSEM
sys.path.append('/Volumes/MacintoshHD/Users/thibautastic/Dropbox/PhD_UBC/telluricpy')
import telluricpy
import matplotlib.pyplot as plt
%matplotlib inline
import copy
In [3]:
# Define the area of interest
bw, be = 557100, 557580
bs, bn = 7133340, 7133960
bb, bt = 0,480
In [4]:
hSizeI,vSizeI = 25., 10.
nrCcoreI = [12, 6, 6, 6, 5, 4, 3, 2, 1]
hPadI = simpeg.Utils.meshTensor([(hSizeI,9,1.5)])
hxI = np.concatenate((hPadI[::-1],np.ones(((be-bw)/hSizeI,))*hSizeI,hPadI))
hyI = np.concatenate((hPadI[::-1],np.ones(((bn-bs)/hSizeI,))*hSizeI,hPadI))
airPadI = simpeg.Utils.meshTensor([(vSizeI,12,1.5)])
vCoreI = np.concatenate([ np.ones(i)*s for i, s in zip(nrCcoreI,(simpeg.Utils.meshTensor([(vSizeI,1),(vSizeI,8,1.3)])))])[::-1]
botPadI = simpeg.Utils.meshTensor([(vCoreI[0],7,-1.5)])
hzI = np.concatenate((botPadI,vCoreI,airPadI))
# Calculate the x0 point
x0I = np.array([bw-np.sum(hPadI),bs-np.sum(hPadI),bt-np.sum(vCoreI)-np.sum(botPadI)])
# Make the mesh
meshInv = simpeg.Mesh.TensorMesh([hxI,hyI,hzI],x0I)
meshFor=copy.deepcopy(meshInv)
In [9]:
print 'the number of cells:',meshFor.nC
print meshFor
In [14]:
# Save the mesh and reload it in VTK
#meshFor.writeVTK('nsmesh_FineHKPK1.vtr',{'id':np.arange(meshFor.nC)})
nsvtr = telluricpy.vtkTools.io.readVTRFile('nsmesh_HPVK1_inv.vtr')
In [15]:
#Extract Topographie
topoSurf = telluricpy.vtkTools.polydata.normFilter(telluricpy.vtkTools.io.readVTPFile('../Geological_model/CDED_Lake_Coarse.vtp'))
activeMod = telluricpy.vtkTools.extraction.extractDataSetWithPolygon(nsvtr,topoSurf)
In [16]:
#telluricpy.vtkTools.io.writeVTUFile('activeModel.vtu',activeMod)
In [17]:
# Get active indieces
activeInd = telluricpy.vtkTools.dataset.getDataArray(activeMod,'id')
In [19]:
# Make the conductivity dictionary
# Note: using the background value for the till, since the extraction gets the ind's below the till surface
geoStructFileDict = {'Till':1e-4,
'PK1':5e-2,
'HK1':1e-3,
'VK':5e-3}
In [20]:
# Loop through
extP = '../Geological_model/'
geoStructIndDict = {}
for key, val in geoStructFileDict.iteritems():
geoPoly = telluricpy.vtkTools.polydata.normFilter(telluricpy.vtkTools.io.readVTPFile(extP+key+'.vtp'))
modStruct = telluricpy.vtkTools.extraction.extractDataSetWithPolygon(activeMod,geoPoly,extBoundaryCells=True,extInside=True,extractBounds=True)
geoStructIndDict[key] = telluricpy.vtkTools.dataset.getDataArray(modStruct,'id')
In [21]:
# Make the physical prop
sigma = np.ones(meshFor.nC)*1e-8
sigma[activeInd] = 1e-3 # 1e-4 is the background and 1e-3 is the till value
# Add the structure
for key in ['Till','PK1','HK1','VK']:
sigma[geoStructIndDict[key]] = geoStructFileDict[key]
In [22]:
#Visualisation
fig = plt.figure(figsize=(10,8))
ax = plt.subplot(111)
model = sigma.reshape(meshFor.vnC,order='F')
a = ax.pcolormesh(meshFor.gridCC[:,0].reshape(meshFor.vnC,order='F')[:,20,:],meshFor.gridCC[:,2].reshape(meshFor.vnC,order='F')[:,20,:],np.log10(model[:,20,:]),edgecolor='k',cmap='viridis')
ax.set_xlim([bw,be])
ax.set_ylim([0,bt])
#ax.grid(which="major")
cb = plt.colorbar(a)
ax.set_aspect("equal")
cb.set_label("Log conductivity (S/m)",fontsize=20)
ax.set_xlabel("Easting (m)",fontsize=20)
ax.set_ylabel("Elevation (m)",fontsize=20)
plt.tight_layout()
In [23]:
# Save the model
#meshFor.writeVTK('nsmesh_FineHKPK1Fin.vtr',{'S/m':sigma})
In [39]:
# Set up the forward modeling
freq = np.logspace(5,2,16)
np.save('MTfrequencies',freq)
In [40]:
# Find the locations on the surface of the model.
# Get the outer shell of the model
import vtk
actModVTP = telluricpy.vtkTools.polydata.normFilter(telluricpy.vtkTools.extraction.geometryFilt(activeMod))
polyBox = vtk.vtkCubeSource()
polyBox.SetBounds(bw-5.,be+5,bs-5.,bn+5.,bb-5.,bt+5)
polyBox.Update()
# Exract the topo of the model
modTopoVTP = telluricpy.vtkTools.extraction.extractDataSetWithPolygon(actModVTP,telluricpy.vtkTools.polydata.normFilter(polyBox.GetOutput()),extractBounds=False)
telluricpy.vtkTools.io.writeVTPFile('topoSurf.vtp',actModVTP)
In [52]:
# Make the rxLocations file
x,y = np.meshgrid(np.arange(bw+10.,be,20),np.arange(bs+10.,bn,20))
xy = np.hstack((x.reshape(-1,1),y.reshape(-1,1)))
# Find the location array
locArr = telluricpy.modelTools.surfaceIntersect.findZofXYOnPolydata(xy,actModVTP) #modTopoVTP)
np.save('MTlocations',locArr)
In [53]:
telluricpy.vtkTools.io.writeVTPFile('MTloc.vtp',telluricpy.dataFiles.XYZtools.makeCylinderPtsVTP(locArr,5,10,10))
In [60]:
# Running the forward modelling on the Cluster.
# Define the forward run in findDiam_MTforward.py
In [71]:
%matplotlib qt
sys.path.append('/home/gudni/Dropbox/code/python/MTview/')
import interactivePlotFunctions as iPf
In [ ]:
In [72]:
# Load the data
mtData = np.load('MTdataStArr_nsmesh_HKPK1.npy')
mtData
Out[72]:
In [110]:
iPf.MTinteractiveMap([mtData])
Out[110]:
In [104]:
# Looking at the data shows that data below 100Hz is affected by the boundary conditions,
# which makes sense for very conductive conditions as we have.
# Invert data in the 1e5-1e2 range.
Run the inversion on the cluster using the inv3d/run1/findDiam_inversion.py
In [106]:
drecAll = np.load('MTdataStArr_nsmesh_0.npy')
In [109]:
np.unique(drecAll['freq'])[10::]
Out[109]:
In [93]:
# Build the Inversion mesh
# Design the tensors
hSizeI,vSizeI = 25., 10.
nrCcoreI = [12, 6, 6, 6, 5, 4, 3, 2, 1]
hPadI = simpeg.Utils.meshTensor([(hSizeI,9,1.5)])
hxI = np.concatenate((hPadI[::-1],np.ones(((be-bw)/hSizeI,))*hSizeI,hPadI))
hyI = np.concatenate((hPadI[::-1],np.ones(((bn-bs)/hSizeI,))*hSizeI,hPadI))
airPadI = simpeg.Utils.meshTensor([(vSizeI,12,1.5)])
vCoreI = np.concatenate([ np.ones(i)*s for i, s in zip(nrCcoreI,(simpeg.Utils.meshTensor([(vSizeI,1),(vSizeI,8,1.3)])))])[::-1]
botPadI = simpeg.Utils.meshTensor([(vCoreI[0],7,-1.5)])
hzI = np.concatenate((botPadI,vCoreI,airPadI))
# Calculate the x0 point
x0I = np.array([bw-np.sum(hPadI),bs-np.sum(hPadI),bt-np.sum(vCoreI)-np.sum(botPadI)])
# Make the mesh
meshInv = simpeg.Mesh.TensorMesh([hxI,hyI,hzI],x0I)
In [94]:
meshInv.writeVTK('nsmesh_HPVK1_inv.vtr',{'id':np.arange(meshInv.nC)})
In [101]:
nsInvvtr = telluricpy.vtkTools.io.readVTRFile('nsmesh_HPVK1_inv.vtr')
activeModInv = telluricpy.vtkTools.extraction.extractDataSetWithPolygon(nsInvvtr,topoSurf,extBoundaryCells=True)
In [102]:
sigma = np.ones(meshInv.nC)*1e-8
indAct = telluricpy.vtkTools.dataset.getDataArray(activeModInv,'id')
sigma[indAct] = 1e-4
In [103]:
meshInv.writeVTK('nsmesh_HPVK1_inv.vtr',{'id':np.arange(meshInv.nC),'S/m':sigma})
In [108]:
from pymatsolver import MumpsSolver
In [106]:
pymatsolver.AvailableSolvers
Out[106]:
In [116]:
NSEM.Utils.skindepth(1000,100000)
Out[116]:
In [117]:
np.unique(mtData['freq'])[5::]
Out[117]:
In [ ]: