start of Casing analytic and test

2026-07-04 04:45:06 +08:00 · 2015-05-06 06:47:53 -07:00
parent 96b0d6a38c
commit 1e52c365bc
4 changed files with 152 additions and 1 deletions
@@ -0,0 +1,91 @@
+from SimPEG import Utils, np
+from scipy.constants import mu_0, epsilon_0
+from simpegEM.Utils.EMUtils import k
+
+def getKc(freq,sigma,a,b,mu=mu_0,eps=epsilon_0):
+    a = float(a)
+    b = float(b)
+    return 2. * np.sqrt(b / a) * np.exp(-1j*k(freq,sigma,mu,eps)*(b-a))
+
+def _r2(xyz):
+    return np.sum(xyz**2,1)
+
+def _getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu=mu_0,eps=epsilon_0,moment=1.):
+    Kc1 = getKc(freq,sigma[1],a,b,mu,eps)
+
+    nobs = obsloc.shape[0]
+    dxyz = obsloc - np.c_[np.ones(nobs)]*np.r_[srcloc]
+    
+    r2 = _r2(dxyz[:,:2])
+    sqrtr2z2 = np.sqrt(r2 + dxyz[:,2]**2)
+    k2 = k(freq,sigma[2],mu,eps) 
+    
+    return Kc1 * moment / (4.*np.pi) *np.exp(-1j*k2*sqrtr2z2) / sqrtr2z2
+
+
+def _getCasingHertzMagDipoleDeriv_r(srcloc,obsloc,freq,sigma,a,b,mu=mu_0,eps=epsilon_0,moment=1.):
+    HertzZ = _getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment)
+    
+    nobs = obsloc.shape[0]
+    dxyz = obsloc - np.c_[np.ones(nobs)]*np.r_[srcloc]
+    
+    r2 = _r2(dxyz[:,:2])
+    sqrtr2z2 = np.sqrt(r2 + dxyz[:,2]**2)
+    k2 = k(freq,sigma[2],mu,eps) 
+    
+    return -HertzZ * np.sqrt(r2) / sqrtr2z2 * (1j*k2 + 1./ sqrtr2z2)
+
+
+def _getCasingHertzMagDipoleDeriv_z(srcloc,obsloc,freq,sigma,a,b,mu=mu_0,eps=epsilon_0,moment=1.):
+    HertzZ = _getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment)
+    
+    nobs = obsloc.shape[0]
+    dxyz = obsloc - np.c_[np.ones(nobs)]*np.r_[srcloc]
+    
+    r2z2 = _r2(dxyz)
+    sqrtr2z2 = np.sqrt(r2z2)
+    k2 = k(freq,sigma[2],mu,eps) 
+    
+    return -HertzZ*dxyz[:,2] /sqrtr2z2 * (1j*k2 + 1./sqrtr2z2)
+
+def _getCasingHertzMagDipole2Deriv_z_r(srcloc,obsloc,freq,sigma,a,b,mu=mu_0,eps=epsilon_0,moment=1.):
+    HertzZ = _getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment)
+    dHertzZdr = _getCasingHertzMagDipoleDeriv_r(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment)
+    
+    nobs = obsloc.shape[0]
+    dxyz = obsloc - np.c_[np.ones(nobs)]*np.r_[srcloc]
+
+    r2 = _r2(dxyz[:,:2])
+    r = np.sqrt(r2)
+    z = dxyz[:,2]
+    sqrtr2z2 = np.sqrt(r2 + z**2)
+    k2 = k(freq,sigma[2],mu,eps) 
+    
+    return dHertzZdr*(-z/sqrtr2z2)*(1j*k2+1./sqrtr2z2) + HertzZ*(z*r/sqrtr2z2**3)*(1j*k2 + 2./sqrtr2z2)
+
+def _getCasingHertzMagDipole2Deriv_z_z(srcloc,obsloc,freq,sigma,a,b,mu=mu_0,eps=epsilon_0,moment=1.):
+    HertzZ = _getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment)
+    dHertzZdz = _getCasingHertzMagDipoleDeriv_z(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment)
+    
+    nobs = obsloc.shape[0]
+    dxyz = obsloc - np.c_[np.ones(nobs)]*np.r_[srcloc]
+
+    r2 = _r2(dxyz[:,:2])
+    r = np.sqrt(r2)
+    z = dxyz[:,2]
+    sqrtr2z2 = np.sqrt(r2 + z**2)
+    k2 = k(freq,sigma[2],mu,eps) 
+    
+    return (dHertzZdz*z + HertzZ)/sqrtr2z2*(-1j*k2 - 1./sqrtr2z2) + HertzZ*z/sqrtr2z2**3*(1j*k2*z + 2.*z/sqrtr2z2)
+    
+def getCasingEphiMagDipole(srcloc,obsloc,freq,sigma,a,b,mu=mu_0,eps=epsilon_0,moment=1):
+    return 1j * omega(freq) * mu * _getCasingHertzMagDipoleDeriv_r(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment)
+
+def getCasingHrMagDipole(srcloc,obsloc,freq,sigma,a,b,mu=mu_0,eps=epsilon_0,moment=1):
+    return _getCasingHertzMagDipole2Deriv_z_r(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment)
+
+def getCasingHzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu=mu_0,eps=epsilon_0,moment=1):
+    d2HertzZdz2 = _getCasingHertzMagDipole2Deriv_z_z(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment)
+    k2 = k(freq,sigma[2],mu,eps)
+    HertzZ(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment)
+    return d2HertzZdz2 + k2**2 * HertzZ
@@ -1,2 +1,3 @@
 from TDEM import hzAnalyticDipoleT
 from FDEM import hzAnalyticDipoleF
+from FDEMcasing import *
@@ -0,0 +1,59 @@
+from SimPEG import Tests, Utils, np
+import simpegEM.Analytics.FDEMcasing as Casing
+import unittest
+
+
+n = 50.
+freq = 1.
+a = 5e-2
+b = a + 1e-2
+sigma = np.r_[1., 5.5e6, 1e-1]
+srcloc = np.r_[0., 0., 0.]
+xobs = np.random.rand(n)+10.
+yobs = np.zeros(n)
+zobs = np.random.randn(n)
+plotit = False
+
+def CasingMagDipoleDeriv_r(x):
+    obsloc = np.vstack([x, yobs, zobs]).T 
+    
+    f = Casing._getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b)
+    g = Utils.sdiag(Casing._getCasingHertzMagDipoleDeriv_r(srcloc,obsloc,freq,sigma,a,b))
+
+    return  f,g 
+
+def CasingMagDipoleDeriv_z(z):
+    obsloc = np.vstack([xobs, yobs, z]).T 
+    
+    f = Casing._getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b)
+    g = Utils.sdiag(Casing._getCasingHertzMagDipoleDeriv_z(srcloc,obsloc,freq,sigma,a,b))
+
+    return  f,g 
+
+def CasingMagDipole2Deriv_z_r(x):
+    obsloc = np.vstack([x, yobs, zobs]).T
+    
+    f = Casing._getCasingHertzMagDipoleDeriv_z(srcloc,obsloc,freq,sigma,a,b)
+    g = Utils.sdiag(Casing._getCasingHertzMagDipole2Deriv_z_r(srcloc,obsloc,freq,sigma,a,b))
+    
+    return f,g
+
+def CasingMagDipole2Deriv_z_z(z):
+    obsloc = np.vstack([xobs, yobs, z]).T 
+    
+    f = Casing._getCasingHertzMagDipoleDeriv_z(srcloc,obsloc,freq,sigma,a,b)
+    g = Utils.sdiag(Casing._getCasingHertzMagDipole2Deriv_z_z(srcloc,obsloc,freq,sigma,a,b))
+    
+    return f,g
+
+
+
+class Casing_DerivTest(unittest.TestCase):
+	Tests.checkDerivative(CasingMagDipoleDeriv_r,np.ones(n)*10+np.random.randn(n),plotIt=False)
+	Tests.checkDerivative(CasingMagDipoleDeriv_z,np.random.randn(n),plotIt=False)
+	Tests.checkDerivative(CasingMagDipole2Deriv_z_r,np.ones(n)*10+np.random.randn(n),plotIt=False)
+	Tests.checkDerivative(CasingMagDipole2Deriv_z_z,np.random.randn(n),plotIt=False)
+
+
+if __name__ == '__main__':
+    unittest.main()
@@ -11,7 +11,7 @@ class FieldsTest(unittest.TestCase):
        srcLoc = np.r_[0,0,0.]
        rxList0 = EM.FDEM.RxFDEM(XYZ, 'exi')
        Src0 = EM.FDEM.SrcFDEM_MagDipole(srcLoc, 3., [rxList0])
-        rxList1 = EM.FDEM.RxFDEM(XYZ, 'bxi')
+        rxList1 = EM.FDEM.RxFDEM(XYZ, 'bxi') 
        Src1 = EM.FDEM.SrcFDEM_MagDipole(srcLoc, 3., [rxList1])
        rxList2 = EM.FDEM.RxFDEM(XYZ, 'bxi')
        Src2 = EM.FDEM.SrcFDEM_MagDipole(srcLoc, 2., [rxList2])