From 6165e619edc8cfdefcfc3048603be7887dcf81f0 Mon Sep 17 00:00:00 2001 From: GudniRos Date: Tue, 10 May 2016 16:07:45 -0700 Subject: [PATCH] Moving Problem1D/3D folders to NSEM file Looking at making Jvec more general. --- SimPEG/NSEM/FieldsNSEM.py | 130 ++++++++- SimPEG/NSEM/NSEM.py | 438 +++++++++++++++++++++++++++++- SimPEG/NSEM/Problem1D/Probs.py | 291 -------------------- SimPEG/NSEM/Problem1D/__init__.py | 1 - SimPEG/NSEM/Problem2D/Probs.py | 0 SimPEG/NSEM/Problem2D/__init__.py | 1 - SimPEG/NSEM/Problem3D/Probs.py | 138 ---------- SimPEG/NSEM/Problem3D/__init__.py | 1 - 8 files changed, 559 insertions(+), 441 deletions(-) delete mode 100644 SimPEG/NSEM/Problem1D/Probs.py delete mode 100644 SimPEG/NSEM/Problem1D/__init__.py delete mode 100644 SimPEG/NSEM/Problem2D/Probs.py delete mode 100644 SimPEG/NSEM/Problem2D/__init__.py delete mode 100644 SimPEG/NSEM/Problem3D/Probs.py delete mode 100644 SimPEG/NSEM/Problem3D/__init__.py diff --git a/SimPEG/NSEM/FieldsNSEM.py b/SimPEG/NSEM/FieldsNSEM.py index 72bf3f6b..a13a0e0f 100644 --- a/SimPEG/NSEM/FieldsNSEM.py +++ b/SimPEG/NSEM/FieldsNSEM.py @@ -4,6 +4,7 @@ import sys from numpy.lib import recfunctions as recFunc from SimPEG.EM.Utils import omega + ############## ### Fields ### ############## @@ -12,8 +13,10 @@ class BaseNSEMFields(Problem.Fields): knownFields = {} dtype = complex - -class Fields1D_e(BaseNSEMFields): +########### +# 1D Fields +########### +class Fields1D_ePrimSec(BaseNSEMFields): """ Fields storage for the 1D NSEM solution. """ @@ -44,6 +47,118 @@ class Fields1D_e(BaseNSEMFields): def _e(self, eSolution, srcList): return self._ePrimary(eSolution,srcList) + self._eSecondary(eSolution,srcList) + def _eDeriv_u(self, src, du_dm_v, adjoint = False): + + + return Utils.Identity()*du_dm_v + + def _eDeriv_m(self, src, v, adjoint = False): + # assuming primary does not depend on the model + return Utils.Zero() + + def _bPrimary(self, eSolution, srcList): + bPrimary = np.zeros([self.survey.mesh.nE,eSolution.shape[1]], dtype = complex) + for i, src in enumerate(srcList): + bp = src.bPrimary(self.survey.prob) + if bp is not None: + bPrimary[:,i] += bp[:,-1] + return bPrimary + + def _bSecondary(self, eSolution, srcList): + C = self.mesh.nodalGrad + b = (C * eSolution) + for i, src in enumerate(srcList): + b[:,i] *= - 1./(1j*omega(src.freq)) + # There is no magnetic source in the MT problem + # S_m, _ = src.eval(self.survey.prob) + # if S_m is not None: + # b[:,i] += 1./(1j*omega(src.freq)) * S_m + return b + + def _b(self, eSolution, srcList): + return self._bPrimary(eSolution, srcList) + self._bSecondary(eSolution, srcList) + + def _bSecondaryDeriv_u(self, src, v, adjoint = False): + C = self.mesh.nodalGrad + if adjoint: + return - 1./(1j*omega(src.freq)) * (C.T * v) + return - 1./(1j*omega(src.freq)) * (C * v) + + def _bSecondaryDeriv_m(self, src, v, adjoint = False): + # Doesn't depend on m + # _, S_eDeriv = src.evalDeriv(self.survey.prob, adjoint) + # S_eDeriv = S_eDeriv(v) + # if S_eDeriv is not None: + # return 1./(1j * omega(src.freq)) * S_eDeriv + return None + + def _bDeriv_u(self, src, v, adjoint=False): + # Primary does not depend on u + return self._bSecondaryDeriv_u(src, v, adjoint) + + def _bDeriv_m(self, src, v, adjoint=False): + # Assuming the primary does not depend on the model + return self._bSecondaryDeriv_m(src, v, adjoint) + + def _fDeriv_u(self, src, v, adjoint=False): + """ + Derivative of the fields object wrt u. + + :param NSEMsrc src: NSEM source + :param numpy.ndarray v: random vector of f_sol.size + This function stacks the fields derivatives appropriately + + return a vector of size (nreEle+nrbEle) + """ + + de_du = v #Utils.spdiag(np.ones((self.nF,))) + db_du = self._bDeriv_u(src, v, adjoint) + # Return the stack + # This doesn't work... + return np.vstack((de_du,db_du)) + + def _fDeriv_m(self, src, v, adjoint=False): + """ + Derivative of the fields object wrt m. + + This function stacks the fields derivatives appropriately + """ + return None + + +class Fields1D_eTotal(BaseNSEMFields): + """ + Fields storage for the 1D NSEM solution solved with for a total domain formulation. + + Used in conjuction with Problem1D_eTotal. + """ + knownFields = {'e_1dSolution':'F'} + aliasFields = { + 'e_1d' : ['e_1dSolution','F','_e'], + 'e_1dPrimary' : ['e_1dSolution','F','_ePrimary'], + 'e_1dSecondary' : ['e_1dSolution','F','_eSecondary'], + 'b_1d' : ['e_1dSolution','E','_b'], + 'b_1dPrimary' : ['e_1dSolution','E','_bPrimary'], + 'b_1dSecondary' : ['e_1dSolution','E','_bSecondary'] + } + + def __init__(self,mesh,survey,**kwargs): + BaseNSEMFields.__init__(self,mesh,survey,**kwargs) + + def _ePrimary(self, eSolution, srcList): + ePrimary = np.zeros_like(eSolution) + for i, src in enumerate(srcList): + ep = src.ePrimary(self.survey.prob) + if ep is not None: + ePrimary[:,i] = ep[:,-1] + return ePrimary + + def _eSecondary(self, eSolution, srcList): + return eSolution + + def _e(self, eSolution, srcList): + return self._ePrimary(eSolution,srcList) + self._eSecondary(eSolution,srcList) + def _eDeriv_u(self, src, v, adjoint = False): return v @@ -120,7 +235,16 @@ class Fields1D_e(BaseNSEMFields): """ return None -class Fields3D_e(BaseNSEMFields): + +########### +# 2D Fields +########### + + +########### +# 3D Fields +########### +class Fields3D_ePrimSec(BaseNSEMFields): """ Fields storage for the 3D NSEM solution. Labels polarizations by px and py. diff --git a/SimPEG/NSEM/NSEM.py b/SimPEG/NSEM/NSEM.py index e9915f5d..d0bd8b77 100644 --- a/SimPEG/NSEM/NSEM.py +++ b/SimPEG/NSEM/NSEM.py @@ -1,8 +1,9 @@ +from SimPEG.EM.Utils import omega, mu_0 from SimPEG import SolverLU as SimpegSolver, PropMaps, Utils, mkvc, sp, np from SimPEG.EM.FDEM.FDEM import BaseFDEMProblem from SurveyNSEM import Survey, Data -from FieldsNSEM import BaseNSEMFields - +from FieldsNSEM import Fields1D_ePrimSec, Fields3D_ePrimSec +from SimPEG.NSEM.Utils.MT1Danalytic import getEHfields class BaseNSEMProblem(BaseFDEMProblem): """ @@ -56,9 +57,9 @@ class BaseNSEMProblem(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. - f_src = f[src,:] + u_src = f[src,:] # u should be a vector by definition. Need to fix this... # 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, f_src, v) # Size: nE,2 (u_px,u_py) in the columns. + dA_dm = self.getADeriv_m(freq, u_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 ) @@ -80,7 +81,7 @@ class BaseNSEMProblem(BaseFDEMProblem): :param numpy.ndarray m (nC, 1) - conductive model :param numpy.ndarray v (nD, 1) - vector - :param NSEMfields object u (optional) - NSEM fields object, if not given it is calculated + :param NSEMfields object f (optional) - NSEM fields object, if not given it is calculated :rtype: NSEMdata object :return: Data sensitivities wrt m """ @@ -103,7 +104,7 @@ class BaseNSEMProblem(BaseFDEMProblem): for src in self.survey.getSrcByFreq(freq): ftype = self._fieldType + 'Solution' - f_src = f[src, :] + f_src = f[src, :] # Need to fix this... for rx in src.rxList: # Get the adjoint evalDeriv @@ -130,3 +131,428 @@ class BaseNSEMProblem(BaseFDEMProblem): # Clean the factorization, clear memory. ATinv.clean() return Jtv + +################################### +## 1D problems +################################### + +class Problem1D_ePrimSec(BaseNSEMProblem): + """ + A NSEM problem soving a e formulation and primary/secondary fields decomposion. + + 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^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 ). + + + """ + + # From FDEMproblem: Used to project the fields. Currently not used for NSEMproblem. + _solutionType = 'e_1dSolution' + _formulation = 'EF' + fieldsPair = Fields1D_ePrimSec + + # Initiate properties + _sigmaPrimary = None + + + def __init__(self, mesh, **kwargs): + BaseNSEMProblem.__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): + """ + A background model, use for the calculation of the primary fields. + + """ + return self._sigmaPrimary + + @sigmaPrimary.setter + def sigmaPrimary(self, val): + # Note: TODO add logic for val, make sure it is the correct size. + self._sigmaPrimary = val + + def getA(self, freq): + """ + Function to get the A matrix. + + :param float freq: Frequency + :rtype: scipy.sparse.csr_matrix + :return: A + """ + + # 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 + MeMui = self.MeMui + MfSigma = self.MfSigma + C = self.mesh.nodalGrad + # Make A + A = C.T*MeMui*C + 1j*omega(freq)*MfSigma + # Either return full or only the inner part of A + return A + + def getADeriv_m(self, freq, u, v, adjoint=False): + """ + The derivative of A wrt sigma + """ + + dsig_dm = self.curModel.sigmaDeriv + MeMui = self.MeMui + # + u_src = u['e_1dSolution'] + dMfSigma_dm = self.mesh.getFaceInnerProductDeriv(self.curModel.sigma)(u_src) * self.curModel.sigmaDeriv + if adjoint: + 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) * ( dMfSigma_dm * v ) + + def getRHS(self, freq): + """ + Function to return the right hand side for the system. + :param float freq: Frequency + :rtype: numpy.ndarray (nF, 1), numpy.ndarray (nF, 1) + :return: RHS for 1 polarizations, primary fields + """ + + # Get sources for the frequncy(polarizations) + Src = self.survey.getSrcByFreq(freq)[0] + S_e = Src.S_e(self) + return -1j * omega(freq) * S_e + + def getRHSDeriv_m(self, freq, v, adjoint=False): + """ + The derivative of the RHS wrt sigma + """ + + Src = self.survey.getSrcByFreq(freq)[0] + S_eDeriv = Src.S_eDeriv_m(self, v, adjoint) + return -1j * omega(freq) * S_eDeriv + + def fields(self, m): + ''' + Function to calculate all the fields for the model m. + + :param np.ndarray (nC,) m: Conductivity model + ''' + # Set the current model + self.curModel = m + # Make the fields object + F = self.fieldsPair(self.mesh, self.survey) + # Loop over the frequencies + 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 = self.getRHS(freq) + Ainv = self.Solver(A, **self.solverOpts) + e_s = Ainv * rhs + + # Store the fields + Src = self.survey.getSrcByFreq(freq)[0] + # NOTE: only store the e_solution(secondary), all other components calculated in the fields object + F[Src, 'e_1dSolution'] = e_s[:,-1] # Only storing the yx polarization as 1d + + # Note curl e = -iwb so b = -curl e /iw + # b = -( self.mesh.nodalGrad * e )/( 1j*omega(freq) ) + # F[Src, 'b_1d'] = b[:,1] + if self.verbose: + print 'Ran for {:f} seconds'.format(time.time()-startTime) + 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 Problem1D_eTotal(BaseNSEMProblem): + """ + A NSEM problem solving a e formulation and a Total bondary domain decompostion. + + Solves the equation: + + Math: + Have to do this... + Not implement correctly....... + """ + + # From FDEMproblem: Used to project the fields. Currently not used for NSEMproblem. + _solutionType = 'e_1dSolution' + _formulation = 'EF' + fieldsPair = Fields1D_eTotal + + def __init__(self, mesh, **kwargs): + BaseNSEMProblem.__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): + """ + Function to get the A matrix. + + :param float freq: Frequency + :param logic full: Return full A or the inner part + :rtype: scipy.sparse.csr_matrix + :return: A + """ + + 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 + # MeMui = self.MfMui + # MfSigma = self.MfSigma + C = self.mesh.nodalGrad + # Make A + 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_m(self, freq, u, v, adjoint=False): + raise NotImplementedError('getADeriv is not implemented') + + def getRHS(self, freq): + """ + 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 + """ + # Get sources for the frequency + # NOTE: Need to use the source information, doesn't really apply in 1D + src = self.survey.getSrcByFreq(freq) + # Get the full A + A = self.getA(freq,full=True) + # Define the outer part of the solution matrix + Aio = A[1:-1,[0,-1]] + Ed, Eu, Hd, Hu = getEHfields(self.mesh,self.curModel.sigma,freq,self.mesh.vectorNx) + Etot = (Ed + Eu) + sourceAmp = 1.0 + Etot = ((Etot/Etot[-1])*sourceAmp) # Scale the fields to be equal to sourceAmp at the top + ## Note: The analytic solution is derived with e^iwt + eBC = np.r_[Etot[0],Etot[-1]] + # The right hand side + + return -Aio*eBC, eBC + + def getRHSderiv_m(self, freq, backSigma, u, v, adjoint=False): + raise NotImplementedError('getRHSDeriv not implemented yet') + return None + + def fields(self, m): + ''' + 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 + # RHS, CalcFields = self.getRHS(freq,m_back), self.calcFields + + F = Fields1D_eTotal(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_o = self.getRHS(freq) + Ainv = self.Solver(A, **self.solverOpts) + e_i = Ainv * rhs + e = mkvc(np.r_[e_o[0], e_i, e_o[1]],2) + # Store the fields + Src = self.survey.getSrcByFreq(freq) + # NOTE: only store e fields + F[Src, 'e_1dSolution'] = e[:,0] + if self.verbose: + print 'Ran for {:f} seconds'.format(time.time()-startTime) + sys.stdout.flush() + return F + + +################################### +## 3D problems +################################### +class Problem3D_ePrimSec(BaseNSEMProblem): + """ + A NSEM problem solving a e formulation and a primary/secondary fields decompostion. + + 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). + + """ + + # From FDEMproblem: Used to project the fields. Currently not used for NSEMproblem. + _solutionType = [ 'e_pxSolution', 'e_pySolution'] # Forces order on the object + _formulation = 'EB' + fieldsPair = Fields3D_ePrimSec + + # Initiate properties + _sigmaPrimary = None + + def __init__(self, mesh, **kwargs): + BaseNSEMProblem.__init__(self, mesh, **kwargs) + + @property + def sigmaPrimary(self): + """ + A background model, use for the calculation of the primary fields. + + """ + return self._sigmaPrimary + @sigmaPrimary.setter + def sigmaPrimary(self, val): + # Note: TODO add logic for val, make sure it is the correct size. + self._sigmaPrimary = val + + def getA(self, freq): + """ + Function to get the A system. + + :param float freq: Frequency + :rtype: scipy.sparse.csr_matrix + :return: A + """ + Mmui = self.MfMui + Msig = self.MeSigma + C = self.mesh.edgeCurl + + 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. + + """ + # Fix u to be a matrix nE,2 + # This considers both polarizations and returns a nE,2 matrix for each polarization + if adjoint: + dMe_dsigV = sp.hstack(( self.MeSigmaDeriv( u['e_pxSolution'] ).T, self.MeSigmaDeriv(u['e_pySolution'] ).T ))*v + else: + # 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 + + + def getRHS(self, freq): + """ + 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 + """ + + # Get sources for the frequncy(polarizations) + Src = self.survey.getSrcByFreq(freq)[0] + S_e = Src.S_e(self) + return -1j * omega(freq) * S_e + + def getRHSDeriv_m(self, freq, v, adjoint=False): + """ + The derivative of the RHS with respect to sigma + """ + + Src = self.survey.getSrcByFreq(freq)[0] + S_eDeriv = Src.S_eDeriv_m(self, v, adjoint) + return -1j * omega(freq) * S_eDeriv + + def fields(self, m): + ''' + Function to calculate all the fields for the model m. + + :param np.ndarray (nC,) m: Conductivity model + ''' + # Set the current model + self.curModel = m + + F = self.fieldsPair(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 = self.getRHS(freq) + # Solve the system + Ainv = self.Solver(A, **self.solverOpts) + e_s = Ainv * rhs + + # Store the fields + Src = self.survey.getSrcByFreq(freq)[0] + # Store the fields + # Use self._solutionType + F[Src, 'e_pxSolution'] = e_s[:,0] + F[Src, 'e_pySolution'] = e_s[:,1] + # Note curl e = -iwb so b = -curl/iw + + if self.verbose: + print 'Ran for {:f} seconds'.format(time.time()-startTime) + sys.stdout.flush() + Ainv.clean() + return F diff --git a/SimPEG/NSEM/Problem1D/Probs.py b/SimPEG/NSEM/Problem1D/Probs.py deleted file mode 100644 index 1ac6771d..00000000 --- a/SimPEG/NSEM/Problem1D/Probs.py +++ /dev/null @@ -1,291 +0,0 @@ -from SimPEG.EM.Utils import omega -from SimPEG import mkvc -from scipy.constants import mu_0 -from SimPEG.NSEM.NSEM import BaseNSEMProblem -from SimPEG.NSEM.SurveyNSEM import Survey, Data -from SimPEG.NSEM.FieldsNSEM import Fields1D_e -from SimPEG.NSEM.Utils.MT1Danalytic import getEHfields -import numpy as np -import multiprocessing, sys, time - - -class eForm_psField(BaseNSEMProblem): - """ - A NSEM problem soving a e formulation and primary/secondary fields decomposion. - - 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^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 ). - - - """ - # From FDEMproblem: Used to project the fields. Currently not used for NSEMproblem. - _fieldType = 'e_1d' - _eqLocs = 'EF' - _sigmaPrimary = None - - - def __init__(self, mesh, **kwargs): - BaseNSEMProblem.__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): - """ - A background model, use for the calculation of the primary fields. - - """ - return self._sigmaPrimary - - @sigmaPrimary.setter - def sigmaPrimary(self, val): - # Note: TODO add logic for val, make sure it is the correct size. - self._sigmaPrimary = val - - def getA(self, freq): - """ - Function to get the A matrix. - - :param float freq: Frequency - :rtype: scipy.sparse.csr_matrix - :return: A - """ - - # 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 - MeMui = self.MeMui - MfSigma = self.MfSigma - C = self.mesh.nodalGrad - # Make A - A = C.T*MeMui*C + 1j*omega(freq)*MfSigma - # Either return full or only the inner part of A - return A - - def getADeriv_m(self, freq, u, v, adjoint=False): - """ - The derivative of A wrt sigma - """ - - dsig_dm = self.curModel.sigmaDeriv - MeMui = self.MeMui - # - u_src = u['e_1dSolution'] - dMfSigma_dm = self.mesh.getFaceInnerProductDeriv(self.curModel.sigma)(u_src) * self.curModel.sigmaDeriv - if adjoint: - 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) * ( dMfSigma_dm * v ) - - def getRHS(self, freq): - """ - Function to return the right hand side for the system. - :param float freq: Frequency - :rtype: numpy.ndarray (nF, 1), numpy.ndarray (nF, 1) - :return: RHS for 1 polarizations, primary fields - """ - - # Get sources for the frequncy(polarizations) - Src = self.survey.getSrcByFreq(freq)[0] - S_e = Src.S_e(self) - return -1j * omega(freq) * S_e - - def getRHSDeriv_m(self, freq, v, adjoint=False): - """ - The derivative of the RHS wrt sigma - """ - - Src = self.survey.getSrcByFreq(freq)[0] - S_eDeriv = Src.S_eDeriv_m(self, v, adjoint) - return -1j * omega(freq) * S_eDeriv - - def fields(self, m): - ''' - Function to calculate all the fields for the model m. - - :param np.ndarray (nC,) m: Conductivity model - ''' - # Set the current model - self.curModel = m - - F = Fields1D_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 = self.getRHS(freq) - Ainv = self.Solver(A, **self.solverOpts) - e_s = Ainv * rhs - - # Store the fields - Src = self.survey.getSrcByFreq(freq)[0] - # NOTE: only store the e_solution(secondary), all other components calculated in the fields object - F[Src, 'e_1dSolution'] = e_s[:,-1] # Only storing the yx polarization as 1d - - # Note curl e = -iwb so b = -curl e /iw - # b = -( self.mesh.nodalGrad * e )/( 1j*omega(freq) ) - # F[Src, 'b_1d'] = b[:,1] - if self.verbose: - print 'Ran for {:f} seconds'.format(time.time()-startTime) - 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(BaseNSEMProblem): - """ - A NSEM problem solving a e formulation and a Total bondary domain decompostion. - - Solves the equation: - - Math: - - - """ - - # From FDEMproblem: Used to project the fields. Currently not used for NSEMproblem. - _fieldType = 'e' - _eqLocs = 'EF' - - - def __init__(self, mesh, **kwargs): - BaseNSEMProblem.__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): - """ - Function to get the A matrix. - - :param float freq: Frequency - :param logic full: Return full A or the inner part - :rtype: scipy.sparse.csr_matrix - :return: A - """ - - 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 - # MeMui = self.MfMui - # MfSigma = self.MfSigma - C = self.mesh.nodalGrad - # Make A - 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_m(self, freq, u, v, adjoint=False): - raise NotImplementedError('getADeriv is not implemented') - - def getRHS(self, freq): - """ - 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 - """ - # Get sources for the frequency - # NOTE: Need to use the source information, doesn't really apply in 1D - src = self.survey.getSrcByFreq(freq) - # Get the full A - A = self.getA(freq,full=True) - # Define the outer part of the solution matrix - Aio = A[1:-1,[0,-1]] - Ed, Eu, Hd, Hu = getEHfields(self.mesh,self.curModel.sigma,freq,self.mesh.vectorNx) - Etot = (Ed + Eu) - sourceAmp = 1.0 - Etot = ((Etot/Etot[-1])*sourceAmp) # Scale the fields to be equal to sourceAmp at the top - ## Note: The analytic solution is derived with e^iwt - eBC = np.r_[Etot[0],Etot[-1]] - # The right hand side - - return -Aio*eBC, eBC - - def getRHSderiv_m(self, freq, backSigma, u, v, adjoint=False): - raise NotImplementedError('getRHSDeriv not implemented yet') - return None - - def fields(self, m): - ''' - 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 - # RHS, CalcFields = self.getRHS(freq,m_back), self.calcFields - - F = Fields1D_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_o = self.getRHS(freq) - Ainv = self.Solver(A, **self.solverOpts) - e_i = Ainv * rhs - e = mkvc(np.r_[e_o[0], e_i, e_o[1]],2) - # Store the fields - Src = self.survey.getSrcByFreq(freq) - # NOTE: only store e fields - F[Src, 'e_1dSolution'] = e[:,0] - if self.verbose: - print 'Ran for {:f} seconds'.format(time.time()-startTime) - sys.stdout.flush() - return F diff --git a/SimPEG/NSEM/Problem1D/__init__.py b/SimPEG/NSEM/Problem1D/__init__.py deleted file mode 100644 index b0a77250..00000000 --- a/SimPEG/NSEM/Problem1D/__init__.py +++ /dev/null @@ -1 +0,0 @@ -from Probs import eForm_TotalField, eForm_psField \ No newline at end of file diff --git a/SimPEG/NSEM/Problem2D/Probs.py b/SimPEG/NSEM/Problem2D/Probs.py deleted file mode 100644 index e69de29b..00000000 diff --git a/SimPEG/NSEM/Problem2D/__init__.py b/SimPEG/NSEM/Problem2D/__init__.py deleted file mode 100644 index fc80254b..00000000 --- a/SimPEG/NSEM/Problem2D/__init__.py +++ /dev/null @@ -1 +0,0 @@ -pass \ No newline at end of file diff --git a/SimPEG/NSEM/Problem3D/Probs.py b/SimPEG/NSEM/Problem3D/Probs.py deleted file mode 100644 index cb217abe..00000000 --- a/SimPEG/NSEM/Problem3D/Probs.py +++ /dev/null @@ -1,138 +0,0 @@ -from SimPEG import Survey, Problem, Utils, Models, np, sp, mkvc, SolverLU as SimpegSolver -from SimPEG.EM.Utils import omega -from scipy.constants import mu_0 -from SimPEG.NSEM.NSEM import BaseNSEMProblem -from SimPEG.NSEM.SurveyNSEM import Survey, Data -from SimPEG.NSEM.FieldsNSEM import Fields3D_e -import multiprocessing, sys, time - - - -class eForm_ps(BaseNSEMProblem): - """ - A NSEM problem solving a e formulation and a primary/secondary fields decompostion. - - 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). - - """ - - # From FDEMproblem: Used to project the fields. Currently not used for NSEMproblem. - _fieldType = 'e' - _eqLocs = 'FE' - fieldsPair = Fields3D_e - _sigmaPrimary = None - - def __init__(self, mesh, **kwargs): - BaseNSEMProblem.__init__(self, mesh, **kwargs) - - @property - def sigmaPrimary(self): - """ - A background model, use for the calculation of the primary fields. - - """ - return self._sigmaPrimary - @sigmaPrimary.setter - def sigmaPrimary(self, val): - # Note: TODO add logic for val, make sure it is the correct size. - self._sigmaPrimary = val - - def getA(self, freq): - """ - Function to get the A system. - - :param float freq: Frequency - :rtype: scipy.sparse.csr_matrix - :return: A - """ - Mmui = self.MfMui - Msig = self.MeSigma - C = self.mesh.edgeCurl - - 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. - - """ - - # This considers both polarizations and returns a nE,2 matrix for each polarization - if adjoint: - dMe_dsigV = sp.hstack(( self.MeSigmaDeriv( u['e_pxSolution'] ).T, self.MeSigmaDeriv(u['e_pySolution'] ).T ))*v - else: - # 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 - - - def getRHS(self, freq): - """ - 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 - """ - - # Get sources for the frequncy(polarizations) - Src = self.survey.getSrcByFreq(freq)[0] - S_e = Src.S_e(self) - return -1j * omega(freq) * S_e - - def getRHSDeriv_m(self, freq, v, adjoint=False): - """ - The derivative of the RHS with respect to sigma - """ - - Src = self.survey.getSrcByFreq(freq)[0] - S_eDeriv = Src.S_eDeriv_m(self, v, adjoint) - return -1j * omega(freq) * S_eDeriv - - def fields(self, m): - ''' - Function to calculate all the fields for the model m. - - :param np.ndarray (nC,) m: Conductivity model - ''' - # Set the current model - self.curModel = m - - 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 = self.getRHS(freq) - # Solve the system - Ainv = self.Solver(A, **self.solverOpts) - e_s = Ainv * rhs - - # Store the fields - Src = self.survey.getSrcByFreq(freq)[0] - # Store the fieldss - F[Src, 'e_pxSolution'] = e_s[:,0] - F[Src, 'e_pySolution'] = e_s[:,1] - # Note curl e = -iwb so b = -curl/iw - - if self.verbose: - print 'Ran for {:f} seconds'.format(time.time()-startTime) - sys.stdout.flush() - Ainv.clean() - return F - diff --git a/SimPEG/NSEM/Problem3D/__init__.py b/SimPEG/NSEM/Problem3D/__init__.py deleted file mode 100644 index 27e761a9..00000000 --- a/SimPEG/NSEM/Problem3D/__init__.py +++ /dev/null @@ -1 +0,0 @@ -from Probs import eForm_ps \ No newline at end of file