Merge branch 'em/dev' into em/primsec

# Conflicts:
#	SimPEG/EM/FDEM/SrcFDEM.py
#	tests/em/fdem/inverse/adjoint/test_FDEM_adjointEB.py

SimPEG/EM/Utils/EMUtils.py

@@ -13,37 +13,4 @@ def k(freq, sigma, mu=mu_0, eps=epsilon_0):
     beta = w * np.sqrt( mu*eps/2 * ( np.sqrt(1. + (sigma / (eps*w))**2 ) - 1) ) 
     return alp - 1j*beta
 
-# Constitutive relations
-def e_from_j(prob,j):
-    eqLocs = prob._eqLocs
-    if eqLocs is 'FE':
-        MSigmaI = prob.MeSigmaI
-    elif eqLocs is 'EF':
-        MSigmaI = prob.MfRho
-    return MSigmaI*j
-
-def j_from_e(prob,e):
-    eqLocs = prob._eqLocs
-    if eqLocs is 'FE':
-        MSigma = prob.MeSigma
-    elif eqLocs is 'EF':
-        MSigma = prob.MfRhoI
-    return MSigma*e
-
-def b_from_h(prob,h):
-    eqLocs = prob._eqLocs
-    if eqLocs is 'FE':
-        MMu = prob.MfMuiI
-    elif eqLocs is 'EF':
-        MMu = prob.MeMu
-    return MMu*h
-
-def h_from_b(prob,b):
-    eqLocs = prob._eqLocs
-    if eqLocs is 'FE':
-        MMuI = prob.MfMui
-    elif eqLocs is 'EF':
-        MMuI = prob.MeMuI
-    return MMuI*b 
-
 

SimPEG/EM/Utils/__init__.py

@@ -1,5 +1,2 @@
-# import Sources
-# import Ana
-# import Solver
-from EMUtils import omega, e_from_j, j_from_e, b_from_h, h_from_b
+from EMUtils import omega, k
 from AnalyticUtils import MagneticDipoleFields, MagneticDipoleVectorPotential, MagneticLoopVectorPotential

SimPEG/EM/Utils/testingUtils.py

@@ -4,19 +4,28 @@ from SimPEG import EM
 import sys
 from scipy.constants import mu_0
 
-def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False):
-    cs = 5.
-    ncx, ncy, ncz = 6, 6, 6
-    npad = 3
+FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
+CONDUCTIVITY = 1e1
+MU = mu_0
+freq = 5e-1
+
+
+def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False):
+    cs = 10.
+    ncx, ncy, ncz = 0, 0, 0
+    npad = 8
     hx = [(cs,npad,-1.3), (cs,ncx), (cs,npad,1.3)]
     hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
     hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)]
     mesh = Mesh.TensorMesh([hx,hy,hz],['C','C','C'])
 
-    mapping = Maps.ExpMap(mesh)
+    if useMu is True:
+        mapping = [('sigma', Maps.ExpMap(mesh)), ('mu', Maps.IdentityMap(mesh))] 
+    else:
+        mapping = Maps.ExpMap(mesh)
 
-    x = np.array([np.linspace(-30,-15,3),np.linspace(15,30,3)]) #don't sample right by the source
-    XYZ = Utils.ndgrid(x,x,np.r_[0.])
+    x = np.array([np.linspace(-5.*cs,-2.*cs,3),np.linspace(5.*cs,2.*cs,3)]) + cs/4. #don't sample right by the source, slightly off alignment from either staggered grid
+    XYZ = Utils.ndgrid(x,x,np.linspace(-2.*cs,2.*cs,5))
     Rx0 = EM.FDEM.Rx(XYZ, comp)
 
     Src = []
@@ -32,15 +41,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))
         elif SrcType is 'PrimSecSigma':
             if fdemType is 'e' or fdemType is 'b':
@@ -77,6 +86,48 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False):
         from pymatsolver import MumpsSolver
         prb.Solver = MumpsSolver
     except ImportError, e:
-        pass
+        prb.Solver = SolverLU
 
-    return prb
+    return prb
+
+def crossCheckTest(SrcList, fdemType1, fdemType2, comp, addrandoms = False, useMu=False, TOL=1e-5, verbose=False):
+
+    l2norm = lambda r: np.sqrt(r.dot(r))
+
+    prb1 = getFDEMProblem(fdemType1, comp, SrcList, freq, useMu, verbose)
+    mesh = prb1.mesh
+    print 'Cross Checking Forward: %s, %s formulations - %s' % (fdemType1, fdemType2, comp)
+    
+    logsig = np.log(np.ones(mesh.nC)*CONDUCTIVITY)
+    mu = np.ones(mesh.nC)*MU
+
+    if addrandoms is True:
+        logsig  += np.random.randn(mesh.nC)*np.log(CONDUCTIVITY)*1e-1
+        mu += np.random.randn(mesh.nC)*MU*1e-1
+
+    if useMu is True:
+        m = np.r_[logsig, mu]
+    else:
+        m = logsig
+
+    survey1 = prb1.survey
+    d1 = survey1.dpred(m)
+
+    if verbose:
+        print '  Problem 1 solved'
+
+
+    prb2 = getFDEMProblem(fdemType2, comp, SrcList, freq, useMu, verbose)
+
+    survey2 = prb2.survey
+    d2 = survey2.dpred(m)
+
+    if verbose:
+        print '  Problem 2 solved'
+
+    r = d2-d1
+    l2r = l2norm(r)
+
+    tol = np.max([TOL*(10**int(np.log10(0.5* (l2norm(d1) + l2norm(d2)) ))),FLR])
+    print l2norm(d1), l2norm(d2),  l2r , tol, l2r < tol
+    return l2r < tol