workking nodal discretizations

SimPEG/EM/Static/DC/BoundaryUtils.py

@@ -0,0 +1,158 @@
+import numpy as np
+
+def getxBCyBC_CC(mesh, alpha, beta, gamma):
+# def getxBCyBC(mesh, alpha, beta, gamma):
+    """
+    This is a subfunction generating mixed-boundary condition:
+
+    .. math::
+
+        \nabla \cdot \vec{j} = -\nabla \cdot \vec{j}_s = q
+
+        \rho \vec{j} = -\nabla \phi \phi
+
+        \alpha \phi + \beta \frac{\partial \phi}{\partial r} = \gamma \ at \ r = \partial \Omega
+
+        xBC = f_1(\alpha, \beta, \gamma)
+        yBC = f(\alpha, \beta, \gamma)
+
+    Computes xBC and yBC for cell-centered discretizations
+    """
+    if mesh.dim == 1: #1D
+        if (len(alpha) != 2 or len(beta) != 2 or len(gamma) != 2):
+            raise Exception("Lenght of list, alpha should be 2")
+        fCCxm,fCCxp = mesh.cellBoundaryInd
+        nBC = fCCxm.sum()+fCCxp.sum()
+        h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp]
+
+        alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0]
+        alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1]
+
+        h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp]
+
+        a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm)
+        b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
+        a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp)
+        b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp)
+
+        xBC_xm = 0.5*a_xm
+        xBC_xp = 0.5*a_xp/b_xp
+        yBC_xm = 0.5*(1.-b_xm)
+        yBC_xp = 0.5*(1.-1./b_xp)
+
+        xBC = np.r_[xBC_xm, xBC_xp]
+        yBC = np.r_[yBC_xm, yBC_xp]
+
+    elif mesh.dim == 2: #2D
+        if (len(alpha) != 4 or len(beta) != 4 or len(gamma) != 4):
+            raise Exception("Lenght of list, alpha should be 4")
+
+        fCCxm,fCCxp,fCCym,fCCyp = mesh.cellBoundaryInd
+        fxm,fxp,fym,fyp = mesh.faceBoundaryInd
+        nBC = fCCxm.sum()+fCCxp.sum()+fCCxm.sum()+fCCxp.sum()
+        h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp]
+        h_ym, h_yp = mesh.gridCC[fCCym], mesh.gridCC[fCCyp]
+
+        alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0]
+        alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1]
+        alpha_ym, beta_ym, gamma_ym = alpha[2], beta[2], gamma[2]
+        alpha_yp, beta_yp, gamma_yp = alpha[3], beta[3], gamma[3]
+
+        h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0]
+        h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1]
+
+        a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm)
+        b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
+        a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp)
+        b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp)
+
+        a_ym = gamma_ym/(0.5*alpha_ym-beta_ym/h_ym)
+        b_ym = (0.5*alpha_ym+beta_ym/h_ym)/(0.5*alpha_ym-beta_ym/h_ym)
+        a_yp = gamma_yp/(0.5*alpha_yp-beta_yp/h_yp)
+        b_yp = (0.5*alpha_yp+beta_yp/h_yp)/(0.5*alpha_yp-beta_yp/h_yp)
+
+        xBC_xm = 0.5*a_xm
+        xBC_xp = 0.5*a_xp/b_xp
+        yBC_xm = 0.5*(1.-b_xm)
+        yBC_xp = 0.5*(1.-1./b_xp)
+        xBC_ym = 0.5*a_ym
+        xBC_yp = 0.5*a_yp/b_yp
+        yBC_ym = 0.5*(1.-b_ym)
+        yBC_yp = 0.5*(1.-1./b_yp)
+
+        sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm], np.arange(mesh.nFx)[fxp]])
+        sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym], np.arange(mesh.nFy)[fyp]])
+
+        xBC_x = np.r_[xBC_xm, xBC_xp][sortindsfx]
+        xBC_y = np.r_[xBC_ym, xBC_yp][sortindsfy]
+        yBC_x = np.r_[yBC_xm, yBC_xp][sortindsfx]
+        yBC_y = np.r_[yBC_ym, yBC_yp][sortindsfy]
+
+        xBC = np.r_[xBC_x, xBC_y]
+        yBC = np.r_[yBC_x, yBC_y]
+
+    elif mesh.dim == 3: #3D
+        if (len(alpha) != 6 or len(beta) != 6 or len(gamma) != 6):
+            raise Exception("Lenght of list, alpha should be 6")
+        fCCxm,fCCxp,fCCym,fCCyp,fCCzm,fCCzp = mesh.cellBoundaryInd
+        fxm,fxp,fym,fyp,fzm,fzp = mesh.faceBoundaryInd
+        nBC = fCCxm.sum()+fCCxp.sum()+fCCxm.sum()+fCCxp.sum()
+        h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp]
+        h_ym, h_yp = mesh.gridCC[fCCym], mesh.gridCC[fCCyp]
+        h_zm, h_zp = mesh.gridCC[fCCzm], mesh.gridCC[fCCzp]
+
+        alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0]
+        alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1]
+        alpha_ym, beta_ym, gamma_ym = alpha[2], beta[2], gamma[2]
+        alpha_yp, beta_yp, gamma_yp = alpha[3], beta[3], gamma[3]
+        alpha_zm, beta_zm, gamma_zm = alpha[2], beta[2], gamma[2]
+        alpha_zp, beta_zp, gamma_zp = alpha[3], beta[3], gamma[3]
+
+        h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0]
+        h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1]
+        h_zm, h_zp = mesh.gridCC[fCCzm,2], mesh.gridCC[fCCzp,2]
+
+        a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm)
+        b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
+        a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp)
+        b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp)
+
+        a_ym = gamma_ym/(0.5*alpha_ym-beta_ym/h_ym)
+        b_ym = (0.5*alpha_ym+beta_ym/h_ym)/(0.5*alpha_ym-beta_ym/h_ym)
+        a_yp = gamma_yp/(0.5*alpha_yp-beta_yp/h_yp)
+        b_yp = (0.5*alpha_yp+beta_yp/h_yp)/(0.5*alpha_yp-beta_yp/h_yp)
+
+        a_zm = gamma_zm/(0.5*alpha_zm-beta_zm/h_zm)
+        b_zm = (0.5*alpha_zm+beta_zm/h_zm)/(0.5*alpha_zm-beta_zm/h_zm)
+        a_zp = gamma_zp/(0.5*alpha_zp-beta_zp/h_zp)
+        b_zp = (0.5*alpha_zp+beta_zp/h_zp)/(0.5*alpha_zp-beta_zp/h_zp)
+
+        xBC_xm = 0.5*a_xm
+        xBC_xp = 0.5*a_xp/b_xp
+        yBC_xm = 0.5*(1.-b_xm)
+        yBC_xp = 0.5*(1.-1./b_xp)
+        xBC_ym = 0.5*a_ym
+        xBC_yp = 0.5*a_yp/b_yp
+        yBC_ym = 0.5*(1.-b_ym)
+        yBC_yp = 0.5*(1.-1./b_yp)
+        xBC_zm = 0.5*a_zm
+        xBC_zp = 0.5*a_zp/b_zp
+        yBC_zm = 0.5*(1.-b_zm)
+        yBC_zp = 0.5*(1.-1./b_zp)
+
+        sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm], np.arange(mesh.nFx)[fxp]])
+        sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym], np.arange(mesh.nFy)[fyp]])
+        sortindsfz = np.argsort(np.r_[np.arange(mesh.nFz)[fzm], np.arange(mesh.nFz)[fzp]])
+
+        xBC_x = np.r_[xBC_xm, xBC_xp][sortindsfx]
+        xBC_y = np.r_[xBC_ym, xBC_yp][sortindsfy]
+        xBC_z = np.r_[xBC_zm, xBC_zp][sortindsfz]
+
+        yBC_x = np.r_[yBC_xm, yBC_xp][sortindsfx]
+        yBC_y = np.r_[yBC_ym, yBC_yp][sortindsfy]
+        yBC_z = np.r_[yBC_zm, yBC_zp][sortindsfz]
+
+        xBC = np.r_[xBC_x, xBC_y, xBC_z]
+        yBC = np.r_[yBC_x, yBC_y, yBC_z]
+
+    return xBC, yBC

SimPEG/EM/Static/DC/ProblemDC.py

@@ -10,9 +10,14 @@ class BaseDCProblem(BaseEMProblem):
 
     surveyPair = Survey
     fieldsPair = Fields
+    Ainv = None
 
     def fields(self, m):
         self.curModel = m
+
+        if not self.Ainv == None:
+            self.Ainv.clean()
+
         f = self.fieldsPair(self.mesh, self.survey)
         A = self.getA()
         self.Ainv = self.Solver(A, **self.solverOpts)
@@ -32,13 +37,12 @@ class BaseDCProblem(BaseEMProblem):
         Jv = self.dataPair(self.survey) #same size as the data
 
         A = self.getA()
-        Ainv = self.Solver(A, **self.solverOpts)
 
         for src in self.survey.srcList:
             u_src = f[src, self._solutionType] # solution vector
             dA_dm_v = self.getADeriv(u_src, v)
             dRHS_dm_v = self.getRHSDeriv(src, v)
-            du_dm_v = Ainv * ( - dA_dm_v + dRHS_dm_v )
+            du_dm_v = self.Ainv * ( - dA_dm_v + dRHS_dm_v )
 
             for rx in src.rxList:
                 df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
@@ -57,8 +61,8 @@ class BaseDCProblem(BaseEMProblem):
             v = self.dataPair(self.survey, v)
 
         Jtv = np.zeros(m.size)
-        AT = self.getA().T
-        ATinv = self.Solver(AT, **self.solverOpts)
+        AT = self.getA()
+
 
         for src in self.survey.srcList:
             u_src = f[src, self._solutionType]
@@ -67,7 +71,7 @@ class BaseDCProblem(BaseEMProblem):
                 df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
                 df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True)
 
-                ATinvdf_duT = ATinv * df_duT
+                ATinvdf_duT = self.Ainv * df_duT
 
                 dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True)
                 dRHS_dmT = self.getRHSDeriv(src, ATinvdf_duT, adjoint=True)
@@ -128,13 +132,18 @@ class Problem3D_N(BaseDCProblem):
         #     return V.T * A
         return A
 
-    def getADeriv(self, u, v, adoint=False):
+    def getADeriv(self, u, v, adjoint=False):
         """
 
         Product of the derivative of our system matrix with respect to the model and a vector
 
         """
-        return Div*self.MfRhoIDeriv(Div.T*u)
+        MeSigma = self.MeSigma
+        Grad = self.mesh.nodalGrad
+        if not adjoint:
+            return Grad.T*(self.MeSigmaDeriv(Grad*u)*v)
+        elif adjoint:
+            return self.MeSigmaDeriv(Grad*u).T * (Grad*v)
 
 
     def getRHS(self):
@@ -196,7 +205,7 @@ class Problem3D_CC(BaseDCProblem):
         if adjoint:
             # if self._makeASymmetric is True:
             #     v = V * v
-            return( MfRhoIDeriv( D.T * u ).T) * ( D.T * v)
+            return(MfRhoIDeriv( D.T * u ).T) * ( D.T * v)
 
         # I think we should deprecate this for DC problem.
         # if self._makeASymmetric is True:

SimPEG/EM/Static/DC/SrcDC.py

@@ -31,9 +31,9 @@ class Dipole(BaseSrc):
             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
+            inds = closestPoints(prob.mesh, self.loc)
+            q = np.zeros(prob.mesh.nN)
+            q[inds] = self.current * np.r_[1., -1.]
         return q
 
     # def bc_contribution
@@ -52,9 +52,9 @@ class Pole(BaseSrc):
             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
+            inds = closestPoints(prob.mesh, self.loc)
+            q = np.zeros(prob.mesh.nN)
+            q[inds] = self.current * np.r_[1.]
         return q
 
     # def bc_contribution

SimPEG/EM/Static/DC/__init__.py

@@ -1,4 +1,4 @@
-from ProblemDC import Problem3D_CC
+from ProblemDC import Problem3D_CC, Problem3D_N
 from SurveyDC import Survey
 import SrcDC as Src #Pole
 import RxDC as Rx

tests/em/static/test_DC.py

@@ -3,7 +3,7 @@ from SimPEG import *
 import SimPEG.EM.Static.DC as DC
 
 
-class DCProblemTests(unittest.TestCase):
+class DCProblemTestsCC(unittest.TestCase):
 
     def setUp(self):
 
@@ -63,5 +63,64 @@ class DCProblemTests(unittest.TestCase):
         passed = Tests.checkDerivative(derChk, self.m0, plotIt=False)
         self.assertTrue(passed)
 
+class DCProblemTestsN(unittest.TestCase):
+
+    def setUp(self):
+
+        aSpacing=2.5
+        nElecs=10
+
+        surveySize = nElecs*aSpacing - aSpacing
+        cs = surveySize/nElecs/4
+
+        mesh = Mesh.TensorMesh([
+                [(cs,10, -1.3),(cs,surveySize/cs),(cs,10, 1.3)],
+                [(cs,3, -1.3),(cs,3,1.3)],
+        #         [(cs,5, -1.3),(cs,10)]
+            ],'CN')
+
+        srcList = DC.Utils.WennerSrcList(nElecs, aSpacing, in2D=True)
+        survey = DC.Survey(srcList)
+        problem = DC.Problem3D_N(mesh, mapping=[('rho', Maps.IdentityMap(mesh))])
+        problem.pair(survey)
+
+        mSynth = np.ones(mesh.nC)
+        survey.makeSyntheticData(mSynth)
+
+        # Now set up the problem to do some minimization
+        dmis = DataMisfit.l2_DataMisfit(survey)
+        reg = Regularization.Tikhonov(mesh)
+        opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
+        invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
+        inv = Inversion.BaseInversion(invProb)
+
+        self.inv = inv
+        self.reg = reg
+        self.p =     problem
+        self.mesh = mesh
+        self.m0 = mSynth
+        self.survey = survey
+        self.dmis = dmis
+
+    def test_misfit(self):
+        derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
+        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False)
+        self.assertTrue(passed)
+
+    def test_adjoint(self):
+        # Adjoint Test
+        u = np.random.rand(self.mesh.nC*self.survey.nSrc)
+        v = np.random.rand(self.mesh.nC)
+        w = np.random.rand(self.survey.dobs.shape[0])
+        wtJv = w.dot(self.p.Jvec(self.m0, v))
+        vtJtw = v.dot(self.p.Jtvec(self.m0, w))
+        passed = np.abs(wtJv - vtJtw) < 1e-10
+        print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
+        self.assertTrue(passed)
+
+    def test_dataObj(self):
+        derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
+        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False)
+        self.assertTrue(passed)
 if __name__ == '__main__':
     unittest.main()