removed MT1d example from this pr

if mumps import fails use LU

SimPEG/Examples/EM_FDEM_SusEffects.py

@@ -0,0 +1,154 @@
+from SimPEG import *
+from SimPEG import EM
+from scipy.constants import mu_0
+
+def run(plotIt=True):
+    """
+        EM: FDEM: Effects of susceptibility
+        ===================================
+
+        When airborne freqeuncy domain EM (AFEM) survey is flown over
+        the earth including significantly susceptible bodies (magnetite-rich rocks),
+        negative data is often observed in the real part of the lowest frequency
+        (e.g. Dighem system 900 Hz). This phenomenon mostly based upon magnetization
+        occurs due to a susceptible body when the magnetic field is applied.
+
+        To clarify what is happening in the earth when we are exciting the earth with
+        a loop source in the frequency domain we run three forward modelling:
+
+            - F[:math:`\sigma`, :math:`\mu`]: Anomalous conductivity and susceptibility
+            - F[:math:`\sigma`, :math:`\mu_0`]: Anomalous conductivity
+            - F[:math:`\sigma_{air}`, :math:`\mu_0`]: primary field
+
+        We plot vector magnetic fields in the earth. For secondary fields we provide
+        F[:math:`\sigma`, :math:`\mu`]-F[:math:`\sigma`, :math:`\mu_0`]. Following
+        figure show both real and parts.
+
+    """
+    # Generate Cylindrical mesh
+    cs, ncx, ncz, npad = 5, 25, 24, 20.
+    hx = [(cs,ncx), (cs,npad,1.3)]
+    hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)]
+    mesh = Mesh.CylMesh([hx,1,hz], '00C')
+    sighalf = 1e-3
+    sigma = np.ones(mesh.nC)*1e-8
+    sigmahomo = sigma.copy()
+    mu = np.ones(mesh.nC)*mu_0
+    sigma[mesh.gridCC[:,-1]<0.] = sighalf
+    blkind  = np.logical_and(mesh.gridCC[:,0]<30., (mesh.gridCC[:,2]<0)&(mesh.gridCC[:,2]>-150)&(mesh.gridCC[:,2]<-50))
+    sigma[blkind] = 1e-1
+    mu[blkind] = mu_0*1.1
+    offset = 0.
+    frequency = np.r_[10., 100., 1000.]
+    rx0 = EM.FDEM.Rx(np.array([[8., 0., 30.]]), 'bzr')
+    rx1 = EM.FDEM.Rx(np.array([[8., 0., 30.]]), 'bzi')
+    srcLists = []
+    nfreq = frequency.size
+    for ifreq in range(nfreq):
+        src = EM.FDEM.Src.CircularLoop([rx0, rx1], frequency[ifreq], np.array([[0., 0., 30.]]), radius=5.)
+        srcLists.append(src)
+    survey = EM.FDEM.Survey(srcLists)
+    iMap = Maps.IdentityMap(nP=int(mesh.nC))
+    # Use PhysPropMap
+    maps = [('sigma', iMap), ('mu', iMap)]
+    prob = EM.FDEM.Problem_b(mesh, mapping=maps)
+    prob.Solver = MumpsSolver
+    survey.pair(prob)
+    m = np.r_[sigma, mu]
+    survey0 = EM.FDEM.Survey(srcLists)
+    prob0 = EM.FDEM.Problem_b(mesh, mapping=maps)
+    
+
+    try:
+        from pymatsolver import MumpsSolver
+        prb.Solver = MumpsSolver
+    except ImportError, e:
+            prb.Solver = SolverLU
+
+    survey0.pair(prob0)
+    m = np.r_[sigma, mu]
+    m0 = np.r_[sigma, np.ones(mesh.nC)*mu_0]
+    m00 = np.r_[np.ones(mesh.nC)*1e-8, np.ones(mesh.nC)*mu_0]
+    # Anomalous conductivity and susceptibility
+    F = prob.fields(m)
+    # Only anomalous conductivity
+    F0 = prob.fields(m0)
+    # Primary field
+    F00 = prob.fields(m00)
+
+    if plotIt:
+        import matplotlib.pyplot as plt
+        def vizfields(ifreq=0, primsec="secondary",realimag="real"):
+
+            titles = ["F[$\sigma$, $\mu$]", "F[$\sigma$, $\mu_0$]", "F[$\sigma$, $\mu$]-F[$\sigma$, $\mu_0$]"]
+            actind  = np.logical_and(mesh.gridCC[:,0]<200., (mesh.gridCC[:,2]>-400)&(mesh.gridCC[:,2]<200))
+
+            if primsec=="secondary":
+                bCCprim = (mesh.aveF2CCV*F00[:,'b'][:,ifreq]).reshape(mesh.nC, 2, order='F')
+                bCC = (mesh.aveF2CCV*F[:,'b'][:,ifreq]).reshape(mesh.nC, 2, order='F')-bCCprim
+                bCC0 = (mesh.aveF2CCV*F0[:,'b'][:,ifreq]).reshape(mesh.nC, 2, order='F')-bCCprim
+            elif primsec=="primary":
+                bCC = (mesh.aveF2CCV*F[:,'b'][:,ifreq]).reshape(mesh.nC, 2, order='F')
+                bCC0 = (mesh.aveF2CCV*F0[:,'b'][:,ifreq]).reshape(mesh.nC, 2, order='F')
+
+            XYZ = mesh.gridCC[actind,:]
+            X = XYZ[:,0].reshape((31,43), order='F')
+            Z = XYZ[:,2].reshape((31,43), order='F')
+            bx = bCC[actind,0].reshape((31,43), order='F')
+            bz = bCC[actind,1].reshape((31,43), order='F')
+            bx0 = bCC0[actind,0].reshape((31,43), order='F')
+            bz0 = bCC0[actind,1].reshape((31,43), order='F')
+
+            bxsec = (bCC[actind,0]-bCC0[actind,0]).reshape((31,43), order='F')
+            bzsec = (bCC[actind,1]-bCC0[actind,1]).reshape((31,43), order='F')
+
+            absbreal = np.sqrt(bx.real**2+bz.real**2)
+            absbimag = np.sqrt(bx.imag**2+bz.imag**2)
+            absb0real = np.sqrt(bx0.real**2+bz0.real**2)
+            absb0imag = np.sqrt(bx0.imag**2+bz0.imag**2)
+
+            absbrealsec = np.sqrt(bxsec.real**2+bzsec.real**2)
+            absbimagsec = np.sqrt(bxsec.imag**2+bzsec.imag**2)
+
+            fig = plt.figure(figsize=(15,5))
+            ax1 = plt.subplot(131)
+            ax2 = plt.subplot(132)
+            ax3 = plt.subplot(133)
+            typefield="real"
+            scale=20
+            if realimag=="real":
+                ax1.contourf(X, Z,np.log10(absbreal), 100)
+                ax1.quiver(X, Z,bx.real/absbreal,bz.real/absbreal,scale=scale,width=0.005, alpha = 0.5)
+                ax2.contourf(X, Z,np.log10(absb0real), 100)
+                ax2.quiver(X, Z,bx0.real/absb0real,bz0.real/absb0real,scale=scale,width=0.005, alpha = 0.5)
+                ax3.contourf(X, Z,np.log10(absbrealsec), 100)
+                ax3.quiver(X, Z,bxsec.real/absbrealsec,bzsec.real/absbrealsec,scale=scale,width=0.005, alpha = 0.5)
+            elif realimag=="imag":
+                ax1.contourf(X, Z,np.log10(absbimag), 100)
+                ax1.quiver(X, Z,bx.imag/absbimag,bz.imag/absbimag,scale=scale,width=0.005, alpha = 0.5)
+                ax2.contourf(X, Z,np.log10(absb0imag), 100)
+                ax2.quiver(X, Z,bx0.imag/absb0imag,bz0.imag/absb0imag,scale=scale,width=0.005, alpha = 0.5)
+                ax3.contourf(X, Z,np.log10(absbimagsec), 100)
+                ax3.quiver(X, Z,bxsec.imag/absbimagsec,bzsec.imag/absbimagsec,scale=scale,width=0.005, alpha = 0.5)
+
+            ax = [ax1, ax2, ax3]
+            ax3.text(30, 50, ("Frequency=%5.2f Hz")%(frequency[ifreq]), color="k", fontsize=18)
+            ax2.text(30, 50, primsec, color="k", fontsize=18)
+            ax1.text(30, 50, realimag, color="k", fontsize=18)
+            for i, axtemp in enumerate(ax):
+                axtemp.plot(np.r_[0, 29.75], np.r_[-50, -50], 'w', lw=3)
+
+                axtemp.plot(np.r_[29.5, 29.5], np.r_[-50, -142.5], 'w', lw=3)
+                axtemp.plot(np.r_[0, 29.5], np.r_[-142.5, -142.5], 'w', lw=3)
+                axtemp.plot(np.r_[0, 100.], np.r_[0, 0], 'w', lw=3)
+                axtemp.set_ylim(-200, 100.)
+                axtemp.set_xlim(10, 100.)
+                axtemp.set_title(titles[i])
+            plt.show()
+            return fig, ax
+        fig1, ax1 = vizfields(1, primsec="primary", realimag="real")
+        fig2, ax2 = vizfields(1, primsec="secondary", realimag="real")
+        fig4, ax4 = vizfields(1, primsec="secondary", realimag="imag")
+
+if __name__ == '__main__':
+    run()

SimPEG/Examples/__init__.py

@@ -3,6 +3,7 @@
 ##### AUTOIMPORTS #####
 import EM_FDEM_1D_Inversion
 import EM_FDEM_Analytic_MagDipoleWholespace
+import EM_FDEM_SusEffects
 import EM_TDEM_1D_Inversion
 import FLOW_Richards_1D_Celia1990
 import Forward_BasicDirectCurrent
@@ -14,8 +15,9 @@ import Mesh_QuadTree_Creation
 import Mesh_QuadTree_FaceDiv
 import Mesh_QuadTree_HangingNodes
 import Mesh_Tensor_Creation
+import MT_1D_analytic_nlayer_Earth
 
-__examples__ = ["EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation"]
+__examples__ = ["EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_FDEM_SusEffects", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_analytic_nlayer_Earth"]
 
 ##### AUTOIMPORTS #####
 

SimPEG/Optimization.py

@@ -990,4 +990,18 @@ class ProjectedGNCG(BFGS, Minimize, Remember):
                     cgFlag = 1
                 # End CG Iterations
 
+            # Take a gradient step on the active cells if exist
+            if temp != self.xc.size:
+
+                rhs_a = (Active) * -self.g
+
+                dm_i = max( abs( delx ) )
+                dm_a = max( abs(rhs_a) )
+
+                delx = delx + rhs_a * dm_i / dm_a /10.
+
+            # Only keep gradients going in the right direction on the active set
+            indx = ((self.xc<=self.lower) & (delx < 0)) | ((self.xc>=self.upper) & (delx > 0))
+            delx[indx] = 0.
+
             return delx

SimPEG/PropMaps.py

@@ -34,18 +34,6 @@ class Property(object):
             setattr(self, '_%sMap'%prop.name, val)
         return property(fget=fget, fset=fset, doc=prop.doc)
 
-    def _getDefaultProperty(self):
-        prop = self
-        def fget(self):
-            return getattr(self, '_%sDefault'%prop.name, None)
-        def fset(self, val):
-            if prop.propertyLink is not None:
-                linkName, linkMap = prop.propertyLink
-                assert getattr(self, '%sDefault'%linkName, None) is None, 'Cannot set both sides of a linked property.'
-            assert isinstance(val, np.ndarray) or np.isscalar(val), 'Default must be a scalar or a numpy array.'
-            setattr(self, '_%sDefault'%prop.name, val)
-        return property(fget=fget, fset=fset, doc=prop.doc)
-
     def _getIndexProperty(self):
         prop = self
         def fget(self):
@@ -59,17 +47,12 @@ class Property(object):
         def fget(self):
             mapping = getattr(self, '%sMap'%prop.name)
             if mapping is None and prop.propertyLink is None:
-                return getattr(self, '%sDefault'%prop.name)
+                return prop.defaultVal
 
             if mapping is None and prop.propertyLink is not None:
                 linkName, linkMapClass = prop.propertyLink
                 linkMap = linkMapClass(None)
-                #                                                                   *
-                print linkName, getattr(self.propMap, '_%sDefault'%linkName, None)
-                if getattr(self, '%sMap'%linkName, None) is None and getattr(self.propMap, '_%sDefault'%linkName, None) is not None:
-                    # We have a default
-                    return linkMap * getattr(self, '%sDefault'%linkName, None)
-                elif getattr(self, '%sMap'%linkName, None) is None:
+                if getattr(self, '%sMap'%linkName, None) is None:
                     return prop.defaultVal
                 m = getattr(self, '%s'%linkName)
                 return linkMap * m
@@ -127,12 +110,6 @@ class Property(object):
             return getattr(self.propMap, '_%sMap'%prop.name, None)
         return property(fget=fget)
 
-    def _getModelDefaultProperty(self):
-        prop = self
-        def fget(self):
-            return getattr(self.propMap, '_%sDefault'%prop.name, prop.defaultVal)
-        return property(fget=fget)
-
 
 
 class PropModel(object):
@@ -173,9 +150,8 @@ class _PropMapMetaClass(type):
         for attr in keys:
             if isinstance(attrs[attr], Property):
                 attrs[attr].name = attr
-                attrs[attr + 'Map'    ] = attrs[attr]._getMapProperty()
-                attrs[attr + 'Default'] = attrs[attr]._getDefaultProperty()
-                attrs[attr + 'Index'  ] = attrs[attr]._getIndexProperty()
+                attrs[attr + 'Map'  ] = attrs[attr]._getMapProperty()
+                attrs[attr + 'Index'] = attrs[attr]._getIndexProperty()
                 _properties[attr] = attrs[attr]
                 attrs.pop(attr)
 
@@ -205,12 +181,11 @@ class _PropMapMetaClass(type):
         for attr in _properties:
             prop = _properties[attr]
 
-            attrs[attr            ] = prop._getProperty()
-            attrs[attr + 'Map'    ] = prop._getModelMapProperty()
-            attrs[attr + 'Default'] = prop._getModelDefaultProperty()
-            attrs[attr + 'Proj'   ] = prop._getModelProjProperty()
-            attrs[attr + 'Model'  ] = prop._getModelProperty()
-            attrs[attr + 'Deriv'  ] = prop._getModelDerivProperty()
+            attrs[attr          ] = prop._getProperty()
+            attrs[attr + 'Map'  ] = prop._getModelMapProperty()
+            attrs[attr + 'Proj' ] = prop._getModelProjProperty()
+            attrs[attr + 'Model'] = prop._getModelProperty()
+            attrs[attr + 'Deriv'] = prop._getModelDerivProperty()
 
         return type(name.replace('PropMap', 'PropModel'), (PropModel, ), attrs)
 
@@ -223,8 +198,8 @@ class PropMap(object):
             PropMap takes a multi parameter model and maps it to the equivalent PropModel
         """
         if type(mappings) is dict:
-            assert np.all([k in ['maps', 'slices', 'defaults'] for k in mappings]), 'Dict must only have properties "maps", "slices" and "defaults"'
-            self.setup(mappings['maps'], slices=mappings.get('slices',{}), defaults=mappings.get('defaults',{}))
+            assert np.all([k in ['maps', 'slices'] for k in mappings]), 'Dict must only have properties "maps" and "slices"'
+            self.setup(mappings['maps'], slices=mappings['slices'])
         elif type(mappings) is list:
             self.setup(mappings)
         elif isinstance(mappings, Maps.IdentityMap):
@@ -233,7 +208,7 @@ class PropMap(object):
             raise Exception('mappings must be a dict, a mapping, or a list of tuples.')
 
 
-    def setup(self, maps, slices=None, defaults=None):
+    def setup(self, maps, slices=None):
         """
             Sets up the maps and slices for the PropertyMap
 
@@ -256,13 +231,6 @@ class PropMap(object):
                 s in self._properties and
                 (type(slices[s]) in [slice, list] or isinstance(slices[s], np.ndarray))
                 for s in slices]), 'Slices must be for each property'
-        if defaults is None:
-            defaults = dict()
-        else:
-            assert np.all([
-                s in self._properties and
-                (np.isscalar(defaults[s]) or isinstance(defaults[s], np.ndarray))
-                for s in defaults]), 'Defaults must be for each property'
 
         self.clearMaps()
 
@@ -271,12 +239,7 @@ class PropMap(object):
             setattr(self, '%sMap'%name, mapping)
             setattr(self, '%sIndex'%name, slices.get(name, slice(nP, nP + mapping.nP)))
             nP += mapping.nP
-        self.nP = nP
-
-        for key in defaults:
-            setattr(self, '%sDefault'%key, defaults[key])
-
-
+        self.nP = nP 
 
     @property
     def defaultInvProp(self):

docs/examples/EM_FDEM_SusEffects.rst

@@ -0,0 +1,41 @@
+.. _examples_EM_FDEM_SusEffects:
+
+.. --------------------------------- ..
+..                                   ..
+..    THIS FILE IS AUTO GENEREATED   ..
+..                                   ..
+..    SimPEG/Examples/__init__.py    ..
+..                                   ..
+.. --------------------------------- ..
+
+
+EM: FDEM: Effects of susceptibility
+===================================
+
+When airborne freqeuncy domain EM (AFEM) survey is flown over
+the earth including significantly susceptible bodies (magnetite-rich rocks),
+negative data is often observed in the real part of the lowest frequency
+(e.g. Dighem system 900 Hz). This phenomenon mostly based upon magnetization
+occurs due to a susceptible body when the magnetic field is applied.
+
+To clarify what is happening in the earth when we are exciting the earth with
+a loop source in the frequency domain we run three forward modelling:
+
+    - F[:math:`\sigma`, :math:`\mu`]: Anomalous conductivity and susceptibility
+    - F[:math:`\sigma`, :math:`\mu_0`]: Anomalous conductivity
+    - F[:math:`\sigma_{air}`, :math:`\mu_0`]: primary field
+
+We plot vector magnetic fields in the earth. For secondary fields we provide
+F[:math:`\sigma`, :math:`\mu`]-F[:math:`\sigma`, :math:`\mu_0`]. Following
+figure show both real and parts.
+
+
+
+.. plot::
+
+    from SimPEG import Examples
+    Examples.EM_FDEM_SusEffects.run()
+
+.. literalinclude:: ../../SimPEG/Examples/EM_FDEM_SusEffects.py
+    :language: python
+    :linenos:

tests/base/test_PropMaps.py

@@ -56,22 +56,8 @@ class TestPropMaps(unittest.TestCase):
             assert np.all(m.sigma == np.exp(np.r_[1.,2,3]))
             assert m.sigmaDeriv is not None
 
-            assert m.mu == mu_0
-
             assert m.nP == 3
 
-    def test_defaultOverride(self):
-        expMap = Maps.ExpMap(Mesh.TensorMesh((3,)))
-        PM = MyReciprocalPropMap({'maps':[('sigma', expMap)], 'defaults':{'mu':mu_0*2}})
-        self.assertRaises(Exception, MyReciprocalPropMap, {'maps':[('sigma', expMap)], 'defaults':{'mu':mu_0*2, 'mui':5}}) # Cannot set both sides of the default
-
-        m = PM(np.r_[1.,2,3])
-        assert np.all(m.sigmaModel == np.r_[1,2,3])
-
-        self.assertEqual(m.mu, mu_0 * 2)
-        # self.assertEqual(m.mui, 1/(mu_0 * 2))
-
-
     def test_slices(self):
         expMap = Maps.ExpMap(Mesh.TensorMesh((3,)))
         PM = MyPropMap({'maps':[('sigma', expMap)], 'slices':{'sigma':[2,1,0]}})