no message

Better list comprehension.

SimPEG/DCIP/BaseDC.py

@@ -200,11 +200,11 @@ class ProblemDC_CC(Problem.BaseProblem):
 
         return F
 
-    def Jvec(self, m, v, u=None):
+    def Jvec(self, m, v, f=None):
         """
             :param numpy.array m: model
             :param numpy.array v: vector to multiply
-            :param numpy.array u: fields
+            :param Fields f: fields
             :rtype: numpy.array
             :return: Jv
 
@@ -225,11 +225,10 @@ class ProblemDC_CC(Problem.BaseProblem):
         # Set current model; clear dependent property $\mathbf{A(m)}$
         self.curModel = m
         sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$
-        if u is None:
+        if f is None:
             # Run forward simulation if $u$ not provided
-            u = self.fields(self.curModel)[self.survey.srcList, 'phi_sol']
-        else:
-            u = u[self.survey.srcList, 'phi_sol']
+            f = self.fields(self.curModel)
+        u = f[self.survey.srcList, 'phi_sol']
 
         D = self.mesh.faceDiv
         G = self.mesh.cellGrad
@@ -251,19 +250,18 @@ class ProblemDC_CC(Problem.BaseProblem):
         if self.Ainv is None:
             self.Ainv = self.Solver(dA_du, **self.solverOpts)
 
-        P        = self.survey.getP(self.mesh)
+        P     = self.survey.getP(self.mesh)
         Jv    = - P * mkvc( self.Ainv * dCdm_x_v )
         return Jv
 
-    def Jtvec(self, m, v, u=None):
+    def Jtvec(self, m, v, f=None):
 
         self.curModel = m
         sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$
-        if u is None:
-            # Run forward simulation if $u$ not provided
-            u = self.fields(self.curModel)[self.survey.srcList, 'phi_sol']
-        else:
-            u = u[self.survey.srcList, 'phi_sol']
+        if f is None:
+            # Run forward simulation if $f$ not provided
+            f = self.fields(self.curModel)
+        u = f[self.survey.srcList, 'phi_sol']
 
         shp = u.shape
         P = self.survey.getP(self.mesh)

SimPEG/DCIP/BaseIP.py

@@ -14,12 +14,12 @@ class SurveyIP(SurveyDC):
         Survey.BaseSurvey.__init__(self, **kwargs)
         self._Ps = {}
 
-    def dpred(self, m, u=None):
+    def dpred(self, m, f=None):
         """
             Predicted data.
 
             .. math::
-                d_\\text{pred} = Pu(m)
+                d_\\text{pred} = Pf(m)
         """
 
         return self.prob.forward(m)
@@ -143,10 +143,10 @@ class ProblemIP(Problem.BaseProblem):
         J_x_v    = - P * mkvc( self.Ainv * dCdm_x_v )
         return -J_x_v
 
-    def Jvec(self, m, v, u=None):
+    def Jvec(self, m, v, f=None):
         return self.forward(v)
 
-    def Jtvec(self, m, v, u=None):
+    def Jtvec(self, m, v, f=None):
 
         self.curModel = m
         # sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$

SimPEG/DCIP/DCIPUtils.py

@@ -169,7 +169,7 @@ def readUBC_DC2DModel(fileName):
 
     return model
 
-def plot_pseudoSection(DCsurvey, axs, stype):
+def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None):
     """
         Read list of 2D tx-rx location and plot a speudo-section of apparent
         resistivity.
@@ -179,7 +179,7 @@ def plot_pseudoSection(DCsurvey, axs, stype):
         Input:
         :param d2D, z0
         :switch stype -> Either 'pdp' (pole-dipole) | 'dpdp' (dipole-dipole)
-
+        :switch dtype=-> Either 'appr' (app. res) | 'appc' (app. con) | 'volt' (potential)
         Output:
         :figure scatter plot overlayed on image
 
@@ -221,20 +221,43 @@ def plot_pseudoSection(DCsurvey, axs, stype):
         Cmid = (Tx[0][0] + Tx[1][0])/2
         Pmid = (Rx[0][:,0] + Rx[1][:,0])/2
 
-        # Compute pant leg of apparent rho
-        if stype == 'pdp':
-            leg =  data * 2*np.pi  * MA * ( MA + MN ) / MN
+        # Change output for dtype
+        if dtype == 'volt':
 
-            leg = np.log10(abs(1/leg))
+            rho = np.hstack([rho,data])
 
-        elif stype == 'dpdp':
-            leg = data * 2*np.pi / ( 1/MA - 1/MB - 1/NB + 1/NA )
+        else:
 
+            # Compute pant leg of apparent rho
+            if stype == 'pdp':
+
+                leg =  data * 2*np.pi  * MA * ( MA + MN ) / MN
+
+            elif stype == 'dpdp':
+
+                leg = data * 2*np.pi / ( 1/MA - 1/MB - 1/NB + 1/NA )
+
+            else:
+                print """dtype must be 'pdp'(pole-dipole) | 'dpdp' (dipole-dipole) """
+                break
+
+
+            if dtype == 'appc':
+
+                leg = np.log10(abs(1./leg))
+                rho = np.hstack([rho,leg])
+
+            elif dtype == 'appr':
+
+                leg = np.log10(abs(leg))
+                rho = np.hstack([rho,leg])
+
+            else:
+                print """dtype must be 'appr' | 'appc' | 'volt' """
+                break
 
         midx = np.hstack([midx, ( Cmid + Pmid )/2 ])
         midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + z0 ])
-        rho = np.hstack([rho,leg])
-
 
     ax = axs
 
@@ -242,26 +265,38 @@ def plot_pseudoSection(DCsurvey, axs, stype):
     grid_x, grid_z = np.mgrid[np.min(midx):np.max(midx), np.min(midz):np.max(midz)]
     grid_rho = griddata(np.c_[midx,midz], rho.T, (grid_x, grid_z), method='linear')
 
+    if clim == None:
+        vmin, vmax = rho.min(), rho.max()
+    else:
+        vmin, vmax = clim[0], clim[1]
 
-    plt.imshow(grid_rho.T, extent = (np.min(midx),np.max(midx),np.min(midz),np.max(midz)), origin='lower', alpha=0.8, vmin = np.min(rho), vmax = np.max(rho))
-    cbar = plt.colorbar(format = '%.2f',fraction=0.04,orientation="horizontal")
+    grid_rho = np.ma.masked_where(np.isnan(grid_rho), grid_rho)
+    ph = plt.pcolormesh(grid_x[:,0],grid_z[0,:],grid_rho.T, clim=(vmin, vmax))
+    cbar = plt.colorbar(format="$10^{%.1f}$",fraction=0.04,orientation="horizontal")
 
     cmin,cmax = cbar.get_clim()
     ticks = np.linspace(cmin,cmax,3)
     cbar.set_ticks(ticks)
+    cbar.ax.tick_params(labelsize=10)
+    
+    if dtype == 'appc':
+        cbar.set_label("App.Cond",size=12)
+    elif dtype == 'appr':
+        cbar.set_label("App.Res.",size=12)
+    elif dtype == 'volt':
+        cbar.set_label("Potential (V)",size=12)
 
     # Plot apparent resistivity
-    plt.scatter(midx,midz,s=50,c=rho.T)
+    ax.scatter(midx,midz,s=10,c=rho.T, vmin =vmin, vmax = vmax, clim=(vmin, vmax))
 
-    ax.set_xticklabels([])
+    #ax.set_xticklabels([])
+    #ax.set_yticklabels([])
 
-    ax.set_ylabel('Z')
-    ax.yaxis.tick_right()
-    ax.yaxis.set_label_position('right')
     plt.gca().set_aspect('equal', adjustable='box')
 
 
-    return ax
+
+    return ph
 
 def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
     """
@@ -361,16 +396,6 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
                 srcClass = DC.SrcDipole([rxClass], M[ii,:],M[ii,:])
             SrcList.append(srcClass)
 
-#==============================================================================
-#     elif re.match(stype,'dpdp'):
-#
-#         for ii in range(0, int(nstn)-2):
-#
-#             indx = np.min([ii+n+1,nstn])
-#             Tx.append(np.c_[M[ii,:],N[ii,:]])
-#             Rx.append(np.c_[M[ii+2:indx,:],N[ii+2:indx,:]])
-#==============================================================================
-
     elif stype == 'gradient':
 
         # Gradient survey only requires Tx at end of line and creates a square
@@ -513,22 +538,22 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype, stype):
 
     fid.close()
 
-def convertObs_DC3D_to_2D(DCsurvey,lineID):
+def convertObs_DC3D_to_2D(DCsurvey,lineID, flag = 'local'):
     """
-        Read DC survey and data and change
-        coordinate system to distance along line assuming
-        all data is acquired along line.
-        First transmitter pole is assumed to be at the origin
+        Read DC survey and projects the coordinate system
+        according to the flag = 'Xloc' | 'Yloc' | 'local' (default)
+        In the 'local' system, station coordinates are referenced
+        to distance from the first srcLoc[0].loc[0]
 
-        Assumes flat topo for now...
+        The Z value is preserved, but Y coordinates zeroed.
 
         Input:
-        :param Tx, Rx
+        :param survey3D
 
         Output:
-        :figure Tx2d, Rx2d
+        :figure survey2D
 
-        Edited Feb 17th, 2016
+        Edited April 6th, 2016
 
         @author: dominiquef
 
@@ -570,25 +595,39 @@ def convertObs_DC3D_to_2D(DCsurvey,lineID):
             Rx = DCsurvey.srcList[indx[ii]].rxList[0].locs
             nrx = Rx[0].shape[0]
 
-            # Find A electrode along line
-            vec, r = r_unit(x0,Tx[ii][0,0:2])
-            A = stn_id(vecTx,vec,r)
+            if flag == 'local':
+                # Find A electrode along line
+                vec, r = r_unit(x0,Tx[ii][0,0:2])
+                A = stn_id(vecTx,vec,r)
 
-            # Find B electrode along line
-            vec, r = r_unit(x0,Tx[ii][1,0:2])
-            B = stn_id(vecTx,vec,r)
+                # Find B electrode along line
+                vec, r = r_unit(x0,Tx[ii][1,0:2])
+                B = stn_id(vecTx,vec,r)
 
-            M = np.zeros(nrx)
-            N = np.zeros(nrx)
-            for kk in range(nrx):
+                M = np.zeros(nrx)
+                N = np.zeros(nrx)
+                for kk in range(nrx):
 
-                # Find all M electrodes along line
-                vec, r = r_unit(x0,Rx[0][kk,0:2])
-                M[kk] = stn_id(vecTx,vec,r)
+                    # Find all M electrodes along line
+                    vec, r = r_unit(x0,Rx[0][kk,0:2])
+                    M[kk] = stn_id(vecTx,vec,r)
 
-                # Find all N electrodes along line
-                vec, r = r_unit(x0,Rx[1][kk,0:2])
-                N[kk] = stn_id(vecTx,vec,r)
+                    # Find all N electrodes along line
+                    vec, r = r_unit(x0,Rx[1][kk,0:2])
+                    N[kk] = stn_id(vecTx,vec,r)
+            elif flag == 'Yloc':
+                """ Flip the XY axis locs"""
+                A = Tx[ii][0,1]
+                B = Tx[ii][1,1]
+                M = Rx[0][:,1]
+                N = Rx[1][:,1]
+
+            elif flag == 'Xloc':
+                """ Copy the rx-tx locs"""
+                A = Tx[ii][0,0]
+                B = Tx[ii][1,0]
+                M = Rx[0][:,0]
+                N = Rx[1][:,0]
 
             Rx = DC.RxDipole(np.c_[M,np.zeros(nrx),Rx[0][:,2]],np.c_[N,np.zeros(nrx),Rx[1][:,2]])
 
@@ -604,16 +643,16 @@ def convertObs_DC3D_to_2D(DCsurvey,lineID):
 
 def readUBC_DC3Dobs(fileName):
     """
-        Read UBC GIF DCIP 3D observation file and generate arrays for tx-rx location
+        Read UBC GIF DCIP 3D observation file and generate survey
 
         Input:
         :param fileName, path to the UBC GIF 3D obs file
 
         Output:
-        :param rx, tx, d, wd
+        :param DCIPsurvey
         :return
 
-        Created on Mon December 7th, 2015
+        Created on Mon April 6th, 2015
 
         @author: dominiquef
 
@@ -688,6 +727,7 @@ def readUBC_DC3Dobs(fileName):
 
 def readUBC_DC2Dobs(fileName):
     """
+        ------- NEEDS TO BE UPDATED ------
         Read UBC GIF 2D observation file and generate arrays for tx-rx location
 
         Input:
@@ -735,6 +775,73 @@ def readUBC_DC2Dobs(fileName):
 
     return tx, rx, d, wd
 
+def readUBC_DC2Dpre(fileName):
+    """
+        Read UBC GIF DCIP 2D observation file and generate arrays for tx-rx location
+
+        Input:
+        :param fileName, path to the UBC GIF 3D obs file
+
+        Output:
+        DCsurvey
+        :return
+
+        Created on Mon March 9th, 2016 << Doug's 70th Birthday !! >>
+
+        @author: dominiquef
+
+    """
+
+    # Load file
+    obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!')
+
+    # Pre-allocate
+    srcLists = []
+    Rx = []
+    d = []
+    zflag = True # Flag for z value provided
+
+    for ii in range(obsfile.shape[0]):
+
+        if not obsfile[ii]:
+            continue
+
+        # First line is transmitter with number of receivers
+
+
+        temp = (np.fromstring(obsfile[ii], dtype=float,sep=' ').T)
+
+
+        # Check if z value is provided, if False -> nan
+        if len(temp)==5:
+            tx = np.r_[temp[0],np.nan,np.nan,temp[1],np.nan,np.nan]
+            zflag = False
+
+        else:
+            tx = np.r_[temp[0],np.nan,temp[1],temp[2],np.nan,temp[3]]
+
+
+        if zflag:
+            rx = np.c_[temp[4],np.nan,temp[5],temp[6],np.nan,temp[7]]
+
+
+        else:
+            rx = np.c_[temp[2],np.nan,np.nan,temp[3],np.nan,np.nan]
+            # Check if there is data with the location
+
+        d.append(temp[-1])
+
+
+        Rx = DC.RxDipole(rx[:,:3],rx[:,3:])
+        srcLists.append( DC.SrcDipole( [Rx], tx[:3],tx[3:]) )
+
+    # Create survey class
+    survey = DC.SurveyDC(srcLists)
+
+    survey.dobs = np.asarray(d)
+
+    return {'DCsurvey':survey}
+
 def readUBC_DC2DMesh(fileName):
     """
         Read UBC GIF 2DTensor mesh and generate 2D Tensor mesh in simpeg
@@ -928,7 +1035,6 @@ def getSrc_locs(DCsurvey):
 
     srcMat = np.zeros((DCsurvey.nSrc,2,3))
     for ii in range(DCsurvey.nSrc):
-        print np.asarray(DCsurvey.srcList[ii].loc).shape
         srcMat[ii,:,:] =  np.asarray(DCsurvey.srcList[ii].loc)
 
     return srcMat

SimPEG/DataMisfit.py

@@ -22,11 +22,11 @@ class BaseDataMisfit(object):
         Utils.setKwargs(self,**kwargs)
 
     @Utils.timeIt
-    def eval(self, m, u=None):
-        """eval(m, u=None)
+    def eval(self, m, f=None):
+        """eval(m, f=None)
 
             :param numpy.array m: geophysical model
-            :param numpy.array u: fields
+            :param Fields f: fields
             :rtype: float
             :return: data misfit
 
@@ -34,11 +34,11 @@ class BaseDataMisfit(object):
         raise NotImplementedError('This method should be overwritten.')
 
     @Utils.timeIt
-    def evalDeriv(self, m, u=None):
-        """evalDeriv(m, u=None)
+    def evalDeriv(self, m, f=None):
+        """evalDeriv(m, f=None)
 
             :param numpy.array m: geophysical model
-            :param numpy.array u: fields
+            :param Fields f: fields
             :rtype: numpy.array
             :return: data misfit derivative
 
@@ -47,12 +47,12 @@ class BaseDataMisfit(object):
 
 
     @Utils.timeIt
-    def eval2Deriv(self, m, v, u=None):
-        """eval2Deriv(m, v, u=None)
+    def eval2Deriv(self, m, v, f=None):
+        """eval2Deriv(m, v, f=None)
 
             :param numpy.array m: geophysical model
             :param numpy.array v: vector to multiply
-            :param numpy.array u: fields
+            :param Fields f: fields
             :rtype: numpy.array
             :return: data misfit derivative
 
@@ -89,7 +89,7 @@ class l2_DataMisfit(BaseDataMisfit):
         """
 
         if getattr(self, '_Wd', None) is None:
-            
+
             survey = self.survey
 
             if getattr(survey,'std', None) is None:
@@ -108,24 +108,20 @@ class l2_DataMisfit(BaseDataMisfit):
         self._Wd = value
 
     @Utils.timeIt
-    def eval(self, m, u=None):
-        "eval(m, u=None)"
-        prob   = self.prob
-        survey = self.survey
-        R = self.Wd * survey.residual(m, u=u)
+    def eval(self, m, f=None):
+        "eval(m, f=None)"
+        if f is None: f = self.prob.fields(m)
+        R = self.Wd * self.survey.residual(m, f)
         return 0.5*np.vdot(R, R)
 
     @Utils.timeIt
-    def evalDeriv(self, m, u=None):
-        "evalDeriv(m, u=None)"
-        prob   = self.prob
-        survey = self.survey
-        if u is None: u = prob.fields(m)
-        return prob.Jtvec(m, self.Wd * (self.Wd * survey.residual(m, u=u)), u=u)
+    def evalDeriv(self, m, f=None):
+        "evalDeriv(m, f=None)"
+        if f is None: f = self.prob.fields(m)
+        return self.prob.Jtvec(m, self.Wd * (self.Wd * self.survey.residual(m, f=f)), f=f)
 
     @Utils.timeIt
-    def eval2Deriv(self, m, v, u=None):
-        "eval2Deriv(m, v, u=None)"
-        prob   = self.prob
-        if u is None: u = prob.fields(m)
-        return prob.Jtvec_approx(m, self.Wd * (self.Wd * prob.Jvec_approx(m, v, u=u)), u=u)
+    def eval2Deriv(self, m, v, f=None):
+        "eval2Deriv(m, v, f=None)"
+        if f is None: f = self.prob.fields(m)
+        return self.prob.Jtvec_approx(m, self.Wd * (self.Wd * self.prob.Jvec_approx(m, v, f=f)), f=f)

SimPEG/Directives.py

@@ -123,10 +123,10 @@ class BetaEstimate_ByEig(InversionDirective):
         if self.debug: print 'Calculating the beta0 parameter.'
 
         m = self.invProb.curModel
-        u = self.invProb.getFields(m, store=True, deleteWarmstart=False)
+        f = self.invProb.getFields(m, store=True, deleteWarmstart=False)
 
         x0 = np.random.rand(*m.shape)
-        t = x0.dot(self.dmisfit.eval2Deriv(m,x0,u=u))
+        t = x0.dot(self.dmisfit.eval2Deriv(m,x0,f=f))
         b = x0.dot(self.reg.eval2Deriv(m, v=x0))
         self.beta0 = self.beta0_ratio*(t/b)
 
@@ -216,7 +216,7 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
         # 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.smoothModel == True:
+        if self.reg.mrefInSmooth == True:
             mref = self.reg.mref
         else:
             mref = 0
@@ -249,7 +249,7 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
         # 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.smoothModel == True:
+        if self.reg.mrefInSmooth == True:
             mref = self.reg.mref
         else:
             mref = 0
@@ -271,7 +271,6 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
         np.savez('{:s}-{:03d}'.format(self.fileName,self.opt.iter), 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)
 
 
-
 # class UpdateReferenceModel(Parameter):
 
 #     mref0 = None
@@ -283,3 +282,63 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
 #             mref = self.mref0
 #         self.m_prev = self.invProb.m_current
 #         return mref
+
+class update_IRLS(InversionDirective):
+
+    eps_min = None
+    factor = None
+    gamma = None
+    phi_m_last = None
+    phi_d_last = None
+    
+    def initialize(self):
+
+        # Scale the regularization for changes in norm
+        if getattr(self, 'phi_m_last', None) is not None:
+            self.reg.gamma = 1.
+            phim_new = self.reg.eval(self.invProb.curModel)
+            self.gamma = self.phi_m_last / phim_new
+
+        self.reg.curModel = self.invProb.curModel
+        self.reg.gamma = self.gamma
+        
+        if getattr(self, 'phi_d_last', None) is None:
+            self.phi_d_last = self.invProb.phi_d
+            
+    def endIter(self):
+        # Cool the threshold parameter
+        if getattr(self, 'factor', None) is not None:
+            eps = self.reg.eps / self.factor
+
+            if getattr(self, 'eps_min', None) is not None:
+                self.reg.eps = np.max([self.eps_min,eps])
+            else:
+                self.reg.eps = eps
+
+        # Get phi_m at the end of current iteration
+        self.phi_m_last = self.invProb.phi_m_last
+        
+         # Update the model used for the IRLS weights
+        self.reg.curModel = self.invProb.curModel
+
+        # Update the pre-conditioner
+        diagA = np.sum(self.prob.G**2.,axis=0) + self.invProb.beta*(self.reg.W.T*self.reg.W).diagonal() * (self.reg.mapping * np.ones(self.reg.curModel.size))**2.
+        PC     = Utils.sdiag(diagA**-1.)
+        self.opt.approxHinv = PC
+         
+        # Temporarely set gamma to 1.
+        self.reg.gamma = 1.
+        
+        # Compute change in model objective function and update scaling
+        phim_new = self.reg.eval(self.invProb.curModel)
+        self.reg.gamma = self.phi_m_last / phim_new
+        
+        # TO DO: Re-scale beta if too much change in misfit
+        self.invProb.beta = self.invProb.beta * self.phi_d_last / self.invProb.phi_d
+        
+#==============================================================================
+#          import pylab as plt
+#          plt.figure()
+#          ax = plt.subplot(221)
+#          self.prob.mesh.plotSlice(self.invProb.curModel, ax = ax, normal = 'Z', ind=-5, clim = (0, 0.005))
+#==============================================================================

SimPEG/EM/Base.py

@@ -2,14 +2,14 @@ from SimPEG import Survey, Problem, Utils, Models, Maps, PropMaps, np, sp, Solve
 from scipy.constants import mu_0
 
 class EMPropMap(Maps.PropMap):
-    """ 
+    """
         Property Map for EM Problems. The electrical conductivity (\\(\\sigma\\)) is the default inversion property, and the default value of the magnetic permeability is that of free space (\\(\\mu = 4\\pi\\times 10^{-7} \\) H/m)
     """
 
     sigma = Maps.Property("Electrical Conductivity", defaultInvProp = True, propertyLink=('rho',Maps.ReciprocalMap))
     mu = Maps.Property("Inverse Magnetic Permeability", defaultVal = mu_0, propertyLink=('mui',Maps.ReciprocalMap))
 
-    rho = Maps.Property("Electrical Resistivity", propertyLink=('sigma', Maps.ReciprocalMap)) 
+    rho = Maps.Property("Electrical Resistivity", propertyLink=('sigma', Maps.ReciprocalMap))
     mui = Maps.Property("Inverse Magnetic Permeability", defaultVal = 1./mu_0, propertyLink=('mu', Maps.ReciprocalMap))
 
 
@@ -21,7 +21,7 @@ class BaseEMProblem(Problem.BaseProblem):
 
     surveyPair = Survey.BaseSurvey
     dataPair = Survey.Data
-    
+
     PropMap = EMPropMap
 
     Solver = SimpegSolver
@@ -51,7 +51,7 @@ class BaseEMProblem(Problem.BaseProblem):
         if self.mapping.muMap is not None or self.mapping.muiMap is not None:
             toDelete += ['_MeMu', '_MeMuI','_MfMui','_MfMuiI']
         return toDelete
-    
+
     @property
     def Me(self):
         """
@@ -71,7 +71,7 @@ class BaseEMProblem(Problem.BaseProblem):
         return self._Mf
 
 
-    # ----- Magnetic Permeability ----- # 
+    # ----- Magnetic Permeability ----- #
     @property
     def MfMui(self):
         """
@@ -109,7 +109,7 @@ class BaseEMProblem(Problem.BaseProblem):
         return self._MeMuI
 
 
-    # ----- Electrical Conductivity ----- # 
+    # ----- Electrical Conductivity ----- #
     #TODO: hardcoded to sigma as the model
     @property
     def MeSigma(self):
@@ -120,18 +120,18 @@ class BaseEMProblem(Problem.BaseProblem):
             self._MeSigma = self.mesh.getEdgeInnerProduct(self.curModel.sigma)
         return self._MeSigma
 
-    # TODO: This should take a vector 
+    # TODO: This should take a vector
     def MeSigmaDeriv(self, u):
         """
             Derivative of MeSigma with respect to the model
-        """ 
+        """
         return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma)(u) * self.curModel.sigmaDeriv
-    
+
 
     @property
     def MeSigmaI(self):
         """
-            Inverse of the edge inner product matrix for \\(\\sigma\\). 
+            Inverse of the edge inner product matrix for \\(\\sigma\\).
         """
         if getattr(self, '_MeSigmaI', None) is None:
             self._MeSigmaI = self.mesh.getEdgeInnerProduct(self.curModel.sigma, invMat=True)
@@ -140,8 +140,8 @@ class BaseEMProblem(Problem.BaseProblem):
     # TODO: This should take a vector
     def MeSigmaIDeriv(self, u):
         """
-            Derivative of :code:`MeSigma` with respect to the model 
-        """ 
+            Derivative of :code:`MeSigma` with respect to the model
+        """
         # TODO: only works for diagonal tensors. getEdgeInnerProductDeriv, invMat=True should be implemented in SimPEG
 
         dMeSigmaI_dI = -self.MeSigmaI**2
@@ -163,7 +163,7 @@ class BaseEMProblem(Problem.BaseProblem):
     # TODO: This should take a vector
     def MfRhoDeriv(self,u):
         """
-            Derivative of :code:`MfRho` with respect to the model. 
+            Derivative of :code:`MfRho` with respect to the model.
         """
         return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) * (-Utils.sdiag(self.curModel.rho**2) * self.curModel.sigmaDeriv)
         # self.curModel.rhoDeriv
@@ -181,6 +181,29 @@ class BaseEMProblem(Problem.BaseProblem):
     # TODO: This should take a vector
     def MfRhoIDeriv(self,u):
         """
-            Derivative of :code:`MfRhoI` with respect to the model. 
+            Derivative of :code:`MfRhoI` with respect to the model.
         """
         return self.mesh.getFaceInnerProductDeriv(self.curModel.rho, invMat=True)(u) * self.curModel.rhoDeriv
+
+class BaseEMSurvey(Survey.BaseSurvey):
+
+    def __init__(self, srcList, **kwargs):
+        # Sort these by frequency
+        self.srcList = srcList
+        Survey.BaseSurvey.__init__(self, **kwargs)
+
+    def eval(self, u):
+        """
+        Project fields to receiver locations
+        :param Fields u: fields object
+        :rtype: numpy.ndarray
+        :return: data
+        """
+        data = Survey.Data(self)
+        for src in self.srcList:
+            for rx in src.rxList:
+                data[src, rx] = rx.eval(src, self.mesh, u)
+        return data
+
+    def evalDeriv(self, u):
+        raise Exception('Use Receivers to project fields deriv.')

SimPEG/EM/FDEM/FDEM.py

@@ -18,9 +18,9 @@ class BaseFDEMProblem(BaseEMProblem):
             {\mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{s_e}}
 
         if using the E-B formulation (:code:`Problem_e`
-        or :code:`Problem_b`). Note that in this case, :math:`\mathbf{s_e}` is an integrated quantity. 
+        or :code:`Problem_b`). Note that in this case, :math:`\mathbf{s_e}` is an integrated quantity.
 
-        If we write Maxwell's equations in terms of 
+        If we write Maxwell's equations in terms of
         \\\(\\\mathbf{h}\\\) and current density \\\(\\\mathbf{j}\\\)
 
         .. math ::
@@ -28,7 +28,7 @@ class BaseFDEMProblem(BaseEMProblem):
             \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{j} + i \omega \mathbf{M_{\mu}^e} \mathbf{h} = \mathbf{s_m} \\\\
             \mathbf{C} \mathbf{h} - \mathbf{j} = \mathbf{s_e}
 
-        if using the H-J formulation (:code:`Problem_j` or :code:`Problem_h`). Note that here, :math:`\mathbf{s_m}` is an integrated quantity. 
+        if using the H-J formulation (:code:`Problem_j` or :code:`Problem_h`). Note that here, :math:`\mathbf{s_m}` is an integrated quantity.
 
         The problem performs the elimination so that we are solving the system for \\\(\\\mathbf{e},\\\mathbf{b},\\\mathbf{j} \\\) or \\\(\\\mathbf{h}\\\)
     """
@@ -36,76 +36,76 @@ class BaseFDEMProblem(BaseEMProblem):
     surveyPair = SurveyFDEM
     fieldsPair = Fields
 
-    def fields(self, m=None):
+    def fields(self, m):
         """
         Solve the forward problem for the fields.
-        
+
         :param numpy.array m: inversion model (nP,)
         :rtype numpy.array:
-        :return F: forward solution
+        :return f: forward solution
         """
 
         self.curModel = m
-        F = self.fieldsPair(self.mesh, self.survey)
+        f = self.fieldsPair(self.mesh, self.survey)
 
         for freq in self.survey.freqs:
             A = self.getA(freq)
             rhs = self.getRHS(freq)
             Ainv = self.Solver(A, **self.solverOpts)
-            sol = Ainv * rhs
+            u = Ainv * rhs
             Srcs = self.survey.getSrcByFreq(freq)
-            F[Srcs, self._solutionType] = sol
+            f[Srcs, self._solutionType] = u
             Ainv.clean()
-        return F
+        return f
 
-    def Jvec(self, m, v, u=None):
+    def Jvec(self, m, v, f=None):
         """
         Sensitivity times a vector.
 
         :param numpy.array m: inversion model (nP,)
         :param numpy.array v: vector which we take sensitivity product with (nP,)
-        :param SimPEG.EM.FDEM.Fields u: fields object 
+        :param SimPEG.EM.FDEM.Fields u: fields object
         :rtype numpy.array:
-        :return: Jv (ndata,) 
+        :return: Jv (ndata,)
         """
 
-        if u is None:
-           u = self.fields(m)
+        if f is None:
+           f = self.fields(m)
 
         self.curModel = m
 
         Jv = self.dataPair(self.survey)
 
         for freq in self.survey.freqs:
-            A = self.getA(freq) 
-            Ainv = self.Solver(A, **self.solverOpts)
+            A = self.getA(freq)
+            Ainv = self.Solver(A, **self.solverOpts) # create the concept of Ainv (actually a solve)
 
             for src in self.survey.getSrcByFreq(freq):
-                u_src = u[src, self._solutionType]
+                u_src = f[src, self._solutionType]
                 dA_dm_v = self.getADeriv(freq, u_src, v)
-                dRHS_dm_v = self.getRHSDeriv(freq, src, v) 
+                dRHS_dm_v = self.getRHSDeriv(freq, src, v)
                 du_dm_v = Ainv * ( - dA_dm_v + dRHS_dm_v )
-                
+
                 for rx in src.rxList:
-                    df_dmFun = getattr(u, '_%sDeriv'%rx.projField, None)
+                    df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
                     df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
-                    Jv[src, rx] = rx.evalDeriv(src, self.mesh, u, df_dm_v)
+                    Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
             Ainv.clean()
         return Utils.mkvc(Jv)
 
-    def Jtvec(self, m, v, u=None):
+    def Jtvec(self, m, v, f=None):
         """
         Sensitivity transpose times a vector
 
         :param numpy.array m: inversion model (nP,)
         :param numpy.array v: vector which we take adjoint product with (nP,)
-        :param SimPEG.EM.FDEM.Fields u: fields object 
+        :param SimPEG.EM.FDEM.Fields u: fields object
         :rtype numpy.array:
-        :return: Jv (ndata,) 
+        :return: Jv (ndata,)
         """
 
-        if u is None:
-            u = self.fields(m)
+        if f is None:
+            f = self.fields(m)
 
         self.curModel = m
 
@@ -120,12 +120,12 @@ class BaseFDEMProblem(BaseEMProblem):
             ATinv = self.Solver(AT, **self.solverOpts)
 
             for src in self.survey.getSrcByFreq(freq):
-                u_src = u[src, self._solutionType]
+                u_src = f[src, self._solutionType]
 
                 for rx in src.rxList:
-                    PTv = rx.evalDeriv(src, self.mesh, u, v[src, rx], adjoint=True) # wrt u, need possibility wrt m
+                    PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m
 
-                    df_duTFun = getattr(u, '_%sDeriv'%rx.projField, None)
+                    df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
                     df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True)
 
                     ATinvdf_duT = ATinv * df_duT
@@ -144,7 +144,7 @@ class BaseFDEMProblem(BaseEMProblem):
                         Jtv += - np.array(df_dmT, dtype=complex).real
                     else:
                         raise Exception('Must be real or imag')
-            
+
             ATinv.clean()
 
         return Utils.mkvc(Jtv)
@@ -154,23 +154,23 @@ class BaseFDEMProblem(BaseEMProblem):
         Evaluates the sources for a given frequency and puts them in matrix form
 
         :param float freq: Frequency
-        :rtype: (numpy.ndarray, numpy.ndarray)  
-        :return: S_m, S_e (nE or nF, nSrc)
+        :rtype: (numpy.ndarray, numpy.ndarray)
+        :return: s_m, s_e (nE or nF, nSrc)
         """
         Srcs = self.survey.getSrcByFreq(freq)
         if self._formulation is 'EB':
-            S_m = np.zeros((self.mesh.nF,len(Srcs)), dtype=complex)
-            S_e = np.zeros((self.mesh.nE,len(Srcs)), dtype=complex)
+            s_m = np.zeros((self.mesh.nF,len(Srcs)), dtype=complex)
+            s_e = np.zeros((self.mesh.nE,len(Srcs)), dtype=complex)
         elif self._formulation is 'HJ':
-            S_m = np.zeros((self.mesh.nE,len(Srcs)), dtype=complex)
-            S_e = np.zeros((self.mesh.nF,len(Srcs)), dtype=complex)
+            s_m = np.zeros((self.mesh.nE,len(Srcs)), dtype=complex)
+            s_e = np.zeros((self.mesh.nF,len(Srcs)), dtype=complex)
 
         for i, src in enumerate(Srcs):
             smi, sei = src.eval(self)
-            S_m[:,i] = S_m[:,i] + smi
-            S_e[:,i] = S_e[:,i] + sei
+            s_m[:,i] = s_m[:,i] + smi
+            s_e[:,i] = s_e[:,i] + sei
 
-        return S_m, S_e
+        return s_m, s_e
 
 
 ##########################################################################################
@@ -207,7 +207,7 @@ class Problem_e(BaseFDEMProblem):
     def getA(self, freq):
         """
         System matrix
-        
+
         .. math ::
             \mathbf{A} = \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{C} + i \omega \mathbf{M^e_{\sigma}}
 
@@ -230,12 +230,12 @@ class Problem_e(BaseFDEMProblem):
         .. math ::
             \\frac{\mathbf{A}(\mathbf{m}) \mathbf{v}}{d \mathbf{m}} = i \omega \\frac{d \mathbf{M^e_{\sigma}}\mathbf{v} }{d\mathbf{m}}
 
-        :param float freq: frequency 
-        :param numpy.ndarray u: solution vector (nE,) 
+        :param float freq: frequency
+        :param numpy.ndarray u: solution vector (nE,)
         :param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint
         :param bool adjoint: adjoint?
         :rtype: numpy.ndarray
-        :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,) 
+        :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
         """
 
         dsig_dm = self.curModel.sigmaDeriv
@@ -248,25 +248,25 @@ class Problem_e(BaseFDEMProblem):
 
     def getRHS(self, freq):
         """
-        Right hand side for the system 
+        Right hand side for the system
 
         .. math ::
             \mathbf{RHS} = \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f}\mathbf{s_m} -i\omega\mathbf{M_e}\mathbf{s_e}
 
         :param float freq: Frequency
-        :rtype: numpy.ndarray 
+        :rtype: numpy.ndarray
         :return: RHS (nE, nSrc)
         """
 
-        S_m, S_e = self.getSourceTerm(freq)
+        s_m, s_e = self.getSourceTerm(freq)
         C = self.mesh.edgeCurl
         MfMui = self.MfMui
 
-        return C.T * (MfMui * S_m) -1j * omega(freq) * S_e
+        return C.T * (MfMui * s_m) -1j * omega(freq) * s_e
 
     def getRHSDeriv(self, freq, src, v, adjoint=False):
         """
-        Derivative of the right hand side with respect to the model 
+        Derivative of the right hand side with respect to the model
 
         :param float freq: frequency
         :param SimPEG.EM.FDEM.Src src: FDEM source
@@ -278,14 +278,14 @@ class Problem_e(BaseFDEMProblem):
 
         C = self.mesh.edgeCurl
         MfMui = self.MfMui
-        S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint)
+        s_mDeriv, s_eDeriv = src.evalDeriv(self, adjoint=adjoint)
 
         if adjoint:
             dRHS = MfMui * (C * v)
-            return S_mDeriv(dRHS) - 1j * omega(freq) * S_eDeriv(v)
+            return s_mDeriv(dRHS) - 1j * omega(freq) * s_eDeriv(v)
 
         else:
-            return C.T * (MfMui * S_mDeriv(v)) -1j * omega(freq) * S_eDeriv(v)
+            return C.T * (MfMui * s_mDeriv(v)) -1j * omega(freq) * s_eDeriv(v)
 
 
 class Problem_b(BaseFDEMProblem):
@@ -346,12 +346,12 @@ class Problem_b(BaseFDEMProblem):
         .. math ::
             \\frac{\mathbf{A}(\mathbf{m}) \mathbf{v}}{d \mathbf{m}} = \mathbf{C} \\frac{\mathbf{M^e_{\sigma}} \mathbf{v}}{d\mathbf{m}}
 
-        :param float freq: frequency 
-        :param numpy.ndarray u: solution vector (nF,) 
+        :param float freq: frequency
+        :param numpy.ndarray u: solution vector (nF,)
         :param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint
         :param bool adjoint: adjoint?
         :rtype: numpy.ndarray
-        :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)  
+        :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
         """
 
         MfMui = self.MfMui
@@ -373,21 +373,21 @@ class Problem_b(BaseFDEMProblem):
 
     def getRHS(self, freq):
         """
-        Right hand side for the system 
+        Right hand side for the system
 
         .. math ::
             \mathbf{RHS} = \mathbf{s_m} + \mathbf{M^e_{\sigma}}^{-1}\mathbf{s_e}
 
         :param float freq: Frequency
-        :rtype: numpy.ndarray 
+        :rtype: numpy.ndarray
         :return: RHS (nE, nSrc)
         """
 
-        S_m, S_e = self.getSourceTerm(freq)
+        s_m, s_e = self.getSourceTerm(freq)
         C = self.mesh.edgeCurl
         MeSigmaI = self.MeSigmaI
 
-        RHS = S_m + C * ( MeSigmaI * S_e )
+        RHS = s_m + C * ( MeSigmaI * s_e )
 
         if self._makeASymmetric is True:
             MfMui = self.MfMui
@@ -408,21 +408,21 @@ class Problem_b(BaseFDEMProblem):
         """
 
         C = self.mesh.edgeCurl
-        S_m, S_e = src.eval(self)
+        s_m, s_e = src.eval(self)
         MfMui = self.MfMui
 
         if self._makeASymmetric and adjoint:
             v = self.MfMui * v
 
-        MeSigmaIDeriv = self.MeSigmaIDeriv(S_e)
-        S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint)
+        MeSigmaIDeriv = self.MeSigmaIDeriv(s_e)
+        s_mDeriv, s_eDeriv = src.evalDeriv(self, adjoint=adjoint)
 
         if not adjoint:
             RHSderiv = C * (MeSigmaIDeriv * v)
-            SrcDeriv = S_mDeriv(v) + C * (self.MeSigmaI * S_eDeriv(v))
+            SrcDeriv = s_mDeriv(v) + C * (self.MeSigmaI * s_eDeriv(v))
         elif adjoint:
             RHSderiv = MeSigmaIDeriv.T * (C.T * v)
-            SrcDeriv = S_mDeriv(v) + self.MeSigmaI.T * (C.T * S_eDeriv(v))
+            SrcDeriv = s_mDeriv(v) + self.MeSigmaI.T * (C.T * s_eDeriv(v))
 
         if self._makeASymmetric is True and not adjoint:
             return MfMui.T * (SrcDeriv + RHSderiv)
@@ -497,12 +497,12 @@ class Problem_j(BaseFDEMProblem):
 
             \\frac{\mathbf{A(\sigma)} \mathbf{v}}{d \mathbf{m}} = \mathbf{C}  \mathbf{M^e_{mu^{-1}}} \mathbf{C^{\\top}} \\frac{d \mathbf{M^f_{\sigma^{-1}}}\mathbf{v} }{d \mathbf{m}}
 
-        :param float freq: frequency 
-        :param numpy.ndarray u: solution vector (nF,) 
+        :param float freq: frequency
+        :param numpy.ndarray u: solution vector (nF,)
         :param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint
         :param bool adjoint: adjoint?
         :rtype: numpy.ndarray
-        :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)  
+        :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
         """
 
         MeMuI = self.MeMuI
@@ -522,7 +522,7 @@ class Problem_j(BaseFDEMProblem):
 
     def getRHS(self, freq):
         """
-        Right hand side for the system 
+        Right hand side for the system
 
         .. math ::
 
@@ -533,11 +533,11 @@ class Problem_j(BaseFDEMProblem):
         :return: RHS
         """
 
-        S_m, S_e = self.getSourceTerm(freq)
+        s_m, s_e = self.getSourceTerm(freq)
         C = self.mesh.edgeCurl
         MeMuI = self.MeMuI
 
-        RHS = C * (MeMuI * S_m) - 1j * omega(freq) * S_e
+        RHS = C * (MeMuI * s_m) - 1j * omega(freq) * s_e
         if self._makeASymmetric is True:
             MfRho = self.MfRho
             return MfRho.T*RHS
@@ -546,7 +546,7 @@ class Problem_j(BaseFDEMProblem):
 
     def getRHSDeriv(self, freq, src, v, adjoint=False):
         """
-        Derivative of the right hand side with respect to the model 
+        Derivative of the right hand side with respect to the model
 
         :param float freq: frequency
         :param SimPEG.EM.FDEM.Src src: FDEM source
@@ -558,16 +558,16 @@ class Problem_j(BaseFDEMProblem):
 
         C = self.mesh.edgeCurl
         MeMuI = self.MeMuI
-        S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint)
+        s_mDeriv, s_eDeriv = src.evalDeriv(self, adjoint=adjoint)
 
         if adjoint:
             if self._makeASymmetric:
                 MfRho = self.MfRho
                 v = MfRho*v
-            return S_mDeriv(MeMuI.T * (C.T * v)) - 1j * omega(freq) * S_eDeriv(v)
+            return s_mDeriv(MeMuI.T * (C.T * v)) - 1j * omega(freq) * s_eDeriv(v)
 
         else:
-            RHSDeriv = C * (MeMuI * S_mDeriv(v)) - 1j * omega(freq) * S_eDeriv(v)
+            RHSDeriv = C * (MeMuI * s_mDeriv(v)) - 1j * omega(freq) * s_eDeriv(v)
 
             if self._makeASymmetric:
                 MfRho = self.MfRho
@@ -626,12 +626,12 @@ class Problem_h(BaseFDEMProblem):
         .. math::
             \\frac{\mathbf{A}(\mathbf{m}) \mathbf{v}}{d \mathbf{m}} = \mathbf{C}^{\\top}\\frac{d \mathbf{M^f_{\\rho}}\mathbf{v} }{d\mathbf{m}}
 
-        :param float freq: frequency 
-        :param numpy.ndarray u: solution vector (nE,) 
+        :param float freq: frequency
+        :param numpy.ndarray u: solution vector (nE,)
         :param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint
         :param bool adjoint: adjoint?
         :rtype: numpy.ndarray
-        :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)  
+        :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,)
         """
 
         MeMu = self.MeMu
@@ -644,26 +644,26 @@ class Problem_h(BaseFDEMProblem):
 
     def getRHS(self, freq):
         """
-        Right hand side for the system 
+        Right hand side for the system
 
         .. math ::
 
             \mathbf{RHS} = \mathbf{M^e} \mathbf{s_m} + \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{s_e}
 
         :param float freq: Frequency
-        :rtype: numpy.ndarray 
+        :rtype: numpy.ndarray
         :return: RHS (nE, nSrc)
         """
 
-        S_m, S_e = self.getSourceTerm(freq)
+        s_m, s_e = self.getSourceTerm(freq)
         C = self.mesh.edgeCurl
         MfRho  = self.MfRho
 
-        return S_m + C.T * ( MfRho * S_e )
+        return s_m + C.T * ( MfRho * s_e )
 
     def getRHSDeriv(self, freq, src, v, adjoint=False):
         """
-        Derivative of the right hand side with respect to the model 
+        Derivative of the right hand side with respect to the model
 
         :param float freq: frequency
         :param SimPEG.EM.FDEM.Src src: FDEM source
@@ -673,17 +673,17 @@ class Problem_h(BaseFDEMProblem):
         :return: product of rhs deriv with a vector
         """
 
-        _, S_e = src.eval(self)
+        _, s_e = src.eval(self)
         C = self.mesh.edgeCurl
         MfRho  = self.MfRho
 
-        MfRhoDeriv = self.MfRhoDeriv(S_e)
+        MfRhoDeriv = self.MfRhoDeriv(s_e)
         if not adjoint:
             RHSDeriv = C.T * (MfRhoDeriv * v)
         elif adjoint:
             RHSDeriv = MfRhoDeriv.T * (C * v)
 
-        S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint)
+        s_mDeriv, s_eDeriv = src.evalDeriv(self, adjoint=adjoint)
 
-        return RHSDeriv + S_mDeriv(v) + C.T * (MfRho * S_eDeriv(v))
+        return RHSDeriv + s_mDeriv(v) + C.T * (MfRho * s_eDeriv(v))
 

SimPEG/EM/FDEM/FieldsFDEM.py

@@ -8,7 +8,7 @@ from SimPEG.Utils import Zero, Identity, sdiag
 
 class Fields(SimPEG.Problem.Fields):
     """
-    
+
     Fancy Field Storage for a FDEM survey. Only one field type is stored for
     each problem, the rest are computed. The fields obejct acts like an array and is indexed by
 
@@ -34,56 +34,56 @@ class Fields(SimPEG.Problem.Fields):
 
     def _e(self, solution, srcList):
         """
-        Total electric field is sum of primary and secondary 
-        
+        Total electric field is sum of primary and secondary
+
         :param numpy.ndarray solution: field we solved for
         :param list srcList: list of sources
         :rtype: numpy.ndarray
         :return: total electric field
         """
-        if getattr(self, '_ePrimary', None) is None or getattr(self, '_eSecondary', None) is None: 
+        if getattr(self, '_ePrimary', None) is None or getattr(self, '_eSecondary', None) is None:
             raise NotImplementedError ('Getting e from %s is not implemented' %self.knownFields.keys()[0])
 
         return self._ePrimary(solution,srcList) + self._eSecondary(solution,srcList)
 
     def _b(self, solution, srcList):
         """
-        Total magnetic flux density is sum of primary and secondary 
-        
+        Total magnetic flux density is sum of primary and secondary
+
         :param numpy.ndarray solution: field we solved for
         :param list srcList: list of sources
         :rtype: numpy.ndarray
-        :return: total magnetic flux density 
+        :return: total magnetic flux density
         """
-        if getattr(self, '_bPrimary', None) is None or getattr(self, '_bSecondary', None) is None: 
+        if getattr(self, '_bPrimary', None) is None or getattr(self, '_bSecondary', None) is None:
             raise NotImplementedError ('Getting b from %s is not implemented' %self.knownFields.keys()[0])
 
         return self._bPrimary(solution, srcList) + self._bSecondary(solution, srcList)
 
     def _h(self, solution, srcList):
         """
-        Total magnetic field is sum of primary and secondary 
-        
+        Total magnetic field is sum of primary and secondary
+
         :param numpy.ndarray solution: field we solved for
         :param list srcList: list of sources
         :rtype: numpy.ndarray
         :return: total magnetic field
         """
-        if getattr(self, '_hPrimary', None) is None or getattr(self, '_hSecondary', None) is None: 
+        if getattr(self, '_hPrimary', None) is None or getattr(self, '_hSecondary', None) is None:
             raise NotImplementedError ('Getting h from %s is not implemented' %self.knownFields.keys()[0])
 
         return self._hPrimary(solution, srcList) + self._hSecondary(solution, srcList)
 
     def _j(self, solution, srcList):
         """
-        Total current density is sum of primary and secondary 
-        
+        Total current density is sum of primary and secondary
+
         :param numpy.ndarray solution: field we solved for
         :param list srcList: list of sources
         :rtype: numpy.ndarray
-        :return: total current density 
+        :return: total current density
         """
-        if getattr(self, '_jPrimary', None) is None or getattr(self, '_jSecondary', None) is None: 
+        if getattr(self, '_jPrimary', None) is None or getattr(self, '_jSecondary', None) is None:
             raise NotImplementedError ('Getting j from %s is not implemented' %self.knownFields.keys()[0])
 
         return self._jPrimary(solution, srcList) + self._jSecondary(solution, srcList)
@@ -99,7 +99,7 @@ class Fields(SimPEG.Problem.Fields):
         :rtype: numpy.ndarray
         :return: derivative times a vector (or tuple for adjoint)
         """
-        if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None: 
+        if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None:
             raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0])
 
         if adjoint:
@@ -117,12 +117,12 @@ class Fields(SimPEG.Problem.Fields):
         :rtype: numpy.ndarray
         :return: derivative times a vector (or tuple for adjoint)
         """
-        if getattr(self, '_bDeriv_u', None) is None or getattr(self, '_bDeriv_m', None) is None: 
+        if getattr(self, '_bDeriv_u', None) is None or getattr(self, '_bDeriv_m', None) is None:
             raise NotImplementedError ('Getting bDerivs from %s is not implemented' %self.knownFields.keys()[0])
 
         if adjoint:
             return self._bDeriv_u(src, v, adjoint), self._bDeriv_m(src, v, adjoint)
-        return np.array(self._bDeriv_u(src, du_dm_v, adjoint) + self._bDeriv_m(src, v, adjoint), dtype = complex) 
+        return np.array(self._bDeriv_u(src, du_dm_v, adjoint) + self._bDeriv_m(src, v, adjoint), dtype = complex)
 
     def _hDeriv(self, src, du_dm_v, v, adjoint = False):
         """
@@ -135,10 +135,10 @@ class Fields(SimPEG.Problem.Fields):
         :rtype: numpy.ndarray
         :return: derivative times a vector (or tuple for adjoint)
         """
-        if getattr(self, '_hDeriv_u', None) is None or getattr(self, '_hDeriv_m', None) is None: 
+        if getattr(self, '_hDeriv_u', None) is None or getattr(self, '_hDeriv_m', None) is None:
             raise NotImplementedError ('Getting hDerivs from %s is not implemented' %self.knownFields.keys()[0])
 
-        if adjoint: 
+        if adjoint:
             return self._hDeriv_u(src, v, adjoint), self._hDeriv_m(src, v, adjoint)
         return np.array(self._hDeriv_u(src, du_dm_v, adjoint) + self._hDeriv_m(src, v, adjoint), dtype = complex)
 
@@ -153,7 +153,7 @@ class Fields(SimPEG.Problem.Fields):
         :rtype: numpy.ndarray
         :return: derivative times a vector (or tuple for adjoint)
         """
-        if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None: 
+        if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None:
             raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0])
 
         if adjoint:
@@ -162,10 +162,10 @@ class Fields(SimPEG.Problem.Fields):
 
 class Fields_e(Fields):
     """
-    Fields object for Problem_e. 
+    Fields object for Problem_e.
 
     :param Mesh mesh: mesh
-    :param Survey survey: survey 
+    :param Survey survey: survey
     """
 
     knownFields = {'eSolution':'E'}
@@ -233,9 +233,9 @@ class Fields_e(Fields):
 
     def _eDeriv_u(self, src, v, adjoint = False):
         """
-        Partial derivative of the total electric field with respect to the thing we 
+        Partial derivative of the total electric field with respect to the thing we
         solved for.
-        
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -247,8 +247,8 @@ class Fields_e(Fields):
 
     def _eDeriv_m(self, src, v, adjoint = False):
         """
-        Partial derivative of the total electric field with respect to the inversion model. Here, we assume that the primary does not depend on the model. Note that this also includes derivative contributions from the sources. 
-        
+        Partial derivative of the total electric field with respect to the inversion model. Here, we assume that the primary does not depend on the model. Note that this also includes derivative contributions from the sources.
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -289,14 +289,14 @@ class Fields_e(Fields):
         b = (C * eSolution)
         for i, src in enumerate(srcList):
             b[:,i] *= - 1./(1j*omega(src.freq))
-            S_m, _ = src.eval(self.prob)
-            b[:,i] = b[:,i]+ 1./(1j*omega(src.freq)) * S_m
+            s_m, _ = src.eval(self.prob)
+            b[:,i] = b[:,i]+ 1./(1j*omega(src.freq)) * s_m
         return b
 
     def _bDeriv_u(self, src, du_dm_v, adjoint = False):
         """
         Derivative of the magnetic flux density with respect to the thing we solved for
-        
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray du_dm_v: vector to take product with
         :param bool adjoint: adjoint?
@@ -312,8 +312,8 @@ class Fields_e(Fields):
 
     def _bDeriv_m(self, src, v, adjoint = False):
         """
-        Derivative of the magnetic flux density with respect to the inversion model. 
-        
+        Derivative of the magnetic flux density with respect to the inversion model.
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -321,8 +321,8 @@ class Fields_e(Fields):
         :return: product of the magnetic flux density derivative with respect to the inversion model with a vector
         """
 
-        S_mDeriv, _ = src.evalDeriv(self.prob, v, adjoint)
-        return 1./(1j * omega(src.freq)) * S_mDeriv
+        s_mDeriv, _ = src.evalDeriv(self.prob, v, adjoint)
+        return 1./(1j * omega(src.freq)) * s_mDeriv
 
     def _j(self,  eSolution, srcList):
         """
@@ -341,7 +341,7 @@ class Fields_e(Fields):
     def _jDeriv_u(self, src, du_dm_v, adjoint = False):
         """
         Derivative of the current density with respect to the thing we solved for
-        
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray du_dm_v: vector to take product with
         :param bool adjoint: adjoint?
@@ -351,15 +351,15 @@ class Fields_e(Fields):
         n = int(self._aveE2CCV.shape[0] / self._nC) # number of components (instead of checking if cyl or not)
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
 
-        if adjoint: 
+        if adjoint:
             return  self._eDeriv_u(src, self._MeSigma.T * (self._aveE2CCV.T * (VI.T * du_dm_v) ), adjoint = adjoint)
         return VI * (self._aveE2CCV * (self._MeSigma * (self._eDeriv_u(src, du_dm_v, adjoint=adjoint) ) ) )
-        
+
 
     def _jDeriv_m(self, src, v, adjoint = False):
         """
-        Derivative of the current density with respect to the inversion model. 
-        
+        Derivative of the current density with respect to the inversion model.
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -373,7 +373,7 @@ class Fields_e(Fields):
         if adjoint:
             return self._MeSigmaDeriv(e).T * (self._aveE2CCV.T * (VI.T * v)) + self._eDeriv_m(src, self._aveE2CCV.T * (VI.T * v), adjoint=adjoint)
         return VI * (self._aveE2CCV * ( self._eDeriv_m(src, v, adjoint=adjoint) + self._MeSigmaDeriv(e) * v))
-        
+
 
 
     def _h(self, eSolution, srcList):
@@ -393,7 +393,7 @@ class Fields_e(Fields):
     def _hDeriv_u(self, src, du_dm_v, adjoint = False):
         """
         Derivative of the magnetic field with respect to the thing we solved for
-        
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray du_dm_v: vector to take product with
         :param bool adjoint: adjoint?
@@ -409,8 +409,8 @@ class Fields_e(Fields):
 
     def _hDeriv_m(self, src, v, adjoint = False):
         """
-        Derivative of the magnetic field with respect to the inversion model. 
-        
+        Derivative of the magnetic field with respect to the inversion model.
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -428,10 +428,10 @@ class Fields_e(Fields):
 
 class Fields_b(Fields):
     """
-    Fields object for Problem_b. 
+    Fields object for Problem_b.
 
     :param Mesh mesh: mesh
-    :param Survey survey: survey 
+    :param Survey survey: survey
     """
 
     knownFields = {'bSolution':'F'}
@@ -506,9 +506,9 @@ class Fields_b(Fields):
 
     def _bDeriv_u(self, src, du_dm_v, adjoint=False):
         """
-        Partial derivative of the total magnetic flux density with respect to the thing we 
+        Partial derivative of the total magnetic flux density with respect to the thing we
         solved for.
-        
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray du_dm_v: vector to take product with
         :param bool adjoint: adjoint?
@@ -520,8 +520,8 @@ class Fields_b(Fields):
 
     def _bDeriv_m(self, src, v, adjoint=False):
         """
-        Partial derivative of the total magnetic flux density with respect to the inversion model. Here, we assume that the primary does not depend on the model. Note that this also includes derivative contributions from the sources. 
-        
+        Partial derivative of the total magnetic flux density with respect to the inversion model. Here, we assume that the primary does not depend on the model. Note that this also includes derivative contributions from the sources.
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -560,15 +560,15 @@ class Fields_b(Fields):
 
         e = ( self._edgeCurl.T * ( self._MfMui * bSolution))
         for i,src in enumerate(srcList):
-            _,S_e = src.eval(self.prob)
-            e[:,i] = e[:,i] + - S_e
+            _,s_e = src.eval(self.prob)
+            e[:,i] = e[:,i] + - s_e
 
         return self._MeSigmaI * e
 
     def _eDeriv_u(self, src, du_dm_v, adjoint=False):
         """
         Derivative of the electric field with respect to the thing we solved for
-        
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -583,8 +583,8 @@ class Fields_b(Fields):
 
     def _eDeriv_m(self, src, v, adjoint=False):
         """
-        Derivative of the electric field with respect to the inversion model 
-        
+        Derivative of the electric field with respect to the inversion model
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -593,15 +593,15 @@ class Fields_b(Fields):
         """
 
         bSolution = Utils.mkvc(self[src, 'bSolution'])
-        _,S_e = src.eval(self.prob)
+        _,s_e = src.eval(self.prob)
 
-        w = -S_e + self._edgeCurl.T * (self._MfMui * bSolution)
-        _, S_eDeriv = src.evalDeriv(self.prob, v, adjoint)
+        w = -s_e + self._edgeCurl.T * (self._MfMui * bSolution)
+        _, s_eDeriv = src.evalDeriv(self.prob, v, adjoint)
 
 
         if adjoint:
-            return self._MeSigmaIDeriv(w).T * v - self._MeSigmaI.T * S_eDeriv
-        return  self._MeSigmaIDeriv(w) * v - self._MeSigmaI * S_eDeriv
+            return self._MeSigmaIDeriv(w).T * v - self._MeSigmaI.T * s_eDeriv
+        return  self._MeSigmaIDeriv(w) * v - self._MeSigmaI * s_eDeriv
 
     def _j(self, bSolution, srcList):
         """
@@ -617,13 +617,13 @@ class Fields_b(Fields):
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
 
         return VI * (self._aveE2CCV * ( self._MeSigma * self._e(bSolution,srcList ) ) )
-        
+
 
     def _jDeriv_u(self, src, du_dm_v, adjoint=False):
         """
-        Partial derivative of the current density with respect to the thing we 
+        Partial derivative of the current density with respect to the thing we
         solved for.
-        
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray du_dm_v: vector to take product with
         :param bool adjoint: adjoint?
@@ -639,8 +639,8 @@ class Fields_b(Fields):
 
     def _jDeriv_m(self, src, v, adjoint=False):
         """
-        Derivative of the current density with respect to the inversion model 
-        
+        Derivative of the current density with respect to the inversion model
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -664,9 +664,9 @@ class Fields_b(Fields):
 
     def _hDeriv_u(self, src, du_dm_v, adjoint=False):
         """
-        Partial derivative of the magnetic field with respect to the thing we 
+        Partial derivative of the magnetic field with respect to the thing we
         solved for.
-        
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray du_dm_v: vector to take product with
         :param bool adjoint: adjoint?
@@ -682,8 +682,8 @@ class Fields_b(Fields):
 
     def _hDeriv_m(self, src, v, adjoint=False):
         """
-        Derivative of the magnetic field with respect to the inversion model 
-        
+        Derivative of the magnetic field with respect to the inversion model
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -695,10 +695,10 @@ class Fields_b(Fields):
 
 class Fields_j(Fields):
     """
-    Fields object for Problem_j. 
+    Fields object for Problem_j.
 
     :param Mesh mesh: mesh
-    :param Survey survey: survey 
+    :param Survey survey: survey
     """
 
     knownFields = {'jSolution':'F'}
@@ -769,12 +769,12 @@ class Fields_j(Fields):
 
     def _j(self, jSolution, srcList):
         """
-        Total current density is sum of primary and secondary 
-        
+        Total current density is sum of primary and secondary
+
         :param numpy.ndarray jSolution: field we solved for
         :param list srcList: list of sources
         :rtype: numpy.ndarray
-        :return: total current density 
+        :return: total current density
         """
 
         return self._jPrimary(jSolution, srcList) + self._jSecondary(jSolution, srcList)
@@ -782,9 +782,9 @@ class Fields_j(Fields):
 
     def _jDeriv_u(self, src, du_dm_v, adjoint=False):
         """
-        Partial derivative of the total current density with respect to the thing we 
+        Partial derivative of the total current density with respect to the thing we
         solved for.
-        
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -797,8 +797,8 @@ class Fields_j(Fields):
 
     def _jDeriv_m(self, src, v, adjoint=False):
         """
-        Partial derivative of the total current density with respect to the inversion model. Here, we assume that the primary does not depend on the model. Note that this also includes derivative contributions from the sources. 
-        
+        Partial derivative of the total current density with respect to the inversion model. Here, we assume that the primary does not depend on the model. Note that this also includes derivative contributions from the sources.
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -837,15 +837,15 @@ class Fields_j(Fields):
         h = (self._edgeCurl.T * (self._MfRho * jSolution) )
         for i, src in enumerate(srcList):
             h[:,i] *= -1./(1j*omega(src.freq))
-            S_m,_ = src.eval(self.prob)
-            h[:,i] = h[:,i] + 1./(1j*omega(src.freq)) * (S_m)
+            s_m,_ = src.eval(self.prob)
+            h[:,i] = h[:,i] + 1./(1j*omega(src.freq)) * (s_m)
         return self._MeMuI * h
 
 
     def _hDeriv_u(self, src, du_dm_v, adjoint=False):
         """
         Derivative of the magnetic field with respect to the thing we solved for
-        
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray du_dm_v: vector to take product with
         :param bool adjoint: adjoint?
@@ -856,13 +856,13 @@ class Fields_j(Fields):
         if adjoint:
             return  -1./(1j*omega(src.freq)) * self._MfRho.T * (self._edgeCurl * ( self._MeMuI.T * du_dm_v))
         return  -1./(1j*omega(src.freq)) * self._MeMuI * (self._edgeCurl.T * (self._MfRho * du_dm_v) )
-        
+
 
 
     def _hDeriv_m(self, src, v, adjoint=False):
         """
-        Derivative of the magnetic field with respect to the inversion model 
-        
+        Derivative of the magnetic field with respect to the inversion model
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -875,19 +875,19 @@ class Fields_j(Fields):
         C = self._edgeCurl
         MfRho = self._MfRho
         MfRhoDeriv = self._MfRhoDeriv
-        S_mDeriv,_ = src.evalDeriv(self.prob, adjoint = adjoint)
+        s_mDeriv,_ = src.evalDeriv(self.prob, adjoint = adjoint)
 
         if not adjoint:
             hDeriv_m =  -1./(1j*omega(src.freq)) * MeMuI * (C.T * (MfRhoDeriv(jSolution)*v ) )
-            S_mDeriv = S_mDeriv(v)
-            hDeriv_m = hDeriv_m + 1./(1j*omega(src.freq)) * MeMuI * ( S_mDeriv)
+            s_mDeriv = s_mDeriv(v)
+            hDeriv_m = hDeriv_m + 1./(1j*omega(src.freq)) * MeMuI * ( s_mDeriv)
 
         elif adjoint:
             hDeriv_m =  -1./(1j*omega(src.freq)) * MfRhoDeriv(jSolution).T * ( C * (MeMuI.T * v ) )
 
-            S_mDeriv = S_mDeriv(MeMuI.T * v)
-            hDeriv_m = hDeriv_m + 1./(1j*omega(src.freq)) * S_mDeriv
-        
+            s_mDeriv = s_mDeriv(MeMuI.T * v)
+            hDeriv_m = hDeriv_m + 1./(1j*omega(src.freq)) * s_mDeriv
+
         return hDeriv_m
 
     def _e(self, jSolution, srcList):
@@ -901,12 +901,12 @@ class Fields_j(Fields):
         """
         n = int(self._aveF2CCV.shape[0] / self._nC) # number of components
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
-        return VI * (self._aveF2CCV * (self._MfRho *  self._j(jSolution, srcList))) 
+        return VI * (self._aveF2CCV * (self._MfRho *  self._j(jSolution, srcList)))
 
     def _eDeriv_u(self, src, du_dm_v, adjoint=False):
         """
         Derivative of the electric field with respect to the thing we solved for
-        
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray du_dm_v: vector to take product with
         :param bool adjoint: adjoint?
@@ -921,8 +921,8 @@ class Fields_j(Fields):
 
     def _eDeriv_m(self, src, v, adjoint=False):
         """
-        Derivative of the electric field with respect to the inversion model 
-        
+        Derivative of the electric field with respect to the inversion model
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -943,17 +943,17 @@ class Fields_j(Fields):
         :param numpy.ndarray hSolution: field we solved for
         :param list srcList: list of sources
         :rtype: numpy.ndarray
-        :return: secondary magnetic flux density 
+        :return: secondary magnetic flux density
         """
         n = int(self._aveE2CCV.shape[0] / self._nC) # number of components
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
 
-        return VI * (self._aveE2CCV * ( self._MeMu * self._h(jSolution,srcList)) ) 
+        return VI * (self._aveE2CCV * ( self._MeMu * self._h(jSolution,srcList)) )
 
     def _bDeriv_u(self, src, du_dm_v, adjoint=False):
         """
         Derivative of the magnetic flux density with respect to the thing we solved for
-        
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray du_dm_v: vector to take product with
         :param bool adjoint: adjoint?
@@ -969,8 +969,8 @@ class Fields_j(Fields):
 
     def _bDeriv_m(self, src, v, adjoint=False):
         """
-        Derivative of the magnetic flux density with respect to the inversion model 
-        
+        Derivative of the magnetic flux density with respect to the inversion model
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -980,20 +980,20 @@ class Fields_j(Fields):
         jSolution = self[src,'jSolution']
         n = int(self._aveE2CCV.shape[0] / self._nC) # number of components
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
-        S_mDeriv,_ = src.evalDeriv(self.prob, adjoint = adjoint)
+        s_mDeriv,_ = src.evalDeriv(self.prob, adjoint = adjoint)
 
         if adjoint:
             v = self._aveE2CCV.T * ( VI.T * v)
-            return 1./(1j * omega(src.freq)) * ( S_mDeriv(v) - self._MfRhoDeriv(jSolution).T * (self._edgeCurl * v ))
-        return 1./(1j * omega(src.freq)) * VI * (self._aveE2CCV * ( S_mDeriv(v) - self._edgeCurl.T * ( self._MfRhoDeriv(jSolution) * v ) ) )
+            return 1./(1j * omega(src.freq)) * ( s_mDeriv(v) - self._MfRhoDeriv(jSolution).T * (self._edgeCurl * v ))
+        return 1./(1j * omega(src.freq)) * VI * (self._aveE2CCV * ( s_mDeriv(v) - self._edgeCurl.T * ( self._MfRhoDeriv(jSolution) * v ) ) )
 
 
 class Fields_h(Fields):
     """
-    Fields object for Problem_h. 
+    Fields object for Problem_h.
 
     :param Mesh mesh: mesh
-    :param Survey survey: survey 
+    :param Survey survey: survey
     """
 
     knownFields = {'hSolution':'E'}
@@ -1065,9 +1065,9 @@ class Fields_h(Fields):
 
     def _hDeriv_u(self, src, du_dm_v, adjoint=False):
         """
-        Partial derivative of the total magnetic field with respect to the thing we 
+        Partial derivative of the total magnetic field with respect to the thing we
         solved for.
-        
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray du_dm_v: vector to take product with
         :param bool adjoint: adjoint?
@@ -1079,8 +1079,8 @@ class Fields_h(Fields):
 
     def _hDeriv_m(self, src, v, adjoint=False):
         """
-        Partial derivative of the total magnetic field with respect to the inversion model. Here, we assume that the primary does not depend on the model. Note that this also includes derivative contributions from the sources. 
-        
+        Partial derivative of the total magnetic field with respect to the inversion model. Here, we assume that the primary does not depend on the model. Note that this also includes derivative contributions from the sources.
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -1119,14 +1119,14 @@ class Fields_h(Fields):
 
         j = self._edgeCurl*hSolution
         for i, src in enumerate(srcList):
-            _,S_e = src.eval(self.prob)
-            j[:,i] = j[:,i]+ -S_e
+            _,s_e = src.eval(self.prob)
+            j[:,i] = j[:,i]+ -s_e
         return j
 
     def _jDeriv_u(self, src, du_dm_v, adjoint=False):
         """
         Derivative of the current density with respect to the thing we solved for
-        
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray du_dm_v: vector to take product with
         :param bool adjoint: adjoint?
@@ -1142,8 +1142,8 @@ class Fields_h(Fields):
 
     def _jDeriv_m(self, src, v, adjoint=False):
         """
-        Derivative of the current density with respect to the inversion model. 
-        
+        Derivative of the current density with respect to the inversion model.
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -1151,9 +1151,9 @@ class Fields_h(Fields):
         :return: product of the current density derivative with respect to the inversion model with a vector
         """
 
-        _,S_eDeriv = src.evalDeriv(self.prob, v, adjoint)
-        return -S_eDeriv
-    
+        _,s_eDeriv = src.evalDeriv(self.prob, v, adjoint)
+        return -s_eDeriv
+
     def _e(self, hSolution, srcList):
         """
         Electric field from hSolution
@@ -1165,12 +1165,12 @@ class Fields_h(Fields):
         """
         n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
-        return VI * (self._aveF2CCV * (self._MfRho *  self._j(hSolution, srcList))) 
+        return VI * (self._aveF2CCV * (self._MfRho *  self._j(hSolution, srcList)))
 
     def _eDeriv_u(self, src, du_dm_v, adjoint=False):
         """
         Derivative of the electric field with respect to the thing we solved for
-        
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray du_dm_v: vector to take product with
         :param bool adjoint: adjoint?
@@ -1181,12 +1181,12 @@ class Fields_h(Fields):
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
         if adjoint:
             return self._edgeCurl.T * ( self._MfRho.T * ( self._aveF2CCV.T * ( VI.T * du_dm_v ) ) )
-        return VI * (self._aveF2CCV * (self._MfRho *  self._edgeCurl * du_dm_v )) 
+        return VI * (self._aveF2CCV * (self._MfRho *  self._edgeCurl * du_dm_v ))
 
     def _eDeriv_m(self, src, v, adjoint=False):
         """
-        Derivative of the electric field with respect to the inversion model. 
-        
+        Derivative of the electric field with respect to the inversion model.
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
@@ -1196,7 +1196,7 @@ class Fields_h(Fields):
         hSolution = Utils.mkvc(self[src,'hSolution'])
         n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
-        if adjoint: 
+        if adjoint:
             return ( self._MfRhoDeriv(self._edgeCurl * hSolution).T * ( self._aveF2CCV.T * (VI.T * v) ) )
         return VI * (self._aveF2CCV * (self._MfRhoDeriv(self._edgeCurl * hSolution) * v ))
 
@@ -1207,10 +1207,10 @@ class Fields_h(Fields):
         :param numpy.ndarray hSolution: field we solved for
         :param list srcList: list of sources
         :rtype: numpy.ndarray
-        :return: magnetic flux density 
+        :return: magnetic flux density
         """
         h = self._h(hSolution, srcList)
-        n = int(self._aveE2CCV.shape[0] / self._nC) #number of components 
+        n = int(self._aveE2CCV.shape[0] / self._nC) #number of components
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
 
         return VI * (self._aveE2CCV * (self._MeMu * h))
@@ -1218,14 +1218,14 @@ class Fields_h(Fields):
     def _bDeriv_u(self, src, du_dm_v, adjoint=False):
         """
         Derivative of the magnetic flux density with respect to the thing we solved for
-        
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray du_dm_v: vector to take product with
         :param bool adjoint: adjoint?
         :rtype: numpy.ndarray
         :return: product of the derivative of the magnetic flux density with respect to the field we solved for with a vector
         """
-        n = int(self._aveE2CCV.shape[0] / self._nC) #number of components 
+        n = int(self._aveE2CCV.shape[0] / self._nC) #number of components
         VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
         if adjoint:
             return self._MeMu.T * (self._aveE2CCV.T * ( VI.T * du_dm_v ))
@@ -1233,8 +1233,8 @@ class Fields_h(Fields):
 
     def _bDeriv_m(self, src, v, adjoint=False):
         """
-        Derivative of the magnetic flux density with respect to the inversion model. 
-        
+        Derivative of the magnetic flux density with respect to the inversion model.
+
         :param SimPEG.EM.FDEM.Src src: source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?

SimPEG/EM/FDEM/SrcFDEM.py

@@ -15,22 +15,22 @@ class BaseSrc(Survey.BaseSrc):
     def eval(self, prob):
         """
         Evaluate the source terms.
-        - :math:`S_m` : magnetic source term
-        - :math:`S_e` : electric source term
+        - :math:`s_m` : magnetic source term
+        - :math:`s_e` : electric source term
 
         :param Problem prob: FDEM Problem
         :rtype: (numpy.ndarray, numpy.ndarray)
         :return: tuple with magnetic source term and electric source term
         """
-        S_m = self.S_m(prob)
-        S_e = self.S_e(prob)
-        return S_m, S_e
+        s_m = self.s_m(prob)
+        s_e = self.s_e(prob)
+        return s_m, s_e
 
     def evalDeriv(self, prob, v=None, adjoint=False):
         """
         Derivatives of the source terms with respect to the inversion model
-        - :code:`S_mDeriv` : derivative of the magnetic source term
-        - :code:`S_eDeriv` : derivative of the electric source term
+        - :code:`s_mDeriv` : derivative of the magnetic source term
+        - :code:`s_eDeriv` : derivative of the electric source term
 
         :param Problem prob: FDEM Problem
         :param numpy.ndarray v: vector to take product with
@@ -39,9 +39,9 @@ class BaseSrc(Survey.BaseSrc):
         :return: tuple with magnetic source term and electric source term derivatives times a vector
         """
         if v is not None:
-            return self.S_mDeriv(prob, v, adjoint), self.S_eDeriv(prob, v, adjoint)
+            return self.s_mDeriv(prob, v, adjoint), self.s_eDeriv(prob, v, adjoint)
         else:
-            return lambda v: self.S_mDeriv(prob, v, adjoint), lambda v: self.S_eDeriv(prob, v, adjoint)
+            return lambda v: self.s_mDeriv(prob, v, adjoint), lambda v: self.s_eDeriv(prob, v, adjoint)
 
     def bPrimary(self, prob):
         """
@@ -83,7 +83,7 @@ class BaseSrc(Survey.BaseSrc):
         """
         return Zero()
 
-    def S_m(self, prob):
+    def s_m(self, prob):
         """
         Magnetic source term
 
@@ -93,7 +93,7 @@ class BaseSrc(Survey.BaseSrc):
         """
         return Zero()
 
-    def S_e(self, prob):
+    def s_e(self, prob):
         """
         Electric source term
 
@@ -103,7 +103,7 @@ class BaseSrc(Survey.BaseSrc):
         """
         return Zero()
 
-    def S_mDeriv(self, prob, v, adjoint = False):
+    def s_mDeriv(self, prob, v, adjoint = False):
         """
         Derivative of magnetic source term with respect to the inversion model
 
@@ -116,7 +116,7 @@ class BaseSrc(Survey.BaseSrc):
 
         return Zero()
 
-    def S_eDeriv(self, prob, v, adjoint = False):
+    def s_eDeriv(self, prob, v, adjoint = False):
         """
         Derivative of electric source term with respect to the inversion model
 
@@ -131,22 +131,22 @@ class BaseSrc(Survey.BaseSrc):
 
 class RawVec_e(BaseSrc):
     """
-    RawVec electric source. It is defined by the user provided vector S_e
+    RawVec electric source. It is defined by the user provided vector s_e
 
     :param list rxList: receiver list
     :param float freq: frequency
-    :param numpy.array S_e: electric source term
+    :param numpy.array s_e: electric source term
     :param bool integrate: Integrate the source term (multiply by Me) [True]
     """
 
-    def __init__(self, rxList, freq, S_e, integrate=True): #, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None):
-        self._S_e = np.array(S_e, dtype=complex)
+    def __init__(self, rxList, freq, s_e, integrate=True): #, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None):
+        self._s_e = np.array(s_e, dtype=complex)
         self.freq = float(freq)
         self.integrate = integrate
 
         BaseSrc.__init__(self, rxList)
 
-    def S_e(self, prob):
+    def s_e(self, prob):
         """
         Electric source term
 
@@ -155,28 +155,28 @@ class RawVec_e(BaseSrc):
         :return: electric source term on mesh
         """
         if prob._formulation is 'EB' and self.integrate is True:
-            return prob.Me * self._S_e
-        return self._S_e
+            return prob.Me * self._s_e
+        return self._s_e
 
 
 class RawVec_m(BaseSrc):
     """
-    RawVec magnetic source. It is defined by the user provided vector S_m
+    RawVec magnetic source. It is defined by the user provided vector s_m
 
     :param float freq: frequency
     :param rxList: receiver list
-    :param numpy.array S_m: magnetic source term
+    :param numpy.array s_m: magnetic source term
     :param bool integrate: Integrate the source term (multiply by Me) [True]
     """
 
-    def __init__(self, rxList, freq, S_m, integrate=True):  #ePrimary=Zero(), bPrimary=Zero(), hPrimary=Zero(), jPrimary=Zero()):
-        self._S_m = np.array(S_m, dtype=complex)
+    def __init__(self, rxList, freq, s_m, integrate=True):  #ePrimary=Zero(), bPrimary=Zero(), hPrimary=Zero(), jPrimary=Zero()):
+        self._s_m = np.array(s_m, dtype=complex)
         self.freq = float(freq)
         self.integrate = integrate
 
         BaseSrc.__init__(self, rxList)
 
-    def S_m(self, prob):
+    def s_m(self, prob):
         """
         Magnetic source term
 
@@ -185,28 +185,28 @@ class RawVec_m(BaseSrc):
         :return: magnetic source term on mesh
         """
         if prob._formulation is 'HJ' and self.integrate is True:
-            return prob.Me * self._S_m
-        return self._S_m
+            return prob.Me * self._s_m
+        return self._s_m
 
 
 class RawVec(BaseSrc):
     """
-    RawVec source. It is defined by the user provided vectors S_m, S_e
+    RawVec source. It is defined by the user provided vectors s_m, s_e
 
     :param rxList: receiver list
     :param float freq: frequency
-    :param numpy.array S_m: magnetic source term
-    :param numpy.array S_e: electric source term
+    :param numpy.array s_m: magnetic source term
+    :param numpy.array s_e: electric source term
     :param bool integrate: Integrate the source term (multiply by Me) [True]
     """
-    def __init__(self, rxList, freq, S_m, S_e, integrate=True):
-        self._S_m = np.array(S_m, dtype=complex)
-        self._S_e = np.array(S_e, dtype=complex)
+    def __init__(self, rxList, freq, s_m, s_e, integrate=True):
+        self._s_m = np.array(s_m, dtype=complex)
+        self._s_e = np.array(s_e, dtype=complex)
         self.freq = float(freq)
         self.integrate = integrate
         BaseSrc.__init__(self, rxList)
 
-    def S_m(self, prob):
+    def s_m(self, prob):
         """
         Magnetic source term
 
@@ -215,10 +215,10 @@ class RawVec(BaseSrc):
         :return: magnetic source term on mesh
         """
         if prob._formulation is 'HJ' and self.integrate is True:
-            return prob.Me * self._S_m
-        return self._S_m
+            return prob.Me * self._s_m
+        return self._s_m
 
-    def S_e(self, prob):
+    def s_e(self, prob):
         """
         Electric source term
 
@@ -227,8 +227,8 @@ class RawVec(BaseSrc):
         :return: electric source term on mesh
         """
         if prob._formulation is 'EB' and self.integrate is True:
-            return prob.Me * self._S_e
-        return self._S_e
+            return prob.Me * self._s_e
+        return self._s_e
 
 
 class MagDipole(BaseSrc):
@@ -335,9 +335,9 @@ class MagDipole(BaseSrc):
         :return: primary magnetic field
         """
         b = self.bPrimary(prob)
-        return 1./self.mu * b 
+        return 1./self.mu * b
 
-    def S_m(self, prob):
+    def s_m(self, prob):
         """
         The magnetic source term
 
@@ -348,10 +348,10 @@ class MagDipole(BaseSrc):
 
         b_p = self.bPrimary(prob)
         if prob._formulation is 'HJ':
-            b_p = prob.Me * b_p 
+            b_p = prob.Me * b_p
         return -1j*omega(self.freq)*b_p
 
-    def S_e(self, prob):
+    def s_e(self, prob):
         """
         The electric source term
 
@@ -453,7 +453,7 @@ class MagDipole_Bfield(BaseSrc):
         b = self.bPrimary(prob)
         return 1/self.mu * b
 
-    def S_m(self, prob):
+    def s_m(self, prob):
         """
         The magnetic source term
 
@@ -466,7 +466,7 @@ class MagDipole_Bfield(BaseSrc):
             b = prob.Me * b
         return -1j*omega(self.freq)*b
 
-    def S_e(self, prob):
+    def s_e(self, prob):
         """
         The electric source term
 
@@ -565,7 +565,7 @@ class CircularLoop(BaseSrc):
         b = self.bPrimary(prob)
         return 1./self.mu*b
 
-    def S_m(self, prob):
+    def s_m(self, prob):
         """
         The magnetic source term
 
@@ -578,7 +578,7 @@ class CircularLoop(BaseSrc):
             b =  prob.Me *  b
         return -1j*omega(self.freq)*b
 
-    def S_e(self, prob):
+    def s_e(self, prob):
         """
         The electric source term
 
@@ -604,6 +604,6 @@ class CircularLoop(BaseSrc):
 
             return -C.T * (MMui_s * self.bPrimary(prob))
 
-            
+
 
 

SimPEG/EM/FDEM/SurveyFDEM.py

@@ -1,5 +1,6 @@
 import SimPEG
 from SimPEG.EM.Utils import *
+from SimPEG.EM.Base import BaseEMSurvey
 from scipy.constants import mu_0
 from SimPEG.Utils import Zero, Identity
 import SrcFDEM as Src
@@ -66,7 +67,7 @@ class Rx(SimPEG.Survey.BaseRx):
         """Grid Location projection (e.g. Ex Fy ...)"""
         return u._GLoc(self.rxType[0]) + self.knownRxTypes[self.rxType][1]
 
-    def eval(self, src, mesh, u):
+    def eval(self, src, mesh, f):
         """
         Project fields to recievers to get data.
 
@@ -79,27 +80,27 @@ class Rx(SimPEG.Survey.BaseRx):
         # projGLoc = u._GLoc(self.knownRxTypes[self.rxType][0])
         # projGLoc += self.knownRxTypes[self.rxType][1]
 
-        P = self.getP(mesh, self.projGLoc(u))
-        u_part_complex = u[src, self.projField]
+        P = self.getP(mesh, self.projGLoc(f))
+        f_part_complex = f[src, self.projField]
         # get the real or imag component
         real_or_imag = self.projComp
-        u_part = getattr(u_part_complex, real_or_imag)
-        
-        return P*u_part
+        f_part = getattr(f_part_complex, real_or_imag)
 
-    def evalDeriv(self, src, mesh, u, v, adjoint=False):
+        return P*f_part
+
+    def evalDeriv(self, src, mesh, f, v, adjoint=False):
         """
         Derivative of projected fields with respect to the inversion model times a vector.
 
         :param Source src: FDEM source
         :param Mesh mesh: mesh used
-        :param Fields u: fields object
+        :param Fields f: fields object
         :param numpy.ndarray v: vector to multiply
         :rtype: numpy.ndarray
         :return: fields projected to recievers
         """
 
-        P = self.getP(mesh, self.projGLoc(u))
+        P = self.getP(mesh, self.projGLoc(f))
 
         if not adjoint:
             Pv_complex = P * v
@@ -123,7 +124,7 @@ class Rx(SimPEG.Survey.BaseRx):
 # Survey
 ####################################################
 
-class Survey(SimPEG.Survey.BaseSurvey):
+class Survey(BaseEMSurvey):
     """
     Frequency domain electromagnetic survey
 
@@ -131,12 +132,12 @@ class Survey(SimPEG.Survey.BaseSurvey):
     """
 
     srcPair = Src.BaseSrc
-    rxPaair = Rx 
+    rxPair = Rx
 
     def __init__(self, srcList, **kwargs):
         # Sort these by frequency
         self.srcList = srcList
-        SimPEG.Survey.BaseSurvey.__init__(self, **kwargs)
+        BaseEMSurvey.__init__(self, srcList, **kwargs)
 
         _freqDict = {}
         for src in srcList:
@@ -171,24 +172,8 @@ class Survey(SimPEG.Survey.BaseSurvey):
         Returns the sources associated with a specific frequency.
         :param float freq: frequency for which we look up sources
         :rtype: dictionary
-        :return: sources at the sepcified frequency 
+        :return: sources at the sepcified frequency
         """
         assert freq in self._freqDict, "The requested frequency is not in this survey."
         return self._freqDict[freq]
 
-    def eval(self, u):
-        """
-        Project fields to receiver locations
-        :param Fields u: fields object
-        :rtype: numpy.ndarray
-        :return: data
-        """
-        data = SimPEG.Survey.Data(self)
-        for src in self.srcList:
-            for rx in src.rxList:
-                data[src, rx] = rx.eval(src, self.mesh, u)
-        return data
-
-    def evalDeriv(self, u):
-        raise Exception('Use Receivers to project fields deriv.')
-

SimPEG/EM/TDEM/BaseTDEM.py

@@ -108,11 +108,11 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
         Ainv.clean()
         return F
 
-    def Jvec(self, m, v, u=None):
+    def Jvec(self, m, v, f=None):
         """
             :param numpy.array m: Conductivity model
             :param numpy.ndarray v: vector (model object)
-            :param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m
+            :param simpegEM.TDEM.FieldsTDEM f: Fields resulting from m
             :rtype: numpy.ndarray
             :return: w (data object)
 
@@ -125,15 +125,15 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
         """
         if self.verbose: print '%s\nCalculating J(v)\n%s'%('*'*50,'*'*50)
         self.curModel = m
-        if u is None:
-            u = self.fields(m)
-        p = self.Gvec(m, v, u)
+        if f is None:
+            f = self.fields(m)
+        p = self.Gvec(m, v, f)
         y = self.solveAh(m, p)
-        Jv = self.survey.evalDeriv(u, v=y)
+        Jv = self.survey.evalDeriv(f, v=y)
         if self.verbose: print '%s\nDone calculating J(v)\n%s'%('*'*50,'*'*50)
         return - mkvc(Jv)
 
-    def Jtvec(self, m, v, u=None):
+    def Jtvec(self, m, v, f=None):
         """
             :param numpy.array m: Conductivity model
             :param numpy.ndarray,SimPEG.Survey.Data v: vector (data object)
@@ -150,15 +150,15 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
         """
         if self.verbose: print '%s\nCalculating J^T(v)\n%s'%('*'*50,'*'*50)
         self.curModel = m
-        if u is None:
-            u = self.fields(m)
+        if f is None:
+            f = self.fields(m)
 
         if not isinstance(v, self.dataPair):
             v = self.dataPair(self.survey, v)
 
-        p = self.survey.evalDeriv(u, v=v, adjoint=True)
+        p = self.survey.evalDeriv(f, v=v, adjoint=True)
         y = self.solveAht(m, p)
-        w = self.Gtvec(m, y, u)
+        w = self.Gtvec(m, y, f)
         if self.verbose: print '%s\nDone calculating J^T(v)\n%s'%('*'*50,'*'*50)
         return - mkvc(w)
 

SimPEG/Examples/DC_Forward_PseudoSection.py

@@ -2,19 +2,27 @@ from SimPEG import Mesh, Utils, np, sp
 import SimPEG.DCIP as DC
 import time
 
-def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', plotIt=True):
+def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', plotIt=True):
     """
         DC Forward Simulation
         =====================
 
-        Forward model conductive spheres in a half-space and plot a pseudo-section
-
+        Forward model two conductive spheres in a half-space and plot a 
+        pseudo-section. Assumes an infinite line source and measures along the
+        center of the spheres.
+        
+        INPUT:
+        loc     = Location of spheres [[x1,y1,z1],[x2,y2,z2]]
+        radi    = Radius of spheres [r1,r2]
+        param   = Conductivity of background and two spheres [m0,m1,m2]
+        stype   = survey type "pdp" (pole dipole) or "dpdp" (dipole dipole)
+        dtype   = Data type "appr" (app res) | "appc" (app cond) | "volt" (potential)
         Created by @fourndo on Mon Feb 01 19:28:06 2016
 
     """
 
     assert stype in ['pdp', 'dpdp'], "Source type (stype) must be pdp or dpdp (pole dipole or dipole dipole)"
-
+    assert dtype in ['appr', 'appc', 'volt'], "Data type (dtype) must be appr (app res) or appc (app cond) or volt (potential)"
 
     if loc is None:
         loc = np.c_[[-50.,0.,-50.],[50.,0.,-50.]]
@@ -27,7 +35,6 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', plotIt=True):
 
 
     # First we need to create a mesh and a model.
-
     # This is our mesh
     dx    = 5.
 
@@ -52,14 +59,10 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', plotIt=True):
     # Get index of the center
     indy = int(mesh.nCy/2)
 
-
     # Plot the model for reference
     # Define core mesh extent
     xlim = 200
-    zlim = 125
-
-    # Specify the survey type: "pdp" | "dpdp"
-
+    zlim = 100
 
     # Then specify the end points of the survey. Let's keep it simple for now and survey above the anomalies, top of the mesh
     ends = [(-175,0),(175,0)]
@@ -77,12 +80,13 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', plotIt=True):
     dl_len = np.sqrt( np.sum((locs[0,:] - locs[1,:])**2) )
     dl_x = ( Tx[-1][0,1] - Tx[0][0,0] ) / dl_len
     dl_y = ( Tx[-1][1,1] - Tx[0][1,0]  ) / dl_len
-    azm =  np.arctan(dl_y/dl_x)
+    #azm =  np.arctan(dl_y/dl_x)
 
     #Set boundary conditions
     mesh.setCellGradBC('neumann')
 
-    # Define the differential operators needed for the DC problem
+    # Define the linear system needed for the DC problem. We assume an infitite
+    # line source for simplicity.
     Div = mesh.faceDiv
     Grad = mesh.cellGrad
     Msig = Utils.sdiag(1./(mesh.aveF2CC.T*(1./model)))
@@ -145,16 +149,23 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', plotIt=True):
     print 'Forward completed'
 
     # Let's just convert the 3D format into 2D (distance along line) and plot
-    # [Tx2d, Rx2d] = DC.convertObs_DC3D_to_2D(survey, np.ones(survey.nSrc))
-    survey2D = DC.convertObs_DC3D_to_2D(survey, np.ones(survey.nSrc))
+    survey2D = DC.convertObs_DC3D_to_2D(survey, np.ones(survey.nSrc) , 'Xloc')
     survey2D.dobs =np.hstack(data)
-    # Here is an example for the first tx-rx array
+
     if plotIt:
         import matplotlib.pyplot as plt
-        fig = plt.figure()
+        fig = plt.figure(figsize=(7,7))
         ax = plt.subplot(2,1,1, aspect='equal')
-        mesh.plotSlice(np.log10(model), ax =ax, normal = 'Y', ind = indy,grid=True)
-        ax.set_title('E-W section at '+str(mesh.vectorCCy[indy])+' m')
+        # Plot the location of the spheres for reference
+        circle1=plt.Circle((loc[0,0],loc[2,0]),radi[0],color='w',fill=False, lw=3)
+        circle2=plt.Circle((loc[0,1],loc[2,1]),radi[1],color='k',fill=False, lw=3)
+        ax.add_artist(circle1)
+        ax.add_artist(circle2)
+
+        dat = mesh.plotSlice(np.log10(model), ax =ax, normal = 'Y',
+                             ind = indy,grid=True, clim = np.log10([sig.min(),sig.max()]))
+
+        ax.set_title('3-D model')
         plt.gca().set_aspect('equal', adjustable='box')
 
         plt.scatter(Tx[0][0,:],Tx[0][2,:],s=40,c='g', marker='v')
@@ -163,22 +174,34 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', plotIt=True):
         plt.ylim([-zlim,mesh.vectorNz[-1]+dx])
 
 
-        ax = plt.subplot(2,1,2, aspect='equal')
+        pos =  ax.get_position()
+        ax.set_position([pos.x0 , pos.y0 + 0.025 ,  pos.width, pos.height])
+        pos =  ax.get_position()
+        cbarax = fig.add_axes([pos.x0 , pos.y0 + 0.025 ,  pos.width, pos.height * 0.04])  ## the parameters are the specified position you set
+        cb = fig.colorbar(dat[0],cax=cbarax, orientation="horizontal",
+                          ax = ax, ticks=np.linspace(np.log10(sig.min()),
+                         np.log10(sig.max()), 3), format="$10^{%.1f}$")
+        cb.set_label("Conductivity (S/m)",size=12)
+        cb.ax.tick_params(labelsize=12)
+
+        # Second plot for the predicted apparent resistivity data
+        ax2 = plt.subplot(2,1,2, aspect='equal')
 
         # Plot the location of the spheres for reference
-        circle1=plt.Circle((loc[0,0]-Tx[0][0,0],loc[2,0]),radi[0],color='w',fill=False, lw=3)
-        circle2=plt.Circle((loc[0,1]-Tx[0][0,0],loc[2,1]),radi[1],color='k',fill=False, lw=3)
-        ax.add_artist(circle1)
-        ax.add_artist(circle2)
+        circle1=plt.Circle((loc[0,0],loc[2,0]),radi[0],color='w',fill=False, lw=3)
+        circle2=plt.Circle((loc[0,1],loc[2,1]),radi[1],color='k',fill=False, lw=3)
+        ax2.add_artist(circle1)
+        ax2.add_artist(circle2)
 
         # Add the speudo section
-        DC.plot_pseudoSection(survey2D,ax,stype)
+        dat = DC.plot_pseudoSection(survey2D,ax2,stype=stype, dtype = dtype)
 
         # plt.scatter(Tx2d[0][:],Tx[0][2,:],s=40,c='g', marker='v')
         # plt.scatter(Rx2d[0][:],Rx[0][:,2::3],s=40,c='y')
         # plt.plot(np.r_[Tx2d[0][0],Rx2d[-1][-1,-1]],np.ones(2)*mesh.vectorNz[-1], color='k')
-        plt.ylim([-zlim,mesh.vectorNz[-1]+dx])
+        ax2.set_title('Apparent Conductivity data')
 
+        plt.ylim([-zlim,mesh.vectorNz[-1]+dx])
         plt.show()
 
         return fig, ax

SimPEG/Examples/EM_FDEM_1D_Inversion.py

@@ -48,8 +48,7 @@ def run(plotIt=True):
     freqs = np.logspace(1,3,10)
     srcLoc = np.array([0., 0., 10.])
 
-    srcList = []
-    [srcList.append(EM.FDEM.Src.MagDipole([bzi],freq, srcLoc,orientation='Z')) for freq in freqs]
+    srcList = [EM.FDEM.Src.MagDipole([bzi],freq, srcLoc,orientation='Z') for freq in freqs]
 
     survey = EM.FDEM.Survey(srcList)
     prb = EM.FDEM.Problem_b(mesh, mapping=mapping)

SimPEG/Examples/EM_Schenkel_Morrison_Casing.py

@@ -0,0 +1,275 @@
+from SimPEG import *
+from SimPEG.EM import FDEM, Analytics, mu_0
+import time
+
+try:
+    from pymatsolver import MumpsSolver
+    solver = MumpsSolver
+except Exception:
+    solver = SolverLU
+    pass
+
+def run(plotIt=True):
+    """
+        EM: Schenkel and Morrison Casing Model
+        ======================================
+
+        Here we create and run a FDEM forward simulation to calculate the vertical
+        current inside a steel-cased. The model is based on the Schenkel and
+        Morrison Casing Model, and the results are used in a 2016 SEG abstract by
+        Yang et al.
+
+        - Schenkel, C.J., and H.F. Morrison, 1990, Effects of well casing on potential field measurements using downhole current sources: Geophysical prospecting, 38, 663-686.
+
+
+        The model consists of:
+        - Air: Conductivity 1e-8 S/m, above z = 0
+        - Background: conductivity 1e-2 S/m, below z = 0
+        - Casing: conductivity 1e6 S/m
+            - 300m long
+            - radius of 0.1m
+            - thickness of 6e-3m
+
+        Inside the casing, we take the same conductivity as the background.
+
+        We are using an EM code to simulate DC, so we use frequency low enough
+        that the skin depth inside the casing is longer than the casing length (f
+        = 1e-6 Hz). The plot produced is of the current inside the casing.
+
+        These results are shown in the SEG abstract by Yang et al., 2016: 3D DC
+        resistivity modeling of steel casing for reservoir monitoring using
+        equivalent resistor network. The solver used to produce these results and
+        achieve the CPU time of ~30s is Mumps, which was installed using pymatsolver_
+
+        .. _pymatsolver: https://github.com/rowanc1/pymatsolver
+
+        This example is on figshare: https://dx.doi.org/10.6084/m9.figshare.3126961.v1
+
+        If you would use this example for a code comparison, or build upon it, a
+        citation would be much appreciated!
+
+    """
+
+    if plotIt:
+        import matplotlib.pylab as plt
+
+    # ------------------ MODEL ------------------
+    sigmaair = 1e-8 # air
+    sigmaback = 1e-2 # background
+    sigmacasing = 1e6 # casing
+    sigmainside = sigmaback # inside the casing
+
+
+    casing_t = 0.006  # 1cm thickness
+    casing_l = 300   # length of the casing
+
+    casing_r = 0.1
+    casing_a = casing_r - casing_t/2. # inner radius
+    casing_b = casing_r + casing_t/2. # outer radius
+    casing_z = np.r_[-casing_l,0.]
+
+
+    # ------------------ SURVEY PARAMETERS ------------------
+    freqs = np.r_[1e-6] #[1e-1, 1, 5] # frequencies
+    dsz = -300 # down-hole z source location
+    src_loc = np.r_[0.,0.,dsz]
+    inf_loc = np.r_[0.,0.,1e4]
+
+    print 'Skin Depth: ', [(500./np.sqrt(sigmaback*_)) for _ in freqs]
+
+
+    # ------------------ MESH ------------------
+    # fine cells near well bore
+    csx1, csx2 = 2e-3, 60.
+    pfx1, pfx2 = 1.3, 1.3
+    ncx1 = np.ceil(casing_b/csx1+2)
+
+    # pad nicely to second cell size
+    npadx1 = np.floor(np.log(csx2/csx1) / np.log(pfx1))
+    hx1a,hx1b = Utils.meshTensor([(csx1,ncx1)]),Utils.meshTensor([(csx1,npadx1,pfx1)])
+    dx1 = sum(hx1a)+sum(hx1b)
+    dx1 = np.floor(dx1/csx2)
+    hx1b *= (dx1*csx2 - sum(hx1a))/sum(hx1b)
+
+    # second chunk of mesh
+    dx2 = 300. # uniform mesh out to here
+    ncx2 = np.ceil((dx2 - dx1)/csx2)
+    npadx2 = 45
+    hx2a, hx2b = Utils.meshTensor([(csx2,ncx2)]), Utils.meshTensor([(csx2,npadx2,pfx2)])
+    hx = np.hstack([hx1a,hx1b,hx2a,hx2b])
+
+    # z-direction
+    csz = 0.05
+    nza = 10
+    ncz, npadzu, npadzd = np.int(np.ceil(np.diff(casing_z)[0]/csz))+10, 68, 68 # cell size, number of core cells, number of padding cells in the x- direction
+    hz = Utils.meshTensor([(csz,npadzd,-1.3), (csz,ncz), (csz,npadzu,1.3)]) # vector of cell widths in the z-direction
+
+    # Mesh
+    mesh = Mesh.CylMesh([hx,1.,hz], [0.,0.,-np.sum(hz[:npadzu+ncz-nza])])
+
+    print 'Mesh Extent xmax: %f,: zmin: %f, zmax: %f'%(mesh.vectorCCx.max(), mesh.vectorCCz.min(), mesh.vectorCCz.max())
+    print 'Number of cells', mesh.nC
+
+    if plotIt is True:
+        fig, ax = plt.subplots(1, 1, figsize=(6, 4))
+        ax.set_title('Simulation Mesh')
+        mesh.plotGrid(ax=ax)
+        plt.show()
+
+    # Put the model on the mesh
+    sigWholespace = sigmaback*np.ones((mesh.nC))
+
+    sigBack = sigWholespace.copy()
+    sigBack[mesh.gridCC[:,2] > 0.] = sigmaair
+
+    sigCasing = sigBack.copy()
+    iCasingZ = (mesh.gridCC[:,2] <= casing_z[1]) & (mesh.gridCC[:,2] >= casing_z[0])
+    iCasingX = (mesh.gridCC[:,0] >= casing_a) & (mesh.gridCC[:,0] <= casing_b)
+    iCasing = iCasingX & iCasingZ
+    sigCasing[iCasing] = sigmacasing
+
+
+    if plotIt is True:
+
+        # plotting parameters
+        xlim = np.r_[0., 0.2]
+        zlim = np.r_[-350., 10.]
+        clim_sig = np.r_[-8,6]
+
+        # plot models
+        fig, ax = plt.subplots(1,1,figsize=(4,4))
+
+        f = plt.colorbar(mesh.plotImage(np.log10(sigCasing),ax=ax)[0], ax=ax)
+        ax.grid(which='both')
+        ax.set_title('Log_10 (Sigma)')
+        ax.set_xlim(xlim)
+        ax.set_ylim(zlim)
+        f.set_clim(clim_sig)
+
+        plt.show()
+
+
+    # -------------- Sources --------------------
+    # Define Custom Current Sources
+
+    # surface source
+    sg_x = np.zeros(mesh.vnF[0],dtype=complex)
+    sg_y = np.zeros(mesh.vnF[1],dtype=complex)
+    sg_z = np.zeros(mesh.vnF[2],dtype=complex)
+
+    nza = 2 # put the wire two cells above the surface
+    ncin = 2
+
+    # vertically directed wire
+    sgv_indx = (mesh.gridFz[:,0] > casing_a) & (mesh.gridFz[:,0] < casing_a + csx1)  # hook it up to casing at the surface
+    sgv_indz = (mesh.gridFz[:,2] <= +csz*nza) & (mesh.gridFz[:,2] >= -csz*2)
+    sgv_ind = sgv_indx & sgv_indz
+    sg_z[sgv_ind] = -1.
+
+    # horizontally directed wire
+    sgh_indx = (mesh.gridFx[:,0] > casing_a)    & (mesh.gridFx[:,0] <= inf_loc[2])
+    sgh_indz = (mesh.gridFx[:,2] > csz*(nza-0.5)) & (mesh.gridFx[:,2] < csz*(nza+0.5))
+    sgh_ind = sgh_indx & sgh_indz
+    sg_x[sgh_ind] = -1.
+
+    sgv2_indx = (mesh.gridFz[:,0] >= mesh.gridFx[sgh_ind,0].max()) & (mesh.gridFz[:,0] <= inf_loc[2]*1.2)  # hook it up to casing at the surface
+    sgv2_indz = (mesh.gridFz[:,2] <= +csz*nza) & (mesh.gridFz[:,2] >= -csz*2)
+    sgv2_ind = sgv2_indx & sgv2_indz
+    sg_z[sgv2_ind] = 1.
+
+    # assemble the source
+    sg = np.hstack([sg_x,sg_y,sg_z])
+    sg_p = [FDEM.Src.RawVec_e([],_,sg/mesh.area) for _ in freqs]
+
+    # downhole source
+    dg_x = np.zeros(mesh.vnF[0],dtype=complex)
+    dg_y = np.zeros(mesh.vnF[1],dtype=complex)
+    dg_z = np.zeros(mesh.vnF[2],dtype=complex)
+
+    # vertically directed wire
+    dgv_indx = (mesh.gridFz[:,0] < csx1)  # go through the center of the well
+    dgv_indz = (mesh.gridFz[:,2] <= +csz*nza) & (mesh.gridFz[:,2] > dsz + csz/2.)
+    dgv_ind = dgv_indx & dgv_indz
+    dg_z[dgv_ind] = -1.
+
+    # couple to the casing downhole
+    dgh_indx = mesh.gridFx[:,0] < casing_a + csx1
+    dgh_indz = (mesh.gridFx[:,2] < dsz + csz)  & (mesh.gridFx[:,2] >= dsz)
+    dgh_ind = dgh_indx & dgh_indz
+    dg_x[dgh_ind] = 1.
+
+    # horizontal part at surface
+    dgh2_indx = mesh.gridFx[:,0] <= inf_loc[2]*1.2
+    dgh2_indz = sgh_indz.copy()
+    dgh2_ind = dgh2_indx & dgh2_indz
+    dg_x[dgh2_ind] = -1.
+
+    # vertical part at surface
+    dgv2_ind = sgv2_ind.copy()
+    dg_z[dgv2_ind] = 1.
+
+    # assemble the source
+    dg = np.hstack([dg_x,dg_y,dg_z])
+    dg_p = [FDEM.Src.RawVec_e([],_,dg/mesh.area) for _ in freqs]
+
+    # ------------ Problem and Survey ---------------
+    survey = FDEM.Survey(sg_p + dg_p)
+    mapping = [('sigma', Maps.IdentityMap(mesh))]
+    problem = FDEM.Problem_h(mesh, mapping=mapping)
+    problem.pair(survey)
+
+    # ------------- Solve ---------------------------
+    t0 = time.time()
+    fieldsCasing = problem.fields(sigCasing)
+    print 'Time to solve 2 sources', time.time() - t0
+
+    # Plot current
+
+    # current density
+    jn0 = fieldsCasing[dg_p,'j']
+    jn1 = fieldsCasing[sg_p,'j']
+
+    # current
+    in0 = [mesh.area*fieldsCasing[dg_p,'j'][:,i] for i in range(len(freqs))]
+    in1 = [mesh.area*fieldsCasing[sg_p,'j'][:,i] for i in range(len(freqs))]
+
+    in0 = np.vstack(in0).T
+    in1 = np.vstack(in1).T
+
+    # integrate to get z-current inside casing
+    inds_inx = (mesh.gridFz[:,0] >= casing_a) & (mesh.gridFz[:,0] <= casing_b)
+    inds_inz = (mesh.gridFz[:,2] >= dsz )     & (mesh.gridFz[:,2] <= 0)
+    inds_fz = inds_inx & inds_inz
+
+    indsx = [False]*mesh.nFx
+    inds = list(indsx) + list(inds_fz)
+
+    in0_in = in0[np.r_[inds]]
+    in1_in = in1[np.r_[inds]]
+    z_in = mesh.gridFz[inds_fz,2]
+
+    in0_in = in0_in.reshape([in0_in.shape[0]/3,3])
+    in1_in = in1_in.reshape([in1_in.shape[0]/3,3])
+    z_in = z_in.reshape([z_in.shape[0]/3,3])
+
+    I0 = in0_in.sum(1).real
+    I1 = in1_in.sum(1).real
+    z_in = z_in[:,0]
+
+    if plotIt is True:
+        fig, ax = plt.subplots(1,2,figsize=(12,4))
+
+        ax[0].plot(z_in,np.absolute(I0), z_in,np.absolute(I1))
+        ax[0].legend(['top casing', 'bottom casing'],loc='best')
+        ax[0].set_title('Magnitude of Vertical Current in Casing')
+
+        ax[1].semilogy(z_in,np.absolute(I0), z_in,np.absolute(I1))
+        ax[1].legend(['top casing', 'bottom casing'],loc='best')
+        ax[1].set_title('Magnitude of Vertical Current in Casing')
+        ax[1].set_ylim([1e-2, 1.])
+
+        plt.show()
+
+if __name__ == '__main__':
+    run()
+

SimPEG/Examples/MT_1D_ForwardAndInversion.py

@@ -100,7 +100,7 @@ def run(plotIt=True):
     # Regularization - with a regularization mesh
     regMesh = simpeg.Mesh.TensorMesh([m1d.hx[problem.mapping.sigmaMap.maps[-1].indActive]],m1d.x0)
     reg = simpeg.Regularization.Tikhonov(regMesh)
-    reg.smoothModel = True
+    reg.mrefInSmooth = True
     reg.alpha_s = 1e-7
     reg.alpha_x = 1.
     # Inversion problem

SimPEG/Examples/__init__.py

@@ -5,6 +5,7 @@ import DC_Analytic_Dipole
 import DC_Forward_PseudoSection
 import EM_FDEM_1D_Inversion
 import EM_FDEM_Analytic_MagDipoleWholespace
+import EM_Schenkel_Morrison_Casing
 import EM_TDEM_1D_Inversion
 import FLOW_Richards_1D_Celia1990
 import Forward_BasicDirectCurrent
@@ -19,7 +20,7 @@ import Mesh_Tensor_Creation
 import MT_1D_ForwardAndInversion
 import MT_3D_Foward
 
-__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_Linear", "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"]
+__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", "Forward_BasicDirectCurrent", "Inversion_Linear", "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"]
 
 ##### AUTOIMPORTS #####
 

SimPEG/FLOW/Richards/RichardsProblem.py

@@ -45,19 +45,19 @@ class RichardsSurvey(Survey.BaseSurvey):
 
     @Utils.count
     @Utils.requires('prob')
-    def dpred(self, m, u=None):
+    def dpred(self, m, f=None):
         """
             Create the projected data from a model.
-            The field, u, (if provided) will be used for the predicted data
+            The field, f, (if provided) will be used for the predicted data
             instead of recalculating the fields (which may be expensive!).
 
             .. math::
-                d_\\text{pred} = P(u(m), m)
+                d_\\text{pred} = P(f(m), m)
 
             Where P is a projection of the fields onto the data space.
         """
-        if u is None: u = self.prob.fields(m)
-        return Utils.mkvc(self.eval(u, m))
+        if f is None: f = self.prob.fields(m)
+        return Utils.mkvc(self.eval(f, m))
 
     @Utils.requires('prob')
     def eval(self, U, m):
@@ -233,16 +233,16 @@ class RichardsProblem(Problem.BaseTimeProblem):
         return r, J
 
     @Utils.timeIt
-    def Jfull(self, m, u=None):
-        if u is None:
-            u = self.fields(m)
+    def Jfull(self, m, f=None):
+        if f is None:
+            f = self.fields(m)
 
-        nn = len(u)-1
+        nn = len(f)-1
         Asubs, Adiags, Bs = range(nn), range(nn), range(nn)
         for ii in range(nn):
             dt = self.timeSteps[ii]
-            bc = self.getBoundaryConditions(ii, u[ii])
-            Asubs[ii], Adiags[ii], Bs[ii] = self.diagsJacobian(m, u[ii], u[ii+1], dt, bc)
+            bc = self.getBoundaryConditions(ii, f[ii])
+            Asubs[ii], Adiags[ii], Bs[ii] = self.diagsJacobian(m, f[ii], f[ii+1], dt, bc)
         Ad = sp.block_diag(Adiags)
         zRight = Utils.spzeros((len(Asubs)-1)*Asubs[0].shape[0],Adiags[0].shape[1])
         zTop = Utils.spzeros(Adiags[0].shape[0], len(Adiags)*Adiags[0].shape[1])
@@ -251,7 +251,7 @@ class RichardsProblem(Problem.BaseTimeProblem):
         B = np.array(sp.vstack(Bs).todense())
 
         Ainv = self.Solver(A, **self.solverOpts)
-        P = self.survey.evalDeriv(u, m)
+        P = self.survey.evalDeriv(f, m)
         AinvB = Ainv * B
         z = np.zeros((self.mesh.nC, B.shape[1]))
         zAinvB = np.vstack((z, AinvB))
@@ -259,41 +259,41 @@ class RichardsProblem(Problem.BaseTimeProblem):
         return J
 
     @Utils.timeIt
-    def Jvec(self, m, v, u=None):
-        if u is None:
-            u = self.fields(m)
+    def Jvec(self, m, v, f=None):
+        if f is None:
+            f = self.fields(m)
 
-        JvC = range(len(u)-1) # Cell to hold each row of the long vector.
+        JvC = range(len(f)-1) # Cell to hold each row of the long vector.
 
         # This is done via forward substitution.
-        bc = self.getBoundaryConditions(0, u[0])
-        temp, Adiag, B = self.diagsJacobian(m, u[0], u[1], self.timeSteps[0], bc)
+        bc = self.getBoundaryConditions(0, f[0])
+        temp, Adiag, B = self.diagsJacobian(m, f[0], f[1], self.timeSteps[0], bc)
         Adiaginv = self.Solver(Adiag, **self.solverOpts)
         JvC[0] = Adiaginv * (B*v)
 
-        for ii in range(1,len(u)-1):
-            bc = self.getBoundaryConditions(ii, u[ii])
-            Asub, Adiag, B = self.diagsJacobian(m, u[ii], u[ii+1], self.timeSteps[ii], bc)
+        for ii in range(1,len(f)-1):
+            bc = self.getBoundaryConditions(ii, f[ii])
+            Asub, Adiag, B = self.diagsJacobian(m, f[ii], f[ii+1], self.timeSteps[ii], bc)
             Adiaginv = self.Solver(Adiag, **self.solverOpts)
             JvC[ii] = Adiaginv * (B*v - Asub*JvC[ii-1])
 
-        P = self.survey.evalDeriv(u, m)
+        P = self.survey.evalDeriv(f, m)
         return P * np.concatenate([np.zeros(self.mesh.nC)] + JvC)
 
     @Utils.timeIt
-    def Jtvec(self, m, v, u=None):
-        if u is None:
-            u = self.field(m)
+    def Jtvec(self, m, v, f=None):
+        if f is None:
+            f = self.field(m)
 
-        P = self.survey.evalDeriv(u, m)
+        P = self.survey.evalDeriv(f, m)
         PTv = P.T*v
 
         # This is done via backward substitution.
         minus = 0
         BJtv = 0
-        for ii in range(len(u)-1,0,-1):
-            bc = self.getBoundaryConditions(ii-1, u[ii-1])
-            Asub, Adiag, B = self.diagsJacobian(m, u[ii-1], u[ii], self.timeSteps[ii-1], bc)
+        for ii in range(len(f)-1,0,-1):
+            bc = self.getBoundaryConditions(ii-1, f[ii-1])
+            Asub, Adiag, B = self.diagsJacobian(m, f[ii-1], f[ii], self.timeSteps[ii-1], bc)
             #select the correct part of v
             vpart = range((ii)*Adiag.shape[0], (ii+1)*Adiag.shape[0])
             AdiaginvT = self.Solver(Adiag.T, **self.solverOpts)

SimPEG/InvProblem.py

@@ -82,23 +82,23 @@ class BaseInvProblem(object):
         self._warmstart = value
 
     def getFields(self, m, store=False, deleteWarmstart=True):
-        u = None
+        f = None
 
         for mtest, u_ofmtest in self.warmstart:
             if m is mtest:
-                u = u_ofmtest
+                f = u_ofmtest
                 if self.debug: print 'InvProb is Warm Starting!'
                 break
 
-        if u is None:
-            u = self.prob.fields(m)
+        if f is None:
+            f = self.prob.fields(m)
 
         if deleteWarmstart:
             self.warmstart = []
         if store:
-            self.warmstart += [(m,u)]
+            self.warmstart += [(m,f)]
 
-        return u
+        return f
 
     @Utils.timeIt
     def evalFunction(self, m, return_g=True, return_H=True):
@@ -109,21 +109,21 @@ class BaseInvProblem(object):
         gc.collect()
 
         # Store fields if doing a line-search
-        u = self.getFields(m, store=(return_g==False and return_H==False))
+        f = self.getFields(m, store=(return_g==False and return_H==False))
 
-        phi_d = self.dmisfit.eval(m, u=u)
+        phi_d = self.dmisfit.eval(m, f=f)
         phi_m = self.reg.eval(m)
 
-        self.dpred = self.survey.dpred(m, u=u)  # This is a cheap matrix vector calculation.
+        self.dpred = self.survey.dpred(m, f=f)  # This is a cheap matrix vector calculation.
 
         self.phi_d, self.phi_d_last  = phi_d, self.phi_d
         self.phi_m, self.phi_m_last  = phi_m, self.phi_m
 
-        f = phi_d + self.beta * phi_m
+        phi = phi_d + self.beta * phi_m
 
-        out = (f,)
+        out = (phi,)
         if return_g:
-            phi_dDeriv = self.dmisfit.evalDeriv(m, u=u)
+            phi_dDeriv = self.dmisfit.evalDeriv(m, f=f)
             phi_mDeriv = self.reg.evalDeriv(m)
 
             g = phi_dDeriv + self.beta * phi_mDeriv
@@ -131,7 +131,7 @@ class BaseInvProblem(object):
 
         if return_H:
             def H_fun(v):
-                phi_d2Deriv = self.dmisfit.eval2Deriv(m, v, u=u)
+                phi_d2Deriv = self.dmisfit.eval2Deriv(m, v, f=f)
                 phi_m2Deriv = self.reg.eval2Deriv(m, v=v)
 
                 return phi_d2Deriv + self.beta * phi_m2Deriv

SimPEG/MT/BaseMT.py

@@ -27,7 +27,7 @@ class BaseMTProblem(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, u=None):
+    def Jvec(self, m, v, f=None):
         """
         Function to calculate the data sensitivities dD/dm times a vector.
 
@@ -39,8 +39,8 @@ class BaseMTProblem(BaseFDEMProblem):
         """
 
         # Calculate the fields
-        if u is None:
-           u = self.fields(m)
+        if f is None:
+           f= self.fields(m)
         # Set current model
         self.curModel = m
         # Initiate the Jv object
@@ -56,9 +56,9 @@ class BaseMTProblem(BaseFDEMProblem):
                 # We need fDeriv_m = df/du*du/dm + df/dm
                 # Construct du/dm, it requires a solve
                 # NOTE: need to account for the 2 polarizations in the derivatives.
-                u_src = u[src,:]
+                f_src = f[src,:]
                 # dA_dm and dRHS_dm should be of size nE,2, so that we can multiply by dA_duI. The 2 columns are each of the polarizations.
-                dA_dm = self.getADeriv_m(freq, u_src, v) # Size: nE,2 (u_px,u_py) in the columns.
+                dA_dm = self.getADeriv_m(freq, f_src, v) # Size: nE,2 (u_px,u_py) in the columns.
                 dRHS_dm = self.getRHSDeriv_m(freq, v) # Size: nE,2 (u_px,u_py) in the columns.
                 if dRHS_dm is None:
                     du_dm = dA_duI * ( -dA_dm )
@@ -68,13 +68,13 @@ class BaseMTProblem(BaseFDEMProblem):
                 for rx in src.rxList:
                     # Get the projection derivative
                     # v should be of size 2*nE (for 2 polarizations)
-                    PDeriv_u = lambda t: rx.evalDeriv(src, self.mesh, u, t) # wrt u, we don't have have PDeriv wrt m
+                    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, u=None):
+    def Jtvec(self, m, v, f=None):
         """
         Function to calculate the transpose of the data sensitivities (dD/dm)^T times a vector.
 
@@ -85,8 +85,8 @@ class BaseMTProblem(BaseFDEMProblem):
             :return: Data sensitivities wrt m
         """
 
-        if u is None:
-            u = self.fields(m)
+        if f is None:
+            f = self.fields(m)
 
         self.curModel = m
 
@@ -103,15 +103,15 @@ class BaseMTProblem(BaseFDEMProblem):
 
             for src in self.survey.getSrcByFreq(freq):
                 ftype = self._fieldType + 'Solution'
-                u_src = u[src, :]
+                f_src = f[src, :]
 
                 for rx in src.rxList:
                     # Get the adjoint evalDeriv
                     # PTv needs to be nE,
-                    PTv = rx.evalDeriv(src, self.mesh, u, mkvc(v[src, rx],2), adjoint=True) # wrt u, need possibility wrt m
+                    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, u_src, mkvc(dA_duIT), adjoint=True)
+                    dA_dmT = self.getADeriv_m(freq, f_src, mkvc(dA_duIT), adjoint=True)
                     dRHS_dmT = self.getRHSDeriv_m(freq, mkvc(dA_duIT), adjoint=True)
                     # Make du_dmT
                     if dRHS_dmT is None:
@@ -129,4 +129,4 @@ class BaseMTProblem(BaseFDEMProblem):
                         raise Exception('Must be real or imag')
             # Clean the factorization, clear memory.
             ATinv.clean()
-        return Jtv
+        return Jtv

SimPEG/MT/SurveyMT.py

@@ -427,15 +427,15 @@ class Survey(SimPEGsurvey.BaseSurvey):
         assert freq in self._freqDict, "The requested frequency is not in this survey."
         return self._freqDict[freq]
 
-    def eval(self, u):
+    def eval(self, f):
         data = Data(self)
         for src in self.srcList:
             sys.stdout.flush()
             for rx in src.rxList:
-                data[src, rx] = rx.eval(src, self.mesh, u)
+                data[src, rx] = rx.eval(src, self.mesh, f)
         return data
 
-    def evalDeriv(self, u):
+    def evalDeriv(self, f):
         raise Exception('Use Transmitters to project fields deriv.')
 
 #################

SimPEG/Mesh/DiffOperators.py

@@ -307,24 +307,28 @@ class DiffOperators(object):
         return BC
     _cellGradBC_list = 'neumann'
 
+    def _cellGradStencil(self):
+        BC = self.setCellGradBC(self._cellGradBC_list)
+        n = self.vnC
+        if(self.dim == 1):
+            G = ddxCellGrad(n[0], BC[0])
+        elif(self.dim == 2):
+            G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC[0]))
+            G2 = sp.kron(ddxCellGrad(n[1], BC[1]), speye(n[0]))
+            G = sp.vstack((G1, G2), format="csr")
+        elif(self.dim == 3):
+            G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC[0]))
+            G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC[1]), speye(n[0]))
+            G3 = kron3(ddxCellGrad(n[2], BC[2]), speye(n[1]), speye(n[0]))
+            G = sp.vstack((G1, G2, G3), format="csr")
+        return G
+
     def cellGrad():
         doc = "The cell centered Gradient, takes you to cell faces."
 
         def fget(self):
             if(self._cellGrad is None):
-                BC = self.setCellGradBC(self._cellGradBC_list)
-                n = self.vnC
-                if(self.dim == 1):
-                    G = ddxCellGrad(n[0], BC[0])
-                elif(self.dim == 2):
-                    G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC[0]))
-                    G2 = sp.kron(ddxCellGrad(n[1], BC[1]), speye(n[0]))
-                    G = sp.vstack((G1, G2), format="csr")
-                elif(self.dim == 3):
-                    G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC[0]))
-                    G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC[1]), speye(n[0]))
-                    G3 = kron3(ddxCellGrad(n[2], BC[2]), speye(n[1]), speye(n[0]))
-                    G = sp.vstack((G1, G2, G3), format="csr")
+                G = self._cellGradStencil()
                 # Compute areas of cell faces & volumes
                 S = self.area
                 V = self.aveCC2F*self.vol  # Average volume between adjacent cells
@@ -361,19 +365,24 @@ class DiffOperators(object):
     _cellGradBC = None
     cellGradBC = property(**cellGradBC())
 
+    def _cellGradxStencil(self):
+        BC = ['neumann', 'neumann']
+        n = self.vnC
+        if(self.dim == 1):
+            G1 = ddxCellGrad(n[0], BC)
+        elif(self.dim == 2):
+            G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC))
+        elif(self.dim == 3):
+            G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC))
+        return G1
+
+
     def cellGradx():
         doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
 
         def fget(self):
             if getattr(self, '_cellGradx', None) is None:
-                BC = ['neumann', 'neumann']
-                n = self.vnC
-                if(self.dim == 1):
-                    G1 = ddxCellGrad(n[0], BC)
-                elif(self.dim == 2):
-                    G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC))
-                elif(self.dim == 3):
-                    G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC))
+                G1 = self._cellGradxStencil()
                 # Compute areas of cell faces & volumes
                 V = self.aveCC2F*self.vol
                 L = self.r(self.area/V, 'F','Fx', 'V')
@@ -382,17 +391,22 @@ class DiffOperators(object):
         return locals()
     cellGradx = property(**cellGradx())
 
+    def _cellGradyStencil(self):
+        if self.dim < 2: return None
+        BC = ['neumann', 'neumann']
+        n = self.vnC
+        if(self.dim == 2):
+            G2 = sp.kron(ddxCellGrad(n[1], BC), speye(n[0]))
+        elif(self.dim == 3):
+            G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC), speye(n[0]))
+        return G2
+
     def cellGrady():
         doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
         def fget(self):
             if self.dim < 2: return None
             if getattr(self, '_cellGrady', None) is None:
-                BC = ['neumann', 'neumann']
-                n = self.vnC
-                if(self.dim == 2):
-                    G2 = sp.kron(ddxCellGrad(n[1], BC), speye(n[0]))
-                elif(self.dim == 3):
-                    G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC), speye(n[0]))
+                G2 = self._cellGradyStencil()
                 # Compute areas of cell faces & volumes
                 V = self.aveCC2F*self.vol
                 L = self.r(self.area/V, 'F','Fy', 'V')
@@ -401,14 +415,19 @@ class DiffOperators(object):
         return locals()
     cellGrady = property(**cellGrady())
 
+    def _cellGradzStencil(self):
+        if self.dim < 3: return None
+        BC = ['neumann', 'neumann']
+        n = self.vnC
+        G3 = kron3(ddxCellGrad(n[2], BC), speye(n[1]), speye(n[0]))
+        return G3
+
     def cellGradz():
         doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
         def fget(self):
             if self.dim < 3: return None
             if getattr(self, '_cellGradz', None) is None:
-                BC = ['neumann', 'neumann']
-                n = self.vnC
-                G3 = kron3(ddxCellGrad(n[2], BC), speye(n[1]), speye(n[0]))
+                G3 = self._cellGradzStencil()
                 # Compute areas of cell faces & volumes
                 V = self.aveCC2F*self.vol
                 L = self.r(self.area/V, 'F','Fz', 'V')
@@ -565,7 +584,67 @@ class DiffOperators(object):
 
         return Pbc, Pin, Pout
 
+    def getBCProjWF_simple(self, discretization='CC'):
+        """
 
+        The weak form boundary condition projection matrices
+        when mixed boundary condition is used
+
+
+        """
+
+        if discretization is not 'CC':
+            raise NotImplementedError('Boundary conditions only implemented for CC discretization.')
+
+        def projBC(n):
+            ij = ([0,n], [0,1])
+            vals = [0,0]
+            vals[0] = 1
+            vals[1] = 1
+            return sp.csr_matrix((vals, ij), shape=(n+1,2))
+
+        def projDirichlet(n, bc):
+            bc = checkBC(bc)
+            ij = ([0,n], [0,1])
+            vals = [0,0]
+            if(bc[0] == 'dirichlet'):
+                vals[0] = -1
+            if(bc[1] == 'dirichlet'):
+                vals[1] = 1
+            return sp.csr_matrix((vals, ij), shape=(n+1,2))
+
+        BC = [['dirichlet','dirichlet'],['dirichlet','dirichlet'],['dirichlet','dirichlet']]
+        n = self.vnC
+        indF = self.faceBoundaryInd
+        if(self.dim == 1):
+            Pbc = projDirichlet(n[0], BC[0])
+            B = projBC(n[0])
+            indF = indF[0] | indF[1]
+            Pbc = Pbc*sdiag(self.area[indF])
+
+        elif(self.dim == 2):
+            Pbc1 = sp.kron(speye(n[1]), projDirichlet(n[0], BC[0]))
+            Pbc2 = sp.kron(projDirichlet(n[1], BC[1]), speye(n[0]))
+            Pbc = sp.block_diag((Pbc1, Pbc2), format="csr")
+            B1 = sp.kron(speye(n[1]), projBC(n[0]))
+            B2 = sp.kron(projBC(n[1]), speye(n[0]))
+            B = sp.block_diag((B1, B2), format="csr")
+            indF = np.r_[(indF[0] | indF[1]), (indF[2] | indF[3])]
+            Pbc = Pbc*sdiag(self.area[indF])
+
+        elif(self.dim == 3):
+            Pbc1 = kron3(speye(n[2]), speye(n[1]), projDirichlet(n[0], BC[0]))
+            Pbc2 = kron3(speye(n[2]), projDirichlet(n[1], BC[1]), speye(n[0]))
+            Pbc3 = kron3(projDirichlet(n[2], BC[2]), speye(n[1]), speye(n[0]))
+            Pbc = sp.block_diag((Pbc1, Pbc2, Pbc3), format="csr")
+            B1 = kron3(speye(n[2]), speye(n[1]), projBC(n[0]))
+            B2 = kron3(speye(n[2]), projBC(n[1]), speye(n[0]))
+            B3 = kron3(projBC(n[2]), speye(n[1]), speye(n[0]))
+            B = sp.block_diag((B1, B2, B3), format="csr")
+            indF = np.r_[(indF[0] | indF[1]), (indF[2] | indF[3]), (indF[4] | indF[5])]
+            Pbc = Pbc*sdiag(self.area[indF])
+
+        return Pbc, B.T
     # --------------- Averaging ---------------------
 
     @property

SimPEG/Mesh/MeshIO.py

@@ -21,7 +21,7 @@ class TensorMeshIO(object):
                 if '*' in seg:
                     st = seg
                     sp = seg.split('*')
-                    re = np.array(sp[0],dtype=int)*(' ' + sp[1])
+                    re = int(sp[0])*(' ' + sp[1])
                     line = line.replace(st,re.strip())
             return np.array(line.split(),dtype=float)
 

SimPEG/Optimization.py

@@ -888,6 +888,8 @@ class ProjectedGNCG(BFGS, Minimize, Remember):
     maxIterCG = 5
     tolCG = 1e-1
 
+    stepOffBoundsFact = 0.1 # perturbation of the inactive set off the bounds
+
     lower = -np.inf
     upper = np.inf
 
@@ -990,4 +992,19 @@ class ProjectedGNCG(BFGS, Minimize, Remember):
                     cgFlag = 1
                 # End CG Iterations
 
+            # Take a gradient step on the active cells if exist
+            if temp != self.xc.size:
+
+                rhs_a = (Active) * -self.g
+
+                dm_i = max( abs( delx ) )
+                dm_a = max( abs(rhs_a) )
+
+                # perturb inactive set off of bounds so that they are included in the step
+                delx = delx + self.stepOffBoundsFact * (rhs_a * dm_i / dm_a)
+
+            # Only keep gradients going in the right direction on the active set
+            indx = ((self.xc<=self.lower) & (delx < 0)) | ((self.xc>=self.upper) & (delx > 0))
+            delx[indx] = 0.
+
             return delx

SimPEG/Problem.py

@@ -88,28 +88,28 @@ class BaseProblem(object):
         return self.survey is not None
 
     @Utils.timeIt
-    def Jvec(self, m, v, u=None):
-        """Jvec(m, v, u=None)
+    def Jvec(self, m, v, f=None):
+        """Jvec(m, v, f=None)
 
             Effect of J(m) on a vector v.
 
             :param numpy.array m: model
             :param numpy.array v: vector to multiply
-            :param numpy.array u: fields
+            :param Fields f: fields
             :rtype: numpy.array
             :return: Jv
         """
         raise NotImplementedError('J is not yet implemented.')
 
     @Utils.timeIt
-    def Jtvec(self, m, v, u=None):
-        """Jtvec(m, v, u=None)
+    def Jtvec(self, m, v, f=None):
+        """Jtvec(m, v, f=None)
 
             Effect of transpose of J(m) on a vector v.
 
             :param numpy.array m: model
             :param numpy.array v: vector to multiply
-            :param numpy.array u: fields
+            :param Fields f: fields
             :rtype: numpy.array
             :return: JTv
         """
@@ -117,32 +117,32 @@ class BaseProblem(object):
 
 
     @Utils.timeIt
-    def Jvec_approx(self, m, v, u=None):
-        """Jvec_approx(m, v, u=None)
+    def Jvec_approx(self, m, v, f=None):
+        """Jvec_approx(m, v, f=None)
 
             Approximate effect of J(m) on a vector v
 
             :param numpy.array m: model
             :param numpy.array v: vector to multiply
-            :param numpy.array u: fields
+            :param Fields f: fields
             :rtype: numpy.array
             :return: approxJv
         """
-        return self.Jvec(m, v, u)
+        return self.Jvec(m, v, f)
 
     @Utils.timeIt
-    def Jtvec_approx(self, m, v, u=None):
-        """Jtvec_approx(m, v, u=None)
+    def Jtvec_approx(self, m, v, f=None):
+        """Jtvec_approx(m, v, f=None)
 
             Approximate effect of transpose of J(m) on a vector v.
 
             :param numpy.array m: model
             :param numpy.array v: vector to multiply
-            :param numpy.array u: fields
+            :param Fields f: fields
             :rtype: numpy.array
             :return: JTv
         """
-        return self.Jtvec(m, v, u)
+        return self.Jtvec(m, v, f)
 
     def fields(self, m):
         """
@@ -224,9 +224,9 @@ class LinearProblem(BaseProblem):
     def fields(self, m):
         return self.G.dot(m)
 
-    def Jvec(self, m, v, u=None):
+    def Jvec(self, m, v, f=None):
         return self.G.dot(v)
 
-    def Jtvec(self, m, v, u=None):
+    def Jtvec(self, m, v, f=None):
         return self.G.T.dot(v)
 

SimPEG/Regularization.py

@@ -1,5 +1,289 @@
 import Utils, Maps, Mesh, numpy as np, scipy.sparse as sp
 
+class RegularizationMesh(object):
+    """
+    **Regularization Mesh**
+
+    This contains the operators used in the regularization. Note that these
+    are not necessarily true differential operators, but are constructed from
+    a SimPEG Mesh.
+
+    :param Mesh mesh: problem mesh
+    :param numpy.array indActive: bool array, size nC, that is True where we have active cells. Used to reduce the operators so we regularize only on active cells
+    """
+
+    def __init__(self, mesh, indActive=None):
+        self.mesh = mesh
+        assert indActive is None or indActive.dtype == 'bool', 'indActive needs to be None or a bool'
+        self.indActive = indActive
+
+    @property
+    def vol(self):
+        """
+        reduced volume vector
+        :rtype: numpy.array
+        :return: reduced cell volume
+        """
+        if getattr(self, '_vol', None) is None:
+            self._vol = self._Pac.T * self.mesh.vol
+        return self._vol
+
+    @property
+    def nC(self):
+        """
+        reduced number of cells
+        :rtype: int
+        :return: number of cells being regularized
+        """
+        if getattr(self, '_nC', None) is None:
+            if self.indActive is None:
+                self._nC = self.mesh.nC
+            else:
+                self._nC = sum(self.indActive)
+        return self._nC
+
+    @property
+    def dim(self):
+        """
+        dimension of regularization mesh (1D, 2D, 3D)
+        :rtype: int
+        :return: dimension
+        """
+        if getattr(self, '_dim', None) is None:
+            self._dim = self.mesh.dim
+        return self._dim
+
+
+    @property
+    def _Pac(self):
+        """
+        projection matrix that takes from the reduced space of active cells to full modelling space (ie. nC x nindActive)
+        :rtype: scipy.sparse.csr_matrix
+        :return: active cell projection matrix
+        """
+        if getattr(self, '__Pac', None) is None:
+            if self.indActive is None:
+                self.__Pac = Utils.speye(self.mesh.nC)
+            else:
+                self.__Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
+        return self.__Pac
+
+    @property
+    def _Pafx(self):
+        """
+        projection matrix that takes from the reduced space of active x-faces to full modelling space (ie. nFx x nindActive_Fx )
+        :rtype: scipy.sparse.csr_matrix
+        :return: active face-x projection matrix
+        """
+        if getattr(self, '__Pafx', None) is None:
+            if self.indActive is None:
+                self.__Pafx = Utils.speye(self.mesh.nFx)
+            else:
+                indActive_Fx = (self.mesh.aveFx2CC.T * self.indActive) == 1
+                self.__Pafx = Utils.speye(self.mesh.nFx)[:,indActive_Fx]
+        return self.__Pafx
+
+    @property
+    def _Pafy(self):
+        """
+        projection matrix that takes from the reduced space of active y-faces to full modelling space (ie. nFy x nindActive_Fy )
+        :rtype: scipy.sparse.csr_matrix
+        :return: active face-y projection matrix
+        """
+        if getattr(self, '__Pafy', None) is None:
+            if self.indActive is None:
+                self.__Pafy = Utils.speye(self.mesh.nFy)
+            else:
+                indActive_Fy = (self.mesh.aveFy2CC.T * self.indActive) == 1
+                self.__Pafy = Utils.speye(self.mesh.nFy)[:,indActive_Fy]
+        return self.__Pafy
+
+    @property
+    def _Pafz(self):
+        """
+        projection matrix that takes from the reduced space of active z-faces to full modelling space (ie. nFz x nindActive_Fz )
+        :rtype: scipy.sparse.csr_matrix
+        :return: active face-z projection matrix
+        """
+        if getattr(self, '__Pafz', None) is None:
+            if self.indActive is None:
+                self.__Pafz = Utils.speye(self.mesh.nFz)
+            else:
+                indActive_Fz = (self.mesh.aveFz2CC.T * self.indActive) == 1
+                self.__Pafz = Utils.speye(self.mesh.nFz)[:,indActive_Fz]
+        return self.__Pafz
+
+    @property
+    def aveFx2CC(self):
+        """
+        averaging from active cell centers to active x-faces
+        :rtype: scipy.sparse.csr_matrix
+        :return: averaging from active cell centers to active x-faces
+        """
+        if getattr(self, '_aveFx2CC', None) is None:
+            self._aveFx2CC =  self._Pac.T * self.mesh.aveFx2CC * self._Pafx
+        return self._aveFx2CC
+
+    @property
+    def aveCC2Fx(self):
+        """
+        averaging from active x-faces to active cell centers
+        :rtype: scipy.sparse.csr_matrix
+        :return: averaging matrix from active x-faces to active cell centers
+        """
+        if getattr(self, '_aveCC2Fx', None) is None:
+            self._aveCC2Fx =  Utils.sdiag(1./(self.aveFx2CC.T).sum(1)) * self.aveFx2CC.T
+        return self._aveCC2Fx
+
+    @property
+    def aveFy2CC(self):
+        """
+        averaging from active cell centers to active y-faces
+        :rtype: scipy.sparse.csr_matrix
+        :return: averaging from active cell centers to active y-faces
+        """
+        if getattr(self, '_aveFy2CC', None) is None:
+            self._aveFy2CC = self._Pac.T * self.mesh.aveFy2CC * self._Pafy
+        return self._aveFy2CC
+
+    @property
+    def aveCC2Fy(self):
+        """
+        averaging from active y-faces to active cell centers
+        :rtype: scipy.sparse.csr_matrix
+        :return: averaging matrix from active y-faces to active cell centers
+        """
+        if getattr(self, '_aveCC2Fy', None) is None:
+            self._aveCC2Fy =  Utils.sdiag(1./(self.aveFy2CC.T).sum(1)) * self.aveFy2CC.T
+        return self._aveCC2Fy
+
+    @property
+    def aveFz2CC(self):
+        """
+        averaging from active cell centers to active z-faces
+        :rtype: scipy.sparse.csr_matrix
+        :return: averaging from active cell centers to active z-faces
+        """
+        if getattr(self, '_aveFz2CC', None) is None:
+            self._aveFz2CC = self._Pac.T * self.mesh.aveFz2CC * self._Pafz
+        return self._aveFz2CC
+
+    @property
+    def aveCC2Fz(self):
+        """
+        averaging from active z-faces to active cell centers
+        :rtype: scipy.sparse.csr_matrix
+        :return: averaging matrix from active z-faces to active cell centers
+        """
+        if getattr(self, '_aveCC2Fz', None) is None:
+            self._aveCC2Fz =  Utils.sdiag(1./(self.aveFz2CC.T).sum(1)) * self.aveFz2CC.T
+        return self._aveCC2Fz
+
+    @property
+    def cellDiffx(self):
+        """
+        cell centered difference in the x-direction
+        :rtype: scipy.sparse.csr_matrix
+        :return: differencing matrix for active cells in the x-direction
+        """
+        if getattr(self, '_cellDiffx', None) is None:
+            self._cellDiffx = self._Pafx.T * self.mesh.cellGradx * self._Pac
+        return self._cellDiffx
+
+    @property
+    def cellDiffy(self):
+        """
+        cell centered difference in the y-direction
+        :rtype: scipy.sparse.csr_matrix
+        :return: differencing matrix for active cells in the y-direction
+        """
+        if getattr(self, '_cellDiffy', None) is None:
+            self._cellDiffy = self._Pafy.T * self.mesh.cellGrady * self._Pac
+        return self._cellDiffy
+
+    @property
+    def cellDiffz(self):
+        """
+        cell centered difference in the z-direction
+        :rtype: scipy.sparse.csr_matrix
+        :return: differencing matrix for active cells in the z-direction
+        """
+        if getattr(self, '_cellDiffz', None) is None:
+            self._cellDiffz = self._Pafz.T * self.mesh.cellGradz * self._Pac
+        return self._cellDiffz
+
+    @property
+    def faceDiffx(self):
+        """
+        x-face differences
+        :rtype: scipy.sparse.csr_matrix
+        :return: differencing matrix for active faces in the x-direction
+        """
+        if getattr(self, '_faceDiffx', None) is None:
+            self._faceDiffx = self._Pac.T * self.mesh.faceDivx * self._Pafx
+        return self._faceDiffx
+
+    @property
+    def faceDiffy(self):
+        """
+        y-face differences
+        :rtype: scipy.sparse.csr_matrix
+        :return: differencing matrix for active faces in the y-direction
+        """
+        if getattr(self, '_faceDiffy', None) is None:
+            self._faceDiffy = self._Pac.T * self.mesh.faceDivy * self._Pafy
+        return self._faceDiffy
+
+    @property
+    def faceDiffz(self):
+        """
+        z-face differences
+        :rtype: scipy.sparse.csr_matrix
+        :return: differencing matrix for active faces in the z-direction
+        """
+        if getattr(self, '_faceDiffz', None) is None:
+            self._faceDiffz = self._Pac.T * self.mesh.faceDivz * self._Pafz
+        return self._faceDiffz
+
+    @property
+    def cellDiffxStencil(self):
+        """
+        cell centered difference stencil (no cell lengths include) in the x-direction
+        :rtype: scipy.sparse.csr_matrix
+        :return: differencing matrix for active cells in the x-direction
+        """
+        if getattr(self, '_cellDiffxStencil', None) is None:
+
+            self._cellDiffxStencil = self._Pafx.T * self.mesh._cellGradxStencil() * self._Pac
+        return self._cellDiffxStencil
+
+    @property
+    def cellDiffyStencil(self):
+        """
+        cell centered difference stencil (no cell lengths include) in the y-direction
+        :rtype: scipy.sparse.csr_matrix
+        :return: differencing matrix for active cells in the y-direction
+        """
+        if self.dim < 2: return None
+        if getattr(self, '_cellDiffyStencil', None) is None:
+
+            self._cellDiffyStencil = self._Pafy.T * self.mesh._cellGradyStencil() * self._Pac
+        return self._cellDiffyStencil
+
+    @property
+    def cellDiffzStencil(self):
+        """
+        cell centered difference stencil (no cell lengths include) in the y-direction
+        :rtype: scipy.sparse.csr_matrix
+        :return: differencing matrix for active cells in the y-direction
+        """
+        if self.dim < 3: return None
+        if getattr(self, '_cellDiffzStencil', None) is None:
+
+            self._cellDiffzStencil = self._Pafz.T * self.mesh._cellGradzStencil() * self._Pac
+        return self._cellDiffzStencil
+
+
 class BaseRegularization(object):
     """
     **Base Regularization Class**
@@ -18,12 +302,16 @@ class BaseRegularization(object):
 
     mapping = None    #: A SimPEG.Map instance.
     mesh    = None    #: A SimPEG.Mesh instance.
-    mref = None       #: Reference model.
+    mref    = None    #: Reference model.
 
     def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
         Utils.setKwargs(self, **kwargs)
-        self.mesh = mesh
         assert isinstance(mesh, Mesh.BaseMesh), "mesh must be a SimPEG.Mesh object."
+        if indActive is not None and indActive.dtype != 'bool':
+            tmp = indActive
+            indActive = np.zeros(mesh.nC, dtype=bool)
+            indActive[tmp] = True
+        self.regmesh = RegularizationMesh(mesh,indActive)
         self.mapping = mapping or self.mapPair(mesh)
         self.mapping._assertMatchesPair(self.mapPair)
         self.indActive = indActive
@@ -55,8 +343,8 @@ class BaseRegularization(object):
     @property
     def W(self):
         """Full regularization weighting matrix W."""
-        return sp.identity(self.mapping.nP)
-
+        return sp.identity(self.regmesh.nC)
+        # self.regmesh._Pac.T * sp.identity(self.regmesh.nC) * self.regmesh._Pac # or do we want sp.identity(self.mesh.nC) or even just Utils.Identity() ?
 
     @Utils.timeIt
     def eval(self, m):
@@ -112,11 +400,10 @@ class BaseRegularization(object):
 
         return mD.T * ( self.W.T * ( self.W * ( mD * v) ) )
 
-
 class Tikhonov(BaseRegularization):
     """
     """
-    smoothModel = True  #: SMOOTH and SMOOTH_MOD_DIF options
+    mrefInSmooth = True  #: SMOOTH and SMOOTH_MOD_DIF options
     alpha_s  = Utils.dependentProperty('_alpha_s', 1e-6, ['_W', '_Ws'], "Smallness weight")
     alpha_x  = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction")
     alpha_y  = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction")
@@ -126,98 +413,58 @@ class Tikhonov(BaseRegularization):
     alpha_zz = Utils.dependentProperty('_alpha_zz', 0.0, ['_W', '_Wzz'], "Weight for the second derivative in the z direction")
 
     def __init__(self, mesh, mapping=None, indActive = None, **kwargs):
-        BaseRegularization.__init__(self, mesh, mapping=mapping, **kwargs)
-        self.indActive = indActive
+        BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
 
     @property
     def Ws(self):
         """Regularization matrix Ws"""
         if getattr(self,'_Ws', None) is None:
-            self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5)
-            if self.indActive is not None:
-                Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
-                self._Ws = Pac.T * self._Ws * Pac
+            self._Ws = Utils.sdiag((self.regmesh.vol*self.alpha_s)**0.5)
         return self._Ws
 
     @property
     def Wx(self):
         """Regularization matrix Wx"""
         if getattr(self, '_Wx', None) is None:
-            Ave_x_vol = self.mesh.aveF2CC[:,:self.mesh.nFx].T*self.mesh.vol
-            self._Wx = Utils.sdiag((Ave_x_vol*self.alpha_x)**0.5)*self.mesh.cellGradx
-
-            if self.indActive is not None:
-                indActive_Fx = (self.mesh.aveFx2CC.T * self.indActive) == 1
-                Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
-                Pafx = Utils.speye(self.mesh.nFx)[:,indActive_Fx]
-                self._Wx = Pafx.T*self._Wx*Pac
-
+            Ave_x_vol = self.regmesh.aveCC2Fx * self.regmesh.vol
+            self._Wx = Utils.sdiag((Ave_x_vol*self.alpha_x)**0.5)*self.regmesh.cellDiffx
         return self._Wx
 
     @property
     def Wy(self):
         """Regularization matrix Wy"""
         if getattr(self, '_Wy', None) is None:
-            Ave_y_vol = self.mesh.aveF2CC[:,self.mesh.nFx:np.sum(self.mesh.vnF[:2])].T*self.mesh.vol
-            self._Wy = Utils.sdiag((Ave_y_vol*self.alpha_y)**0.5)*self.mesh.cellGrady
-
-            if self.indActive is not None:
-                indActive_Fy = (self.mesh.aveFy2CC.T * self.indActive) == 1
-                Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
-                Pafy = Utils.speye(self.mesh.nFy)[:,indActive_Fy]
-                self._Wy = Pafy.T*self._Wy*Pac
-
+            Ave_y_vol = self.regmesh.aveCC2Fy * self.regmesh.vol
+            self._Wy = Utils.sdiag((Ave_y_vol*self.alpha_y)**0.5)*self.regmesh.cellDiffy
         return self._Wy
 
     @property
     def Wz(self):
         """Regularization matrix Wz"""
         if getattr(self, '_Wz', None) is None:
-            Ave_z_vol = self.mesh.aveF2CC[:,np.sum(self.mesh.vnF[:2]):].T*self.mesh.vol
-            self._Wz = Utils.sdiag((Ave_z_vol*self.alpha_z)**0.5)*self.mesh.cellGradz
-
-            if self.indActive is not None:
-                indActive_Fz = (self.mesh.aveFz2CC.T * self.indActive) == 1
-                Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
-                Pafz = Utils.speye(self.mesh.nFz)[:,indActive_Fz]
-                self._Wz = Pafz.T*self._Wz*Pac
-
+            Ave_z_vol = self.regmesh.aveCC2Fz * self.regmesh.vol
+            self._Wz = Utils.sdiag((Ave_z_vol*self.alpha_z)**0.5)*self.regmesh.cellDiffz
         return self._Wz
 
     @property
     def Wxx(self):
         """Regularization matrix Wxx"""
         if getattr(self, '_Wxx', None) is None:
-            self._Wxx = Utils.sdiag((self.mesh.vol*self.alpha_xx)**0.5)*self.mesh.faceDivx*self.mesh.cellGradx
-
-            if self.indActive is not None:
-                Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
-                self._Wxx = Pac.T*self._Wxx*Pac
-
+            self._Wxx = Utils.sdiag((self.regmesh.vol*self.alpha_xx)**0.5)*self.regmesh.faceDiffx*self.regmesh.cellDiffx
         return self._Wxx
 
     @property
     def Wyy(self):
         """Regularization matrix Wyy"""
         if getattr(self, '_Wyy', None) is None:
-            self._Wyy = Utils.sdiag((self.mesh.vol*self.alpha_yy)**0.5)*self.mesh.faceDivy*self.mesh.cellGrady
-
-            if self.indActive is not None:
-                Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
-                self._Wyy = Pac.T*self._Wyy*Pac
-
+            self._Wyy = Utils.sdiag((self.regmesh.vol*self.alpha_yy)**0.5)*self.regmesh.faceDiffy*self.regmesh.cellDiffy
         return self._Wyy
 
     @property
     def Wzz(self):
         """Regularization matrix Wzz"""
         if getattr(self, '_Wzz', None) is None:
-            self._Wzz = Utils.sdiag((self.mesh.vol*self.alpha_zz)**0.5)*self.mesh.faceDivz*self.mesh.cellGradz
-
-            if self.indActive is not None:
-                Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
-                self._Wzz = Pac.T*self._Wzz*Pac
-
+            self._Wzz = Utils.sdiag((self.regmesh.vol*self.alpha_zz)**0.5)*self.regmesh.faceDiffz*self.regmesh.cellDiffz
         return self._Wzz
 
     @property
@@ -225,9 +472,9 @@ class Tikhonov(BaseRegularization):
         """Full smoothness regularization matrix W"""
         if getattr(self, '_Wsmooth', None) is None:
             wlist = (self.Wx, self.Wxx)
-            if self.mesh.dim > 1:
+            if self.regmesh.dim > 1:
                 wlist += (self.Wy, self.Wyy)
-            if self.mesh.dim > 2:
+            if self.regmesh.dim > 2:
                 wlist += (self.Wz, self.Wzz)
             self._Wsmooth = sp.vstack(wlist)
         return self._Wsmooth
@@ -242,11 +489,11 @@ class Tikhonov(BaseRegularization):
 
     @Utils.timeIt
     def eval(self, m):
-        if self.smoothModel == True:
+        if self.mrefInSmooth == True:
             r1 = self.Wsmooth * ( self.mapping * (m) )
             r2 = self.Ws * ( self.mapping * (m - self.mref) )
             return 0.5*(r1.dot(r1)+r2.dot(r2))
-        elif self.smoothModel == False:
+        elif self.mrefInSmooth == False:
             r = self.W * ( self.mapping * (m - self.mref) )
             return 0.5*r.dot(r)
 
@@ -268,7 +515,7 @@ class Tikhonov(BaseRegularization):
             R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
 
         """
-        if self.smoothModel == True:
+        if self.mrefInSmooth == True:
             mD1 = self.mapping.deriv(m)
             mD2 = self.mapping.deriv(m - self.mref)
             r1 = self.Wsmooth * ( self.mapping * (m))
@@ -276,249 +523,223 @@ class Tikhonov(BaseRegularization):
             out1 = mD1.T * ( self.Wsmooth.T * r1 )
             out2 = mD2.T * ( self.Ws.T * r2 )
             out = out1+out2
-        elif self.smoothModel == False:
+        elif self.mrefInSmooth == False:
             mD = self.mapping.deriv(m - self.mref)
             r = self.W * ( self.mapping * (m - self.mref) )
             out = mD.T * ( self.W.T * r )
         return out
 
-# <<<<<<< HEAD
 
-# class Simple(BaseRegularization):
-#     """
-#         Only for tensor mesh
-#     """
+class Simple(BaseRegularization):
+    """
+        Only for tensor mesh
+    """
 
-#     smoothModel = True  #: SMOOTH and SMOOTH_MOD_DIF options
-#     alpha_s  = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Ws'], "Smallness weight")
-#     alpha_x  = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction")
-#     alpha_y  = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction")
-#     alpha_z  = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'], "Weight for the first derivative in the z direction")
-#     alpha_xx = Utils.dependentProperty('_alpha_xx', 0.0, ['_W', '_Wxx'], "Weight for the second derivative in the x direction")
-#     alpha_yy = Utils.dependentProperty('_alpha_yy', 0.0, ['_W', '_Wyy'], "Weight for the second derivative in the y direction")
-#     alpha_zz = Utils.dependentProperty('_alpha_zz', 0.0, ['_W', '_Wzz'], "Weight for the second derivative in the z direction")
+    mrefInSmooth = True  #: SMOOTH and SMOOTH_MOD_DIF options
+    alpha_s     = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Ws'], "Smallness weight")
+    alpha_x     = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction")
+    alpha_y     = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction")
+    alpha_z     = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'], "Weight for the first derivative in the z direction")
+    wght        = 1.
+    
+    def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
+        BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
+        
+        if isinstance(self.wght,float):
+            self.wght = np.ones(self.regmesh.nC) * self.wght
 
-#     def __init__(self, mesh, mapping=None, **kwargs):
-#         BaseRegularization.__init__(self, mesh, mapping=mapping, **kwargs)
+    @property
+    def Ws(self):
+        """Regularization matrix Ws"""
+        if getattr(self,'_Ws', None) is None:
+                self._Ws = Utils.sdiag((self.regmesh.vol*self.alpha_s*self.wght)**0.5)
+        return self._Ws
+
+    @property
+    def Wx(self):
+        """Regularization matrix Wx"""
+        if getattr(self, '_Wx', None) is None:
+            self._Wx = Utils.sdiag((self.regmesh.aveCC2Fx * self.regmesh.vol*self.alpha_x*(self.regmesh.aveCC2Fx*self.wght))**0.5)*self.regmesh.cellDiffxStencil
+        return self._Wx
+
+    @property
+    def Wy(self):
+        """Regularization matrix Wy"""
+        if getattr(self, '_Wy', None) is None:
+            self._Wy = Utils.sdiag((self.regmesh.aveCC2Fy * self.regmesh.vol * self.alpha_y*(self.regmesh.aveCC2Fy*self.wght))**0.5)*self.regmesh.cellDiffyStencil
+        return self._Wy
+
+    @property
+    def Wz(self):
+        """Regularization matrix Wz"""
+        if getattr(self, '_Wz', None) is None:
+            self._Wz = Utils.sdiag((self.regmesh.aveCC2Fz * self.regmesh.vol*self.alpha_z*(self.regmesh.aveCC2Fz*self.wght))**0.5)*self.regmesh.cellDiffzStencil
+        return self._Wz
+
+    @property
+    def Wsmooth(self):
+        """Full smoothness regularization matrix W"""
+        if getattr(self, '_Wsmooth', None) is None:
+            wlist = (self.Wx,)
+            if self.regmesh.dim > 1:
+                wlist += (self.Wy,)
+            if self.regmesh.dim > 2:
+                wlist += (self.Wz,)
+            self._Wsmooth = sp.vstack(wlist)
+        return self._Wsmooth
+
+    @property
+    def W(self):
+        """Full regularization matrix W"""
+        if getattr(self, '_W', None) is None:
+            wlist = (self.Ws, self.Wsmooth)
+            self._W = sp.vstack(wlist)
+        return self._W
+
+
+    @Utils.timeIt
+    def eval(self, m):
+        if self.mrefInSmooth == True:
+            r1 = self.Wsmooth * ( self.mapping * (m) )
+            r2 = self.Ws * ( self.mapping * (m - self.mref) )
+            return 0.5*(r1.dot(r1)+r2.dot(r2))
+        elif self.mrefInSmooth == False:
+            r = self.W * ( self.mapping * (m - self.mref) )
+            return 0.5*r.dot(r)
+        return phim
 
 
 
-#     @property
-#     def Ws(self):
-#         """Regularization matrix Ws"""
-#         if getattr(self,'_Ws', None) is None:
-#                 self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5)
-#         return self._Ws
+    @Utils.timeIt
+    def evalDeriv(self, m):
+        """
 
-#     @property
-#     def Wx(self):
-#         """Regularization matrix Wx"""
-#         if getattr(self, '_Wx', None) is None:
-#             self._Wx = Utils.sdiag((self.mesh.vol*self.alpha_x)**0.5)*self.mesh.unitCellGradx
-#         return self._Wx
+        The regularization is:
 
-#     @property
-#     def Wy(self):
-#         """Regularization matrix Wy"""
-#         if getattr(self, '_Wy', None) is None:
-#             self._Wy = Utils.sdiag((self.mesh.vol*self.alpha_y)**0.5)*self.mesh.unitCellGrady
-#         return self._Wy
+        .. math::
 
-#     @property
-#     def Wz(self):
-#         """Regularization matrix Wz"""
-#         if getattr(self, '_Wz', None) is None:
-#             self._Wz = Utils.sdiag((self.mesh.vol*self.alpha_z)**0.5)*self.mesh.unitCellGradz
-#         return self._Wz
+            R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})}
 
-#     @property
-#     def Wxx(self):
-#         """Regularization matrix Wxx"""
-#         if getattr(self, '_Wxx', None) is None:
-#             self._Wxx = Utils.sdiag((self.mesh.vol*self.alpha_xx)**0.5)*self.mesh.faceDivx*self.mesh.cellGradx
-#         return self._Wxx
+        So the derivative is straight forward:
 
-#     @property
-#     def Wyy(self):
-#         """Regularization matrix Wyy"""
-#         if getattr(self, '_Wyy', None) is None:
-#             self._Wyy = Utils.sdiag((self.mesh.vol*self.alpha_yy)**0.5)*self.mesh.faceDivy*self.mesh.cellGrady
-#         return self._Wyy
+        .. math::
 
-#     @property
-#     def Wzz(self):
-#         """Regularization matrix Wzz"""
-#         if getattr(self, '_Wzz', None) is None:
-#             self._Wzz = Utils.sdiag((self.mesh.vol*self.alpha_zz)**0.5)*self.mesh.faceDivz*self.mesh.cellGradz
-#         return self._Wzz
+            R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
 
-#     @property
-#     def Wsmooth(self):
-#         """Full smoothness regularization matrix W"""
-#         if getattr(self, '_Wsmooth', None) is None:
-#             wlist = (self.Wx, self.Wxx)
-#             if self.mesh.dim > 1:
-#                 wlist += (self.Wy, self.Wyy)
-#             if self.mesh.dim > 2:
-#                 wlist += (self.Wz, self.Wzz)
-#             self._Wsmooth = sp.vstack(wlist)
-#         return self._Wsmooth
-
-#     @property
-#     def W(self):
-#         """Full regularization matrix W"""
-#         if getattr(self, '_W', None) is None:
-#             wlist = (self.Ws, self.Wsmooth)
-#             self._W = sp.vstack(wlist)
-#         return self._W
-
-#     @Utils.timeIt
-#     def eval(self, m):
-#         if self.smoothModel == True:
-#             r1 = self.Wsmooth * ( self.mapping * (m) )
-#             r2 = self.Ws * ( self.mapping * (m - self.mref) )
-#             return 0.5*(r1.dot(r1)+r2.dot(r2))
-#         elif self.smoothModel == False:
-#             r = self.W * ( self.mapping * (m - self.mref) )
-#             return 0.5*r.dot(r)
+        """
+        if self.mrefInSmooth == True:
+            mD1 = self.mapping.deriv(m)
+            mD2 = self.mapping.deriv(m - self.mref)
+            r1 = self.Wsmooth * ( self.mapping * (m))
+            r2 = self.Ws * ( self.mapping * (m - self.mref) )
+            out1 = mD1.T * ( self.Wsmooth.T * r1 )
+            out2 = mD2.T * ( self.Ws.T * r2 )
+            out = out1+out2
+        elif self.mrefInSmooth == False:
+            mD = self.mapping.deriv(m - self.mref)
+            r = self.W * ( self.mapping * (m - self.mref) )
+            out = mD.T * ( self.W.T * r )
+        return out
 
 
-#     @Utils.timeIt
-#     def evalDeriv(self, m):
-#         """
+class Sparse(Simple):
 
-#         The regularization is:
+    # set default values
+    eps      = 1e-1
+    curModel = None # use a model to compute the weights
+    gamma    = 1.
+    p        = 0.
+    qx       = 2.
+    qy       = 2.
+    qz       = 2.
+    wght     = 1.
 
-#         .. math::
+    def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
+        Simple.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
+        
+        if isinstance(self.wght,float):
+            self.wght = np.ones(self.regmesh.nC) * self.wght
 
-#             R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})}
+    @property
+    def Ws(self):
+        """Regularization matrix Ws"""
+        if getattr(self, 'curModel', None) is None:
+            self.Rs = Utils.speye(self.regmesh.nC)
 
-#         So the derivative is straight forward:
-
-#         .. math::
-
-#             R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
-
-#         """
-#         if self.smoothModel == True:
-#             mD1 = self.mapping.deriv(m)
-#             mD2 = self.mapping.deriv(m - self.mref)
-#             r1 = self.Wsmooth * ( self.mapping * (m))
-#             r2 = self.Ws * ( self.mapping * (m - self.mref) )
-#             out1 = mD1.T * ( self.Wsmooth.T * r1 )
-#             out2 = mD2.T * ( self.Ws.T * r2 )
-#             out = out1+out2
-#         elif self.smoothModel == False:
-#             mD = self.mapping.deriv(m - self.mref)
-#             r = self.W * ( self.mapping * (m - self.mref) )
-#             out = mD.T * ( self.W.T * r )
-#         return out
-
-# class SparseRegularization(Simple):
-
-#     eps = 1e-1
-
-#     m = None
-#     gamma = 1.
-#     p = 0.
-#     qx = 2.
-#     qy = 2.
-#     qz = 2.
-
-#     def __init__(self, mesh, mapping=None, **kwargs):
-#         Simple.__init__(self, mesh, mapping=mapping, **kwargs)
+        else:
+            f_m = self.curModel - self.reg.mref
+            self.rs = self.R(f_m , self.p)
+            #print "Min rs: " + str(np.max(self.rs)) + "Max rs: " + str(np.min(self.rs))
+            self.Rs = Utils.sdiag( self.rs )
+            
+        return Utils.sdiag((self.regmesh.vol*self.alpha_s*self.gamma*self.wght)**0.5)*self.Rs
 
 
-#     @property
-#     def Wsmooth(self):
-#         """Full smoothness regularization matrix W"""
-#         if getattr(self, '_Wsmooth', None) is None:
-#             wlist = (self.Wx, self.Wxx)
-#             if self.mesh.dim > 1:
-#                 wlist += (self.Wy, self.Wyy)
-#             if self.mesh.dim > 2:
-#                 wlist += (self.Wz, self.Wzz)
-#             self._Wsmooth = sp.vstack(wlist)
-#         return self._Wsmooth
+    @property
+    def Wx(self):
+        """Regularization matrix Wx"""
+            
+        if getattr(self, 'curModel', None) is None:
+            self.Rx = Utils.speye(self.regmesh.cellDiffxStencil.shape[0])
 
-#     @property
-#     def W(self):
-#         """Full regularization matrix W"""
-#         if getattr(self, '_W', None) is None:
-#             wlist = (self.Ws, self.Wsmooth)
-#             self._W = sp.vstack(wlist)
-#         return self._W
+        else:
+            f_m = self.regmesh.cellDiffxStencil * self.curModel
+            self.rx = self.R( f_m , self.qx)
+            self.Rx = Utils.sdiag( self.rx )
 
-#     @property
-#     def Ws(self):
-#         """Regularization matrix Ws"""
-#         if getattr(self, 'm', None) is None:
-#             self.Rs = Utils.speye(self.mesh.nC)
+        return Utils.sdiag(( (self.regmesh.aveCC2Fx * self.regmesh.vol) *self.alpha_x*self.gamma*(self.regmesh.aveCC2Fx*self.wght))**0.5)*self.Rx*self.regmesh.cellDiffxStencil
 
-#         else:
-#             f_m = self.m
-#             self.rs = self.R(f_m , self.p, self.eps)
-#             #print "Min rs: " + str(np.max(self.rs)) + "Max rs: " + str(np.min(self.rs))
-#             self.Rs = Utils.sdiag( self.rs )
+    @property
+    def Wy(self):
+        """Regularization matrix Wy"""
 
-#         self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s*self.gamma)**0.5)*self.Rs
+        if getattr(self, 'curModel', None) is None:
+            self.Ry = Utils.speye(self.regmesh.cellDiffyStencil.shape[0])
 
-#         return self._Ws
+        else:
+            f_m = self.regmesh.cellDiffyStencil * self.curModel
+            self.ry = self.R( f_m , self.qy)
+            self.Ry = Utils.sdiag( self.ry )
+            
+        return Utils.sdiag(((self.regmesh.aveCC2Fy * self.regmesh.vol)*self.alpha_y*self.gamma*(self.regmesh.aveCC2Fy*self.wght))**0.5)*self.Ry*self.regmesh.cellDiffyStencil
 
-#     @property
-#     def Wx(self):
-#         """Regularization matrix Wx"""
+    @property
+    def Wz(self):
+        """Regularization matrix Wz"""
 
-#         if getattr(self, 'm', None) is None:
-#             self.Rx = Utils.speye(self.mesh.unitCellGradx.shape[0])
+        if getattr(self, 'curModel', None) is None:
+            self.Rz = Utils.speye(self.regmesh.cellDiffzStencil.shape[0])
 
-#         else:
-#             f_m = self.mesh.unitCellGradx * self.m
-#             self.rx = self.R( f_m , self.qx, self.eps)
-#             self.Rx = Utils.sdiag( self.rx )
+        else:
+            f_m = self.regmesh.cellDiffzStencil * self.curModel
+            self.rz = self.R( f_m , self.qz)
+            self.Rz = Utils.sdiag( self.rz )
 
-#         if getattr(self, '_Wx', None) is None:
-#             self._Wx = Utils.sdiag((self.mesh.vol*self.alpha_x*self.gamma)**0.5)*self.Rx*self.mesh.unitCellGradx
-#         return self._Wx
+        return Utils.sdiag(((self.regmesh.aveCC2Fz * self.regmesh.vol)*self.alpha_z*self.gamma*(self.regmesh.aveCC2Fz*self.wght))**0.5)*self.Rz*self.regmesh.cellDiffzStencil
 
-#     @property
-#     def Wy(self):
-#         """Regularization matrix Wy"""
+    @property
+    def Wsmooth(self):
+        """Full smoothness regularization matrix W"""
+        #if getattr(self, '_Wsmooth', None) is None:
+        wlist = (self.Wx,)
+        if self.regmesh.dim > 1:
+            wlist += (self.Wy,)
+        if self.regmesh.dim > 2:
+            wlist += (self.Wz,)
+        #self._Wsmooth = sp.vstack(wlist)
+        return sp.vstack(wlist)
 
-#         if getattr(self, 'm', None) is None:
-#             self.Ry = Utils.speye(self.mesh.unitCellGrady.shape[0])
+    @property
+    def W(self):
+        """Full regularization matrix W"""
+        #if getattr(self, '_W', None) is None:
+        wlist = (self.Ws, self.Wsmooth)
+        #self._W = sp.vstack(wlist)
+        return sp.vstack(wlist)
+        
+    def R(self, f_m , exponent):
 
-#         else:
-#             f_m = self.mesh.unitCellGrady * self.m
-#             self.ry = self.R( f_m , self.qy, self.eps)
-#             self.Ry = Utils.sdiag( self.ry )
+        eta = (self.eps**(1-exponent/2.))**0.5
+        r = eta / (f_m**2.+self.eps**2.)**((1-exponent/2.)/2.)
 
-#         if getattr(self, '_Wy', None) is None:
-#             self._Wy = Utils.sdiag((self.mesh.vol*self.alpha_y*self.gamma)**0.5)*self.Ry*self.mesh.unitCellGrady
-#         return self._Wy
-
-#     @property
-#     def Wz(self):
-#         """Regularization matrix Wz"""
-
-#         if getattr(self, 'm', None) is None:
-#             self.Rz = Utils.speye(self.mesh.unitCellGradz.shape[0])
-
-#         else:
-#             f_m = self.mesh.unitCellGradz * self.m
-#             self.rz = self.R( f_m , self.qz, self.eps)
-#             self.Rz = Utils.sdiag( self.rz )
-
-#         if getattr(self, '_Wz', None) is None:
-#             self._Wz = Utils.sdiag((self.mesh.vol*self.alpha_z*self.gamma)**0.5)*self.Rz*self.mesh.unitCellGradz
-#         return self._Wz
-
-
-#     def R(self, f_m , p, dec):
-
-#         eta = (self.eps**(1-p/2.))**0.5
-#         r = eta / (f_m**2.+self.eps**2.)**((1-p/2.)/2.)
-
-#         return r
-# =======
-# >>>>>>> 834de582844e8e1eac95819fbe03eed55dbeb001
+        return r

SimPEG/Survey.py

@@ -295,38 +295,38 @@ class BaseSurvey(object):
 
     @Utils.count
     @Utils.requires('prob')
-    def dpred(self, m, u=None):
-        """dpred(m, u=None)
+    def dpred(self, m, f=None):
+        """dpred(m, f=None)
 
             Create the projected data from a model.
-            The field, u, (if provided) will be used for the predicted data
+            The fields, f, (if provided) will be used for the predicted data
             instead of recalculating the fields (which may be expensive!).
 
             .. math::
 
-                d_\\text{pred} = P(u(m))
+                d_\\text{pred} = P(f(m))
 
             Where P is a projection of the fields onto the data space.
         """
-        if u is None: u = self.prob.fields(m)
-        return Utils.mkvc(self.eval(u))
+        if f is None: f = self.prob.fields(m)
+        return Utils.mkvc(self.eval(f))
 
 
     @Utils.count
-    def eval(self, u):
-        """eval(u)
+    def eval(self, f):
+        """eval(f)
 
             This function projects the fields onto the data space.
 
             .. math::
 
-                d_\\text{pred} = \mathbf{P} u(m)
+                d_\\text{pred} = \mathbf{P} f(m)
         """
         raise NotImplemented('eval is not yet implemented.')
 
     @Utils.count
-    def evalDeriv(self, u):
-        """evalDeriv(u)
+    def evalDeriv(self, f):
+        """evalDeriv(f)
 
             This function s the derivative of projects the fields onto the data space.
 
@@ -337,11 +337,11 @@ class BaseSurvey(object):
         raise NotImplemented('eval is not yet implemented.')
 
     @Utils.count
-    def residual(self, m, u=None):
-        """residual(m, u=None)
+    def residual(self, m, f=None):
+        """residual(m, f=None)
 
             :param numpy.array m: geophysical model
-            :param numpy.array u: fields
+            :param numpy.array f: fields
             :rtype: numpy.array
             :return: data residual
 
@@ -352,14 +352,14 @@ class BaseSurvey(object):
                 \mu_\\text{data} = \mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}
 
         """
-        return Utils.mkvc(self.dpred(m, u=u) - self.dobs)
+        return Utils.mkvc(self.dpred(m, f=f) - self.dobs)
 
     @property
     def isSynthetic(self):
         "Check if the data is synthetic."
         return self.mtrue is not None
 
-    def makeSyntheticData(self, m, std=0.05, u=None, force=False):
+    def makeSyntheticData(self, m, std=0.05, f=None, force=False):
         """
             Make synthetic data given a model, and a standard deviation.
 
@@ -372,16 +372,16 @@ class BaseSurvey(object):
         if getattr(self, 'dobs', None) is not None and not force:
             raise Exception('Survey already has dobs. You can use force=True to override this exception.')
         self.mtrue = m
-        self.dtrue = self.dpred(m, u=u)
+        self.dtrue = self.dpred(m, f=f)
         noise = std*abs(self.dtrue)*np.random.randn(*self.dtrue.shape)
         self.dobs = self.dtrue+noise
         self.std = self.dobs*0 + std
         return self.dobs
 
 class LinearSurvey(BaseSurvey):
-    def eval(self, u):
-        return u
-    
+    def eval(self, f):
+        return f
+
     @property
     def nD(self):
         return self.prob.G.shape[0]

docs/examples/EM_Schenkel_Morrison_Casing.rst

@@ -0,0 +1,58 @@
+.. _examples_EM_Schenkel_Morrison_Casing:
+
+.. --------------------------------- ..
+..                                   ..
+..    THIS FILE IS AUTO GENEREATED   ..
+..                                   ..
+..    SimPEG/Examples/__init__.py    ..
+..                                   ..
+.. --------------------------------- ..
+
+
+EM: Schenkel and Morrison Casing Model
+======================================
+
+Here we create and run a FDEM forward simulation to calculate the vertical
+current inside a steel-cased. The model is based on the Schenkel and
+Morrison Casing Model, and the results are used in a 2016 SEG abstract by
+Yang et al.
+
+- Schenkel, C.J., and H.F. Morrison, 1990, Effects of well casing on potential field measurements using downhole current sources: Geophysical prospecting, 38, 663-686.
+
+
+The model consists of:
+- Air: Conductivity 1e-8 S/m, above z = 0
+- Background: conductivity 1e-2 S/m, below z = 0
+- Casing: conductivity 1e6 S/m
+    - 300m long
+    - radius of 0.1m
+    - thickness of 6e-3m
+
+Inside the casing, we take the same conductivity as the background.
+
+We are using an EM code to simulate DC, so we use frequency low enough
+that the skin depth inside the casing is longer than the casing length (f
+= 1e-6 Hz). The plot produced is of the current inside the casing.
+
+These results are shown in the SEG abstract by Yang et al., 2016: 3D DC
+resistivity modeling of steel casing for reservoir monitoring using
+equivalent resistor network. The solver used to produce these results and
+achieve the CPU time of ~30s is Mumps, which was installed using pymatsolver_
+
+.. _pymatsolver: https://github.com/rowanc1/pymatsolver
+
+This example is on figshare: https://dx.doi.org/10.6084/m9.figshare.3126961.v1
+
+If you would use this example for a code comparison, or build upon it, a
+citation would be much appreciated!
+
+
+
+.. plot::
+
+    from SimPEG import Examples
+    Examples.EM_Schenkel_Morrison_Casing.run()
+
+.. literalinclude:: ../../SimPEG/Examples/EM_Schenkel_Morrison_Casing.py
+    :language: python
+    :linenos:

tests/base/test_regularization.py

@@ -5,6 +5,8 @@ from scipy.sparse.linalg import dsolve
 import inspect
 
 TOL = 1e-20
+testReg = True
+testRegMesh = True
 
 class RegularizationTests(unittest.TestCase):
 
@@ -16,44 +18,82 @@ class RegularizationTests(unittest.TestCase):
         mesh3 = Mesh.TensorMesh([hx, hy, hz])
         self.meshlist = [mesh1,mesh2, mesh3]
 
-    def test_regularization(self):
-        for R in dir(Regularization):
-            r = getattr(Regularization, R)
-            if not inspect.isclass(r): continue
-            if not issubclass(r, Regularization.BaseRegularization):
-                continue
+    if testReg:
+        def test_regularization(self):
+            for R in dir(Regularization):
+                r = getattr(Regularization, R)
+                if not inspect.isclass(r): continue
+                if not issubclass(r, Regularization.BaseRegularization):
+                    continue
+
+                for i, mesh in enumerate(self.meshlist):
+
+                    print 'Testing %iD'%mesh.dim
+
+                    mapping = r.mapPair(mesh)
+                    reg = r(mesh, mapping=mapping)
+                    m = np.random.rand(mapping.nP)
+                    reg.mref = np.ones_like(m)*np.mean(m)
+
+                    print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref)
+                    passed = reg.eval(reg.mref) < TOL
+                    self.assertTrue(passed)
+
+                    print 'Check:', R
+                    passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False)
+                    self.assertTrue(passed)
+
+                    print 'Check 2 Deriv:', R
+                    passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False)
+                    self.assertTrue(passed)
+
+        def test_regularization_ActiveCells(self):
+            for R in dir(Regularization):
+                r = getattr(Regularization, R)
+                if not inspect.isclass(r): continue
+                if not issubclass(r, Regularization.BaseRegularization):
+                    continue
+
+                for i, mesh in enumerate(self.meshlist):
+
+                    print 'Testing Active Cells %iD'%(mesh.dim)
+
+                    if mesh.dim == 1:
+                        indActive = Utils.mkvc(mesh.gridCC <= 0.8)
+                    elif mesh.dim == 2:
+                        indActive = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5)
+                    elif mesh.dim == 3:
+                        indActive = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5 * 2*np.sin(2*np.pi*mesh.gridCC[:,1])+0.5)
+
+                    mapping = Maps.IdentityMap(nP=indActive.nonzero()[0].size)
+
+                    for indAct in [indActive, indActive.nonzero()[0]]: # test both bool and integers
+                        reg = r(mesh, mapping=mapping, indActive=indAct)
+                        m = np.random.rand(mesh.nC)[indAct]
+                        reg.mref = np.ones_like(m)*np.mean(m)
+
+                    print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref)
+                    passed = reg.eval(reg.mref) < TOL
+                    self.assertTrue(passed)
+
+                    print 'Check:', R
+                    passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False)
+                    self.assertTrue(passed)
+
+                    print 'Check 2 Deriv:', R
+                    passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False)
+                    self.assertTrue(passed)
+
+    if testRegMesh:
+        def test_regularizationMesh(self):
 
             for i, mesh in enumerate(self.meshlist):
 
                 print 'Testing %iD'%mesh.dim
 
-                mapping = r.mapPair(mesh)
-                reg = r(mesh, mapping=mapping)
-                m = np.random.rand(mapping.nP)
-                reg.mref = np.ones_like(m)*np.mean(m)
-
-                print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref)
-                passed = reg.eval(reg.mref) < TOL
-                self.assertTrue(passed)
-
-                print 'Check:', R
-                passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False)
-                self.assertTrue(passed)
-
-                print 'Check 2 Deriv:', R
-                passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False)
-                self.assertTrue(passed)
-
-    def test_regularization_ActiveCells(self):
-        for R in dir(Regularization):
-            r = getattr(Regularization, R)
-            if not inspect.isclass(r): continue
-            if not issubclass(r, Regularization.BaseRegularization):
-                continue
-
-            for i, mesh in enumerate(self.meshlist):
-
-                print 'Testing Active Cells %iD'%(mesh.dim)
+                # mapping = r.mapPair(mesh)
+                # reg = r(mesh, mapping=mapping)
+                # m = np.random.rand(mapping.nP)
 
                 if mesh.dim == 1:
                     indAct = Utils.mkvc(mesh.gridCC <= 0.8)
@@ -62,23 +102,9 @@ class RegularizationTests(unittest.TestCase):
                 elif mesh.dim == 3:
                     indAct = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5 * 2*np.sin(2*np.pi*mesh.gridCC[:,1])+0.5)
 
-                mapping = Maps.IdentityMap(nP=indAct.nonzero()[0].size)
+                regmesh = Regularization.RegularizationMesh(mesh, indActive=indAct)
 
-                reg = r(mesh, mapping=mapping, indActive=indAct)
-                m = np.random.rand(mesh.nC)[indAct]
-                reg.mref = np.ones_like(m)*np.mean(m)
-
-                print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref)
-                passed = reg.eval(reg.mref) < TOL
-                self.assertTrue(passed)
-
-                print 'Check:', R
-                passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False)
-                self.assertTrue(passed)
-
-                print 'Check 2 Deriv:', R
-                passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False)
-                self.assertTrue(passed)
+                assert (regmesh.vol == mesh.vol[indAct]).all()
 
 
 if __name__ == '__main__':

tests/flow/test_Richards.py

@@ -116,8 +116,8 @@ class RichardsTests1D(unittest.TestCase):
         v = np.random.rand(self.survey.nD)
         z = np.random.rand(self.M.nC)
         Hs = self.prob.fields(self.Ks)
-        vJz = v.dot(self.prob.Jvec(self.Ks,z,u=Hs))
-        zJv = z.dot(self.prob.Jtvec(self.Ks,v,u=Hs))
+        vJz = v.dot(self.prob.Jvec(self.Ks,z,f=Hs))
+        zJv = z.dot(self.prob.Jtvec(self.Ks,v,f=Hs))
         tol = TOL*(10**int(np.log10(np.abs(zJv))))
         passed = np.abs(vJz - zJv) < tol
         print 'Richards Adjoint Test - PressureHead'
@@ -188,8 +188,8 @@ class RichardsTests2D(unittest.TestCase):
         v = np.random.rand(self.survey.nD)
         z = np.random.rand(self.M.nC)
         Hs = self.prob.fields(self.Ks)
-        vJz = v.dot(self.prob.Jvec(self.Ks,z,u=Hs))
-        zJv = z.dot(self.prob.Jtvec(self.Ks,v,u=Hs))
+        vJz = v.dot(self.prob.Jvec(self.Ks,z,f=Hs))
+        zJv = z.dot(self.prob.Jtvec(self.Ks,v,f=Hs))
         tol = TOL*(10**int(np.log10(np.abs(zJv))))
         passed = np.abs(vJz - zJv) < tol
         print '2D: Richards Adjoint Test - PressureHead'
@@ -260,8 +260,8 @@ class RichardsTests3D(unittest.TestCase):
         v = np.random.rand(self.survey.nD)
         z = np.random.rand(self.M.nC)
         Hs = self.prob.fields(self.Ks)
-        vJz = v.dot(self.prob.Jvec(self.Ks,z,u=Hs))
-        zJv = z.dot(self.prob.Jtvec(self.Ks,v,u=Hs))
+        vJz = v.dot(self.prob.Jvec(self.Ks,z,f=Hs))
+        zJv = z.dot(self.prob.Jtvec(self.Ks,v,f=Hs))
         tol = TOL*(10**int(np.log10(np.abs(zJv))))
         passed = np.abs(vJz - zJv) < tol
         print '3D: Richards Adjoint Test - PressureHead'

tests/mesh/test_Mixed_boundaryPoisson.py

@@ -0,0 +1,396 @@
+import numpy as np
+import scipy.sparse as sp
+import unittest
+import matplotlib.pyplot as plt
+from SimPEG import *
+
+MESHTYPES = ['uniformTensorMesh']
+
+def getxBCyBC(mesh, alpha, beta, gamma):
+# def getxBCyBC(mesh, alpha, beta, gamma):
+    """
+
+    """
+    if mesh.dim == 1: #1D
+        if (len(alpha) != 2 or len(beta) != 2 or len(gamma) != 2):
+            raise Exception("Lenght of list, alpha should be 2")
+        fCCxm,fCCxp = mesh.cellBoundaryInd
+        nBC = fCCxm.sum()+fCCxp.sum()
+        h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp]
+
+        alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0]
+        alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1]
+
+        h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp]
+
+        a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm)
+        b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
+        a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp)
+        b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp)
+
+        xBC_xm = 0.5*a_xm
+        xBC_xp = 0.5*a_xp/b_xp
+        yBC_xm = 0.5*(1.-b_xm)
+        yBC_xp = 0.5*(1.-1./b_xp)
+
+        xBC = np.r_[xBC_xm, xBC_xp]
+        yBC = np.r_[yBC_xm, yBC_xp]
+
+    elif mesh.dim == 2: #2D
+        if (len(alpha) != 4 or len(beta) != 4 or len(gamma) != 4):
+            raise Exception("Lenght of list, alpha should be 4")
+
+        fCCxm,fCCxp,fCCym,fCCyp = mesh.cellBoundaryInd
+        fxm,fxp,fym,fyp = mesh.faceBoundaryInd
+        nBC = fCCxm.sum()+fCCxp.sum()+fCCxm.sum()+fCCxp.sum()
+        h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp]
+        h_ym, h_yp = mesh.gridCC[fCCym], mesh.gridCC[fCCyp]
+
+        alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0]
+        alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1]
+        alpha_ym, beta_ym, gamma_ym = alpha[2], beta[2], gamma[2]
+        alpha_yp, beta_yp, gamma_yp = alpha[3], beta[3], gamma[3]
+
+        h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0]
+        h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1]
+
+        a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm)
+        b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
+        a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp)
+        b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp)
+
+        a_ym = gamma_ym/(0.5*alpha_ym-beta_ym/h_ym)
+        b_ym = (0.5*alpha_ym+beta_ym/h_ym)/(0.5*alpha_ym-beta_ym/h_ym)
+        a_yp = gamma_yp/(0.5*alpha_yp-beta_yp/h_yp)
+        b_yp = (0.5*alpha_yp+beta_yp/h_yp)/(0.5*alpha_yp-beta_yp/h_yp)
+
+        xBC_xm = 0.5*a_xm
+        xBC_xp = 0.5*a_xp/b_xp
+        yBC_xm = 0.5*(1.-b_xm)
+        yBC_xp = 0.5*(1.-1./b_xp)
+        xBC_ym = 0.5*a_ym
+        xBC_yp = 0.5*a_yp/b_yp
+        yBC_ym = 0.5*(1.-b_ym)
+        yBC_yp = 0.5*(1.-1./b_yp)
+
+        sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm], np.arange(mesh.nFx)[fxp]])
+        sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym], np.arange(mesh.nFy)[fyp]])
+
+        xBC_x = np.r_[xBC_xm, xBC_xp][sortindsfx]
+        xBC_y = np.r_[xBC_ym, xBC_yp][sortindsfy]
+        yBC_x = np.r_[yBC_xm, yBC_xp][sortindsfx]
+        yBC_y = np.r_[yBC_ym, yBC_yp][sortindsfy]
+
+        xBC = np.r_[xBC_x, xBC_y]
+        yBC = np.r_[yBC_x, yBC_y]
+
+    elif mesh.dim == 3: #3D
+        if (len(alpha) != 6 or len(beta) != 6 or len(gamma) != 6):
+            raise Exception("Lenght of list, alpha should be 6")
+        fCCxm,fCCxp,fCCym,fCCyp,fCCzm,fCCzp = mesh.cellBoundaryInd
+        fxm,fxp,fym,fyp,fzm,fzp = mesh.faceBoundaryInd
+        nBC = fCCxm.sum()+fCCxp.sum()+fCCxm.sum()+fCCxp.sum()
+        h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp]
+        h_ym, h_yp = mesh.gridCC[fCCym], mesh.gridCC[fCCyp]
+        h_zm, h_zp = mesh.gridCC[fCCzm], mesh.gridCC[fCCzp]
+
+        alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0]
+        alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1]
+        alpha_ym, beta_ym, gamma_ym = alpha[2], beta[2], gamma[2]
+        alpha_yp, beta_yp, gamma_yp = alpha[3], beta[3], gamma[3]
+        alpha_zm, beta_zm, gamma_zm = alpha[2], beta[2], gamma[2]
+        alpha_zp, beta_zp, gamma_zp = alpha[3], beta[3], gamma[3]
+
+        h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0]
+        h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1]
+        h_zm, h_zp = mesh.gridCC[fCCzm,2], mesh.gridCC[fCCzp,2]
+
+        a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm)
+        b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
+        a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp)
+        b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp)
+
+        a_ym = gamma_ym/(0.5*alpha_ym-beta_ym/h_ym)
+        b_ym = (0.5*alpha_ym+beta_ym/h_ym)/(0.5*alpha_ym-beta_ym/h_ym)
+        a_yp = gamma_yp/(0.5*alpha_yp-beta_yp/h_yp)
+        b_yp = (0.5*alpha_yp+beta_yp/h_yp)/(0.5*alpha_yp-beta_yp/h_yp)
+
+        a_zm = gamma_zm/(0.5*alpha_zm-beta_zm/h_zm)
+        b_zm = (0.5*alpha_zm+beta_zm/h_zm)/(0.5*alpha_zm-beta_zm/h_zm)
+        a_zp = gamma_zp/(0.5*alpha_zp-beta_zp/h_zp)
+        b_zp = (0.5*alpha_zp+beta_zp/h_zp)/(0.5*alpha_zp-beta_zp/h_zp)
+
+        xBC_xm = 0.5*a_xm
+        xBC_xp = 0.5*a_xp/b_xp
+        yBC_xm = 0.5*(1.-b_xm)
+        yBC_xp = 0.5*(1.-1./b_xp)
+        xBC_ym = 0.5*a_ym
+        xBC_yp = 0.5*a_yp/b_yp
+        yBC_ym = 0.5*(1.-b_ym)
+        yBC_yp = 0.5*(1.-1./b_yp)
+        xBC_zm = 0.5*a_zm
+        xBC_zp = 0.5*a_zp/b_zp
+        yBC_zm = 0.5*(1.-b_zm)
+        yBC_zp = 0.5*(1.-1./b_zp)
+
+        sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm], np.arange(mesh.nFx)[fxp]])
+        sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym], np.arange(mesh.nFy)[fyp]])
+        sortindsfz = np.argsort(np.r_[np.arange(mesh.nFz)[fzm], np.arange(mesh.nFz)[fzp]])
+
+        xBC_x = np.r_[xBC_xm, xBC_xp][sortindsfx]
+        xBC_y = np.r_[xBC_ym, xBC_yp][sortindsfy]
+        xBC_z = np.r_[xBC_zm, xBC_zp][sortindsfz]
+
+        yBC_x = np.r_[yBC_xm, yBC_xp][sortindsfx]
+        yBC_y = np.r_[yBC_ym, yBC_yp][sortindsfy]
+        yBC_z = np.r_[yBC_zm, yBC_zp][sortindsfz]
+
+        xBC = np.r_[xBC_x, xBC_y, xBC_z]
+        yBC = np.r_[yBC_x, yBC_y, yBC_z]
+
+    return xBC, yBC
+
+class Test1D_InhomogeneousMixed(Tests.OrderTest):
+    name = "1D - Mixed"
+    meshTypes = MESHTYPES
+    meshDimension = 1
+    expectedOrders = 2
+    meshSizes = [4, 8, 16, 32]
+
+    def getError(self):
+        #Test function
+        phi_fun = lambda x: np.cos(np.pi*x)
+        j_fun = lambda x: np.pi*np.sin(np.pi*x)
+        phi_deriv = lambda x: -j_fun(x)
+        q_fun = lambda x: (np.pi**2)*np.cos(np.pi*x)
+
+        xc_ana = phi_fun(self.M.gridCC)
+        q_ana = q_fun(self.M.gridCC)
+        j_ana = j_fun(self.M.gridFx)
+
+        # Get boundary locations
+        vecN = self.M.vectorNx
+        vecC = self.M.vectorCCx
+
+        # Setup Mixed B.C (alpha, beta, gamma)
+        alpha_xm, alpha_xp = 1., 1.
+        beta_xm, beta_xp =  1., 1.
+        alpha = np.r_[alpha_xm, alpha_xp]
+        beta = np.r_[beta_xm, beta_xp]
+        vecN = self.M.vectorNx
+        vecC = self.M.vectorCCx
+        phi_bc = phi_fun(vecN[[0,-1]])
+        phi_deriv_bc = phi_deriv(vecN[[0,-1]])
+        gamma = alpha*phi_bc + beta*phi_deriv_bc
+        x_BC, y_BC = getxBCyBC(self.M, alpha, beta, gamma)
+
+
+        sigma = np.ones(self.M.nC)
+        Mfrho = self.M.getFaceInnerProduct(1./sigma)
+        MfrhoI = self.M.getFaceInnerProduct(1./sigma, invMat=True)
+        V = Utils.sdiag(self.M.vol)
+        Div = V*self.M.faceDiv
+        P_BC, B = self.M.getBCProjWF_simple()
+        q = q_fun(self.M.gridCC)
+        M = B*self.M.aveCC2F
+        G = Div.T - P_BC*Utils.sdiag(y_BC)*M
+        # Mrhoj = D.T V phi + P_BC*Utils.sdiag(y_BC)*M phi - P_BC*x_BC
+        rhs = V*q + Div*MfrhoI*P_BC*x_BC
+        A = Div*MfrhoI*G
+
+        if self.myTest == 'xc':
+            #TODO: fix the null space
+            Ainv = Solver(A)
+            xc = Ainv*rhs
+            err = np.linalg.norm((xc-xc_ana), np.inf)
+        else:
+            NotImplementedError
+        return err
+
+
+    def test_order(self):
+        print "==== Testing Mixed boudary conduction for CC-problem ===="
+        self.name = "1D"
+        self.myTest = 'xc'
+        self.orderTest()
+
+class Test2D_InhomogeneousMixed(Tests.OrderTest):
+    name = "2D - Mixed"
+    meshTypes = MESHTYPES
+    meshDimension = 2
+    expectedOrders = 2
+    meshSizes = [4, 8, 16, 32]
+
+    def getError(self):
+        #Test function
+        phi_fun = lambda x: np.cos(np.pi*x[:,0])*np.cos(np.pi*x[:,1])
+        j_funX = lambda x: +np.pi*np.sin(np.pi*x[:,0])*np.cos(np.pi*x[:,1])
+        j_funY = lambda x: +np.pi*np.cos(np.pi*x[:,0])*np.sin(np.pi*x[:,1])
+        phideriv_funX = lambda x: -j_funX(x)
+        phideriv_funY = lambda x: -j_funY(x)
+        q_fun = lambda x: +2*(np.pi**2)*phi_fun(x)
+
+        xc_ana = phi_fun(self.M.gridCC)
+        q_ana = q_fun(self.M.gridCC)
+        jX_ana = j_funX(self.M.gridFx)
+        jY_ana = j_funY(self.M.gridFy)
+        j_ana = np.r_[jX_ana,jY_ana]
+
+        # Get boundary locations
+        fxm,fxp,fym,fyp = self.M.faceBoundaryInd
+        gBFxm = self.M.gridFx[fxm,:]
+        gBFxp = self.M.gridFx[fxp,:]
+        gBFym = self.M.gridFy[fym,:]
+        gBFyp = self.M.gridFy[fyp,:]
+
+        # Setup Mixed B.C (alpha, beta, gamma)
+        alpha_xm, alpha_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
+        beta_xm, beta_xp =  np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
+        alpha_ym, alpha_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
+        beta_ym, beta_yp =  np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
+
+        phi_bc_xm, phi_bc_xp = phi_fun(gBFxm), phi_fun(gBFxp)
+        phi_bc_ym, phi_bc_yp = phi_fun(gBFym), phi_fun(gBFyp)
+
+        phiderivX_bc_xm, phiderivX_bc_xp = phideriv_funX(gBFxm), phideriv_funX(gBFxp)
+        phiderivY_bc_ym, phiderivY_bc_yp = phideriv_funY(gBFym), phideriv_funY(gBFyp)
+
+        gamma_fun = lambda alpha, beta, phi, phi_deriv: alpha*phi + beta*phi_deriv
+        gamma_xm = gamma_fun(alpha_xm, beta_xm, phi_bc_xm, phiderivX_bc_xm)
+        gamma_xp = gamma_fun(alpha_xp, beta_xp, phi_bc_xp, phiderivX_bc_xp)
+        gamma_ym = gamma_fun(alpha_ym, beta_ym, phi_bc_ym, phiderivY_bc_ym)
+        gamma_yp = gamma_fun(alpha_yp, beta_yp, phi_bc_yp, phiderivY_bc_yp)
+
+        alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp]
+        beta = [beta_xm, beta_xp, beta_ym, beta_yp]
+        gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp]
+
+        x_BC, y_BC = getxBCyBC(self.M, alpha, beta, gamma)
+
+
+        sigma = np.ones(self.M.nC)
+        Mfrho = self.M.getFaceInnerProduct(1./sigma)
+        MfrhoI = self.M.getFaceInnerProduct(1./sigma, invMat=True)
+        V = Utils.sdiag(self.M.vol)
+        Div = V*self.M.faceDiv
+        P_BC, B = self.M.getBCProjWF_simple()
+        q = q_fun(self.M.gridCC)
+        M = B*self.M.aveCC2F
+        G = Div.T - P_BC*Utils.sdiag(y_BC)*M
+        rhs = V*q + Div*MfrhoI*P_BC*x_BC
+        A = Div*MfrhoI*G
+
+        if self.myTest == 'xc':
+            Ainv = Solver(A)
+            xc = Ainv*rhs
+            err = np.linalg.norm((xc-xc_ana), np.inf)
+        else:
+            NotImplementedError
+        return err
+
+
+    def test_order(self):
+        print "==== Testing Mixed boudary conduction for CC-problem ===="
+        self.name = "2D"
+        self.myTest = 'xc'
+        self.orderTest()
+
+class Test3D_InhomogeneousMixed(Tests.OrderTest):
+    name = "3D - Mixed"
+    meshTypes = MESHTYPES
+    meshDimension = 3
+    expectedOrders = 2
+    meshSizes = [4, 8, 16]
+
+    def getError(self):
+        #Test function
+        phi_fun = lambda x: np.cos(np.pi*x[:,0])*np.cos(np.pi*x[:,1])*np.cos(np.pi*x[:,2])
+        j_funX = lambda x: +np.pi*np.sin(np.pi*x[:,0])*np.cos(np.pi*x[:,1])*np.cos(np.pi*x[:,2])
+        j_funY = lambda x: +np.pi*np.cos(np.pi*x[:,0])*np.sin(np.pi*x[:,1])*np.cos(np.pi*x[:,2])
+        j_funZ = lambda x: +np.pi*np.cos(np.pi*x[:,0])*np.cos(np.pi*x[:,1])*np.sin(np.pi*x[:,2])
+
+        phideriv_funX = lambda x: -j_funX(x)
+        phideriv_funY = lambda x: -j_funY(x)
+        phideriv_funZ = lambda x: -j_funZ(x)
+
+        q_fun = lambda x: 3*(np.pi**2)*phi_fun(x)
+
+        xc_ana = phi_fun(self.M.gridCC)
+        q_ana = q_fun(self.M.gridCC)
+        jX_ana = j_funX(self.M.gridFx)
+        jY_ana = j_funY(self.M.gridFy)
+        j_ana = np.r_[jX_ana,jY_ana,jY_ana]
+
+        # Get boundary locations
+        fxm,fxp,fym,fyp,fzm,fzp = self.M.faceBoundaryInd
+        gBFxm = self.M.gridFx[fxm,:]
+        gBFxp = self.M.gridFx[fxp,:]
+        gBFym = self.M.gridFy[fym,:]
+        gBFyp = self.M.gridFy[fyp,:]
+        gBFzm = self.M.gridFz[fzm,:]
+        gBFzp = self.M.gridFz[fzp,:]
+
+        # Setup Mixed B.C (alpha, beta, gamma)
+        alpha_xm, alpha_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
+        beta_xm, beta_xp =  np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
+        alpha_ym, alpha_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
+        beta_ym, beta_yp =  np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
+        alpha_zm, alpha_zp = np.ones_like(gBFzm[:,1]), np.ones_like(gBFzp[:,1])
+        beta_zm, beta_zp =  np.ones_like(gBFzm[:,1]), np.ones_like(gBFzp[:,1])
+
+
+        phi_bc_xm, phi_bc_xp = phi_fun(gBFxm), phi_fun(gBFxp)
+        phi_bc_ym, phi_bc_yp = phi_fun(gBFym), phi_fun(gBFyp)
+        phi_bc_zm, phi_bc_zp = phi_fun(gBFzm), phi_fun(gBFzp)
+
+        phiderivX_bc_xm, phiderivX_bc_xp = phideriv_funX(gBFxm), phideriv_funX(gBFxp)
+        phiderivY_bc_ym, phiderivY_bc_yp = phideriv_funY(gBFym), phideriv_funY(gBFyp)
+        phiderivY_bc_zm, phiderivY_bc_zp = phideriv_funY(gBFzm), phideriv_funY(gBFzp)
+
+        gamma_fun = lambda alpha, beta, phi, phi_deriv: alpha*phi + beta*phi_deriv
+        gamma_xm = gamma_fun(alpha_xm, beta_xm, phi_bc_xm, phiderivX_bc_xm)
+        gamma_xp = gamma_fun(alpha_xp, beta_xp, phi_bc_xp, phiderivX_bc_xp)
+        gamma_ym = gamma_fun(alpha_ym, beta_ym, phi_bc_ym, phiderivY_bc_ym)
+        gamma_yp = gamma_fun(alpha_yp, beta_yp, phi_bc_yp, phiderivY_bc_yp)
+        gamma_zm = gamma_fun(alpha_zm, beta_zm, phi_bc_zm, phiderivY_bc_zm)
+        gamma_zp = gamma_fun(alpha_zp, beta_zp, phi_bc_zp, phiderivY_bc_zp)
+
+        alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp, alpha_zm, alpha_zp]
+        beta = [beta_xm, beta_xp, beta_ym, beta_yp, beta_zm, beta_zp]
+        gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp, gamma_zm, gamma_zp]
+
+        x_BC, y_BC = getxBCyBC(self.M, alpha, beta, gamma)
+
+
+        sigma = np.ones(self.M.nC)
+        Mfrho = self.M.getFaceInnerProduct(1./sigma)
+        MfrhoI = self.M.getFaceInnerProduct(1./sigma, invMat=True)
+        V = Utils.sdiag(self.M.vol)
+        Div = V*self.M.faceDiv
+        P_BC, B = self.M.getBCProjWF_simple()
+        q = q_fun(self.M.gridCC)
+        M = B*self.M.aveCC2F
+        G = Div.T - P_BC*Utils.sdiag(y_BC)*M
+        rhs = V*q + Div*MfrhoI*P_BC*x_BC
+        A = Div*MfrhoI*G
+
+        if self.myTest == 'xc':
+            #TODO: fix the null space
+            Ainv = Solver(A)
+            xc = Ainv*rhs
+            err = np.linalg.norm((xc-xc_ana), np.inf)
+        else:
+            NotImplementedError
+        return err
+
+
+    def test_order(self):
+        print "==== Testing Mixed boudary conduction for CC-problem ===="
+        self.name = "3D"
+        self.myTest = 'xc'
+        self.orderTest()
+
+
+
+if __name__ == '__main__':
+    unittest.main()