diff --git a/SimPEG/MT/BaseMT.py b/SimPEG/MT/BaseMT.py index ddcff33f..d6440bdb 100644 --- a/SimPEG/MT/BaseMT.py +++ b/SimPEG/MT/BaseMT.py @@ -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: diff --git a/SimPEG/MT/FieldsMT.py b/SimPEG/MT/FieldsMT.py index 87ea109d..86ac7343 100644 --- a/SimPEG/MT/FieldsMT.py +++ b/SimPEG/MT/FieldsMT.py @@ -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. diff --git a/SimPEG/MT/Problem1D/Probs.py b/SimPEG/MT/Problem1D/Probs.py index 54d7c91b..6232b2e9 100644 --- a/SimPEG/MT/Problem1D/Probs.py +++ b/SimPEG/MT/Problem1D/Probs.py @@ -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. diff --git a/SimPEG/MT/Problem3D/Probs.py b/SimPEG/MT/Problem3D/Probs.py index 768a2d61..512362b8 100644 --- a/SimPEG/MT/Problem3D/Probs.py +++ b/SimPEG/MT/Problem3D/Probs.py @@ -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 - diff --git a/SimPEG/MT/Problem3D/__init__.py b/SimPEG/MT/Problem3D/__init__.py index 58c374ca..27e761a9 100644 --- a/SimPEG/MT/Problem3D/__init__.py +++ b/SimPEG/MT/Problem3D/__init__.py @@ -1 +1 @@ -from Probs import eForm_ps, eForm_Tp \ No newline at end of file +from Probs import eForm_ps \ No newline at end of file diff --git a/SimPEG/MT/Sources/__init__.py b/SimPEG/MT/Sources/__init__.py deleted file mode 100644 index 7b4a53bc..00000000 --- a/SimPEG/MT/Sources/__init__.py +++ /dev/null @@ -1 +0,0 @@ -from backgroundModelSources import * \ No newline at end of file diff --git a/SimPEG/MT/SrcMT.py b/SimPEG/MT/SrcMT.py index 148ecab4..70698c88 100644 --- a/SimPEG/MT/SrcMT.py +++ b/SimPEG/MT/SrcMT.py @@ -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).' diff --git a/SimPEG/MT/SurveyMT.py b/SimPEG/MT/SurveyMT.py index 83c1f81a..f7d4b4e8 100644 --- a/SimPEG/MT/SurveyMT.py +++ b/SimPEG/MT/SurveyMT.py @@ -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] diff --git a/SimPEG/MT/Sources/backgroundModelSources.py b/SimPEG/MT/Utils/sourceUtils.py similarity index 100% rename from SimPEG/MT/Sources/backgroundModelSources.py rename to SimPEG/MT/Utils/sourceUtils.py diff --git a/SimPEG/MT/__init__.py b/SimPEG/MT/__init__.py index dfbd9a15..b9f5536b 100644 --- a/SimPEG/MT/__init__.py +++ b/SimPEG/MT/__init__.py @@ -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