diff --git a/simpegMT/FieldsMT.py b/simpegMT/FieldsMT.py index 622460a3..3dea0c11 100644 --- a/simpegMT/FieldsMT.py +++ b/simpegMT/FieldsMT.py @@ -45,7 +45,7 @@ class FieldsMT_1D(FieldsMT): return self._ePrimary(eSolution,srcList) + self._eSecondary(eSolution,srcList) def _eDeriv_u(self, src, v, adjoint = False): - return None + return v def _eDeriv_m(self, src, v, adjoint = False): # assuming primary does not depend on the model @@ -124,23 +124,199 @@ class FieldsMT_3D(FieldsMT): """ Fields storage for the 3D MT solution. """ - knownFields = {'e_px':'E','e_py':'E','b_px':'F','b_py':'F'} - aliasFields = { } - # 'e_1d' : ['e_1dSolution','F','_e'], - # 'e_1dPrimary' : ['e_1dSolution','F','_ePrimary'], - # 'e_1dSecondary' : ['e_1dSolution','F','_eSecondary'], - # 'b_1d' : ['e_1dSolution','E','_b'], - # 'b_1dPrimary' : ['e_1dSolution','E','_bPrimary'], - # 'b_1dSecondary' : ['e_1dSolution','E','_bSecondary'] - # } + # Define the known the alias fields + # Assume that the solution of e on the E. + ## NOTE: Need to make this more general, to allow for other solutions formats. + knownFields = {'e_pxSolution':'E','e_pySoluiton':'E'} + aliasFields = { + 'e_px' : ['e_pxSolution','E','_e_px'], + 'e_pxPrimary' : ['e_pxSolution','E','_e_pxPrimary'], + 'e_pxSecondary' : ['e_pxSolution','E','_e_pxSecondary'], + 'b_px' : ['e_pxSolution','F','_b_px'], + 'b_pxPrimary' : ['e_pxSolution','F','_b_pxPrimary'], + 'b_pxSecondary' : ['e_pxSolution','F','_b_pxSecondary'] + } + def __init__(self,mesh,survey,**kwargs): + FieldsMT.__init__(self,mesh,survey,**kwargs) - # knownFields = {'e_pxSolution':'E','e_pySoluiton':'E'} - # aliasFields = { - # 'e_px' : ['e_pxSolution','E','_epx'], - # 'e_pxPrimary' : ['e_pxSolution','E','_epxPrimary'], - # 'e_pxSecondary' : ['e_pxSolution','E','_epxSecondary'], - # 'b_px' : ['e_pxSolution','F','_bpx'], - # 'b_pxPrimary' : ['e_pxSolution','F','_bpxPrimary'], - # 'b_pxSecondary' : ['e_pxSolution','F','_bpxSecondary'] - # } \ No newline at end of file + def _e_pxPrimary(self, e_pxSolution, srcList): + e_pxPrimary = np.zeros_like(e_pxSolution) + for i, src in enumerate(srcList): + ep = src.ePrimary(self.survey.prob) + if ep is not None: + e_pxPrimary[:,i] = ep[:,0] + return e_pxPrimary + + def _e_pyPrimary(self, e_pySolution, srcList): + e_pyPrimary = np.zeros_like(e_pySolution) + for i, src in enumerate(srcList): + ep = src.ePrimary(self.survey.prob) + if ep is not None: + e_pyPrimary[:,i] = ep[:,1] + return e_pyPrimary + + def _e_pxSecondary(self, e_pxSolution, srcList): + return e_pxSolution + + def _e_pySecondary(self, e_pySolution, srcList): + return e_pySolution + + def _e_px(self, e_pxSolution, srcList): + return self._e_pxPrimary(e_pxSolution,srcList) + self._e_pxSecondary(e_pxSolution,srcList) + + def _e_py(self, e_pySolution, srcList): + return self._e_pyPrimary(e_pySolution,srcList) + self._e_pySecondary(e_pySolution,srcList) + + def _e_pxDeriv_u(self, src, v, adjoint = False): + return v + + def _e_pyDeriv_u(self, src, v, adjoint = False): + return v + + def _e_pxDeriv_m(self, src, v, adjoint = False): + # assuming primary does not depend on the model + return None + def _e_pyDeriv_m(self, src, v, adjoint = False): + # assuming primary does not depend on the model + return None + + def _b_pxPrimary(self, e_pxSolution, srcList): + b_pxPrimary = np.zeros([self.survey.mesh.nE,e_pxSolution.shape[1]], dtype = complex) + for i, src in enumerate(srcList): + bp = src.bPrimary(self.survey.prob) + if bp is not None: + b_pxPrimary[:,i] += bp[:,-1] + return b_pxPrimary + + def _b_pyPrimary(self, e_pySolution, srcList): + b_pyPrimary = np.zeros([self.survey.mesh.nE,e_pySolution.shape[1]], dtype = complex) + for i, src in enumerate(srcList): + bp = src.bPrimary(self.survey.prob) + if bp is not None: + b_pyPrimary[:,i] += bp[:,-1] + return b_pyPrimary + + def _b_pxSecondary(self, e_pxSolution, srcList): + C = self.mesh.edgeCurl + b = (C * e_pxSolution) + for i, src in enumerate(srcList): + b[:,i] *= - 1./(1j*omega(src.freq)) + # There is no magnetic source in the MT problem + # S_m, _ = src.eval(self.survey.prob) + # if S_m is not None: + # b[:,i] += 1./(1j*omega(src.freq)) * S_m + return b + + def _b_pySecondary(self, e_pySolution, srcList): + C = self.mesh.edgeCurl + b = (C * e_pySolution) + for i, src in enumerate(srcList): + b[:,i] *= - 1./(1j*omega(src.freq)) + # There is no magnetic source in the MT problem + # S_m, _ = src.eval(self.survey.prob) + # if S_m is not None: + # b[:,i] += 1./(1j*omega(src.freq)) * S_m + return b + + def _b_px(self, eSolution, srcList): + return self._bPrimary(eSolution, srcList) + self._bSecondary(eSolution, srcList) + + def _b_pxSecondaryDeriv_u(self, src, v, adjoint = False): + C = self.mesh.edgeCurl + if adjoint: + return - 1./(1j*omega(src.freq)) * (C.T * v) + return - 1./(1j*omega(src.freq)) * (C * v) + + def _b_py(self, eSolution, srcList): + return self._b_pyPrimary(eSolution, srcList) + self._b_pySecondary(eSolution, srcList) + + def _b_pxSecondaryDeriv_u(self, src, v, adjoint = False): + C = self.mesh.edgeCurl + if adjoint: + return - 1./(1j*omega(src.freq)) * (C.T * v) + return - 1./(1j*omega(src.freq)) * (C * v) + + def _b_pySecondaryDeriv_m(self, src, v, adjoint = False): + # Doesn't depend on m + # _, S_eDeriv = src.evalDeriv(self.survey.prob, adjoint) + # S_eDeriv = S_eDeriv(v) + # if S_eDeriv is not None: + # return 1./(1j * omega(src.freq)) * S_eDeriv + return None + + def _b_pxSecondaryDeriv_m(self, src, v, adjoint = False): + # Doesn't depend on m + # _, S_eDeriv = src.evalDeriv(self.survey.prob, adjoint) + # S_eDeriv = S_eDeriv(v) + # if S_eDeriv is not None: + # return 1./(1j * omega(src.freq)) * S_eDeriv + return None + + def _b_pxDeriv_u(self, src, v, adjoint=False): + # Primary does not depend on u + return self._b_pxSecondaryDeriv_u(src, v, adjoint) + + def _b_pyDeriv_u(self, src, v, adjoint=False): + # Primary does not depend on u + return self._b_pySecondaryDeriv_u(src, v, adjoint) + + def _b_pxDeriv_m(self, src, v, adjoint=False): + # Assuming the primary does not depend on the model + return self._b_pxSecondaryDeriv_m(src, v, adjoint) + + def _b_pyDeriv_m(self, src, v, adjoint=False): + # Assuming the primary does not depend on the model + return self._b_pySecondaryDeriv_m(src, v, adjoint) + + def _f_pxDeriv_u(self, src, v, adjoint=False): + """ + Derivative of the fields object wrt u. + + :param MTsrc src: MT source + :param numpy.ndarray v: random vector of f_sol.size + This function stacks the fields derivatives appropriately + + return a vector of size (nreEle+nrbEle) + """ + + de_du = v #Utils.spdiag(np.ones((self.nF,))) + db_du = self._b_pxDeriv_u(src, v, adjoint) + # Return the stack + # This doesn't work... + return np.vstack((de_du,db_du)) + + def _f_pyDeriv_u(self, src, v, adjoint=False): + """ + Derivative of the fields object wrt u. + + :param MTsrc src: MT source + :param numpy.ndarray v: random vector of f_sol.size + This function stacks the fields derivatives appropriately + + return a vector of size (nreEle+nrbEle) + """ + + de_du = v #Utils.spdiag(np.ones((self.nF,))) + db_du = self._b_pyDeriv_u(src, v, adjoint) + # Return the stack + # This doesn't work... + return np.vstack((de_du,db_du)) + + def _f_pxDeriv_m(self, src, v, adjoint=False): + """ + Derivative of the fields object wrt m. + + This function stacks the fields derivatives appropriately + """ + # The fields have no dependance to the model. + return None + + def _f_pyDeriv_m(self, src, v, adjoint=False): + """ + Derivative of the fields object wrt m. + + This function stacks the fields derivatives appropriately + """ + # The fields have no dependance to the model. + return None \ No newline at end of file diff --git a/simpegMT/ProblemMT3D/Problems.py b/simpegMT/ProblemMT3D/Problems.py index 6574a78f..4804be33 100644 --- a/simpegMT/ProblemMT3D/Problems.py +++ b/simpegMT/ProblemMT3D/Problems.py @@ -104,16 +104,13 @@ class eForm_ps(BaseMTProblem): # Store the fields Src = self.survey.getSrcByFreq(freq)[0] - # Calculate total e - - e = Src.ePrimary(self) + e_s # Store the fieldss - F[Src, 'e_px'] = e[:,0] - F[Src, 'e_py'] = e[:,1] + F[Src, 'e_pxSolution'] = e_s[:,0] + F[Src, 'e_pySolution'] = e_s[:,1] # Note curl e = -iwb so b = -curl/iw - b = -( self.mesh.edgeCurl * e )/( 1j*omega(freq) ) - F[Src, 'b_px'] = b[:,0] - F[Src, 'b_py'] = b[:,1] + # b = -( self.mesh.edgeCurl * e )/( 1j*omega(freq) ) + # F[Src, 'b_px'] = b[:,0] + # F[Src, 'b_py'] = b[:,1] if self.verbose: print 'Ran for {:f} seconds'.format(time.time()-startTime) sys.stdout.flush() diff --git a/simpegMT/SurveyMT.py b/simpegMT/SurveyMT.py index 766e0472..8e877c4f 100644 --- a/simpegMT/SurveyMT.py +++ b/simpegMT/SurveyMT.py @@ -175,7 +175,57 @@ class RxMT(Survey.BaseRx): elif self.projType is 'Z2D': raise NotImplementedError('Has not be implement for 2D impedance tensor') elif self.projType is 'Z3D': - raise NotImplementedError('Has not be implement for full 3D impedance tensor') + if self.locs.ndim == 3: + eFLocs = self.locs[:,:,0] + bFLocs = self.locs[:,:,1] + else: + eFLocs = self.locs + bFLocs = self.locs + # Get the projection + Pex = mesh.getInterpolationMat(eFLocs,'Ex') + Pey = mesh.getInterpolationMat(eFLocs,'Ey') + Pbx = mesh.getInterpolationMat(bFLocs,'Fx') + Pby = mesh.getInterpolationMat(bFLocs,'Fy') + # Get the fields at location + # px: x-polaration and py: y-polaration. + ex_px = Pex*f[src,'e_px'] + ey_px = Pey*f[src,'e_px'] + ex_py = Pex*f[src,'e_py'] + ey_py = Pey*f[src,'e_py'] + hx_px = Pbx*f[src,'b_px']/mu_0 + hy_px = Pby*f[src,'b_px']/mu_0 + hx_py = Pbx*f[src,'b_py']/mu_0 + hy_py = Pby*f[src,'b_py']/mu_0 + # Derivatives as lambda functions + ex_px_u = lambda vec: Pex*f._e_pxDeriv_u(src,vec) + ey_px_u = lambda vec: Pey*f._e_pxDeriv_u(src,vec) + ex_py_u = lambda vec: Pex*f._e_pyDeriv_u(src,vec) + ey_py_u = lambda vec: Pey*f._e_pyDeriv_u(src,vec) + hx_px_u = lambda vec: Pbx*f._b_pxDeriv_u(src,vec)/mu_0 + hy_px_u = lambda vec: Pby*f._b_pxDeriv_u(src,vec)/mu_0 + hx_py_u = lambda vec: Pbx*f._b_pyDeriv_u(src,vec)/mu_0 + hy_py_u = lambda vec: Pby*f._b_pyDeriv_u(src,vec)/mu_0 + + # Define the components of the derivative + Hd = (hx_px*hy_py - hx_py*hy_px) + Hd_uV = hx_px_u(hy_py)*v - hx_py*hy_px_u(v) + hx_px*hy_py_u(v) - hx_py_u(hy_px) + + if 'zxx' in self.rxType: + Zij = ( ex_px*hy_py - ex_py*hy_px)/Hd + ZijN_uV = ex_px_u(hy_py)*v - ex_py*hy_px_u(v) + ex_px*hy_py_u(v) - ex_py_u(hy_px)*v + elif 'zxy' in self.rxType: + Zij = (-ex_px*hx_py + ex_py*hx_px)/Hd + ZijN_uV = -ex_px_u(hx_py)*v + ex_py*hx_px_u(v) - ex_px*hx_py_u(v) + ex_py_u(hx_px)*v + elif 'zyx' in self.rxType: + Zij = ( ey_px*hy_py - ey_py*hy_px)/Hd + ZijN_uV = ey_px_u(hy_py)*v - ey_py*hy_px_u(v) + ey_px*hy_py_u(v) - ey_py_u(hy_px)*v + elif 'zyy' in self.rxType: + Zij = (-ey_px*hx_py + ey_py*hx_px)/Hd + ZijN_uV = -ey_px_u(hx_py)*v + ey_py*hx_px_u(v) - ey_px*hx_py_u(v) +ey_py_u(hx_px)*v + + # Calculate the complex derivative + PDeriv_complex = ZijN_uV/Hd + Zij * (Hd_uV/Hd) + # Extract the real number for the real/imag components. Pv = np.array(getattr(PDeriv_complex, real_or_imag)) elif adjoint: @@ -192,7 +242,56 @@ class RxMT(Survey.BaseRx): elif self.projType is 'Z2D': raise NotImplementedError('Has not be implement for 2D impedance tensor') elif self.projType is 'Z3D': - raise NotImplementedError('Has not be implement for full 3D impedance tensor') + if self.locs.ndim == 3: + eFLocs = self.locs[:,:,0] + bFLocs = self.locs[:,:,1] + else: + eFLocs = self.locs + bFLocs = self.locs + # Get the projection + Pex = mesh.getInterpolationMat(eFLocs,'Ex') + Pey = mesh.getInterpolationMat(eFLocs,'Ey') + Pbx = mesh.getInterpolationMat(bFLocs,'Fx') + Pby = mesh.getInterpolationMat(bFLocs,'Fy') + # Get the fields at location + # px: x-polaration and py: y-polaration. + ex_px = Pex*f[src,'e_px'] + ey_px = Pey*f[src,'e_px'] + ex_py = Pex*f[src,'e_py'] + ey_py = Pey*f[src,'e_py'] + hx_px = Pbx*f[src,'b_px']/mu_0 + hy_px = Pby*f[src,'b_px']/mu_0 + hx_py = Pbx*f[src,'b_py']/mu_0 + hy_py = Pby*f[src,'b_py']/mu_0 + # Derivatives as lambda functions + ex_px_u = lambda vec: Pex*f._e_pxDeriv_u(src,vec) + ey_px_u = lambda vec: Pey*f._e_pxDeriv_u(src,vec) + ex_py_u = lambda vec: Pex*f._e_pyDeriv_u(src,vec) + ey_py_u = lambda vec: Pey*f._e_pyDeriv_u(src,vec) + hx_px_u = lambda vec: Pbx*f._b_pxDeriv_u(src,vec)/mu_0 + hy_px_u = lambda vec: Pby*f._b_pxDeriv_u(src,vec)/mu_0 + hx_py_u = lambda vec: Pbx*f._b_pyDeriv_u(src,vec)/mu_0 + hy_py_u = lambda vec: Pby*f._b_pyDeriv_u(src,vec)/mu_0 + + # Define the components of the derivative + Hd = (hx_px*hy_py - hx_py*hy_px) + Hd_uV = hx_px_u(hy_py)*v - hx_py*hy_px_u(v) + hx_px*hy_py_u(v) - hx_py_u(hy_px) + # Need to fix this to reflect the adjoint + if 'zxx' in self.rxType: + Zij = ( ex_px*hy_py - ex_py*hy_px)/Hd + ZijN_uV = ex_px_u(hy_py)*v - ex_py*hy_px_u(v) + ex_px*hy_py_u(v) - ex_py_u(hy_px)*v + elif 'zxy' in self.rxType: + Zij = (-ex_px*hx_py + ex_py*hx_px)/Hd + ZijN_uV = -ex_px_u(hx_py)*v + ex_py*hx_px_u(v) - ex_px*hx_py_u(v) + ex_py_u(hx_px)*v + elif 'zyx' in self.rxType: + Zij = ( ey_px*hy_py - ey_py*hy_px)/Hd + ZijN_uV = ey_px_u(hy_py)*v - ey_py*hy_px_u(v) + ey_px*hy_py_u(v) - ey_py_u(hy_px)*v + elif 'zyy' in self.rxType: + Zij = (-ey_px*hx_py + ey_py*hx_px)/Hd + ZijN_uV = -ey_px_u(hx_py)*v + ey_py*hx_px_u(v) - ey_px*hx_py_u(v) +ey_py_u(hx_px)*v + + # Calculate the complex derivative + PDeriv_real = ZijN_uV/Hd + (Hd_uV/Hd)*Zij^T # Extract the data if real_or_imag == 'imag': Pv = 1j*PDeriv_real