wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-06-27 19:48:52 +08:00
Code Issues Packages Projects Releases Wiki Activity

Updated doc strings, moved functions around and cleaned unused text.

Browse Source
This commit is contained in:
GudniRos
2016-02-04 00:23:07 -08:00
parent 01bf46c8bf
commit 345e39ed43
10 changed files with 103 additions and 176 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/MT/BaseMT.py
+28 -14
View File
@@ -5,6 +5,9 @@ from FieldsMT import BaseMTFields
class BaseMTProblem(BaseFDEMProblem):
"""
Base class for all Natural source problems.
"""
def __init__(self, mesh, **kwargs):
BaseFDEMProblem.__init__(self, mesh, **kwargs)
@@ -23,20 +26,21 @@ class BaseMTProblem(BaseFDEMProblem):
# Use the forward and devs from BaseFDEMProblem
# Might need to add more stuff here.
def Jvec(self, m, v, u=None):
## NEED to clean up the Jvec and Jtvec to use Zero and Identities for None components.
def Jvec(self, m, v, f=None):
"""
Function to calculate the data sensitivities dD/dm times a vector.
:param numpy.ndarray (nC, 1) - conductive model
:param numpy.ndarray (nC, 1) - random vector
:param numpy.ndarray m (nC, 1) - conductive model
:param numpy.ndarray v (nC, 1) - random vector
:param MTfields object (optional) - MT fields object, if not given it is calculated
:rtype: MTdata object
:return: Data sensitivities wrt m
"""
# Calculate the fields
if u is None:
u = self.fields(m)
if f is None:
f = self.fields(m)
# Set current model
self.curModel = m
# Initiate the Jv object
@@ -52,9 +56,9 @@ class BaseMTProblem(BaseFDEMProblem):
# We need fDeriv_m = df/du*du/dm + df/dm
# Construct du/dm, it requires a solve
# NOTE: need to account for the 2 polarizations in the derivatives.
u_src = u[src,:]
f_src = f[src,:]
# dA_dm and dRHS_dm should be of size nE,2, so that we can multiply by dA_duI. The 2 columns are each of the polarizations.
dA_dm = self.getADeriv_m(freq, u_src, v) # Size: nE,2 (u_px,u_py) in the columns.
dA_dm = self.getADeriv_m(freq, f_src, v) # Size: nE,2 (u_px,u_py) in the columns.
dRHS_dm = self.getRHSDeriv_m(freq, v) # Size: nE,2 (u_px,u_py) in the columns.
if dRHS_dm is None:
du_dm = dA_duI * ( -dA_dm )
@@ -64,15 +68,25 @@ class BaseMTProblem(BaseFDEMProblem):
for rx in src.rxList:
# Get the projection derivative
# v should be of size 2*nE (for 2 polarizations)
PDeriv_u = lambda t: rx.projectFieldsDeriv(src, self.mesh, u, t) # wrt u, we don't have have PDeriv wrt m
PDeriv_u = lambda t: rx.projectFieldsDeriv(src, self.mesh, f, t) # wrt u, we don't have have PDeriv wrt m
Jv[src, rx] = PDeriv_u(mkvc(du_dm))
dA_duI.clean()
# Return the vectorized sensitivities
return mkvc(Jv)
def Jtvec(self, m, v, u=None):
if u is None:
u = self.fields(m)
def Jtvec(self, m, v, f=None):
"""
Function to calculate the transpose of the data sensitivities (dD/dm)^T times a vector.
:param numpy.ndarray m (nC, 1) - conductive model
:param numpy.ndarray v (nD, 1) - vector
:param MTfields object f (optional) - MT fields object, if not given it is calculated
:rtype: MTdata object
:return: Data sensitivities wrt m
"""
if f is None:
f = self.fields(m)
self.curModel = m
@@ -89,15 +103,15 @@ class BaseMTProblem(BaseFDEMProblem):
for src in self.survey.getSrcByFreq(freq):
ftype = self._fieldType + 'Solution'
u_src = u[src, :]
f_src = f[src, :]
for rx in src.rxList:
# Get the adjoint projectFieldsDeriv
# PTv needs to be nE,
PTv = rx.projectFieldsDeriv(src, self.mesh, u, mkvc(v[src, rx],2), adjoint=True) # wrt u, need possibility wrt m
PTv = rx.projectFieldsDeriv(src, self.mesh, f, mkvc(v[src, rx],2), adjoint=True) # wrt u, need possibility wrt m
# Get the
dA_duIT = ATinv * PTv
dA_dmT = self.getADeriv_m(freq, u_src, mkvc(dA_duIT), adjoint=True)
dA_dmT = self.getADeriv_m(freq, f_src, mkvc(dA_duIT), adjoint=True)
dRHS_dmT = self.getRHSDeriv_m(freq, mkvc(dA_duIT), adjoint=True)
# Make du_dmT
if dRHS_dmT is None:
SimPEG/MT/FieldsMT.py
+4 -1
View File
@@ -122,7 +122,10 @@ class Fields1D_e(BaseMTFields):
class Fields3D_e(BaseMTFields):
"""
Fields storage for the 3D MT solution.
Fields storage for the 3D MT solution. Labels polarizations by px and py.
:param SimPEG object mesh: The solution mesh
:param SimPEG object survey: A survey object
"""
# Define the known the alias fields
# Assume that the solution of e on the E.
SimPEG/MT/Problem1D/Probs.py
+24 -1
View File
@@ -13,7 +13,21 @@ class eForm_psField(BaseMTProblem):
"""
A MT problem soving a e formulation and primary/secondary fields decomposion.
Solves the equation
By eliminating the magnetic flux density using
.. math ::
\mathbf{b} = \\frac{1}{i \omega}\\left(-\mathbf{C} \mathbf{e} \\right)
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}
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 ).
"""
# From FDEMproblem: Used to project the fields. Currently not used for MTproblem.
@@ -29,6 +43,10 @@ class eForm_psField(BaseMTProblem):
@property
def sigmaPrimary(self):
"""
A background model, use for the calculation of the primary fields.
"""
return self._sigmaPrimary
@sigmaPrimary.setter
def sigmaPrimary(self, val):
@@ -128,6 +146,11 @@ class eForm_psField(BaseMTProblem):
sys.stdout.flush()
return F
# Note this is not fully functional.
# Missing:
# Fields class corresponding to the fields
# Update Jvec and Jtvec to include all the derivatives components
# Other things ...
class eForm_TotalField(BaseMTProblem):
"""
A MT problem solving a e formulation and a Total bondary domain decompostion.
SimPEG/MT/Problem3D/Probs.py
+24 -151
View File
@@ -12,9 +12,20 @@ class eForm_ps(BaseMTProblem):
"""
A MT problem solving a e formulation and a primary/secondary fields decompostion.
Solves the equation:
By eliminating the magnetic flux density using
.. math ::
\mathbf{b} = \\frac{1}{i \omega}\\left(-\mathbf{C} \mathbf{e} \\right)
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}
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 as a 1D model).
"""
@@ -29,6 +40,10 @@ class eForm_ps(BaseMTProblem):
@property
def sigmaPrimary(self):
"""
A background model, use for the calculation of the primary fields.
"""
return self._sigmaPrimary
@sigmaPrimary.setter
def sigmaPrimary(self, val):
@@ -37,7 +52,7 @@ class eForm_ps(BaseMTProblem):
def getA(self, freq):
"""
Function to get the A matrix.
Function to get the A system.
:param float freq: Frequency
:rtype: scipy.sparse.csr_matrix
@@ -50,16 +65,10 @@ class eForm_ps(BaseMTProblem):
return C.T*Mmui*C + 1j*omega(freq)*Msig
def getADeriv_m(self, freq, u, v, adjoint=False):
"""
Calculate the derivative of A wrt m.
# Nee to account for both the polarizations
# dMe_dsig = (self.MeSigmaDeriv( u['e_pxSolution'] ) + self.MeSigmaDeriv( u['e_pySolution'] ))
# dMe_dsig = (self.MeSigmaDeriv( u['e_pxSolution'] + u['e_pySolution'] ))
# # dMe_dsig = self.MeSigmaDeriv( u )
# if adjoint:
# return 1j * omega(freq) * ( dMe_dsig.T * v ) # As in simpegEM
# return 1j * omega(freq) * ( dMe_dsig * v ) # As in simpegEM
"""
# This considers both polarizations and returns a nE,2 matrix for each polarization
if adjoint:
@@ -68,19 +77,12 @@ class eForm_ps(BaseMTProblem):
# Need a nE,2 matrix to be returned
dMe_dsigV = np.hstack(( mkvc(self.MeSigmaDeriv( u['e_pxSolution'] )*v,2), mkvc( self.MeSigmaDeriv(u['e_pySolution'] )*v,2) ))
return 1j * omega(freq) * dMe_dsigV
# Stacking them
# if adjoint:
# dMe_dsig = sp.vstack((self.MeSigmaDeriv( u['e_pxSolution'] ), self.MeSigmaDeriv( u['e_pxSolution'] ) )).T
# # self.MeSigmaDeriv(u['e_pySolution'] ).T*v,2) ))
# else:
# dMe_dsig = sp.vstack((self.MeSigmaDeriv( u['e_pxSolution'] ), self.MeSigmaDeriv( u['e_pxSolution'] ) ))
# return 1j * omega(freq) * dMe_dsig*v
def getRHS(self, freq):
"""
Function to return the right hand side for the system.
Function to return the right hand side for the system.
:param float freq: Frequency
:rtype: numpy.ndarray (nE, 2), numpy.ndarray (nE, 2)
:return: RHS for both polarizations, primary fields
@@ -117,6 +119,7 @@ class eForm_ps(BaseMTProblem):
sys.stdout.flush()
A = self.getA(freq)
rhs = self.getRHS(freq)
# Solve the system
Ainv = self.Solver(A, **self.solverOpts)
e_s = Ainv * rhs
@@ -126,140 +129,10 @@ class eForm_ps(BaseMTProblem):
F[Src, 'e_pxSolution'] = e_s[:,0]
F[Src, 'e_pySolution'] = e_s[:,1]
# 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]
if self.verbose:
print 'Ran for {:f} seconds'.format(time.time()-startTime)
sys.stdout.flush()
Ainv.clean()
return F
class eForm_Tp(BaseMTProblem):
"""
A MT problem solving a e formulation and a total/primary fields decompostion.
Solves the equation
"""
_fieldType = 'e'
_eqLocs = 'FE'
fieldsPair = Fields3D_e
# Set new properties
# Background model
@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 MeSigmaBack(self):
#TODO: hardcoded to sigma as the model
if getattr(self, '_MeSigmaBack', None) is None:
sigma = self.curModel
sigmaBG = self.backModel
self._MeSigmaBack = self.mesh.getEdgeInnerProduct(sigmaBG)
return self._MeSigmaBack
def __init__(self, mesh, **kwargs):
BaseMTProblem.__init__(self, mesh, **kwargs)
def getA(self, freq):
"""
Function to get the A matrix.
:param float freq: Frequency
:rtype: scipy.sparse.csr_matrix
:return: A
"""
mui = self.MfMui
sig = self.MeSigma
C = self.mesh.edgeCurl
return C.T*mui*C + 1j*omega(freq)*sig
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 getRHS(self, freq, backSigma):
"""
Function to return the right hand side for the system.
:param float freq: Frequency
:param numpy.ndarray (nC,) backSigma: Background conductivity model
:rtype: numpy.ndarray (nE, 2)
:return: one RHS for both polarizations
"""
# Get sources for the frequency
src = self.survey.getSources(freq)
# Make sure that there is 2 polarizations.
# assert len()
# Get the background electric fields
from SimPEG.MT.Sources import homo1DModelSource
eBG_bp = homo1DModelSource(self.mesh,freq,backSigma)
MeBack = self.MeSigmaBack
# Set up the A system
mui = self.MfMui
C = self.mesh.edgeCurl
Abg = C.T*mui*C + 1j*omega(freq)*MeBack
return Abg*eBG_bp, eBG_bp
def getRHSderiv(self, freq, backSigma, u, v, adjoint=False):
raise NotImplementedError('getRHSDeriv not implemented yet')
return None
def fields(self, m, m_back):
'''
Function to calculate all the fields for the model m.
:param np.ndarray (nC,) m: Conductivity model
:param np.ndarray (nC,) m_back: Background conductivity model
'''
self.curModel = m
self.backModel = m_back
# RHS, CalcFields = self.getRHS(freq,m_back), self.calcFields
F = Fields3D_e(self.mesh, self.survey)
for freq in self.survey.freqs:
if self.verbose:
startTime = time.time()
print 'Starting work for {:.3e}'.format(freq)
sys.stdout.flush()
A = self.getA(freq)
rhs, e_p = self.getRHS(freq,m_back)
Ainv = self.Solver(A, **self.solverOpts)
e_s = Ainv * rhs
e = e_s
# Store the fields
Src = self.survey.getSources(freq)
# Store the fieldss
F[Src, 'e_px'] = e[:,0]
F[Src, 'e_py'] = e[:,1]
# 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]
if self.verbose:
print 'Ran for {:f} seconds'.format(time.time()-startTime)
sys.stdout.flush()
return F
SimPEG/MT/Problem3D/__init__.py
+1 -1
View File
@@ -1 +1 @@
from Probs import eForm_ps, eForm_Tp
from Probs import eForm_ps
SimPEG/MT/Sources/__init__.py
-1
View File
@@ -1 +0,0 @@
from backgroundModelSources import *
SimPEG/MT/SrcMT.py
+6 -4
View File
@@ -3,7 +3,7 @@ from SimPEG.EM.FDEM.SrcFDEM import BaseSrc as FDEMBaseSrc
from SimPEG.EM.Utils import omega
from scipy.constants import mu_0
from numpy.lib import recfunctions as recFunc
from Sources import homo1DModelSource
from Utils.sourceUtils import homo1DModelSource
from Utils import rec2ndarr
import sys
@@ -31,7 +31,9 @@ class BaseMTSrc(FDEMBaseSrc):
# 1D sources
class polxy_1DhomotD(BaseMTSrc):
"""
MT source for both polarizations (x and y) for the total Domain. It calculates fields calculated based on conditions on the boundary of the domain.
MT source for both polarizations (x and y) for the total Domain.
It calculates fields calculated based on conditions on the boundary of the domain.
"""
def __init__(self, rxList, freq):
BaseMTSrc.__init__(self, rxList, freq)
@@ -42,8 +44,8 @@ class polxy_1DhomotD(BaseMTSrc):
# Need to implement such that it works for all dims.
class polxy_1Dprimary(BaseMTSrc):
"""
MT source for both polarizations (x and y) given a 1D primary models. It assigns fields calculated from the 1D model
as fields in the full space of the problem.
MT source for both polarizations (x and y) given a 1D primary models.
It assigns fields calculated from the 1D model as fields in the full space of the problem.
"""
def __init__(self, rxList, freq):
# assert mkvc(self.mesh.hz.shape,1) == mkvc(sigma1d.shape,1),'The number of values in the 1D background model does not match the number of vertical cells (hz).'
SimPEG/MT/SurveyMT.py
+16 -2
View File
@@ -3,7 +3,6 @@ from SimPEG.EM.FDEM.SrcFDEM import BaseSrc as FDEMBaseSrc
from SimPEG.EM.Utils import omega
from scipy.constants import mu_0
from numpy.lib import recfunctions as recFunc
from Sources import homo1DModelSource
from Utils import rec2ndarr
import SrcMT
import sys
@@ -12,6 +11,15 @@ import sys
### Receivers ###
#################
class Rx(SimPEGsurvey.BaseRx):
"""
Class that defines natural source receivers.
See knownRxTypes for types of allowed receivers.
:param ndArray locs: Locations of the receivers
:param str rxType: The type of receiver
"""
knownRxTypes = {
# 3D impedance
@@ -81,9 +89,14 @@ class Rx(SimPEGsurvey.BaseRx):
def projectFields(self, src, mesh, f):
'''
Project the fields and return the correct data.
Project the fields to natural source data.
:param SrcMT src: The source of the fields to project
:param SimPEG.Mesh mesh:
:param FieldsMT f: Natural source fields object to project
'''
## NOTE: Assumes that e is on t
if self.projType is 'Z1D':
Pex = mesh.getInterpolationMat(self.locs[:,-1],'Fx')
Pbx = mesh.getInterpolationMat(self.locs[:,-1],'Ex')
@@ -94,6 +107,7 @@ class Rx(SimPEGsurvey.BaseRx):
f_part_complex = -ex/bx
# elif self.projType is 'Z2D':
elif self.projType is 'Z3D':
## NOTE: Assumes that e is on edges and b on the faces. Need to generalize that or use a prop of fields to determine that.
if self.locs.ndim == 3:
eFLocs = self.locs[:,:,0]
bFLocs = self.locs[:,:,1]
SimPEG/MT/Sources/backgroundModelSources.py → SimPEG/MT/Utils/sourceUtils.py
View File
SimPEG/MT/__init__.py
-1
View File
@@ -1,5 +1,4 @@
import Utils
import Sources
from SurveyMT import Rx, Survey, Data
from FieldsMT import Fields1D_e, Fields3D_e
import Problem1D, Problem2D, Problem3D