mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-03 02:00:37 +08:00
Merge branch 'FDEMrefactor' of https://github.com/simpeg/simpegmt into FDEMrefactor
This commit is contained in:
Lindsey Heagy
@@ -42,3 +42,4 @@ nosetests.xml
|
||||
*.sublime-workspace
|
||||
docs/_build/
|
||||
*.ipynb_checkpoints
|
||||
notebooks/scipy2015/001-Inversion_NoStopping.npy
|
||||
|
||||
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
@@ -31,26 +31,30 @@
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"## Setup the modelling\n",
|
||||
"# Setting up 1D mesh and conductivity models to forward model data.\n",
|
||||
"\n",
|
||||
"# Frequency\n",
|
||||
"nFreq = 31\n",
|
||||
"freqs = np.logspace(3,-3,nFreq)\n",
|
||||
"# Set mesh parameters\n",
|
||||
"ct = 10\n",
|
||||
"air = simpeg.Utils.meshTensor([(ct,25,1.3)])\n",
|
||||
"core = np.concatenate( ( np.kron(simpeg.Utils.meshTensor([(ct,5,-1.2)]),np.ones((3,))) , simpeg.Utils.meshTensor([(ct,5)]) ) )\n",
|
||||
"bot = simpeg.Utils.meshTensor([(core[0],25,-1.3)])\n",
|
||||
"ct = 20\n",
|
||||
"air = simpeg.Utils.meshTensor([(ct,16,1.4)])\n",
|
||||
"core = np.concatenate( ( np.kron(simpeg.Utils.meshTensor([(ct,10,-1.3)]),np.ones((5,))) , simpeg.Utils.meshTensor([(ct,5)]) ) )\n",
|
||||
"bot = simpeg.Utils.meshTensor([(core[0],10,-1.4)])\n",
|
||||
"x0 = -np.array([np.sum(np.concatenate((core,bot)))])\n",
|
||||
"# Make the model\n",
|
||||
"m1d = simpeg.Mesh.TensorMesh([np.concatenate((bot,core,air))], x0=x0)\n",
|
||||
"\n",
|
||||
"# Setup model varibles\n",
|
||||
"active = m1d.vectorCCx<0.\n",
|
||||
"layer1 = (m1d.vectorCCx<-200.) & (m1d.vectorCCx>=-600.)\n",
|
||||
"layer2 = (m1d.vectorCCx<-2000.) & (m1d.vectorCCx>=-4000.)\n",
|
||||
"layer1 = (m1d.vectorCCx<-500.) & (m1d.vectorCCx>=-800.)\n",
|
||||
"layer2 = (m1d.vectorCCx<-3500.) & (m1d.vectorCCx>=-5000.)\n",
|
||||
"# Set the conductivity values\n",
|
||||
"sig_half = 2e-3\n",
|
||||
"sig_air = 1e-8\n",
|
||||
"sig_layer1 = 1\n",
|
||||
"sig_layer2 = .1\n",
|
||||
"sig_layer1 = .2\n",
|
||||
"sig_layer2 = .2\n",
|
||||
"# Make the true model\n",
|
||||
"sigma_true = np.ones(m1d.nCx)*sig_air\n",
|
||||
"sigma_true[active] = sig_half\n",
|
||||
@@ -63,8 +67,7 @@
|
||||
"sigma_0[active] = sig_half\n",
|
||||
"m_0 = np.log(sigma_0[active])\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"# Set the mapping# Set the mapping\n",
|
||||
"# Set the mapping\n",
|
||||
"actMap = simpeg.Maps.ActiveCells(m1d, active, np.log(1e-8), nC=m1d.nCx)\n",
|
||||
"mappingExpAct = simpeg.Maps.ExpMap(m1d) * actMap"
|
||||
]
|
||||
|
||||
Binary file not shown.
Binary file not shown.
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
@@ -11,6 +11,175 @@
|
||||
"## Forward model a data"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import SimPEG as simpeg\n",
|
||||
"import simpegMT as simpegmt\n",
|
||||
"import cPickle as pickle"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"## Setup the forward modeling\n",
|
||||
"# Read the model\n",
|
||||
"modelname = \"simpegTDmodel.con\"\n",
|
||||
"sigma = np.loadtxt(modelname)\n",
|
||||
"# Make the mesh.\n",
|
||||
"mTensor = simpeg.Utils.meshTensor\n",
|
||||
"cSize = [50,20]\n",
|
||||
"# Cells constant size mesh\n",
|
||||
"hx = mTensor([(cSize[0],50)])\n",
|
||||
"hy = mTensor([(cSize[0],50)])\n",
|
||||
"hz = mTensor([(cSize[1],48)])\n",
|
||||
"x0 = np.array([-1250,-1250,- 30*20])\n",
|
||||
"mesh3dCons = simpeg.Mesh.TensorMesh([hx,hy,hz],x0)\n",
|
||||
"# With padding\n",
|
||||
"hPad = mTensor([(cSize[0],5,1.5)])\n",
|
||||
"aPad = mTensor([(cSize[1],13,1.3)])\n",
|
||||
"bPad = mTensor([(cSize[1],5,-1.5)])\n",
|
||||
"hxPad = np.hstack((hPad[::-1],mTensor([(cSize[0],40)]),hPad))\n",
|
||||
"hyPad = np.hstack((hPad[::-1],mTensor([(cSize[0],40)]),hPad))\n",
|
||||
"hzPad = np.hstack((bPad,mTensor([(cSize[1],30)]),aPad))\n",
|
||||
"x0Pad = np.array([-(np.sum(hPad)+1000),-(np.sum(hPad)+1000),-(np.sum(bPad)+600)])\n",
|
||||
"mesh3d = simpeg.Mesh.TensorMesh([hxPad,hyPad,hzPad],x0Pad)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Load the model to the uniform cell mesh\n",
|
||||
"modelUniCell = simpeg.Utils.meshutils.readUBCTensorModel(modelname,mesh3dCons)\n",
|
||||
"# Save as a vtk file\n",
|
||||
"simpeg.Utils.meshutils.writeVTRFile('modelTDuniMesh.vtr',mesh3dCons,{'S/m':modelUniCell})"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 5,
|
||||
"metadata": {
|
||||
"collapsed": false,
|
||||
"scrolled": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Load the model to the mesh with padding cells\n",
|
||||
"modelTD = simpeg.Utils.meshutils.readUBCTensorModel(modelname,mesh3d)\n",
|
||||
"# Save as a vtk file\n",
|
||||
"simpeg.Utils.meshutils.writeVTRFile('modelTDpaddedMesh.vtr',mesh3d,{'S/m':modelTD})"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 6,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Define the data locations\n",
|
||||
"xG,yG = np.meshgrid(np.linspace(-700,700,8),np.linspace(-700,700,8))\n",
|
||||
"zG = np.zeros_like(xG)\n",
|
||||
"locs = np.hstack((simpeg.mkvc(xG.ravel(),2),simpeg.mkvc(yG.ravel(),2),simpeg.mkvc(zG.ravel(),2)))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 7,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"ename": "SolverException",
|
||||
"evalue": "Mumps Exception [-13] - An error occurred in a Fortran ALLOCATE statement. The size that the package requested is available in INFO(2). If INFO(2) is negative, then the size that the package requested is obtained by multiplying the absolute value of INFO(2) by 1 million.",
|
||||
"output_type": "error",
|
||||
"traceback": [
|
||||
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
|
||||
"\u001b[1;31mSolverException\u001b[0m Traceback (most recent call last)",
|
||||
"\u001b[1;32m<ipython-input-7-10aaec78a1f3>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m 20\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 21\u001b[0m \u001b[1;31m# Forward model the data\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 22\u001b[1;33m \u001b[0mfields\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mproblem\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfields\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmodelTD\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 23\u001b[0m \u001b[0mmtData\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0msurvey\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mprojectFields\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfields\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpegmt/simpegMT/ProblemMT3D/Problems.pyc\u001b[0m in \u001b[0;36mfields\u001b[1;34m(self, m)\u001b[0m\n\u001b[0;32m 100\u001b[0m \u001b[0mA\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetA\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfreq\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 101\u001b[0m \u001b[0mrhs\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetRHS\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfreq\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 102\u001b[1;33m \u001b[0mAinv\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mSolver\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mA\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msolverOpts\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 103\u001b[0m \u001b[0me_s\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mAinv\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mrhs\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 104\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/pymatsolver/pymatsolver/Mumps/__init__.pyc\u001b[0m in \u001b[0;36m__init__\u001b[1;34m(self, A, symmetric, fromPointer)\u001b[0m\n\u001b[0;32m 96\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 97\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0mfromPointer\u001b[0m \u001b[1;32mis\u001b[0m \u001b[0mNone\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 98\u001b[1;33m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfactor\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 99\u001b[0m \u001b[1;32melif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfromPointer\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0m_Pointer\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 100\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mpointer\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mfromPointer\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/pymatsolver/pymatsolver/Mumps/__init__.pyc\u001b[0m in \u001b[0;36mfactor\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m 132\u001b[0m self.A.indptr+1)\n\u001b[0;32m 133\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0mierr\u001b[0m \u001b[1;33m<\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 134\u001b[1;33m \u001b[1;32mraise\u001b[0m \u001b[0mSolverException\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"Mumps Exception [%d] - %s\"\u001b[0m \u001b[1;33m%\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mierr\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0m_mumpsErrors\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mierr\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 135\u001b[0m \u001b[1;32melif\u001b[0m \u001b[0mierr\u001b[0m \u001b[1;33m>\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 136\u001b[0m \u001b[1;32mprint\u001b[0m \u001b[1;34m\"Mumps Warning [%d] - %s\"\u001b[0m \u001b[1;33m%\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mierr\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0m_mumpsErrors\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mierr\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[1;31mSolverException\u001b[0m: Mumps Exception [-13] - An error occurred in a Fortran ALLOCATE statement. The size that the package requested is available in INFO(2). If INFO(2) is negative, then the size that the package requested is obtained by multiplying the absolute value of INFO(2) by 1 million."
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# Make the receiver list\n",
|
||||
"rxList = []\n",
|
||||
"for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi']:\n",
|
||||
" rxList.append(simpegmt.SurveyMT.RxMT(locs,rxType)) \n",
|
||||
"# Source list\n",
|
||||
"srcList =[]\n",
|
||||
"freqs = np.logspace(3,0,13)\n",
|
||||
"for freq in freqs:\n",
|
||||
" srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freq))\n",
|
||||
"# Survey MT\n",
|
||||
"survey = simpegmt.SurveyMT.SurveyMT(srcList)\n",
|
||||
"\n",
|
||||
"# Setup the problem object\n",
|
||||
"sigma1d = mesh3d.r(modelTD,'CC','CC','M')[0,0,:] # Use the edge column as a background model\n",
|
||||
"problem = simpegmt.ProblemMT3D.eForm_ps(mesh3d,sigmaPrimary = sigma1d)\n",
|
||||
"problem.verbose = False\n",
|
||||
"from pymatsolver import MumpsSolver\n",
|
||||
"problem.Solver = MumpsSolver\n",
|
||||
"problem.pair(survey)\n",
|
||||
"\n",
|
||||
"# Forward model the data\n",
|
||||
"fields = problem.fields(modelTD)\n",
|
||||
"mtData = survey.projectFields(fields)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
" 50*50*48"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"np.sum(modTD<1e-7)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"45000/50**2"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
@@ -18,9 +187,7 @@
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Define the model\n"
|
||||
]
|
||||
"source": []
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
|
||||
@@ -0,0 +1,71 @@
|
||||
## Script to run 3D forward model
|
||||
## Forward model a data
|
||||
import numpy as np, sys, os, time, gzip, cPickle as pickle
|
||||
sys.path.append('/tera_raid/gudni/gitCodes/simpegmt')
|
||||
sys.path.append('/tera_raid/gudni/gitCodes/simpegem')
|
||||
sys.path.append('/tera_raid/gudni/gitCodes/simpeg')
|
||||
sys.path.append('/tera_raid/gudni')
|
||||
from pymatsolver import MumpsSolver
|
||||
import simpegMT as simpegmt, SimPEG as simpeg
|
||||
import numpy as np, scipy
|
||||
## Setup the forward modeling
|
||||
# Read the model
|
||||
modelname = "simpegTDmodel.con"
|
||||
sigma = np.loadtxt(modelname)
|
||||
# Make the mesh.
|
||||
mTensor = simpeg.Utils.meshTensor
|
||||
cSize = [50,20]
|
||||
# Cells constant size mesh
|
||||
hx = mTensor([(cSize[0],50)])
|
||||
hy = mTensor([(cSize[0],50)])
|
||||
hz = mTensor([(cSize[1],48)])
|
||||
x0 = np.array([-1250,-1250,- 30*20])
|
||||
mesh3dCons = simpeg.Mesh.TensorMesh([hx,hy,hz],x0)
|
||||
# With padding
|
||||
hPad = mTensor([(cSize[0],5,1.5)])
|
||||
aPad = mTensor([(cSize[1],13,1.3)])
|
||||
bPad = mTensor([(cSize[1],5,-1.5)])
|
||||
hxPad = np.hstack((hPad[::-1],mTensor([(cSize[0],40)]),hPad))
|
||||
hyPad = np.hstack((hPad[::-1],mTensor([(cSize[0],40)]),hPad))
|
||||
hzPad = np.hstack((bPad,mTensor([(cSize[1],30)]),aPad))
|
||||
x0Pad = np.array([-(np.sum(hPad)+1000),-(np.sum(hPad)+1000),-(np.sum(bPad)+600)])
|
||||
mesh3d = simpeg.Mesh.TensorMesh([hxPad,hyPad,hzPad],x0Pad)
|
||||
|
||||
# Load the model to the uniform cell mesh
|
||||
modelUniCell = simpeg.Utils.meshutils.readUBCTensorModel(modelname,mesh3dCons)
|
||||
|
||||
# Load the model to the mesh with padding cells
|
||||
modelTD = simpeg.Utils.meshutils.readUBCTensorModel(modelname,mesh3d)
|
||||
|
||||
|
||||
# Define the data locations
|
||||
xG,yG = np.meshgrid(np.linspace(-700,700,8),np.linspace(-700,700,8))
|
||||
zG = np.zeros_like(xG)
|
||||
locs = np.hstack((simpeg.mkvc(xG.ravel(),2),simpeg.mkvc(yG.ravel(),2),simpeg.mkvc(zG.ravel(),2)))
|
||||
|
||||
# Make the receiver list
|
||||
rxList = []
|
||||
for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi']:
|
||||
rxList.append(simpegmt.SurveyMT.RxMT(locs,rxType))
|
||||
# Source list
|
||||
srcList =[]
|
||||
freqs = np.logspace(3,0,13)
|
||||
for freq in freqs:
|
||||
srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freq))
|
||||
# Survey MT
|
||||
survey = simpegmt.SurveyMT.SurveyMT(srcList)
|
||||
|
||||
# Setup the problem object
|
||||
sigma1d = mesh3d.r(modelTD,'CC','CC','M')[0,0,:] # Use the edge column as a background model
|
||||
problem = simpegmt.ProblemMT3D.eForm_ps(mesh3d,sigmaPrimary = sigma1d)
|
||||
problem.verbose = False
|
||||
from pymatsolver import MumpsSolver
|
||||
problem.Solver = MumpsSolver
|
||||
problem.pair(survey)
|
||||
|
||||
# Forward model the data
|
||||
fields = problem.fields(modelTD)
|
||||
mtData = survey.projectFields(fields)
|
||||
# Save the data
|
||||
np.save('seogiModel_MTdata.npy',simpeg.mkvc(mtData,1))
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
@@ -23,40 +23,19 @@ class eForm_ps(BaseMTProblem):
|
||||
_fieldType = 'e'
|
||||
_eqLocs = 'FE'
|
||||
fieldsPair = FieldsMT_3D
|
||||
# Need to add the src ....
|
||||
|
||||
|
||||
# Set new properties
|
||||
# Background model
|
||||
# Shouldn't need the commented block.
|
||||
# @property
|
||||
# def backModel(self):
|
||||
# """
|
||||
# Sets the model, and removes dependent mass matrices.
|
||||
# """
|
||||
# return getattr(self, '_backModel', None)
|
||||
|
||||
# @backModel.setter
|
||||
# def backModel(self, value):
|
||||
# if value is self.backModel:
|
||||
# return # it is the same!
|
||||
# self._backModel = Models.Model(value, self.mapping)
|
||||
# for prop in self.deleteTheseOnModelUpdate:
|
||||
# if hasattr(self, prop):
|
||||
# delattr(self, prop)
|
||||
|
||||
# @property
|
||||
# def MeDeltaSigma(self):
|
||||
# #TODO: hardcoded to sigma as the model
|
||||
# if getattr(self, '_MeDeltaSigma', None) is None:
|
||||
# sigma = self.curModel
|
||||
# sigmaBG = self.backModel
|
||||
# self._MeDeltaSigma = self.mesh.getEdgeInnerProduct(sigma - sigmaBG)
|
||||
# return self._MeDeltaSigma
|
||||
_sigmaPrimary = None
|
||||
|
||||
def __init__(self, mesh, **kwargs):
|
||||
BaseMTProblem.__init__(self, mesh, **kwargs)
|
||||
|
||||
@property
|
||||
def sigmaPrimary(self):
|
||||
return self._sigmaPrimary
|
||||
@sigmaPrimary.setter
|
||||
def sigmaPrimary(self, val):
|
||||
# Note: TODO add logic for val, make sure it is the correct size.
|
||||
self._sigmaPrimary = val
|
||||
|
||||
def getA(self, freq):
|
||||
"""
|
||||
Function to get the A matrix.
|
||||
|
||||
@@ -36,7 +36,7 @@ def setupSurvey(sigmaHalf,tD=True):
|
||||
srcList.append(simpegmt.SurveyMT.srcMT_polxy_1DhomotD(rxList,freq))
|
||||
else:
|
||||
for freq in freqs:
|
||||
srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freq,sigma))
|
||||
srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freq))
|
||||
|
||||
survey = simpegmt.SurveyMT.SurveyMT(srcList)
|
||||
return survey, sigma, m1d
|
||||
@@ -63,7 +63,7 @@ def appRes_TotalFieldNorm(sigmaHalf):
|
||||
|
||||
# Make the survey
|
||||
survey, sigma, mesh = setupSurvey(sigmaHalf)
|
||||
problem = simpegmt.ProblemMT1D.eForm_TotalField(mesh,sigma)
|
||||
problem = simpegmt.ProblemMT1D.eForm_TotalField(mesh)
|
||||
problem.pair(survey)
|
||||
|
||||
# Get the fields
|
||||
@@ -99,7 +99,7 @@ def appRes_psFieldNorm(sigmaHalf):
|
||||
|
||||
# Make the survey
|
||||
survey, sigma, mesh = setupSurvey(sigmaHalf,False)
|
||||
problem = simpegmt.ProblemMT1D.eForm_psField(mesh)
|
||||
problem = simpegmt.ProblemMT1D.eForm_psField(mesh, sigmaPrimary = sigma)
|
||||
problem.pair(survey)
|
||||
|
||||
# Get the fields
|
||||
@@ -117,7 +117,7 @@ def appPhs_psFieldNorm(sigmaHalf):
|
||||
|
||||
# Make the survey
|
||||
survey, sigma, mesh = setupSurvey(sigmaHalf,False)
|
||||
problem = simpegmt.ProblemMT1D.eForm_psField(mesh)
|
||||
problem = simpegmt.ProblemMT1D.eForm_psField(mesh, sigmaPrimary = sigma)
|
||||
problem.pair(survey)
|
||||
|
||||
# Get the fields
|
||||
|
||||
@@ -109,7 +109,7 @@ def dataMis_AnalyticPrimarySecondary(sigmaHalf):
|
||||
# Make the survey
|
||||
# Primary secondary
|
||||
surveyPS, sigmaPS, mesh = setupSurvey(sigmaHalf,False)
|
||||
problemPS = simpegmt.ProblemMT1D.eForm_psField(mesh,sigma)
|
||||
problemPS = simpegmt.ProblemMT1D.eForm_psField(mesh,sigmaPS)
|
||||
problemPS.pair(surveyPS)
|
||||
# Analytic data
|
||||
dataAna = calculateAnalyticSolution(surveyPS.srcList,mesh,sigma)
|
||||
|
||||
@@ -18,7 +18,7 @@ def getInputs():
|
||||
# M = simpeg.Mesh.TensorMesh([[(100,5,-1.5),(100.,10),(100,5,1.5)],[(100,5,-1.5),(100.,10),(100,5,1.5)],[(100,5,1.6),(100.,10),(100,3,2)]], x0=['C','C',-3529.5360])
|
||||
M = simpeg.Mesh.TensorMesh([[(1000,6,-1.5),(1000.,6),(1000,6,1.5)],[(1000,6,-1.5),(1000.,2),(1000,6,1.5)],[(1000,10,-1.3),(1000.,2),(1000,10,1.3)]], x0=['C','C','C'])# Setup the model
|
||||
# Set the frequencies
|
||||
freqs = np.logspace(3,-3,7)
|
||||
freqs = np.logspace(1,-3,5)
|
||||
elev = 0
|
||||
|
||||
## Setup the the survey object
|
||||
@@ -73,15 +73,18 @@ def runSimpegMTfwd_eForm_ps(inputsProblem):
|
||||
srcList =[]
|
||||
sigma1d = M.r(sigBG,'CC','CC','M')[0,0,:]
|
||||
for freq in freqs:
|
||||
srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freq,sigma1d))
|
||||
srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freq))
|
||||
# Survey MT
|
||||
survey = simpegmt.SurveyMT.SurveyMT(srcList)
|
||||
|
||||
## Setup the problem object
|
||||
problem = simpegmt.ProblemMT3D.eForm_ps(M)
|
||||
problem = simpegmt.ProblemMT3D.eForm_ps(M,sigmaPrimary=sigma1d)
|
||||
problem.verbose = False
|
||||
from pymatsolver import MumpsSolver
|
||||
problem.Solver = MumpsSolver
|
||||
try:
|
||||
from pymatsolver import MumpsSolver
|
||||
problem.Solver = MumpsSolver
|
||||
except:
|
||||
pass
|
||||
problem.pair(survey)
|
||||
|
||||
fields = problem.fields(sig)
|
||||
@@ -106,19 +109,19 @@ def appResPhsHalfspace_eFrom_ps_Norm(sigmaHalf,appR=True):
|
||||
# Calculate the app phs
|
||||
app_rpxy, app_rpyx = np.array(getAppResPhs(data))
|
||||
if appR:
|
||||
return np.linalg.norm(np.abs(app_rpxy[0,:] - np.ones(survey.nFreq)/sigmaHalf) * sigmaHalf)
|
||||
return np.all(np.abs(app_rpxy[0,:] - np.ones(survey.nFreq)/sigmaHalf) * sigmaHalf < .35)
|
||||
else:
|
||||
return np.linalg.norm(np.abs(app_rpxy[1,:] + np.ones(survey.nFreq)*135) / 135)
|
||||
return np.all(np.abs(app_rpxy[1,:] + np.ones(survey.nFreq)*135) / 135 < .35)
|
||||
|
||||
class TestAnalytics(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
pass
|
||||
def test_appRes2en1(self):self.assertLess(appResPhsHalfspace_eFrom_ps_Norm(2e-1), TOLr)
|
||||
def test_appRes2en2(self):self.assertLess(appResPhsHalfspace_eFrom_ps_Norm(2e-2), TOLr)
|
||||
def test_appRes2en3(self):self.assertLess(appResPhsHalfspace_eFrom_ps_Norm(2e-3), TOLr)
|
||||
def test_appRes2en1(self):self.assertLess(appResPhsHalfspace_eFrom_ps_Norm(2e-1,False), TOLr)
|
||||
def test_appRes2en2(self):self.assertLess(appResPhsHalfspace_eFrom_ps_Norm(2e-2,False), TOLr)
|
||||
def test_appRes2en3(self):self.assertLess(appResPhsHalfspace_eFrom_ps_Norm(2e-3,False), TOLr)
|
||||
# def test_appRes2en1(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(2e-1))
|
||||
def test_appRes1en2(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-2))
|
||||
def test_appRes1en3(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-3))
|
||||
# def test_appRes2en1(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(2e-1,False))
|
||||
def test_appPhs1en2(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-2,False))
|
||||
def test_appPhs1en3(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-3,False))
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -65,7 +65,6 @@ def plotMT1DModelData(problem,models,symList=None):
|
||||
|
||||
# if not symList:
|
||||
# symList = ['x']*len(models)
|
||||
sys.path.append('/home/gudni/Dropbox/code/python/MTview')
|
||||
import plotDataTypes as pDt
|
||||
# Loop through the models.
|
||||
modelList = [problem.survey.mtrue]
|
||||
@@ -122,4 +121,8 @@ def plotMT1DModelData(problem,models,symList=None):
|
||||
for ax in [axM,axR,axP]:
|
||||
ax.xaxis.set_tick_params(labelsize=fontSize)
|
||||
ax.yaxis.set_tick_params(labelsize=fontSize)
|
||||
return fig
|
||||
return fig
|
||||
|
||||
def printTime():
|
||||
import time
|
||||
print time.strftime("%a, %d %b %Y %H:%M:%S +0000", time.localtime())
|
||||
@@ -7,7 +7,12 @@ from simpegMT.Utils.dataUtils import rec2ndarr
|
||||
|
||||
# Import modules
|
||||
import numpy as np
|
||||
import os, osr, sys, re
|
||||
import os, sys, re
|
||||
try:
|
||||
import osr
|
||||
except ImportError as e:
|
||||
print 'Could not import osr, missing the gdal package'
|
||||
pass
|
||||
|
||||
class EDIimporter:
|
||||
"""
|
||||
|
||||
@@ -0,0 +1,416 @@
|
||||
from matplotlib import pyplot as plt, colors, numpy as np
|
||||
|
||||
|
||||
def rec2nd(structArray):
|
||||
""" Converts a structured/record array to ndarray to do operations on."""
|
||||
return structArray.view((np.float,len(structArray.dtype.names)))
|
||||
|
||||
def plotIsoFreqNSimpedance(ax,freq,array,flag,par='abs',colorbar=True,colorNorm='SymLog',cLevel=True,contour=True):
|
||||
|
||||
indUniFreq = np.where(freq==array['freq'])
|
||||
|
||||
|
||||
x, y = array['x'][indUniFreq],array['y'][indUniFreq]
|
||||
if par == 'abs':
|
||||
zPlot = np.abs(array[flag][indUniFreq])
|
||||
cmap = plt.get_cmap('OrRd_r')#seismic')
|
||||
level = np.logspace(0,-5,31)
|
||||
clevel = np.logspace(0,-4,5)
|
||||
plotNorm = colors.LogNorm()
|
||||
elif par == 'real':
|
||||
zPlot = np.real(array[flag][indUniFreq])
|
||||
cmap = plt.get_cmap('RdYlBu')
|
||||
if cLevel:
|
||||
level = np.concatenate((-np.logspace(0,-10,31),np.logspace(-10,0,31)))
|
||||
clevel = np.concatenate((-np.logspace(0,-8,5),np.logspace(-8,0,5)))
|
||||
else:
|
||||
level = np.linspace(zPlot.min(),zPlot.max(),100)
|
||||
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
|
||||
if colorNorm=='SymLog':
|
||||
plotNorm = colors.SymLogNorm(1e-10,linscale=2)
|
||||
else:
|
||||
plotNorm = colors.Normalize()
|
||||
elif par == 'imag':
|
||||
zPlot = np.imag(array[flag][indUniFreq])
|
||||
cmap = plt.get_cmap('RdYlBu')
|
||||
level = np.concatenate((-np.logspace(0,-10,31),np.logspace(-10,0,31)))
|
||||
clevel = np.concatenate((-np.logspace(0,-8,5),np.logspace(-8,0,5)))
|
||||
plotNorm = colors.SymLogNorm(1e-10,linscale=2)
|
||||
if cLevel:
|
||||
level = np.concatenate((-np.logspace(0,-10,31),np.logspace(-10,0,31)))
|
||||
clevel = np.concatenate((-np.logspace(0,-8,5),np.logspace(-8,0,5)))
|
||||
else:
|
||||
level = np.linspace(zPlot.min(),zPlot.max(),100)
|
||||
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
|
||||
if colorNorm=='SymLog':
|
||||
plotNorm = colors.SymLogNorm(1e-10,linscale=2)
|
||||
elif colorNorm=='Lin':
|
||||
plotNorm = colors.Normalize()
|
||||
if contour:
|
||||
cs = ax.tricontourf(x,y,zPlot,levels=level,cmap=cmap,norm=plotNorm)#,extend='both')
|
||||
else:
|
||||
uniX,uniY = np.unique(x),np.unique(y)
|
||||
X,Y = np.meshgrid(np.append(uniX-25,uniX[-1]+25),np.append(uniY-25,uniY[-1]+25))
|
||||
cs = ax.pcolor(X,Y,np.reshape(zPlot,(len(uniY),len(uniX))),cmap=cmap,norm=plotNorm)
|
||||
if colorbar:
|
||||
plt.colorbar(cs,cax=ax.cax,ticks=clevel,format='%1.2e')
|
||||
ax.set_title(flag+' '+par,fontsize=8)
|
||||
return cs
|
||||
|
||||
def plotIsoFreqNSDiff(ax,freq,arrayList,flag,par='abs',colorbar=True,cLevel=True,mask=None,contourLine=True,useLog=False):
|
||||
|
||||
indUniFreq0 = np.where(freq==arrayList[0]['freq'])
|
||||
indUniFreq1 = np.where(freq==arrayList[1]['freq'])
|
||||
seicmap = plt.get_cmap('RdYlBu')#seismic')
|
||||
x, y = arrayList[0]['x'][indUniFreq0],arrayList[0]['y'][indUniFreq0]
|
||||
if par == 'abs':
|
||||
if useLog:
|
||||
zPlot = (np.log10(np.abs(arrayList[0][flag][indUniFreq0])) - np.log10(np.abs(arrayList[1][flag][indUniFreq1])))/np.log10(np.abs(arrayList[1][flag][indUniFreq1]))
|
||||
else:
|
||||
zPlot = (np.abs(arrayList[0][flag][indUniFreq0]) - np.abs(arrayList[1][flag][indUniFreq1]))/np.abs(arrayList[1][flag][indUniFreq1])
|
||||
if mask:
|
||||
maskInd = np.logical_or(np.abs(arrayList[0][flag][indUniFreq0])< 1e-3,np.abs(arrayList[1][flag][indUniFreq1]) < 1e-3)
|
||||
zPlot = np.ma.array(zPlot)
|
||||
zPlot[maskInd] = mask
|
||||
if cLevel:
|
||||
level = np.arange(-200,201,10)
|
||||
clevel = np.arange(-200,201,25)
|
||||
else:
|
||||
level = np.linspace(zPlot.min(),zPlot.max(),100)
|
||||
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
|
||||
elif par == 'real':
|
||||
if useLog:
|
||||
zPlot = (np.log10(np.real(arrayList[0][flag][indUniFreq0])) -np.log10(np.real(arrayList[1][flag][indUniFreq1])))/np.log10(np.abs((np.real(arrayList[1][flag][indUniFreq1]))))
|
||||
else:
|
||||
zPlot = (np.real(arrayList[0][flag][indUniFreq0]) -np.real(arrayList[1][flag][indUniFreq1]))/np.abs((np.real(arrayList[1][flag][indUniFreq1])))
|
||||
if mask:
|
||||
maskInd = np.logical_or(np.abs(np.real(arrayList[0][flag][indUniFreq0])) < 1e-3,np.abs(np.real(arrayList[1][flag][indUniFreq1])) < 1e-3)
|
||||
zPlot = np.ma.array(zPlot)
|
||||
zPlot[maskInd] = mask
|
||||
if cLevel:
|
||||
level = np.arange(-200,201,10)
|
||||
clevel = np.arange(-200,201,25)
|
||||
else:
|
||||
level = np.linspace(zPlot.min(),zPlot.max(),100)
|
||||
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
|
||||
elif par == 'imag':
|
||||
if useLog:
|
||||
zPlot = (np.log10(np.imag(arrayList[0][flag][indUniFreq0])) -np.log10(np.imag(arrayList[1][flag][indUniFreq1])))/np.log10(np.abs((np.imag(arrayList[1][flag][indUniFreq1]))))
|
||||
else:
|
||||
zPlot = (np.imag(arrayList[0][flag][indUniFreq0]) -np.imag(arrayList[1][flag][indUniFreq1]))/np.abs((np.imag(arrayList[1][flag][indUniFreq1])))
|
||||
if mask:
|
||||
maskInd = np.logical_or(np.abs(np.imag(arrayList[0][flag][indUniFreq0])) < 1e-3,np.abs(np.imag(arrayList[1][flag][indUniFreq1])) < 1e-3)
|
||||
zPlot = np.ma.array(zPlot)
|
||||
zPlot[maskInd] = mask
|
||||
if cLevel:
|
||||
level = np.arange(-200,201,10)
|
||||
clevel = np.arange(-200,201,25)
|
||||
else:
|
||||
level = np.linspace(zPlot.min(),zPlot.max(),100)
|
||||
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
|
||||
cs = ax.tricontourf(x,y,zPlot*100,levels=level*100,cmap=seicmap,extend='both') #,norm=colors.SymLogNorm(1e-2,linscale=2))
|
||||
if contourLine:
|
||||
csl = ax.tricontour(x,y,zPlot*100,levels=clevel*100,colors='k')
|
||||
plt.clabel(csl, fontsize=7, inline=1,fmt='%1.1e',inline_spacing=10)
|
||||
if colorbar:
|
||||
cb = plt.colorbar(cs,cax=ax.cax,ticks=clevel*100,format='%1.1e')
|
||||
for t in cb.ax.get_yticklabels():
|
||||
t.set_rotation(60)
|
||||
t.set_fontsize(8)
|
||||
|
||||
ax.set_title(flag+' '+par,fontsize=8)
|
||||
|
||||
def plotIsoFreqNStipper(ax,freq,array,flag,par='abs',colorbar=True,colorNorm='SymLog',cLevel=True,contour=True):
|
||||
|
||||
indUniFreq = np.where(freq==array['freq'])
|
||||
|
||||
x, y = array['x'][indUniFreq],array['y'][indUniFreq]
|
||||
if par == 'abs':
|
||||
cmap = plt.get_cmap('OrRd_r')#seismic')
|
||||
zPlot = np.abs(array[flag][indUniFreq])
|
||||
if cLevel:
|
||||
level = np.logspace(-4,0,33)
|
||||
clevel = np.logspace(-4,0,5)
|
||||
else:
|
||||
level = np.linspace(zPlot.min(),zPlot.max(),100)
|
||||
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
|
||||
if colorNorm=='SymLog':
|
||||
plotNorm = colors.LogNorm()
|
||||
else:
|
||||
plotNorm = colors.Normalize()
|
||||
elif par == 'real':
|
||||
cmap = plt.get_cmap('RdYlBu')
|
||||
zPlot = np.real(array[flag][indUniFreq])
|
||||
if cLevel:
|
||||
level = np.concatenate((-np.logspace(0,-4,33),np.logspace(-4,0,33)))
|
||||
clevel = np.concatenate((-np.logspace(0,-4,5),np.logspace(-4,0,5)))
|
||||
else:
|
||||
level = np.linspace(zPlot.min(),zPlot.max(),100)
|
||||
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
|
||||
if colorNorm=='SymLog':
|
||||
plotNorm = colors.SymLogNorm(1e-4,linscale=2)
|
||||
else:
|
||||
plotNorm = colors.Normalize()
|
||||
elif par == 'imag':
|
||||
cmap = plt.get_cmap('RdYlBu')
|
||||
zPlot = np.imag(array[flag][indUniFreq])
|
||||
if cLevel:
|
||||
level = np.concatenate((-np.logspace(0,-4,33),np.logspace(-4,0,33)))
|
||||
clevel = np.concatenate((-np.logspace(0,-4,5),np.logspace(-4,0,5)))
|
||||
else:
|
||||
level = np.linspace(zPlot.min(),zPlot.max(),100)
|
||||
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
|
||||
if colorNorm=='SymLog':
|
||||
plotNorm = colors.SymLogNorm(1e-4,linscale=2)
|
||||
else:
|
||||
plotNorm = colors.Normalize()
|
||||
if contour:
|
||||
cs = ax.tricontourf(x,y,zPlot,levels=level,cmap=cmap,norm=plotNorm)#,extend='both')
|
||||
else:
|
||||
uniX,uniY = np.unique(x),np.unique(y)
|
||||
X,Y = np.meshgrid(np.append(uniX-25,uniX[-1]+25),np.append(uniY-25,uniY[-1]+25))
|
||||
cs = ax.pcolor(X,Y,np.reshape(zPlot,(len(uniY),len(uniX))),levels=level,cmap=cmap,norm=plotNorm,edgecolors='k', linewidths=0.5)
|
||||
if colorbar:
|
||||
plt.colorbar(cs,cax=ax.cax,ticks=clevel,format='%1.2e')
|
||||
ax.set_title(flag+' '+par,fontsize=8)
|
||||
|
||||
def plotIsoStaImpedance(ax,loc,array,flag,par='abs',pSym='s',pColor=None):
|
||||
|
||||
appResFact = 1/(8*np.pi**2*10**(-7))
|
||||
treshold = 1.0 # 1 meter
|
||||
indUniSta = np.sqrt(np.sum((rec2nd(array[['x','y']])-loc)**2,axis=1)) < treshold
|
||||
freq = array['freq'][indUniSta]
|
||||
|
||||
if par == 'abs':
|
||||
zPlot = np.abs(array[flag][indUniSta])
|
||||
elif par == 'real':
|
||||
zPlot = np.real(array[flag][indUniSta])
|
||||
elif par == 'imag':
|
||||
zPlot = np.imag(array[flag][indUniSta])
|
||||
elif par == 'res':
|
||||
zPlot = (appResFact/freq)*np.abs(array[flag][indUniSta])**2
|
||||
elif par == 'phs':
|
||||
zPlot = np.arctan2(array[flag][indUniSta].imag,array[flag][indUniSta].real)*(180/np.pi)
|
||||
|
||||
if not pColor:
|
||||
if 'xx' in flag:
|
||||
lab = 'XX'
|
||||
pColor = 'g'
|
||||
elif 'xy' in flag:
|
||||
lab = 'XY'
|
||||
pColor = 'r'
|
||||
elif 'yx' in flag:
|
||||
lab = 'YX'
|
||||
pColor = 'b'
|
||||
elif 'yy' in flag:
|
||||
lab = 'YY'
|
||||
pColor = 'y'
|
||||
|
||||
ax.plot(freq,zPlot,color=pColor,marker=pSym,label=flag)
|
||||
|
||||
|
||||
def plotPsudoSectNSimpedance(ax,sectDict,array,flag,par='abs',colorbar=True,colorNorm='None',cLevel=None,contour=True):
|
||||
|
||||
indSect = np.where(sectDict.values()[0]==array[sectDict.keys()[0]])
|
||||
|
||||
# Define the plot axes
|
||||
if 'x' in sectDict.keys()[0]:
|
||||
x = array['y'][indSect]
|
||||
else:
|
||||
x = array['x'][indSect]
|
||||
y = array['freq'][indSect]
|
||||
|
||||
if par == 'abs':
|
||||
zPlot = np.abs(array[flag][indSect])
|
||||
cmap = plt.get_cmap('OrRd_r')#seismic')
|
||||
if cLevel:
|
||||
level = np.logspace(0,-5,31,endpoint=True)
|
||||
clevel = np.logspace(0,-4,5,endpoint=True)
|
||||
else:
|
||||
level = np.linspace(zPlot.min(),zPlot.max(),100,endpoint=True)
|
||||
clevel = np.linspace(zPlot.min(),zPlot.max(),10,endpoint=True)
|
||||
|
||||
elif par == 'ares':
|
||||
zPlot = np.abs(array[flag][indSect])**2/(8*np.pi**2*10**(-7)*array['freq'][indSect])
|
||||
cmap = plt.get_cmap('RdYlBu')#seismic)
|
||||
if cLevel:
|
||||
zMax = np.log10(cLevel[1])
|
||||
zMin = np.log10(cLevel[0])
|
||||
else:
|
||||
zMax = (np.ceil(np.log10(np.abs(zPlot).max())))
|
||||
zMin = (np.floor(np.log10(np.abs(zPlot).min())))
|
||||
level = np.logspace(zMin,zMax,(zMax-zMin)*8+1,endpoint=True)
|
||||
clevel = np.logspace(zMin,zMax,(zMax-zMin)*2+1,endpoint=True)
|
||||
plotNorm = colors.LogNorm()
|
||||
|
||||
elif par == 'aphs':
|
||||
zPlot = np.arctan2(array[flag][indSect].imag,array[flag][indSect].real)*(180/np.pi)
|
||||
cmap = plt.get_cmap('RdYlBu')#seismic)
|
||||
if cLevel:
|
||||
zMax = cLevel[1]
|
||||
zMin = cLevel[0]
|
||||
else:
|
||||
zMax = (np.ceil(zPlot).max())
|
||||
zMin = (np.floor(zPlot).min())
|
||||
level = np.arange(zMin,zMax+.1,1)
|
||||
clevel = np.arange(zMin,zMax+.1,10)
|
||||
plotNorm = colors.Normalize()
|
||||
|
||||
elif par == 'real':
|
||||
zPlot = np.real(array[flag][indSect])
|
||||
cmap = plt.get_cmap('Spectral') #('RdYlBu')
|
||||
if cLevel:
|
||||
zMax = np.log10(cLevel[1])
|
||||
zMin = np.log10(cLevel[0])
|
||||
else:
|
||||
zMax = (np.ceil(np.log10(np.abs(zPlot).max())))
|
||||
zMin = (np.floor(np.log10(np.abs(zPlot).min())))
|
||||
level = np.concatenate((-np.logspace(zMax,zMin-.125,(zMax-zMin)*8+1,endpoint=True),np.logspace(zMin-.125,zMax,(zMax-zMin)*8+1,endpoint=True)))
|
||||
clevel = np.concatenate((-np.logspace(zMax,zMin,(zMax-zMin)*1+1,endpoint=True),np.logspace(zMin,zMax,(zMax-zMin)*1+1,endpoint=True)))
|
||||
plotNorm = colors.SymLogNorm(np.abs(level).min(),linscale=0.1)
|
||||
elif par == 'imag':
|
||||
zPlot = np.imag(array[flag][indSect])
|
||||
cmap = plt.get_cmap('Spectral') #('RdYlBu')
|
||||
|
||||
if cLevel:
|
||||
zMax = np.log10(cLevel[1])
|
||||
zMin = np.log10(cLevel[0])
|
||||
else:
|
||||
zMax = (np.ceil(np.log10(np.abs(zPlot).max())))
|
||||
zMin = (np.floor(np.log10(np.abs(zPlot).min())))
|
||||
level = np.concatenate((-np.logspace(zMax,zMin-.125,(zMax-zMin)*8+1,endpoint=True),np.logspace(zMin-.125,zMax,(zMax-zMin)*8+1,endpoint=True)))
|
||||
clevel = np.concatenate((-np.logspace(zMax,zMin,(zMax-zMin)*1+1,endpoint=True),np.logspace(zMin,zMax,(zMax-zMin)*1+1,endpoint=True)))
|
||||
plotNorm = colors.SymLogNorm(np.abs(level).min(),linscale=0.1)
|
||||
|
||||
if colorNorm=='SymLog':
|
||||
plotNorm = colors.SymLogNorm(np.abs(level).min(),linscale=0.1)
|
||||
elif colorNorm=='Lin':
|
||||
plotNorm = colors.Normalize()
|
||||
elif colorNorm=='Log':
|
||||
plotNorm = colors.LogNorm()
|
||||
if contour:
|
||||
cs = ax.tricontourf(x,y,zPlot,levels=level,cmap=cmap,norm=plotNorm)#,extend='both')
|
||||
else:
|
||||
uniX,uniY = np.unique(x),np.unique(y)
|
||||
X,Y = np.meshgrid(np.append(uniX-25,uniX[-1]+25),np.append(uniY-25,uniY[-1]+25))
|
||||
cs = ax.pcolor(X,Y,np.reshape(zPlot,(len(uniY),len(uniX))),cmap=cmap,norm=plotNorm)
|
||||
if colorbar:
|
||||
csB = plt.colorbar(cs,cax=ax.cax,ticks=clevel,format='%1.2e')
|
||||
# csB.on_mappable_changed(cs)
|
||||
ax.set_title(flag+' '+par,fontsize=8)
|
||||
return cs, csB
|
||||
return cs,None
|
||||
|
||||
|
||||
def plotPsudoSectNSDiff(ax,sectDict,arrayList,flag,par='abs',colorbar=True,colorNorm='SymLog',cLevel=None,contour=True,mask=None,useLog=False):
|
||||
|
||||
def sortInArr(arr):
|
||||
return np.sort(arr,order=['freq','x','y','z'])
|
||||
# Find the index for the slice
|
||||
indSect0 = np.where(sectDict.values()[0]==arrayList[0][sectDict.keys()[0]])
|
||||
indSect1 = np.where(sectDict.values()[0]==arrayList[1][sectDict.keys()[0]])
|
||||
# Extract and sort the mats
|
||||
arr0 = sortInArr(arrayList[0][indSect0])
|
||||
arr1 = sortInArr(arrayList[1][indSect1])
|
||||
|
||||
# Define the plot axes
|
||||
if 'x' in sectDict.keys()[0]:
|
||||
x0 = arr0['y']
|
||||
x1 = arr1['y']
|
||||
else:
|
||||
x0 = arr0['x']
|
||||
x1 = arr1['x']
|
||||
y0 = arr0['freq']
|
||||
y1 = arr1['freq']
|
||||
|
||||
|
||||
if par == 'abs':
|
||||
if useLog:
|
||||
zPlot = (np.log10(np.abs(arr0[flag])) - np.log10(np.abs(arr1[flag])))/np.log10(np.abs(arr1[flag]))
|
||||
else:
|
||||
zPlot = (np.abs(arr0[flag]) - np.abs(arr1[flag]))/np.abs(arr1[flag])
|
||||
if mask:
|
||||
maskInd = np.logical_or(np.abs(arr0[flag])< 1e-3,np.abs(arr1[flag]) < 1e-3)
|
||||
zPlot = np.ma.array(zPlot)
|
||||
zPlot[maskInd] = mask
|
||||
cmap = plt.get_cmap('RdYlBu')#seismic)
|
||||
elif par == 'ares':
|
||||
arF = 1/(8*np.pi**2*10**(-7))
|
||||
if useLog:
|
||||
zPlot = (np.log10((arF/arr0['freq'])*np.abs(arr0[flag])**2) - np.log10((arF/arr1['freq'])*np.abs(arr1[flag])**2))/np.log10((arF/arr1['freq'])*np.abs(arr1[flag])**2)
|
||||
else:
|
||||
zPlot = ((arF/arr0['freq'])*np.abs(arr0[flag])**2 - (arF/arr1['freq'])*np.abs(arr1[flag])**2)/((arF/arr1['freq'])*np.abs(arr1[flag])**2)
|
||||
if mask:
|
||||
maskInd = np.logical_or(np.abs(arr0[flag])< 1e-3,np.abs(arr1[flag]) < 1e-3)
|
||||
zPlot = np.ma.array(zPlot)
|
||||
zPlot[maskInd] = mask
|
||||
cmap = plt.get_cmap('Spectral')#seismic)
|
||||
|
||||
elif par == 'aphs':
|
||||
if useLog:
|
||||
zPlot = (np.log10(np.arctan2(arr0[flag].imag,arr0[flag].real)*(180/np.pi)) - np.log10(np.arctan2(arr1[flag].imag,arr1[flag].real)*(180/np.pi)) )/np.log10(np.arctan2(arr1[flag].imag,arr1[flag].real)*(180/np.pi))
|
||||
else:
|
||||
zPlot = ( np.arctan2(arr0[flag].imag,arr0[flag].real)*(180/np.pi) - np.arctan2(arr1[flag].imag,arr1[flag].real)*(180/np.pi) )/(np.arctan2(arr1[flag].imag,arr1[flag].real)*(180/np.pi))
|
||||
if mask:
|
||||
maskInd = np.logical_or(np.abs(arr0[flag])< 1e-3,np.abs(arr1[flag]) < 1e-3)
|
||||
zPlot = np.ma.array(zPlot)
|
||||
zPlot[maskInd] = mask
|
||||
cmap = plt.get_cmap('Spectral')#seismic)
|
||||
elif par == 'real':
|
||||
if useLog:
|
||||
zPlot = (np.log10(arr0[flag].real) - np.log10(arr1[flag].real))/np.log10(arr1[flag].real)
|
||||
else:
|
||||
zPlot = (arr0[flag].real - arr1[flag].real)/arr1[flag].real
|
||||
if mask:
|
||||
maskInd = np.logical_or(arr0[flag].real< 1e-3,arr1[flag].real < 1e-3)
|
||||
zPlot = np.ma.array(zPlot)
|
||||
zPlot[maskInd] = mask
|
||||
cmap = plt.get_cmap('Spectral') #('Spectral')
|
||||
|
||||
elif par == 'imag':
|
||||
if useLog:
|
||||
zPlot = (np.log10(arr0[flag].imag) - np.log10(arr1[flag].imag))/np.log10(arr1[flag].imag)
|
||||
else:
|
||||
zPlot = (arr0[flag].imag - arr1[flag].imag)/arr1[flag].imag
|
||||
if mask:
|
||||
maskInd = np.logical_or(arr0[flag].imag< 1e-3,arr1[flag].imag < 1e-3)
|
||||
zPlot = np.ma.array(zPlot)
|
||||
zPlot[maskInd] = mask
|
||||
cmap = plt.get_cmap('Spectral') #('RdYlBu')
|
||||
|
||||
if cLevel:
|
||||
zMax = np.log10(cLevel[1])
|
||||
zMin = np.log10(cLevel[0])
|
||||
else:
|
||||
zMax = (np.ceil(np.log10(np.abs(zPlot).max())))
|
||||
zMin = (np.floor(np.log10(np.abs(zPlot).min())))
|
||||
|
||||
|
||||
if colorNorm=='SymLog':
|
||||
level = np.concatenate((-np.logspace(zMax,zMin-.125,(zMax-zMin)*8+1,endpoint=True),np.logspace(zMin-.125,zMax,(zMax-zMin)*8+1,endpoint=True)))
|
||||
clevel = np.concatenate((-np.logspace(zMax,zMin,(zMax-zMin)*1+1,endpoint=True),np.logspace(zMin,zMax,(zMax-zMin)*1+1,endpoint=True)))
|
||||
plotNorm = colors.SymLogNorm(np.abs(level).min(),linscale=0.1)
|
||||
elif colorNorm=='Lin':
|
||||
if cLevel:
|
||||
level = np.arange(cLevel[0],cLevel[1]+.1,(cLevel[1] - cLevel[0])/50.)
|
||||
clevel = np.arange(cLevel[0],cLevel[1]+.1,(cLevel[1] - cLevel[0])/10.)
|
||||
else:
|
||||
level = np.arange(zPlot.min(),zPlot.max(),(zPlot.max() - zPlot.min())/50.)
|
||||
clevel = np.arange(zPlot.min(),zPlot.max(),(zPlot.max() - zPlot.min())/10.)
|
||||
plotNorm = colors.Normalize()
|
||||
elif colorNorm=='Log':
|
||||
level = np.logspace(zMin-.125,zMax,(zMax-zMin)*8+1,endpoint=True)
|
||||
clevel = np.logspace(zMin,zMax,(zMax-zMin)*2+1,endpoint=True)
|
||||
plotNorm = colors.LogNorm()
|
||||
if contour:
|
||||
cs = ax.tricontourf(x0,y0,zPlot*100,levels=level*100,cmap=cmap,norm=plotNorm,extend='both')#,extend='both')
|
||||
else:
|
||||
uniX,uniY = np.unique(x0),np.unique(y0)
|
||||
X,Y = np.meshgrid(np.append(uniX-25,uniX[-1]+25),np.append(uniY-25,uniY[-1]+25))
|
||||
cs = ax.pcolor(X,Y,np.reshape(zPlot,(len(uniY),len(uniX))),cmap=cmap,norm=plotNorm)
|
||||
if colorbar:
|
||||
csB = plt.colorbar(cs,cax=ax.cax,ticks=clevel*100,format='%1.2e')
|
||||
# csB.on_mappable_changed(cs)
|
||||
ax.set_title(flag+' '+par + ' diff',fontsize=8)
|
||||
return cs, csB
|
||||
return cs,None
|
||||
Reference in New Issue
Block a user