Working 2.5D fwd (nodal discretization)

On going Jvec and Jtvec

SimPEG/EM/Base.py

@@ -134,7 +134,6 @@ def MeSigmaDeriv(self, u):
         """
         return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma)(u) * self.curModel.sigmaDeriv
 
-
     @property
     def MeSigmaI(self):
         """
@@ -157,7 +156,6 @@ def MeSigmaIDeriv(self, u):
         return dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm))
         # return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma, invMat=True)(u)
 
-
     @property
     def MfRho(self):
         """

SimPEG/EM/Static/DC/FieldsDC_2D.py

@@ -105,3 +105,42 @@ def _j(self, phiSolution, srcList):
 
     def _e(self, phiSolution, srcList):
         raise NotImplementedError
+
+class Fields_ky_N(Fields_ky):
+    knownFields = {'phiSolution':'N'}
+    aliasFields = {
+                    'phi': ['phiSolution','N','_phi'],
+                    'j' : ['phiSolution','E','_j'],
+                    'e' : ['phiSolution','E','_e'],
+                  }
+                  # primary - secondary
+                  # CC variables
+
+    def __init__(self, mesh, survey, **kwargs):
+        Fields_ky.__init__(self, mesh, survey, **kwargs)
+
+    def startup(self):
+        self.prob = self.survey.prob
+
+    def _GLoc(self, fieldType):
+        if fieldType == 'phi':
+            return 'N'
+        elif fieldType == 'e' or fieldType == 'j':
+            return 'E'
+        else:
+            raise Exception('Field type must be phi, e, j')
+
+    def _phi(self, phiSolution, src, kyInd):
+        return phiSolution
+
+    def _phiDeriv_u(self, kyInd, src, v, adjoint = False):
+        return Identity()*v
+
+    def _phiDeriv_m(self, kyInd, src, v, adjoint = False):
+        return Zero()
+
+    def _j(self, phiSolution, srcList):
+        raise NotImplementedError
+
+    def _e(self, phiSolution, srcList):
+        raise NotImplementedError

SimPEG/EM/Static/DC/ProblemDC.py

@@ -267,8 +267,6 @@ def getA(self):
         # Handling Null space of A
         A[0,0] = A[0,0] + 1.
 
-        # if self._makeASymmetric is True:
-        #     return V.T * A
         return A
 
     def getADeriv(self, u, v, adjoint=False):
@@ -293,8 +291,6 @@ def getRHS(self):
         """
 
         RHS = self.getSourceTerm()
-        # if self._makeASymmetric is True:
-        #     return self.Vol.T * RHS
         return RHS
 
     def getRHSDeriv(self, src, v, adjoint=False):

SimPEG/EM/Static/DC/ProblemDC_2D.py

@@ -1,7 +1,7 @@
 from SimPEG import Problem, Utils
 from SimPEG.EM.Base import BaseEMProblem
 from SurveyDC import Survey, Survey_ky
-from FieldsDC_2D import Fields_ky, Fields_ky_CC
+from FieldsDC_2D import Fields_ky, Fields_ky_CC, Fields_ky_N
 from SimPEG.Utils import sdiag
 import numpy as np
 from SimPEG.Utils import Zero
@@ -90,14 +90,10 @@ def Jtvec(self, m, v, f=None):
         dky = np.r_[dky[0], dky]
         y = 0.
 
-
         for src in self.survey.srcList:
-
             for rx in src.rxList:
-
                 Jtv_temp1 = np.zeros(m.size)
                 Jtv_temp0 = np.zeros(m.size)
-
                 for iky in range(self.nky):
                     u_src = f[src, self._solutionType, iky]
                     ky = self.kys[iky]
@@ -152,7 +148,7 @@ class Problem2D_CC(BaseDCProblem_2D):
 
     _solutionType = 'phiSolution'
     _formulation  = 'HJ' # CC potentials means J is on faces
-    fieldsPair    = Fields_ky_CC
+    fieldsPair    = Fields_ky_N
 
     def __init__(self, mesh, **kwargs):
         BaseDCProblem_2D.__init__(self, mesh, **kwargs)
@@ -269,3 +265,85 @@ def setBC(self):
         P_BC, B = self.mesh.getBCProjWF_simple()
         M = B*self.mesh.aveCC2F
         self.Grad = self.Div.T - P_BC*Utils.sdiag(y_BC)*M
+
+class Problem2D_N(BaseDCProblem_2D):
+
+    _solutionType = 'phiSolution'
+    _formulation  = 'EB' # CC potentials means J is on faces
+    fieldsPair    = Fields_ky_N
+
+    def __init__(self, mesh, **kwargs):
+        BaseDCProblem_2D.__init__(self, mesh, **kwargs)
+        # self.setBC()
+
+    @property
+    def MnSigma(self):
+        """
+            Node inner product matrix for \\(\\sigma\\). Used in the E-B formulation
+        """
+        # TODO: only works isotropic sigma
+        sigma = self.curModel.sigma
+        vol = self.mesh.vol
+        MnSigma = Utils.sdiag(self.mesh.aveN2CC.T*(Utils.sdiag(vol)*sigma))
+
+        return MnSigma
+
+    def MnSigmaDeriv(self, u):
+        """
+            Derivative of MnSigma with respect to the model
+        """
+        sigma = self.curModel.sigma
+        sigmaderiv = self.curModel.sigmaDeriv
+        vol = self.mesh.vol
+        return Utils.sdiag(u)*self.mesh.aveN2CC.T*Utils.sdiag(vol) * self.curModel.sigmaDeriv
+
+    def getA(self, ky):
+        """
+
+        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
+        MeSigma = self.MeSigma
+        MnSigma = self.MnSigma
+        Grad = self.mesh.nodalGrad
+        # Get conductivity sigma
+        sigma = self.curModel.sigma
+        A = Grad.T * MeSigma * Grad + ky**2*MnSigma
+
+        # Handling Null space of A
+        A[0,0] = A[0,0] + 1.
+        return A
+
+    def getADeriv(self, ky, u, v, adjoint= False):
+
+        MeSigma = self.MeSigma
+        Grad = self.mesh.nodalGrad
+        sigma = self.curModel.sigma
+        vol = self.mesh.vol
+
+        if adjoint:
+            return Grad.T*(self.MeSigmaDeriv(Grad*u)*v) + ky**2*self.MnSigmaDeriv(u)*v
+        return self.MeSigmaDeriv(Grad*u).T * (Grad*v) + ky**2*self.MnSigmaDeriv(u)*v
+
+    def getRHS(self, ky):
+        """
+        RHS for the DC problem
+
+        q
+        """
+
+        RHS = self.getSourceTerm(ky)
+        return RHS
+
+    def getRHSDeriv(self, ky, src, v, adjoint=False):
+        """
+        Derivative of the right hand side with respect to the model
+        """
+        # TODO: add qDeriv for RHS depending on m
+        # qDeriv = src.evalDeriv(self, ky, adjoint=adjoint)
+        # return qDeriv
+        return Zero()

SimPEG/EM/Static/DC/SrcDC.py

@@ -47,8 +47,7 @@ def eval(self, prob):
             q = np.zeros(prob.mesh.nC)
             q[inds] = self.current * np.r_[1.]
         elif prob._formulation == 'EB':
-            inds = closestPoints(prob.mesh, self.loc)
-            q = np.zeros(prob.mesh.nN)
-            q[inds] = self.current * np.r_[1.]
+            q = prob.mesh.getInterpolationMat(self.loc, locType='N').todense()
+            q = self.current * mkvc(q)
         return q
 

SimPEG/EM/Static/DC/__init__.py

@@ -1,5 +1,5 @@
 from ProblemDC import Problem3D_CC, Problem3D_N
-from ProblemDC_2D import Problem2D_CC
+from ProblemDC_2D import Problem2D_CC, Problem2D_N
 from SurveyDC import Survey, Survey_ky
 import SrcDC as Src #Pole
 import RxDC as Rx