From 31e7bcb3f0790b6e6561dc4541d037534ebddbc3 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 22 Jun 2016 11:24:46 -0600 Subject: [PATCH] Analytic test --- simpegPF/MagAnalytics.py | 5 ++- simpegPF/Tests/test_forward_PFproblem.py | 1 - simpegPF/Tests/test_magnetics_analytics.py | 45 ++++++++++++++++++++ simpegPF/Tests/test_sensitivity_PFproblem.py | 16 +++---- 4 files changed, 55 insertions(+), 12 deletions(-) create mode 100644 simpegPF/Tests/test_magnetics_analytics.py diff --git a/simpegPF/MagAnalytics.py b/simpegPF/MagAnalytics.py index b4a69af7..64123fcc 100644 --- a/simpegPF/MagAnalytics.py +++ b/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. diff --git a/simpegPF/Tests/test_forward_PFproblem.py b/simpegPF/Tests/test_forward_PFproblem.py index 5f7d13a5..32e24f59 100644 --- a/simpegPF/Tests/test_forward_PFproblem.py +++ b/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() diff --git a/simpegPF/Tests/test_magnetics_analytics.py b/simpegPF/Tests/test_magnetics_analytics.py new file mode 100644 index 00000000..a0cd47cc --- /dev/null +++ b/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() diff --git a/simpegPF/Tests/test_sensitivity_PFproblem.py b/simpegPF/Tests/test_sensitivity_PFproblem.py index b2a79549..8e27a0b4 100644 --- a/simpegPF/Tests/test_sensitivity_PFproblem.py +++ b/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