diff --git a/SimPEG/EM/TDEM/BaseTDEM.py b/SimPEG/EM/TDEM/BaseTDEM.py index e36d76b8..4d4b39a7 100644 --- a/SimPEG/EM/TDEM/BaseTDEM.py +++ b/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) diff --git a/SimPEG/EM/TDEM/SurveyTDEM.py b/SimPEG/EM/TDEM/SurveyTDEM.py index 7f6e5c04..6a65c0fe 100644 --- a/SimPEG/EM/TDEM/SurveyTDEM.py +++ b/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): diff --git a/SimPEG/Examples/EM_FDEM_Analytic_MagDipoleWholespace.py b/SimPEG/Examples/EM_FDEM_Analytic_MagDipoleWholespace.py new file mode 100644 index 00000000..791f8c63 --- /dev/null +++ b/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() diff --git a/SimPEG/Examples/EM_FDEM_1D_Inversion.py b/SimPEG/Examples/EM_TDEM_1D_Inversion.py similarity index 97% rename from SimPEG/Examples/EM_FDEM_1D_Inversion.py rename to SimPEG/Examples/EM_TDEM_1D_Inversion.py index aba70f4b..d4d80e55 100644 --- a/SimPEG/Examples/EM_FDEM_1D_Inversion.py +++ b/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. """ diff --git a/SimPEG/Examples/__init__.py b/SimPEG/Examples/__init__.py index 8431e4ba..7ef15451 100644 --- a/SimPEG/Examples/__init__.py +++ b/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() diff --git a/SimPEG/Maps.py b/SimPEG/Maps.py index 0c86d476..d14400f8 100644 --- a/SimPEG/Maps.py +++ b/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. diff --git a/SimPEG/Problem.py b/SimPEG/Problem.py index da03d402..9289a95e 100644 --- a/SimPEG/Problem.py +++ b/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' diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 97b5e7a2..6c5b5f97 100644 --- a/SimPEG/Regularization.py +++ b/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 diff --git a/docs/examples/EM_FDEM_Analytic_MagDipoleWholespace.rst b/docs/examples/EM_FDEM_Analytic_MagDipoleWholespace.rst new file mode 100644 index 00000000..0e1b2d93 --- /dev/null +++ b/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: diff --git a/docs/examples/EM_FDEM_1D_Inversion.rst b/docs/examples/EM_TDEM_1D_Inversion.rst similarity index 64% rename from docs/examples/EM_FDEM_1D_Inversion.rst rename to docs/examples/EM_TDEM_1D_Inversion.rst index acbc8cdc..53f6f9ef 100644 --- a/docs/examples/EM_FDEM_1D_Inversion.rst +++ b/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: