diff --git a/simpegMT/FieldsMT.py b/simpegMT/FieldsMT.py index 9487d122..2fc2e1c8 100644 --- a/simpegMT/FieldsMT.py +++ b/simpegMT/FieldsMT.py @@ -2,11 +2,19 @@ from SimPEG import Survey, Utils, Problem, np, sp, mkvc from scipy.constants import mu_0 import sys from numpy.lib import recfunctions as recFunc +from simpegEM.Utils.EMUtils import omega ############## ### Fields ### ############## class FieldsMT(Problem.Fields): - """Fancy Field Storage for a MT survey.""" + """Field Storage for a MT survey.""" knownFields = {'b_px': 'F','b_py': 'F', 'e_px': 'E','e_py': 'E','b_1d':'E','e_1d':'F'} dtype = complex + + + def _b_1dDeriv_u(self,src,v,adjoint=False): + """ + The derivative of b_1d wrt u + """ + return -( self.mesh.nodalGrad * v)/( 1j*omega(src.freq) ) \ No newline at end of file diff --git a/simpegMT/ProblemMT1D/Problems.py b/simpegMT/ProblemMT1D/Problems.py index 74c62292..51982255 100644 --- a/simpegMT/ProblemMT1D/Problems.py +++ b/simpegMT/ProblemMT1D/Problems.py @@ -10,15 +10,125 @@ import numpy as np import multiprocessing, sys, time -# class eForm_ps(BaseMTProblem): +class eForm_psField(BaseMTProblem): + """ + A MT problem soving a e formulation and primary/secondary fields decomposion. + Solves the equation + + """ + # From FDEMproblem: Used to project the fields. Currently not used for MTproblem. + _fieldType = 'e' + _eqLocs = 'EF' + + + def __init__(self, mesh, **kwargs): + BaseMTProblem.__init__(self, mesh, **kwargs) + + def getA(self, freq,): + """ + 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 + """ + + 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 + C = self.mesh.nodalGrad + # Make A + A = C.T*Mmui*C + 1j*omega(freq)*Msig + # 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 + dMf_dsig = self.mesh.getFaceInnerProductDeriv(self.curModel.sigma)(u) * self.curModel.sigmaDeriv + if adjoint: + return 1j * omega(freq) * ( dsig_dm.T * ( dMf_dsig.T * v ) ) + + return 1j * omega(freq) * ( dMf_dsig * ( dsig_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, u, v, adjoint=False): + """ + The derivative of the RHS wrt sigma + """ + + Src = self.survey.getSrcByFreq(freq)[0] + S_eDeriv = Src.S_eDeriv(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 = FieldsMT(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] + # Calculate total e + + e = Src.ePrimary(self) + e_s + + # Store the fields + # NOTE: only store + F[Src, 'e_1d'] = e[:,1] # Only storing the yx polarization as 1d + # F[Src, 'e_py'] = 0*e[:,0] + # Note curl e = -iwb so b = -curl e /iw + b = -( self.mesh.nodalGrad * e )/( 1j*omega(freq) ) + # F[Src, 'b_px'] = 0*b[:,0] + F[Src, 'b_1d'] = b[:,1] + if self.verbose: + print 'Ran for {:f} seconds'.format(time.time()-startTime) + sys.stdout.flush() + return F class eForm_TotalField(BaseMTProblem): """ - A MT problem solving a e formulation and a primary/secondary fields decompostion. + A MT problem solving a e formulation and a Total bondary domain decompostion. Solves the equation: + Math: + """ diff --git a/simpegMT/ProblemMT1D/__init__.py b/simpegMT/ProblemMT1D/__init__.py index 767c81df..d1908621 100644 --- a/simpegMT/ProblemMT1D/__init__.py +++ b/simpegMT/ProblemMT1D/__init__.py @@ -1 +1 @@ -from Problems import eForm_TotalField \ No newline at end of file +from Problems import eForm_TotalField, eForm_psField \ No newline at end of file diff --git a/simpegMT/ProblemMT3D/Problems.py b/simpegMT/ProblemMT3D/Problems.py index b067ebba..a2b61edd 100644 --- a/simpegMT/ProblemMT3D/Problems.py +++ b/simpegMT/ProblemMT3D/Problems.py @@ -74,7 +74,7 @@ class eForm_ps(BaseMTProblem): def getADeriv(self, freq, u, v, adjoint=False): dsig_dm = self.curModel.sigmaDeriv - dMe_dsig = self.MeSimgaDeriv( v=u) + dMe_dsig = self.MeSigmaDeriv( v=u) if adjoint: return 1j * omega(freq) * ( dsig_dm.T * ( dMe_dsig.T * v ) ) diff --git a/simpegMT/Sources/backgroundModelSources.py b/simpegMT/Sources/backgroundModelSources.py index 261cdbe5..ebd6cd37 100644 --- a/simpegMT/Sources/backgroundModelSources.py +++ b/simpegMT/Sources/backgroundModelSources.py @@ -22,26 +22,43 @@ def homo1DModelSource(mesh,freq,sigma_1d): mesh1d = simpeg.Mesh.TensorMesh([mesh.hz],np.array([mesh.x0[2]])) # # Note: Everything is using e^iwt e0_1d = get1DEfields(mesh1d,sigma_1d,freq) - # Setup x (east) polarization (_x) - ex_px = np.zeros(mesh.vnEx,dtype=complex) - ey_px = np.zeros((mesh.nEy,1),dtype=complex) - ez_px = np.zeros((mesh.nEz,1),dtype=complex) - # Assign the source to ex_x - for i in np.arange(mesh.vnEx[0]): - for j in np.arange(mesh.vnEx[1]): - ex_px[i,j,:] = -e0_1d - eBG_px = np.vstack((simpeg.Utils.mkvc(ex_px,2),ey_px,ez_px)) - # Setup y (north) polarization (_py) - ex_py = np.zeros((mesh.nEx,1), dtype='complex128') - ey_py = np.zeros(mesh.vnEy, dtype='complex128') - ez_py = np.zeros((mesh.nEz,1), dtype='complex128') - # Assign the source to ey_py - - for i in np.arange(mesh.vnEy[0]): - for j in np.arange(mesh.vnEy[1]): - ey_py[i,j,:] = e0_1d - # ey_py[1:-1,1:-1,1:-1] = 0 - eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py)) + if mesh.dim == 1: + eBG_px = -simpeg.mkvc(e0_1d,2) + eBG_py = simpeg.mkvc(e0_1d,2) + elif mesh.dim == 2: + ex_px = np.zeros(mesh.vnEx,dtype=complex) + ey_px = np.zeros((mesh.nEy,1),dtype=complex) + for i in np.arange(mesh.vnEx[0]): + ex_px[i,:] = -e0_1d + eBG_px = np.vstack((simpeg.Utils.mkvc(ex_px,2),ey_px)) + # Setup y (north) polarization (_py) + ex_py = np.zeros((mesh.nEx,1), dtype='complex128') + ey_py = np.zeros(mesh.vnEy, dtype='complex128') + # Assign the source to ey_py + for i in np.arange(mesh.vnEy[0]): + ey_py[i,:] = e0_1d + # ey_py[1:-1,1:-1,1:-1] = 0 + eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py)) + elif mesh.dim == 3: + # Setup x (east) polarization (_x) + ex_px = np.zeros(mesh.vnEx,dtype=complex) + ey_px = np.zeros((mesh.nEy,1),dtype=complex) + ez_px = np.zeros((mesh.nEz,1),dtype=complex) + # Assign the source to ex_x + for i in np.arange(mesh.vnEx[0]): + for j in np.arange(mesh.vnEx[1]): + ex_px[i,j,:] = -e0_1d + eBG_px = np.vstack((simpeg.Utils.mkvc(ex_px,2),ey_px,ez_px)) + # Setup y (north) polarization (_py) + ex_py = np.zeros((mesh.nEx,1), dtype='complex128') + ey_py = np.zeros(mesh.vnEy, dtype='complex128') + ez_py = np.zeros((mesh.nEz,1), dtype='complex128') + # Assign the source to ey_py + for i in np.arange(mesh.vnEy[0]): + for j in np.arange(mesh.vnEy[1]): + ey_py[i,j,:] = e0_1d + # ey_py[1:-1,1:-1,1:-1] = 0 + eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py)) # Return the electric fields eBG_bp = np.hstack((eBG_px,eBG_py)) diff --git a/simpegMT/SurveyMT.py b/simpegMT/SurveyMT.py index dfe174a4..6b21bfda 100644 --- a/simpegMT/SurveyMT.py +++ b/simpegMT/SurveyMT.py @@ -94,6 +94,7 @@ class RxMT(Survey.BaseRx): ex = Pex*mkvc(u[src,'e_1d'],2) bx = Pbx*mkvc(u[src,'b_1d'],2)/mu_0 f_part_complex = ex/bx + # elif self.projType is 'Z2D': elif self.projType is 'Z3D': # Get the projection Pex = mesh.getInterpolationMat(self.locs,'Ex') @@ -124,25 +125,28 @@ class RxMT(Survey.BaseRx): # Get the real or imag component real_or_imag = self.projComp f_part = getattr(f_part_complex, real_or_imag) + # print f_part return f_part - def projectFieldsDeriv(self, src, mesh, u, v, adjoint=False): - P = self.getP(mesh) + def projectFieldsDeriv(self, src, mesh, f, v, adjoint=False): + """ + The derivative of the projection wrt u + """ + real_or_imag = self.projComp if not adjoint: - Pv_complex = P * v - real_or_imag = self.projComp - Pv = getattr(Pv_complex, real_or_imag) + if self.projType is 'Z1D': + Pex = mesh.getInterpolationMat(self.locs,'Fx') + Pbx = mesh.getInterpolationMat(self.locs,'Ex') + # ex = Pex*mkvc(f[src,'e_1d'],2) + # bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0 + deriv_complex = Utils.sdiag(1/(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))*(Pex*v) - Utils.sdiag(1/(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)).T*Utils.sdiag(1/(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))*(Pbx*f._b_1dDeriv_u(src,v)/mu_0) + # elif self.projType is 'Z2D + elif self.projType is 'Z3D': + pass + Pv = getattr(deriv_complex, real_or_imag) elif adjoint: - Pv_real = P.T * v - - real_or_imag = self.projComp - if real_or_imag == 'imag': - Pv = 1j*Pv_real - elif real_or_imag == 'real': - Pv = Pv_real.astype(complex) - else: - raise NotImplementedError('must be real or imag') + raise NotImplementedError('must be real or imag') return Pv @@ -187,7 +191,7 @@ class srcMT_polxy_1Dprimary(srcMT): as fields in the full space of the problem. """ def __init__(self, rxList, freq, sigma1d): - 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).' + # 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).' self.sigma1d = sigma1d srcMT.__init__(self, rxList, freq) @@ -201,7 +205,7 @@ class srcMT_polxy_1Dprimary(srcMT): def bPrimary(self,problem): # Project ePrimary to bPrimary # Satisfies the primary(background) field conditions - bBG_bp = (- self.mesh.edgeCurl * self.ePrimary )/( 1j*omega(freq) ) + bBG_bp = (- problem.mesh.edgeCurl * self.ePrimary )/( 1j*omega(freq) ) return bBG_bp def S_e(self,problem): @@ -213,12 +217,25 @@ class srcMT_polxy_1Dprimary(srcMT): sigma_p = Map_sigma_p._transform(self.sigma1d) # Make mass matrix # Note: M(sig) - M(sig_p) = M(sig - sig_p) - Mesigma = problem.MeSigma - Mesigma_p = problem.mesh.getEdgeInnerProduct(sigma_p) + # Need to deal with the edge/face discrepencies between 1d/2d/3d + if problem.mesh.dim == 1: + Mesigma = problem.mesh.getFaceInnerProduct(problem.curModel.sigma) + Mesigma_p = problem.mesh.getFaceInnerProduct(sigma_p) + if problem.mesh.dim == 2: + pass + if problem.mesh.dim == 3: + Mesigma = problem.MeSigma + Mesigma_p = problem.mesh.getEdgeInnerProduct(sigma_p) return (Mesigma - Mesigma_p) * e_p def S_eDeriv(self, problem, v, adjoint = False): - MesigmaDeriv = problem.MeSigmaDeriv(self.ePrimary(problem)) + # Need to deal with + if problem.mesh.dim == 1: + pass + if problem.mesh.dim == 2: + pass + if problem.mesh.dim == 3: + MesigmaDeriv = problem.MeSigmaDeriv(self.ePrimary(problem)) if adjoint: return MesigmaDeriv.T * v else: diff --git a/simpegMT/Tests/test_Problem1D_againstAnalyticHalfspace.py b/simpegMT/Tests/test_Problem1D_againstAnalyticHalfspace.py index 8aa54880..b8a50980 100644 --- a/simpegMT/Tests/test_Problem1D_againstAnalyticHalfspace.py +++ b/simpegMT/Tests/test_Problem1D_againstAnalyticHalfspace.py @@ -8,7 +8,7 @@ TOLr = 5e-2 TOLp = 5e-2 -def setupSurvey(sigmaHalf): +def setupSurvey(sigmaHalf,tD=True): # Frequency nFreq = 33 @@ -31,8 +31,13 @@ def setupSurvey(sigmaHalf): rxList.append(simpegmt.SurveyMT.RxMT(simpeg.mkvc(np.array([0.0]),2).T,rxType)) # Source list srcList =[] - for freq in freqs: - srcList.append(simpegmt.SurveyMT.srcMT_polxy_1DhomotD(rxList,freq)) + if tD: + for freq in freqs: + srcList.append(simpegmt.SurveyMT.srcMT_polxy_1DhomotD(rxList,freq)) + else: + for freq in freqs: + srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freq,sigma)) + survey = simpegmt.SurveyMT.SurveyMT(srcList) return survey, sigma, m1d @@ -90,23 +95,68 @@ def appPhs_TotalFieldNorm(sigmaHalf): return np.linalg.norm(np.abs(app_p - np.ones(survey.nFreq)*135)/ 135) +def appRes_psFieldNorm(sigmaHalf): + + # Make the survey + survey, sigma, mesh = setupSurvey(sigmaHalf,False) + problem = simpegmt.ProblemMT1D.eForm_psField(mesh) + problem.pair(survey) + + # Get the fields + fields = problem.fields(sigma) + + # Project the data + data = survey.projectFields(fields) + + # Calculate the app res and phs + app_r = np.array(getAppResPhs(data))[:,0] + + return np.linalg.norm(np.abs(app_r - np.ones(survey.nFreq)/sigmaHalf)*sigmaHalf) + +def appPhs_psFieldNorm(sigmaHalf): + + # Make the survey + survey, sigma, mesh = setupSurvey(sigmaHalf,False) + problem = simpegmt.ProblemMT1D.eForm_psField(mesh) + problem.pair(survey) + + # Get the fields + fields = problem.fields(sigma) + + # Project the data + data = survey.projectFields(fields) + + # Calculate the app phs + app_p = np.array(getAppResPhs(data))[:,1] + + return np.linalg.norm(np.abs(app_p - np.ones(survey.nFreq)*135)/ 135) + class TestAnalytics(unittest.TestCase): def setUp(self): pass - def test_appRes2en1(self):self.assertLess(appRes_TotalFieldNorm(2e-1), TOLr) - def test_appRes2en2(self):self.assertLess(appRes_TotalFieldNorm(2e-2), TOLr) - def test_appRes2en3(self):self.assertLess(appRes_TotalFieldNorm(2e-3), TOLr) - def test_appRes2en4(self):self.assertLess(appRes_TotalFieldNorm(2e-4), TOLr) - def test_appRes2en5(self):self.assertLess(appRes_TotalFieldNorm(2e-5), TOLr) - def test_appRes2en6(self):self.assertLess(appRes_TotalFieldNorm(2e-6), TOLr) - def test_appPhs2en1(self):self.assertLess(appPhs_TotalFieldNorm(2e-1), TOLp) - def test_appPhs2en2(self):self.assertLess(appPhs_TotalFieldNorm(2e-2), TOLp) - def test_appPhs2en3(self):self.assertLess(appPhs_TotalFieldNorm(2e-3), TOLp) - def test_appPhs2en4(self):self.assertLess(appPhs_TotalFieldNorm(2e-4), TOLp) - def test_appPhs2en5(self):self.assertLess(appPhs_TotalFieldNorm(2e-5), TOLp) - def test_appPhs2en6(self):self.assertLess(appPhs_TotalFieldNorm(2e-6), TOLp) + # Total Fields + # def test_appRes2en1(self):self.assertLess(appRes_TotalFieldNorm(2e-1), TOLr) + # def test_appPhs2en1(self):self.assertLess(appPhs_TotalFieldNorm(2e-1), TOLp) + def test_appRes2en2(self):self.assertLess(appRes_TotalFieldNorm(2e-2), TOLr) + def test_appPhs2en2(self):self.assertLess(appPhs_TotalFieldNorm(2e-2), TOLp) + + # def test_appRes2en3(self):self.assertLess(appRes_TotalFieldNorm(2e-3), TOLr) + # def test_appPhs2en3(self):self.assertLess(appPhs_TotalFieldNorm(2e-3), TOLp) + + # def test_appRes2en4(self):self.assertLess(appRes_TotalFieldNorm(2e-4), TOLr) + # def test_appPhs2en4(self):self.assertLess(appPhs_TotalFieldNorm(2e-4), TOLp) + + # def test_appRes2en5(self):self.assertLess(appRes_TotalFieldNorm(2e-5), TOLr) + # def test_appPhs2en5(self):self.assertLess(appPhs_TotalFieldNorm(2e-5), TOLp) + + # def test_appRes2en6(self):self.assertLess(appRes_TotalFieldNorm(2e-6), TOLr) + # def test_appPhs2en6(self):self.assertLess(appPhs_TotalFieldNorm(2e-6), TOLp) + + # Primary/secondary + def test_appRes2en2_ps(self):self.assertLess(appRes_psFieldNorm(2e-2), TOLr) + def test_appPhs2en2_ps(self):self.assertLess(appPhs_psFieldNorm(2e-2), TOLp) if __name__ == '__main__': unittest.main() \ No newline at end of file