Analytic test

simpegPF/MagAnalytics.py

@@ -3,6 +3,7 @@ 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
@@ -13,6 +14,7 @@ def spheremodel(mesh, x0, y0, z0, r):
     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
@@ -29,7 +31,6 @@ def MagSphereAnaFun(x, y, z, R, x0, y0, z0, mu1, mu2, H0, flag='total'):
 
 
     """
-    print H0
 
     if (~np.size(x)==np.size(y)==np.size(z)):
         print "Specify same size of x, y, z"
@@ -180,6 +181,7 @@ def MagSphereAnaFunA(x, y, z, R, xc, yc, zc, chi, Bo, flag):
 
     return Bx, By, Bz
 
+
 def IDTtoxyz(Inc, Dec, Btot):
     """
         Convert from Inclination, Declination, Total  intensity of earth field to x, y, z
@@ -190,6 +192,7 @@ def IDTtoxyz(Inc, Dec, Btot):
 
     return np.r_[Bx, By, Bz]
 
+
 def MagSphereFreeSpace(x, y, z, R, xc, yc, zc, chi, Bo):
     """
         Computing boundary condition using Congrous sphere method.

simpegPF/Tests/test_forward_PFproblem.py

@@ -25,7 +25,6 @@ class MagFwdProblemTests(unittest.TestCase):
         self.M = M
         self.chi = chi
 
-
     def test_ana_forward(self):
 
         survey = PF.BaseMag.BaseMagSurvey()

simpegPF/Tests/test_magnetics_analytics.py

@@ -0,0 +1,45 @@
+import unittest
+from SimPEG import *
+import matplotlib.pyplot as plt
+import simpegPF as PF
+from scipy.constants import mu_0
+
+
+class MagFwdProblemTests(unittest.TestCase):
+
+    def test_ana_forward(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]
+
+        # 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()

simpegPF/Tests/test_sensitivity_PFproblem.py

@@ -10,9 +10,9 @@ 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)]
+        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
@@ -37,9 +37,7 @@ class MagSensProblemTests(unittest.TestCase):
         rxLoc = np.c_[Utils.mkvc(X), Utils.mkvc(Y), Utils.mkvc(Z)]
         survey.rxLoc = rxLoc
 
-
-
-        prob = PF.Magnetics.MagneticsDiffSecondary(M, mapping=model)
+        prob = PF.Magnetics.Problem3D_DiffSecondary(M, mapping=model)
         prob.pair(survey)
         dpre = survey.dpred(chi)
 
@@ -53,11 +51,10 @@ class MagSensProblemTests(unittest.TestCase):
         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
@@ -83,7 +80,6 @@ class MagSensProblemTests(unittest.TestCase):
         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)
 
 
@@ -118,7 +114,7 @@ class MagSensProblemTests(unittest.TestCase):
             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  )
+            dCdm_A = Div * (Utils.sdiag( Div.T * u ) * dMfMuI * dmudm)
 
             return dCdm_A*v