diff --git a/simpegMT/BaseMT.py b/simpegMT/BaseMT.py index 6105d106..7eb2573f 100644 --- a/simpegMT/BaseMT.py +++ b/simpegMT/BaseMT.py @@ -14,33 +14,33 @@ class BaseMTProblem(BaseFDEMProblem): dataPair = DataMT fieldsPair = FieldsMT - # Pickleing support methods - def __getstate__(self): - ''' - Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object. + # # Pickleing support methods + # def __getstate__(self): + # ''' + # Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object. - Used when doing: - pickle.dump(pickleFile,object) - ''' - odict = self.__dict__.copy() - # Remove fields that are not needed - del odict['hook'] - del odict['setKwargs'] - # Return the dict - return odict + # Used when doing: + # pickle.dump(pickleFile,object) + # ''' + # odict = self.__dict__.copy() + # # Remove fields that are not needed + # del odict['hook'] + # del odict['setKwargs'] + # # Return the dict + # return odict - def __setstate__(self,odict): - ''' - Function that sets a pickle dictionary in to an object. + # def __setstate__(self,odict): + # ''' + # Function that sets a pickle dictionary in to an object. - Used when doing: - object = pickle.load(pickleFile) - ''' - # Update the dict - self.__dict__.update(odict) - # Re-hook the methods to the object - Utils.codeutils.hook(self,Utils.codeutils.hook) - Utils.codeutils.hook(self,Utils.codeutils.setKwargs) + # Used when doing: + # object = pickle.load(pickleFile) + # ''' + # # Update the dict + # self.__dict__.update(odict) + # # Re-hook the methods to the object + # Utils.codeutils.hook(self,Utils.codeutils.hook) + # Utils.codeutils.hook(self,Utils.codeutils.setKwargs) # Set the solver Solver = SimpegSolver diff --git a/simpegMT/Sources/backgroundModelSources.py b/simpegMT/Sources/backgroundModelSources.py index e152bbd8..fc18783e 100644 --- a/simpegMT/Sources/backgroundModelSources.py +++ b/simpegMT/Sources/backgroundModelSources.py @@ -20,6 +20,7 @@ def homo1DModelSource(mesh,freq,sigma_1d): mesh1d = simpeg.Mesh.TensorMesh([mesh.hy],np.array([mesh.x0[1]])) elif mesh.dim == 3: mesh1d = simpeg.Mesh.TensorMesh([mesh.hz],np.array([mesh.x0[2]])) + # # Note: Everything is using e^iwt e0_1d = get1DEfields(mesh1d,sigma_1d,freq) if mesh.dim == 1: @@ -62,4 +63,116 @@ def homo1DModelSource(mesh,freq,sigma_1d): # Return the electric fields eBG_bp = np.hstack((eBG_px,eBG_py)) - return eBG_bp \ No newline at end of file + return eBG_bp + +def analytic1DModelSource(mesh,freq,sigma_1d): + ''' + Function that calculates and return background fields + + :param Simpeg mesh object mesh: Holds information on the discretization + :param float freq: The frequency to solve at + :param np.array sigma_1d: Background model of conductivity to base the calculations on, 1d model. + :rtype: numpy.ndarray (mesh.nE,2) + :return: eBG_bp, E fields for the background model at both polarizations. + + ''' + # import + from simpegMT.Utils import getEHfields + # Get a 1d solution for a halfspace background + if mesh.dim == 1: + mesh1d = mesh + elif mesh.dim == 2: + mesh1d = simpeg.Mesh.TensorMesh([mesh.hy],np.array([mesh.x0[1]])) + elif mesh.dim == 3: + mesh1d = simpeg.Mesh.TensorMesh([mesh.hz],np.array([mesh.x0[2]])) + + # # Note: Everything is using e^iwt + Eu, Ed, _, _ = getEHfields(mesh1d,sigma_1d,freq,mesh.vectorNz) + # Make the fields into a dictionary of location and the fields + e0_1d = Eu+Ed + E1dFieldDict = dict(zip(mesh.vectorNz,e0_1d)) + if mesh.dim == 1: + eBG_px = simpeg.mkvc(e0_1d,2) + eBG_py = -simpeg.mkvc(e0_1d,2) # added a minus to make the results in the correct quadrents. + elif mesh.dim == 2: + ex_px = np.zeros(mesh.vnEx,dtype=complex) + ey_px = np.zeros((mesh.nEy,1),dtype=complex) + for i in np.arange(mesh.vnEx[0]): + ex_px[i,:] = -e0_1d + eBG_px = np.vstack((simpeg.Utils.mkvc(ex_px,2),ey_px)) + # Setup y (north) polarization (_py) + ex_py = np.zeros((mesh.nEx,1), dtype='complex128') + ey_py = np.zeros(mesh.vnEy, dtype='complex128') + # Assign the source to ey_py + for i in np.arange(mesh.vnEy[0]): + ey_py[i,:] = e0_1d + # ey_py[1:-1,1:-1,1:-1] = 0 + eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py)) + elif mesh.dim == 3: + # Setup x (east) polarization (_x) + ex_px = -np.array([E1dFieldDict[i] for i in mesh.gridEx[:,2]]).reshape(-1,1) + ey_px = np.zeros((mesh.nEy,1),dtype=complex) + ez_px = np.zeros((mesh.nEz,1),dtype=complex) + # Construct the full fields + eBG_px = np.vstack((ex_px,ey_px,ez_px)) + # Setup y (north) polarization (_py) + ex_py = np.zeros((mesh.nEx,1), dtype='complex128') + ey_py = np.array([E1dFieldDict[i] for i in mesh.gridEy[:,2]]).reshape(-1,1) + ez_py = np.zeros((mesh.nEz,1), dtype='complex128') + # Construct the full fields + eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py)) + + # Return the electric fields + eBG_bp = np.hstack((eBG_px,eBG_py)) + return eBG_bp + +# def homo3DModelSource(mesh,model,freq): +# ''' +# Function that estimates 1D analytic background fields from a 3D model. + +# :param Simpeg mesh object mesh: Holds information on the discretization +# :param float freq: The frequency to solve at +# :param np.array sigma_1d: Background model of conductivity to base the calculations on, 1d model. +# :rtype: numpy.ndarray (mesh.nE,2) +# :return: eBG_bp, E fields for the background model at both polarizations. + +# ''' + +# if mesh.dim < 3: +# raise IOError('Input mesh has to have 3 dimensions.') + + +# # Get the locations +# a = mesh.gridCC[:,0:2].copy() +# unixy = np.unique(a.view(a.dtype.descr * a.shape[1])).view(float).reshape(-1,2) +# uniz = np.unique(mesh.gridCC[:,2]) +# # # Note: Everything is using e^iwt +# # Need to loop thourgh the xy locations, assess the model and calculate the fields at the phusdo cell centers. +# # Then interpolate the cc fields to the edges. + +# e0_1d = get1DEfields(mesh1d,sigma_1d,freq) + +# elif mesh.dim == 3: +# # Setup x (east) polarization (_x) +# ex_px = np.zeros(mesh.vnEx,dtype=complex) +# ey_px = np.zeros((mesh.nEy,1),dtype=complex) +# ez_px = np.zeros((mesh.nEz,1),dtype=complex) +# # Assign the source to ex_x +# for i in np.arange(mesh.vnEx[0]): +# for j in np.arange(mesh.vnEx[1]): +# ex_px[i,j,:] = -e0_1d +# eBG_px = np.vstack((simpeg.Utils.mkvc(ex_px,2),ey_px,ez_px)) +# # Setup y (north) polarization (_py) +# ex_py = np.zeros((mesh.nEx,1), dtype='complex128') +# ey_py = np.zeros(mesh.vnEy, dtype='complex128') +# ez_py = np.zeros((mesh.nEz,1), dtype='complex128') +# # Assign the source to ey_py +# for i in np.arange(mesh.vnEy[0]): +# for j in np.arange(mesh.vnEy[1]): +# ey_py[i,j,:] = e0_1d +# # ey_py[1:-1,1:-1,1:-1] = 0 +# eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py)) + +# # Return the electric fields +# eBG_bp = np.hstack((eBG_px,eBG_py)) +# return eBG_bp \ No newline at end of file diff --git a/simpegMT/SurveyMT.py b/simpegMT/SurveyMT.py index 7547ca5d..c8a759cc 100644 --- a/simpegMT/SurveyMT.py +++ b/simpegMT/SurveyMT.py @@ -430,6 +430,82 @@ class srcMT_polxy_1Dprimary(srcMT): # v should be nC size return MsigmaDeriv * v +class srcMT_polxy_3Dprimary(srcMT): + """ + MT source for both polarizations (x and y) given a 3D primary model. It assigns fields calculated from the 1D model + as fields in the full space of the problem. + """ + def __init__(self, rxList, freq): + # assert mkvc(self.mesh.hz.shape,1) == mkvc(sigma1d.shape,1),'The number of values in the 1D background model does not match the number of vertical cells (hz).' + self.sigmaPrimary = None + srcMT.__init__(self, rxList, freq) + # Hidden property of the ePrimary + self._ePrimary = None + + def ePrimary(self,problem): + # Get primary fields for both polarizations + self.sigmaPrimary = problem._sigmaPrimary + + if self._ePrimary is None: + self._ePrimary = homo3DModelSource(problem.mesh,self.sigmaPrimary,self.freq) + return self._ePrimary + + def bPrimary(self,problem): + # Project ePrimary to bPrimary + # Satisfies the primary(background) field conditions + if problem.mesh.dim == 1: + C = problem.mesh.nodalGrad + elif problem.mesh.dim == 3: + C = problem.mesh.edgeCurl + bBG_bp = (- C * self.ePrimary(problem) )*(1/( 1j*omega(self.freq) )) + return bBG_bp + + def S_e(self,problem): + """ + Get the electrical field source + """ + e_p = self.ePrimary(problem) + Map_sigma_p = Maps.Vertical1DMap(problem.mesh) + sigma_p = Map_sigma_p._transform(self.sigma1d) + # Make mass matrix + # Note: M(sig) - M(sig_p) = M(sig - sig_p) + # Need to deal with the edge/face discrepencies between 1d/2d/3d + if problem.mesh.dim == 1: + Mesigma = problem.mesh.getFaceInnerProduct(problem.curModel.sigma) + Mesigma_p = problem.mesh.getFaceInnerProduct(sigma_p) + if problem.mesh.dim == 2: + pass + if problem.mesh.dim == 3: + Mesigma = problem.MeSigma + Mesigma_p = problem.mesh.getEdgeInnerProduct(sigma_p) + return (Mesigma - Mesigma_p) * e_p + + def S_eDeriv_m(self, problem, v, adjoint = False): + ''' + Get the derivative of S_e wrt to sigma (m) + ''' + # Need to deal with + if problem.mesh.dim == 1: + # Need to use the faceInnerProduct + MsigmaDeriv = problem.mesh.getFaceInnerProductDeriv(problem.curModel.sigma)(self.ePrimary(problem)[:,1]) * problem.curModel.sigmaDeriv + # MsigmaDeriv = ( MsigmaDeriv * MsigmaDeriv.T)**2 + if problem.mesh.dim == 2: + pass + if problem.mesh.dim == 3: + # Need to take the derivative of both u_px and u_py + ePri = self.ePrimary(problem) + # MsigmaDeriv = problem.MeSigmaDeriv(ePri[:,0]) + problem.MeSigmaDeriv(ePri[:,1]) + # MsigmaDeriv = problem.MeSigmaDeriv(np.sum(ePri,axis=1)) + if adjoint: + return sp.hstack(( problem.MeSigmaDeriv(ePri[:,0]).T, problem.MeSigmaDeriv(ePri[:,1]).T ))*v + else: + return np.hstack(( mkvc(problem.MeSigmaDeriv(ePri[:,0]) * v,2), mkvc(problem.MeSigmaDeriv(ePri[:,1])*v,2) )) + if adjoint: + # + return MsigmaDeriv.T * v + else: + # v should be nC size + return MsigmaDeriv * v ############## ### Survey ###