start of DC refactor in EM/Static

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)
@@ -70,6 +71,12 @@ class BaseEMProblem(Problem.BaseProblem):
             self._Mf = self.mesh.getFaceInnerProduct()
         return self._Mf
 
+    @property
+    def Vol(self):
+        if getattr(self, '_Vol', None) is None:
+            self._Vol = Utils.sdiag(self.mesh.vol)
+        return self._Vol
+
 
     # ----- Magnetic Permeability ----- #
     @property
@@ -192,7 +199,7 @@ class BaseEMSurvey(Survey.BaseSurvey):
         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 +209,8 @@ class BaseEMSurvey(Survey.BaseSurvey):
         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

@@ -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

SimPEG/EM/FDEM/SrcFDEM.py

@@ -14,7 +14,6 @@ class BaseSrc(Survey.BaseSrc):
 
     def eval(self, prob):
         """
-        Evaluate the source terms.
         - :math:`s_m` : magnetic source term
         - :math:`s_e` : electric source term
 

SimPEG/EM/FDEM/SurveyFDEM.py

@@ -63,9 +63,9 @@ class Rx(SimPEG.Survey.BaseRx):
         """Component projection (real/imag)"""
         return self.knownRxTypes[self.rxType][2]
 
-    def projGLoc(self, u):
+    def projGLoc(self, f):
         """Grid Location projection (e.g. Ex Fy ...)"""
-        return u._GLoc(self.rxType[0]) + self.knownRxTypes[self.rxType][1]
+        return f._GLoc(self.rxType[0]) + self.knownRxTypes[self.rxType][1]
 
     def eval(self, src, mesh, f):
         """

SimPEG/EM/Static/DC/FieldsDC.py

@@ -0,0 +1,72 @@
+import SimPEG
+from SimPEG.Utils import Identity, Zero
+
+class Fields(SimPEG.Problem.Fields):
+    knownFields = {}
+    dtype = float
+
+    def _phiDeriv(self, src, du_dm_v, v, adjoint=False):
+        if getattr(self, '_phiDeriv_u', None) is None or getattr(self, '_phiDeriv_m', None) is None:
+            raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %self.knownFields.keys()[0])
+
+        if adjoint:
+            return self._phiDeriv_u(src, v, adjoint), self._phiDeriv_m(src, v, adjoint)
+        return np.array(self._phiDeriv_u(src, du_dm_v, adjoint) + self._phiDeriv_m(src, v, adjoint), dtype = complex)
+
+    def _eDeriv(self, src, du_dm_v, v, adjoint=False):
+        if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None:
+            raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0])
+
+        if adjoint:
+            return self._eDeriv_u(src, v, adjoint), self._eDeriv_m(src, v, adjoint)
+        return np.array(self._eDeriv_u(src, du_dm_v, adjoint) + self._eDeriv_m(src, v, adjoint), dtype = complex)
+
+    def _jDeriv(self, src, du_dm_v, v, adjoint=False):
+        if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None:
+            raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0])
+
+        if adjoint:
+            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_CC(Fields):
+    knownFields = {'phiSolution':'CC'}
+    aliasFields = {
+                    'phi': ['phiSolution','CC','_phi'],
+                    'j' : ['phiSolution','F','_j'],
+                    'e' : ['phiSolution','F','_e'],
+                  }
+                  # primary - secondary
+                  # CC variables
+
+    def __init__(self, mesh, survey, **kwargs):
+        Fields.__init__(self, mesh, survey, **kwargs)
+
+    def startup(self):
+        self.prob = self.survey.prob
+
+    def _GLoc(self, fieldType):
+        if fieldType == 'phi':
+            return 'CC'
+        elif fieldType == 'e' or fieldType == 'j':
+            return 'F'
+        else:
+            raise Exception('Field type must be phi, e, j')
+
+    def _phi(self, phiSolution, srcList):
+        return phiSolution
+
+    def _phiDeriv_u():
+        return Identity()
+
+    def _phiDeriv_m():
+        return Zero()
+
+    def _j(self, phiSolution, srcList):
+        raise NotImplementedError
+
+    def _e(self, phiSolution, srcList):
+        raise NotImplementedError
+
+

SimPEG/EM/Static/DC/ProblemDC.py

@@ -0,0 +1,79 @@
+from SimPEG import Problem
+from SimPEG.EM.Base import BaseEMProblem
+from SurveyDC import Survey
+from FieldsDC import Fields, Fields_CC
+
+class BaseDCProblem(BaseEMProblem):
+
+    surveyPair = Survey
+    fieldsPair = Fields
+
+    def fields(self, m):
+        self.curModel = m
+        f = self.fieldsPair(self.mesh, self.survey)
+        A = self.getA()
+        self.Ainv = self.Solver(A, **self.solverOpts)
+        RHS = self.getRHS()
+        u = self.Ainv * RHS
+        Srcs = self.survey.srcList
+        f[Srcs, self._solutionType] = u
+        return f
+
+    def Jvec(self, m, v, f=None):
+        raise NotImplementedError
+
+    def Jtvec(self, m, v, f=None):
+        raise NotImplementedError
+
+    def getSourceTerm(self):
+        """
+        takes concept of source and turns it into a matrix
+        """
+        raise NotImplementedError
+
+
+
+class Problem3D_CC(BaseDCProblem):
+
+    _solutionType = 'phiSolution'
+    _formulation  = 'HJ' # CC potentials means J is on faces
+    fieldsPair    = Fields_CC
+
+    def __init__(self, mesh, **kwargs):
+        BaseDCProblem.__init__(self, mesh, **kwargs)
+
+
+    def getA(self):
+        """
+
+        Make the A matrix for the cell centered DC resistivity problem
+
+        A = D MfRhoI D^\\top V
+
+        """
+
+        # TODO: this won't work for full anisotropy
+
+        D = self.mesh.faceDiv
+        MfRhoI = self.MfRhoI
+        V = self.Vol
+        A = D * ( MfRhoI * ( D.T * V ) )
+
+        if self._makeASymmetric is True:
+            return V.T * A
+        return A
+
+
+    def getRHS(self):
+        """
+        RHS for the DC problem
+
+        q
+        """
+
+        RHS = self.getSourceTerm()
+        if self._makeASymmetric is True:
+            return self.Vol.T * RHS
+        return RHS
+
+

SimPEG/EM/Static/DC/RxDC.py

@@ -0,0 +1,66 @@
+import SimPEG
+# from SimPEG.EM.Base import BaseEMSurvey
+from SimPEG.Utils import Zero, closestPoints
+
+class BaseRx(SimPEG.Survey.BaseRx):
+    loc = None
+    rxType = None
+
+    knownRxTypes = {
+                    'phi':['phi',None],
+                    'ex':['e','x'],
+                    'ey':['e','y'],
+                    'ez':['e','z'],
+                    'jx':['j','x'],
+                    'jy':['j','y'],
+                    'jz':['j','z'],
+                    }
+
+    def __init__(self, **kwargs):
+        SimPEG.Survey.BaseRx.__init__(locs, rxType, **kwargs)
+
+    @property
+    def projField(self):
+        """Field Type projection (e.g. e b ...)"""
+        return self.knownRxTypes[self.rxType][0]
+
+    def projGLoc(self, f):
+        """Grid Location projection (e.g. Ex Fy ...)"""
+        comp = self.knownRxTypes[self.rxType][1]
+        if comp is not None:
+            return f._GLoc(self.rxType[0]) + comp
+        return f._GLoc(self.rxType[0])
+
+    def eval(self, src, mesh, f):
+        P = self.getP(self.prob.mesh)
+        return P*f[src, self.projField]
+
+# DC.Rx.Dipole(locs)
+class Dipole(BaseRx):
+
+    def __init__(self, locsM, locsN, rxType = 'phi', **kwargs):
+        assert locsM.shape == locsN.shape, 'locsM and locsN need to be the same size'
+        self.locs = [locsM, locsN]
+        BaseRx.__init__(self)
+
+    @property
+    def nD(self):
+        """Number of data in the receiver."""
+        return self.locs[0].shape[0]
+
+    def getP(self,mesh):
+        if mesh in self._Ps:
+            return self._Ps[mesh]
+
+        P0 = mesh.getInterpolationMat(self.locs[0], self.projGLoc)
+        P1 = mesh.getInterpolationMat(self.locs[1], self.projGLoc)
+        P = P0 - P1
+
+        if self.storeProjections:
+            self._Ps[mesh] = P
+        return P
+
+
+# class Pole(BaseRx):
+
+

SimPEG/EM/Static/DC/SrcDC.py

@@ -0,0 +1,60 @@
+import SimPEG
+# from SimPEG.EM.Base import BaseEMSurvey
+from SimPEG.Utils import Zero, closestPoints
+
+class BaseSrc(SimPEG.Survey.BaseSrc):
+
+    current = 1
+    loc = None
+
+    def __init__(self, rxList, **kwargs):
+        SimPEG.Survey.BaseSrc.__init__(self, rxList, **kwargs)
+
+    def eval(self, prob):
+        raise NotImplementedError
+
+    def evalDeriv(self, prob):
+        Zero()
+
+
+class Dipole(BaseSrc):
+
+    def __init__(self, rxList, locA, locB, **kwargs):
+        assert locA.shape == locB.shape, 'Shape of locA and locB should be the same'
+        self.loc = [locA, locB]
+        BaseSrc.__init__(self, rxList, **kwargs)
+
+    def eval(self, prob):
+        if prob._formulation == 'HJ':
+            inds = closestPoints(prob.mesh, self.loc)
+            q = np.zeros(prob.mesh.nC)
+            q[inds] = self.current * np.r_[1., -1.]
+        elif prob._formulation == 'EB':
+            # TODO: there is probably a faster way to do this
+            # Utils.cellNodes , Utils.cellFaces, Utils.cellEdges
+            raise NotImplementedError
+        return q
+
+    # def bc_contribution
+
+
+# How to treat boundary conditions here
+
+class Pole(BaseSrc):
+
+    def __init__(self, rxList, loc, **kwargs):
+        BaseSrc.__init__(self, rxList, loc=loc, **kwargs)
+
+    def eval(self, prob):
+        if prob._formulation == 'HJ':
+            inds = closestPoints(prob.mesh, self.loc)
+            q = np.zeros(prob.mesh.nC)
+            q[inds] = self.current * np.r_[1.]
+        elif prob._formulation == 'EB':
+            # TODO: there is probably a faster way to do this
+            # Utils.cellNodes , Utils.cellFaces, Utils.cellEdges
+            raise NotImplementedError
+        return q
+
+    # def bc_contribution
+

SimPEG/EM/Static/DC/SurveyDC.py

@@ -0,0 +1,18 @@
+import SimPEG
+from SimPEG.EM.Base import BaseEMSurvey
+from SimPEG import sp
+from SimPEG.Utils import Zero, Identity
+from RxDC import BaseRx
+from SrcDC import BaseSrc
+
+class Survey(BaseEMSurvey):
+    rxPair  = BaseRx
+    srcPair = BaseSrc
+
+    def __init__(self, srcList, **kwargs):
+        self.srcList = srcList
+        BaseEMSurvey.__init__(self, srcList, **kwargs)
+
+
+
+

SimPEG/EM/Static/DC/__init__.py

@@ -0,0 +1,5 @@
+from ProblemDC import Problem3D_CC
+from SurveyDC import Survey
+import SrcDC as Src #Pole
+import RxDC as Rx
+from FieldsDC import Fields_CC