From 9e69455d6e1116e67fe3ad889e1693effb8b02ac Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 8 Dec 2015 23:35:06 -0800 Subject: [PATCH] b from h and h from b --- SimPEG/EM/FDEM/FieldsFDEM.py | 108 ++++++++++++++++++++- SimPEG/EM/Utils/testingUtils.py | 8 +- tests/em/fdem/forward/test_FDEM_forward.py | 80 +++++++++++++++ 3 files changed, 187 insertions(+), 9 deletions(-) diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py index 22838e17..fb05b8c6 100644 --- a/SimPEG/EM/FDEM/FieldsFDEM.py +++ b/SimPEG/EM/FDEM/FieldsFDEM.py @@ -34,6 +34,7 @@ class Fields_e(Fields): self._aveF2CCV = self.survey.prob.mesh.aveF2CCV self._sigma = self.survey.prob.curModel.sigma self._sigmaDeriv = self.survey.prob.curModel.sigmaDeriv + self._mui = self.survey.prob.curModel.mui self._nC = self.survey.prob.mesh.nC def _GLoc(self,fieldType): @@ -108,12 +109,25 @@ class Fields_e(Fields): sigma = self._sigma aveE2CCV = self._aveE2CCV n = int(aveE2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy + # Sigma = sdiag(np.kron(np.ones(n), sigma)) + Sigma = self.prob.MeSigma VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol)) - Sigma = sdiag(np.kron(np.ones(n), sigma)) e = self._e(eSolution, srcList) - return Sigma * (aveE2CCV * e) + return VI * (aveE2CCV * (Sigma *e) ) + + def _h(self, eolution, srcList): + b = self._b(eSolution, srcList) + Mui = self.survey.prob.MfMui + aveF2CCV = self._aveF2CCV + n = int(aveF2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy + # Mui = sdiag(sp.kron(np.ones(n), mui)) + VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol)) + + + return VI * (aveF2CCV * (Mui * b)) + def _jDeriv_u(self, src, v, adjoint=False): raise NotImplementedError @@ -149,6 +163,8 @@ class Fields_b(Fields): 'e' : ['bSolution','E','_e'], 'ePrimary' : ['bSolution','E','_ePrimary'], 'eSecondary' : ['bSolution','E','_eSecondary'], + 'j' : ['bSolution','C','_j'], + 'h' : ['bSolution','C','_h'], } def __init__(self,mesh,survey,**kwargs): @@ -161,6 +177,13 @@ class Fields_b(Fields): self._MfMui = self.survey.prob.MfMui self._MeSigmaIDeriv = self.survey.prob.MeSigmaIDeriv self._Me = self.survey.prob.Me + self._aveF2CCV = self.survey.prob.mesh.aveF2CCV + self._aveE2CCV = self.survey.prob.mesh.aveE2CCV + self._sigma = self.survey.prob.curModel.sigma + self._mui = self.survey.prob.curModel.mui + self._nC = self.survey.prob.mesh.nC + + def _GLoc(self,fieldType): if fieldType == 'e': @@ -245,6 +268,29 @@ class Fields_b(Fields): # assuming primary doesn't depend on model return self._eSecondaryDeriv_m(src, v, adjoint) + def _j(self, bSolution, srcList): + sigma = self._sigma + aveE2CCV = self._aveE2CCV + n = int(aveE2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy + # Sigma = sdiag(np.kron(np.ones(n), sigma)) + Sigma = self.prob.MeSigma + VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol)) + + e = self._e(bSolution, srcList) + + return VI * (aveE2CCV * (Sigma *e) ) + + def _h(self, bSolution, srcList): + b = self._b(bSolution, srcList) + Mui = self.survey.prob.MfMui + aveF2CCV = self._aveF2CCV + n = int(aveF2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy + # Mui = sdiag(sp.kron(np.ones(n), mui)) + VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol)) + + + return VI * (aveF2CCV * (Mui * b)) + class Fields_j(Fields): knownFields = {'jSolution':'F'} @@ -256,6 +302,7 @@ class Fields_j(Fields): 'hPrimary' : ['jSolution','E','_hPrimary'], 'hSecondary' : ['jSolution','E','_hSecondary'], 'e' : ['jSolution','C','_e'], + 'b' : ['jSolution','C','_b'], } def __init__(self,mesh,survey,**kwargs): @@ -269,7 +316,9 @@ class Fields_j(Fields): self._MfRhoDeriv = self.survey.prob.MfRhoDeriv self._Me = self.survey.prob.Me self._rho = self.survey.prob.curModel.rho + self._mu = self.survey.prob.curModel.mui self._aveF2CCV = self.survey.prob.mesh.aveF2CCV + self._aveE2CCV = self.survey.prob.mesh.aveE2CCV self._nC = self.survey.prob.mesh.nC def _GLoc(self,fieldType): @@ -361,11 +410,14 @@ class Fields_j(Fields): rho = self._rho aveF2CCV = self._aveF2CCV n = int(aveF2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy - Rho = sdiag(np.kron(np.ones(n), rho)) + + Rho = self.prob.MfRho + VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol)) + j = self._j(jSolution, srcList) - return 1./self.survey.prob.mesh.dim * Rho * (aveF2CCV * j) + return VI * (aveF2CCV * (Rho * j)) def _eDeriv_u(self, src, v, adjoint=False): raise NotImplementedError @@ -373,6 +425,16 @@ class Fields_j(Fields): def _eDeriv_m(self, src, v, adjoint=False): raise NotImplementedError + def _b(self, jSolution, srcList): + h = self._h(jSolution, srcList) + Mu = self.prob.MeMu + aveE2CCV = self._aveE2CCV + n = int(aveE2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy + # Mu = sdiag(sp.kron(np.ones(n), mu)) + VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol)) + + return VI * (aveE2CCV * (Mu * h)) + class Fields_h(Fields): knownFields = {'hSolution':'E'} @@ -382,7 +444,9 @@ class Fields_h(Fields): 'hSecondary' : ['hSolution','E','_hSecondary'], 'j' : ['hSolution','F','_j'], 'jPrimary' : ['hSolution','F','_jPrimary'], - 'jSecondary' : ['hSolution','F','_jSecondary'] + 'jSecondary' : ['hSolution','F','_jSecondary'], + 'e' : ['hSolution','C','_e'], + 'b' : ['hSolution','C','_b'], } def __init__(self,mesh,survey,**kwargs): @@ -393,6 +457,11 @@ class Fields_h(Fields): self._edgeCurl = self.survey.prob.mesh.edgeCurl self._MeMuI = self.survey.prob.MeMuI self._MfRho = self.survey.prob.MfRho + self._rho = self.survey.prob.curModel.rho + self._mu = self.survey.prob.curModel.mui + self._aveF2CCV = self.survey.prob.mesh.aveF2CCV + self._aveE2CCV = self.survey.prob.mesh.aveE2CCV + self._nC = self.survey.prob.mesh.nC def _GLoc(self,fieldType): if fieldType == 'h': @@ -458,3 +527,32 @@ class Fields_h(Fields): def _jDeriv_m(self, src, v, adjoint=False): # assuming the primary does not depend on the model return self._jSecondaryDeriv_m(src,v,adjoint) + + def _e(self, hSolution, srcList): + rho = self._rho + aveF2CCV = self._aveF2CCV + n = int(aveF2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy + + Rho = self.prob.MfRho + VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol)) + + + j = self._j(hSolution, srcList) + + return VI * (aveF2CCV * (Rho * j)) + + def _eDeriv_u(self, src, v, adjoint=False): + raise NotImplementedError + + def _eDeriv_m(self, src, v, adjoint=False): + raise NotImplementedError + + def _b(self, hSolution, srcList): + h = self._h(hSolution, srcList) + Mu = self.prob.MeMu + aveE2CCV = self._aveE2CCV + n = int(aveE2CCV.shape[0] / self._nC) #TODO: This is a bit sloppy + # Mu = sdiag(sp.kron(np.ones(n), mu)) + VI = sdiag(1./np.kron(np.ones(n), self.prob.mesh.vol)) + + return VI * (aveE2CCV * (Mu * h)) diff --git a/SimPEG/EM/Utils/testingUtils.py b/SimPEG/EM/Utils/testingUtils.py index 59971192..9979417d 100644 --- a/SimPEG/EM/Utils/testingUtils.py +++ b/SimPEG/EM/Utils/testingUtils.py @@ -32,15 +32,15 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False): if fdemType is 'e' or fdemType is 'b': S_m = np.zeros(mesh.nF) S_e = np.zeros(mesh.nE) - S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1. - S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1. + S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1e-3 + S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1e-3 Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e)) elif fdemType is 'h' or fdemType is 'j': S_m = np.zeros(mesh.nE) S_e = np.zeros(mesh.nF) - S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1. - S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1. + S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1e-3 + S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1e-3 Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e)) if verbose: diff --git a/tests/em/fdem/forward/test_FDEM_forward.py b/tests/em/fdem/forward/test_FDEM_forward.py index 21c0d114..da22c950 100644 --- a/tests/em/fdem/forward/test_FDEM_forward.py +++ b/tests/em/fdem/forward/test_FDEM_forward.py @@ -8,6 +8,7 @@ from SimPEG.EM.Utils.testingUtils import getFDEMProblem testEB = True testHJ = True testEJ = True +testBH = True verbose = False TOLEBHJ = 1e-5 @@ -145,5 +146,84 @@ class FDEM_CrossCheck(unittest.TestCase): def test_EJ_CrossCheck_ezi_Jform(self): self.assertTrue(crossCheckTest('e', 'j', 'ezi', TOL=TOLEJHB)) + def test_EJ_CrossCheck_bxr_Jform(self): + self.assertTrue(crossCheckTest('e', 'j', 'bxr', TOL=TOLEJHB)) + def test_EJ_CrossCheck_byr_Jform(self): + self.assertTrue(crossCheckTest('e', 'j', 'byr', TOL=TOLEJHB)) + def test_EJ_CrossCheck_bzr_Jform(self): + self.assertTrue(crossCheckTest('e', 'j', 'bzr', TOL=TOLEJHB)) + def test_EJ_CrossCheck_bxi_Jform(self): + self.assertTrue(crossCheckTest('e', 'j', 'bxi', TOL=TOLEJHB)) + def test_EJ_CrossCheck_byi_Jform(self): + self.assertTrue(crossCheckTest('e', 'j', 'byi', TOL=TOLEJHB)) + def test_EJ_CrossCheck_bzi_Jform(self): + self.assertTrue(crossCheckTest('e', 'j', 'bzi', TOL=TOLEJHB)) + + def test_EJ_CrossCheck_hxr_Jform(self): + self.assertTrue(crossCheckTest('e', 'j', 'hxr', TOL=TOLEJHB)) + def test_EJ_CrossCheck_hyr_Jform(self): + self.assertTrue(crossCheckTest('e', 'j', 'hyr', TOL=TOLEJHB)) + def test_EJ_CrossCheck_hzr_Jform(self): + self.assertTrue(crossCheckTest('e', 'j', 'hzr', TOL=TOLEJHB)) + def test_EJ_CrossCheck_hxi_Jform(self): + self.assertTrue(crossCheckTest('e', 'j', 'hxi', TOL=TOLEJHB)) + def test_EJ_CrossCheck_hyi_Jform(self): + self.assertTrue(crossCheckTest('e', 'j', 'hyi', TOL=TOLEJHB)) + def test_EJ_CrossCheck_hzi_Jform(self): + self.assertTrue(crossCheckTest('e', 'j', 'hzi', TOL=TOLEJHB)) + + if testBH: + def test_BH_CrossCheck_jxr_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'jxr', TOL=TOLEJHB)) + def test_BH_CrossCheck_jyr_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'jyr', TOL=TOLEJHB)) + def test_BH_CrossCheck_jzr_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'jzr', TOL=TOLEJHB)) + def test_BH_CrossCheck_jxi_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'jxi', TOL=TOLEJHB)) + def test_BH_CrossCheck_jyi_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'jyi', TOL=TOLEJHB)) + def test_BH_CrossCheck_jzi_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'jzi', TOL=TOLEJHB)) + + def test_BH_CrossCheck_exr_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'exr', TOL=TOLEJHB)) + def test_BH_CrossCheck_eyr_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'eyr', TOL=TOLEJHB)) + def test_BH_CrossCheck_ezr_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'ezr', TOL=TOLEJHB)) + def test_BH_CrossCheck_exi_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'exi', TOL=TOLEJHB)) + def test_BH_CrossCheck_eyi_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'eyi', TOL=TOLEJHB)) + def test_BH_CrossCheck_ezi_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'ezi', TOL=TOLEJHB)) + + def test_BH_CrossCheck_bxr_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'bxr', TOL=TOLEJHB)) + def test_BH_CrossCheck_byr_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'byr', TOL=TOLEJHB)) + def test_BH_CrossCheck_bzr_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'bzr', TOL=TOLEJHB)) + def test_BH_CrossCheck_bxi_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'bxi', TOL=TOLEJHB)) + def test_BH_CrossCheck_byi_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'byi', TOL=TOLEJHB)) + def test_BH_CrossCheck_bzi_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'bzi', TOL=TOLEJHB)) + + def test_BH_CrossCheck_hxr_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'hxr', TOL=TOLEJHB)) + def test_BH_CrossCheck_hyr_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'hyr', TOL=TOLEJHB)) + def test_BH_CrossCheck_hzr_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'hzr', TOL=TOLEJHB)) + def test_BH_CrossCheck_hxi_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'hxi', TOL=TOLEJHB)) + def test_BH_CrossCheck_hyi_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'hyi', TOL=TOLEJHB)) + def test_BH_CrossCheck_hzi_Jform(self): + self.assertTrue(crossCheckTest('b', 'h', 'hzi', TOL=TOLEJHB)) + if __name__ == '__main__': unittest.main() \ No newline at end of file