Start of JH implementation, solving for h. Derivatives for ProblemFDEM_h NOT working yet

simpegEM/FDEM/FDEM.py

@@ -338,6 +338,7 @@ class ProblemFDEM_b(BaseFDEMProblem):
         raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
 
 
+
 ##########################################################################################
 ################################ H-J Formulation #########################################
 ##########################################################################################
@@ -424,6 +425,7 @@ class ProblemFDEM_j(BaseFDEMProblem):
             rhs[i] = np.concatenate((SRCx, SRCy, SRCz))
 
         a = np.concatenate(rhs).reshape((self.mesh.nF, len(Txs)), order='F')
+        a = Utils.mkvc(a)
         MeMui = self.MeMui
         C = self.mesh.edgeCurl
 
@@ -463,3 +465,153 @@ class ProblemFDEM_j(BaseFDEMProblem):
             else:
                 return -(1./(1j*omega(freq))) * dsig_dm.T * ( dsigi_dsig.T * ( dMf_dsigi.T * ( C * ( MeMui.T * v ) ) ) )       
         raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
+
+
+
+
+# Solving for h! - NOTE: WE ARE GOING TO NEED dRHS / dm !
+class ProblemFDEM_h(BaseFDEMProblem):
+    """
+        Using the H-J formulation of Maxwell's equations
+
+        .. math::
+            \\nabla \\times \\sigma^{-1} \\vec{J} + i\\omega\\mu\\vec{H} = 0
+            \\nabla \\times \\vec{H} - \\vec{J} = \\vec{J_s}
+
+        Since \(\\vec{J}\) is a flux and \(\\vec{H}\) is a field, we discretize \(\\vec{J}\) on faces and \(\\vec{H}\) on edges.
+
+        For this implementation, we solve for J using \( \\vec{J} =  \\nabla \\times \\vec{H} - \\vec{J_s} \)
+
+        .. math::
+            \\nabla \\times \\sigma^{-1} \\nabla \\times \\vec{H} + i\\omega\\mu\\vec{H} = \\nabla \\times \\sigma^{-1} \\vec{J_s}
+
+        .. note::
+            This implementation does not yet work with full anisotropy!!
+
+    """
+
+    solType = 'h'
+    storeTheseFields = ['j','h']
+
+    def __init__(self, model, **kwargs):
+        BaseFDEMProblem.__init__(self, model, **kwargs)
+
+    def getA(self, freq):
+        """
+            :param float freq: Frequency
+            :rtype: scipy.sparse.csr_matrix
+            :return: A
+        """
+
+        MeMu = self.MeMu
+        MfSigi = self.MfSigmai
+        C = self.mesh.edgeCurl
+
+        return C.T * MfSigi * C + 1j*omega(freq)*MeMu
+
+    def getADeriv(self, freq, u, v, adjoint=False):
+
+        MeMu = self.MeMu
+        C = self.mesh.edgeCurl
+        sig = self.curModel.transform
+        sigi = 1/sig
+        dsig_dm = self.curModel.transformDeriv
+        dsigi_dsig = -Utils.sdiag(sigi)**2
+
+        dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(C*u)
+
+        if adjoint:
+            return (dsig_dm.T * (dsigi_dsig.T * (dMf_dsigi.T * (C * v))))
+        return  (C.T * (dMf_dsigi * (dsigi_dsig * (dsig_dm * v))))
+
+
+    def getjs(self,freq):
+        """
+            :param float freq: Frequency
+            :rtype: numpy.ndarray (nE, nTx)
+            :return: j_s
+        """
+        Txs = self.survey.getTransmitters(freq)
+        rhs = range(len(Txs))
+        for i, tx in enumerate(Txs):
+            if tx.txType == 'VMD':
+                src = Sources.MagneticDipoleVectorPotential
+                SRCx = src(tx.loc, self.mesh.gridFx, 'x')
+                SRCy = src(tx.loc, self.mesh.gridFy, 'y')
+                SRCz = src(tx.loc, self.mesh.gridFz, 'z')
+
+            elif tx.txType == 'CircularLoop':
+                src = Sources.MagneticLoopVectorPotential                
+                SRCx = src(tx.loc, self.mesh.gridFx, 'x', tx.radius)
+                SRCy = src(tx.loc, self.mesh.gridFy, 'y', tx.radius)
+                SRCz = src(tx.loc, self.mesh.gridFz, 'z', tx.radius)
+            else:
+                raise NotImplemented('%s txType is not implemented' % tx.txType)
+            rhs[i] = np.concatenate((SRCx, SRCy, SRCz))
+
+        a = np.concatenate(rhs).reshape((self.mesh.nF, len(Txs)), order='F')
+        a = Utils.mkvc(a)
+
+        MeMui = self.MeMui
+        C = self.mesh.edgeCurl
+
+        return C*MeMui*C.T*a
+
+    def getRHS(self, freq):
+        """
+            :param float freq: Frequency
+            :rtype: numpy.ndarray (nE, nTx)
+            :return: RHS
+        """
+        MfSigi = self.MfSigmai
+        C = self.mesh.edgeCurl
+        j_s = self.getjs(freq)
+        return C.T*MfSigi*j_s
+
+    def getRHSDeriv(self, freq, v, adjoint=False):
+        """
+            :param float freq: Frequency
+            :rtype: numpy.ndarray (nE, nTx)
+            :return: RHSDeriv
+        """
+        C = self.mesh.edgeCurl
+        sig = self.curModel.transform
+        sigi = 1/sig
+        j_s = self.getjs(freq)
+        dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(j_s)
+        dsig_dm = self.curModel.transformDeriv
+        dsigi_dsig = -Utils.sdiag(sigi)**2 # only works for diagonal matrices 
+
+        if adjoint:
+            return dsig_dm.T * dsigi_dsig.T * dMf_dsigi.T * C * v
+        return C.T * dMf_dsigi * dsigi_dsig * dsig_dm * v
+        
+    def calcFields(self, sol, freq, fieldType, adjoint=False):
+        h = sol
+        if fieldType == 'j':
+            NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
+            C = self.mesh.edgeCurl
+            j_s = self.getjs(freq)
+            if adjoint: 
+                NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
+                # return C.T*h # TODO: THIS IS WRONG 
+            NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
+            return C*h # - j_s
+        elif fieldType == 'h':
+            return h
+        raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
+
+    def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False):
+        h = sol
+        if fieldType == 'j':
+            NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
+            C = self.mesh.edgeCurl
+            drhs = self.getRHSDeriv(freq,v,adjoint)
+            if adjoint:
+                NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
+                # return -C.T * drhs
+            NotImplementedError('fieldType "%s" is not implemented.' % fieldType)
+            # return -C * drhs
+        elif fieldType == 'h':
+            return -self.getRHSDeriv(freq,v,adjoint)
+        raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType)

simpegEM/FDEM/__init__.py

@@ -1,2 +1,2 @@
 from SurveyFDEM import *
-from FDEM import ProblemFDEM_e, ProblemFDEM_b, ProblemFDEM_j
+from FDEM import ProblemFDEM_e, ProblemFDEM_b, ProblemFDEM_j, ProblemFDEM_h

simpegEM/Tests/test_FDEM.py

@@ -8,6 +8,7 @@ TOL = 1e-4
 FLR = 1e-15 # "zero", so if residual below this --> pass regardless of order
 CONDUCTIVITY = 1e1
 MU = mu_0
+freq = 1
 addrandoms = True # important to addrandoms if testing HJ formulation with VMD source! (or else jz ~ 0)
 
 def getProblem(fdemType, comp):
@@ -24,7 +25,7 @@ def getProblem(fdemType, comp):
     x = np.linspace(-30,30,6)
     XYZ = Utils.ndgrid(x,x,np.r_[0])
     Rx0 = EM.FDEM.RxFDEM(XYZ, comp)
-    Tx0 = EM.FDEM.TxFDEM(np.r_[0.,0.,0.], 'VMD', 1, [Rx0])
+    Tx0 = EM.FDEM.TxFDEM(np.r_[0.,0.,0.], 'VMD', freq, [Rx0])
 
     survey = EM.FDEM.SurveyFDEM([Tx0])
 
@@ -36,6 +37,8 @@ def getProblem(fdemType, comp):
         prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping)
     elif fdemType == 'j':
         prb = EM.FDEM.ProblemFDEM_j(mesh, mapping=mapping)
+    elif fdemType == 'h':
+        prb = EM.FDEM.ProblemFDEM_h(mesh, mapping=mapping)
     else:
         raise NotImplementedError()
     prb.pair(survey)
@@ -83,7 +86,7 @@ def derivTest(fdemType, comp):
         x0  = x0 + np.random.randn(prb.mesh.nC)*CONDUCTIVITY*1e-1 
         mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-1
 
-    prb.mu = mu 
+    prb.mu = mu
     survey = prb.survey
     def fun(x):
         return survey.dpred(x), lambda x: prb.Jvec(x0, x)
@@ -221,6 +224,32 @@ class FDEM_DerivTests(unittest.TestCase):
     def test_Jvec_hzi_Jform(self):
         self.assertTrue(derivTest('j', 'hzi'))
 
+    # def test_Jvec_hxr_Hform(self):
+    #     self.assertTrue(derivTest('h', 'hxr'))
+    # def test_Jvec_hyr_Hform(self):
+    #     self.assertTrue(derivTest('h', 'hyr'))
+    # def test_Jvec_hzr_Hform(self):
+    #     self.assertTrue(derivTest('h', 'hzr'))
+    # def test_Jvec_hxi_Hform(self):
+    #     self.assertTrue(derivTest('h', 'hxi'))
+    # def test_Jvec_hyi_Hform(self):
+    #     self.assertTrue(derivTest('h', 'hyi'))
+    # def test_Jvec_hzi_Hform(self):
+    #     self.assertTrue(derivTest('h', 'hzi'))
+
+    # def test_Jvec_hxr_Hform(self):
+    #     self.assertTrue(derivTest('h', 'jxr'))
+    # def test_Jvec_hyr_Hform(self):
+    #     self.assertTrue(derivTest('h', 'jyr'))
+    # def test_Jvec_hzr_Hform(self):
+    #     self.assertTrue(derivTest('h', 'jzr'))
+    # def test_Jvec_hxi_Hform(self):
+    #     self.assertTrue(derivTest('h', 'jxi'))
+    # def test_Jvec_hyi_Hform(self):
+    #     self.assertTrue(derivTest('h', 'jyi'))
+    # def test_Jvec_hzi_Hform(self):
+    #     self.assertTrue(derivTest('h', 'jzi'))
+
     def test_Jtvec_adjointTest_jxr_Jform(self):
         self.assertTrue(adjointTest('j', 'jxr'))
     def test_Jtvec_adjointTest_jyr_Jform(self):
@@ -247,6 +276,33 @@ class FDEM_DerivTests(unittest.TestCase):
     def test_Jtvec_adjointTest_hzi_Jform(self):
         self.assertTrue(adjointTest('j', 'hzi'))
 
+    # def test_Jtvec_adjointTest_hxr_Hform(self):
+    #     self.assertTrue(adjointTest('h', 'hxr'))
+    # def test_Jtvec_adjointTest_hyr_Hform(self):
+    #     self.assertTrue(adjointTest('h', 'hyr'))
+    # def test_Jtvec_adjointTest_hzr_Hform(self):
+    #     self.assertTrue(adjointTest('h', 'hzr'))
+    # def test_Jtvec_adjointTest_hxi_Hform(self):
+    #     self.assertTrue(adjointTest('h', 'hxi'))
+    # def test_Jtvec_adjointTest_hyi_Hform(self):
+    #     self.assertTrue(adjointTest('h', 'hyi'))
+    # def test_Jtvec_adjointTest_hzi_Hform(self):
+    #     self.assertTrue(adjointTest('h', 'hzi'))
+
+    # def test_Jtvec_adjointTest_hxr_Hform(self):
+    #     self.assertTrue(adjointTest('h', 'jxr'))
+    # def test_Jtvec_adjointTest_hyr_Hform(self):
+    #     self.assertTrue(adjointTest('h', 'jyr'))
+    # def test_Jtvec_adjointTest_hzr_Hform(self):
+    #     self.assertTrue(adjointTest('h', 'jzr'))
+    # def test_Jtvec_adjointTest_hxi_Hform(self):
+    #     self.assertTrue(adjointTest('h', 'jxi'))
+    # def test_Jtvec_adjointTest_hyi_Hform(self):
+    #     self.assertTrue(adjointTest('h', 'jyi'))
+    # def test_Jtvec_adjointTest_hzi_Hform(self):
+    #     self.assertTrue(adjointTest('h', 'jzi'))
+
+
 
 if __name__ == '__main__':
     unittest.main()