mag work, forward problem running

simpegPF/BaseMag.py

@@ -9,20 +9,69 @@ class BaseMagData(Data.BaseData):
     rxType = None #: receiver type
 
     def __init__(self, **kwargs):
-        Data.BaseData.__init__(**kwargs)
-
+        Data.BaseData.__init__(self, **kwargs)
 
     #TODO: change to inc, dec, intensity
-    def setBackgroundField(self, x, y, z):
+    def setBackgroundField(self, x=1., y=0., z=0.):
         # Primary field in x-direction (background)
-        self.B0 = np.r_[1.,0,0]
+        self.B0 = np.r_[x,y,z]
+
+    @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, B):
+        """
+            This function projects the fields onto the data space.
+
+
+            .. math::
+                d_\\text{pred} = \mathbf{P} u(m)
+        """
+
+        bfx = self.Qfx*B
+        bfy = self.Qfy*B
+        bfz = self.Qfz*B
+
+        return np.sqrt(bfx**2 + bfy**2 + bfz**2)
+
+    def projectFieldsAsVector(self, B):
+
+        bfx = self.Qfx*B
+        bfy = self.Qfy*B
+        bfz = self.Qfz*B
+
+        return np.c_[bfx, bfy, bfz]
+
+class MagDataBx(object):
+    """docstring for MagDataBx"""
+    def __init__(self, **kwargs):
+        Data.BaseData.__init__(self, **kwargs)
+
+    def projectFields(self, B):
+        bfx = self.Qfx*B
+        return bfx
 
 
 class BaseMagModel(Model.BaseModel):
     """BaseMagModel"""
 
     def __init__(self, mesh, **kwargs):
-        Model.BaseModel.__init__(mesh, **kwargs)
+        Model.BaseModel.__init__(self, mesh)
 
     def transform(self, m, asMu=True):
         if asMu:

simpegPF/MagAnalytics.py

@@ -1,5 +1,5 @@
 from scipy.constants import mu_0
-from SimPEG.Utils.sputils import kron3, speye, sdiag
+from SimPEG.Utils.matutils import kron3, speye, sdiag
 import matplotlib.pyplot as plt
 from SimPEG import *
 import numpy as np

simpegPF/Magnetics.py

@@ -1,6 +1,7 @@
-from SimPEG import Problem, Utils
+from SimPEG import Mesh, Problem, Utils, np, sp
 import BaseMag
 from scipy.constants import mu_0
+from MagAnalytics import spheremodel, CongruousMagBC
 
 
 
@@ -11,27 +12,27 @@ class MagneticsDiffSecondary(Problem.BaseProblem):
     modelPair = BaseMag.BaseMagModel
 
     def __init__(self, mesh, model, **kwargs):
-        Problem.BaseProblem.__init__(mesh, model, **kwargs)
+        Problem.BaseProblem.__init__(self, mesh, model, **kwargs)
 
-        self._Pbc, self._Pin, self._Pout = /
+        Pbc, Pin, self._Pout = \
             self.mesh.getBCProjWF('neumann', discretization='CC')
 
         Dface = self.mesh.faceDiv
         Mc = Utils.sdiag(self.mesh.vol)
-        self._Div = Mc*Dface*self._Pin.T*self._Pin
+        self._Div = Mc*Dface*Pin.T*Pin
 
     @property
     def MfMuI(self): return self._MfMuI
 
     @property
-    def Mfmu0(self): return self._Mfmu0
+    def MfMu0(self): return self._MfMu0
 
     def makeMassMatrices(self, m):
         mu = self.model.transform(m)
-        Mfmui = self.mesh.getFaceInnerProduct(1./mu)
+        MfMui = self.mesh.getFaceInnerProduct(1./mu)
         #TODO: this will break if tensor mu
-        self._MfmuI = Utils.sdiag(1./Mfmui.diagonal())
-        self._Mfmu0 = self.mesh.getFaceInnerProduct(1/mu_0)
+        self._MfMuI = Utils.sdiag(1./MfMui.diagonal())
+        self._MfMu0 = self.mesh.getFaceInnerProduct(1/mu_0)
 
     def getRHS(self, m):
         b0 = self.data.B0
@@ -45,7 +46,7 @@ class MagneticsDiffSecondary(Problem.BaseProblem):
         chi = self.model.transform(m, asMu=False)
         Bbc = CongruousMagBC(self.mesh, self.data.B0, chi)
 
-        return -self._Div*self.MfmuI*self.Mfmu0*B0 + self._Div*B0 - Mc*Dface*self._Pout.T*Bbc
+        return -self._Div*self.MfMuI*self.MfMu0*B0 + self._Div*B0 - Mc*Dface*self._Pout.T*Bbc
 
     def getA(self, m):
         """
@@ -60,15 +61,74 @@ class MagneticsDiffSecondary(Problem.BaseProblem):
 
 
         """
-        return -self._Div*self.MfmuI*self._Div.T
+        return -self._Div*self.MfMuI*self._Div.T
 
 
     def fields(self, m):
         self.makeMassMatrices(m)
+        #TODO: change to pos def A
         A = self.getA(m)
         rhs = self.getRHS(m)
 
+        m1 = sp.linalg.interface.aslinearoperator(Utils.sdiag(-1/A.diagonal()))
+        phi, info = sp.linalg.bicgstab(A, rhs, tol=1e-6, maxiter=1000, M=m1)
+
+        #TODO: make onPair function call
+        b0 = self.data.B0
+        B0 = np.r_[b0[0]*np.ones(self.mesh.nFx),
+                   b0[1]*np.ones(self.mesh.nFy),
+                   b0[2]*np.ones(self.mesh.nFz)]
+
+        B = self.MfMuI*self.MfMu0*B0-B0-self.MfMuI*self._Div.T*phi
+
+        return B
+
         # F = self.getInitialFields()
         # return self.forward(m, self.getRHS, self.calcFields, F=F)
 
 
+if __name__ == '__main__':
+    import matplotlib.pyplot as plt
+    hxind = ((5,25,1.3),(41, 12.5),(5,25,1.3))
+    hyind = ((5,25,1.3),(41, 12.5),(5,25,1.3))
+    hzind = ((5,25,1.3),(40, 12.5),(1,25,1.3))
+    hx, hy, hz = Utils.meshTensors(hxind, hyind, hzind)
+    mesh = Mesh.TensorMesh([hx, hy, hz], [-hx.sum()/2,-hy.sum()/2,-hz.sum()/2])
+
+    chibkg = 0.
+    chiblk = 0.01
+    chi = np.ones(mesh.nC)*chibkg
+    sph_ind = spheremodel(mesh, 0., 0., 0., 100)
+    chi[sph_ind] = chiblk
+    model = BaseMag.BaseMagModel(mesh)
+    # mu = (1.+chi)*mu_0
+
+    data = BaseMag.BaseMagData()
+    data.setBackgroundField(x=1., y=1., z=0.)
+    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)]
+    data.rxLoc = rxLoc
+
+    prob = MagneticsDiffSecondary(mesh, model)
+
+    prob.pair(data)
+
+    B = prob.fields(chi)
+    mesh.plotSlice(B, 'F', view='vec', showIt=True)
+
+    dpred = data.dpred(chi, u=B)
+
+    plt.pcolor(X, Y, dpred.reshape(X.shape, order='F'))
+
+    plt.show()
+
+
+
+
+
+
+
+