Merge branch 'dev' into ref/dev

Conflicts:
	docs/examples/DC_Forward_PseudoSection.rst

SimPEG/EM/Analytics/DC.py

@@ -0,0 +1,118 @@
+import numpy as np
+from scipy.constants import mu_0, pi
+from scipy import special
+
+def DCAnalyticHalf(txloc, rxlocs, sigma, earth_type="wholespace"):
+    """
+        Analytic solution for electric potential from a postive pole
+
+        :param array txloc: a xyz location of A (+) electrode (np.r_[xa, ya, za])
+        :param list rxlocs: xyz locations of M (+) and N (-) electrodes [M, N]
+
+            e.g.
+            rxlocs = [M, N]
+            M: xyz locations of M (+) electrode (np.c_[xmlocs, ymlocs, zmlocs])
+            N: xyz locations of N (-) electrode (np.c_[xnlocs, ynlocs, znlocs])
+
+        :param float or complex sigma: values of conductivity
+        :param string earth_type: values of conductivity ("wholsespace" or "halfspace")
+
+    """
+    M = rxlocs[0]
+    N = rxlocs[1]
+
+    rM = np.sqrt( (M[:,0]-txloc[0])**2 + (M[:,1]-txloc[1])**2 + (M[:,2]-txloc[1])**2 )
+    rN = np.sqrt( (N[:,0]-txloc[0])**2 + (N[:,1]-txloc[1])**2 + (N[:,2]-txloc[1])**2 )
+
+    phiM = 1./(4*np.pi*rM*sigma)
+    phiN = 1./(4*np.pi*rN*sigma)
+    phi = phiM - phiN
+
+    if earth_type == "halfspace":
+        phi *= 2
+
+    return phi
+
+deg2rad  = lambda deg: deg/180.*np.pi
+rad2deg  = lambda rad: rad*180./np.pi
+
+def DCAnalyticSphere(txloc, rxloc, xc, radius, sigma, sigma1, \
+                 field_type = "secondary", order=12, halfspace=False):
+# def DCSpherePointCurrent(txloc, rxloc, xc, radius, rho, rho1, \
+#                  field_type = "secondary", order=12):
+    """
+
+        Parameters:
+
+            :param array txloc: A (+) current electrode location (x,y,z)
+            :param array xc: x center of depressed sphere
+            :param array rxloc: M(+) electrode locations / (Nx3 array, # of electrodes)
+
+            :param float radius: radius (float): radius of the sphere (m)
+            :param float rho: resistivity of the background (ohm-m)
+            :param float rho1: resistivity of the sphere
+            :param string field_type: : "secondary", "total", "primary"
+                  (default="secondary")
+                  "secondary": secondary potential only due to sphere
+                  "primary": primary potential from the point source
+                  "total": "secondary"+"primary"
+            :param float order: maximum order of Legendre polynomial (default=12)
+
+        Written by Seogi Kang (skang@eos.ubc.ca)
+        Ph.D. Candidate of University of British Columbia, Canada
+
+    """
+
+    Pleg = []
+    # Compute Legendre Polynomial
+    for i in range(order):
+        Pleg.append(special.legendre(i, monic=0))
+
+
+    rho = 1./sigma
+    rho1 = 1./sigma1
+
+    # Center of the sphere should be aligned in txloc in y-direction
+    yc = txloc[1]
+    xyz = np.c_[rxloc[:,0]-xc, rxloc[:,1]-yc, rxloc[:,2]]
+    r = np.sqrt( (xyz**2).sum(axis=1) )
+
+    x0 = abs(txloc[0]-xc)
+
+    costheta = xyz[:,0]/r * (txloc[0]-xc)/x0
+    phi = np.zeros_like(r)
+    R = (r**2+x0**2.-2.*r*x0*costheta)**0.5
+    # primary potential in a whole space
+    prim = rho*1./(4*np.pi*R)
+
+    if field_type =="primary":
+        return prim
+
+    sphind = r < radius
+    out = np.zeros_like(r)
+    for n in range(order):
+        An, Bn = AnBnfun(n, radius, x0, rho, rho1)
+        dumout = An*r[~sphind]**(-n-1.)*Pleg[n](costheta[~sphind])
+        out[~sphind] += dumout
+        dumin = Bn*r[sphind]**(n)*Pleg[n](costheta[sphind])
+        out[sphind] += dumin
+
+    out[~sphind] += prim[~sphind]
+
+    if halfspace:
+        scale = 2
+    else:
+        scale = 1
+
+    if field_type == "secondary":
+        return scale*(out-prim)
+    elif field_type == "total":
+        return scale*out
+
+def AnBnfun(n, radius, x0, rho, rho1, I=1.):
+    const = I*rho/(4*np.pi)
+    bunmo = n*rho + (n+1)*rho1
+    An = const * radius**(2*n+1) / x0 ** (n+1.) * n * \
+        (rho1-rho) / bunmo
+    Bn = const * 1. / x0 ** (n+1.) * (2*n+1) * (rho1) / bunmo
+    return An, Bn

SimPEG/EM/Analytics/__init__.py

@@ -1,3 +1,4 @@
 from TDEM import hzAnalyticDipoleT
 from FDEM import hzAnalyticDipoleF
 from FDEMcasing import *
+from DC import DCAnalyticHalf, DCAnalyticSphere

SimPEG/EM/Base.py

@@ -1,6 +1,7 @@
 from SimPEG import Survey, Problem, Utils, Models, Maps, PropMaps, np, sp, Solver as SimpegSolver
 from scipy.constants import mu_0
 
+
 class EMPropMap(Maps.PropMap):
     """
         Property Map for EM Problems. The electrical conductivity (\\(\\sigma\\)) is the default inversion property, and the default value of the magnetic permeability is that of free space (\\(\\mu = 4\\pi\\times 10^{-7} \\) H/m)
@@ -61,6 +62,15 @@ def Me(self):
             self._Me = self.mesh.getEdgeInnerProduct()
         return self._Me
 
+    @property
+    def MeI(self):
+        """
+            Edge inner product matrix
+        """
+        if getattr(self, '_MeI', None) is None:
+            self._MeI = self.mesh.getEdgeInnerProduct(invMat=True)
+        return self._MeI
+
     @property
     def Mf(self):
         """
@@ -70,6 +80,20 @@ def Mf(self):
             self._Mf = self.mesh.getFaceInnerProduct()
         return self._Mf
 
+    @property
+    def MfI(self):
+        """
+            Face inner product matrix
+        """
+        if getattr(self, '_MfI', None) is None:
+            self._MfI = self.mesh.getFaceInnerProduct(invMat=True)
+        return self._MfI
+
+    @property
+    def Vol(self):
+        if getattr(self, '_Vol', None) is None:
+            self._Vol = Utils.sdiag(self.mesh.vol)
+        return self._Vol
 
     # ----- Magnetic Permeability ----- #
     @property
@@ -127,7 +151,6 @@ def MeSigmaDeriv(self, u):
         """
         return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma)(u) * self.curModel.sigmaDeriv
 
-
     @property
     def MeSigmaI(self):
         """
@@ -146,10 +169,7 @@ def MeSigmaIDeriv(self, u):
 
         dMeSigmaI_dI = -self.MeSigmaI**2
         dMe_dsig = self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma)(u)
-        dsig_dm = self.curModel.sigmaDeriv
-        return dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm))
-        # return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma, invMat=True)(u)
-
+        return dMeSigmaI_dI * ( dMe_dsig * self.curModel.sigmaDeriv )
 
     @property
     def MfRho(self):
@@ -165,8 +185,7 @@ def MfRhoDeriv(self,u):
         """
             Derivative of :code:`MfRho` with respect to the model.
         """
-        return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) * (-Utils.sdiag(self.curModel.rho**2) * self.curModel.sigmaDeriv)
-        # self.curModel.rhoDeriv
+        return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) * self.curModel.rhoDeriv
 
     @property
     def MfRhoI(self):
@@ -183,7 +202,10 @@ def MfRhoIDeriv(self,u):
         """
             Derivative of :code:`MfRhoI` with respect to the model.
         """
-        return self.mesh.getFaceInnerProductDeriv(self.curModel.rho, invMat=True)(u) * self.curModel.rhoDeriv
+
+        dMfRhoI_dI = -self.MfRhoI**2
+        dMf_drho = self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u)
+        return dMfRhoI_dI * ( dMf_drho * self.curModel.rhoDeriv )
 
 class BaseEMSurvey(Survey.BaseSurvey):
 
@@ -192,7 +214,7 @@ def __init__(self, srcList, **kwargs):
         self.srcList = srcList
         Survey.BaseSurvey.__init__(self, **kwargs)
 
-    def eval(self, u):
+    def eval(self, f):
         """
         Project fields to receiver locations
         :param Fields u: fields object
@@ -202,8 +224,8 @@ def eval(self, u):
         data = Survey.Data(self)
         for src in self.srcList:
             for rx in src.rxList:
-                data[src, rx] = rx.eval(src, self.mesh, u)
+                data[src, rx] = rx.eval(src, self.mesh, f)
         return data
 
-    def evalDeriv(self, u):
+    def evalDeriv(self, f):
         raise Exception('Use Receivers to project fields deriv.')

SimPEG/EM/FDEM/FieldsFDEM.py

@@ -160,9 +160,9 @@ def _jDeriv(self, src, du_dm_v, v, adjoint = False):
             return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint)
         return np.array(self._jDeriv_u(src, du_dm_v, adjoint) + self._jDeriv_m(src, v, adjoint), dtype = complex)
 
-class Fields_e(Fields):
+class Fields3D_e(Fields):
     """
-    Fields object for Problem_e.
+    Fields object for Problem3D_e.
 
     :param Mesh mesh: mesh
     :param Survey survey: survey
@@ -181,7 +181,7 @@ class Fields_e(Fields):
                   }
 
     def __init__(self, mesh, survey, **kwargs):
-        Fields.__init__(self,mesh,survey,**kwargs)
+        Fields.__init__(self, mesh, survey, **kwargs)
 
     def startup(self):
         self.prob = self.survey.prob
@@ -426,9 +426,9 @@ def _hDeriv_m(self, src, v, adjoint = False):
 
 
 
-class Fields_b(Fields):
+class Fields3D_b(Fields):
     """
-    Fields object for Problem_b.
+    Fields object for Problem3D_b.
 
     :param Mesh mesh: mesh
     :param Survey survey: survey
@@ -693,9 +693,9 @@ def _hDeriv_m(self, src, v, adjoint=False):
         return Zero()
 
 
-class Fields_j(Fields):
+class Fields3D_j(Fields):
     """
-    Fields object for Problem_j.
+    Fields object for Problem3D_j.
 
     :param Mesh mesh: mesh
     :param Survey survey: survey
@@ -988,9 +988,9 @@ def _bDeriv_m(self, src, v, adjoint=False):
         return 1./(1j * omega(src.freq)) * VI * (self._aveE2CCV * ( s_mDeriv(v) - self._edgeCurl.T * ( self._MfRhoDeriv(jSolution) * v ) ) )
 
 
-class Fields_h(Fields):
+class Fields3D_h(Fields):
     """
-    Fields object for Problem_h.
+    Fields object for Problem3D_h.
 
     :param Mesh mesh: mesh
     :param Survey survey: survey