Fixed directive to save iteration dictionary.

Fixing MT namespace to NSEM.

.travis.yml

@@ -35,7 +35,7 @@ before_install:
 
 # Install packages
 install:
-  - conda install --yes pip python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib cython ipython nose vtk
+  - conda install --yes pip python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib cython ipython ipywidgets nose vtk
   - pip install nose-cov python-coveralls
 
   - git clone https://github.com/rowanc1/pymatsolver.git

SimPEG/Directives.py

@@ -213,42 +213,33 @@ class SaveOutputEveryIteration(_SaveEveryIteration):
         f.close()
 
 class SaveOutputDictEveryIteration(_SaveEveryIteration):
-    """SaveOutputDictEveryIteration"""
+    """
+    Saves inversion parameters at every iteraion.
+
+
+    """
 
     def initialize(self):
         print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '###-%s.npz'"%self.fileName
 
     def endIter(self):
-        # Save the data.
-        ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.reg.mref) )
-        phi_ms = 0.5*ms.dot(ms)
-        if self.reg.mrefInSmooth == True:
-            mref = self.reg.mref
-        else:
-            mref = 0
-        mx = self.reg.Wx * ( self.reg.mapping * (self.invProb.curModel - mref) )
-        phi_mx = 0.5 * mx.dot(mx)
-        if self.prob.mesh.dim >= 2:
-            my = self.reg.Wy * ( self.reg.mapping * (self.invProb.curModel - mref) )
-            phi_my = 0.5 * my.dot(my)
-        else:
-            phi_my = 'NaN'
-        if self.prob.mesh.dim==3:
-            mz = self.reg.Wz * ( self.reg.mapping * (self.invProb.curModel - mref) )
-            phi_mz = 0.5 * mz.dot(mz)
-        else:
-            phi_mz = 'NaN'
 
+        # Initialize the output dict
+        outDict = {}
+        # Save the data.
+        outDict['iter'] = self.opt.iter
+        outDict['beta'] = self.invProb.beta
+        outDict['phi_d'] = self.invProb.phi_d
+        outDict['phi_ms'] = self.reg._evalSmall(self.invProb.curModel)
+        outDict['phi_mx'] = self.reg._evalSmoothx(self.invProb.curModel)
+        outDict['phi_my'] = self.reg._evalSmoothy(self.invProb.curModel) if self.prob.mesh.dim >= 2 else 'NaN'
+        outDict['phi_mz'] = self.reg._evalSmoothz(self.invProb.curModel) if self.prob.mesh.dim==3 else 'NaN'
+        outDict['f'] = self.opt.f
+        outDict['m'] = self.invProb.curModel
+        outDict['dpred'] = self.invProb.dpred
 
         # Save the file as a npz
-        np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred)
+        np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), outDict)
-
-#         mref = getattr(self, 'm_prev', None)
-#         if mref is None:
-#             if self.debug: print 'UpdateReferenceModel is using mref0'
-#             mref = self.mref0
-#         self.m_prev = self.invProb.m_current
-#         return mref
 
 class Update_IRLS(InversionDirective):
 

SimPEG/Examples/MT_1D_ForwardAndInversion.py

@@ -1,9 +1,11 @@
 import SimPEG as simpeg
 import numpy as np
-import SimPEG.MT as MT
+from SimPEG import NSEM
 from scipy.constants import mu_0
 import matplotlib.pyplot as plt
 
+np.random.seed(1983)
+
 def run(plotIt=True):
     """
         MT: 1D: Inversion
@@ -17,13 +19,13 @@ def run(plotIt=True):
     ## Setup the forward modeling
     # Setting up 1D mesh and conductivity models to forward model data.
     # Frequency
-    nFreq = 31
+    nFreq = 26
-    freqs = np.logspace(3,-3,nFreq)
+    freqs = np.logspace(2,-3,nFreq)
     # Set mesh parameters
-    ct = 20
+    ct = 10
-    air = simpeg.Utils.meshTensor([(ct,16,1.4)])
+    air = simpeg.Utils.meshTensor([(ct,25,1.4)])
     core = np.concatenate( (  np.kron(simpeg.Utils.meshTensor([(ct,10,-1.3)]),np.ones((5,))) , simpeg.Utils.meshTensor([(ct,5)]) ) )
-    bot = simpeg.Utils.meshTensor([(core[0],10,-1.4)])
+    bot = simpeg.Utils.meshTensor([(core[0],25,-1.4)])
     x0 = -np.array([np.sum(np.concatenate((core,bot)))])
     # Make the model
     m1d = simpeg.Mesh.TensorMesh([np.concatenate((bot,core,air))], x0=x0)
@@ -33,7 +35,7 @@ def run(plotIt=True):
     layer1 = (m1d.vectorCCx<-500.) & (m1d.vectorCCx>=-800.)
     layer2 = (m1d.vectorCCx<-3500.) & (m1d.vectorCCx>=-5000.)
     # Set the conductivity values
-    sig_half = 2e-3
+    sig_half = 1e-2
     sig_air = 1e-8
     sig_layer1 = .2
     sig_layer2 = .2
@@ -50,38 +52,38 @@ def run(plotIt=True):
     m_0 = np.log(sigma_0[active])
 
     # Set the mapping
-    actMap = simpeg.Maps.ActiveCells(m1d, active, np.log(1e-8), nC=m1d.nCx)
+    actMap = simpeg.Maps.InjectActiveCells(m1d, active, np.log(1e-8), nC=m1d.nCx)
     mappingExpAct = simpeg.Maps.ExpMap(m1d) * actMap
 
     ## Setup the layout of the survey, set the sources and the connected receivers
     # Receivers
     rxList = []
     for rxType in ['z1dr','z1di']:
-        rxList.append(MT.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType))
+        rxList.append(NSEM.Rx(simpeg.mkvc(np.array([-0.5]),2).T,rxType))
     # Source list
     srcList =[]
     for freq in freqs:
-            srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
+            srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq))
     # Make the survey
-    survey = MT.Survey(srcList)
+    survey = NSEM.Survey(srcList)
     survey.mtrue = m_true
 
     ## Set the problem
-    problem = MT.Problem1D.eForm_psField(m1d,sigmaPrimary=sigma_0,mapping=mappingExpAct)
+    problem = NSEM.Problem1D_ePrimSec(m1d,sigmaPrimary=sigma_0,mapping=mappingExpAct)
     problem.pair(survey)
 
     ## Forward model data
     # Project the data
     survey.dtrue = survey.dpred(m_true)
-    survey.dobs = survey.dtrue + 0.025*abs(survey.dtrue)*np.random.randn(*survey.dtrue.shape)
+    survey.dobs = survey.dtrue + 0.01*abs(survey.dtrue)*np.random.randn(*survey.dtrue.shape)
 
     if plotIt:
-        fig = MT.Utils.dataUtils.plotMT1DModelData(problem)
+        fig = NSEM.Utils.dataUtils.plotMT1DModelData(problem,[])
         fig.suptitle('Target - smooth true')
 
 
     # Assign uncertainties
-    std = 0.05 # 5% std
+    std = 0.025 # 5% std
     survey.std = np.abs(survey.dobs*std)
     # Assign the data weight
     Wd = 1./survey.std
@@ -90,30 +92,33 @@ def run(plotIt=True):
     # Define a counter
     C =  simpeg.Utils.Counter()
     # Set the optimization
-    opt = simpeg.Optimization.InexactGaussNewton(maxIter = 30)
+    opt = simpeg.Optimization.ProjectedGNCG(maxIter = 25)
     opt.counter = C
-    opt.LSshorten = 0.5
+    opt.lower = np.log(1e-4)
+    opt.upper = np.log(5)
+    opt.LSshorten = 0.1
     opt.remember('xc')
     # Data misfit
     dmis = simpeg.DataMisfit.l2_DataMisfit(survey)
     dmis.Wd = Wd
     # Regularization - with a regularization mesh
-    regMesh = simpeg.Mesh.TensorMesh([m1d.hx[problem.mapping.sigmaMap.maps[-1].indActive]],m1d.x0)
+    regMesh = simpeg.Mesh.TensorMesh([m1d.hx[active]],m1d.x0)
     reg = simpeg.Regularization.Tikhonov(regMesh)
     reg.mrefInSmooth = True
-    reg.alpha_s = 1e-7
+    reg.alpha_s = 1e-1
     reg.alpha_x = 1.
+
     # Inversion problem
     invProb = simpeg.InvProblem.BaseInvProblem(dmis, reg, opt)
     invProb.counter = C
     # Beta cooling
     beta = simpeg.Directives.BetaSchedule()
-    beta.coolingRate = 4
+    beta.coolingRate = 4.
-    betaest = simpeg.Directives.BetaEstimate_ByEig(beta0_ratio=0.75)
+    beta.coolingFactor = 4.
+    betaest = simpeg.Directives.BetaEstimate_ByEig(beta0_ratio=1.)
+    betaest.beta0 = 1.
     targmis = simpeg.Directives.TargetMisfit()
     targmis.target = survey.nD
-    saveModel = simpeg.Directives.SaveModelEveryIteration()
-    saveModel.fileName = 'Inversion_TargMisEqnD_smoothTrue'
     # Create an inversion object
     inv = simpeg.Inversion.BaseInversion(invProb, directiveList=[beta,betaest,targmis])
 
@@ -121,8 +126,9 @@ def run(plotIt=True):
     mopt = inv.run(m_0)
 
     if plotIt:
-        fig = MT.Utils.dataUtils.plotMT1DModelData(problem,[mopt])
+        fig = NSEM.Utils.dataUtils.plotMT1DModelData(problem,[mopt])
         fig.suptitle('Target - smooth true')
+        fig.axes[0].set_ylim([-10000,500])
         plt.show()
 
 if __name__ == '__main__':

SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py

@@ -0,0 +1,428 @@
+from scipy.constants import epsilon_0, mu_0
+import matplotlib.pyplot as plt
+import numpy as np
+from ipywidgets import *
+from SimPEG.EM.Utils import k, omega
+
+"""
+MT1D: n layered earth problem
+*****************************
+
+Author: Thibaut Astic
+Contact: thast@eos.ubc.ca
+Date: January 2016
+
+This code compute the analytic response of a n-layered Earth to a plane wave (Magneto-Tellurics).
+
+We start by looking at Maxwell's equations in the electric
+field \\\(\\\mathbf{E}\\) and the magnetic flux
+\\\(\\\mathbf{H}\\) to write the wave equations
+\\(\\ \nabla ^2  \mathbf{E_x} + k^2 \mathbf{E_x} = 0 \\) &
+\\(\\ \nabla ^2  \mathbf{H_y} + k^2 \mathbf{H_y} = 0 \\)
+
+Then solving the equations in each layer "j" between z_{j-1} and z_j in the form of
+\\(\\ E_{x,j} (z) = U_j e^{i k (z-z_{j-1})} + D_j e^{-i k (z-z_{j-1})} \\)
+\\(\\ H_{y,j} (z) = \frac{1}{Z_j} (D_j e^{-i k (z-z_{j-1})} - U_j e^{i k (z-z_{j-1})}) \\)
+
+With U and D the Up and Down components of the E-field.
+
+The iteration from one layer to another is ensure by:
+
+\\(\\  \left(\begin{matrix} E_{x,j} \\ H_{y,j} \end{matrix} \right) =
+    P_j T_j P^{-1}_J \left(\begin{matrix} E_{x,j+1} \\ H_{y,j+1} \end{matrix} \right) \\)
+
+And the Boundary Condition is set for the E-field in the last layer, with no Up component (=0)
+and only a down component (=1 then normalized by the highest amplitude to ensure numeric stability)
+
+The layer 0 is assumed to be the air layer.
+
+"""
+
+#Define a frquency range for a survey
+frange = lambda minfreq, maxfreq, step: np.logspace(minfreq,maxfreq,num = step, base = 10.)
+
+#Functions to create random physical Perties for a n-layered earth
+thick = lambda minthick, maxthick, nlayer: np.append(np.array([1.2*10.**5]),
+                                                     np.ndarray.round(minthick + (maxthick-minthick)* np.random.rand(nlayer-1,1)
+                                                            ,decimals =1))
+
+sig = lambda minsig, maxsig, nlayer: np.append(np.array([0.]),
+                                               np.ndarray.round(10.**minsig + (10.**maxsig-10.**minsig)* np.random.rand(nlayer,1)
+                                                      ,decimals=3))
+
+mu  = lambda minmu, maxmu, nlayer: np.append(np.array([1.]),
+                                             np.ndarray.round(minmu + (maxmu-minmu)* np.random.rand(nlayer,1)
+                                                    ,decimals=1))
+
+eps = lambda mineps, maxeps, nlayer: np.append(np.array([1.]),
+                                               np.ndarray.round(mineps + (maxeps-mineps)* np.random.rand(nlayer,1)
+                                                                ,decimals=1))
+
+#Evaluate Impedance Z of a layer
+ImpZ = lambda f, mu, k: omega(f)*mu*mu_0/k
+
+#Complex Cole-Cole Conductivity - EM utils
+PCC= lambda siginf,m,t,c,f: siginf*(1.-(m/(1.+(1j*omega(f)*t)**c)))
+
+#Converted thickness array into top of layer array
+top = lambda thick: np.cumsum(thick)
+
+#Propagation Matrix and theirs inverses
+
+#matrix T for transition of Up and Down components accross a layer
+T = lambda h,k: np.matrix([[np.exp(1j*k*h),0.],[0.,np.exp(-1j*k*h)]],dtype='complex_')
+
+Tinv = lambda h,k: np.matrix([[np.exp(-1j*k*h),0.],[0.,np.exp(1j*k*h)]],dtype='complex_')
+
+#transition of Up and Down components accross a layer
+UD_Z = lambda UD,z,zj,k : T((z-zj),k)*UD
+
+
+#matrix P relating Up and Down components with E and H fields
+P = lambda z: np.matrix([[1.,1,],[-1./z,1./z]],dtype='complex_')
+
+Pinv = lambda z: np.matrix([[1.,-z],[1.,z]],dtype='complex_')/2.
+
+
+#Time Variation of E and H
+E_ZT = lambda U,D,f,t : np.exp(1j*omega(f)*t)*(U+D)
+H_ZT = lambda U,D,Z,f,t : (1./Z)*np.exp(1j*omega(f)*t)*(D-U)
+
+#Plot the configuration of the problem
+def PlotConfiguration(thick,sig,eps,mu,ax,widthg,z):
+
+    topn = top(thick)
+    widthn = np.arange(-widthg,widthg+widthg/10.,widthg/10.)
+
+    ax.set_ylim([z.min(),z.max()])
+    ax.set_xlim([-widthg,widthg])
+
+    ax.set_ylabel("Depth (m)", fontsize=16.)
+    ax.yaxis.tick_right()
+    ax.yaxis.set_label_position("right")
+
+    #define filling for the different layers
+    hatches=['/'  , '+', 'x', '|' , '\\', '-'  , 'o' , 'O' , '.' , '*' ]
+
+    #Write the physical properties of air
+    ax.annotate(("Air, $\sigma$ =%1.0f mS/m")%(sig[0]*10**(3)),
+            xy=(-widthg/2., -np.abs(z.max())/2.), xycoords='data',
+            xytext=(-widthg/2., -np.abs(z.max())/2.), textcoords='data',
+            fontsize=14.)
+
+    ax.annotate(("$\epsilon_r$= %1i")%(eps[0]),
+            xy=(-widthg/2., -np.abs(z.max())/3.), xycoords='data',
+            xytext=(-widthg/2., -np.abs(z.max())/3.), textcoords='data',
+            fontsize=14.)
+
+    ax.annotate(("$\mu_r$= %1i")%(mu[0]),
+            xy=(-widthg/2., -np.abs(z.max())/3.), xycoords='data',
+            xytext=(0, -np.abs(z.max())/3.), textcoords='data',
+            fontsize=14.)
+
+    #Write the physical properties of the differents layers up to the (n-1)-th and fill it with pattern
+    for i in range(1,len(topn)-1,1):
+        if topn[i] == topn[i+1]:
+            pass
+        else:
+            ax.annotate(("$\sigma$ =%3.3f mS/m")%(sig[i]*10**(3)),
+                xy=(0., (2.*topn[i]+topn[i+1])/3), xycoords='data',
+                xytext=(0., (2.*topn[i]+topn[i+1])/3), textcoords='data',
+                fontsize=14.)
+
+            ax.annotate(("$\epsilon_r$= %1i")%(eps[i]),
+                xy=(-widthg/1.1, (2.*topn[i]+topn[i+1])/3), xycoords='data',
+                xytext=(-widthg/1.1, (2.*topn[i]+topn[i+1])/3), textcoords='data',
+                fontsize=14.)
+
+            ax.annotate(("$\mu_r$= %1.2f")%(mu[i]),
+                xy=(-widthg/2., (2.*topn[i]+topn[i+1])/3), xycoords='data',
+                xytext=(-widthg/2., (2.*topn[i]+topn[i+1])/3), textcoords='data',
+                fontsize=14.)
+
+            ax.plot(widthn,topn[i]*np.ones_like(widthn),color='black')
+            ax.fill_between(widthn,topn[i],topn[i+1],alpha=0.3,color="none",edgecolor='black', hatch=hatches[(i-1)%10])
+
+    #Write the physical properties of the n-th layer and fill it with pattern
+    ax.plot(widthn,topn[-1]*np.ones_like(widthn),color='black')
+    ax.fill_between(widthn,topn[-1],z.max(),alpha=0.3,color="none",edgecolor='black', hatch=hatches[(len(topn)-2)%10])
+
+    ax.annotate(("$\sigma$ =%3.3f mS/m")%(sig[-1]*10**(3)),
+            xy=(0., (2.*topn[-1]+z.max())/3), xycoords='data',
+            xytext=(0., (2.*topn[-1]+z.max())/3), textcoords='data',
+            fontsize=14.)
+
+    ax.annotate(("$\epsilon_r$= %1i")%(eps[-1]),
+            xy=(-widthg/1.1, (2.*topn[-1]+z.max())/3), xycoords='data',
+            xytext=(-widthg/1.1, (2.*topn[-1]+z.max())/3), textcoords='data',
+            fontsize=14.)
+
+    ax.annotate(("$\mu_r$= %1.2f")%(mu[-1]),
+            xy=(-widthg/2., (2.*topn[-1]+z.max())/3), xycoords='data',
+            xytext=(-widthg/2., (2.*topn[-1]+z.max())/3), textcoords='data',
+            fontsize=14.)
+
+    #plot Trees!
+    ax.annotate("",
+            xy=(widthg/2., -1.*z.max()/5.), xycoords='data',
+            xytext=(widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.2,head_length=1.2',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(widthg/2., -3./4.*z.max()/5.), xycoords='data',
+            xytext=(widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.4,head_length=1.4',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(widthg/2., -1./2.*z.max()/5.), xycoords='data',
+            xytext=(widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(1.2*widthg/2., -1.*z.max()/5.), xycoords='data',
+            xytext=(1.2*widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.2,head_length=1.2',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(1.2*widthg/2., -3./4.*z.max()/5.), xycoords='data',
+            xytext=(1.2*widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.4,head_length=1.4',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(1.2*widthg/2., -1./2.*z.max()/5.), xycoords='data',
+            xytext=(1.2*widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(1.5*widthg/2., -1.*z.max()/5.), xycoords='data',
+            xytext=(1.5*widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.2,head_length=1.2',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(1.5*widthg/2., -3./4.*z.max()/5.), xycoords='data',
+            xytext=(1.5*widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.4,head_length=1.4',color='green',linewidth=2.)
+            )
+
+    ax.annotate("",
+            xy=(1.5*widthg/2., -1./2.*z.max()/5.), xycoords='data',
+            xytext=(1.5*widthg/2., 0.), textcoords='data',
+            arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.)
+            )
+
+
+    ax.invert_yaxis()
+
+    return ax
+
+#Propagate Up and Down component for a certain frequency & evaluate E and H field
+
+def Propagate(f,H,sig,chg,taux,c,mu,eps,n):
+
+    sigcm = np.zeros_like(sig,dtype='complex_')
+
+    for j in range(1,len(sig)):
+        sigcm[j]=PCC(sig[j],chg[j],taux[j],c[j],f)
+
+    K = k(f, sigcm, mu, eps)
+    Z = ImpZ(f,mu,K)
+
+    EH = np.matrix(np.zeros((2,n+1),dtype = 'complex_'),dtype = 'complex_')
+    UD = np.matrix(np.zeros((2,n+1),dtype = 'complex_'),dtype = 'complex_')
+
+    UD[1,-1] = 1.
+
+    for i in range(-2,-(n+2),-1):
+
+        UD[:,i] = Tinv(H[i+1],K[i])*Pinv(Z[i])*P(Z[i+1])*UD[:,i+1]
+        UD = UD/((np.abs(UD[0,:]+UD[1,:])).max())
+
+    for j in range(0,n+1):
+        EH[:,j] = np.matrix([[1.,1,],[-1./Z[j],1./Z[j]]])*UD[:,j]
+
+    return UD, EH, Z ,K
+
+
+#Evaluate the apparent resistivity and phase for a frequency range
+def appres(F,H,sig,chg,taux,c,mu,eps,n):
+
+    Res = np.zeros_like(F)
+    Phase = np.zeros_like(F)
+    App_ImpZ= np.zeros_like(F,dtype='complex_')
+
+    for i in range(0,len(F)):
+
+        UD,EH,Z ,K = Propagate(F[i],H,sig,chg,taux,c,mu,eps,n)
+
+        App_ImpZ[i] = EH[0,1]/EH[1,1]
+
+        Res[i] = np.abs(App_ImpZ[i])**2./(mu_0*omega(F[i]))
+        Phase[i] = np.angle(App_ImpZ[i], deg = True)
+
+    return Res,Phase
+
+#Evaluate Up, Down components, E and H field, for a frequency range,
+#a discretized depth range and a time range (use to calculate envelope)
+def calculateEHzt(F,H,sig,chg,taux,c,mu,eps,n,zsample,tsample):
+
+    topc = top(H)
+
+    layer = np.zeros(len(zsample),dtype=np.int)-1
+
+    Exzt = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
+    Hyzt = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
+    Uz = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
+    Dz = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
+    UDaux = np.matrix(np.zeros((2,len(zsample)),dtype = 'complex_'),dtype = 'complex_')
+
+    for i in range(0,n+1,1):
+        layer = layer+(zsample>=topc[i])*1
+
+    for j in range(0,len(F)):
+
+        UD,EH,Z ,K = Propagate(F[j],H,sig,chg,taux,c,mu,eps,n)
+
+        for p in range(0,len(zsample)):
+
+            UDaux[:,p] = UD_Z(UD[:,layer[p]],zsample[p],topc[layer[p]],K[layer[p]])
+
+            for q in range(0,len(tsample)):
+
+                Exzt[p,q]  = Exzt[p,q] + E_ZT(UDaux[0,p],UDaux[1,p],F[j],tsample[q])/len(F)
+                Hyzt[p,q] = Hyzt[p,q] + H_ZT(UDaux[0,p],UDaux[1,p],Z[layer[p]],F[j],tsample[q])/len(F)
+                Uz[p,q] = Uz[p,q] + UDaux[0,p]*np.exp(1j*omega(F[j])*tsample[q])/len(F)
+                Dz[p,q] = Dz[p,q] + UDaux[1,p]*np.exp(1j*omega(F[j])*tsample[q])/len(F)
+
+    return  Exzt,Hyzt,Uz,Dz,UDaux,layer
+
+
+#Function to Plot Apparent Resistivity and Phase
+def PlotAppRes(F,H,sig,chg,taux,c,mu,eps,n,fenvelope,PlotEnvelope):
+
+    Res, Phase = appres(F,H,sig,chg,taux,c,mu,eps,n)
+
+    fig,ax = plt.subplots(1,2,figsize=(16,10))
+
+    ax[0].scatter(Res,F,color='black')
+    ax[0].set_xscale('Log')
+    ax[0].set_yscale('Log')
+    ax[0].set_xlim([10.**(np.log10(Res.min())-1.),10.**(np.log10(Res.max())+1.)])
+    ax[0].set_ylim([F.min(),F.max()])
+    ax[0].set_xlabel('Apparent Resistivity (Ohm*m)',fontsize=16.,color="black")
+    ax[0].set_ylabel('Frequency (Hz)',fontsize=16.)
+    ax[0].grid(which='major')
+
+    ax0 = ax[0].twiny()
+
+    ax0.set_xlim([0.,90.])
+    ax0.set_ylim([F.min(),F.max()])
+    ax0.scatter(Phase,F,color='purple')
+    ax0.set_xlabel('Phase (Degrees)',fontsize=16.,color="purple")
+
+    zc=np.arange(-(H[1:].max()+10)*n,(H[1:].max()+10)*n,10.)
+
+    ax[0].tick_params(labelsize=16)
+    ax[1].tick_params(labelsize=16)
+    ax0.tick_params(labelsize=16)
+
+    if PlotEnvelope:
+
+        widthn=np.logspace(np.log10(Res.min())-1., np.log10(Res.max())+1., num=100, endpoint=True, base=10.0)
+        fenvelope1n=np.ones(100)*fenvelope
+        ax[0].plot(widthn,fenvelope1n,linestyle='dashed',color='black')
+
+        tc=np.arange(0.,1./fenvelope,0.01/(fenvelope))
+        Exzt,Hyzt,Uz,Dz,UDaux,layer = calculateEHzt(np.array([fenvelope]),H,sig,chg,taux,c,mu,eps,n,zc,tc)
+
+        ax1=ax[1].twiny()
+
+        ax[1].tick_params(labelsize=16)
+        ax1.tick_params(labelsize=16)
+
+        ax[1].set_xlabel('Amplitude Electric Field E (V/m)',color='blue',fontsize=16)
+
+        ax1.set_xlabel('Amplitude Magnetic Field H (A/m)',color='red',fontsize=16)
+
+        ax[1].fill_betweenx(zc,np.squeeze(np.asarray(np.real(Exzt.min(axis=1)))),
+                      np.squeeze(np.asarray(np.real(Exzt.max(axis=1)))),
+                      color='blue', alpha=0.1)
+
+        ax1.fill_betweenx(zc,np.squeeze(np.asarray(np.real(Hyzt.min(axis=1)))),
+                      np.squeeze(np.asarray(np.real(Hyzt.max(axis=1)))),
+                      color='red', alpha=0.1)
+
+        ax[1] = PlotConfiguration(H,sig,eps,mu,ax[1],(1.5*np.abs(Exzt).max()),zc)
+        ax1.set_xlim([-1.5*np.abs(Hyzt).max(),1.5*np.abs(Hyzt).max()])
+        ax1.set_xlim([-1.5*np.abs(Hyzt).max(),1.5*np.abs(Hyzt).max()])
+    else:
+        print 'No envelop (if True, might be slow)'
+        ax[1] = PlotConfiguration(H,sig,eps,mu,ax[1],1.,zc)
+        ax[1].get_xaxis().set_ticks([])
+
+    plt.show()
+
+#Interactive MT for Notebook
+def PlotAppRes3LayersInteract(h1,h2,sigl1,sigl2,sigl3,mul1,mul2,mul3,epsl1,epsl2,epsl3,PlotEnvelope,F_Envelope):
+
+    frangn=frange(-5,5,100.)
+    sig3= np.array([0.,0.001,0.1, 0.001])
+    thick3 = np.array([120000.,50.,50.])
+    eps3=np.array([1.,1.,1.,1])
+    mu3=np.array([1.,1.,1.,1])
+    chg3=np.array([0.,0.1,0.,0.2])
+    chg3_0=np.array([0.,0.1,0.,0.])
+    taux3=np.array([0.,0.1,0.,0.1])
+    c3=np.array([1.,1.,1.,1.])
+
+    sig3[1]=sigl1
+    sig3[1]=10.**sig3[1]
+    sig3[2]=sigl2
+    sig3[2]=10.**sig3[2]
+    sig3[3]=sigl3
+    sig3[3]=10.**sig3[3]
+    mu3[1]=mul1
+    mu3[2]=mul2
+    mu3[3]=mul3
+    eps3[1]=epsl1
+    eps3[2]=epsl2
+    eps3[3]=epsl3
+    thick3[1]=h1
+    thick3[2]=h2
+
+    PlotAppRes(frangn,thick3,sig3,chg3_0,taux3,c3,mu3,eps3,3,F_Envelope,PlotEnvelope)
+
+
+def run(n=3,plotIt=True):
+    # something to make a plot
+
+    F = frange(-5.,5.,20)
+    H = thick(50.,100.,n)
+    sign = sig(-5.,0.,n)
+    mun = mu(1.,2.,n)
+    epsn = eps(1.,9.,n)
+    chg = np.zeros_like(sign)
+    taux = np.zeros_like(sign)
+    c = np.zeros_like(sign)
+
+    Res, Phase = appres(F,H,sign,chg,taux,c,mun,epsn,n)
+
+    if plotIt:
+
+        PlotAppRes(F, H, sign, chg, taux, c, mun, epsn, n, fenvelope=1000., PlotEnvelope=True)
+
+    return Res, Phase
+
+if __name__ == '__main__':
+    run()
+
+
+
+
+

SimPEG/Examples/MT_3D_Foward.py

@@ -2,7 +2,7 @@
 
 # Import
 import SimPEG as simpeg
-from SimPEG import MT
+from SimPEG import NSEM
 import numpy as np
 try:
     from pymatsolver import MumpsSolver as Solver
@@ -37,16 +37,16 @@ def run(plotIt=True, nFreq=1):
     for loc in rx_loc:
         # NOTE: loc has to be a (1,3) np.ndarray otherwise errors accure
         for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi','tzxr','tzxi','tzyr','tzyi']:
-            rxList.append(MT.Rx(simpeg.mkvc(loc,2).T,rxType))
+            rxList.append(NSEM.Rx(simpeg.mkvc(loc,2).T,rxType))
     # Source list
     srcList =[]
     for freq in np.logspace(3,-3,nFreq):
-        srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
+        srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq))
     # Survey MT
-    survey = MT.Survey(srcList)
+    survey = NSEM.Survey(srcList)
 
     ## Setup the problem object
-    problem = MT.Problem3D.eForm_ps(M, sigmaPrimary=sigBG)
+    problem = NSEM.Problem3D_ePrimSec(M, sigmaPrimary=sigBG)
     problem.pair(survey)
     problem.Solver = Solver
 
@@ -55,7 +55,7 @@ def run(plotIt=True, nFreq=1):
     dataVec = survey.eval(fields)
 
     # Make the data
-    mtData = MT.Data(survey,dataVec)
+    mtData = NSEM.Data(survey,dataVec)
     # Add plots
     if plotIt:
         pass

SimPEG/Examples/__init__.py

@@ -1,28 +1,29 @@
 # Run this file to add imports.
 
 ##### AUTOIMPORTS #####
-import DC_Analytic_Dipole
-import DC_Forward_PseudoSection
 import EM_FDEM_1D_Inversion
-import EM_FDEM_Analytic_MagDipoleWholespace
+import Mesh_QuadTree_Creation
-import EM_Schenkel_Morrison_Casing
 import EM_TDEM_1D_Inversion
+import Mesh_QuadTree_FaceDiv
+import Mesh_Tensor_Creation
 import FLOW_Richards_1D_Celia1990
+import DC_Forward_PseudoSection
+import Mesh_Operators_CahnHilliard
+import Mesh_Basic_Types
 import Inversion_IRLS
 import Inversion_Linear
-import Mesh_Basic_ForwardDC
+import EM_Schenkel_Morrison_Casing
-import Mesh_Basic_PlotImage
-import Mesh_Basic_Types
-import Mesh_Operators_CahnHilliard
-import Mesh_QuadTree_Creation
-import Mesh_QuadTree_FaceDiv
-import Mesh_QuadTree_HangingNodes
-import Mesh_Tensor_Creation
-import MT_1D_ForwardAndInversion
 import MT_3D_Foward
+import Mesh_Basic_ForwardDC
+import MT_1D_ForwardAndInversion
 import Utils_surface2ind_topo
+import MT_1D_analytic_nlayer_Earth
+import EM_FDEM_Analytic_MagDipoleWholespace
+import Mesh_Basic_PlotImage
+import DC_Analytic_Dipole
+import Mesh_QuadTree_HangingNodes
 
-__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Inversion_IRLS", "Inversion_Linear", "Mesh_Basic_ForwardDC", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward", "Utils_surface2ind_topo"]
+__examples__ = ["EM_FDEM_1D_Inversion", "Mesh_QuadTree_Creation", "EM_TDEM_1D_Inversion", "Mesh_QuadTree_FaceDiv", "Mesh_Tensor_Creation", "FLOW_Richards_1D_Celia1990", "DC_Forward_PseudoSection", "Mesh_Operators_CahnHilliard", "Mesh_Basic_Types", "Inversion_IRLS", "Inversion_Linear", "EM_Schenkel_Morrison_Casing", "MT_3D_Foward", "Mesh_Basic_ForwardDC", "MT_1D_ForwardAndInversion", "Utils_surface2ind_topo", "MT_1D_analytic_nlayer_Earth", "EM_FDEM_Analytic_MagDipoleWholespace", "Mesh_Basic_PlotImage", "DC_Analytic_Dipole", "Mesh_QuadTree_HangingNodes"]
 
 ##### AUTOIMPORTS #####
 

SimPEG/MT/BaseMT.py

@@ -1,132 +0,0 @@
-from SimPEG import SolverLU as SimpegSolver, PropMaps, Utils, mkvc, sp, np
-from SimPEG.EM.FDEM.ProblemFDEM import BaseFDEMProblem
-from SurveyMT import Survey, Data
-from FieldsMT import BaseMTFields
-
-
-class BaseMTProblem(BaseFDEMProblem):
-    """
-        Base class for all Natural source problems.
-    """
-
-    def __init__(self, mesh, **kwargs):
-        BaseFDEMProblem.__init__(self, mesh, **kwargs)
-        Utils.setKwargs(self, **kwargs)
-    # Set the default pairs of the problem
-    surveyPair = Survey
-    dataPair = Data
-    fieldsPair = BaseMTFields
-
-    # Set the solver
-    Solver = SimpegSolver
-    solverOpts = {}
-
-    verbose = False
-    # Notes:
-    # Use the forward and devs from BaseFDEMProblem
-    # Might need to add more stuff here.
-
-    ## 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 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 f is None:
-           f= self.fields(m)
-        # Set current model
-        self.curModel = m
-        # Initiate the Jv object
-        Jv = self.dataPair(self.survey)
-
-        # Loop all the frequenies
-        for freq in self.survey.freqs:
-            dA_du = self.getA(freq) #
-
-            dA_duI = self.Solver(dA_du, **self.solverOpts)
-
-            for src in self.survey.getSrcByFreq(freq):
-                # 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,:]
-                # 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.
-                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 )
-                else:
-                    du_dm = dA_duI * ( -dA_dm + dRHS_dm )
-                # Calculate the projection derivatives
-                for rx in src.rxList:
-                    # Get the projection derivative
-                    # v should be of size 2*nE (for 2 polarizations)
-                    PDeriv_u = lambda t: rx.evalDeriv(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, 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 u (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
-
-        # Ensure v is a data object.
-        if not isinstance(v, self.dataPair):
-            v = self.dataPair(self.survey, v)
-
-        Jtv = np.zeros(m.size)
-
-        for freq in self.survey.freqs:
-            AT = self.getA(freq).T
-
-            ATinv = self.Solver(AT, **self.solverOpts)
-
-            for src in self.survey.getSrcByFreq(freq):
-                ftype = self._fieldType + 'Solution'
-                f_src = f[src, :]
-
-                for rx in src.rxList:
-                    # Get the adjoint evalDeriv
-                    # PTv needs to be nE,
-                    PTv = rx.evalDeriv(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, 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:
-                        du_dmT = -dA_dmT
-                    else:
-                        du_dmT = -dA_dmT + dRHS_dmT
-                    # Select the correct component
-                    # du_dmT needs to be of size nC,
-                    real_or_imag = rx.projComp
-                    if real_or_imag == 'real':
-                        Jtv +=  du_dmT.real
-                    elif real_or_imag == 'imag':
-                        Jtv +=  -du_dmT.real
-                    else:
-                        raise Exception('Must be real or imag')
-            # Clean the factorization, clear memory.
-            ATinv.clean()
-        return Jtv

SimPEG/MT/Problem1D/Probs.py

@@ -1,291 +0,0 @@
-from SimPEG.EM.Utils import omega
-from SimPEG import mkvc
-from scipy.constants import mu_0
-from SimPEG.MT.BaseMT import BaseMTProblem
-from SimPEG.MT.SurveyMT import Survey, Data
-from SimPEG.MT.FieldsMT import Fields1D_e
-from SimPEG.MT.Utils.MT1Danalytic import getEHfields
-import numpy as np
-import multiprocessing, sys, time
-
-
-class eForm_psField(BaseMTProblem):
-    """
-    A MT 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 MTproblem.
-    _fieldType = 'e_1d'
-    _eqLocs    = 'EF'
-    _sigmaPrimary = None
-
-
-    def __init__(self, mesh, **kwargs):
-        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):
-        """
-        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(BaseMTProblem):
-    """
-    A MT 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 MTproblem.
-    _fieldType = 'e'
-    _eqLocs    = 'EF'
-
-
-    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):
-        """
-            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

SimPEG/MT/Problem1D/__init__.py

@@ -1 +0,0 @@
-from Probs import eForm_TotalField, eForm_psField

SimPEG/MT/Problem2D/__init__.py

@@ -1 +0,0 @@
-pass

SimPEG/MT/Problem3D/Probs.py

@@ -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.MT.BaseMT import BaseMTProblem
-from SimPEG.MT.SurveyMT import Survey, Data
-from SimPEG.MT.FieldsMT import Fields3D_e
-import multiprocessing, sys, time
-
-
-
-class eForm_ps(BaseMTProblem):
-    """
-    A MT 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 MTproblem.
-    _fieldType = 'e'
-    _eqLocs    = 'FE'
-    fieldsPair = Fields3D_e
-    _sigmaPrimary = None
-
-    def __init__(self, mesh, **kwargs):
-        BaseMTProblem.__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
-

SimPEG/MT/Problem3D/__init__.py

@@ -1 +0,0 @@
-from Probs import eForm_ps

SimPEG/MT/Utils/__init__.py

@@ -1,4 +0,0 @@
-from MT1Dsolutions import *  # Add the names of the functions
-from MT1Danalytic import *
-from dataUtils import *
-from ediFilesUtils import *

SimPEG/MT/Utils/srcUtils.py

@@ -1,46 +0,0 @@
-import SimPEG as simpeg, numpy as np
-
-def homo1DModelSource(mesh,freq,m_back):
-    '''
-        Function that calculates and return background fields for a 3D mesh and model.
-        The calculuations use 1D field solution for a vertical slice throught model (south-western most column),
-        which is assigned at the fields everywhere for the respective polarizations.2
-
-        :param Simpeg mesh object mesh: Holds information on the discretization
-        :param float freq: The frequency to solve at
-        :param np.array m_back: Background model of conductivity to base the calculations on.
-        :rtype: numpy.ndarray (mesh.nE,2)
-        :return: eBG_bp, E fields for the background model at both polarizations.
-
-    '''
-
-    # import
-    from SimPEG.MT.Utils import get1DEfields
-    # Get a 1d solution for a halfspace background
-    mesh1d = simpeg.Mesh.TensorMesh([mesh.hz],np.array([mesh.x0[2]]))
-    # Note: Everything is using e^iwt
-    e0_1d = get1DEfields(mesh1d,mesh.r(m_back,'CC','CC','M')[0,0,:],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))
-
-    # Return the electric fields
-    eBG_bp = np.hstack((eBG_px,eBG_py))
-    return eBG_bp

SimPEG/MT/__init__.py

@@ -1,5 +0,0 @@
-import Utils
-from SurveyMT import Rx, Survey, Data
-from FieldsMT import Fields1D_e, Fields3D_e
-import Problem1D, Problem2D, Problem3D
-import SrcMT

SimPEG/Mesh/MeshIO.py

@@ -140,7 +140,6 @@ class TensorMeshIO(object):
                 vtkObj.GetCellData().AddArray(vtkDoubleArr)
             # Set the active scalar
             vtkObj.GetCellData().SetActiveScalars(models.keys()[0])
-        # vtkObj.Update()
 
         # Check the extension of the fileName
         ext = os.path.splitext(fileName)[1]
@@ -157,6 +156,50 @@ class TensorMeshIO(object):
         vtrWriteFilter.SetFileName(fileName)
         vtrWriteFilter.Update()
 
+    def _toVTRObj(mesh,models=None):
+        """
+        Makes and saves a VTK rectilinear file (vtr) for a simpeg Tensor mesh and model.
+
+        Input:
+        :param str, path to the output vtk file
+        :param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK
+        :param models, dictionary of numpy.array - Name('s) and array('s). Match number of cells
+
+        """
+        # Import
+        from vtk import vtkRectilinearGrid as rectGrid, VTK_VERSION
+        from vtk.util.numpy_support import numpy_to_vtk
+
+        # Deal with dimensionalities
+        if mesh.dim >= 1:
+            vX = mesh.vectorNx
+            xD = mesh.nNx
+            yD,zD = 1,1
+            vY, vZ = np.array([0,0])
+        if mesh.dim >= 2:
+            vY = mesh.vectorNy
+            yD = mesh.nNy
+        if mesh.dim == 3:
+            vZ = mesh.vectorNz
+            zD = mesh.nNz
+        # Use rectilinear VTK grid.
+        # Assign the spatial information.
+        vtkObj = rectGrid()
+        vtkObj.SetDimensions(xD,yD,zD)
+        vtkObj.SetXCoordinates(numpy_to_vtk(vX,deep=1))
+        vtkObj.SetYCoordinates(numpy_to_vtk(vY,deep=1))
+        vtkObj.SetZCoordinates(numpy_to_vtk(vZ,deep=1))
+
+        # Assign the model('s) to the object
+        if models is not None:
+            for item in models.iteritems():
+                # Convert numpy array
+                vtkDoubleArr = numpy_to_vtk(item[1],deep=1)
+                vtkDoubleArr.SetName(item[0])
+                vtkObj.GetCellData().AddArray(vtkDoubleArr)
+            # Set the active scalar
+            vtkObj.GetCellData().SetActiveScalars(models.keys()[0])
+        return vtkObj
 
     def readModelUBC(mesh, fileName):
         """

SimPEG/NSEM/FieldsNSEM.py

@@ -4,18 +4,21 @@ import sys
 from numpy.lib import recfunctions as recFunc
 from SimPEG.EM.Utils import omega
 
+
 ##############
 ### Fields ###
 ##############
-class BaseMTFields(Problem.Fields):
+class BaseNSEMFields(Problem.Fields):
-    """Field Storage for a MT survey."""
+    """Field Storage for a NSEM survey."""
     knownFields = {}
     dtype = complex
 
-
+###########
-class Fields1D_e(BaseMTFields):
+# 1D Fields
+###########
+class Fields1D_ePrimSec(BaseNSEMFields):
     """
-    Fields storage for the 1D MT solution.
+    Fields storage for the 1D NSEM solution.
     """
     knownFields = {'e_1dSolution':'F'}
     aliasFields = {
@@ -28,7 +31,119 @@ class Fields1D_e(BaseMTFields):
                   }
 
     def __init__(self,mesh,survey,**kwargs):
-        BaseMTFields.__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, 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)
@@ -99,7 +214,7 @@ class Fields1D_e(BaseMTFields):
         """
         Derivative of the fields object wrt u.
 
-        :param MTsrc src: MT source
+        :param NSEMsrc src: NSEM source
         :param numpy.ndarray v: random vector of f_sol.size
         This function stacks the fields derivatives appropriately
 
@@ -120,9 +235,18 @@ class Fields1D_e(BaseMTFields):
         """
         return None
 
-class Fields3D_e(BaseMTFields):
+
+###########
+# 2D Fields
+###########
+
+
+###########
+# 3D Fields
+###########
+class Fields3D_ePrimSec(BaseNSEMFields):
     """
-    Fields storage for the 3D MT solution. Labels polarizations by px and py.
+    Fields storage for the 3D NSEM solution. Labels polarizations by px and py.
 
         :param SimPEG object mesh: The solution mesh
         :param SimPEG object survey: A survey object
@@ -147,7 +271,7 @@ class Fields3D_e(BaseMTFields):
                   }
 
     def __init__(self,mesh,survey,**kwargs):
-        BaseMTFields.__init__(self,mesh,survey,**kwargs)
+        BaseNSEMFields.__init__(self,mesh,survey,**kwargs)
 
     def _e_pxPrimary(self, e_pxSolution, srcList):
         e_pxPrimary = np.zeros_like(e_pxSolution)
@@ -228,7 +352,7 @@ class Fields3D_e(BaseMTFields):
         b = (C * e_pxSolution)
         for i, src in enumerate(srcList):
             b[:,i] *= - 1./(1j*omega(src.freq))
-            # There is no magnetic source in the MT problem
+            # There is no magnetic source in the NSEM problem
             # S_m, _ = src.eval(self.survey.prob)
             # if S_m is not None:
             #     b[:,i] += 1./(1j*omega(src.freq)) * S_m
@@ -239,7 +363,7 @@ class Fields3D_e(BaseMTFields):
         b = (C * e_pySolution)
         for i, src in enumerate(srcList):
             b[:,i] *= - 1./(1j*omega(src.freq))
-            # There is no magnetic source in the MT problem
+            # There is no magnetic source in the NSEM problem
             # S_m, _ = src.eval(self.survey.prob)
             # if S_m is not None:
             #     b[:,i] += 1./(1j*omega(src.freq)) * S_m
@@ -302,7 +426,7 @@ class Fields3D_e(BaseMTFields):
         """
         Derivative of the fields object wrt u.
 
-        :param MTsrc src: MT source
+        :param NSEMsrc src: NSEM source
         :param numpy.ndarray v: random vector of f_sol.size
         This function stacks the fields derivatives appropriately
 
@@ -319,7 +443,7 @@ class Fields3D_e(BaseMTFields):
         """
         Derivative of the fields object wrt u.
 
-        :param MTsrc src: MT source
+        :param NSEMsrc src: NSEM source
         :param numpy.ndarray v: random vector of f_sol.size
         This function stacks the fields derivatives appropriately
 

SimPEG/NSEM/ProblemNSEM.py

@@ -0,0 +1,560 @@
+from SimPEG.EM.Utils.EMUtils import omega, mu_0
+from SimPEG import SolverLU as SimpegSolver, PropMaps, Utils, mkvc, sp, np
+from SimPEG.EM.FDEM.ProblemFDEM import BaseFDEMProblem
+from SurveyNSEM import Survey, Data
+from FieldsNSEM import BaseNSEMFields, Fields1D_ePrimSec, Fields3D_ePrimSec
+from SimPEG.NSEM.Utils.MT1Danalytic import getEHfields
+import time, sys
+
+class BaseNSEMProblem(BaseFDEMProblem):
+    """
+        Base class for all Natural source problems.
+    """
+
+    def __init__(self, mesh, **kwargs):
+        BaseFDEMProblem.__init__(self, mesh, **kwargs)
+        Utils.setKwargs(self, **kwargs)
+    # Set the default pairs of the problem
+    surveyPair = Survey
+    dataPair = Data
+    fieldsPair = BaseNSEMFields
+
+    # Set the solver
+    Solver = SimpegSolver
+    solverOpts = {}
+
+    verbose = False
+    # Notes:
+    # Use the forward and devs from BaseFDEMProblem
+    # Might need to add more stuff here.
+
+    ## 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 m (nC, 1) - conductive model
+            :param numpy.ndarray v (nC, 1) - random vector
+            :param NSEMfields object (optional) - NSEM fields object, if not given it is calculated
+            :rtype: NSEMdata object
+            :return: Data sensitivities wrt m
+        """
+
+        # Calculate the fields
+        if f is None:
+           f= self.fields(m)
+        # Set current model
+        self.curModel = m
+        # Initiate the Jv object
+        Jv = self.dataPair(self.survey)
+
+        # Loop all the frequenies
+        for freq in self.survey.freqs:
+            dA_du = self.getA(freq) #
+
+            dA_duI = self.Solver(dA_du, **self.solverOpts)
+
+            for src in self.survey.getSrcByFreq(freq):
+                # 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 = 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, 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 )
+                else:
+                    du_dm = dA_duI * ( -dA_dm + dRHS_dm )
+                # Calculate the projection derivatives
+                for rx in src.rxList:
+                    # Get the projection derivative
+                    # v should be of size 2*nE (for 2 polarizations)
+                    PDeriv_u = lambda t: rx.evalDeriv(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, 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 NSEMfields object f (optional) - NSEM fields object, if not given it is calculated
+            :rtype: NSEMdata object
+            :return: Data sensitivities wrt m
+        """
+
+        if f is None:
+            f = self.fields(m)
+
+        self.curModel = m
+
+        # Ensure v is a data object.
+        if not isinstance(v, self.dataPair):
+            v = self.dataPair(self.survey, v)
+
+        Jtv = np.zeros(m.size)
+
+        for freq in self.survey.freqs:
+            AT = self.getA(freq).T
+
+            ATinv = self.Solver(AT, **self.solverOpts)
+
+            for src in self.survey.getSrcByFreq(freq):
+                ftype = self._solutionType
+                f_src = f[src, :] # Need to fix this...
+
+                for rx in src.rxList:
+                    # Get the adjoint evalDeriv
+                    # PTv needs to be nE,
+                    PTv = rx.evalDeriv(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, 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:
+                        du_dmT = -dA_dmT
+                    else:
+                        du_dmT = -dA_dmT + dRHS_dmT
+                    # Select the correct component
+                    # du_dmT needs to be of size nC,
+                    real_or_imag = rx.projComp
+                    if real_or_imag == 'real':
+                        Jtv +=  du_dmT.real
+                    elif real_or_imag == 'imag':
+                        Jtv +=  -du_dmT.real
+                    else:
+                        raise Exception('Must be real or imag')
+            # 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._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

SimPEG/NSEM/SrcNSEM.py

@@ -11,9 +11,9 @@ import sys
 ###   Sources ###
 #################
 
-class BaseMTSrc(FDEMBaseSrc):
+class BaseNSEMSrc(FDEMBaseSrc):
     '''
-    Sources for the MT problem.
+    Sources for the NSEM problem.
     Use the SimPEG BaseSrc, since the source fields share properties with the transmitters.
 
     :param float freq: The frequency of the source
@@ -29,28 +29,28 @@ class BaseMTSrc(FDEMBaseSrc):
         FDEMBaseSrc.__init__(self, rxList)
 
 # 1D sources
-class polxy_1DhomotD(BaseMTSrc):
+class polxy_1DhomotD(BaseNSEMSrc):
     """
-    MT source for both polarizations (x and y) for the total Domain.
+    NSEM 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)
+        BaseNSEMSrc.__init__(self, rxList, freq)
 
 
     # TODO: need to add the  primary fields calc and source terms into the problem.
 
 # Need to implement such that it works for all dims.
-class polxy_1Dprimary(BaseMTSrc):
+class polxy_1Dprimary(BaseNSEMSrc):
     """
-    MT source for both polarizations (x and y) given a 1D primary models.
+    NSEM 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).'
         self.sigma1d = None
-        BaseMTSrc.__init__(self, rxList, freq)
+        BaseNSEMSrc.__init__(self, rxList, freq)
         # Hidden property of the ePrimary
         self._ePrimary = None
 
@@ -86,7 +86,7 @@ class polxy_1Dprimary(BaseMTSrc):
         Get the electrical field source
         """
         e_p = self.ePrimary(problem)
-        Map_sigma_p = Maps.Vertical1DMap(problem.mesh)
+        Map_sigma_p = Maps.SurjectVertical1D(problem.mesh)
         sigma_p = Map_sigma_p._transform(self.sigma1d)
         # Make mass matrix
         # Note: M(sig) - M(sig_p) = M(sig - sig_p)
@@ -128,15 +128,15 @@ class polxy_1Dprimary(BaseMTSrc):
             # v should be nC size
             return MsigmaDeriv * v
 
-class polxy_3Dprimary(BaseMTSrc):
+class polxy_3Dprimary(BaseNSEMSrc):
     """
-    MT source for both polarizations (x and y) given a 3D primary model. It assigns fields calculated from the 1D model
+    NSEM source for both polarizations (x and y) given a 3D primary model. 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).'
         self.sigmaPrimary = None
-        BaseMTSrc.__init__(self, rxList, freq)
+        BaseNSEMSrc.__init__(self, rxList, freq)
         # Hidden property of the ePrimary
         self._ePrimary = None
 
@@ -163,7 +163,7 @@ class polxy_3Dprimary(BaseMTSrc):
         Get the electrical field source
         """
         e_p = self.ePrimary(problem)
-        Map_sigma_p = Maps.Vertical1DMap(problem.mesh)
+        Map_sigma_p = Maps.SurjectVertical1D(problem.mesh)
         sigma_p = Map_sigma_p._transform(self.sigma1d)
         # Make mass matrix
         # Note: M(sig) - M(sig_p) = M(sig - sig_p)

SimPEG/NSEM/SurveyNSEM.py

@@ -4,7 +4,7 @@ from SimPEG.EM.Utils import omega
 from scipy.constants import mu_0
 from numpy.lib import recfunctions as recFunc
 from Utils import rec2ndarr
-import SrcMT
+import SrcNSEM
 import sys
 
 #################
@@ -63,9 +63,9 @@ class Rx(SimPEGsurvey.BaseRx):
         '''
         Project the fields to natural source data.
 
-            :param SrcMT src: The source of the fields to project
+            :param SrcNSEM src: The source of the fields to project
             :param SimPEG.Mesh mesh:
-            :param FieldsMT f: Natural source fields object to project
+            :param FieldsNSEM f: Natural source fields object to project
         '''
 
         ## NOTE: Assumes that e is on t
@@ -143,9 +143,9 @@ class Rx(SimPEGsurvey.BaseRx):
         """
         The derivative of the projection wrt u
 
-        :param MTsrc src: MT source
+        :param NSEMsrc src: NSEM source
         :param TensorMesh mesh: Mesh defining the topology of the problem
-        :param MTfields f: MT fields object of the source
+        :param NSEMfields f: NSEM fields object of the source
         :param numpy.ndarray v: Random vector of size
         """
 
@@ -390,12 +390,12 @@ class Rx(SimPEGsurvey.BaseRx):
 #################
 class Survey(SimPEGsurvey.BaseSurvey):
     """
-        Survey class for MT. Contains all the sources associated with the survey.
+        Survey class for NSEM. Contains all the sources associated with the survey.
 
         :param list srcList: List of sources associated with the survey
 
     """
-    srcPair = SrcMT.BaseMTSrc
+    srcPair = SrcNSEM.BaseNSEMSrc
 
     def __init__(self, srcList, **kwargs):
         # Sort these by frequency
@@ -443,7 +443,7 @@ class Survey(SimPEGsurvey.BaseSurvey):
 #################
 class Data(SimPEGsurvey.Data):
     '''
-    Data class for MTdata. Stores the data vector indexed by the survey.
+    Data class for NSEMdata. Stores the data vector indexed by the survey.
 
     :param SimPEG survey object survey:
     :param v vector of the data in order matching of the survey
@@ -461,7 +461,7 @@ class Data(SimPEGsurvey.Data):
 
     def toRecArray(self,returnType='RealImag'):
         '''
-        Function that returns a numpy.recarray for a SimpegMT impedance data object.
+        Function that returns a numpy.recarray for a SimpegNSEM impedance data object.
 
         :param str returnType: Switches between returning a rec array where the impedance is split to real and imaginary ('RealImag') or is a complex ('Complex')
 
@@ -483,7 +483,7 @@ class Data(SimPEGsurvey.Data):
                 locs = np.hstack((np.array([[0.0]]),locs))
             tArrRec = np.concatenate((src.freq*np.ones((locs.shape[0],1)),locs,np.nan*np.ones((locs.shape[0],12))),axis=1).view(dtRI)
             # np.array([(src.freq,rx.locs[0,0],rx.locs[0,1],rx.locs[0,2],np.nan ,np.nan ,np.nan ,np.nan ,np.nan ,np.nan ,np.nan ,np.nan ) for rx in src.rxList],dtype=dtRI)
-            # Get the type and the value for the DataMT object as a list
+            # Get the type and the value for the DataNSEM object as a list
             typeList = [[rx.rxType.replace('z1d','zyx'),self[src,rx]] for rx in src.rxList]
             # Insert the values to the temp array
             for nr,(key,val) in enumerate(typeList):
@@ -517,17 +517,17 @@ class Data(SimPEGsurvey.Data):
     @classmethod
     def fromRecArray(cls, recArray, srcType='primary'):
         """
-        Class method that reads in a numpy record array to MTdata object.
+        Class method that reads in a numpy record array to NSEMdata object.
 
         Only imports the impedance data.
 
         """
         if srcType=='primary':
-            src = SrcMT.polxy_1Dprimary
+            src = SrcNSEM.polxy_1Dprimary
         elif srcType=='total':
-            src = SrcMT.polxy_1DhomotD
+            src = SrcNSEM.polxy_1DhomotD
         else:
-            raise NotImplementedError('{:s} is not a valid source type for MTdata')
+            raise NotImplementedError('{:s} is not a valid source type for NSEMdata')
 
         # Find all the frequencies in recArray
         uniFreq = np.unique(recArray['freq'])

SimPEG/NSEM/Utils/MT1Danalytic.py

@@ -3,7 +3,7 @@
 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):
+def getEHfields(m1d,sigma,freq,zd,scaleUD=True,scaleValue=1):
     '''Analytic solution for MT 1D layered earth. Returns E and H fields.
 
     :param SimPEG.mesh, object m1d: Mesh object with the 1D spatial information.
@@ -12,7 +12,7 @@ def getEHfields(m1d,sigma,freq,zd,scaleUD=True):
     :param numpy array, vector zd: location to calculate EH fields at
     :param bollean, scaleUD: scales the output to be 1 at the top, increases numeracal stability.
 
-    Assumes a halfspace with the same conductive as the last cell below.
+    Assumes a halfspace with the same conductive as the deepest cell.
 
     '''
     # Note add an error check for the mesh and sigma are the same size.
@@ -29,7 +29,7 @@ def getEHfields(m1d,sigma,freq,zd,scaleUD=True):
 
     # Initiate the propagation matrix, in the order down up.
     UDp = np.zeros((2,m1d.nC+1),dtype=complex)
-    UDp[1,0] = 1. # Set the wave amplitude as 1 into the half-space at the bottom of the mesh
+    UDp[1,0] = scaleValue # Set the wave amplitude as 1 into the half-space at the bottom of the mesh
     # Loop over all the layers, starting at the bottom layer
     for lnr, h in enumerate(m1d.hx): # lnr-number of layer, h-thickness of the layer
         # Calculate
@@ -38,9 +38,9 @@ def getEHfields(m1d,sigma,freq,zd,scaleUD=True):
         # Build the propagation matrix
 
         # Convert fields to down/up going components in layer below current layer
-        Pj1 = np.array([[1,1],[yp1,-yp1]])
+        Pj1 = np.array([[1,1],[yp1,-yp1]],dtype=complex)
         # Convert fields to down/up going components in current layer
-        Pjinv = 1./2*np.array([[1,zp],[1,-zp]])
+        Pjinv = 1./2*np.array([[1,zp],[1,-zp]],dtype=complex)
         # Propagate down and up components through the current layer
         elamh = np.array([[np.exp(-1j*k[lnr+1]*h),0],[0,np.exp(1j*k[lnr+1]*h)]])
 
@@ -48,7 +48,14 @@ def getEHfields(m1d,sigma,freq,zd,scaleUD=True):
         UDp[:,lnr+1] = elamh.dot(Pjinv.dot(Pj1)).dot(UDp[:,lnr])
 
         if scaleUD:
-             UDp[:,lnr+1::-1] = UDp[:,lnr+1::-1]/UDp[1,lnr+1]
+            # Scale the values such that 1 at the top
+            scaleVal = UDp[:,lnr+1::-1]/UDp[1,lnr+1]
+            if np.any(np.isnan(scaleVal)):
+                # If there is a nan (thickness very great), rebuild the move up cell
+                scaleVal = np.zeros_like(UDp[:,lnr+1::-1],dtype=complex)
+                scaleVal[1,0] = scaleValue
+
+            UDp[:,lnr+1::-1] = scaleVal
 
     # Calculate the fields
     Ed = np.empty((zd.size,),dtype=complex)

SimPEG/NSEM/Utils/__init__.py

@@ -0,0 +1,5 @@
+from MT1Dsolutions import get1DEfields  # Add the names of the functions
+from MT1Danalytic import getEHfields, getImpedance
+from dataUtils import *
+from ediFilesUtils import *
+from testUtils import *

SimPEG/NSEM/Utils/dataUtils.py

@@ -5,25 +5,25 @@ import numpy.lib.recfunctions as recFunc
 from scipy.constants import mu_0
 from scipy import interpolate as sciint
 
-def getAppRes(MTdata):
+def getAppRes(NSEMdata):
     # Make impedance
     zList = []
-    for src in MTdata.survey.srcList:
+    for src in NSEMdata.survey.srcList:
         zc = [src.freq]
         for rx in src.rxList:
             if 'i' in rx.rxType:
                 m=1j
             else:
                 m = 1
-            zc.append(m*MTdata[src,rx])
+            zc.append(m*NSEMdata[src,rx])
         zList.append(zc)
     return [appResPhs(zList[i][0],np.sum(zList[i][1:3])) for i in np.arange(len(zList))]
 
-def rotateData(MTdata,rotAngle):
+def rotateData(NSEMdata,rotAngle):
     '''
     Function that rotates clockwist by rotAngle (- negative for a counter-clockwise rotation)
     '''
-    recData = MTdata.toRecArray('Complex')
+    recData = NSEMdata.toRecArray('Complex')
     impData = rec2ndarr(recData[['zxx','zxy','zyx','zyy']],complex)
     # Make the rotation matrix
     # c,s,zxx,zxy,zyx,zyy = sympy.symbols('c,s,zxx,zxy,zyx,zyy')
@@ -40,8 +40,8 @@ def rotateData(MTdata,rotAngle):
     for nr,comp in enumerate(['zxx','zxy','zyx','zyy']):
         outRec[comp] = rotData[:,nr]
 
-    from SimPEG import MT
+    from SimPEG import NSEM
-    return MT.Data.fromRecArray(outRec)
+    return NSEM.Data.fromRecArray(outRec)
 
 
 def appResPhs(freq,z):
@@ -57,10 +57,10 @@ def rec2ndarr(x,dt=float):
     return x.view((dt, len(x.dtype.names)))
 
 def makeAnalyticSolution(mesh,model,elev,freqs):
-    from SimPEG import MT
+    from SimPEG import NSEM
     data1D = []
     for freq in freqs:
-        anaEd, anaEu, anaHd, anaHu = MT.Utils.MT1Danalytic.getEHfields(mesh,model,freq,elev)
+        anaEd, anaEu, anaHd, anaHu = NSEM.Utils.MT1Danalytic.getEHfields(mesh,model,freq,elev)
         anaE = anaEd+anaEu
         anaH = anaHd+anaHu
 
@@ -71,7 +71,7 @@ def makeAnalyticSolution(mesh,model,elev,freqs):
     return dataRec
 
 def plotMT1DModelData(problem,models,symList=None):
-    from SimPEG import MT
+    from SimPEG import NSEM
     # Setup the figure
     fontSize = 15
 
@@ -79,7 +79,7 @@ def plotMT1DModelData(problem,models,symList=None):
     axM = fig.add_axes([0.075,.1,.25,.875])
     axM.set_xlabel('Resistivity [Ohm*m]',fontsize=fontSize)
     axM.set_xlim(1e-1,1e5)
-    axM.set_ylim(-10000,5000)
+    # axM.set_ylim(-10000,5000)
     axM.set_ylabel('Depth [km]',fontsize=fontSize)
     axR = fig.add_axes([0.42,.575,.5,.4])
     axR.set_xscale('log')
@@ -132,38 +132,94 @@ def plotMT1DModelData(problem,models,symList=None):
     freq = simpeg.mkvc(data1D['freq'],2)
     res, phs = appResPhs(freq,allData)
 
-    stdCol = 'gray'
+    if False:
-    axRtw = axR.twinx()
+        stdCol = 'gray'
-    axRtw.set_ylabel('Std of log10',color=stdCol)
+        axRtw = axR.twinx()
-    [(t.set_color(stdCol), t.set_rotation(-45)) for t in axRtw.get_yticklabels()]
+        axRtw.set_ylabel('Std of log10',color=stdCol)
-    axPtw = axP.twinx()
+        [(t.set_color(stdCol), t.set_rotation(-45)) for t in axRtw.get_yticklabels()]
-    axPtw.set_ylabel('Std ',color=stdCol)
+        axPtw = axP.twinx()
-    [t.set_color(stdCol) for t in axPtw.get_yticklabels()]
+        axPtw.set_ylabel('Std ',color=stdCol)
-    axRtw.plot(freq, np.std(np.log10(res),1),'--',color=stdCol)
+        [t.set_color(stdCol) for t in axPtw.get_yticklabels()]
-    axPtw.plot(freq, np.std(phs,1),'--',color=stdCol)
+        axRtw.plot(freq, np.std(np.log10(res),1),'--',color=stdCol)
+        axPtw.plot(freq, np.std(phs,1),'--',color=stdCol)
 
     # Fix labels and ticks
 
-    yMtick = [l/1000 for l in axM.get_yticks().tolist()]
+    # yMtick = [l/1000 for l in axM.get_yticks().tolist()]
-    axM.set_yticklabels(yMtick)
+    # axM.set_yticklabels(yMtick)
     [ l.set_rotation(90) for l in axM.get_yticklabels()]
     [ l.set_rotation(90) for l in axR.get_yticklabels()]
-    [(t.set_color(stdCol), t.set_rotation(-45)) for t in axRtw.get_yticklabels()]
+    # [(t.set_color(stdCol), t.set_rotation(-45)) for t in axRtw.get_yticklabels()]
-    [t.set_color(stdCol) for t in axPtw.get_yticklabels()]
+    # [t.set_color(stdCol) for t in axPtw.get_yticklabels()]
     for ax in [axM,axR,axP]:
         ax.xaxis.set_tick_params(labelsize=fontSize)
         ax.yaxis.set_tick_params(labelsize=fontSize)
     return fig
 
+def plotImpAppRes(dataArrays,plotLoc,textStr=[]):
+    ''' Plots amplitude impedance and phase'''
+    # fig = plt.figure(1,(7, 7))
+    import plotDataTypes as pDt
+    # axes = ImageGrid(fig, (0.05,0.05,0.875,0.875),nrows_ncols = (2, 2),axes_pad = 0.25,add_all=True,share_all=True,label_mode = "L")
+    # Make the figure and axes
+    fig,axT=plt.subplots(2,2,sharex=True)
+    axes = axT.ravel()
+    fig.set_size_inches((13.5,7.0))
+    fig.suptitle('{:s}\nStation at: {:.1f}x ; {:.1f}y'.format(textStr,plotLoc[0],plotLoc[1]))
+    # Have to deal with axes
+    # Set log
+    for ax in axes.ravel():
+        ax.set_xscale('log')
+
+    axes[0].invert_xaxis()
+    axes[0].set_yscale('log')
+    axes[2].set_yscale('log')
+    # Set labels
+    axes[2].set_xlabel('Frequency [Hz]')
+    axes[3].set_xlabel('Frequency [Hz]')
+    axes[0].set_ylabel('Apperent resistivity [Ohm m]')
+    axes[1].set_ylabel('Apperent phase [degrees]')
+    axes[1].set_ylim(-180,180)
+    axes[2].set_ylabel('Impedance amplitude [V/A]')
+    axes[3].set_ylim(-180,180)
+    axes[3].set_ylabel('Impedance angle [degrees]')
+
+
+    # Plot the data
+    for nr,dataArray in enumerate(dataArrays):
+        if nr==1:
+            parSym = '*'
+        else:
+            parSym = 's'
+        # app res
+        pDt.plotIsoStaImpedance(axes[0],plotLoc,dataArray,'zxy',par='res',pSym=parSym)
+        pDt.plotIsoStaImpedance(axes[0],plotLoc,dataArray,'zyx',par='res',pSym=parSym)
+        # app phs
+        pDt.plotIsoStaImpedance(axes[1],plotLoc,dataArray,'zxy',par='phs',pSym=parSym)
+        pDt.plotIsoStaImpedance(axes[1],plotLoc,dataArray,'zyx',par='phs',pSym=parSym)
+        # imp abs
+        pDt.plotIsoStaImpedance(axes[2],plotLoc,dataArray,'zxx',par='abs',pSym=parSym)
+        pDt.plotIsoStaImpedance(axes[2],plotLoc,dataArray,'zxy',par='abs',pSym=parSym)
+        pDt.plotIsoStaImpedance(axes[2],plotLoc,dataArray,'zyx',par='abs',pSym=parSym)
+        pDt.plotIsoStaImpedance(axes[2],plotLoc,dataArray,'zyy',par='abs',pSym=parSym)
+            # imp abs
+        pDt.plotIsoStaImpedance(axes[3],plotLoc,dataArray,'zxx',par='phs',pSym=parSym)
+        pDt.plotIsoStaImpedance(axes[3],plotLoc,dataArray,'zxy',par='phs',pSym=parSym)
+        pDt.plotIsoStaImpedance(axes[3],plotLoc,dataArray,'zyx',par='phs',pSym=parSym)
+        pDt.plotIsoStaImpedance(axes[3],plotLoc,dataArray,'zyy',par='phs',pSym=parSym)
+
+    return fig,axes
+
+
 def printTime():
     import time
     print time.strftime("%a, %d %b %Y %H:%M:%S +0000", time.localtime())
 
-def convert3Dto1Dobject(MTdata,rxType3D='zyx'):
+def convert3Dto1Dobject(NSEMdata,rxType3D='zyx'):
-    from SimPEG import MT
+    from SimPEG import NSEM
     # Find the unique locations
     # Need to find the locations
-    recDataTemp = MTdata.toRecArray()
+    recDataTemp = NSEMdata.toRecArray()
     # Check if survey.std has been assigned.
     ## NEED TO: write this...
     # Calculte and add the DET of the tensor to the recArray
@@ -185,24 +241,24 @@ def convert3Dto1Dobject(MTdata,rxType3D='zyx'):
         # Make the receiver list
         rx1DList = []
         for rxType in ['z1dr','z1di']:
-            rx1DList.append(MT.Rx(simpeg.mkvc(loc,2).T,rxType))
+            rx1DList.append(NSEM.Rx(simpeg.mkvc(loc,2).T,rxType))
         # Source list
         locrecData = recData[np.sqrt(np.sum( (rec2ndarr(recData[['x','y','z']]).data - loc )**2,axis=1)) < 1e-5]
         dat1DList = []
         src1DList = []
         for freq in locrecData['freq']:
-            src1DList.append(MT.SrcMT.src_polxy_1Dprimary(rx1DList,freq))
+            src1DList.append(NSEM.SrcNSEM.src_polxy_1Dprimary(rx1DList,freq))
             for comp  in ['r','i']:
                 dat1DList.append( corr * locrecData[rxType3D+comp][locrecData['freq']== freq].data )
 
         # Make the survey
-        sur1D = MT.Survey(src1DList)
+        sur1D = NSEM.Survey(src1DList)
 
         # Make the data
         dataVec = np.hstack(dat1DList)
-        dat1D = MT.Data(sur1D,dataVec)
+        dat1D = NSEM.Data(sur1D,dataVec)
         sur1D.dobs = dataVec
-        # Need to take MTdata.survey.std and split it as well.
+        # Need to take NSEMdata.survey.std and split it as well.
         std=0.05
         sur1D.std =  np.abs(sur1D.dobs*std) #+ 0.01*np.linalg.norm(sur1D.dobs)
         mtData1DList.append(dat1D)
@@ -210,29 +266,29 @@ def convert3Dto1Dobject(MTdata,rxType3D='zyx'):
     # Return the the list of data.
     return mtData1DList
 
-def resampleMTdataAtFreq(MTdata,freqs):
+def resampleNSEMdataAtFreq(NSEMdata,freqs):
     """
-    Function to resample MTdata at set of frequencies
+    Function to resample NSEMdata at set of frequencies
 
     """
-    from SimPEG import MT
+    from SimPEG import NSEM
     # Make a rec array
-    MTrec = MTdata.toRecArray().data
+    NSEMrec = NSEMdata.toRecArray().data
 
     # Find unique locations
-    uniLoc = np.unique(MTrec[['x','y','z']])
+    uniLoc = np.unique(NSEMrec[['x','y','z']])
-    uniFreq = MTdata.survey.freqs
+    uniFreq = NSEMdata.survey.freqs
     # Get the comps
-    dNames = MTrec.dtype
+    dNames = NSEMrec.dtype
 
     # Loop over all the locations and interpolate
     for loc in uniLoc:
         # Find the index of the station
-        ind = np.sqrt(np.sum((rec2ndarr(MTrec[['x','y','z']]) - rec2ndarr(loc))**2,axis=1)) < 1. # Find dist of 1 m accuracy
+        ind = np.sqrt(np.sum((rec2ndarr(NSEMrec[['x','y','z']]) - rec2ndarr(loc))**2,axis=1)) < 1. # Find dist of 1 m accuracy
         # Make a temporary recArray and interpolate all the components
         tArrRec = np.concatenate((simpeg.mkvc(freqs,2),np.ones((len(freqs),1))*rec2ndarr(loc),np.nan*np.ones((len(freqs),12))),axis=1).view(dNames)
         for comp in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi','tzxr','tzxi','tzyr','tzyi']:
-            int1d = sciint.interp1d(MTrec[ind]['freq'],MTrec[ind][comp],bounds_error=False)
+            int1d = sciint.interp1d(NSEMrec[ind]['freq'],NSEMrec[ind][comp],bounds_error=False)
             tArrRec[comp] = simpeg.mkvc(int1d(freqs),2)
 
         # Join together
@@ -241,5 +297,5 @@ def resampleMTdataAtFreq(MTdata,freqs):
         except NameError as e:
             outRecArr = tArrRec
 
-    # Make the MTdata and return
+    # Make the NSEMdata and return
-    return MT.Data.fromRecArray(outRecArr)
+    return NSEM.Data.fromRecArray(outRecArr)

SimPEG/NSEM/Utils/ediFilesUtils.py

@@ -2,7 +2,7 @@
 from SimPEG import mkvc
 from scipy.constants import mu_0
 from numpy.lib import recfunctions as recFunc
-from SimPEG.MT.Utils.dataUtils import rec2ndarr
+from SimPEG.NSEM.Utils.dataUtils import rec2ndarr
 
 # Import modules
 import numpy as np

SimPEG/NSEM/Utils/sourceUtils.py

@@ -12,7 +12,7 @@ def homo1DModelSource(mesh,freq,sigma_1d):
 
     '''
     # import
-    from SimPEG.MT.Utils import get1DEfields
+    from SimPEG.NSEM.Utils import get1DEfields
     # Get a 1d solution for a halfspace background
     if mesh.dim == 1:
         mesh1d = mesh
@@ -77,7 +77,7 @@ def analytic1DModelSource(mesh,freq,sigma_1d):
 
     '''
     # import
-    from SimPEG.MT.Utils import getEHfields
+    from SimPEG.NSEM.Utils import getEHfields
     # Get a 1d solution for a halfspace background
     if mesh.dim == 1:
         mesh1d = mesh

SimPEG/NSEM/Utils/testUtils.py

@@ -0,0 +1,198 @@
+import unittest
+import sys
+from scipy.constants import mu_0
+import SimPEG as simpeg
+
+from SimPEG.Utils import meshTensor
+import numpy as np
+
+np.random.seed(1100)
+# Define the tolerances
+TOLr = 5e-2
+TOLp = 5e-2
+
+
+def getAppResPhs(NSEMdata):
+    # Make impedance
+    from SimPEG.NSEM.Utils import appResPhs
+    zList = []
+    for src in NSEMdata.survey.srcList:
+        zc = [src.freq]
+        for rx in src.rxList:
+            if 'i' in rx.rxType:
+                m=1j
+            else:
+                m = 1
+            zc.append(m*NSEMdata[src,rx])
+        zList.append(zc)
+    return [appResPhs(zList[i][0],np.sum(zList[i][1:3])) for i in np.arange(len(zList))]
+
+
+def setup1DSurvey(sigmaHalf,tD=True,structure=False):
+    from SimPEG import NSEM
+    # Frequency
+    nFreq = 33
+    freqs = np.logspace(3,-3,nFreq)
+    # Make the mesh
+    ct = 5
+    air = meshTensor([(ct,25,1.3)])
+    # coreT0 = meshTensor([(ct,15,1.2)])
+    # coreT1 = np.kron(meshTensor([(coreT0[-1],15,1.3)]),np.ones((7,)))
+    core = np.concatenate( (  np.kron(meshTensor([(ct,15,-1.2)]),np.ones((10,))) , meshTensor([(ct,20)]) ) )
+    bot = meshTensor([(core[0],20,-1.3)])
+    x0 = -np.array([np.sum(np.concatenate((core,bot)))])
+    m1d = simpeg.Mesh.TensorMesh([np.concatenate((bot,core,air))], x0=x0)
+    # Make the model
+    sigma = np.zeros(m1d.nC) + sigmaHalf
+    sigma[m1d.gridCC > 0 ] = 1e-8
+    sigmaBack = sigma.copy()
+    # Add structure
+    if structure:
+        shallow = (m1d.gridCC < -200) * (m1d.gridCC > -600)
+        deep = (m1d.gridCC < -3000) * (m1d.gridCC > -5000)
+        sigma[shallow] = 1
+        sigma[deep] = 0.1
+
+    rxList = []
+    for rxType in ['z1dr','z1di']:
+        rxList.append(NSEM.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType))
+    # Source list
+    srcList =[]
+    if tD:
+        for freq in freqs:
+            srcList.append(NSEM.SrcNSEM.polxy_1DhomotD(rxList,freq))
+    else:
+        for freq in freqs:
+            srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq))
+
+    survey = NSEM.Survey(srcList)
+    return survey, sigma, m1d
+
+
+def setupSimpegNSEM_ePrimSec(inputSetup,comp='Imp',singleFreq=False,expMap=True):
+    from SimPEG import NSEM
+
+    M,freqs,sig,sigBG,rx_loc = inputSetup
+    # Make a receiver list
+    rxList = []
+    if comp == 'All':
+        for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi','tzxr','tzxi','tzyr','tzyi']:
+            rxList.append(NSEM.Rx(rx_loc,rxType))
+    elif comp == 'Imp':
+        for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi']:
+            rxList.append(NSEM.Rx(rx_loc,rxType))
+    elif comp == 'Tip':
+        for rxType in ['tzxr','tzxi','tzyr','tzyi']:
+            rxList.append(NSEM.Rx(rx_loc,rxType))
+    else:
+        rxList.append(NSEM.Rx(rx_loc,comp))
+    # Source list
+    srcList =[]
+
+    if singleFreq:
+        srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,singleFreq))
+    else:
+        for freq in freqs:
+            srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq))
+    # Survey NSEM
+    survey = NSEM.Survey(srcList)
+
+    ## Setup the problem object
+    sigma1d = M.r(sigBG,'CC','CC','M')[0,0,:]
+    if expMap:
+        problem = NSEM.Problem3D_ePrimSec(M,sigmaPrimary= np.log(sigma1d) )
+        problem.mapping = simpeg.Maps.ExpMap(problem.mesh)
+        problem.curModel = np.log(sig)
+    else:
+        problem = NSEM.Problem3D_ePrimSec(M,sigmaPrimary= sigma1d)
+        problem.curModel = sig
+    problem.pair(survey)
+    problem.verbose = False
+    try:
+        from pymatsolver import MumpsSolver
+        problem.Solver = MumpsSolver
+    except:
+        pass
+
+    return (survey, problem)
+
+def getInputs():
+    """
+    Function that returns Mesh, freqs, rx_loc, elev.
+    """
+    # Make a mesh
+    # M = simpeg.Mesh.TensorMesh([[(100,5,-1.5),(100.,10),(100,5,1.5)],[(100,5,-1.5),(100.,10),(100,5,1.5)],[(100,5,1.6),(100.,10),(100,3,2)]], x0=['C','C',-3529.5360])
+    # M = simpeg.Mesh.TensorMesh([[(1000,6,-1.5),(1000.,6),(1000,6,1.5)],[(1000,6,-1.5),(1000.,2),(1000,6,1.5)],[(1000,6,-1.3),(1000.,6),(1000,6,1.3)]], x0=['C','C','C'])# Setup the model
+    M = simpeg.Mesh.TensorMesh([[(200,6,-1.5),(200.,4),(200,6,1.5)],[(200,6,-1.5),(200.,4),(200,6,1.5)],[(200,8,-1.5),(200.,8),(200,8,1.5)]], x0=['C','C','C'])# Setup the model
+    # Set the frequencies
+    freqs = np.logspace(1,-3,5)
+    elev = 0
+
+    ## Setup the the survey object
+    # Receiver locations
+    rx_x, rx_y = np.meshgrid(np.arange(-350,350,200),np.arange(-350,350,200))
+    rx_loc = np.hstack((simpeg.Utils.mkvc(rx_x,2),simpeg.Utils.mkvc(rx_y,2),elev+np.zeros((np.prod(rx_x.shape),1))))
+
+    return M, freqs, rx_loc, elev
+
+def random(conds):
+    ''' Returns a halfspace model based on the inputs'''
+    M, freqs, rx_loc, elev = getInputs()
+
+    # Backround
+    sigBG = np.ones(M.nC)*conds
+    # Add randomness to the model (10% of the value).
+    sig = np.exp( np.log(sigBG) + np.random.randn(M.nC)*(conds)*1e-1 )
+
+    return (M, freqs, sig, sigBG, rx_loc)
+
+def halfSpace(conds):
+    ''' Returns a halfspace model based on the inputs'''
+    M, freqs, rx_loc, elev = getInputs()
+
+    # Model
+    ccM = M.gridCC
+    # conds = [1e-2]
+    groundInd = ccM[:,2] < elev
+    sig = np.zeros(M.nC) + 1e-8
+    sig[groundInd] = conds
+    # Set the background, not the same as the model
+    sigBG = np.zeros(M.nC) + 1e-8
+    sigBG[groundInd] = conds
+
+    return (M, freqs, sig, sigBG, rx_loc)
+
+def blockInhalfSpace(conds):
+    ''' Returns a halfspace model based on the inputs'''
+    M, freqs, rx_loc, elev = getInputs()
+
+    # Model
+    ccM = M.gridCC
+    # conds = [1e-2]
+    groundInd = ccM[:,2] < elev
+    sig = simpeg.Utils.ModelBuilder.defineBlock(M.gridCC,np.array([-1000,-1000,-1500]),np.array([1000,1000,-1000]),conds)
+    sig[~groundInd] = 1e-8
+    # Set the background, not the same as the model
+    sigBG = np.zeros(M.nC) + 1e-8
+    sigBG[groundInd] = conds[1]
+
+    return (M, freqs, sig, sigBG, rx_loc)
+
+def twoLayer(conds):
+    ''' Returns a 2 layer model based on the conductivity values given'''
+    M, freqs, rx_loc, elev = getInputs()
+
+    # Model
+    ccM = M.gridCC
+    groundInd = ccM[:,2] < elev
+    botInd = ccM[:,2] < -3000
+    sig = np.zeros(M.nC) + 1e-8
+    sig[groundInd] = conds[1]
+    sig[botInd] = conds[0]
+    # Set the background, not the same as the model
+    sigBG = np.zeros(M.nC) + 1e-8
+    sigBG[groundInd] = conds[1]
+
+
+    return (M, freqs, sig, sigBG, rx_loc)
+

SimPEG/NSEM/__init__.py

@@ -0,0 +1,5 @@
+import Utils
+from SurveyNSEM import Rx, Survey, Data
+from FieldsNSEM import Fields1D_ePrimSec, Fields3D_ePrimSec
+from ProblemNSEM import Problem1D_ePrimSec, Problem3D_ePrimSec
+import SrcNSEM

docs/examples/MT_1D_analytic_nlayer_Earth.rst

@@ -0,0 +1,21 @@
+.. _examples_MT_1D_analytic_nlayer_Earth:
+
+.. --------------------------------- ..
+..                                   ..
+..    THIS FILE IS AUTO GENEREATED   ..
+..                                   ..
+..    SimPEG/Examples/__init__.py    ..
+..                                   ..
+.. --------------------------------- ..
+
+MT 1D analytic nlayer Earth
+===========================
+
+.. plot::
+
+    from SimPEG import Examples
+    Examples.MT_1D_analytic_nlayer_Earth.run()
+
+.. literalinclude:: ../../SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py
+    :language: python
+    :linenos:

tests/mt/forward/test_AnalyticFunctionVsAppResPhs.py

@@ -1,13 +1,10 @@
 import unittest
 from SimPEG import *
-from SimPEG import MT
+from SimPEG import NSEM
 
 TOL = 1e-6
 
-def appResPhs(freq,z):
+
-    app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
-    app_phs = np.arctan2(-z.imag,z.real)*(180/np.pi)
-    return app_res, app_phs
 
 def appResNorm(sigmaHalf):
     nFreq = 26
@@ -20,12 +17,12 @@ def appResNorm(sigmaHalf):
     freqs = np.logspace(4,-4,nFreq)
     Z = []
     for freq in freqs:
-        Ed, Eu, Hd, Hu = MT.Utils.getEHfields(m1d,sigma,freq,np.array([200]))
+        Ed, Eu, Hd, Hu = NSEM.Utils.getEHfields(m1d,sigma,freq,np.array([200]))
         Z.append((Ed + Eu)/(Hd + Hu))
 
     Zarr = np.concatenate(Z)
 
-    app_r, app_p = appResPhs(freqs,Zarr)
+    app_r, app_p = NSEM.Utils.appResPhs(freqs,Zarr)
 
     return np.linalg.norm(np.abs(app_r - np.ones(nFreq)/sigmaHalf)) / np.log10(sigmaHalf)
 

tests/mt/forward/test_Problem1D_AnalyticVsNumeric.py

@@ -1,6 +1,6 @@
 import unittest
 import SimPEG as simpeg
-from SimPEG import MT
+from SimPEG import NSEM
 from SimPEG.Utils import meshTensor
 import numpy as np
 # Define the tolerances
@@ -8,69 +8,30 @@ TOLr = 5e-2
 TOLp = 5e-2
 
 
-def setupSurvey(sigmaHalf,tD=True):
+def getAppResPhs(NSEMdata):
-
-    # Frequency
-    nFreq = 33
-    freqs = np.logspace(3,-3,nFreq)
-    # Make the mesh
-    ct = 5
-    air = meshTensor([(ct,25,1.3)])
-    # coreT0 = meshTensor([(ct,15,1.2)])
-    # coreT1 = np.kron(meshTensor([(coreT0[-1],15,1.3)]),np.ones((7,)))
-    core = np.concatenate( (  np.kron(meshTensor([(ct,15,-1.2)]),np.ones((10,))) , meshTensor([(ct,20)]) ) )
-    bot = meshTensor([(core[0],15,-1.3)])
-    x0 = -np.array([np.sum(np.concatenate((core,bot)))])
-    m1d = simpeg.Mesh.TensorMesh([np.concatenate((bot,core,air))], x0=x0)
-    # Make the model
-    sigma = np.zeros(m1d.nC) + sigmaHalf
-    sigma[m1d.gridCC > 0 ] = 1e-8
-    sigmaBack = sigma.copy()
-    # Add structure
-    shallow = (m1d.gridCC < -200) * (m1d.gridCC > -600)
-    deep = (m1d.gridCC < -3000) * (m1d.gridCC > -5000)
-    sigma[shallow] = 1
-    sigma[deep] = 0.1
-
-    rxList = []
-    for rxType in ['z1dr','z1di']:
-        rxList.append(MT.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType))
-    # Source list
-    srcList =[]
-    if tD:
-        for freq in freqs:
-            srcList.append(MT.SrcMT.polxy_1DhomotD(rxList,freq))
-    else:
-        for freq in freqs:
-            srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
-
-    survey = MT.Survey(srcList)
-    return survey, sigma, m1d
-
-def getAppResPhs(MTdata):
     # Make impedance
     def appResPhs(freq,z):
         app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
         app_phs = np.arctan2(z.imag,z.real)*(180/np.pi)
         return app_res, app_phs
     zList = []
-    for src in MTdata.survey.srcList:
+    for src in NSEMdata.survey.srcList:
         zc = [src.freq]
         for rx in src.rxList:
             if 'i' in rx.rxType:
                 m=1j
             else:
                 m = 1
-            zc.append(m*MTdata[src,rx])
+            zc.append(m*NSEMdata[src,rx])
         zList.append(zc)
     return [appResPhs(zList[i][0],np.sum(zList[i][1:3])) for i in np.arange(len(zList))]
 
 def calculateAnalyticSolution(srcList,mesh,model):
-    surveyAna = MT.Survey(srcList)
+    surveyAna = NSEM.Survey(srcList)
-    data1D = MT.Data(surveyAna)
+    data1D = NSEM.Data(surveyAna)
     for src in surveyAna.srcList:
         elev = src.rxList[0].locs[0]
-        anaEd, anaEu, anaHd, anaHu = MT.Utils.MT1Danalytic.getEHfields(mesh,model,src.freq,elev)
+        anaEd, anaEu, anaHd, anaHu = NSEM.Utils.MT1Danalytic.getEHfields(mesh,model,src.freq,elev)
         anaE = anaEd+anaEu
         anaH = anaHd+anaHu
         # Scale the solution
@@ -86,12 +47,12 @@ def dataMis_AnalyticTotalDomain(sigmaHalf):
     # Make the survey
 
     # Total domain solution
-    surveyTD, sigma, mesh = setupSurvey(sigmaHalf)
+    surveyTD, sigma, mesh = NSEM.Utils.testUtils.setup1DSurvey(sigmaHalf)
-    problemTD = MT.Problem1D.eForm_TotalField(mesh)
+    problemTD = NSEM.Problem1D_eTotal(mesh) # This not fully implemented
     problemTD.pair(surveyTD)
     # Analytic data
     dataAnaObj = calculateAnalyticSolution(surveyTD.srcList,mesh,sigma)
-    # dataTDObj = MT.DataMT.DataMT(surveyTD, surveyTD.dpred(sigma))
+    # dataTDObj = NSEM.DataNSEM.DataNSEM(surveyTD, surveyTD.dpred(sigma))
     dataTD = surveyTD.dpred(sigma)
     dataAna = simpeg.mkvc(dataAnaObj)
     return np.all((dataTD - dataAna)/dataAna < 2.)
@@ -108,16 +69,16 @@ def dataMis_AnalyticPrimarySecondary(sigmaHalf):
 
     # Make the survey
     # Primary secondary
-    surveyPS, sigmaPS, mesh = setupSurvey(sigmaHalf,tD=False)
+    survey, sigma, mesh = NSEM.Utils.testUtils.setup1DSurvey(sigmaHalf,False,structure=True)
-    problemPS = MT.Problem1D.eForm_psField(mesh)
-    problemPS.sigmaPrimary = sigmaPS
-    problemPS.pair(surveyPS)
     # Analytic data
-    dataAnaObj = calculateAnalyticSolution(surveyPS.srcList,mesh,sigmaPS)
+    problem = NSEM.Problem1D_ePrimSec(mesh, sigmaPrimary = sigma)
+    problem.pair(survey)
 
-    dataPS = surveyPS.dpred(sigmaPS)
+    dataAnaObj = calculateAnalyticSolution(survey.srcList,mesh,sigma)
+
+    data = survey.dpred(sigma)
     dataAna = simpeg.mkvc(dataAnaObj)
-    return np.all((dataPS - dataAna)/dataAna < 2.)
+    return np.all((data - dataAna)/dataAna < 2.)
 
 
 

tests/mt/forward/test_Problem1D_VsAnalyticHalfspace.py

@@ -0,0 +1,102 @@
+import unittest
+import SimPEG as simpeg
+from SimPEG import NSEM
+from SimPEG.Utils import meshTensor
+import numpy as np
+# Define the tolerances
+TOLr = 5e-1
+TOLp = 5e-1
+
+
+
+def appRes_TotalFieldNorm(sigmaHalf):
+
+    # Make the survey
+    survey, sigma, mesh = NSEM.Utils.testUtils.setup1DSurvey(sigmaHalf)
+    problem = NSEM.Problem1D_eTotal(mesh)
+    problem.pair(survey)
+
+    # Get the fields
+    fields = problem.fields(sigma)
+
+    # Project the data
+    data = survey.eval(fields)
+
+    # Calculate the app res and phs
+    app_r = np.array(NSEM.Utils.testUtils.getAppResPhs(data))[:,0]
+
+    return np.linalg.norm(np.abs(np.log(app_r) - np.log(np.ones(survey.nFreq)/sigmaHalf))*np.log(sigmaHalf))
+
+def appPhs_TotalFieldNorm(sigmaHalf):
+
+    # Make the survey
+    survey, sigma, mesh = NSEM.Utils.testUtils.setup1DSurvey(sigmaHalf)
+    problem = NSEM.Problem1D_eTotal(mesh)
+    problem.pair(survey)
+
+    # Get the fields
+    fields = problem.fields(sigma)
+
+    # Project the data
+    data = survey.eval(fields)
+
+    # Calculate the app  phs
+    app_p = np.array(NSEM.Utils.testUtils.getAppResPhs(data))[:,1]
+
+    return np.linalg.norm(np.abs(app_p - np.ones(survey.nFreq)*45)/ 45)
+
+def appRes_psFieldNorm(sigmaHalf):
+
+    # Make the survey
+    survey, sigma, mesh = NSEM.Utils.testUtils.setup1DSurvey(sigmaHalf,False)
+    problem = NSEM.Problem1D_ePrimSec(mesh, sigmaPrimary = sigma)
+    problem.pair(survey)
+
+    # Get the fields
+    fields = problem.fields(sigma)
+
+    # Project the data
+    data = survey.eval(fields)
+
+    # Calculate the app res and phs
+    app_r = np.array(NSEM.Utils.testUtils.getAppResPhs(data))[:,0]
+
+    return np.linalg.norm(np.abs(np.log(app_r) - np.log(np.ones(survey.nFreq)/sigmaHalf))*np.log(sigmaHalf))
+
+def appPhs_psFieldNorm(sigmaHalf):
+
+    # Make the survey
+    survey, sigma, mesh = NSEM.Utils.testUtils.setup1DSurvey(sigmaHalf,False)
+    problem = NSEM.Problem1D_ePrimSec(mesh, sigmaPrimary = sigma)
+    problem.pair(survey)
+
+    # Get the fields
+    fields = problem.fields(sigma)
+
+    # Project the data
+    data = survey.eval(fields)
+
+    # Calculate the app  phs
+    app_p = np.array(NSEM.Utils.testUtils.getAppResPhs(data))[:,1]
+
+    return np.linalg.norm(np.abs(app_p - np.ones(survey.nFreq)*45)/ 45)
+
+class TestAnalytics(unittest.TestCase):
+
+    def setUp(self):
+        pass
+    # Total Fields
+    # def test_appRes2en1(self):self.assertLess(appRes_TotalFieldNorm(2e-1), TOLr)
+    # def test_appPhs2en1(self):self.assertLess(appPhs_TotalFieldNorm(2e-1), TOLp)
+
+    # Primary/secondary
+    def test_appRes1en0_ps(self):self.assertLess(appRes_psFieldNorm(1e-0), TOLr)
+    def test_appPhs1en0_ps(self):self.assertLess(appPhs_psFieldNorm(1e-0), TOLp)
+    def test_appRes2en1_ps(self):self.assertLess(appRes_psFieldNorm(2e-1), TOLr)
+    def test_appPhs2en1_ps(self):self.assertLess(appPhs_psFieldNorm(2e-1), TOLp)
+    def test_appRes2en3_ps(self):self.assertLess(appRes_psFieldNorm(2e-3), TOLr)
+    def test_appPhs2en3_ps(self):self.assertLess(appPhs_psFieldNorm(2e-3), TOLp)
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/mt/forward/test_Problem3D_VsAnalyticSolution.py

@@ -0,0 +1,54 @@
+# Test functions
+from glob import glob
+import numpy as np, sys, os, time, scipy, subprocess
+import SimPEG as simpeg
+import unittest
+from SimPEG import NSEM
+from SimPEG.Utils import meshTensor
+from scipy.constants import mu_0
+
+np.random.seed(1100)
+
+TOLr = 1
+TOLp = 2
+FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
+CONDUCTIVITY = 1e1
+MU = mu_0
+freq = [1e-1, 2e-1]
+addrandoms = True
+
+def appResPhsHalfspace_eFrom_ps_Norm(sigmaHalf,appR=True,expMap=False):
+    if appR:
+        label = 'resistivity'
+    else:
+        label = 'phase'
+    print 'Apperent {:s} test of eFormulation primary/secondary at {:g}\n\n'.format(label,sigmaHalf)
+
+    # Calculate the app  phs
+    survey, problem = NSEM.Utils.testUtils.setupSimpegNSEM_ePrimSec(NSEM.Utils.testUtils.halfSpace(sigmaHalf),expMap=expMap)
+    data = problem.dataPair(survey,survey.dpred(problem.curModel))
+    recData = data.toRecArray('Complex')
+    app_rpxy, app_rpyx = NSEM.Utils.appResPhs(recData['freq'],recData['zxy'])[0], NSEM.Utils.appResPhs(recData['freq'],recData['zyx'])[0]
+    if appR:
+        return np.linalg.norm( np.abs(np.log10(app_rpxy[0]) - np.log10(1./sigmaHalf)) * np.log10(sigmaHalf ))
+    else:
+        return np.linalg.norm( np.abs(app_rpxy[1] + 135) / 135 )
+
+
+class TestAnalytics(unittest.TestCase):
+
+    def setUp(self):
+        # Make the survey and the problem
+        pass
+
+    # # Test apparent resistivity and phase
+    def test_appRes1en2(self):self.assertLess(appResPhsHalfspace_eFrom_ps_Norm(1e-2),TOLr)
+    def test_appPhs1en2(self):self.assertLess(appResPhsHalfspace_eFrom_ps_Norm(1e-2,False),TOLp)
+
+    def test_appRes1en1(self):self.assertLess(appResPhsHalfspace_eFrom_ps_Norm(1e-1),TOLr)
+    def test_appPhs1en1(self):self.assertLess(appResPhsHalfspace_eFrom_ps_Norm(1e-1,False),TOLp)
+
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/mt/inversion/__init__.py

@@ -0,0 +1,12 @@
+import os
+import glob
+import unittest
+
+if __name__ == '__main__':
+    test_file_strings = glob.glob('test_*.py')
+    module_strings = [str[0:len(str)-3] for str in test_file_strings]
+    suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str
+              in module_strings]
+    testSuite = unittest.TestSuite(suites)
+
+    unittest.TextTestRunner(verbosity=2).run(testSuite)

tests/mt/inversion/test_Problem3D_Adjoint.py

@@ -0,0 +1,58 @@
+# Test functions
+from glob import glob
+import numpy as np, sys, os, time, scipy, subprocess
+import SimPEG as simpeg
+import unittest
+from SimPEG import NSEM
+from SimPEG.Utils import meshTensor
+from scipy.constants import mu_0
+
+
+TOLr = 5e-2
+TOL = 1e-4
+FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
+CONDUCTIVITY = 1e1
+MU = mu_0
+freq = [1e-1, 2e-1]
+addrandoms = True
+
+
+
+def JvecAdjointTest(inputSetup,comp='All',freq=False):
+    (M, freqs, sig, sigBG, rx_loc) = inputSetup
+    survey, problem = NSEM.Utils.testUtils.setupSimpegNSEM_ePrimSec(inputSetup,comp='All',singleFreq=freq)
+    print 'Adjoint test of eForm primary/secondary for {:s} comp at {:s}\n'.format(comp,str(survey.freqs))
+
+    m  = sig
+    u = problem.fields(m)
+
+    v = np.random.rand(survey.nD,)
+    # print problem.PropMap.PropModel.nP
+    w = np.random.rand(problem.mesh.nC,)
+
+    vJw = v.ravel().dot(problem.Jvec(m, w, u))
+    wJtv = w.ravel().dot(problem.Jtvec(m, v, u))
+    tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR])
+    print ' vJw   wJtv  vJw - wJtv     tol    abs(vJw - wJtv) < tol'
+    print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol
+    return np.abs(vJw - wJtv) < tol
+
+
+class NSEM_AdjointTests(unittest.TestCase):
+
+    def setUp(self):
+        pass
+
+    # Test the adjoint of Jvec and Jtvec
+    # def test_JvecAdjoint_zxxr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxxr',.1))
+    # def test_JvecAdjoint_zxxi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxxi',.1))
+    # def test_JvecAdjoint_zxyr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxyr',.1))
+    # def test_JvecAdjoint_zxyi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxyi',.1))
+    # def test_JvecAdjoint_zyxr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyxr',.1))
+    # def test_JvecAdjoint_zyxi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyxi',.1))
+    # def test_JvecAdjoint_zyyr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyyr',.1))
+    # def test_JvecAdjoint_zyyi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyyi',.1))
+    def test_JvecAdjoint_All(self):self.assertTrue(JvecAdjointTest(NSEM.Utils.testUtils.random(1e-2),'All',.1))
+
+if __name__ == '__main__':
+    unittest.main()

tests/mt/inversion/test_Problem3D_Derivs.py

@@ -0,0 +1,83 @@
+# Test functions
+from glob import glob
+import numpy as np, sys, os, time, scipy, subprocess
+import SimPEG as simpeg
+import unittest
+from SimPEG import NSEM
+from SimPEG.Utils import meshTensor
+from scipy.constants import mu_0
+
+np.random.seed(1100)
+
+TOLr = 5e-2
+TOL = 1e-4
+FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
+CONDUCTIVITY = 1e1
+MU = mu_0
+freq = [1e-1, 2e-1]
+addrandoms = True
+
+
+# Test the Jvec derivative
+def DerivJvecTest(inputSetup,comp='All',freq=False,expMap=True):
+    (M, freqs, sig, sigBG, rx_loc) = inputSetup
+    survey, problem = NSEM.Utils.testUtils.setupSimpegNSEM_ePrimSec(inputSetup,comp=comp,singleFreq=freq,expMap=expMap)
+    print 'Derivative test of Jvec for eForm primary/secondary for {:s} comp at {:s}\n'.format(comp,survey.freqs)
+    # problem.mapping = simpeg.Maps.ExpMap(problem.mesh)
+    # problem.sigmaPrimary = np.log(sigBG)
+    x0 = np.log(sigBG)
+    # cond = sig[0]
+    # x0 = np.log(np.ones(problem.mesh.nC)*cond)
+    # problem.sigmaPrimary = x0
+    # if True:
+    #     x0  = x0 + np.random.randn(problem.mesh.nC)*cond*1e-1
+    survey = problem.survey
+    def fun(x):
+        return survey.dpred(x), lambda x: problem.Jvec(x0, x)
+    return simpeg.Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=FLR)
+
+def DerivProjfieldsTest(inputSetup,comp='All',freq=False):
+
+    survey, problem = NSEM.Utils.testUtils.setupSimpegNSEM_ePrimSec(inputSetup,comp,freq)
+    print 'Derivative test of data projection for eFormulation primary/secondary\n\n'
+    # problem.mapping = simpeg.Maps.ExpMap(problem.mesh)
+    # Initate things for the derivs Test
+    src = survey.srcList[0]
+    rx = src.rxList[0]
+
+    u0x = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j
+    u0y = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j
+    u0 = np.vstack((simpeg.mkvc(u0x,2),simpeg.mkvc(u0y,2)))
+    f0 = problem.fieldsPair(survey.mesh,survey)
+    # u0 = np.hstack((simpeg.mkvc(u0_px,2),simpeg.mkvc(u0_py,2)))
+    f0[src,'e_pxSolution'] =  u0[:len(u0)/2]#u0x
+    f0[src,'e_pySolution'] = u0[len(u0)/2::]#u0y
+
+    def fun(u):
+        f = problem.fieldsPair(survey.mesh,survey)
+        f[src,'e_pxSolution'] = u[:len(u)/2]
+        f[src,'e_pySolution'] = u[len(u)/2::]
+        return rx.eval(src,survey.mesh,f), lambda t: rx.evalDeriv(src,survey.mesh,f0,simpeg.mkvc(t,2))
+
+    return simpeg.Tests.checkDerivative(fun, u0, num=3, plotIt=False, eps=FLR)
+
+
+
+class NSEM_DerivTests(unittest.TestCase):
+
+    def setUp(self):
+        pass
+
+    # Do a derivative test of Jvec
+    # def test_derivJvec_zxxr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxxr',.1))
+    # def test_derivJvec_zxxi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxxi',.1))
+    # def test_derivJvec_zxyr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxyr',.1))
+    # def test_derivJvec_zxyi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxyi',.1))
+    # def test_derivJvec_zyxr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyxr',.1))
+    # def test_derivJvec_zyxi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyxi',.1))
+    # def test_derivJvec_zyyr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyyr',.1))
+    # def test_derivJvec_zyyi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyyi',.1))
+    def test_derivJvec_All(self):self.assertTrue(DerivJvecTest(NSEM.Utils.testUtils.random(1e-2),'All',.1))
+
+if __name__ == '__main__':
+    unittest.main()

tests/mt/test_Problem1D_againstAnalyticHalfspace.py

@@ -1,162 +0,0 @@
-import unittest
-import SimPEG as simpeg
-from SimPEG import MT
-from SimPEG.Utils import meshTensor
-import numpy as np
-# Define the tolerances
-TOLr = 5e-2
-TOLp = 5e-2
-
-
-def setupSurvey(sigmaHalf,tD=True):
-
-    # Frequency
-    nFreq = 33
-    freqs = np.logspace(3,-3,nFreq)
-    # Make the mesh
-    ct = 5
-    air = meshTensor([(ct,25,1.3)])
-    # coreT0 = meshTensor([(ct,15,1.2)])
-    # coreT1 = np.kron(meshTensor([(coreT0[-1],15,1.3)]),np.ones((7,)))
-    core = np.concatenate( (  np.kron(meshTensor([(ct,15,-1.2)]),np.ones((10,))) , meshTensor([(ct,20)]) ) )
-    bot = meshTensor([(core[0],10,-1.3)])
-    x0 = -np.array([np.sum(np.concatenate((core,bot)))])
-    m1d = simpeg.Mesh.TensorMesh([np.concatenate((bot,core,air))], x0=x0)
-    # Make the model
-    sigma = np.zeros(m1d.nC) + sigmaHalf
-    sigma[m1d.gridCC > 0 ] = 1e-8
-
-    rxList = []
-    for rxType in ['z1dr','z1di']:
-        rxList.append(MT.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType))
-    # Source list
-    srcList =[]
-    if tD:
-        for freq in freqs:
-            srcList.append(MT.SrcMT.polxy_1DhomotD(rxList,freq))
-    else:
-        for freq in freqs:
-            srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
-
-    survey = MT.Survey(srcList)
-    return survey, sigma, m1d
-
-def getAppResPhs(MTdata):
-    # Make impedance
-    def appResPhs(freq,z):
-        app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
-        app_phs = np.arctan2(z.imag,z.real)*(180/np.pi)
-        return app_res, app_phs
-    zList = []
-    for src in MTdata.survey.srcList:
-        zc = [src.freq]
-        for rx in src.rxList:
-            if 'i' in rx.rxType:
-                m=1j
-            else:
-                m = 1
-            zc.append(m*MTdata[src,rx])
-        zList.append(zc)
-    return [appResPhs(zList[i][0],np.sum(zList[i][1:3])) for i in np.arange(len(zList))]
-
-def appRes_TotalFieldNorm(sigmaHalf):
-
-    # Make the survey
-    survey, sigma, mesh = setupSurvey(sigmaHalf)
-    problem = MT.Problem1D.eForm_TotalField(mesh)
-    problem.pair(survey)
-
-    # Get the fields
-    fields = problem.fields(sigma)
-
-    # Project the data
-    data = survey.eval(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_TotalFieldNorm(sigmaHalf):
-
-    # Make the survey
-    survey, sigma, mesh = setupSurvey(sigmaHalf)
-    problem = MT.Problem1D.eForm_TotalField(mesh)
-    problem.pair(survey)
-
-    # Get the fields
-    fields = problem.fields(sigma)
-
-    # Project the data
-    data = survey.eval(fields)
-
-    # Calculate the app  phs
-    app_p = np.array(getAppResPhs(data))[:,1]
-
-    return np.linalg.norm(np.abs(app_p - np.ones(survey.nFreq)*45)/ 45)
-
-def appRes_psFieldNorm(sigmaHalf):
-
-    # Make the survey
-    survey, sigma, mesh = setupSurvey(sigmaHalf,False)
-    problem = MT.Problem1D.eForm_psField(mesh, sigmaPrimary = sigma)
-    problem.pair(survey)
-
-    # Get the fields
-    fields = problem.fields(sigma)
-
-    # Project the data
-    data = survey.eval(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 = MT.Problem1D.eForm_psField(mesh, sigmaPrimary = sigma)
-    problem.pair(survey)
-
-    # Get the fields
-    fields = problem.fields(sigma)
-
-    # Project the data
-    data = survey.eval(fields)
-
-    # Calculate the app  phs
-    app_p = np.array(getAppResPhs(data))[:,1]
-
-    return np.linalg.norm(np.abs(app_p - np.ones(survey.nFreq)*45)/ 45)
-
-class TestAnalytics(unittest.TestCase):
-
-    def setUp(self):
-        pass
-    # 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()

tests/mt/test_Problem3D_againstAnalytic.py

@@ -1,268 +0,0 @@
-# Test functions
-from glob import glob
-import numpy as np, sys, os, time, scipy, subprocess
-import SimPEG as simpeg
-import unittest
-from SimPEG import MT
-from SimPEG.Utils import meshTensor
-from scipy.constants import mu_0
-
-TOLr = 5e-2
-TOL = 1e-4
-FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
-CONDUCTIVITY = 1e1
-MU = mu_0
-freq = [1e-1, 2e-1]
-addrandoms = True
-
-
-def getInputs():
-    """
-    Function that returns Mesh, freqs, rx_loc, elev.
-    """
-    # Make a mesh
-    # M = simpeg.Mesh.TensorMesh([[(100,5,-1.5),(100.,10),(100,5,1.5)],[(100,5,-1.5),(100.,10),(100,5,1.5)],[(100,5,1.6),(100.,10),(100,3,2)]], x0=['C','C',-3529.5360])
-    # M = simpeg.Mesh.TensorMesh([[(1000,6,-1.5),(1000.,6),(1000,6,1.5)],[(1000,6,-1.5),(1000.,2),(1000,6,1.5)],[(1000,6,-1.3),(1000.,6),(1000,6,1.3)]], x0=['C','C','C'])# Setup the model
-    M = simpeg.Mesh.TensorMesh([[(1000,6,-1.5),(1000.,4),(1000,6,1.5)],[(1000,6,-1.5),(1000.,4),(1000,6,1.5)],[(500,8,-1.3),(500.,8),(500,8,1.3)]], x0=['C','C','C'])# Setup the model
-    # Set the frequencies
-    freqs = np.logspace(1,-3,5)
-    elev = 0
-
-    ## Setup the the survey object
-    # Receiver locations
-    rx_x, rx_y = np.meshgrid(np.arange(-1000,1001,500),np.arange(-1000,1001,500))
-    rx_loc = np.hstack((simpeg.Utils.mkvc(rx_x,2),simpeg.Utils.mkvc(rx_y,2),elev+np.zeros((np.prod(rx_x.shape),1))))
-
-    return M, freqs, rx_loc, elev
-
-def random(conds):
-    ''' Returns a halfspace model based on the inputs'''
-    M, freqs, rx_loc, elev = getInputs()
-
-    # Backround
-    sigBG = np.ones(M.nC)*conds
-    # Add randomness to the model (10% of the value).
-    sig = np.exp( np.log(sigBG) + np.random.randn(M.nC)*(conds)*1e-1 )
-
-    return (M, freqs, sig, sigBG, rx_loc)
-
-def halfSpace(conds):
-    ''' Returns a halfspace model based on the inputs'''
-    M, freqs, rx_loc, elev = getInputs()
-
-    # Model
-    ccM = M.gridCC
-    # conds = [1e-2]
-    groundInd = ccM[:,2] < elev
-    sig = np.zeros(M.nC) + 1e-8
-    sig[groundInd] = conds
-    # Set the background, not the same as the model
-    sigBG = np.zeros(M.nC) + 1e-8
-    sigBG[groundInd] = conds
-
-    return (M, freqs, sig, sigBG, rx_loc)
-
-def blockInhalfSpace(conds):
-    ''' Returns a halfspace model based on the inputs'''
-    M, freqs, rx_loc, elev = getInputs()
-
-    # Model
-    ccM = M.gridCC
-    # conds = [1e-2]
-    groundInd = ccM[:,2] < elev
-    sig = simpeg.Utils.ModelBuilder.defineBlock(M.gridCC,np.array([-1000,-1000,-1500]),np.array([1000,1000,-1000]),conds)
-    sig[~groundInd] = 1e-8
-    # Set the background, not the same as the model
-    sigBG = np.zeros(M.nC) + 1e-8
-    sigBG[groundInd] = conds[1]
-
-    return (M, freqs, sig, sigBG, rx_loc)
-
-def twoLayer(conds):
-    ''' Returns a 2 layer model based on the conductivity values given'''
-    M, freqs, rx_loc, elev = getInputs()
-
-    # Model
-    ccM = M.gridCC
-    groundInd = ccM[:,2] < elev
-    botInd = ccM[:,2] < -3000
-    sig = np.zeros(M.nC) + 1e-8
-    sig[groundInd] = conds[1]
-    sig[botInd] = conds[0]
-    # Set the background, not the same as the model
-    sigBG = np.zeros(M.nC) + 1e-8
-    sigBG[groundInd] = conds[1]
-
-
-    return (M, freqs, sig, sigBG, rx_loc)
-
-
-
-def setupSimpegMTfwd_eForm_ps(inputSetup,comp='Imp',singleFreq=False,expMap=True):
-    M,freqs,sig,sigBG,rx_loc = inputSetup
-    # Make a receiver list
-    rxList = []
-    if comp == 'All':
-        for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi','tzxr','tzxi','tzyr','tzyi']:
-            rxList.append(MT.Rx(rx_loc,rxType))
-    elif comp == 'Imp':
-        for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi']:
-            rxList.append(MT.Rx(rx_loc,rxType))
-    elif comp == 'Tip':
-        for rxType in ['tzxr','tzxi','tzyr','tzyi']:
-            rxList.append(MT.Rx(rx_loc,rxType))
-    else:
-        rxList.append(MT.Rx(rx_loc,comp))
-    # Source list
-    srcList =[]
-
-    if singleFreq:
-        srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,singleFreq))
-    else:
-        for freq in freqs:
-            srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
-    # Survey MT
-    survey = MT.Survey(srcList)
-
-    ## Setup the problem object
-    sigma1d = M.r(sigBG,'CC','CC','M')[0,0,:]
-    if expMap:
-        problem = MT.Problem3D.eForm_ps(M,sigmaPrimary= np.log(sigma1d) )
-        problem.mapping = simpeg.Maps.ExpMap(problem.mesh)
-        problem.curModel = np.log(sig)
-    else:
-        problem = MT.Problem3D.eForm_ps(M,sigmaPrimary= sigma1d)
-        problem.curModel = sig
-    problem.pair(survey)
-    problem.verbose = False
-    try:
-        from pymatsolver import MumpsSolver
-        problem.Solver = MumpsSolver
-    except:
-        pass
-
-    return (survey, problem)
-
-def getAppResPhs(MTdata):
-    # Make impedance
-    def appResPhs(freq,z):
-        app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
-        app_phs = np.arctan2(z.imag,z.real)*(180/np.pi)
-        return app_res, app_phs
-    recData = MTdata.toRecArray('Complex')
-    return appResPhs(recData['freq'],recData['zxy']), appResPhs(recData['freq'],recData['zyx'])
-
-def JvecAdjointTest(inputSetup,comp='All',freq=False):
-    (M, freqs, sig, sigBG, rx_loc) = inputSetup
-    survey, problem = setupSimpegMTfwd_eForm_ps(inputSetup,comp='All',singleFreq=freq)
-    print 'Adjoint test of eForm primary/secondary for {:s} comp at {:s}\n'.format(comp,str(survey.freqs))
-
-    m  = sig
-    u = problem.fields(m)
-
-    v = np.random.rand(survey.nD,)
-    # print problem.PropMap.PropModel.nP
-    w = np.random.rand(problem.mesh.nC,)
-
-    vJw = v.ravel().dot(problem.Jvec(m, w, u))
-    wJtv = w.ravel().dot(problem.Jtvec(m, v, u))
-    tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR])
-    print ' vJw   wJtv  vJw - wJtv     tol    abs(vJw - wJtv) < tol'
-    print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol
-    return np.abs(vJw - wJtv) < tol
-
-# Test the Jvec derivative
-def DerivJvecTest(inputSetup,comp='All',freq=False,expMap=True):
-    (M, freqs, sig, sigBG, rx_loc) = inputSetup
-    survey, problem = setupSimpegMTfwd_eForm_ps(inputSetup,comp=comp,singleFreq=freq,expMap=expMap)
-    print 'Derivative test of Jvec for eForm primary/secondary for {:s} comp at {:s}\n'.format(comp,survey.freqs)
-    # problem.mapping = simpeg.Maps.ExpMap(problem.mesh)
-    # problem.sigmaPrimary = np.log(sigBG)
-    x0 = np.log(sigBG)
-    # cond = sig[0]
-    # x0 = np.log(np.ones(problem.mesh.nC)*cond)
-    # problem.sigmaPrimary = x0
-    # if True:
-    #     x0  = x0 + np.random.randn(problem.mesh.nC)*cond*1e-1
-    survey = problem.survey
-    def fun(x):
-        return survey.dpred(x), lambda x: problem.Jvec(x0, x)
-    return simpeg.Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=FLR)
-
-def DerivProjfieldsTest(inputSetup,comp='All',freq=False):
-
-    survey, problem = setupSimpegMTfwd_eForm_ps(inputSetup,comp,freq)
-    print 'Derivative test of data projection for eFormulation primary/secondary\n\n'
-    # problem.mapping = simpeg.Maps.ExpMap(problem.mesh)
-    # Initate things for the derivs Test
-    src = survey.srcList[0]
-    rx = src.rxList[0]
-
-    u0x = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j
-    u0y = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j
-    u0 = np.vstack((simpeg.mkvc(u0x,2),simpeg.mkvc(u0y,2)))
-    f0 = problem.fieldsPair(survey.mesh,survey)
-    # u0 = np.hstack((simpeg.mkvc(u0_px,2),simpeg.mkvc(u0_py,2)))
-    f0[src,'e_pxSolution'] =  u0[:len(u0)/2]#u0x
-    f0[src,'e_pySolution'] = u0[len(u0)/2::]#u0y
-
-    def fun(u):
-        f = problem.fieldsPair(survey.mesh,survey)
-        f[src,'e_pxSolution'] = u[:len(u)/2]
-        f[src,'e_pySolution'] = u[len(u)/2::]
-        return rx.eval(src,survey.mesh,f), lambda t: rx.evalDeriv(src,survey.mesh,f0,simpeg.mkvc(t,2))
-
-    return simpeg.Tests.checkDerivative(fun, u0, num=3, plotIt=False, eps=FLR)
-
-def appResPhsHalfspace_eFrom_ps_Norm(sigmaHalf,appR=True,expMap=False):
-    if appR:
-        label = 'resistivity'
-    else:
-        label = 'phase'
-    # Make the survey and the problem
-    survey, problem = setupSimpegMTfwd_eForm_ps(halfSpace(sigmaHalf),expMap=expMap)
-    print 'Apperent {:s} test of eFormulation primary/secondary at {:g}\n\n'.format(label,sigmaHalf)
-
-    data = problem.dataPair(survey,survey.dpred(problem.curModel))
-    # Calculate the app  phs
-    app_rpxy, app_rpyx = np.array(getAppResPhs(data))
-    if appR:
-        return np.all(np.abs(app_rpxy[0,:] - 1./sigmaHalf) * sigmaHalf < .4)
-    else:
-        return np.all(np.abs(app_rpxy[1,:] + 135) / 135 < .4)
-
-class TestAnalytics(unittest.TestCase):
-
-    def setUp(self):
-        pass
-    # # Test apparent resistivity and phase
-    def test_appRes1en2(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-2))
-    def test_appPhs1en2(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-2,False))
-
-    def test_appRes1en3(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-3))
-    def test_appPhs1en3(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-3,False))
-
-    # Do a derivative test of Jvec
-    # def test_derivJvec_zxxr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxxr',.1))
-    # def test_derivJvec_zxxi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxxi',.1))
-    # def test_derivJvec_zxyr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxyr',.1))
-    # def test_derivJvec_zxyi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxyi',.1))
-    # def test_derivJvec_zyxr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyxr',.1))
-    # def test_derivJvec_zyxi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyxi',.1))
-    # def test_derivJvec_zyyr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyyr',.1))
-    # def test_derivJvec_zyyi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyyi',.1))
-    def test_derivJvec_All(self):self.assertTrue(DerivJvecTest(random(1e-2),'All',.1))
-
-    # Test the adjoint of Jvec and Jtvec
-    # def test_JvecAdjoint_zxxr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxxr',.1))
-    # def test_JvecAdjoint_zxxi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxxi',.1))
-    # def test_JvecAdjoint_zxyr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxyr',.1))
-    # def test_JvecAdjoint_zxyi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxyi',.1))
-    # def test_JvecAdjoint_zyxr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyxr',.1))
-    # def test_JvecAdjoint_zyxi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyxi',.1))
-    # def test_JvecAdjoint_zyyr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyyr',.1))
-    # def test_JvecAdjoint_zyyi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyyi',.1))
-    def test_JvecAdjoint_All(self):self.assertTrue(JvecAdjointTest(random(1e-2),'All',.1))
-
-if __name__ == '__main__':
-    unittest.main()