Allow to forward G data without storing the kernel. Add figure-axes options to plot data function

Conflicts:
	.gitignore
	SimPEG/Examples/__init__.py

.travis.yml

@@ -18,6 +18,7 @@ env:
     - TEST_DIR="tests/mesh tests/base tests/utils"
     - TEST_DIR=tests/em/fdem/inverse/derivs
     - TEST_DIR=tests/em/tdem
+    - TEST_DIR=tests/pf
     - TEST_DIR=tests/dcip
     - TEST_DIR=tests/flow
     - TEST_DIR=tests/mt

SimPEG/Directives.py

@@ -300,19 +300,16 @@ class Update_IRLS(InversionDirective):
 
             # Either use the supplied epsilon, or fix base on distribution of
             # model values
-            if getattr(self, 'reg.eps', None) is None:
+            if getattr(self, 'eps', None) is None:
                 self.reg.eps_p = np.percentile(np.abs(self.invProb.curModel),self.prctile)
-            else:                 
+            else:
                 self.reg.eps_p = self.eps[0]
-                
-            if getattr(self, 'reg.eps', None) is None:
+
+            if getattr(self, 'eps', None) is None:
                 self.reg.eps_q = np.percentile(np.abs(self.reg.regmesh.cellDiffxStencil*(self.reg.mapping * self.invProb.curModel)),self.prctile)
-            else:                 
+            else:
                 self.reg.eps_q = self.eps[1]
-            
-            print "L[p qx qy qz]-norm : " + str(self.reg.norms)
-            print "eps_p: " + str(self.reg.eps_p) + " eps_q: " + str(self.reg.eps_q)
-            
+
             self.reg.norms = self.norms
             self.coolingFactor = 1.
             self.coolingRate = 1
@@ -323,6 +320,9 @@ class Update_IRLS(InversionDirective):
             self.reg.l2model = self.invProb.curModel
             self.reg.curModel = self.invProb.curModel
 
+            print "L[p qx qy qz]-norm : " + str(self.reg.norms)
+            print "eps_p: " + str(self.reg.eps_p) + " eps_q: " + str(self.reg.eps_q)
+
             if getattr(self, 'f_old', None) is None:
                 self.f_old = self.reg.eval(self.invProb.curModel)#self.invProb.evalFunction(self.invProb.curModel, return_g=False, return_H=False)
 

SimPEG/Examples/Mesh_Basic_Types.py

@@ -1,4 +1,5 @@
-from SimPEG import *
+from SimPEG import Mesh, Utils, np
+
 
 def run(plotIt=True):
     """
@@ -8,15 +9,15 @@ def run(plotIt=True):
         Here we show SimPEG used to create three different types of meshes.
 
     """
-    sz = [16,16]
+    sz = [16, 16]
     tM = Mesh.TensorMesh(sz)
     qM = Mesh.TreeMesh(sz)
     qM.refine(lambda cell: 4 if np.sqrt(((np.r_[cell.center]-0.5)**2).sum()) < 0.4 else 3)
-    rM = Mesh.CurvilinearMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate'))
+    rM = Mesh.CurvilinearMesh(Utils.meshutils.exampleLrmGrid(sz, 'rotate'))
 
     if plotIt:
         import matplotlib.pyplot as plt
-        fig, axes = plt.subplots(1,3,figsize=(14,5))
+        fig, axes = plt.subplots(1, 3, figsize=(14, 5))
         opts = {}
         tM.plotGrid(ax=axes[0], **opts)
         axes[0].set_title('TensorMesh')

SimPEG/Examples/PF_Magnetics_Analytics.py

@@ -0,0 +1,65 @@
+from SimPEG import Mesh, np, PF
+
+
+def run(plotIt=True):
+    """
+        PF: Magnetics: Analytics
+        ========================
+
+        Comparing the magnetics field in Vancouver to Seoul
+
+    """
+
+    xr = np.linspace(-300, 300, 41)
+    yr = np.linspace(-300, 300, 41)
+    X, Y = np.meshgrid(xr, yr)
+    Z = np.ones((np.size(xr), np.size(yr)))*150
+
+    # Bz component in Korea
+    inckr = -8. + 3./60
+    deckr = 54. + 9./60
+    btotkr = 50898.6
+    Bokr = PF.MagAnalytics.IDTtoxyz(inckr, deckr, btotkr)
+
+    bx, by, bz = PF.MagAnalytics.MagSphereAnaFunA(
+        X, Y, Z, 100., 0., 0., 0., 0.01, Bokr, 'secondary'
+    )
+    Bzkr = np.reshape(bz, (np.size(xr), np.size(yr)), order='F')
+
+    # Bz component in Canada
+    incca = 16. + 49./60
+    decca = 70. + 19./60
+    btotca = 54692.1
+    Boca = PF.MagAnalytics.IDTtoxyz(incca, decca, btotca)
+
+    bx, by, bz = PF.MagAnalytics.MagSphereAnaFunA(
+        X, Y, Z, 100., 0., 0., 0., 0.01, Boca, 'secondary'
+    )
+    Bzca = np.reshape(bz, (np.size(xr), np.size(yr)), order='F')
+
+    if plotIt:
+        import matplotlib.pyplot as plt
+        from mpl_toolkits.axes_grid1 import make_axes_locatable
+        fig = plt.figure(figsize=(14, 5))
+
+        ax1 = plt.subplot(121)
+        dat1 = plt.imshow(Bzkr, extent=[min(xr), max(xr), min(yr), max(yr)])
+        divider = make_axes_locatable(ax1)
+        cax1 = divider.append_axes("right", size="5%", pad=0.05)
+        ax1.set_xlabel('East-West (m)')
+        ax1.set_ylabel('South-North (m)')
+        plt.colorbar(dat1, cax=cax1)
+        ax1.set_title('$B_z$ field at Seoul, South Korea')
+
+        ax2 = plt.subplot(122)
+        dat2 = plt.imshow(Bzca, extent=[min(xr), max(xr), min(yr), max(yr)])
+        divider = make_axes_locatable(ax2)
+        cax2 = divider.append_axes("right", size="5%", pad=0.05)
+        ax2.set_xlabel('East-West (m)')
+        ax2.set_ylabel('South-North (m)')
+        plt.colorbar(dat2, cax=cax2)
+        ax2.set_title('$B_z$ field at Vancouver, Canada')
+        plt.show()
+
+if __name__ == '__main__':
+    run()

SimPEG/PF/BaseGrav.py

@@ -0,0 +1,53 @@
+from SimPEG import Maps, Survey, Utils, np, sp
+from scipy.constants import mu_0
+import re
+
+
+class LinearSurvey(Survey.BaseSurvey):
+    """Base Magnetics Survey"""
+
+    rxLoc = None #: receiver locations
+    rxType = None #: receiver type
+
+    def __init__(self, srcField, **kwargs):
+        self.srcField = srcField
+        Survey.BaseSurvey.__init__(self, **kwargs)
+
+    def eval(self, u):
+        return u
+
+    @property
+    def nD(self):
+        return self.prob.G.shape[0]
+
+    @property
+    def nRx(self):
+        return self.srcField.rxList[0].locs.shape[0]
+    # def setBackgroundField(self, SrcField):
+
+    #     if getattr(self, 'B0', None) is None:
+    #         self._B0 = SrcField.param[0] * dipazm_2_xyz(SrcField.param[1],SrcField.param[2])
+
+    #     return self._B0
+
+
+class SrcField(Survey.BaseSrc):
+    """ Define the inducing field """
+
+    param = None #: Inducing field param (Amp, Incl, Decl)
+
+    def __init__(self, rxList, **kwargs):
+        super(SrcField, self).__init__(rxList, **kwargs)
+
+
+class RxObs(Survey.BaseRx):
+    """A station location must have be located in 3-D"""
+    def __init__(self, locsXYZ, **kwargs):
+        locs = locsXYZ
+        assert locsXYZ.shape[1] == 3, 'locs must in 3-D (x,y,z).'
+        super(RxObs, self).__init__(locs, 'tmi', storeProjections=False, **kwargs)
+
+    @property
+    def nD(self):
+        """Number of data in the receiver."""
+        return self.locs[0].shape[0]

SimPEG/PF/BaseMag.py

@@ -0,0 +1,194 @@
+from SimPEG import Maps, Survey, Utils, np, sp
+from scipy.constants import mu_0
+import re
+
+
+class BaseMagSurvey(Survey.BaseSurvey):
+    """Base Magnetics Survey"""
+
+    rxLoc = None   #: receiver locations
+    rxType = None   #: receiver type
+
+    def __init__(self, **kwargs):
+        Survey.BaseSurvey.__init__(self, **kwargs)
+
+    def setBackgroundField(self, Inc, Dec, Btot):
+
+        Bx = Btot*np.cos(Inc/180.*np.pi)*np.sin(Dec/180.*np.pi)
+        By = Btot*np.cos(Inc/180.*np.pi)*np.cos(Dec/180.*np.pi)
+        Bz = -Btot*np.sin(Inc/180.*np.pi)
+
+        self.B0 = np.r_[Bx, By, Bz]
+
+    @property
+    def Qfx(self):
+        if getattr(self, '_Qfx', None) is None:
+            self._Qfx = self.prob.mesh.getInterpolationMat(self.rxLoc, 'Fx')
+        return self._Qfx
+
+    @property
+    def Qfy(self):
+        if getattr(self, '_Qfy', None) is None:
+            self._Qfy = self.prob.mesh.getInterpolationMat(self.rxLoc, 'Fy')
+        return self._Qfy
+
+    @property
+    def Qfz(self):
+        if getattr(self, '_Qfz', None) is None:
+            self._Qfz = self.prob.mesh.getInterpolationMat(self.rxLoc, 'Fz')
+        return self._Qfz
+
+    def projectFields(self, u):
+        """
+            This function projects the fields onto the data space.
+
+            Especially, here for we use total magnetic intensity (TMI) data,
+            which is common in practice.
+
+            First we project our B on to data location
+
+            .. math::
+
+                \mathbf{B}_{rec} = \mathbf{P} \mathbf{B}
+
+            then we take the dot product between B and b_0
+
+            .. math ::
+
+                \\text{TMI} = \\vec{B}_s \cdot \hat{B}_0
+
+        """
+        # TODO: There can be some different tyes of data like |B| or B
+
+        bfx = self.Qfx*u['B']
+        bfy = self.Qfy*u['B']
+        bfz = self.Qfz*u['B']
+
+        # Generate unit vector
+        B0 = self.prob.survey.B0
+        Bot = np.sqrt(B0[0]**2+B0[1]**2+B0[2]**2)
+        box = B0[0]/Bot
+        boy = B0[1]/Bot
+        boz = B0[2]/Bot
+
+        # return bfx*box + bfx*boy + bfx*boz
+        return bfx*box + bfy*boy + bfz*boz
+
+    @Utils.count
+    def projectFieldsDeriv(self, B):
+        """
+            This function projects the fields onto the data space.
+
+            .. math::
+
+                \\frac{\partial d_\\text{pred}}{\partial \mathbf{B}} = \mathbf{P}
+
+            Especially, this function is for TMI data type
+
+        """
+        # Generate unit vector
+        B0 = self.prob.survey.B0
+        Bot = np.sqrt(B0[0]**2+B0[1]**2+B0[2]**2)
+        box = B0[0]/Bot
+        boy = B0[1]/Bot
+        boz = B0[2]/Bot
+
+        return self.Qfx*box+self.Qfy*boy+self.Qfz*boz
+
+    def projectFieldsAsVector(self, B):
+
+        bfx = self.Qfx*B
+        bfy = self.Qfy*B
+        bfz = self.Qfz*B
+
+        return np.r_[bfx, bfy, bfz]
+
+
+class LinearSurvey(Survey.BaseSurvey):
+    """Base Magnetics Survey"""
+
+    rxLoc = None #: receiver locations
+    rxType = None #: receiver type
+
+    def __init__(self, srcField, **kwargs):
+        self.srcField = srcField
+        Survey.BaseSurvey.__init__(self, **kwargs)
+
+    def eval(self, u):
+        return u
+
+    @property
+    def nD(self):
+        return self.prob.G.shape[0]
+
+    @property
+    def nRx(self):
+        return self.srcField.rxList[0].locs.shape[0]
+    # def setBackgroundField(self, SrcField):
+
+    #     if getattr(self, 'B0', None) is None:
+    #         self._B0 = SrcField.param[0] * dipazm_2_xyz(SrcField.param[1],SrcField.param[2])
+
+    #     return self._B0
+
+
+class SrcField(Survey.BaseSrc):
+    """ Define the inducing field """
+
+    param = None #: Inducing field param (Amp, Incl, Decl)
+
+    def __init__(self, rxList, **kwargs):
+        super(SrcField, self).__init__(rxList, **kwargs)
+
+
+class RxObs(Survey.BaseRx):
+    """A station location must have be located in 3-D"""
+    def __init__(self, locsXYZ, **kwargs):
+        locs = locsXYZ
+        assert locsXYZ.shape[1] == 3, 'locs must in 3-D (x,y,z).'
+        super(RxObs, self).__init__(locs, 'tmi', storeProjections=False, **kwargs)
+
+    @property
+    def nD(self):
+        """Number of data in the receiver."""
+        return self.locs[0].shape[0]
+
+
+class MagSurveyBx(object):
+    """docstring for MagSurveyBx"""
+    def __init__(self, **kwargs):
+        Survey.BaseData.__init__(self, **kwargs)
+
+    def projectFields(self, B):
+        bfx = self.Qfx*B
+        return bfx
+
+
+class BaseMagMap(Maps.IdentityMap):
+    """BaseMagMap"""
+
+    def __init__(self, mesh, **kwargs):
+        Maps.IdentityMap.__init__(self, mesh)
+
+    def _transform(self, m):
+
+        return mu_0*(1 + m)
+
+    def deriv(self, m):
+
+        return mu_0*sp.identity(self.nP)
+
+
+class WeightMap(Maps.IdentityMap):
+    """Weighted Map for distributed parameters"""
+
+    def __init__(self, nP, weight, **kwargs):
+        Maps.IdentityMap.__init__(self, nP)
+        self.mesh = None
+        self.weight = weight
+
+    def _transform(self, m):
+        return m*self.weight
+
+    def deriv(self, m):
+        return Utils.sdiag(self.weight)

SimPEG/PF/Gravity.py

@@ -0,0 +1,576 @@
+from SimPEG import *
+import BaseGrav as GRAV
+import re
+
+
+class GravityIntegral(Problem.BaseProblem):
+
+    # surveyPair = Survey.LinearSurvey
+    forwardOnly = False #: Determine if the forward matrix is stored (defaut:yes)
+    actInd = None #: Active cell indices provided
+    rtype       = 'z'
+
+    def __init__(self, mesh, mapping=None, **kwargs):
+        Problem.BaseProblem.__init__(self, mesh, mapping=mapping, **kwargs)
+
+    def fwr_op(self):
+        # Add forward function
+        # kappa = self.curModel.kappa TODO
+        rho = self.mapping*self.curModel
+
+        if self.forwardOnly:
+
+            if getattr(self, 'actInd', None) is not None:
+
+                if self.actInd.dtype=='bool':
+                    inds = np.asarray([inds for inds, elem in enumerate(self.actInd, 1) if elem], dtype = int) - 1
+                else:
+                    inds = self.actInd
+
+            else:
+
+                inds = np.asarray(range(self.mesh.nC))
+
+            nC = len(inds)
+
+            # Create active cell projector
+            P = sp.csr_matrix(
+                (np.ones(nC), (inds, range(nC))),
+                shape=(self.mesh.nC, nC)
+            )
+
+            # Create vectors of nodal location (lower and upper corners for each cell)
+            xn = self.mesh.vectorNx
+            yn = self.mesh.vectorNy
+            zn = self.mesh.vectorNz
+
+            yn2, xn2, zn2 = np.meshgrid(yn[1:], xn[1:], zn[1:])
+            yn1, xn1, zn1 = np.meshgrid(yn[0:-1], xn[0:-1], zn[0:-1])
+
+            Yn = P.T*np.c_[mkvc(yn1), mkvc(yn2)]
+            Xn = P.T*np.c_[mkvc(xn1), mkvc(xn2)]
+            Zn = P.T*np.c_[mkvc(zn1), mkvc(zn2)]
+
+            rxLoc = self.survey.srcField.rxList[0].locs
+            ndata = rxLoc.shape[0]
+
+
+            # Pre-allocate space and create magnetization matrix if required
+            # Pre-allocate space
+            if self.rtype == 'z':
+
+                fwr_d = np.zeros(self.survey.nRx)
+
+            elif self.rtype == 'xyz':
+
+                fwr_d = np.zeros(3*self.survey.nRx)
+
+            else:
+
+                print """Flag must be either 'z' | 'xyz', please revised"""
+                return
+
+
+            # Add counter to dsiplay progress. Good for large problems
+            count = -1;
+            for ii in range(ndata):
+
+
+                tx, ty, tz = get_T_mat(Xn, Yn, Zn, rxLoc[ii, :])
+
+
+                if self.rtype =='z':
+                    fwr_d[ii] =tz.dot(rho)
+
+                elif self.rtype =='xyz':
+                    fwr_d[ii] = tx.dot(rho)
+                    fwr_d[ii+ndata] = ty.dot(rho)
+                    fwr_d[ii+2*ndata] = tz.dot(rho)
+
+
+            # Display progress
+                count = progress(ii,count,ndata)
+
+            print "Done 100% ...forward operator completed!!\n"
+
+            return fwr_d
+
+        else:
+            return self.G.dot(rho)
+
+    def fields(self, m):
+        self.curModel = m
+
+        fields = self.fwr_op()
+
+        return fields
+
+        # return self.G.dot(self.mapping*(m))
+
+    def Jvec(self, m, v, f=None):
+        dmudm = self.mapping.deriv(m)
+        return self.G.dot(dmudm*v)
+
+    def Jtvec(self, m, v, f=None):
+        dmudm = self.mapping.deriv(m)
+        return dmudm.T * (self.G.T.dot(v))
+
+    @property
+    def G(self):
+        if not self.ispaired:
+            raise Exception('Need to pair!')
+
+        if getattr(self, '_G', None) is None:
+            self._G = self.Intrgl_Fwr_Op( 'z' )
+
+        return self._G
+
+    def Intrgl_Fwr_Op(self, flag):
+
+        """
+
+        Gravity forward operator in integral form
+
+        flag        = 'z' | 'xyz'
+
+        Return
+        _G        = Linear forward modeling operation
+
+        Created on March, 15th 2016
+
+        @author: dominiquef
+
+         """
+        # Find non-zero cells
+        # inds = np.nonzero(actv)[0]
+        if getattr(self, 'actInd', None) is not None:
+
+            if self.actInd.dtype=='bool':
+                inds = np.asarray([inds for inds, elem in enumerate(self.actInd, 1) if elem], dtype = int) - 1
+            else:
+                inds = self.actInd
+
+        else:
+
+            inds = np.asarray(range(self.mesh.nC))
+
+        nC = len(inds)
+
+        # Create active cell projector
+        P = sp.csr_matrix(
+            (np.ones(nC), (inds, range(nC))),
+            shape=(self.mesh.nC, nC)
+        )
+
+        # Create vectors of nodal location (lower and upper corners for each cell)
+        xn = self.mesh.vectorNx
+        yn = self.mesh.vectorNy
+        zn = self.mesh.vectorNz
+
+        yn2, xn2, zn2 = np.meshgrid(yn[1:], xn[1:], zn[1:])
+        yn1, xn1, zn1 = np.meshgrid(yn[0:-1], xn[0:-1], zn[0:-1])
+
+        Yn = P.T*np.c_[mkvc(yn1), mkvc(yn2)]
+        Xn = P.T*np.c_[mkvc(xn1), mkvc(xn2)]
+        Zn = P.T*np.c_[mkvc(zn1), mkvc(zn2)]
+
+        rxLoc = self.survey.srcField.rxList[0].locs
+        ndata = rxLoc.shape[0]
+
+        # Pre-allocate space and create magnetization matrix if required
+        # Pre-allocate space
+        if flag == 'z':
+
+            G = np.zeros((ndata, nC))
+
+        elif flag == 'xyz':
+
+            G = np.zeros((int(3*ndata), nC))
+
+        else:
+
+            print """Flag must be either 'z' | 'xyz', please revised"""
+            return
+
+
+        # Loop through all observations and create forward operator (ndata-by-nC)
+        print "Begin calculation of forward operator: " + flag
+
+        # Add counter to dsiplay progress. Good for large problems
+        count = -1;
+        for ii in range(ndata):
+
+            if flag=='z':
+                tt = get_T_mat(Xn, Yn, Zn, rxLoc[ii, :])
+                G[ii, :] = tt
+
+            elif flag == 'xyz':
+                print "Sorry 3-component not implemented yet"
+
+            # Display progress
+            count = progress(ii, count, ndata)
+
+        print "Done 100% ...forward operator completed!!\n"
+
+        return G
+
+
+def get_T_mat(Xn, Yn, Zn, rxLoc):
+    """
+    Load in the active nodes of a tensor mesh and computes the gravity tensor
+    for a given observation location rxLoc[obsx, obsy, obsz]
+
+    INPUT:
+    Xn, Yn, Zn: Node location matrix for the lower and upper most corners of
+                all cells in the mesh shape[nC,2]
+    M
+    OUTPUT:
+    Tx = [Txx Txy Txz]
+    Ty = [Tyx Tyy Tyz]
+    Tz = [Tzx Tzy Tzz]
+
+    where each elements have dimension 1-by-nC.
+    Only the upper half 5 elements have to be computed since symetric.
+    Currently done as for-loops but will eventually be changed to vector
+    indexing, once the topography has been figured out.
+
+    """
+    NewtG=6.6738e-3
+    eps = 1e-10 # add a small value to the locations to avoid /0
+
+    nC = Xn.shape[0]
+
+    # Pre-allocate space for 1D array
+    tx = np.zeros((1,nC))
+    ty = np.zeros((1,nC))
+    tz = np.zeros((1,nC))
+
+    dz = rxLoc[2] - Zn + eps
+
+    dy = Yn - rxLoc[1] + eps
+
+    dx = Xn - rxLoc[0] + eps
+
+    # Compute contribution from each corners
+    for aa in range(2):
+        for bb in range(2):
+            for cc in range(2):
+
+                r = (
+                        dx[:, aa] ** 2 +
+                        dy[:, bb] ** 2 +
+                        dz[:, cc] ** 2
+                    ) ** (0.50)
+
+                tx = tx - NewtG * (-1) ** aa * (-1) ** bb * (-1) ** cc * (
+                    dy[:, bb] * np.log(dz[:, cc] + r) +
+                    dz[:, cc] * np.log(dy[:, bb] + r) -
+                    dx[:, aa] * np.arctan(dy[:, bb] * dz[:, cc] / (dx[:, aa] * r)))
+
+                ty = ty - NewtG * (-1) ** aa * (-1) ** bb * (-1) ** cc * (
+                    dx[:, aa] * np.log(dz[:, cc] + r) +
+                    dz[:, cc] * np.log(dx[:, aa] + r) -
+                    dy[:, bb] * np.arctan(dx[:, aa] * dz[:, cc] / (dy[:, bb] * r)))
+
+                tz = tz - NewtG * (-1) ** aa * (-1) ** bb * (-1) ** cc * (
+                    dx[:, aa] * np.log(dy[:, bb] + r) +
+                    dy[:, bb] * np.log(dx[:, aa] + r) -
+                    dz[:, cc] * np.arctan(dx[:, aa] * dy[:, bb] / (dz[:, cc] * r)))
+
+    return tx,ty,tz
+
+
+def progress(iter, prog, final):
+    """
+    progress(iter,prog,final)
+
+    Function measuring the progress of a process and print to screen the %.
+    Useful to estimate the remaining runtime of a large problem.
+
+    Created on Dec, 20th 2015
+
+    @author: dominiquef
+    """
+    arg = np.floor(float(iter)/float(final)*10.)
+
+    if  arg > prog:
+
+        strg = "Done " + str(arg*10) + " %"
+        print strg
+        prog = arg
+
+    return prog
+
+
+def writeUBCobs(filename, survey, d):
+    """
+    writeUBCobs(filename,survey,d)
+
+    Function writing an observation file in UBC-GRAV3D format.
+
+    INPUT
+    filename    : Name of out file including directory
+    survey
+    flag          : dobs | dpred
+
+    OUTPUT
+    Obsfile
+
+    """
+
+    rxLoc = survey.srcField.rxList[0].locs
+
+    wd = survey.std
+
+    data = np.c_[rxLoc, d, wd]
+
+    with file(filename, 'w') as fid:
+        fid.write('%i\n' % len(d))
+        np.savetxt(fid, data, fmt='%e', delimiter=' ', newline='\n')
+
+    print "Observation file saved to: " + filename
+
+
+def getActiveTopo(mesh, topo, flag):
+    """
+    getActiveTopo(mesh,topo)
+
+    Function creates an active cell model from topography
+
+    INPUT
+    mesh        : Mesh in SimPEG format
+    topo        : Scatter points defining topography [x,y,z]
+
+    OUTPUT
+    actv        : Active cell model
+
+    """
+    import scipy.interpolate as interpolation
+
+    if flag == 'N':
+        Zn = np.zeros((mesh.nNx, mesh.nNy))
+        # wght    = np.zeros((mesh.nNx,mesh.nNy))
+        cx = mesh.vectorNx
+        cy = mesh.vectorNy
+
+    F = interpolation.NearestNDInterpolator(topo[:, 0:2], topo[:, 2])
+    [Y, X] = np.meshgrid(cy, cx)
+
+    Zn = F(X, Y)
+
+    actv = np.zeros((mesh.nCx, mesh.nCy, mesh.nCz))
+
+    if flag == 'N':
+        Nz = mesh.vectorNz[1:]
+
+    for jj in range(mesh.nCy):
+
+        for ii in range(mesh.nCx):
+
+            temp = [kk for kk in range(len(Nz)) if np.all(Zn[ii:(ii+2), jj:(jj+2)] > Nz[kk]) ]
+            actv[ii, jj, temp] = 1
+
+    actv = mkvc(actv == 1)
+
+    inds = np.asarray([inds for inds, elem in enumerate(actv, 1) if elem], dtype = int) - 1
+
+    return inds
+
+def plot_obs_2D(survey,varstr, fig = None):
+    """ Function plot_obs(rxLoc,d,wd)
+    Generate a 2d interpolated plot from scatter points of data
+
+    INPUT
+    rxLoc       : Observation locations [x,y,z]
+    d           : Data vector
+    wd          : Uncertainty vector
+
+    OUTPUT
+    figure()
+
+    Created on Dec, 27th 2015
+
+    @author: dominiquef
+
+    """
+
+    from scipy.interpolate import griddata
+    import pylab as plt
+
+    rxLoc   = survey.srcField.rxList[0].locs
+    d       = survey.dobs
+    wd      = survey.std
+
+    # Create grid of points
+    x = np.linspace(rxLoc[:,0].min(), rxLoc[:,0].max(), 100)
+    y = np.linspace(rxLoc[:,1].min(), rxLoc[:,1].max(), 100)
+
+    X, Y = np.meshgrid(x,y)
+
+    # Interpolate
+    d_grid = griddata(rxLoc[:,0:2],d,(X,Y), method ='linear')
+
+    # Plot result
+    if fig is None:
+        fig = plt.figure()
+
+    ax = plt.subplot()
+    plt.imshow(d_grid, extent=[x.min(), x.max(), y.min(), y.max()],origin = 'lower', cmap='plasma')
+    plt.colorbar(fraction=0.02)
+    plt.contour(X,Y, d_grid,10)
+    plt.scatter(rxLoc[:,0],rxLoc[:,1], c=d, s=20)
+    plt.title(varstr)
+    plt.gca().set_aspect('equal', adjustable='box')
+
+def readUBCgravObs(obs_file):
+
+    """
+    Read UBC grav file format
+
+    INPUT:
+    :param fileName, path to the UBC obs grav file
+
+    OUTPUT:
+    :param survey
+
+    """
+
+    fid = open(obs_file,'r')
+
+    # First line has the number of rows
+    line = fid.readline()
+    ndat = np.array(line.split(),dtype=int)
+
+    # Pre-allocate space for obsx, obsy, obsz, data, uncert
+    line = fid.readline()
+    temp = np.array(line.split(),dtype=float)
+
+    d  = np.zeros(ndat, dtype=float)
+    wd = np.zeros(ndat, dtype=float)
+    locXYZ = np.zeros( (ndat,3), dtype=float)
+
+    for ii in range(ndat):
+
+        temp = np.array(line.split(),dtype=float)
+        locXYZ[ii,:] = temp[:3]
+        d[ii] = temp[3]
+        wd[ii] = temp[4]
+        line = fid.readline()
+
+    rxLoc = GRAV.RxObs(locXYZ)
+    srcField = GRAV.SrcField([rxLoc])
+    survey = GRAV.LinearSurvey(srcField)
+    survey.dobs =  d
+    survey.std =  wd
+    return survey
+
+
+def read_GRAVinv_inp(input_file):
+    """Read input files for forward modeling MAG data with integral form
+    INPUT:
+    input_file: File name containing the forward parameter
+
+    OUTPUT:
+    mshfile
+    obsfile
+    topofile
+    start model
+    ref model
+    weightfile
+    chi_target
+    as, ax ,ay, az
+    upper, lower bounds
+    lp, lqx, lqy, lqz
+
+    # All files should be in the working directory, otherwise the path must
+    # be specified.
+
+    Created on Dec 21th, 2015
+
+    @author: dominiquef
+    """
+
+
+    fid = open(input_file,'r')
+
+    # Line 1
+    line = fid.readline()
+    l_input  = line.split('!')
+    mshfile = l_input[0].rstrip()
+
+    # Line 2
+    line = fid.readline()
+    l_input  = line.split('!')
+    obsfile = l_input[0].rstrip()
+
+    # Line 3
+    line = fid.readline()
+    l_input = re.split('[!\s]',line)
+    if l_input=='null':
+        topofile = []
+
+    else:
+        topofile = l_input[0].rstrip()
+
+
+    # Line 4
+    line = fid.readline()
+    l_input = re.split('[!\s]',line)
+    if l_input[0]=='VALUE':
+        mstart = float(l_input[1])
+
+    else:
+        mstart = l_input[0].rstrip()
+
+    # Line 5
+    line = fid.readline()
+    l_input = re.split('[!\s]',line)
+    if l_input[0]=='VALUE':
+        mref = float(l_input[1])
+
+    else:
+        mref = l_input[0].rstrip()
+
+
+    # Line 7
+    line = fid.readline()
+    l_input = re.split('[!\s]',line)
+    if l_input[0]=='DEFAULT':
+        wgtfile = None
+
+    else:
+        wgtfile = l_input[0].rstrip()
+
+    # Line 8
+    line = fid.readline()
+    l_input = re.split('[!\s]',line)
+    chi = float(l_input[0])
+
+    # Line 9
+    line = fid.readline()
+    l_input = re.split('[!\s]',line)
+    val = np.array(l_input[0:4])
+    alphas = val.astype(np.float)
+
+    # Line 10
+    line = fid.readline()
+    l_input = re.split('[!\s]',line)
+    if l_input[0]=='VALUE':
+        val   = np.array(l_input[1:3])
+        bounds = val.astype(np.float)
+
+    else:
+        bounds = l_input[0].rstrip()
+
+    # Line 11
+    line = fid.readline()
+    l_input = re.split('[!\s]',line)
+    if l_input[0]=='VALUE':
+        val   = np.array(l_input[1:6])
+        lpnorms = val.astype(np.float)
+
+    else:
+        lpnorms = l_input[0].rstrip()
+
+    return mshfile, obsfile, topofile, mstart, mref, wgtfile, chi, alphas, bounds, lpnorms
+

SimPEG/PF/GravityDriver.py

@@ -0,0 +1,295 @@
+import re, os
+from SimPEG import Mesh, np, Utils
+import BaseGrav, Gravity
+
+
+class GravityDriver_Inv(object):
+    """docstring for GravityDriver_Inv"""
+
+    def __init__(self, input_file=None):
+        if input_file is not None:
+            self.basePath = os.path.sep.join(input_file.split(os.path.sep)[:-1])
+            if len(self.basePath) > 0:
+                self.basePath += os.path.sep
+            self.readDriverFile(input_file.split(os.path.sep)[-1])
+
+    def readDriverFile(self, input_file):
+        """
+            Read input files for forward modeling GRAV data with integral form
+            INPUT:
+            input_file: File name containing the forward parameter
+
+            OUTPUT:
+            mshfile
+            obsfile
+            topofile
+            start model
+            ref model
+            active cells model
+            weightfile
+            chi_target
+            as, ax ,ay, az
+            upper, lower bounds
+            lp, lqx, lqy, lqz
+            eps_p, eps_q
+            # All files should be in the working directory, otherwise the path must
+            # be specified.
+
+        """
+
+        fid = open(self.basePath + input_file, 'r')
+
+        # Line 1
+        line = fid.readline()
+        l_input  = line.split('!')
+        mshfile = l_input[0].rstrip()
+
+        # Line 2
+        line = fid.readline()
+        l_input  = line.split('!')
+        obsfile = l_input[0].rstrip()
+
+        # Line 3
+        line = fid.readline()
+        l_input = re.split('[!\s]', line)
+        if l_input=='null':
+            topofile = []
+
+        else:
+            topofile = l_input[0].rstrip()
+
+
+        # Line 4
+        line = fid.readline()
+        l_input = re.split('[!\s]', line)
+        if l_input[0]=='VALUE':
+            mstart = float(l_input[1])
+
+        else:
+            mstart = l_input[0].rstrip()
+
+        # Line 5
+        line = fid.readline()
+        l_input = re.split('[!\s]',line)
+        if l_input[0]=='VALUE':
+            mref = float(l_input[1])
+
+        else:
+            mref = l_input[0].rstrip()
+
+        # Line 6
+        line = fid.readline()
+        l_input = re.split('[!\s]', line)
+        if l_input[0]=='VALUE':
+            staticInput = float(l_input[1])
+
+        elif l_input[0]=='DEFAULT':
+            staticInput = None
+
+        else:
+            staticInput = l_input[0].rstrip()
+
+        # Line 7
+        line = fid.readline()
+        l_input = re.split('[!\s]', line)
+        if l_input=='DEFAULT':
+            wgtfile = []
+
+        else:
+            wgtfile = l_input[0].rstrip()
+
+        # Line 8
+        line = fid.readline()
+        l_input = re.split('[!\s]', line)
+        chi = float(l_input[0])
+
+        # Line 9
+        line = fid.readline()
+        l_input = re.split('[!\s]', line)
+        val = np.array(l_input[0:4])
+        alphas = val.astype(np.float)
+
+        # Line 10
+        line = fid.readline()
+        l_input = re.split('[!\s]', line)
+        if l_input[0]=='VALUE':
+            val   = np.array(l_input[1:3])
+            bounds = val.astype(np.float)
+
+        else:
+            bounds = l_input[0].rstrip()
+
+        # Line 11
+        line = fid.readline()
+        l_input = re.split('[!\s]', line)
+        if l_input[0]=='VALUE':
+            val   = np.array(l_input[1:6])
+            lpnorms = val.astype(np.float)
+
+        else:
+            lpnorms = l_input[0].rstrip()
+
+        # Line 12
+        line = fid.readline()
+        l_input = re.split('[!\s]', line)
+        if l_input[0]=='VALUE':
+            val   = np.array(l_input[1:3])
+            eps = val.astype(np.float)
+
+        else:
+            eps = [None, None]
+
+        self.mshfile  = mshfile
+        self.obsfile  = obsfile
+        self.topofile = topofile
+        self.mstart   = mstart
+        self._mrefInput = mref
+        self._staticInput = staticInput
+        self.wgtfile  = wgtfile
+        self.chi      = chi
+        self.alphas   = alphas
+        self.bounds   = bounds
+        self.lpnorms  = lpnorms
+        self.eps      = eps
+
+    @property
+    def mesh(self):
+        if getattr(self, '_mesh', None) is None:
+            self._mesh = Mesh.TensorMesh.readUBC(self.basePath + self.mshfile)
+        return self._mesh
+
+    @property
+    def survey(self):
+        if getattr(self, '_survey', None) is None:
+            self._survey = self.readGravityObservations(self.basePath + self.obsfile)
+        return self._survey
+
+    @property
+    def activeCells(self):
+        if getattr(self, '_activeCells', None) is None:
+            if self.topofile == 'null':
+                self._activeCells = np.arange(mesh.nC)
+            else:
+                topo = np.genfromtxt(self.basePath + self.topofile, skip_header=1)
+                # Find the active cells
+                active = Utils.surface2ind_topo(self.mesh,topo,'N')
+                inds = np.asarray([inds for inds, elem in enumerate(active, 1) if elem], dtype = int) - 1
+                self._activeCells = inds
+
+        return self._activeCells
+
+    @property
+    def staticCells(self):
+        if getattr(self, '_staticCells', None) is None:
+
+            if getattr(self, '_staticInput', None) is None:
+                # All cells are dynamic: 1's
+                self._dynamicCells = np.arange(len(self.m0))
+                self._staticCells = []
+
+            # Cells with specific value are static: 0's
+            else:
+                if isinstance(self._staticInput, float):
+                    staticCells = self.m0 == self._staticInput
+
+                else:
+                    # Read from file active cells with 0:air, 1:dynamic, -1 static
+                    staticCells = Mesh.TensorMesh.readModelUBC(self.mesh, self.basePath + self._staticInput)
+                    staticCells = staticCells[self.activeCells] == -1
+
+                inds = np.asarray([inds for inds, elem in enumerate(staticCells, 1) if elem], dtype = int) - 1
+                self._staticCells = inds
+
+        return self._staticCells
+
+    @property
+    def dynamicCells(self):
+        if getattr(self, '_dynamicCells', None) is None:
+
+            if getattr(self, '_staticInput', None) is None:
+                # All cells are dynamic: 1's
+                self._dynamicCells = np.arange(len(self.m0))
+
+            # Cells with specific value are static: 0's
+            else:
+                if isinstance(self._staticInput, float):
+                    dynamicCells = self.m0 != self._staticInput
+
+                else:
+                    # Read from file active cells with 0:air, 1:dynamic, -1 static
+                    dynamicCells = Mesh.TensorMesh.readModelUBC(self.mesh, self.basePath + self._staticInput)
+                    dynamicCells = dynamicCells[self.activeCells] == 1
+
+                inds = np.asarray([inds for inds, elem in enumerate(dynamicCells, 1) if elem], dtype = int) - 1
+                self._dynamicCells = inds
+
+        return self._dynamicCells
+
+    @property
+    def nC(self):
+        if getattr(self, '_nC', None) is None:
+            self._nC = len(self.activeCells)
+        return self._nC
+
+    @property
+    def m0(self):
+        if getattr(self, '_m0', None) is None:
+            if isinstance(self.mstart, float):
+                self._m0 = np.ones(self.nC) * self.mstart
+            else:
+
+                self._m0 = Mesh.TensorMesh.readModelUBC(self.mesh, self.basePath + self.mstart)
+                self._m0 = self._m0[self.activeCells]
+
+        return self._m0
+
+    @property
+    def mref(self):
+        if getattr(self, '_mref', None) is None:
+            if isinstance(self._mrefInput, float):
+                self._mref = np.ones(self.nC) * self._mrefInput
+            else:
+                self._mref = Mesh.TensorMesh.readModelUBC(self.mesh, self.basePath + self._mrefInput)
+                self._mref = self._mref[self.activeCells]
+        return self._mref
+
+    def readGravityObservations(self, obs_file):
+        """
+        Read UBC grav file format
+
+        INPUT:
+        :param fileName, path to the UBC obs grav file
+
+        OUTPUT:
+        :param survey
+
+        """
+
+        fid = open(obs_file,'r')
+
+        # First line has the number of rows
+        line = fid.readline()
+        ndat = np.array(line.split(),dtype=int)
+
+        # Pre-allocate space for obsx, obsy, obsz, data, uncert
+        line = fid.readline()
+        temp = np.array(line.split(),dtype=float)
+
+        d  = np.zeros(ndat, dtype=float)
+        wd = np.zeros(ndat, dtype=float)
+        locXYZ = np.zeros( (ndat,3), dtype=float)
+
+        for ii in range(ndat):
+
+            temp = np.array(line.split(),dtype=float)
+            locXYZ[ii,:] = temp[:3]
+            d[ii] = temp[3]
+            wd[ii] = temp[4]
+            line = fid.readline()
+
+        rxLoc = BaseGrav.RxObs(locXYZ)
+        srcField = BaseGrav.SrcField([rxLoc])
+        survey = BaseGrav.LinearSurvey(srcField)
+        survey.dobs =  d
+        survey.std =  wd
+        return survey

SimPEG/PF/MagAnalytics.py

@@ -0,0 +1,278 @@
+from scipy.constants import mu_0
+from SimPEG import *
+from SimPEG.Utils import kron3, speye, sdiag
+import matplotlib.pyplot as plt
+
+
+def spheremodel(mesh, x0, y0, z0, r):
+    """
+        Generate model indicies for sphere
+        - (x0, y0, z0 ): is the center location of sphere
+        - r: is the radius of the sphere
+        - it returns logical indicies of cell-center model
+    """
+    ind = np.sqrt( (mesh.gridCC[:,0]-x0)**2+(mesh.gridCC[:,1]-y0)**2+(mesh.gridCC[:,2]-z0)**2 ) < r
+    return ind
+
+
+def MagSphereAnaFun(x, y, z, R, x0, y0, z0, mu1, mu2, H0, flag='total'):
+    """
+        test
+        Analytic function for Magnetics problem. The set up here is
+        magnetic sphere in whole-space assuming that the inducing field is oriented in the x-direction.
+
+        * (x0,y0,z0)
+        * (x0, y0, z0 ): is the center location of sphere
+        * r: is the radius of the sphere
+
+    .. math::
+
+        \mathbf{H}_0 = H_0\hat{x}
+
+
+    """
+
+    if (~np.size(x)==np.size(y)==np.size(z)):
+        print "Specify same size of x, y, z"
+        return
+    dim = x.shape
+    x = Utils.mkvc(x)
+    y = Utils.mkvc(y)
+    z = Utils.mkvc(z)
+
+    ind = np.sqrt((x-x0)**2+(y-y0)**2+(z-z0)**2 ) < R
+    r = Utils.mkvc(np.sqrt((x-x0)**2+(y-y0)**2+(z-z0)**2 ))
+    Bx = np.zeros(x.size)
+    By = np.zeros(x.size)
+    Bz = np.zeros(x.size)
+
+    # Inside of the sphere
+    rf2 = 3*mu1/(mu2+2*mu1)
+    if flag is 'total' and any(ind):
+        Bx[ind] = mu2*H0*(rf2)
+    elif (flag == 'secondary'):
+        Bx[ind] = mu2*H0*(rf2)-mu1*H0
+
+    By[ind] = 0.
+    Bz[ind] = 0.
+    # Outside of the sphere
+    rf1 = (mu2-mu1)/(mu2+2*mu1)
+    if (flag == 'total'):
+        Bx[~ind] = mu1*(H0+H0/r[~ind]**5*(R**3)*rf1*(2*(x[~ind]-x0)**2-(y[~ind]-y0)**2-(z[~ind]-z0)**2))
+    elif (flag == 'secondary'):
+        Bx[~ind] = mu1*(H0/r[~ind]**5*(R**3)*rf1*(2*(x[~ind]-x0)**2-(y[~ind]-y0)**2-(z[~ind]-z0)**2))
+
+    By[~ind] = mu1*(H0/r[~ind]**5*(R**3)*rf1*(3*(x[~ind]-x0)*(y[~ind]-y0)))
+    Bz[~ind] = mu1*(H0/r[~ind]**5*(R**3)*rf1*(3*(x[~ind]-x0)*(z[~ind]-z0)))
+    return np.reshape(Bx, x.shape, order='F'), np.reshape(By, x.shape, order='F'), np.reshape(Bz, x.shape, order='F')
+
+
+def CongruousMagBC(mesh, Bo, chi):
+    """
+        Computing boundary condition using Congrous sphere method.
+        This is designed for secondary field formulation.
+
+        >> Input
+
+        * mesh:   Mesh class
+        * Bo:     np.array([Box, Boy, Boz]): Primary magnetic flux
+        * chi:    susceptibility at cell volume
+
+        .. math::
+
+            \\vec{B}(r) = \\frac{\mu_0}{4\pi} \\frac{m}{ \| \\vec{r} - \\vec{r}_0\|^3}[3\hat{m}\cdot\hat{r}-\hat{m}]
+
+    """
+
+    ind = chi > 0.
+    V = mesh.vol[ind].sum()
+
+    gamma = 1/V*(chi*mesh.vol).sum() # like a mass!
+
+    Bot = np.sqrt(sum(Bo**2))
+    mx = Bo[0]/Bot
+    my = Bo[1]/Bot
+    mz = Bo[2]/Bot
+
+    mom = 1/mu_0*Bot*gamma*V/(1+gamma/3)
+    xc = sum(chi[ind]*mesh.gridCC[:,0][ind])/sum(chi[ind])
+    yc = sum(chi[ind]*mesh.gridCC[:,1][ind])/sum(chi[ind])
+    zc = sum(chi[ind]*mesh.gridCC[:,2][ind])/sum(chi[ind])
+
+    indxd, indxu, indyd, indyu, indzd, indzu =  mesh.faceBoundaryInd
+
+    const = mu_0/(4*np.pi)*mom
+    rfun = lambda x: np.sqrt((x[:,0]-xc)**2 + (x[:,1]-yc)**2 + (x[:,2]-zc)**2)
+
+    mdotrx = (mx*(mesh.gridFx[(indxd|indxu),0]-xc)/rfun(mesh.gridFx[(indxd|indxu),:]) +
+              my*(mesh.gridFx[(indxd|indxu),1]-yc)/rfun(mesh.gridFx[(indxd|indxu),:]) +
+              mz*(mesh.gridFx[(indxd|indxu),2]-zc)/rfun(mesh.gridFx[(indxd|indxu),:]))
+
+    Bbcx = const/(rfun(mesh.gridFx[(indxd|indxu),:])**3)*(3*mdotrx*(mesh.gridFx[(indxd|indxu),0]-xc)/rfun(mesh.gridFx[(indxd|indxu),:])-mx)
+
+    mdotry = (mx*(mesh.gridFy[(indyd|indyu),0]-xc)/rfun(mesh.gridFy[(indyd|indyu),:]) +
+              my*(mesh.gridFy[(indyd|indyu),1]-yc)/rfun(mesh.gridFy[(indyd|indyu),:]) +
+              mz*(mesh.gridFy[(indyd|indyu),2]-zc)/rfun(mesh.gridFy[(indyd|indyu),:]))
+
+    Bbcy = const/(rfun(mesh.gridFy[(indyd|indyu),:])**3)*(3*mdotry*(mesh.gridFy[(indyd|indyu),1]-yc)/rfun(mesh.gridFy[(indyd|indyu),:])-my)
+
+    mdotrz = (mx*(mesh.gridFz[(indzd|indzu),0]-xc)/rfun(mesh.gridFz[(indzd|indzu),:])  +
+              my*(mesh.gridFz[(indzd|indzu),1]-yc)/rfun(mesh.gridFz[(indzd|indzu),:]) +
+              mz*(mesh.gridFz[(indzd|indzu),2]-zc)/rfun(mesh.gridFz[(indzd|indzu),:]))
+
+    Bbcz = const/(rfun(mesh.gridFz[(indzd|indzu),:])**3)*(3*mdotrz*(mesh.gridFz[(indzd|indzu),2]-zc)/rfun(mesh.gridFz[(indzd|indzu),:])-mz)
+
+    return np.r_[Bbcx, Bbcy, Bbcz], (1/gamma-1/(3+gamma))*1/V
+
+
+def MagSphereAnaFunA(x, y, z, R, xc, yc, zc, chi, Bo, flag):
+    """
+        Computing boundary condition using Congrous sphere method.
+        This is designed for secondary field formulation.
+        >> Input
+        mesh:   Mesh class
+        Bo:     np.array([Box, Boy, Boz]): Primary magnetic flux
+        Chi:    susceptibility at cell volume
+
+        .. math::
+
+            \\vec{B}(r) = \\frac{\mu_0}{4\pi}\\frac{m}{\| \\vec{r}-\\vec{r}_0\|^3}[3\hat{m}\cdot\hat{r}-\hat{m}]
+
+    """
+    if (~np.size(x)==np.size(y)==np.size(z)):
+        print "Specify same size of x, y, z"
+        return
+    dim = x.shape
+    x = Utils.mkvc(x)
+    y = Utils.mkvc(y)
+    z = Utils.mkvc(z)
+
+    Bot = np.sqrt(sum(Bo**2))
+    mx = Bo[0]/Bot
+    my = Bo[1]/Bot
+    mz = Bo[2]/Bot
+
+    ind = np.sqrt((x-xc)**2+(y-yc)**2+(z-zc)**2 ) < R
+
+    Bx = np.zeros(x.size)
+    By = np.zeros(x.size)
+    Bz = np.zeros(x.size)
+
+    # Inside of the sphere
+    rf2 = 3/(chi+3)*(1+chi)
+    if (flag == 'total'):
+        Bx[ind] = Bo[0]*(rf2)
+        By[ind] = Bo[1]*(rf2)
+        Bz[ind] = Bo[2]*(rf2)
+    elif (flag == 'secondary'):
+        Bx[ind] = Bo[0]*(rf2)-Bo[0]
+        By[ind] = Bo[1]*(rf2)-Bo[1]
+        Bz[ind] = Bo[2]*(rf2)-Bo[2]
+
+    r = Utils.mkvc(np.sqrt((x-xc)**2+(y-yc)**2+(z-zc)**2 ))
+    V = 4*np.pi*R**3/3
+    mom = Bot/mu_0*chi/(1+chi/3)*V
+    const = mu_0/(4*np.pi)*mom
+    mdotr = (mx*(x[~ind]-xc)/r[~ind] + my*(y[~ind]-yc)/r[~ind] + mz*(z[~ind]-zc)/r[~ind])
+    Bx[~ind] = const/(r[~ind]**3)*(3*mdotr*(x[~ind]-xc)/r[~ind]-mx)
+    By[~ind] = const/(r[~ind]**3)*(3*mdotr*(y[~ind]-yc)/r[~ind]-my)
+    Bz[~ind] = const/(r[~ind]**3)*(3*mdotr*(z[~ind]-zc)/r[~ind]-mz)
+
+
+    return Bx, By, Bz
+
+
+def IDTtoxyz(Inc, Dec, Btot):
+    """
+        Convert from Inclination, Declination, Total  intensity of earth field to x, y, z
+    """
+    Bx = Btot*np.cos(Inc/180.*np.pi)*np.sin(Dec/180.*np.pi)
+    By = Btot*np.cos(Inc/180.*np.pi)*np.cos(Dec/180.*np.pi)
+    Bz = -Btot*np.sin(Inc/180.*np.pi)
+
+    return np.r_[Bx, By, Bz]
+
+
+def MagSphereFreeSpace(x, y, z, R, xc, yc, zc, chi, Bo):
+    """
+        Computing boundary condition using Congrous sphere method.
+        This is designed for secondary field formulation.
+        >> Input
+        mesh:   Mesh class
+        Bo:     np.array([Box, Boy, Boz]): Primary magnetic flux
+        Chi:    susceptibility at cell volume
+
+        .. math::
+
+            \\vec{B}(r) = \\frac{\mu_0}{4\pi}\\frac{m}{\| \\vec{r}-\\vec{r}_0\|^3}[3\hat{m}\cdot\hat{r}-\hat{m}]
+
+    """
+    if (~np.size(x)==np.size(y)==np.size(z)):
+        print "Specify same size of x, y, z"
+        return
+
+    x = Utils.mkvc(x)
+    y = Utils.mkvc(y)
+    z = Utils.mkvc(z)
+
+    nobs = len(x)
+
+    Bot = np.sqrt(sum(Bo**2))
+
+    mx = np.ones([nobs]) * Bo[0,0] *  R**3 / 3. * chi
+    my = np.ones([nobs]) * Bo[0,1] *  R**3 / 3. * chi
+    mz = np.ones([nobs]) * Bo[0,2] *  R**3 / 3. * chi
+
+    M = np.c_[mx, my, mz]
+
+    rx = (x - xc)
+    ry = (y - yc)
+    rz = (zc - z)
+
+    rvec = np.c_[rx, ry, rz]
+    r = np.sqrt((rx)**2+(ry)**2+(rz)**2 )
+
+    B = -Utils.sdiag(1./r**3)*M + Utils.sdiag((3 * np.sum(M*rvec,axis=1))/r**5)*rvec
+
+    Bx = B[:,0]
+    By = B[:,1]
+    Bz = B[:,2]
+
+    return Bx, By, Bz
+
+if __name__ == '__main__':
+
+    hxind = [(0,25,1.3),(21, 12.5),(0,25,1.3)]
+    hyind = [(0,25,1.3),(21, 12.5),(0,25,1.3)]
+    hzind = [(0,25,1.3),(20, 12.5),(0,25,1.3)]
+    # hx, hy, hz = Utils.meshTensors(hxind, hyind, hzind)
+    M3 = Mesh.TensorMesh([hxind, hyind, hzind], "CCC")
+    indxd, indxu, indyd, indyu, indzd, indzu = M3.faceBoundaryInd
+    mu0 = 4*np.pi*1e-7
+    chibkg = 0.
+    chiblk = 0.01
+    chi = np.ones(M3.nC)*chibkg
+    sph_ind = spheremodel(M3, 0, 0, 0, 100)
+    chi[sph_ind] = chiblk
+    mu = (1.+chi)*mu0
+    Bbc, const = CongruousMagBC(M3, np.array([1., 0., 0.]), chi)
+
+    flag = 'secondary'
+    Box = 1.
+    H0 = Box/mu_0
+    Bbcxx, Bbcxy, Bbcxz  = MagSphereAnaFun(M3.gridFx[(indxd|indxu),0], M3.gridFx[(indxd|indxu),1], M3.gridFx[(indxd|indxu),2], 100, 0., 0., 0., mu_0, mu_0*(1+chiblk), H0, flag)
+    Bbcyx, Bbcyy, Bbcyz  = MagSphereAnaFun(M3.gridFy[(indyd|indyu),0], M3.gridFy[(indyd|indyu),1], M3.gridFy[(indyd|indyu),2], 100, 0., 0., 0., mu_0, mu_0*(1+chiblk), H0, flag)
+    Bbczx, Bbczy, Bbczz  = MagSphereAnaFun(M3.gridFz[(indzd|indzu),0], M3.gridFz[(indzd|indzu),1], M3.gridFz[(indzd|indzu),2], 100, 0., 0., 0., mu_0, mu_0*(1+chiblk), H0, flag)
+    Bbc_ana = np.r_[Bbcxx, Bbcyy, Bbczz]
+
+    # fig, ax = plt.subplots(1,1, figsize = (10, 10))
+    # ax.plot(Bbc_ana)
+    # ax.plot(Bbc)
+    # plt.show()
+    err = np.linalg.norm(Bbc-Bbc_ana)/np.linalg.norm(Bbc_ana)
+
+    if err < 0.1:
+        print 'Mag Boundary computation is valid, err = ', err
+    else:
+        print 'Mag Boundary computation is wrong!!, err = ', err
+    pass

SimPEG/PF/MagneticsDriver.py

@@ -0,0 +1,334 @@
+import re, os
+from SimPEG import Mesh, np, Utils
+import BaseMag, Magnetics
+
+class MagneticsDriver_Inv(object):
+    """docstring for MagneticsDriver_Inv"""
+
+    def __init__(self, input_file=None):
+        if input_file is not None:
+            self.basePath = os.path.sep.join(input_file.split(os.path.sep)[:-1])
+            if len(self.basePath) > 0:
+                self.basePath += os.path.sep
+            self.readDriverFile(input_file.split(os.path.sep)[-1])
+
+
+    def readDriverFile(self, input_file):
+        """
+            Read input files for forward modeling MAG data with integral form
+            INPUT:
+            input_file: File name containing the forward parameter
+
+            OUTPUT:
+            mshfile
+            obsfile
+            topofile
+            start model
+            ref model
+            mag model
+            weightfile
+            chi_target
+            as, ax ,ay, az
+            upper, lower bounds
+            lp, lqx, lqy, lqz
+
+            # All files should be in the working directory, otherwise the path must
+            # be specified.
+
+        """
+
+
+        fid = open(self.basePath + input_file,'r')
+
+        # Line 1
+        line = fid.readline()
+        l_input  = line.split('!')
+        mshfile = l_input[0].rstrip()
+
+        # Line 2
+        line = fid.readline()
+        l_input  = line.split('!')
+        obsfile = l_input[0].rstrip()
+
+        # Line 3
+        line = fid.readline()
+        l_input = re.split('[!\s]',line)
+        if l_input=='null':
+            topofile = []
+
+        else:
+            topofile = l_input[0].rstrip()
+
+
+        # Line 4
+        line = fid.readline()
+        l_input = re.split('[!\s]',line)
+        if l_input[0]=='VALUE':
+            mstart = float(l_input[1])
+
+        else:
+            mstart = l_input[0].rstrip()
+
+        # Line 5
+        line = fid.readline()
+        l_input = re.split('[!\s]',line)
+        if l_input[0]=='VALUE':
+            mref = float(l_input[1])
+
+        else:
+            mref = l_input[0].rstrip()
+
+        # Line 6
+        line = fid.readline()
+        l_input = re.split('[!\s]',line)
+        if l_input[0]=='VALUE':
+            staticInput = float(l_input[1])
+
+        elif l_input[0]=='DEFAULT':
+            staticInput = None
+
+        else:
+            staticInput = l_input[0].rstrip()
+
+
+        # Line 7
+        line = fid.readline()
+        l_input = re.split('[!\s]',line)
+        if l_input=='DEFAULT':
+            magfile = []
+
+        else:
+            magfile = l_input[0].rstrip()
+
+        # Line 8
+        line = fid.readline()
+        l_input = re.split('[!\s]',line)
+        if l_input=='DEFAULT':
+            wgtfile = []
+
+        else:
+            wgtfile = l_input[0].rstrip()
+
+        # Line 9
+        line = fid.readline()
+        l_input = re.split('[!\s]',line)
+        chi = float(l_input[0])
+
+        # Line 10
+        line = fid.readline()
+        l_input = re.split('[!\s]',line)
+        val = np.array(l_input[0:4])
+        alphas = val.astype(np.float)
+
+        # Line 11
+        line = fid.readline()
+        l_input = re.split('[!\s]',line)
+        if l_input[0]=='VALUE':
+            val   = np.array(l_input[1:3])
+            bounds = val.astype(np.float)
+
+        else:
+            bounds = l_input[0].rstrip()
+
+        # Line 12
+        line = fid.readline()
+        l_input = re.split('[!\s]',line)
+        if l_input[0]=='VALUE':
+            val   = np.array(l_input[1:6])
+            lpnorms = val.astype(np.float)
+
+        else:
+            lpnorms = l_input[0].rstrip()
+
+        # Line 13
+        line = fid.readline()
+        l_input = re.split('[!\s]',line)
+        if l_input[0]=='VALUE':
+            val   = np.array(l_input[1:3])
+            eps = val.astype(np.float)
+
+        else:
+            eps = [None,None]
+
+        self.mshfile  = mshfile
+        self.obsfile  = obsfile
+        self.topofile = topofile
+        self.mstart   = mstart
+        self._mrefInput = mref
+        self._staticInput = staticInput
+        self.magfile  = magfile
+        self.wgtfile  = wgtfile
+        self.chi      = chi
+        self.alphas   = alphas
+        self.bounds   = bounds
+        self.lpnorms  = lpnorms
+        self.eps      = eps
+
+    @property
+    def mesh(self):
+        if getattr(self, '_mesh', None) is None:
+            self._mesh = Mesh.TensorMesh.readUBC(self.basePath + self.mshfile)
+        return self._mesh
+
+    @property
+    def survey(self):
+        if getattr(self, '_survey', None) is None:
+            self._survey = self.readMagneticsObservations(self.obsfile)
+        return self._survey
+
+    @property
+    def activeCells(self):
+        if getattr(self, '_activeCells', None) is None:
+            if self.topofile == 'null':
+                self._activeCells = np.arange(self.mesh.nC)
+            else:
+                topo = np.genfromtxt(self.basePath + self.topofile, skip_header=1)
+                # Find the active cells
+                active = Utils.surface2ind_topo(self.mesh,topo,'N')
+                inds = np.asarray([inds for inds, elem in enumerate(active, 1) if elem], dtype = int) - 1
+                self._activeCells = inds
+
+        return self._activeCells
+
+    @property
+    def staticCells(self):
+        if getattr(self, '_staticCells', None) is None:
+
+            if getattr(self, '_staticInput', None) is None:
+                # All cells are dynamic: 1's
+                self._dynamicCells = np.arange(len(self.m0))
+                self._staticCells = []
+
+            # Cells with specific value are static: 0's
+            else:
+                if isinstance(self._staticInput, float):
+                    staticCells = self.m0 == self._staticInput
+
+                else:
+                    # Read from file active cells with 0:air, 1:dynamic, -1 static
+                    staticCells = Mesh.TensorMesh.readModelUBC(self.mesh, self.basePath + self._staticInput)
+                    staticCells = staticCells[self.activeCells] == -1
+
+                inds = np.asarray([inds for inds, elem in enumerate(staticCells, 1) if elem], dtype = int) - 1
+                self._staticCells = inds
+
+        return self._staticCells
+
+    @property
+    def dynamicCells(self):
+        if getattr(self, '_dynamicCells', None) is None:
+
+            if getattr(self, '_staticInput', None) is None:
+                # All cells are dynamic: 1's
+                self._dynamicCells = np.arange(len(self.m0))
+
+            # Cells with specific value are static: 0's
+            else:
+                if isinstance(self._staticInput, float):
+                    dynamicCells = self.m0 != self._staticInput
+
+                else:
+                    # Read from file active cells with 0:air, 1:dynamic, -1 static
+                    dynamicCells = Mesh.TensorMesh.readModelUBC(self.mesh, self.basePath + self._staticInput)
+                    dynamicCells = dynamicCells[self.activeCells] == 1
+
+                inds = np.asarray([inds for inds, elem in enumerate(dynamicCells, 1) if elem], dtype = int) - 1
+                self._dynamicCells = inds
+
+        return self._dynamicCells
+
+    @property
+    def nC(self):
+        if getattr(self, '_nC', None) is None:
+            self._nC = len(self.activeCells)
+        return self._nC
+
+    @property
+    def m0(self):
+        if getattr(self, '_m0', None) is None:
+            if isinstance(self.mstart, float):
+                self._m0 = np.ones(self.nC) * self.mstart
+            else:
+                self._m0 = Mesh.TensorMesh.readModelUBC(self.mesh,self.basePath + self.mstart)
+                self._m0 = self._m0[self.activeCells]
+
+        return self._m0
+
+    @property
+    def mref(self):
+        if getattr(self, '_mref', None) is None:
+            if isinstance(self._mrefInput, float):
+                self._mref = np.ones(self.nC) * self._mrefInput
+            else:
+                self._mref = Mesh.TensorMesh.readModelUBC(self.mesh,self.basePath + self._mrefInput)
+                self._mref = self._mref[self.activeCells]
+        return self._mref
+
+
+    @property
+    def magnetizationModel(self):
+        """
+            magnetization vector
+        """
+
+        if self.magfile == 'DEFAULT':
+            return Magnetics.dipazm_2_xyz(np.ones(self.nC) * self.survey.srcField.param[1], np.ones(self.nC) * self.survey.srcField.param[2])
+
+        else:
+            raise NotImplementedError("this will require you to read in a three column vector model")
+            self._mref = Utils.meshutils.readUBCTensorModel(self.basePath + self._mrefInput, self.mesh)
+            return np.genfromtxt(self.magfile,delimiter=' \n',dtype=np.str,comments='!')
+
+    def readMagneticsObservations(self, obs_file):
+        """
+            Read and write UBC mag file format
+
+            INPUT:
+            :param fileName, path to the UBC obs mag file
+
+            OUTPUT:
+            :param survey
+            :param M, magnetization orentiaton (MI, MD)
+        """
+
+        fid = open(self.basePath + obs_file,'r')
+
+        # First line has the inclination,declination and amplitude of B0
+        line = fid.readline()
+        B = np.array(line.split(),dtype=float)
+
+        # Second line has the magnetization orientation and a flag
+        line = fid.readline()
+        M = np.array(line.split(),dtype=float)
+
+        # Third line has the number of rows
+        line = fid.readline()
+        ndat = np.array(line.split(),dtype=int)
+
+        # Pre-allocate space for obsx, obsy, obsz, data, uncert
+        line = fid.readline()
+        temp = np.array(line.split(),dtype=float)
+
+        d  = np.zeros(ndat, dtype=float)
+        wd = np.zeros(ndat, dtype=float)
+        locXYZ = np.zeros( (ndat,3), dtype=float)
+
+        for ii in range(ndat):
+
+            temp = np.array(line.split(),dtype=float)
+            locXYZ[ii,:] = temp[:3]
+
+            if len(temp) > 3:
+                d[ii] = temp[3]
+
+                if len(temp)==5:
+                    wd[ii] = temp[4]
+
+            line = fid.readline()
+
+        rxLoc = BaseMag.RxObs(locXYZ)
+        srcField = BaseMag.SrcField([rxLoc],param=(B[2],B[0],B[1]))
+        survey = BaseMag.LinearSurvey(srcField)
+        survey.dobs =  d
+        survey.std =  wd
+        return survey

SimPEG/PF/__init__.py

@@ -0,0 +1,7 @@
+import MagAnalytics
+import BaseMag
+import Magnetics
+import BaseGrav
+import Gravity
+import MagneticsDriver
+import GravityDriver

docs/conf.py

@@ -289,4 +289,6 @@ nitpick_ignore = [
     ('py:class', 'SurveyDC'),
     ('py:class', 'BaseMTFields'),
     ('py:class', 'SolverLU'),
+    ('py:class', 'BaseMagSurvey'),
+    ('py:class', 'BaseMagMap'),
 ]

docs/content/examples/PF_Magnetics_Analytics.rst

@@ -0,0 +1,26 @@
+.. _examples_PF_Magnetics_Analytics:
+
+.. --------------------------------- ..
+..                                   ..
+..    THIS FILE IS AUTO GENEREATED   ..
+..                                   ..
+..    SimPEG/Examples/__init__.py    ..
+..                                   ..
+.. --------------------------------- ..
+
+
+PF: Magnetics: Analytics
+========================
+
+Comparing the magnetics field in Vancouver to Seoul
+
+
+
+.. plot::
+
+    from SimPEG import Examples
+    Examples.PF_Magnetics_Analytics.run()
+
+.. literalinclude:: ../../../SimPEG/Examples/PF_Magnetics_Analytics.py
+    :language: python
+    :linenos:

docs/content/pf/index.rst

@@ -0,0 +1,170 @@
+.. math::
+
+    \renewcommand{\div}{\nabla\cdot\,}
+    \newcommand{\grad}{\vec \nabla}
+    \newcommand{\curl}{{\vec \nabla}\times\,}
+    \newcommand {\J}{{\vec J}}
+    \renewcommand{\H}{{\vec H}}
+    \newcommand {\E}{{\vec E}}
+    \newcommand{\dcurl}{{\mathbf C}}
+    \newcommand{\dgrad}{{\mathbf G}}
+    \newcommand{\Acf}{{\mathbf A_c^f}}
+    \newcommand{\Ace}{{\mathbf A_c^e}}
+    \renewcommand{\S}{{\mathbf \Sigma}}
+    \renewcommand{\Div}{{\mathbf {Div}}}
+    \newcommand{\St}{{\mathbf \Sigma_\tau}}
+    \newcommand{\T}{{\mathbf T}}
+    \newcommand{\Tt}{{\mathbf T_\tau}}
+    \newcommand{\diag}{\mathbf{diag}}
+    \newcommand{\M}{{\mathbf M}}
+    \newcommand{\MfMui}{{\M^f_{\mu^{-1}}}}
+    \newcommand{\dMfMuI}{{d_m (\M^f_{\mu^{-1}})^{-1}}}
+    \newcommand{\MeSig}{{\M^e_\sigma}}
+    \newcommand{\MeSigInf}{{\M^e_{\sigma_\infty}}}
+    \newcommand{\MeSigO}{{\M^e_{\sigma_0}}}
+    \newcommand{\Me}{{\M^e}}
+    \newcommand{\Mes}[1]{{\M^e_{#1}}}
+    \newcommand{\Mee}{{\M^e_e}}
+    \newcommand{\Mej}{{\M^e_j}}
+    \newcommand{\BigO}[1]{\mathcal{O}\bigl(#1\bigr)}
+    \newcommand{\bE}{\mathbf{E}}
+    \newcommand{\bH}{\mathbf{H}}
+    \newcommand{\B}{\vec{B}}
+    \newcommand{\D}{\vec{D}}
+    \renewcommand{\H}{\vec{H}}
+    \newcommand{\s}{\vec{s}}
+    \newcommand{\bfJ}{\bf{J}}
+    \newcommand{\vecm}{\vec m}
+    \renewcommand{\Re}{\mathsf{Re}}
+    \renewcommand{\Im}{\mathsf{Im}}
+    \renewcommand {\j}  { {\vec j} }
+    \newcommand {\h}  { {\vec h} }
+    \renewcommand {\b}  { {\vec b} }
+    \newcommand {\e}  { {\vec e} }
+    \newcommand {\c}  { {\vec c} }
+    \renewcommand {\d}  { {\vec d} }
+    \renewcommand {\u}  { {\vec u} }
+    \newcommand{\I}{\vec{I}}
+
+Magnetics
+*********
+
+
+The geomagnetic field can be ranked as the longest studied of all the geophysical properties of the earth. In addition, magnetic survey, has been used broadly in diverse realm e.g., mining, oil and gas industry and environmental engineering. Although, this geophysical application is quite common in geoscience; however, we do not have modular, well-documented and well-tested open-source codes, which perform forward and inverse problems of magnetic survey. Therefore, here we are going to build up magnetic forward and inverse modeling code based on  two common methodologies for forward problem - differential equation and integral equation approaches. \
+
+First, we start with some backgrounds of magnetics, e.g., Maxwell's equations. Based on that secondly, we use differential equation approach to solve forward problem with secondary field formulation. In order to discretzie our system here, we use finite volume approach with weak formulation. Third, we solve inverse problem through Gauss-Newton method.
+
+Backgrounds
+===========
+Maxwell's equations for static case with out current source can be written as
+
+.. math::
+
+    \nabla U = \frac{1}{\mu}\vec{B} \\
+
+    \nabla \cdot \vec{B} = 0
+
+where \\(\\vec{B}\\) is magnetic flux (\\(\T\\)) and \\(\U\\) is magnetic potential and \\(\\mu\\) is permeability. Since we do not have any source term in above equations, boundary condition is going to be the driving force of our system as given below
+
+.. math::
+
+    (\vec{B}\cdot{\vec{n}})_{\partial\Omega} = B_{BC}
+
+where \\(\\vec{n}\\) means the unit normal vector on the boundary surface (\\(\\partial \\Omega\\)). By using seocondary field formulation we can rewrite above equations as
+
+.. math::
+
+    \frac{1}{\mu}\vec{B}_s = (\frac{1}{\mu}_0-\frac{1}{\mu})\vec{B}_0+\nabla\phi_s
+
+    \nabla \cdot \vec{B}_s = 0
+
+    (\vec{B}_s\cdot{\vec{n}})_{\partial\Omega} = B_{sBC}
+
+where \\(\\vec{B}_s\\) is the secondary magnetic flux and \\(\\vec{B}_0\\) is the background or primary magnetic flux. In practice, we consider our earth field, which we can get from International Geomagnetic Reference Field (IGRF) by specifying the time and location, as  \\(\\vec{B}_0\\). And based on this background fields, we compute secondary fields (\\(\\vec{B}_s\\)). Now we introduce the susceptibility as
+
+.. math::
+
+    \chi = \frac{\mu}{\mu_0} - 1 \\
+
+    \mu = \mu_0(1+\chi)
+
+Since most materials in the earth  have lower permeability than \\(\\mu_0\\), usually \\(\\chi\\) is greater than 0.
+
+.. note::
+
+    Actually, this is an assumption, which means we are not sure exactly this is true, although we are sure, it is very rare that we can encounter those materials. Anyway, practical range of the susceptibility is \\(0 < \\chi < 1 \\).
+
+Since we compute secondary field based on the earth field, which can be different from different locations in the world, we can expect different anomalous responses in different locations in the earth. For instance, assume we have two susceptible spheres, which are exactly same. However, anomalous responses in Seoul and Vancouver are going to be different.
+
+.. plot::
+
+    from SimPEG import Examples
+    Examples.PF_Magnetics_Analytics.run()
+
+Since we can measure total fields ( \\(\\vec{B}\\)), and usually have reasonably accurate earth field (\\(\\vec{B}_0\\)), we can compute anomalous fields, \\(\\vec{B}_s\\) from our observed data. If you want to download earth magnetic fields at specific location see this website (`noaa <http://www.ngdc.noaa.gov/geomag-web/>`_).
+
+What is our data?
+-----------------
+
+In applied geophysics, which means in practice, it is common to refer to measurements as "the magnetic anomaly" and we can consider this as our observed data. For further descriptions in `GPG <http://www.eos.ubc.ca/courses/eosc350/content/>`_ materials for magnetic survey. Now we have the simple relation ship between "the magnetic anomaly" and the total field as
+
+.. math::
+
+	\triangle\vec{B} = |\hat{B}_o-\vec{B}_s|-|\hat{B}_o| \approx |\vec{B}_s|cos \theta
+
+where \\(\\theta\\) is the angle between total and anomalous fields, \\(\\hat{B}_o\\) is the unit vector for \\(\\vec{B}_o\\). Equivalently, we can use the vector dot product to show that the anomalous field is approximately equal to the projection of that field onto the direction of the inducing field. Using this approach we would write
+
+.. math::
+
+	\triangle\vec{B} = |\vec{B}_s|cos \theta = |\hat{B}_o||\vec{B}_s|cos \theta = \hat{B}_o \cdot \vec{B}_s
+
+This is important because, in practice we usually use a total field magnetometer (like a proton precession or optically pumped sensor), which can measure only that part of the anomalous field which is in the direction of the earth's main field.
+
+Sphere in a whole space
+-----------------------
+
+
+Forward problem
+===============
+
+Differential equation approach
+------------------------------
+
+    .. math::
+
+        \mathbf{A}\mathbf{u} = \mathbf{rhs}
+
+        \mathbf{A} =  \Div(\MfMui)^{-1}\Div^{T}
+
+        \mathbf{rhs} =  \Div(\MfMui)^{-1}\mathbf{M}^f_{\mu_0^{-1}}\mathbf{B}_0 - \Div\mathbf{B}_0+\diag(v)\mathbf{D} \mathbf{P}_{out}^T \mathbf{B}_{sBC}
+
+        \mathbf{B}_s = (\MfMui)^{-1}\mathbf{M}^f_{\mu_0^{-1}}\mathbf{B}_0-\mathbf{B}_0 -(\MfMui)^{-1}\Div^T \mathbf{u}
+
+
+Mag Differential eq. approach
+=============================
+
+.. autoclass:: SimPEG.PF.Magnetics.Problem3D_DiffSecondary
+    :show-inheritance:
+    :members:
+    :undoc-members:
+    :inherited-members:
+
+.. autoclass:: SimPEG.PF.BaseMag.BaseMagSurvey
+    :show-inheritance:
+    :members:
+    :undoc-members:
+    :inherited-members:
+
+
+Mag Integral eq. approach
+=========================
+
+Mag analytic solutions
+----------------------
+
+.. automodule:: SimPEG.PF.MagAnalytics
+    :show-inheritance:
+    :members:
+    :undoc-members:
+    :inherited-members:

docs/index.rst

@@ -62,6 +62,7 @@ Packages
    content/dc/index
    content/ip/index
    content/mt/index
+   content/pf/index
    content/flow/index
 
 Finite Volume

tests/flow/test_Richards.py

@@ -14,6 +14,7 @@ TOL = 1E-8
 
 np.random.seed(0)
 
+
 class TestModels(unittest.TestCase):
 
     def test_BaseHaverkamp_Theta(self):

tests/mesh/test_MeshIO.py

@@ -1,5 +1,6 @@
 import numpy as np
-import unittest, os
+import unittest
+import os
 import SimPEG as simpeg
 from SimPEG.Mesh import TensorMesh, TreeMesh
 
@@ -8,7 +9,7 @@ class TestTensorMeshIO(unittest.TestCase):
 
     def setUp(self):
         h = np.ones(16)
-        mesh = TensorMesh([h,2*h,3*h])
+        mesh = TensorMesh([h, 2*h, 3*h])
         self.mesh = mesh
 
     def test_UBCfiles(self):
@@ -17,9 +18,9 @@ class TestTensorMeshIO(unittest.TestCase):
         # Make a vector
         vec = np.arange(mesh.nC)
         # Write and read
-        mesh.writeUBC('temp.msh', {'arange.txt':vec})
+        mesh.writeUBC('temp.msh', {'arange.txt': vec})
         meshUBC = TensorMesh.readUBC('temp.msh')
-        vecUBC  = meshUBC.readModelUBC('arange.txt')
+        vecUBC = meshUBC.readModelUBC('arange.txt')
 
         # The mesh
         assert mesh.__str__() == meshUBC.__str__()
@@ -27,12 +28,11 @@ class TestTensorMeshIO(unittest.TestCase):
         assert np.sum(vec - vecUBC) == 0
         assert np.all(np.array(mesh.h) - np.array(meshUBC.h) == 0)
 
-
-        vecUBC  = mesh.readModelUBC('arange.txt')
+        vecUBC = mesh.readModelUBC('arange.txt')
         assert np.sum(vec - vecUBC) == 0
 
         mesh.writeModelUBC('arange2.txt', vec + 1)
-        vec2UBC  = mesh.readModelUBC('arange2.txt')
+        vec2UBC = mesh.readModelUBC('arange2.txt')
         assert np.sum(vec + 1 - vec2UBC) == 0
 
         print 'IO of UBC tensor mesh files is working'
@@ -44,7 +44,7 @@ class TestTensorMeshIO(unittest.TestCase):
         mesh = self.mesh
         vec = np.arange(mesh.nC)
 
-        mesh.writeVTK('temp.vtr', {'arange.txt':vec})
+        mesh.writeVTK('temp.vtr', {'arange.txt': vec})
         meshVTR, models = TensorMesh.readVTK('temp.vtr')
 
         assert mesh.__str__() == meshVTR.__str__()
@@ -62,10 +62,10 @@ class TestOcTreeMeshIO(unittest.TestCase):
 
     def setUp(self):
         h = np.ones(16)
-        mesh = TreeMesh([h,2*h,3*h])
+        mesh = TreeMesh([h, 2*h, 3*h])
         mesh.refine(3)
-        mesh._refineCell([0,0,0,3])
-        mesh._refineCell([0,2,0,3])
+        mesh._refineCell([0, 0, 0, 3])
+        mesh._refineCell([0, 2, 0, 3])
         self.mesh = mesh
 
     def test_UBCfiles(self):
@@ -74,9 +74,9 @@ class TestOcTreeMeshIO(unittest.TestCase):
         # Make a vector
         vec = np.arange(mesh.nC)
         # Write and read
-        mesh.writeUBC('temp.msh', {'arange.txt':vec})
+        mesh.writeUBC('temp.msh', {'arange.txt': vec})
         meshUBC = TreeMesh.readUBC('temp.msh')
-        vecUBC  = meshUBC.readModelUBC('arange.txt')
+        vecUBC = meshUBC.readModelUBC('arange.txt')
 
         # The mesh
         assert mesh.__str__() == meshUBC.__str__()
@@ -90,11 +90,10 @@ class TestOcTreeMeshIO(unittest.TestCase):
     def test_VTUfiles(self):
         mesh = self.mesh
         vec = np.arange(mesh.nC)
-        mesh.writeVTK('temp.vtu',{'arange':vec})
+        mesh.writeVTK('temp.vtu', {'arange': vec})
         print 'Writing of VTU files is working'
         os.remove('temp.vtu')
 
 
-
 if __name__ == '__main__':
     unittest.main()

tests/mesh/test_TreeMesh.py

@@ -6,6 +6,7 @@ import unittest
 
 TOL = 1e-8
 
+
 class TestSimpleQuadTree(unittest.TestCase):
 
     def test_counts(self):

tests/pf/__init__.py

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

tests/pf/assets/magnetics/Mesh_10m.msh

@@ -0,0 +1,5 @@
+47 47 20
+-60.00 -60.00 280.00
+50.00 30.00 20.00 41*10.00 20.00 30.00 50.00 
+50.00 30.00 20.00 41*10.00 20.00 30.00 50.00 
+17*10.00 20.00 30.00 50

tests/pf/assets/magnetics/Obs_loc_TMI.obs

@@ -0,0 +1,619 @@
+90 0  50000
+90 0  1
+616
+           105            155     240.831177    -0.7588417            1
+           115            155     244.046936    -0.7915504            1
+           125            155     247.198578     -0.805861            1
+           135            155     250.256927     -0.798687            1
+           145            155     253.192322    -0.7701055            1
+           155            155     255.975342    -0.7241221            1
+           165            155     258.576965    -0.6681905            1
+           175            155     260.969513    -0.6113366            1
+           185            155     263.126709    -0.5615882            1
+           195            155     265.024445    -0.5237785            1
+           205            155     266.641052    -0.4986475            1
+           215            155     267.957825    -0.4833724            1
+           225            155      268.95929    -0.4730379            1
+           235            155     269.633667    -0.4623284            1
+           245            155       269.9729    -0.4469757            1
+           255            155       269.9729    -0.4247823            1
+           265            155     269.633667      -0.39628            1
+           275            155      268.95929    -0.3650762            1
+           285            155     267.957825    -0.3378048            1
+           295            155     266.641052    -0.3234333            1
+           305            155     265.024445    -0.3316285            1
+           315            155     263.126709    -0.3701661            1
+           325            155     260.969513    -0.4418891            1
+           335            155     258.576965    -0.5423904            1
+           345            155     255.975342     -0.659748            1
+           355            155     253.192322    -0.7770281            1
+           365            155     250.256927    -0.8768104            1
+           375            155     247.198608    -0.9458876            1
+           105            165     242.760101    -0.8195894            1
+           115            165     246.086945    -0.8419799            1
+           125            165     249.347504    -0.8356238            1
+           135            165     252.511536    -0.7963749            1
+           145            165      255.54834    -0.7257527            1
+           155            165      258.42749    -0.6322459            1
+           165            165     261.119019    -0.5301054            1
+           175            165     263.594238    -0.4354313            1
+           185            165     265.825958    -0.3611376            1
+           195            165     267.789246    -0.3131647            1
+           205            165       269.4617    -0.2895989            1
+           215            165     270.823975    -0.2825252            1
+           225            165     271.860046     -0.281212            1
+           235            165     272.557739    -0.2751351            1
+           245            165      272.90863    -0.2561399            1
+           255            165      272.90863    -0.2197196            1
+           265            165     272.557739    -0.1658328            1
+           275            165     271.860046   -0.09951203            1
+           285            165     270.823975   -0.03119349            1
+           295            165       269.4617     0.0236968            1
+           305            165     267.789246    0.04674543            1
+           315            165     265.825958    0.02064035            1
+           325            165     263.594238   -0.06508273            1
+           335            165     261.119019     -0.208611            1
+           345            165      258.42749     -0.393637            1
+           355            165      255.54834    -0.5925498            1
+           365            165     252.511536    -0.7744437            1
+           375            165     249.347534    -0.9143868            1
+           105            175     244.517365    -0.8731715            1
+           115            175     247.945435    -0.8756285            1
+           125            175     251.305176      -0.83312            1
+           135            175      254.56546     -0.739917            1
+           145            175     257.694702    -0.6004066            1
+           155            175     260.661377    -0.4314286            1
+           165            175     263.434845    -0.2595249            1
+           175            175     265.985352    -0.1129547            1
+           185            175     268.284973   -0.01194964            1
+           195            175     270.307983    0.03755872            1
+           205            175     272.031311    0.04395232            1
+           215            175     273.435028     0.0249682            1
+           225            175     274.502655   0.001294248            1
+           235            175     275.221558  -0.008426053            1
+           245            175      275.58313   0.009649532            1
+           255            175      275.58313     0.0634047            1
+           265            175     275.221558      0.154099            1
+           275            175     274.502655     0.2754614            1
+           285            175     273.435028     0.4125105            1
+           295            175     272.031311     0.5409548            1
+           305            175     270.307983     0.6288864            1
+           315            175     268.284973     0.6424912            1
+           325            175     265.985352     0.5560076            1
+           335            175     263.434845     0.3630851            1
+           345            175     260.661377    0.08384954            1
+           355            175     257.694702     -0.237926            1
+           365            175      254.56546     -0.548451            1
+           375            175     251.305176    -0.8018905            1
+           105            185     246.083527    -0.9140503            1
+           115            185     249.601807    -0.8826051            1
+           125            185     253.049957    -0.7811191            1
+           135            185     256.396027    -0.6012919            1
+           145            185     259.607605    -0.3526822            1
+           155            185     262.652374   -0.06704467            1
+           165            185     265.498779     0.2074117            1
+           175            185     268.116394     0.4220746            1
+           185            185     270.476563     0.5463553            1
+           195            185     272.552795     0.5775926            1
+           205            185     274.321503     0.5378997            1
+           215            185     275.762146     0.4625557            1
+           225            185     276.857849     0.3880164            1
+           235            185     277.595642     0.3444396            1
+           245            185     277.966797     0.3530179            1
+           255            185     277.966797     0.4259006            1
+           265            185     277.595642     0.5663478            1
+           275            185     276.857849     0.7675966            1
+           285            185     275.762146      1.010183            1
+           295            185     274.321503      1.259027            1
+           305            185     272.552795      1.463291            1
+           315            185     270.476563      1.563029            1
+           325            185     268.116394      1.504866            1
+           335            185     265.498779      1.263533            1
+           345            185     262.652374     0.8593219            1
+           355            185     259.607605     0.3588987            1
+           365            185     256.396027    -0.1463404            1
+           375            185     253.049988    -0.5744338            1
+           105            195     247.440887    -0.9365266            1
+           115            195     251.037354    -0.8525355            1
+           125            195     254.562134    -0.6607845            1
+           135            195     257.982544    -0.3485086            1
+           145            195     261.265472    0.06623109            1
+           155            195      264.37793     0.5262602            1
+           165            195     267.287598     0.9466072            1
+           175            195     269.963379      1.245725            1
+           185            195     272.375977      1.379536            1
+           195            195     274.498352      1.355977            1
+           205            195     276.306335      1.223446            1
+           215            195     277.778992      1.046257            1
+           225            195     278.899048     0.8838316            1
+           235            195     279.653259     0.7809416            1
+           245            195     280.032623     0.7668838            1
+           255            195     280.032623     0.8582527            1
+           265            195     279.653259      1.060967            1
+           275            195     278.899048      1.369147            1
+           285            195     277.778992      1.760239            1
+           295            195     276.306335      2.187949            1
+           305            195     274.498352      2.577668            1
+           315            195     272.375977      2.831963            1
+           325            195     269.963379      2.852671            1
+           335            195     267.287598      2.577787            1
+           345            195      264.37793       2.01826            1
+           355            195     261.265472       1.27011            1
+           365            195     257.982544     0.4833555            1
+           375            195     254.562134    -0.2003685            1
+           105            205     248.573853    -0.9364405            1
+           115            205     252.235565    -0.7778356            1
+           125            205      255.82428    -0.4572914            1
+           135            205     259.306763     0.0458342            1
+           145            205     262.649292     0.7008426            1
+           155            205     265.818176      1.409817            1
+           165            205      268.78067      2.028682            1
+           175            205     271.505005       2.42519            1
+           185            205     273.961365      2.539679            1
+           195            205     276.122253      2.404351            1
+           205            205     277.963074      2.113149            1
+           215            205     279.462433      1.774462            1
+           225            205     280.602814       1.47808            1
+           235            205     281.370697      1.284939            1
+           245            205     281.756958      1.231368            1
+           255            205     281.756958      1.337238            1
+           265            205     281.370697      1.611238            1
+           275            205     280.602814      2.050138            1
+           285            205     279.462433      2.631077            1
+           295            205     277.963074        3.2984            1
+           305            205     276.122253      3.951254            1
+           315            205     273.961365      4.444233            1
+           325            205     271.505005      4.614739            1
+           335            205      268.78067      4.339799            1
+           345            205     265.818176      3.603591            1
+           355            205     262.649292      2.534426            1
+           365            205     259.306763      1.366716            1
+           375            205      255.82431     0.3338567            1
+           105            215     249.469238    -0.9136291            1
+           115            215     253.182526    -0.6589502            1
+           125            215     256.821808    -0.1698693            1
+           135            215     260.353333     0.5874462            1
+           145            215      263.74292      1.565753            1
+           155            215     266.956482      2.607309            1
+           165            215     269.960632      3.479923            1
+           175            215     272.723358      3.978535            1
+           185            215     275.214325      4.028715            1
+           195            215      277.40564      3.707969            1
+           205            215       279.2724      3.180547            1
+           215            215     280.792847      2.615157            1
+           225            215      281.94931       2.13744            1
+           235            215     282.727997      1.823207            1
+           245            215      283.11969       1.71236            1
+           255            215      283.11969      1.825263            1
+           265            215     282.727997      2.172536            1
+           275            215      281.94931      2.755052            1
+           285            215     280.792847      3.553031            1
+           295            215       279.2724      4.505293            1
+           305            215      277.40564      5.485875            1
+           315            215     275.214325      6.295436            1
+           325            215     272.723358      6.690773            1
+           335            215     269.960632      6.463354            1
+           345            215     266.956482      5.547766            1
+           355            215      263.74292       4.10083            1
+           365            215     260.353333      2.465151            1
+           375            215     256.821838      1.000294            1
+           105            225     250.116516    -0.8743188            1
+           115            225     253.867065    -0.5101706            1
+           125            225     257.542908     0.1755437            1
+           135            225     261.109894      1.236342            1
+           145            225     264.533508      2.606781            1
+           155            225     267.779327      4.051729            1
+           165            225      270.81366      5.220251            1
+           175            225     273.604126      5.813692            1
+           185            225     276.120087      5.748304            1
+           195            225     278.333435      5.170528            1
+           205            225     280.218872      4.338887            1
+           215            225     281.754639      3.495014            1
+           225            225     282.922668      2.801622            1
+           235            225     283.709198      2.344926            1
+           245            225     284.104828       2.16313            1
+           255            225     284.104828      2.273224            1
+           265            225     283.709198      2.685922            1
+           275            225     282.922668       3.40649            1
+           285            225     281.754639      4.420648            1
+           295            225     280.218872      5.665849            1
+           305            225     278.333435      6.995332            1
+           315            225     276.120087      8.157348            1
+           325            225     273.604126       8.82303            1
+           335            225      270.81366      8.684887            1
+           345            225     267.779327      7.610068            1
+           355            225     264.533508      5.774337            1
+           365            225     261.109894      3.638668            1
+           375            225     257.542908      1.709743            1
+           105            235     250.507996    -0.8317488            1
+           115            235     254.281128    -0.3622331            1
+           125            235     257.979034     0.5184553            1
+           135            235     261.567444      1.888687            1
+           145            235     265.011658      3.667592            1
+           155            235     268.276978      5.534542            1
+           165            235      271.32959      7.004851            1
+           175            235      274.13678      7.678418            1
+           185            235     276.667908      7.464011            1
+           195            235     278.894531      6.593227            1
+           205            235     280.791321      5.432199            1
+           215            235     282.336304      4.298311            1
+           225            235     283.511383      3.386243            1
+           235            235     284.302643      2.786252            1
+           245            235     284.700653      2.529866            1
+           255            235     284.700653       2.62728            1
+           265            235     284.302643      3.087025            1
+           275            235     283.511383      3.917431            1
+           285            235     282.336304      5.109731            1
+           295            235     280.791321       6.60235            1
+           305            235     278.894531      8.233635            1
+           315            235     276.667908       9.70905            1
+           325            235      274.13678      10.62564            1
+           335            235      271.32959      10.58569            1
+           345            235     268.276978      9.389627            1
+           355            235     265.011658      7.223036            1
+           365            235     261.567444      4.651613            1
+           375            235     257.979065      2.317601            1
+           105            245     250.639008    -0.8031906            1
+           115            245     254.419678    -0.2567715            1
+           125            245        258.125     0.7720825            1
+           135            245     261.720581      2.386057            1
+           145            245     265.171661      4.496398            1
+           155            245     268.443542      6.708923            1
+           165            245     271.502228      8.428509            1
+           175            245     274.315063      9.166789            1
+           185            245     276.851227      8.817329            1
+           195            245     279.082336      7.692689            1
+           205            245      280.98291      6.255138            1
+           215            245     282.531006      4.884361            1
+           225            245     283.708405      3.796578            1
+           235            245     284.501251        3.0811            1
+           245            245     284.900055      2.761007            1
+           255            245     284.900055      2.838641            1
+           265            245     284.501251      3.318654            1
+           275            245     283.708405      4.209807            1
+           285            245     282.531006      5.506134            1
+           295            245      280.98291      7.146722            1
+           305            245     279.082336      8.961084            1
+           315            245     276.851227      10.62863            1
+           325            245     274.315063      11.70051            1
+           335            245     271.502228      11.72266            1
+           345            245     268.443542      10.45331            1
+           355            245     265.171661      8.084093            1
+           365            245     261.720581       5.24673            1
+           375            245        258.125      2.668323            1
+           105            255     250.507996    -0.8032413            1
+           115            255     254.281128    -0.2308932            1
+           125            255     257.979034     0.8544926            1
+           135            255     261.567444      2.570146            1
+           145            255     265.011658      4.828936            1
+           155            255     268.276978      7.205737            1
+           165            255      271.32959      9.052127            1
+           175            255      274.13678      9.832447            1
+           185            255     276.667908      9.425632            1
+           195            255     278.894531      8.181864            1
+           205            255     280.791321       6.61156            1
+           215            255     282.336304      5.126034            1
+           225            255     283.511383      3.951786            1
+           235            255     284.302643      3.176769            1
+           245            255     284.700653      2.818616            1
+           255            255     284.700653      2.874094            1
+           265            255     284.302643      3.343708            1
+           275            255     283.511383      4.233653            1
+           285            255     282.336304      5.536428            1
+           295            255     280.791321      7.189605            1
+           305            255     278.894531      9.020029            1
+           315            255     276.667908      10.70241            1
+           325            255      274.13678      11.78181            1
+           335            255      271.32959      11.79881            1
+           345            255     268.276978      10.51058            1
+           355            255     265.011658      8.114479            1
+           365            255     261.567444       5.25268            1
+           375            255     257.979065      2.659622            1
+           105            265     250.116516    -0.8364838            1
+           115            265     253.867065    -0.2987034            1
+           125            265     257.542908      0.730639            1
+           135            265     261.109894      2.367063            1
+           145            265     264.533508      4.533685            1
+           155            265     267.779327      6.832248            1
+           165            265      270.81366      8.639468            1
+           175            265     273.604126      9.430596            1
+           185            265     276.120087      9.076918            1
+           195            265     278.333435      7.902471            1
+           205            265     280.218872      6.395443            1
+           215            265     281.754639      4.956879            1
+           225            265     282.922668      3.811898            1
+           235            265     283.709198      3.049655            1
+           245            265     284.104828      2.688498            1
+           255            265     284.104828      2.724281            1
+           265            265     283.709198      3.154537            1
+           275            265     282.922668      3.980581            1
+           285            265     281.754639      5.189228            1
+           295            265     280.218872      6.714512            1
+           305            265     278.333435      8.387216            1
+           315            265     276.120087      9.900153            1
+           325            265     273.604126      10.83422            1
+           335            265      270.81366      10.77838            1
+           345            265     267.779327      9.531084            1
+           355            265     264.533508      7.293797            1
+           365            265     261.109894      4.659332            1
+           375            265     257.542908      2.288398            1
+           105            275     249.469238    -0.8944128            1
+           115            275     253.182526    -0.4425263            1
+           125            275     256.821808     0.4343149            1
+           135            275     260.353333      1.834889            1
+           145            275      263.74292      3.697208            1
+           155            275     266.956482      5.696558            1
+           165            275     269.960632       7.30593            1
+           175            275     272.723358      8.066821            1
+           185            275     275.214325      7.854281            1
+           195            275      277.40564       6.91226            1
+           205            275       279.2724      5.644381            1
+           215            275     280.792847      4.401628            1
+           225            275      281.94931      3.394523            1
+           235            275     282.727997      2.714301            1
+           245            275      283.11969      2.385276            1
+           255            275      283.11969      2.406946            1
+           265            275     282.727997      2.775713            1
+           275            275      281.94931      3.487118            1
+           285            275     280.792847      4.520212            1
+           295            275       279.2724       5.80555            1
+           305            275      277.40564      7.185402            1
+           315            275     275.214325      8.390831            1
+           325            275     272.723358      9.072539            1
+           335            275     269.960632      8.907607            1
+           345            275     266.956482      7.763454            1
+           355            275      263.74292      5.837826            1
+           365            275     260.353333      3.625323            1
+           375            275     256.821838      1.653854            1
+           105            285     248.573853    -0.9591397            1
+           115            285     252.235565    -0.6212285            1
+           125            285      255.82428      0.051057            1
+           135            285     259.306763      1.132654            1
+           145            285     262.649292      2.578346            1
+           155            285     265.818176      4.151248            1
+           165            285      268.78067      5.462206            1
+           175            285     271.505005       6.15344            1
+           185            285     273.961365      6.103346            1
+           195            285     276.122253      5.464153            1
+           205            285     277.963074      4.527912            1
+           215            285     279.462433      3.569985            1
+           225            285     280.602814      2.771323            1
+           235            285     281.370697      2.220876            1
+           245            285     281.756958      1.949892            1
+           255            285     281.756958      1.963643            1
+           265            285     281.370697      2.258935            1
+           275            285     280.602814      2.826297            1
+           285            285     279.462433      3.638017            1
+           295            285     277.963074      4.624463            1
+           305            285     276.122253      5.647346            1
+           315            285     273.961365      6.489984            1
+           325            285     271.505005       6.89047            1
+           335            285      268.78067      6.628548            1
+           345            285     265.818176      5.643015            1
+           355            285     262.649292      4.113403            1
+           365            285     259.306763      2.412046            1
+           375            285      255.82431     0.9149938            1
+           105            295     247.440887     -1.011423            1
+           115            295     251.037354    -0.7896024            1
+           125            295     254.562134    -0.3250723            1
+           135            295     257.982544     0.4339865            1
+           145            295     261.265472      1.458079            1
+           155            295      264.37793      2.590895            1
+           165            295     267.287598      3.574395            1
+           175            295     269.963379      4.158915            1
+           185            295     272.375977      4.235457            1
+           195            295     274.498352      3.878613            1
+           205            295     276.306335       3.27391            1
+           215            295     277.778992      2.615522            1
+           225            295     278.899048      2.044877            1
+           235            295     279.653259      1.641445            1
+           245            295     280.032623      1.439964            1
+           255            295     280.032623      1.450583            1
+           265            295     279.653259      1.671221            1
+           275            295     278.899048      2.089455            1
+           285            295     277.778992      2.674484            1
+           295            295     276.306335       3.36202            1
+           305            295     274.498352      4.039838            1
+           315            295     272.375977      4.548412            1
+           325            295     269.963379      4.711965            1
+           335            295     267.287598      4.401741            1
+           345            295      264.37793       3.60888            1
+           355            295     261.265472      2.480467            1
+           365            295     257.982544      1.271313            1
+           375            295     254.562134     0.2239833            1
+           105            305     246.083527     -1.037815            1
+           115            305     249.601807    -0.9153721            1
+           125            305     253.049957    -0.6277985            1
+           135            305     256.396027    -0.1413393            1
+           145            305     259.607605     0.5270341            1
+           155            305     262.652374      1.283511            1
+           165            305     265.498779      1.970314            1
+           175            305     268.116394      2.427473            1
+           185            305     270.476563       2.56995            1
+           195            305     272.552795      2.421773            1
+           205            305     274.321503      2.085502            1
+           215            305     275.762146      1.684076            1
+           225            305     276.857849      1.317886            1
+           235            305     277.595642      1.051018            1
+           245            305     277.966797     0.9164558            1
+           255            305     277.966797     0.9264848            1
+           265            305     277.595642      1.080114            1
+           275            305     276.857849      1.364457            1
+           285            305     275.762146      1.750175            1
+           295            305     274.321503      2.183568            1
+           305            305     272.552795      2.581238            1
+           315            305     270.476563      2.836312            1
+           325            305     268.116394      2.843435            1
+           335            305     265.498779       2.53956            1
+           345            305     262.652374      1.942758            1
+           355            305     259.607605      1.161767            1
+           365            305     256.396027     0.3578575            1
+           375            305     253.049988    -0.3240296            1
+           105            315     244.517365     -1.033498            1
+           115            315     247.945435    -0.9854098            1
+           125            315     251.305176    -0.8302344            1
+           135            315      254.56546    -0.5464416            1
+           145            315     257.694702    -0.1422417            1
+           155            315     260.661377     0.3302703            1
+           165            315     263.434845     0.7810533            1
+           175            315     265.985352       1.11401            1
+           185            315     268.284973       1.26902            1
+           195            315     270.307983      1.245204            1
+           205            315     272.031311      1.091521            1
+           215            315     273.435028     0.8779288            1
+           225            315     274.502655     0.6689956            1
+           235            315     275.221558     0.5111202            1
+           245            315      275.58313     0.4313801            1
+           255            315      275.58313     0.4413795            1
+           265            315     275.221558     0.5408217            1
+           275            315     274.502655     0.7183567            1
+           285            315     273.435028     0.9495843            1
+           295            315     272.031311      1.194077            1
+           305            315     270.307983      1.395294            1
+           315            315     268.284973      1.488117            1
+           325            315     265.985352      1.416343            1
+           335            315     263.434845      1.155969            1
+           345            315     260.661377     0.7327568            1
+           355            315     257.694702     0.2205752            1
+           365            315      254.56546    -0.2842715            1
+           375            315     251.305176    -0.6998433            1
+           105            325     242.760101     -1.001009            1
+           115            325     246.086945     -1.002393            1
+           125            325     249.347504    -0.9365406            1
+           135            325     252.511536    -0.7895898            1
+           145            325      255.54834    -0.5643467            1
+           155            325      258.42749    -0.2877616            1
+           165            325     261.119019  -0.008496879            1
+           175            325     263.594238     0.2185697            1
+           185            325     265.825958     0.3539491            1
+           195            325     267.789246     0.3879928            1
+           205            325       269.4617     0.3399999            1
+           215            325     270.823975      0.245651            1
+           225            325     271.860046     0.1425535            1
+           235            325     272.557739    0.06089945            1
+           245            325      272.90863    0.02018288            1
+           255            325      272.90863    0.02955809            1
+           265            325     272.557739    0.08900298            1
+           275            325     271.860046     0.1897431            1
+           285            325     270.823975     0.3137744            1
+           295            325       269.4617     0.4335785            1
+           305            325     267.789246     0.5141889            1
+           315            325     265.825958     0.5195781            1
+           325            325     263.594238      0.423505            1
+           335            325     261.119019     0.2213678            1
+           345            325      258.42749   -0.06328931            1
+           355            325      255.54834     -0.383331            1
+           365            325     252.511536    -0.6835498            1
+           375            325     249.347534    -0.9195023            1
+           105            335     240.831177    -0.9470594            1
+           115            335     244.046936    -0.9776354            1
+           125            335     247.198578    -0.9667871            1
+           135            335     250.256927    -0.9053537            1
+           145            335     253.192322    -0.7932306            1
+           155            335     255.975342     -0.643329            1
+           165            335     258.576965    -0.4807753            1
+           175            335     260.969513    -0.3357394            1
+           185            335     263.126709    -0.2327569            1
+           195            335     265.024445    -0.1823665            1
+           205            335     266.641052    -0.1795066            1
+           215            335     267.957825    -0.2082184            1
+           225            335      268.95929    -0.2487825            1
+           235            335     269.633667    -0.2835538            1
+           245            335       269.9729     -0.300115            1
+           255            335       269.9729    -0.2922643            1
+           265            335     269.633667    -0.2600168            1
+           275            335      268.95929     -0.209445            1
+           285            335     267.957825    -0.1524777            1
+           295            335     266.641052    -0.1060946            1
+           305            335     265.024445   -0.08997072            1
+           315            335     263.126709     -0.121989            1
+           325            335     260.969513     -0.212178            1
+           335            335     258.576965    -0.3573267            1
+           345            335     255.975342    -0.5393556            1
+           355            335     253.192322    -0.7294713            1
+           365            335     250.256927    -0.8969562            1
+           375            335     247.198608     -1.018534            1
+           105            345     238.751434    -0.8795042            1
+           115            345     241.847382    -0.9247606            1
+           125            345     244.881653    -0.9446973            1
+           135            345      247.82608    -0.9330818            1
+           145            345      250.65213    -0.8885093            1
+           155            345     253.331451    -0.8165421            1
+           165            345     255.836212    -0.7295913            1
+           175            345     258.139618    -0.6437526            1
+           185            345     260.216461    -0.5736571            1
+           195            345     262.043518    -0.5278839            1
+           205            345     263.599884    -0.5071621            1
+           215            345     264.867584    -0.5057246            1
+           225            345     265.831787    -0.5144912            1
+           235            345     266.481018    -0.5243066            1
+           245            345     266.807617    -0.5281933            1
+           255            345     266.807617    -0.5225207            1
+           265            345     266.481018    -0.5073952            1
+           275            345     265.831787    -0.4866288            1
+           285            345     264.867584    -0.4673505            1
+           295            345     263.599884    -0.4590167            1
+           305            345     262.043518    -0.4715563            1
+           315            345     260.216461    -0.5125949            1
+           325            345     258.139618    -0.5843691            1
+           335            345     255.836212    -0.6815979            1
+           345            345     253.331451    -0.7916309            1
+           355            345      250.65213    -0.8974017            1
+           365            345      247.82608    -0.9822296            1
+           375            345     244.881653     -1.034367            1
+           105            355     236.542786    -0.8053916            1
+           115            355     239.511536    -0.8559309            1
+           125            355     242.421112    -0.8909921            1
+           135            355     245.244568    -0.9062259            1
+           145            355     247.954498    -0.8998395            1
+           155            355     250.523712    -0.8737812            1
+           165            355     252.925537    -0.8338691            1
+           175            355     255.134338    -0.7884543            1
+           185            355     257.125824    -0.7460073            1
+           195            355     258.877808     -0.712674            1
+           205            355     260.370239    -0.6909303            1
+           215            355     261.585815    -0.6797073            1
+           225            355     262.510406    -0.6756343            1
+           235            355     263.132996    -0.6746618            1
+           245            355     263.446106    -0.6734288            1
+           255            355     263.446106    -0.6701413            1
+           265            355     263.132996    -0.6649869            1
+           275            355     262.510406    -0.6601467            1
+           285            355     261.585815    -0.6594277            1
+           295            355     260.370239     -0.667488            1
+           305            355     258.877808    -0.6886132            1
+           315            355     257.125824    -0.7251665            1
+           325            355     255.134338    -0.7761745            1
+           335            355     252.925537    -0.8366687            1
+           345            355     250.523712    -0.8982714            1
+           355            355     247.954498    -0.9510985            1
+           365            355     245.244568     -0.986296            1
+           375            355     242.421143    -0.9982168            1
+           105            365     234.227783    -0.7301244            1
+           115            365     237.063202    -0.7803932            1
+           125            365     239.842102    -0.8211433            1
+           135            365     242.538727    -0.8492868            1
+           145            365     245.126953    -0.8630885            1
+           155            365      247.58078    -0.8628479            1
+           165            365     249.874756    -0.8510816            1
+           175            365     251.984314    -0.8319991            1
+           185            365     253.886353    -0.8103872            1
+           195            365     255.559631    -0.7903531            1
+           205            365     256.985046    -0.7744435            1
+           215            365     258.146057    -0.7634233            1
+           225            365     259.029114     -0.756661            1
+           235            365     259.623718    -0.7528488            1
+           245            365     259.922791    -0.7507278            1
+           255            365     259.922791    -0.7496345            1
+           265            365     259.623718    -0.7497882            1
+           275            365     259.029114    -0.7523054            1
+           285            365     258.146057    -0.7589456            1
+           295            365     256.985046    -0.7716019            1
+           305            365     255.559631    -0.7915938            1
+           315            365     253.886353    -0.8188929            1
+           325            365     251.984314    -0.8515587            1
+           335            365     249.874756    -0.8856623            1
+           345            365      247.58078    -0.9158597            1
+           355            365     245.126953    -0.9365292            1
+           365            365     242.538727    -0.9430979            1
+           375            365     239.842102    -0.9330602            1

tests/pf/assets/magnetics/PYMAG3D_inv.inp

@@ -0,0 +1,13 @@
+Mesh_10m.msh			! Mesh file
+Obs_loc_TMI.obs			! Obsfile
+Gaussian.topo			! Topofile | null
+ModelStart.sus			! Starting model
+VALUE 0				! Reference model
+ActiveCells.dat			! [ Active file (0:Inactive | 1: Active-Dynamic | -1: Active-Static) ] || [ VALUE ### (from starting model) ] || [ DEFAULT (all Active-Dynamic)]
+DEFAULT !..\AzmDip.dat		! Magnetization vector model
+DEFAULT				! Cell based weight file
+1				! target chi factor | DEFAULT=1
+1 1 1 1				! alpha s, x ,y ,z
+VALUE 0 1			! Lower and Upper Bounds for p-component
+VALUE 0 1 1 1 1			! lp-norm for amplitude inversion FILE pqxqyqzr.dat  ! Norms VALUE p, qx, qy, qz, r | FILE m-by-5 matrix
+DEFAULT		! Threshold value for the norm on model and model gradient VALUE eps_p, eps_q | DEFAULT

tests/pf/test_forward_PFproblem.py

@@ -0,0 +1,56 @@
+import unittest
+from SimPEG import Mesh, Utils, np, PF
+import matplotlib.pyplot as plt
+
+
+class MagFwdProblemTests(unittest.TestCase):
+
+    def setUp(self):
+
+        cs = 25.
+        hxind = [(cs, 5, -1.3), (cs/2.0, 41), (cs, 5, 1.3)]
+        hyind = [(cs, 5, -1.3), (cs/2.0, 41), (cs, 5, 1.3)]
+        hzind = [(cs, 5, -1.3), (cs/2.0, 40), (cs, 5, 1.3)]
+        M = Mesh.TensorMesh([hxind, hyind, hzind], 'CCC')
+
+        chibkg = 0.
+        chiblk = 0.01
+        chi = np.ones(M.nC)*chibkg
+        sph_ind = PF.MagAnalytics.spheremodel(M, 0., 0., 0., 100)
+        chi[sph_ind] = chiblk
+        model = PF.BaseMag.BaseMagMap(M)
+        prob = PF.Magnetics.Problem3D_DiffSecondary(M, mapping=model)
+        self.prob = prob
+        self.M = M
+        self.chi = chi
+
+    def test_ana_forward(self):
+
+        survey = PF.BaseMag.BaseMagSurvey()
+
+        Inc = 45.
+        Dec = 45.
+        Btot = 51000
+
+        b0 = PF.MagAnalytics.IDTtoxyz(Inc, Dec, Btot)
+        survey.setBackgroundField(Inc, Dec, Btot)
+        xr = np.linspace(-300, 300, 41)
+        yr = np.linspace(-300, 300, 41)
+        X, Y = np.meshgrid(xr, yr)
+        Z = np.ones((xr.size, yr.size))*150
+        rxLoc = np.c_[Utils.mkvc(X), Utils.mkvc(Y), Utils.mkvc(Z)]
+        survey.rxLoc = rxLoc
+
+        self.prob.pair(survey)
+        u = self.prob.fields(self.chi)
+        B = u['B']
+
+        bxa, bya, bza = PF.MagAnalytics.MagSphereAnaFunA(rxLoc[:, 0], rxLoc[:, 1], rxLoc[:, 2], 100., 0., 0., 0., 0.01, b0, 'secondary')
+
+        dpred = survey.projectFieldsAsVector(B)
+        err = np.linalg.norm(dpred-np.r_[bxa, bya, bza])/np.linalg.norm(np.r_[bxa, bya, bza])
+        self.assertTrue(err < 0.08)
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/pf/test_magnetics_IO.py

@@ -0,0 +1,36 @@
+import unittest
+from SimPEG import Mesh, np, PF
+from scipy.constants import mu_0
+import os
+
+
+class MagSensProblemTests(unittest.TestCase):
+
+    def setUp(self):
+        pass
+
+    def test_magnetics_inversion(self):
+
+        driver = PF.MagneticsDriver.MagneticsDriver_Inv(
+            os.path.sep.join([
+                os.path.split(__file__)[0],
+                'assets', 'magnetics', 'PYMAG3D_inv.inp'
+            ])
+        )
+
+        print driver.mesh
+        print driver.survey
+        print driver.m0
+        print driver.mref
+        print driver.activeCells
+        print driver.staticCells
+        print driver.dynamicCells
+        print driver.chi
+        print driver.nC
+        print driver.alphas
+        print driver.bounds
+        print driver.lpnorms
+        print driver.eps
+
+if __name__ == '__main__':
+    unittest.main()

tests/pf/test_magnetics_analytics.py

@@ -0,0 +1,48 @@
+import unittest
+from SimPEG import Mesh, np, PF
+from scipy.constants import mu_0
+
+
+plotIt = False
+
+
+class TestBoundaryConditionAnalytics(unittest.TestCase):
+
+    def test_ana_boundary_computation(self):
+
+        hxind = [(0, 25, 1.3), (21, 12.5), (0, 25, 1.3)]
+        hyind = [(0, 25, 1.3), (21, 12.5), (0, 25, 1.3)]
+        hzind = [(0, 25, 1.3), (20, 12.5), (0, 25, 1.3)]
+        # hx, hy, hz = Utils.meshTensors(hxind, hyind, hzind)
+        M3 = Mesh.TensorMesh([hxind, hyind, hzind], "CCC")
+        indxd, indxu, indyd, indyu, indzd, indzu = M3.faceBoundaryInd
+        mu0 = 4*np.pi*1e-7
+        chibkg = 0.
+        chiblk = 0.01
+        chi = np.ones(M3.nC)*chibkg
+        sph_ind = PF.MagAnalytics.spheremodel(M3, 0, 0, 0, 100)
+        chi[sph_ind] = chiblk
+        mu = (1.+chi)*mu0
+        Bbc, const = PF.MagAnalytics.CongruousMagBC(M3, np.array([1., 0., 0.]), chi)
+
+        flag = 'secondary'
+        Box = 1.
+        H0 = Box/mu_0
+        Bbcxx, Bbcxy, Bbcxz  = PF.MagAnalytics.MagSphereAnaFun(M3.gridFx[(indxd|indxu),0], M3.gridFx[(indxd|indxu),1], M3.gridFx[(indxd|indxu),2], 100, 0., 0., 0., mu_0, mu_0*(1+chiblk), H0, flag)
+        Bbcyx, Bbcyy, Bbcyz  = PF.MagAnalytics.MagSphereAnaFun(M3.gridFy[(indyd|indyu),0], M3.gridFy[(indyd|indyu),1], M3.gridFy[(indyd|indyu),2], 100, 0., 0., 0., mu_0, mu_0*(1+chiblk), H0, flag)
+        Bbczx, Bbczy, Bbczz  = PF.MagAnalytics.MagSphereAnaFun(M3.gridFz[(indzd|indzu),0], M3.gridFz[(indzd|indzu),1], M3.gridFz[(indzd|indzu),2], 100, 0., 0., 0., mu_0, mu_0*(1+chiblk), H0, flag)
+        Bbc_ana = np.r_[Bbcxx, Bbcyy, Bbczz]
+
+        if plotIt:
+            import matplotlib.pyplot as plt
+            fig, ax = plt.subplots(1,1, figsize = (10, 10))
+            ax.plot(Bbc_ana)
+            ax.plot(Bbc)
+            plt.show()
+        err = np.linalg.norm(Bbc-Bbc_ana) / np.linalg.norm(Bbc_ana)
+
+        assert err < 0.1, 'Mag Boundary computation is wrong!!, err = {}'.format(err)
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/pf/test_sensitivity_PFproblem.py

@@ -0,0 +1,306 @@
+# import unittest
+# from SimPEG import *
+# from simpegPF import BaseMag
+# import matplotlib.pyplot as plt
+# import simpegPF as PF
+# from scipy.constants import mu_0
+
+
+# class MagSensProblemTests(unittest.TestCase):
+
+#     def setUp(self):
+#         cs = 25.
+#         hxind = [(cs, 5, -1.3), (cs, 21), (cs, 5, 1.3)]
+#         hyind = [(cs, 5, -1.3), (cs, 21), (cs, 5, 1.3)]
+#         hzind = [(cs, 5, -1.3), (cs, 20), (cs, 5, 1.3)]
+#         M = Mesh.TensorMesh([hxind, hyind, hzind], 'CCC')
+#         chibkg = 0.001
+#         chiblk = 0.01
+#         chi = np.ones(M.nC)*chibkg
+
+#         Inc = 90.
+#         Dec = 0.
+#         Btot = 51000
+
+#         b0 = PF.MagAnalytics.IDTtoxyz(Inc, Dec, Btot)
+#         sph_ind = PF.MagAnalytics.spheremodel(M, 0., 0., 0., 100)
+#         chi[sph_ind] = chiblk
+#         model = PF.BaseMag.BaseMagMap(M)
+
+#         survey = BaseMag.BaseMagSurvey()
+
+#         survey.setBackgroundField(Inc, Dec, Btot)
+#         xr = np.linspace(-300, 300, 41)
+#         yr = np.linspace(-300, 300, 41)
+#         X, Y = np.meshgrid(xr, yr)
+#         Z = np.ones((xr.size, yr.size))*150
+#         rxLoc = np.c_[Utils.mkvc(X), Utils.mkvc(Y), Utils.mkvc(Z)]
+#         survey.rxLoc = rxLoc
+
+#         prob = PF.Magnetics.Problem3D_DiffSecondary(M, mapping=model)
+#         prob.pair(survey)
+#         dpre = survey.dpred(chi)
+
+#         fields = prob.fields(chi)
+#         self.u = fields['u']
+#         self.B = fields['B']
+
+#         self.survey = survey
+#         self.model = model
+#         self.prob = prob
+#         self.M = M
+#         self.chi = chi
+
+#     def test_mass(self):
+#         print '\n >>Derivative for MfMuI works.'
+#         mu = self.model*self.chi
+
+#         def MfmuI(mu):
+
+#             chi = mu/mu_0-1
+#             self.prob.makeMassMatrices(chi)
+#             vol = self.prob.mesh.vol
+#             aveF2CC = self.prob.mesh.aveF2CC
+#             MfMuI = self.prob.MfMuI.diagonal()
+
+#             return MfMuI
+
+#         def dMfmuI(mu, v):
+
+#             chi = mu/mu_0-1
+#             self.prob.makeMassMatrices(chi)
+#             vol = self.prob.mesh.vol
+#             aveF2CC = self.prob.mesh.aveF2CC
+#             MfMuI = self.prob.MfMuI.diagonal()
+#             dMfMuI = Utils.sdiag(MfMuI**2)*aveF2CC.T*Utils.sdiag(vol*1./mu**2)
+
+#             return dMfMuI*v
+
+#         d_mu = mu*0.8
+#         derChk = lambda m: [MfmuI(m), lambda mx: dMfmuI(self.chi, mx)]
+#         passed = Tests.checkDerivative(derChk, mu, num=4, dx = d_mu, plotIt=False)
+
+#         self.assertTrue(passed)
+
+
+#     def test_dCdm_Av(self):
+#         print '\n >>Derivative for Cm_A.'
+#         Div = self.prob._Div
+#         vol = self.prob.mesh.vol
+#         aveF2CC = self.prob.mesh.aveF2CC
+
+#         def Cm_A(chi):
+#             dmudm = self.model.deriv(chi)
+#             u = self.u
+#             # chi = mu/mu_0-1
+#             self.prob.makeMassMatrices(chi)
+#             mu = self.model*(self.chi)
+#             A = self.prob.getA(self.chi)
+#             MfMuIvec = 1/self.prob.MfMui.diagonal()
+#             dMfMuI = Utils.sdiag(MfMuIvec**2)*aveF2CC.T*Utils.sdiag(vol*1./mu**2)
+
+#             Cm_A =  A*u
+
+#             return Cm_A
+
+#         def dCdm_A(chi, v):
+
+#             dmudm = self.model.deriv(chi)
+#             u = self.u
+#             self.prob.makeMassMatrices(chi)
+#             mu = self.model*self.chi
+#             A = self.prob.getA(self.chi)
+#             MfMuIvec = 1/self.prob.MfMui.diagonal()
+#             dMfMuI = Utils.sdiag(MfMuIvec**2)*aveF2CC.T*Utils.sdiag(vol*1./mu**2)
+
+#             Cm_A =  A*u
+#             dCdm_A = Div * (Utils.sdiag( Div.T * u ) * dMfMuI * dmudm)
+
+#             return dCdm_A*v
+
+#         d_chi = self.chi*0.8
+#         derChk = lambda m: [Cm_A(m), lambda mx: dCdm_A(self.chi, mx)]
+#         passed = Tests.checkDerivative(derChk, self.chi, num=4, dx = d_chi, plotIt=False)
+#         self.assertTrue(passed)
+
+
+#     def test_dCdmu_RHS(self):
+#         print '\n >>Derivative for Cm_RHS.'
+#         u = self.u
+#         Div = self.prob._Div
+#         mu = self.model*self.chi
+#         vol = self.prob.mesh.vol
+#         Mc = Utils.sdiag(vol)
+#         aveF2CC = self.prob.mesh.aveF2CC
+#         B0 = self.prob.getB0()
+#         Dface = self.prob.mesh.faceDiv
+
+#         def Cm_RHS(chi):
+
+#             self.prob.makeMassMatrices(chi)
+#             dmudm = self.model.deriv(chi)
+#             dchidmu = Utils.sdiag(1/(dmudm.diagonal()))
+#             Bbc, Bbc_const = PF.MagAnalytics.CongruousMagBC(self.prob.mesh, self.survey.B0, chi)
+#             MfMuIvec = 1/self.prob.MfMui.diagonal()
+#             dMfMuI = Utils.sdiag(MfMuIvec**2)*aveF2CC.T*Utils.sdiag(vol*1./mu**2)
+#             RHS1 = Div*self.prob.MfMuI*self.prob.MfMu0*B0
+#             RHS2 =  Mc*Dface*self.prob._Pout.T*Bbc
+#             RHS = RHS1 + RHS2 + Div*B0
+
+#             return RHS
+
+
+#         def dCdm_RHS(chi, v):
+
+
+#             self.prob.makeMassMatrices(chi)
+#             dmudm = self.model.deriv(chi)
+#             dmdmu = Utils.sdiag(1/(dmudm.diagonal()))
+#             Bbc, Bbc_const = PF.MagAnalytics.CongruousMagBC(self.prob.mesh, self.survey.B0, chi)
+#             MfMuIvec = 1/self.prob.MfMui.diagonal()
+#             dMfMuI = Utils.sdiag(MfMuIvec**2)*aveF2CC.T*Utils.sdiag(vol*1./mu**2)
+#             dCdm_RHS1 = Div * (Utils.sdiag( self.prob.MfMu0*B0  ) * dMfMuI)
+#             temp1 = (Dface*(self.prob._Pout.T*Bbc_const*Bbc))
+#             dCdm_RHS2v  = (Utils.sdiag(vol)*temp1)*np.inner(vol, v)
+#             dCdm_RHSv =  dCdm_RHS1*(dmudm*v) + dCdm_RHS2v
+
+#             return dCdm_RHSv
+
+#         d_chi = self.chi*0.8
+#         derChk = lambda m: [Cm_RHS(m), lambda mx: dCdm_RHS(self.chi, mx)]
+#         passed = Tests.checkDerivative(derChk, self.chi, num=4, dx = d_chi, plotIt=False)
+#         self.assertTrue(passed)
+
+
+#     # def test_dudm(self):
+#     #     print ">> Derivative test for dudm"
+#     #     u = self.u
+#     #     Div = self.prob._Div
+#     #     mu = self.model*(self.chi)
+#     #     vol = self.prob.mesh.vol
+#     #     Mc = Utils.sdiag(vol)
+#     #     aveF2CC = self.prob.mesh.aveF2CC
+#     #     B0 = self.prob.getB0()
+#     #     Dface = self.prob.mesh.faceDiv
+
+#     #     def ufun(chi):
+#     #         u = self.prob.fields(chi)['u']
+#     #         return u
+
+#     #     def dudm(chi, v):
+
+#     #         chi = mu/mu_0-1
+#     #         self.prob.makeMassMatrices(chi)
+#     #         u = self.u
+#     #         dmudm = self.model.deriv(chi)
+#     #         dmdmu = Utils.sdiag(1/(dmudm.diagonal()))
+#     #         Bbc, Bbc_const = PF.MagAnalytics.CongruousMagBC(self.prob.mesh, self.survey.B0, chi)
+#     #         MfMuIvec = 1/self.prob.MfMui.diagonal()
+#     #         dMfMuI = Utils.sdiag(MfMuIvec**2)*aveF2CC.T*Utils.sdiag(vol*1./mu**2)
+#     #         dCdu = self.prob.getA(chi)
+#     #         dCdm_A = Div * ( Utils.sdiag( Div.T * u )* dMfMuI *dmudm  )
+#     #         dCdm_RHS1 = Div * (Utils.sdiag( self.prob.MfMu0*B0  ) * dMfMuI)
+#     #         temp1 = (Dface*(self.prob._Pout.T*Bbc_const*Bbc))
+#     #         dCdm_RHS2v  = (Utils.sdiag(vol)*temp1)*np.inner(vol, v)
+#     #         dCdm_RHSv =  dCdm_RHS1*(dmudm*v) +  dCdm_RHS2v
+#     #         dCdm_v = dCdm_A*v - dCdm_RHSv
+#     #         m1 = sp.linalg.interface.aslinearoperator(Utils.sdiag(1/dCdu.diagonal()))
+#     #         sol, info = sp.linalg.bicgstab(dCdu, dCdm_v, tol=1e-8, maxiter=1000, M=m1)
+
+#     #         dudm = -sol
+
+#     #         return dudm
+
+#     #     d_chi = 10.0*self.chi #np.random.rand(mesh.nCz)
+#     #     d_sph_ind = PF.MagAnalytics.spheremodel(self.prob.mesh, 0., 0., -50., 50)
+#     #     d_chi[d_sph_ind] = 0.1
+
+#     #     derChk = lambda m: [ufun(m), lambda mx: dudm(self.chi, mx)]
+#     #     # TODO: I am not sure why the order get worse as step decreases .. --;
+#     #     passed = Tests.checkDerivative(derChk, self.chi, num=2, dx = d_chi, plotIt=False)
+#     #     self.assertTrue(passed)
+
+
+#     # def test_dBdm(self):
+#     #     print ">> Derivative test for dBdm"
+#     #     u = self.u
+#     #     Div = self.prob._Div
+#     #     mu = self.model*(self.chi)
+#     #     vol = self.prob.mesh.vol
+#     #     Mc = Utils.sdiag(vol)
+#     #     aveF2CC = self.prob.mesh.aveF2CC
+#     #     B0 = self.prob.getB0()
+#     #     Dface = self.prob.mesh.faceDiv
+
+#     #     def Bfun(chi):
+#     #         B = self.prob.fields(chi)['B']
+#     #         return B
+
+#     #     def dBdm(chi, v):
+
+#     #         chi = mu/mu_0-1
+#     #         self.prob.makeMassMatrices(chi)
+#     #         u = self.u
+#     #         dmudm = self.model.deriv(chi)
+#     #         dmdmu = Utils.sdiag(1/(dmudm.diagonal()))
+#     #         Bbc, Bbc_const = PF.MagAnalytics.CongruousMagBC(self.prob.mesh, self.survey.B0, chi)
+#     #         MfMuIvec = 1/self.prob.MfMui.diagonal()
+#     #         dMfMuI = Utils.sdiag(MfMuIvec**2)*aveF2CC.T*Utils.sdiag(vol*1./mu**2)
+#     #         dCdu = self.prob.getA(chi)
+#     #         dCdm_A = Div * ( Utils.sdiag( Div.T * u )* dMfMuI *dmudm  )
+#     #         dCdm_RHS1 = Div * (Utils.sdiag( self.prob.MfMu0*B0  ) * dMfMuI)
+#     #         temp1 = (Dface*(self.prob._Pout.T*Bbc_const*Bbc))
+#     #         dCdm_RHS2v  = (Utils.sdiag(vol)*temp1)*np.inner(vol, v)
+#     #         dCdm_RHSv =  dCdm_RHS1*(dmudm*v) +  dCdm_RHS2v
+#     #         dCdm_v = dCdm_A*v - dCdm_RHSv
+#     #         m1 = sp.linalg.interface.aslinearoperator(Utils.sdiag(1/dCdu.diagonal()))
+#     #         sol, info = sp.linalg.bicgstab(dCdu, dCdm_v, tol=1e-8, maxiter=1000, M=m1)
+
+#     #         dudm = -sol
+#     #         dBdmv =       (  Utils.sdiag(self.prob.MfMu0*B0)*(dMfMuI * (dmudm*v)) \
+#     #                      - Utils.sdiag(Div.T*u)*(dMfMuI* (dmudm*v)) \
+#     #                      - self.prob.MfMuI*(Div.T* (dudm)) )
+
+#     #         return dBdmv
+
+#     #     d_chi = 10.0*self.chi #np.random.rand(mesh.nCz)
+#     #     d_sph_ind = PF.MagAnalytics.spheremodel(self.prob.mesh, 0., 0., -50., 50)
+#     #     d_chi[d_sph_ind] = 0.1
+
+#     #     derChk = lambda m: [Bfun(m), lambda mx: dBdm(self.chi, mx)]
+#     #     # TODO: I am not sure why the order get worse as step decreases .. --;
+#     #     passed = Tests.checkDerivative(derChk, self.chi, num=2, dx = d_chi, plotIt=False)
+#     #     self.assertTrue(passed)
+
+
+
+#     def test_Jvec(self):
+#         print ">> Derivative test for Jvec"
+
+#         d_chi = 10.0*self.chi #np.random.rand(mesh.nCz)
+#         d_sph_ind = PF.MagAnalytics.spheremodel(self.prob.mesh, 0., 0., -50., 50)
+#         d_chi[d_sph_ind] = 0.1
+
+#         derChk = lambda m: (self.survey.dpred(m), lambda v: self.prob.Jvec(m, v))
+#         # TODO: I am not sure why the order get worse as step decreases .. --;
+#         passed = Tests.checkDerivative(derChk, self.chi, num=2, dx = d_chi, plotIt=False)
+#         self.assertTrue(passed)
+
+#     def test_Jtvec(self):
+#         print ">> Derivative test for Jtvec"
+#         dobs = self.survey.dpred(self.chi)
+
+#         def misfit(m):
+#             dpre = self.survey.dpred(m)
+#             misfit = 0.5*np.linalg.norm(dpre-dobs)**2
+#             residual = dpre-dobs
+#             dmisfit = self.prob.Jtvec(self.chi, residual)
+
+#             return misfit, dmisfit
+
+#         # TODO: I am not sure why the order get worse as step decreases .. --;
+#         passed = Tests.checkDerivative(misfit, self.chi, num=4, plotIt=False)
+#         self.assertTrue(passed)
+
+# if __name__ == '__main__':
+#     unittest.main()