mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-30 06:47:05 +08:00
Fixed bugs. Added notebooks with halfspace and layer examples.
Code is working (returning results with in couple of % for a 1D analytic solution) but test need to be "automated".
This commit is contained in:
Gudni Karl Rosenkjaer
+378
-37
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
+961
-260
File diff suppressed because one or more lines are too long
@@ -0,0 +1,128 @@
|
||||
{
|
||||
"metadata": {
|
||||
"name": "",
|
||||
"signature": "sha256:260d0d5c93203805d78b8fa9e20a9079ce0e65e4b6a3453fedfcda305fffaa80"
|
||||
},
|
||||
"nbformat": 3,
|
||||
"nbformat_minor": 0,
|
||||
"worksheets": [
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"import SimPEG as simpeg\n",
|
||||
"import simpegMT as simpegmt\n"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"output_type": "stream",
|
||||
"stream": "stdout",
|
||||
"text": [
|
||||
"Efficiency Warning: Interpolation will be slow, use setup.py!\n",
|
||||
"\n",
|
||||
" python setup.py build_ext --inplace\n",
|
||||
" \n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"prompt_number": 1
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"#Define the mesh\n",
|
||||
"m1d = simpeg.Mesh.TensorMesh([[(100,5,1.5),(100.,10),(100,5,1.5)]], x0=['C'])\n",
|
||||
"sigma = np.zeros(m1d.nC) + 2e-3\n",
|
||||
"sigma[m1d.gridCC[:]>200] = 1e-8\n"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": 6
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"# Calculate the analytic fields\n",
|
||||
"freqs = np.logspace(4,-1,26)\n",
|
||||
"Z = []\n",
|
||||
"for freq in freqs:\n",
|
||||
" Ed, Eu, Hd, Hu = simpegmt.Utils.getEHfields(m1d,sigma,freq,np.array([200]))\n",
|
||||
" Z.append((Ed + Eu)/(Hd + Hu))"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": 17
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"Zarr = np.concatenate(Z)"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": 19
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"def appResPhs(freq,z):\n",
|
||||
" app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2\n",
|
||||
" app_phs = np.arctan2(-z.imag,z.real)*(180/np.pi)\n",
|
||||
" return app_res, app_phs\n",
|
||||
"app_r, app_p = appResPhs(freqs,Zarr)"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": 20
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"app_r"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"metadata": {},
|
||||
"output_type": "pyout",
|
||||
"prompt_number": 21,
|
||||
"text": [
|
||||
"array([ 499.99998067, 499.99999231, 499.99999694, 499.99999878,\n",
|
||||
" 499.99999951, 499.99999981, 499.99999992, 499.99999997,\n",
|
||||
" 499.99999999, 500. , 500. , 500. ,\n",
|
||||
" 500. , 500. , 500. , 500. ,\n",
|
||||
" 500. , 500. , 500. , 500. ,\n",
|
||||
" 500. , 500. , 500. , 500. ,\n",
|
||||
" 500. , 500. ])"
|
||||
]
|
||||
}
|
||||
],
|
||||
"prompt_number": 21
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": []
|
||||
}
|
||||
],
|
||||
"metadata": {}
|
||||
}
|
||||
]
|
||||
}
|
||||
@@ -26,12 +26,12 @@ for loc in rx_loc:
|
||||
rxList.append(simpegmt.SurveyMT.RxMT(simpeg.mkvc(loc,2).T,rxType))
|
||||
# Source list
|
||||
srcList =[]
|
||||
for freq in np.logspace(3,-1,5):
|
||||
for freq in np.logspace(3,-3,7):
|
||||
srcList.append(simpegmt.SurveyMT.srcMT(freq,rxList))
|
||||
# Survey MT
|
||||
survey = simpegmt.SurveyMT.SurveyMT(srcList)
|
||||
|
||||
## Setup the problem objec
|
||||
## Setup the problem object
|
||||
problem = simpegmt.ProblemMT.MTProblem(M)
|
||||
problem.pair(survey)
|
||||
|
||||
|
||||
@@ -1,4 +1,4 @@
|
||||
from SimPEG import Survey, Problem, Utils, Models, np, sp, Solver as SimpegSolver
|
||||
from SimPEG import Survey, Problem, Utils, Models, np, sp, SolverLU as SimpegSolver
|
||||
from scipy.constants import mu_0
|
||||
from SurveyMT import SurveyMT, FieldsMT
|
||||
|
||||
@@ -120,7 +120,8 @@ class MTProblem(Problem.BaseProblem):
|
||||
# Store the fields
|
||||
F[Src, 'e_px'] = e[:,0]
|
||||
F[Src, 'e_py'] = e[:,1]
|
||||
b = self.mesh.edgeCurl * e
|
||||
# Note curl e = -iwb so b = -curl/iw
|
||||
b = -( self.mesh.edgeCurl * e )/( 1j*omega(freq) )
|
||||
F[Src, 'b_px'] = b[:,0]
|
||||
F[Src, 'b_py'] = b[:,1]
|
||||
return F
|
||||
@@ -138,7 +139,7 @@ class MTProblem(Problem.BaseProblem):
|
||||
sig = self.MeSigma
|
||||
C = self.mesh.edgeCurl
|
||||
|
||||
return C.T*mui*C + 1j*omega(freq)*sig
|
||||
return C.T*mui*C - 1j*omega(freq)*sig
|
||||
|
||||
def getAbg(self, freq):
|
||||
"""
|
||||
@@ -152,7 +153,7 @@ class MTProblem(Problem.BaseProblem):
|
||||
sigBG = self.MeSigmaBG
|
||||
C = self.mesh.edgeCurl
|
||||
|
||||
return C.T*mui*C + 1j*omega(freq)*sigBG
|
||||
return C.T*mui*C - 1j*omega(freq)*sigBG
|
||||
|
||||
def getADeriv(self, freq, u, v, adjoint=False):
|
||||
sig = self.curTModel
|
||||
@@ -181,7 +182,7 @@ class MTProblem(Problem.BaseProblem):
|
||||
eBG_bp = homo1DModelSource(self.mesh,freq,backSigma)
|
||||
Abg = self.getAbg(freq)
|
||||
|
||||
return Abg*eBG_bp
|
||||
return -Abg*eBG_bp
|
||||
|
||||
##################################################################
|
||||
# Inversion stuff
|
||||
|
||||
@@ -16,7 +16,8 @@ def homo1DModelSource(mesh,freq,m_back):
|
||||
from simpegMT.Utils import get1DEfields
|
||||
# Get a 1d solution for a halfspace background
|
||||
mesh1d = simpeg.Mesh.TensorMesh([mesh.hz],np.array([mesh.x0[2]]))
|
||||
e0_1d = get1DEfields(mesh1d,mesh.r(m_back,'CC','CC','M')[0,0,:],freq)
|
||||
# Note: Need to conjugate the source field to comply with orientations
|
||||
e0_1d = get1DEfields(mesh1d,mesh.r(m_back,'CC','CC','M')[0,0,:],freq).conj()
|
||||
# Setup x (east) polarization (_x)
|
||||
ex_px = np.zeros(mesh.vnEx,dtype=complex)
|
||||
ey_px = np.zeros((mesh.nEy,1),dtype=complex)
|
||||
@@ -24,18 +25,19 @@ def homo1DModelSource(mesh,freq,m_back):
|
||||
# Assign the source to ex_x
|
||||
for i in np.arange(mesh.vnEx[0]):
|
||||
for j in np.arange(mesh.vnEx[1]):
|
||||
ex_px[i,j,:] = e0_1d
|
||||
ex_px[1:-1,1:-1,1:-1] = 0
|
||||
eBG_px = np.vstack((simpeg.Utils.mkvc(mesh.r(ex_px,'Ex','Ex','V'),2),ey_px,ez_px))
|
||||
ex_px[i,j,:] = -e0_1d
|
||||
# ex_px[1:-1,1:-1,1:-1] = 0
|
||||
eBG_px = np.vstack((simpeg.Utils.mkvc(ex_px,2),ey_px,ez_px))
|
||||
# Setup y (north) polarization (_py)
|
||||
ex_py = np.zeros((mesh.nEx,1), dtype='complex128')
|
||||
ey_py = np.zeros(mesh.vnEy, dtype='complex128')
|
||||
ez_py = np.zeros((mesh.nEz,1), dtype='complex128')
|
||||
# Assign the source to ey_py
|
||||
|
||||
for i in np.arange(mesh.vnEy[0]):
|
||||
for j in np.arange(mesh.vnEy[1]):
|
||||
ey_py[i,j,:] = e0_1d
|
||||
ey_py[1:-1,1:-1,1:-1] = 0
|
||||
ey_py[i,j,:] = e0_1d
|
||||
# ey_py[1:-1,1:-1,1:-1] = 0
|
||||
eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py))
|
||||
|
||||
# Return the electric fields
|
||||
|
||||
@@ -269,7 +269,7 @@ class DataMT(Survey.Data):
|
||||
from numpy.lib import recfunctions as recFunc
|
||||
dt = dtCP
|
||||
for uniFL in uniFLmarr:
|
||||
mTemp = simpeg.mkvc(rec2ndarr(mArrRec[np.ma.where(mArrRec[['freq','x','y','z']].data == np.array(uniFL))][impList]).sum(axis=0),2).T
|
||||
mTemp = mkvc(rec2ndarr(mArrRec[np.ma.where(mArrRec[['freq','x','y','z']].data == np.array(uniFL))][impList]).sum(axis=0),2).T
|
||||
compBlock = np.sum(mTemp.data.reshape((4,2))*np.array([[1,1j],[1,1j],[1,1j],[1,1j]]),axis=1).copy().view(dt[4::])
|
||||
dataBlock = mkvc(recFunc.merge_arrays((np.array(uniFL),compBlock),flatten=True),2).T
|
||||
try:
|
||||
|
||||
@@ -47,12 +47,19 @@ def getEHfields(m1d,sigma,freq,zd):
|
||||
UDp[:,lnr+1] = elamh.dot(Pjinv.dot(Pj1)).dot(UDp[:,lnr])
|
||||
|
||||
# Calculate the fields
|
||||
Ed = np.zeros((zd.size,),dtype=complex)
|
||||
Eu = np.zeros((zd.size,),dtype=complex)
|
||||
Hd = np.zeros((zd.size,),dtype=complex)
|
||||
Hu = np.zeros((zd.size,),dtype=complex)
|
||||
Ed = np.empty((zd.size,),dtype=complex)
|
||||
Eu = np.empty((zd.size,),dtype=complex)
|
||||
Hd = np.empty((zd.size,),dtype=complex)
|
||||
Hu = np.empty((zd.size,),dtype=complex)
|
||||
|
||||
# Loop over the layers and calculate the fields
|
||||
# In the halfspace below the mesh
|
||||
dup = m1d.vectorNx[0]
|
||||
dind = dup >= zd
|
||||
Ed[dind] = UDp[1,0]*np.exp(-1j*k[0]*(dup-zd[dind]))
|
||||
Eu[dind] = UDp[0,0]*np.exp(1j*k[0]*(dup-zd[dind]))
|
||||
Hd[dind] = (k[0]/(w*mu[0]))*UDp[1,0]*np.exp(-1j*k[0]*(dup-zd[dind]))
|
||||
Hu[dind] = -(k[0]/(w*mu[0]))*UDp[0,0]*np.exp(1j*k[0]*(dup-zd[dind]))
|
||||
for ki,mui,epsi,dlow,dup,Up,Dp in zip(k[1::],mu[1::],eps[1::],m1d.vectorNx[:-1],m1d.vectorNx[1::],UDp[0,1::],UDp[1,1::]):
|
||||
dind = np.logical_and(dup >= zd, zd > dlow)
|
||||
Ed[dind] = Dp*np.exp(-1j*ki*(dup-zd[dind]))
|
||||
|
||||
@@ -20,10 +20,12 @@ def get1DEfields(m1d,sigma,freq,sourceAmp=1.0):
|
||||
Aio = A[1:-1,[0,-1]]
|
||||
|
||||
# Set the boundary conditions
|
||||
Ed_low, Eu_low, Hd_low, Hu_low = getEHfields(m1d,sigma,freq,np.array([m1d.vectorNx[0]]))
|
||||
Etot_low = Ed_low + Eu_low
|
||||
Ed, Eu, Hd, Hu = getEHfields(m1d,sigma,freq,m1d.vectorNx)
|
||||
Etot = Ed + Eu
|
||||
if sourceAmp is not None:
|
||||
Etot = ((Etot/Etot[-1])*sourceAmp) # Scale the fields to be equal to sourceAmp at the top
|
||||
## Note: need to use conjugate of the analytic solution. It is derived with e^iwt
|
||||
bc = np.r_[Etot_low.conj(),sourceAmp]
|
||||
bc = np.r_[Etot[0],Etot[-1]]
|
||||
# The right hand side
|
||||
rhs = -Aio*bc
|
||||
# Solve the system
|
||||
|
||||
Reference in New Issue
Block a user