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

Fixed comments from Lindsey and Rowan

Browse Source
This commit is contained in:
GudniRos
2016-02-04 22:38:19 -08:00
parent b8bd011662
commit 5705fabc5a
2 changed files with 57 additions and 34 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/MT/Problem1D/Probs.py
+51 -25
View File
@@ -23,7 +23,7 @@ class eForm_psField(BaseMTProblem):
we can write Maxwell's equations as a second order system in \\\(\\\mathbf{e}\\\) only:
.. math ::
\\left(\mathbf{C}^T \mathbf{M^f_{\mu^{-1}}} \mathbf{C} + i \omega \mathbf{M^e_\sigma}] \mathbf{e}_{s} =& i \omega \mathbf{M^e_{\delta \sigma}} \mathbf{e}_{p}
\\left(\mathbf{C}^T \mathbf{M^e_{\mu^{-1}}} \mathbf{C} + i \omega \mathbf{M^f_\sigma}] \mathbf{e}_{s} =& i \omega \mathbf{M^f_{\delta \sigma}} \mathbf{e}_{p}
which we solve for \\\(\\\mathbf{e_s}\\\). The total field \\\mathbf{e}\\ = \\\mathbf{e_p}\\ + \\\mathbf{e_s}\\.
The primary field is estimated from a background model (commonly half space ).
@@ -40,6 +40,23 @@ class eForm_psField(BaseMTProblem):
BaseMTProblem.__init__(self, mesh, **kwargs)
self.fieldsPair = Fields1D_e
# self._sigmaPrimary = sigmaPrimary
@property
def MeMui(self):
"""
Edge inner product matrix
"""
if getattr(self, '_MeMui', None) is None:
self._MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
return self._MeMui
@property
def MfSigma(self):
"""
Edge inner product matrix
"""
if getattr(self, '_MfSigma', None) is None:
self._MfSigma = self.mesh.getFaceInnerProduct(self.curModel.sigma)
return self._MfSigma
@property
def sigmaPrimary(self):
@@ -48,6 +65,7 @@ class eForm_psField(BaseMTProblem):
"""
return self._sigmaPrimary
@sigmaPrimary.setter
def sigmaPrimary(self, val):
# Note: TODO add logic for val, make sure it is the correct size.
@@ -62,16 +80,14 @@ class eForm_psField(BaseMTProblem):
:return: A
"""
Mmui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
Msig = self.mesh.getFaceInnerProduct(self.curModel.sigma)
# Note: need to use the code above since in the 1D problem I want
# e to live on Faces(nodes) and h on edges(cells). Might need to rethink this
# Possible that _fieldType and _eqLocs can fix this
# Mmui = self.MfMui
# Msig = self.MeSigma
MeMui = self.MfMui
MfSigma = self.MfSigma
C = self.mesh.nodalGrad
# Make A
A = C.T*Mmui*C + 1j*omega(freq)*Msig
A = C.T*MeMui*C + 1j*omega(freq)*MfSigma
# Either return full or only the inner part of A
return A
@@ -81,15 +97,15 @@ class eForm_psField(BaseMTProblem):
"""
dsig_dm = self.curModel.sigmaDeriv
MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
MeMui = self.MeMui
#
u_src = u['e_1dSolution']
dMf_dsig = self.mesh.getFaceInnerProductDeriv(self.curModel.sigma)(u_src) * self.curModel.sigmaDeriv
dMfSigma_dm = self.mesh.getFaceInnerProductDeriv(self.curModel.sigma)(u_src) * self.curModel.sigmaDeriv
if adjoint:
return 1j * omega(freq) * ( dMf_dsig.T * v )
return 1j * omega(freq) * ( dMfSigma_dm.T * v )
# Note: output has to be nN/nF, not nC/nE.
# v should be nC
return 1j * omega(freq) * ( dMf_dsig * v )
return 1j * omega(freq) * ( dMfSigma_dm * v )
def getRHS(self, freq):
"""
@@ -169,6 +185,23 @@ class eForm_TotalField(BaseMTProblem):
def __init__(self, mesh, **kwargs):
BaseMTProblem.__init__(self, mesh, **kwargs)
@property
def MeMui(self):
"""
Edge inner product matrix
"""
if getattr(self, '_MeMui', None) is None:
self._MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
return self._MeMui
@property
def MfSigma(self):
"""
Edge inner product matrix
"""
if getattr(self, '_MfSigma', None) is None:
self._MfSigma = self.mesh.getFaceInnerProduct(self.curModel.sigma)
return self._MfSigma
def getA(self, freq, full=False):
"""
@@ -180,31 +213,24 @@ class eForm_TotalField(BaseMTProblem):
:return: A
"""
Mmui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
Msig = self.mesh.getFaceInnerProduct(self.curModel.sigma)
MeMui = self.MeMui
MfSigma = self.MfSigma
# Note: need to use the code above since in the 1D problem I want
# e to live on Faces(nodes) and h on edges(cells). Might need to rethink this
# Possible that _fieldType and _eqLocs can fix this
# Mmui = self.MfMui
# Msig = self.MeSigma
# MeMui = self.MfMui
# MfSigma = self.MfSigma
C = self.mesh.nodalGrad
# Make A
A = C.T*Mmui*C + 1j*omega(freq)*Msig
A = C.T*MeMui*C + 1j*omega(freq)*MfSigma
# Either return full or only the inner part of A
if full:
return A
else:
return A[1:-1,1:-1]
def getADeriv(self, freq, u, v, adjoint=False):
sig = self.curTModel
dsig_dm = self.curTModelDeriv
dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, v=u)
if adjoint:
return 1j * omega(freq) * ( dsig_dm.T * ( dMe_dsig.T * v ) )
return 1j * omega(freq) * ( dMe_dsig * ( dsig_dm * v ) )
def getADeriv_m(self, freq, u, v, adjoint=False):
raise NotImplementedError('getADeriv is not implemented')
def getRHS(self, freq):
"""
@@ -230,7 +256,7 @@ class eForm_TotalField(BaseMTProblem):
return -Aio*eBC, eBC
def getRHSderiv(self, freq, backSigma, u, v, adjoint=False):
def getRHSderiv_m(self, freq, backSigma, u, v, adjoint=False):
raise NotImplementedError('getRHSDeriv not implemented yet')
return None
SimPEG/MT/Utils/MT1Danalytic.py
+6 -9
View File
@@ -1,6 +1,7 @@
# Analytic solution of EM fields due to a plane wave
import numpy as np, SimPEG as simpeg
from scipy.constants import mu_0, epsilon_0 as eps_0
def getEHfields(m1d,sigma,freq,zd,scaleUD=True):
'''Analytic solution for MT 1D layered earth. Returns E and H fields.
@@ -17,8 +18,8 @@ def getEHfields(m1d,sigma,freq,zd,scaleUD=True):
# Note add an error check for the mesh and sigma are the same size.
# Constants: Assume constant
mu = 4*np.pi*1e-7*np.ones((m1d.nC+1))
eps = 8.85*1e-12*np.ones((m1d.nC+1))
mu = mu_0*np.ones((m1d.nC+1))
eps = eps_0*np.ones((m1d.nC+1))
# Angular freq
w = 2*np.pi*freq
# Add the halfspace value to the property
@@ -83,10 +84,6 @@ def getImpedance(m1d,sigma,freq):
"""
# 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]
@@ -95,13 +92,13 @@ def getImpedance(m1d,sigma,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)
Zall[0] = (mu_0*om)/np.sqrt(mu_0*eps_0*(om)**2 - 1j*mu_0*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
k = np.sqrt(mu_0*eps_0*om**2 - 1j*mu_0*sigma[nr]*om)
Z = (mu_0*om)/k
Zall[nr+1] = Z *((Zall[nr] + Z*np.tanh(1j*k*hi))/(Z + Zall[nr]*np.tanh(1j*k*hi)))