wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-08-16 11:28:21 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #2 from simpeg/moveAnalytics

Browse Source 
Move analytics
This commit is contained in:
Rowan Cockett
2015-02-12 11:40:30 -08:00
parent 7f817c96ad bae79ecb25
commit 48ada2d63c
6 changed files with 111 additions and 254 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simpegMT/Analytics/MT1Danalytic.py
-101
View File
@@ -1,101 +0,0 @@
# Analytic solution of EM fields due to a plane wave
import numpy as np, SimPEG as simpeg
def getEHfields(m1d,sigma,freq,zd):
'''Analytic solution for MT 1D layered earth. Returns E and H fields.
:param SimPEG.mesh, object m1d: Mesh object with the 1D spatial information.
:param numpy.array, vector sigma: Physical property of conductivity corresponding with the mesh.
:param float, freq: Frequency to calculate data at.
:param numpy array, vector zd: location to calculate EH fields at
Assumes a halfspace with the same conductive as the last cell below.
'''
# Note add an error check for the mesh and sigma are the same size.
# Need make the solution e^-iwt dependent
# Constants: Assume constant
mu = 4*np.pi*1e-7*np.ones((m1d.nC+1))
eps = 8.85*1e-12*np.ones((m1d.nC+1))
# Angular freq
w = 2*np.pi*freq
# Add the halfspace value to the property
sig = np.concatenate((np.array([sigma[0]]),sigma))
# Calculate the wave number
k = np.sqrt(eps*mu*w**2-1j*mu*sig*w)
# Initiate the propagation matrix, in the order down up.
UDp = np.zeros((2,m1d.nC+1),dtype=complex)
UDp[1,0] = 1 # Set the wave amplitude as 1 into the half-space at the bottom of the mesh
# Loop over all the layers, starting at the bottom layer
for lnr, h in enumerate(m1d.hx): # lnr-number of layer, h-thickness of the layer
# Calculate
yp1 = k[lnr]/(w*mu[lnr]) # Admittance of the layer below the current layer
zp = (w*mu[lnr+1])/k[lnr+1] # Impedance in the current layer
# Build the propagation matrix
# Convert fields to down/up going components in layer below current layer
Pj1 = np.array([[1,1],[yp1,-yp1]])
# Convert fields to down/up going components in current layer
Pjinv = 1./2*np.array([[1,zp],[1,-zp]])
# Propagate down and up components through the current layer
elamh = np.array([[np.exp(-1j*k[lnr+1]*h),0],[0,np.exp(1j*k[lnr+1]*h)]])
# The down and up component in current layer.
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)
# Loop over the layers and calculate the fields
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]))
Eu[dind] = Up*np.exp(1j*ki*(dup-zd[dind]))
Hd[dind] = (ki/(w*mui))*Dp*np.exp(-1j*ki*(dup-zd[dind]))
Hu[dind] = -(ki/(w*mui))*Up*np.exp(1j*ki*(dup-zd[dind]))
# Return return the fields
return Ed, Eu, Hd, Hu
def getImpedance(m1d,sigma,freq):
"""Analytic solution for MT 1D layered earth. Returns the impedance at the surface.
:param SimPEG.mesh, object m1d: Mesh object with the 1D spatial information.
:param numpy.array, vector sigma: Physical property corresponding with the mesh.
:param numpy.array, vector freq: Frequencies to calculate data at.
"""
# Define constants
mu0 = 4*np.pi*1e-7
eps0 = 8.85e-12
# Initiate the impedances
Z1d = np.empty(len(freq) , dtype='complex')
h = m1d.hx #vectorNx[:-1]
# Start the process
for nrFr, fr in enumerate(freq):
om = 2*np.pi*fr
Zall = np.empty(len(h)+1,dtype='complex')
# Calculate the impedance for the bottom layer
Zall[0] = (mu0*om)/np.sqrt(mu0*eps0*(om)**2 - 1j*mu0*sigma[0]*om)
for nr,hi in enumerate(h):
# Calculate the wave number
# print nr,sigma[nr]
k = np.sqrt(mu0*eps0*om**2 - 1j*mu0*sigma[nr]*om)
Z = (mu0*om)/k
Zall[nr+1] = Z *((Zall[nr] + Z*np.tanh(1j*k*hi))/(Z + Zall[nr]*np.tanh(1j*k*hi)))
#pdb.set_trace()
Z1d[nrFr] = Zall[-1]
return Z1d
simpegMT/Analytics/MT1Dsolutions.py
-43
View File
@@ -1,43 +0,0 @@
import numpy as np, SimPEG as simpeg
from MT1Danalytic import getEHfields
from scipy.constants import mu_0
def get1DEfields(m1d,sigma,freq,sourceAmp=1.0):
"""Function to get 1D electrical fields"""
# Get the gradient
G = m1d.nodalGrad
# Mass matrices
# Magnetic permeability
Mmu = simpeg.Utils.sdiag(m1d.vol*(1.0/mu_0))
# Conductivity
Msig = m1d.getFaceInnerProduct(sigma)
# Set up the solution matrix
A = G.T*Mmu*G - 1j*2.*np.pi*freq*Msig
# Define the inner part of the solution matrix
Aii = A[1:-1,1:-1]
# Define the outer part of the solution matrix
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
## Note: need to use conjugate of the analytic solution. It is derived with e^iwt
bc = np.r_[Etot_low.conj(),sourceAmp]
# The right hand side
rhs = -Aio*bc
# Solve the system
Aii_inv = simpeg.Solver(Aii)
eii = Aii_inv*rhs
# Assign the boundary conditions
e = np.r_[bc[0],eii,bc[1]]
# Return the electrical fields
return e
if __name__ == '__main__':
hz = [(100.,18)]
M = simpeg.Mesh.TensorMesh([hz],'C')
sig = np.zeros(M.nC) + 1e-8
sig[M.vectorCCx<=0] = sigHalf
simpegMT/Analytics/__init__.py
-1
View File
@@ -1 +0,0 @@
from MT1Dsolutions import * # Add the names of the functions
simpegMT/Utils/MT1Danalytic.py
+77 -76
View File
@@ -3,98 +3,99 @@
import numpy as np, SimPEG as simpeg import numpy as np, SimPEG as simpeg
def getEHfields(m1d,sigma,freq,zd): def getEHfields(m1d,sigma,freq,zd):
'''Analytic solution for MT 1D layered earth. Returns E and H fields. '''Analytic solution for MT 1D layered earth. Returns E and H fields.
:param SimPEG.mesh, object m1d: Mesh object with the 1D spatial information. :param SimPEG.mesh, object m1d: Mesh object with the 1D spatial information.
:param numpy.array, vector sigma: Physical property of conductivity corresponding with the mesh. :param numpy.array, vector sigma: Physical property of conductivity corresponding with the mesh.
:param float, freq: Frequency to calculate data at. :param float, freq: Frequency to calculate data at.
:param numpy array, vector zd: location to calculate EH fields at :param numpy array, vector zd: location to calculate EH fields at
Assumes a halfspace with the same conductive as the last cell below. Assumes a halfspace with the same conductive as the last cell below.
''' '''
# Note add an error check for the mesh. # Note add an error check for the mesh and sigma are the same size.
# Need make the solution e^-iwt dependent
# Constants: Assume constant # Constants: Assume constant
mu = 4*np.pi*1e-7*np.ones((m1d.nC+1)) mu = 4*np.pi*1e-7*np.ones((m1d.nC+1))
eps = 8.85*1e-12*np.ones((m1d.nC+1)) eps = 8.85*1e-12*np.ones((m1d.nC+1))
# Angular freq # Angular freq
w = 2*np.pi*freq w = 2*np.pi*freq
# Add the halfspace value to the property # Add the halfspace value to the property
sig = np.concatenate((np.array([sigma[0]]),sigma)) sig = np.concatenate((np.array([sigma[0]]),sigma))
# Calculate the wave number # Calculate the wave number
k = np.sqrt(eps*mu*w**2-1j*mu*sig*w) k = np.sqrt(eps*mu*w**2-1j*mu*sig*w)
# Initiate the propagation matrix, in the order down up. # Initiate the propagation matrix, in the order down up.
UDp = np.zeros((2,m1d.nC+1),dtype=complex) UDp = np.zeros((2,m1d.nC+1),dtype=complex)
UDp[1,0] = 1 # Set the wave amplitude as 1 into the half-space at the bottom of the mesh UDp[1,0] = 1 # Set the wave amplitude as 1 into the half-space at the bottom of the mesh
# Loop over all the layers, starting at the bottom layer # Loop over all the layers, starting at the bottom layer
for lnr, h in enumerate(m1d.hx): # lnr-number of layer, h-thickness of the layer for lnr, h in enumerate(m1d.hx): # lnr-number of layer, h-thickness of the layer
# Calculate # Calculate
yp1 = k[lnr]/(w*mu[lnr]) # Admittance of the layer below the current layer yp1 = k[lnr]/(w*mu[lnr]) # Admittance of the layer below the current layer
zp = (w*mu[lnr+1])/k[lnr+1] # Impedance in the current layer zp = (w*mu[lnr+1])/k[lnr+1] # Impedance in the current layer
# Build the propagation matrix # Build the propagation matrix
# Convert fields to down/up going components in layer below current layer
Pj1 = np.array([[1,1],[yp1,-yp1]])
# Convert fields to down/up going components in current layer
Pjinv = 1./2*np.array([[1,zp],[1,-zp]])
# Propagate down and up components through the current layer
elamh = np.array([[np.exp(-1j*k[lnr+1]*h),0],[0,np.exp(1j*k[lnr+1]*h)]])
# The down and up component in current layer. # Convert fields to down/up going components in layer below current layer
UDp[:,lnr+1] = elamh.dot(Pjinv.dot(Pj1)).dot(UDp[:,lnr]) Pj1 = np.array([[1,1],[yp1,-yp1]])
# Convert fields to down/up going components in current layer
Pjinv = 1./2*np.array([[1,zp],[1,-zp]])
# Propagate down and up components through the current layer
elamh = np.array([[np.exp(-1j*k[lnr+1]*h),0],[0,np.exp(1j*k[lnr+1]*h)]])
# Calculate the fields # The down and up component in current layer.
Ed = np.zeros((zd.size,),dtype=complex) UDp[:,lnr+1] = elamh.dot(Pjinv.dot(Pj1)).dot(UDp[:,lnr])
Eu = np.zeros((zd.size,),dtype=complex)
Hd = np.zeros((zd.size,),dtype=complex)
Hu = np.zeros((zd.size,),dtype=complex)
# Loop over the layers and calculate the fields # Calculate the fields
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::]): Ed = np.zeros((zd.size,),dtype=complex)
dind = np.logical_and(dup >= zd, zd > dlow) Eu = np.zeros((zd.size,),dtype=complex)
Ed[dind] = Dp*np.exp(-1j*ki*(dup-zd[dind])) Hd = np.zeros((zd.size,),dtype=complex)
Eu[dind] = Up*np.exp(1j*ki*(dup-zd[dind])) Hu = np.zeros((zd.size,),dtype=complex)
Hd[dind] = (ki/(w*mui))*Dp*np.exp(-1j*ki*(dup-zd[dind]))
Hu[dind] = -(ki/(w*mui))*Up*np.exp(1j*ki*(dup-zd[dind]))
# Return return the fields # Loop over the layers and calculate the fields
return Ed, Eu, Hd, Hu 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]))
Eu[dind] = Up*np.exp(1j*ki*(dup-zd[dind]))
Hd[dind] = (ki/(w*mui))*Dp*np.exp(-1j*ki*(dup-zd[dind]))
Hu[dind] = -(ki/(w*mui))*Up*np.exp(1j*ki*(dup-zd[dind]))
# Return return the fields
return Ed, Eu, Hd, Hu
def getImpedance(m1d,sigma,freq): def getImpedance(m1d,sigma,freq):
"""Analytic solution for MT 1D layered earth. Returns the impedance at the surface. """Analytic solution for MT 1D layered earth. Returns the impedance at the surface.
:param SimPEG.mesh, object m1d: Mesh object with the 1D spatial information. :param SimPEG.mesh, object m1d: Mesh object with the 1D spatial information.
:param numpy.array, vector sigma: Physical property corresponding with the mesh. :param numpy.array, vector sigma: Physical property corresponding with the mesh.
:param numpy.array, vector freq: Frequencies to calculate data at. :param numpy.array, vector freq: Frequencies to calculate data at.
""" """
# Define constants # Define constants
mu0 = 4*np.pi*1e-7 mu0 = 4*np.pi*1e-7
eps0 = 8.85e-12 eps0 = 8.85e-12
# Initiate the impedances # Initiate the impedances
Z1d = np.empty(len(freq) , dtype='complex') Z1d = np.empty(len(freq) , dtype='complex')
h = m1d.hx #vectorNx[:-1] h = m1d.hx #vectorNx[:-1]
# Start the process # Start the process
for nrFr, fr in enumerate(freq): for nrFr, fr in enumerate(freq):
om = 2*np.pi*fr om = 2*np.pi*fr
Zall = np.empty(len(h)+1,dtype='complex') Zall = np.empty(len(h)+1,dtype='complex')
# Calculate the impedance for the bottom layer # Calculate the impedance for the bottom layer
Zall[0] = (mu0*om)/np.sqrt(mu0*eps0*(om)**2 - 1j*mu0*sigma[0]*om) Zall[0] = (mu0*om)/np.sqrt(mu0*eps0*(om)**2 - 1j*mu0*sigma[0]*om)
for nr,hi in enumerate(h): for nr,hi in enumerate(h):
# Calculate the wave number # Calculate the wave number
# print nr,sigma[nr] # print nr,sigma[nr]
k = np.sqrt(mu0*eps0*om**2 - 1j*mu0*sigma[nr]*om) k = np.sqrt(mu0*eps0*om**2 - 1j*mu0*sigma[nr]*om)
Z = (mu0*om)/k Z = (mu0*om)/k
Zall[nr+1] = Z *((Zall[nr] + Z*np.tanh(1j*k*hi))/(Z + Zall[nr]*np.tanh(1j*k*hi)))
#pdb.set_trace() Zall[nr+1] = Z *((Zall[nr] + Z*np.tanh(1j*k*hi))/(Z + Zall[nr]*np.tanh(1j*k*hi)))
Z1d[nrFr] = Zall[-1]
return Z1d #pdb.set_trace()
Z1d[nrFr] = Zall[-1]
return Z1d
simpegMT/Utils/MT1Dsolutions.py
+33 -32
View File
@@ -3,40 +3,41 @@ from MT1Danalytic import getEHfields
from scipy.constants import mu_0 from scipy.constants import mu_0
def get1DEfields(m1d,sigma,freq,sourceAmp=1.0): def get1DEfields(m1d,sigma,freq,sourceAmp=1.0):
"""Function to get 1D electrical fields""" """Function to get 1D electrical fields"""
# Get the gradient
G = m1d.nodalGrad
# Mass matrices
# Magnetic permeability
Mmu = simpeg.Utils.sdiag(m1d.vol*(1.0/mu_0))
# Conductivity
Msig = m1d.getFaceInnerProduct(sigma)
# Set up the solution matrix
A = G.T*Mmu*G - 1j*2.*np.pi*freq*Msig
# Define the inner part of the solution matrix
Aii = A[1:-1,1:-1]
# Define the outer part of the solution matrix
Aio = A[1:-1,[0,-1]]
# Set the boundary conditions # Get the gradient
Ed_low, Eu_low, Hd_low, Hu_low = getEHfields(m1d,sigma,freq,np.array([m1d.vectorNx[0]])) G = m1d.nodalGrad
Etot_low = Ed_low + Eu_low # Mass matrices
bc = np.r_[Etot_low,sourceAmp] # Magnetic permeability
# The right hand side Mmu = simpeg.Utils.sdiag(m1d.vol*(1.0/mu_0))
rhs = -Aio*bc # Conductivity
# Solve the system Msig = m1d.getFaceInnerProduct(sigma)
Aii_inv = simpeg.Solver(Aii) # Set up the solution matrix
eii = Aii_inv*rhs A = G.T*Mmu*G - 1j*2.*np.pi*freq*Msig
# Assign the boundary conditions # Define the inner part of the solution matrix
e = np.r_[bc[0],eii,bc[1]] Aii = A[1:-1,1:-1]
# Return the electrical fields # Define the outer part of the solution matrix
return e 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
## Note: need to use conjugate of the analytic solution. It is derived with e^iwt
bc = np.r_[Etot_low.conj(),sourceAmp]
# The right hand side
rhs = -Aio*bc
# Solve the system
Aii_inv = simpeg.Solver(Aii)
eii = Aii_inv*rhs
# Assign the boundary conditions
e = np.r_[bc[0],eii,bc[1]]
# Return the electrical fields
return e
if __name__ == '__main__': if __name__ == '__main__':
hz = [(100.,18)] hz = [(100.,18)]
M = simpeg.Mesh.TensorMesh([hz],'C') M = simpeg.Mesh.TensorMesh([hz],'C')
sig = np.zeros(M.nC) + 1e-8 sig = np.zeros(M.nC) + 1e-8
sig[M.vectorCCx<=0] = sigHalf sig[M.vectorCCx<=0] = sigHalf
simpegMT/Utils/__init__.py
+1 -1
View File
@@ -1,2 +1,2 @@
from MT1Dsolutions import * # Add the names of the functions from MT1Dsolutions import * # Add the names of the functions
from MT1Danalytic import * from MT1Danalytic import *