Merge branch 'feat/GlobalProblem' of https://github.com/simpeg/simpeg into feat/GlobalProblem

SimPEG/EM/TDEM/BaseTDEM.py

@@ -31,7 +31,7 @@ class FieldsTDEM(Problem.TimeFields):
 
 
 class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
-    """docstring for ProblemTDEM1D"""
+    """docstring for BaseTDEMProblem"""
     def __init__(self, mesh, mapping=None, **kwargs):
         BaseTimeProblem.__init__(self, mesh, mapping=mapping, **kwargs)
 
@@ -43,7 +43,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
         # Create a fields storage object
         F = self._FieldsForward_pair(self.mesh, self.survey)
         for src in self.survey.srcList:
-            # Set the initial conditions
+            # Set the initial conditions            
             F[src,:,0] = src.getInitialFields(self.mesh)
         F = self.forward(m, self.getRHS, F=F)
         if self.verbose: print '%s\nDone calculating fields(m)\n%s'%('*'*50,'*'*50)

SimPEG/EM/TDEM/SurveyTDEM.py

@@ -99,14 +99,33 @@ class SrcTDEM_VMD_MVP(SrcTDEM):
 
 
 class SrcTDEM_CircularLoop_MVP(SrcTDEM):
-
-    def __init__(self,rxList,loc,radius):
+    def __init__(self,rxList,loc,radius,waveformType):
         self.loc = loc
         self.radius = radius
-        SrcTDEM.__init__(self,rxList)
+        self.waveformType = waveformType
+        SrcTDEM.__init__(self,rxList)        
 
     def getInitialFields(self, mesh):
         """Circular Loop, magnetic vector potential"""
+        if self.waveformType == "STEPOFF":
+            print ">> Step waveform: Non-zero initial condition"
+            if mesh._meshType is 'CYL':
+                if mesh.isSymmetric:
+                    MVP = MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
+                else:
+                    raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
+            elif mesh._meshType is 'TENSOR':
+                MVP = MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius)
+            else:
+                raise Exception('Unknown mesh for CircularLoop')
+            return {"b": mesh.edgeCurl*MVP}
+        elif self.waveformType == "GENERAL":
+            print ">> General waveform: Zero initial condition"
+            return {"b": np.zeros(mesh.nF)}
+        else:
+            raise NotImplementedError("Only use STEPOFF or GENERAL")
+
+    def getMeS(self, mesh, MfMui):
         if mesh._meshType is 'CYL':
             if mesh.isSymmetric:
                 MVP = MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
@@ -115,9 +134,8 @@ class SrcTDEM_CircularLoop_MVP(SrcTDEM):
         elif mesh._meshType is 'TENSOR':
             MVP = MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius)
         else:
-            raise Exception('Unknown mesh for CircularLoop')
-
-        return {"b": mesh.edgeCurl*MVP}
+            raise Exception('Unknown mesh for CircularLoop')            
+        return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP
 
 
 class SurveyTDEM(Survey.BaseSurvey):

SimPEG/Examples/EM_FDEM_Analytic_MagDipoleWholespace.py

@@ -0,0 +1,43 @@
+from SimPEG import *
+import SimPEG.EM as EM
+
+def run(XYZ=None, sig=1.0, freq=1.0, orientation='Z', plotIt=True):
+    """
+        EM: Magnetic Dipole in a Whole-Space
+        ====================================
+
+        Here we plot the magnetic flux density from a harmonic dipole in a wholespace.
+
+    """
+
+    if XYZ is None:
+        x = np.arange(-100.5,100.5,step = 1.) #(avoid putting measurement points where source is located)
+        y = np.r_[0]
+        z = x
+        XYZ = Utils.ndgrid(x,y,z)
+
+
+    Bx, By, Bz = EM.Analytics.FDEM.MagneticDipoleWholeSpace(XYZ, np.r_[0.,0.,0.], sig, freq, orientation=orientation)
+    absB = np.sqrt(Bx*Bx.conj()+By*By.conj()+Bz*Bz.conj()).real
+
+
+    if plotIt:
+        import matplotlib.pyplot as plt
+        from matplotlib.colors import LogNorm
+        fig, ax = plt.subplots(1,1,figsize=(6,5))
+        bxplt = Bx.reshape(x.size,z.size)
+        bzplt = Bz.reshape(x.size,z.size)
+        pc = ax.pcolor(x,z,absB.reshape(x.size,z.size),norm=LogNorm())
+        ax.streamplot(x,z,bxplt.real,bzplt.real,color='k',density=1)
+        ax.set_xlim([x.min(),x.max()])
+        ax.set_ylim([z.min(),z.max()])
+        ax.set_xlabel('x')
+        ax.set_ylabel('z')
+        cb = plt.colorbar(pc,ax = ax)
+        cb.set_label('|B| (T)')
+        plt.show()
+
+        return fig, ax
+
+if __name__ == '__main__':
+    run()

SimPEG/Examples/EM_TDEM_1D_Inversion.py

@@ -5,10 +5,10 @@ from scipy.constants import mu_0
 
 def run(plotIt=True):
     """
-        EM: FDEM: 1D: Inversion
+        EM: TDEM: 1D: Inversion
         =======================
 
-        Here we will create and run a FDEM 1D inversion.
+        Here we will create and run a TDEM 1D inversion.
 
     """
 

SimPEG/Examples/__init__.py

@@ -1,7 +1,8 @@
 # Run this file to add imports.
 
 ##### AUTOIMPORTS #####
-import EM_FDEM_1D_Inversion
+import EM_FDEM_Analytic_MagDipoleWholespace
+import EM_TDEM_1D_Inversion
 import FLOW_Richards_1D_Celia1990
 import Forward_BasicDirectCurrent
 import Inversion_Linear
@@ -13,7 +14,7 @@ import Mesh_QuadTree_FaceDiv
 import Mesh_QuadTree_HangingNodes
 import Mesh_Tensor_Creation
 
-__examples__ = ["EM_FDEM_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_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"]
 
 ##### AUTOIMPORTS #####
 
@@ -28,16 +29,17 @@ if __name__ == '__main__':
     from SimPEG import Examples
 
     # Create the examples dir in the docs folder.
-    docExamplesDir = os.path.sep.join(os.path.realpath(__file__).split(os.path.sep)[:-3] + ['docs', 'examples'])
+    fName = os.path.realpath(__file__)
+    docExamplesDir = os.path.sep.join(fName.split(os.path.sep)[:-3] + ['docs', 'examples'])
     shutil.rmtree(docExamplesDir)
     os.makedirs(docExamplesDir)
 
     # Get all the python examples in this folder
-    thispath = os.path.sep.join(__file__.split(os.path.sep)[:-1])
+    thispath = os.path.sep.join(fName.split(os.path.sep)[:-1])
     exfiles  = [f[:-3] for f in os.listdir(thispath) if os.path.isfile(os.path.join(thispath, f)) and f.endswith('.py') and not f.startswith('_')]
 
     # Add the imports to the top in the AUTOIMPORTS section
-    f = file(__file__, 'r')
+    f = file(fName, 'r')
     inimports = False
     out = ''
     for line in f:
@@ -52,7 +54,7 @@ if __name__ == '__main__':
                 out += '\n##### AUTOIMPORTS #####\n'
     f.close()
 
-    f = file(__file__, 'w')
+    f = file(fName, 'w')
     f.write(out)
     f.close()
 

SimPEG/Maps.py

@@ -449,6 +449,84 @@ class Mesh2Mesh(IdentityMap):
         return self.P
 
 
+class Mesh2MeshTopo(IdentityMap):
+    """
+        Takes a model on one mesh are translates it to another mesh
+        with consideration of topography
+
+    """
+    tree = None
+    nIterpPts = 6
+
+    def __init__(self, meshes, actinds, **kwargs):
+        Utils.setKwargs(self, **kwargs)
+
+        assert type(meshes) is list, "meshes must be a list of two meshes"
+        assert len(meshes) == 2, "meshes must be a list of two meshes"
+        assert type(actinds) is list, "actinds must be a list of two meshes"
+        assert len(actinds) == 2, "actinds must be a list of two meshes"
+        assert meshes[0].dim == meshes[1].dim, """The two meshes must be the same dimension"""
+
+        self.mesh  = meshes[0]
+        self.mesh2 = meshes[1]
+        self.actind = actinds[0]
+        self.actind2 = actinds[1]
+        self.getP()
+
+        # Old version using SimPEG interpolation
+        # self.P = self.mesh2.getInterpolationMat(self.mesh.gridCC,'CC',zerosOutside=True)
+    
+    from scipy.interpolate import NearestNDInterpolator        
+    def genActiveindfromTopo(mesh, xyztopo):    
+        #TODO: This possibly needs to be improved use vtk(?)
+        if mesh.dim==3:
+            nCxy = mesh.nCx*mesh.nCy
+            Zcc = mesh.gridCC[:,2].reshape((nCxy, mesh.nCz), order='F')
+            Ftopo = NearestNDInterpolator(xyztopo[:,:2], xyztopo[:,2])
+            XY = Utils.ndgrid(mesh.vectorCCx, mesh.vectorCCy)
+            XY.shape
+            topo = Ftopo(XY)
+            actind = []                    
+            for ixy in range(nCxy):
+                actind.append(topo[ixy] <= Zcc[ixy,:])
+        else:
+            raise NotImplementedError("Only 3D is working")
+
+        return Utils.mkvc(np.vstack(actind))
+
+    #Question .. is it only generated once?
+    def getP(self):
+        """
+
+        """
+        if self.tree==None:
+            self.tree = cKDTree(zip(self.mesh.gridCC[self.actind,0], self.mesh.gridCC[self.actind,1], self.mesh.gridCC[self.actind,2]))
+        d, inds = tree.query(zip(self.mesh2.gridCC[self.actind2,0],self.mesh2.gridCC[self.actind2,1],self.mesh2.gridCC[self.actind2,2]), k=self.nIterpPts)
+        w = 1./ d**2
+        w = Utils.sdiag(1./np.sum(w, axis=1)) * w
+        I = Utils.mkvc(np.arange(inds.shape[0]).reshape([-1,1]).repeat(6, axis=1))
+        J = Utils.mkvc(inds)
+        P = sp.coo_matrix( (Utils.mkvc(w),(I, J)), shape=(inds.shape[0], (self.actind).sum()) )
+        self.P = P.tocsc() 
+
+    @property
+    def shape(self):
+        """Number of parameters in the model."""
+        # return (self.mesh.nC, self.mesh2.nC)
+        return (self.actind2.sum(), self.actind.sum())
+
+    @property
+    def nP(self):
+        """Number of parameters in the model."""
+        # return self.mesh2.nC
+        return self.actind2.sum()
+
+    def _transform(self, m):
+        return self.P*m
+
+    def deriv(self, m):
+        return self.P
+
 class ActiveCells(IdentityMap):
     """
         Active model parameters.

SimPEG/Problem.py

@@ -166,6 +166,9 @@ class BaseProblem(object):
 
 class BaseTimeProblem(BaseProblem):
     """Sets up that basic needs of a time domain problem."""
+    
+    waveformType = "STEPOFF"    
+    current = None
 
     @property
     def timeSteps(self):
@@ -192,6 +195,11 @@ class BaseTimeProblem(BaseProblem):
         self._timeSteps = Utils.meshTensor(value)
         del self.timeMesh
 
+    def currentwaveform(self, wave):
+        self._timeSteps = np.diff(wave[:,0])
+        self.current = wave[:,1]
+        self.waveformType = "GENERAL"
+
     @property
     def nT(self):
         "Number of time steps."
@@ -306,7 +314,7 @@ class GlobalProblem(BaseProblem):
             raise Exceptions.PairingException(reason='The meshes are not the the same length as the number of groups')
 
     def getSubProblem(self, ind):
-
+        #This is a core place that we can proceed parallelization
         assert self.ispaired, 'You must be paired to a survey'
         assert type(ind) in [int,long] and ind >= 0 and ind < self.nGroups, 'ind must be an index into the group list'
 

SimPEG/Regularization.py

@@ -24,7 +24,7 @@ class BaseRegularization(object):
         Utils.setKwargs(self, **kwargs)
         self.mesh = mesh
         assert isinstance(mesh, Mesh.BaseMesh), "mesh must be a SimPEG.Mesh object."
-        self.mapping = mapping or Maps.IdentityMap(mesh)
+        self.mapping = mapping or self.mapPair(mesh)
         self.mapping._assertMatchesPair(self.mapPair)
 
     @property

docs/examples/EM_FDEM_Analytic_MagDipoleWholespace.rst

@@ -0,0 +1,26 @@
+.. _examples_EM_FDEM_Analytic_MagDipoleWholespace:
+
+.. --------------------------------- ..
+..                                   ..
+..    THIS FILE IS AUTO GENEREATED   ..
+..                                   ..
+..    SimPEG/Examples/__init__.py    ..
+..                                   ..
+.. --------------------------------- ..
+
+
+EM: Magnetic Dipole in a Whole-Space
+====================================
+
+Here we plot the magnetic flux density from a harmonic dipole in a wholespace.
+
+
+
+.. plot::
+
+    from SimPEG import Examples
+    Examples.EM_FDEM_Analytic_MagDipoleWholespace.run()
+
+.. literalinclude:: ../../SimPEG/Examples/EM_FDEM_Analytic_MagDipoleWholespace.py
+    :language: python
+    :linenos:

docs/examples/EM_TDEM_1D_Inversion.rst

@@ -1,4 +1,4 @@
-.. _examples_EM_FDEM_1D_Inversion:
+.. _examples_EM_TDEM_1D_Inversion:
 
 .. --------------------------------- ..
 ..                                   ..
@@ -9,18 +9,18 @@
 .. --------------------------------- ..
 
 
-EM: FDEM: 1D: Inversion
+EM: TDEM: 1D: Inversion
 =======================
 
-Here we will create and run a FDEM 1D inversion.
+Here we will create and run a TDEM 1D inversion.
 
 
 
 .. plot::
 
     from SimPEG import Examples
-    Examples.EM_FDEM_1D_Inversion.run()
+    Examples.EM_TDEM_1D_Inversion.run()
 
-.. literalinclude:: ../../SimPEG/Examples/EM_FDEM_1D_Inversion.py
+.. literalinclude:: ../../SimPEG/Examples/EM_TDEM_1D_Inversion.py
     :language: python
     :linenos: