diff --git a/docs/examples/FDEM_b.py b/docs/examples/FDEM_b.py index bb3ca167..f0d3c0fc 100644 --- a/docs/examples/FDEM_b.py +++ b/docs/examples/FDEM_b.py @@ -19,9 +19,9 @@ model = Model.LogModel(mesh) x = np.linspace(-10,10,5) XYZ = Utils.ndgrid(x,np.r_[0],np.r_[0]) rxList = EM.FDEM.RxListFDEM(XYZ, 'Ex') -Tx0 = EM.FDEM.TxFDEM(np.r_[0.,0.,0.], 'VMD', 1e2, rxList) +Src0 = EM.FDEM.SrcFDEM(np.r_[0.,0.,0.], 'VMD', 1e2, rxList) -survey = EM.FDEM.SurveyFDEM([Tx0]) +survey = EM.FDEM.SurveyFDEM([Src0]) prb = EM.FDEM.ProblemFDEM_b(model) prb.pair(survey) @@ -31,26 +31,26 @@ sigma = np.ones(mesh.nC)*sig sigma[mesh.gridCC[:,2] > 0] = 1e-8 m = np.log(sigma) -skin = 500*np.sqrt(1/(sig*Tx0.freq)) +skin = 500*np.sqrt(1/(sig*Src0.freq)) print 'The skin depth is: %4.2f m' % skin prb.Solver = Utils.SolverUtils.DSolverWrap(sp.linalg.spsolve, factorize=False, checkAccuracy=True) u = prb.fields(m) -plt.colorbar(mesh.plotImage(np.log10(np.abs(u[Tx0, 'b'].real)), 'Fz')) +plt.colorbar(mesh.plotImage(np.log10(np.abs(u[Src0, 'b'].real)), 'Fz')) -bfz = mesh.r(u[Tx0, 'b'],'F','Fz','M') +bfz = mesh.r(u[Src0, 'b'],'F','Fz','M') x = np.linspace(-55,55,12) XYZ = Utils.ndgrid(x,np.r_[0],np.r_[0]) P = mesh.getInterpolationMat(XYZ, 'Fz') -an = EM.Utils.Ana.FEM.hzAnalyticDipoleF(x, Tx0.freq, sig) +an = EM.Utils.Ana.FEM.hzAnalyticDipoleF(x, Src0.freq, sig) plt.figure(2) -plt.plot(x,np.log10(np.abs(P*np.imag(u[Tx0, 'b'])))) +plt.plot(x,np.log10(np.abs(P*np.imag(u[Src0, 'b'])))) plt.plot(x,np.log10(np.abs(mu_0*np.imag(an))), 'r') plt.xlabel('Distance, m') plt.ylabel('Log10 Response imag($B_z$)') diff --git a/simpegEM/Analytics/FDEM.py b/simpegEM/Analytics/FDEM.py index 8166075a..9abb0a15 100644 --- a/simpegEM/Analytics/FDEM.py +++ b/simpegEM/Analytics/FDEM.py @@ -5,7 +5,8 @@ from scipy.special import erf import matplotlib.pyplot as plt from SimPEG import Utils -def hzAnalyticDipoleF(r, freq, sigma, secondary=True): + +def hzAnalyticDipoleF(r, freq, sigma, secondary=True, mu=mu_0): """ 4.56 in Ward and Hohmann @@ -25,7 +26,7 @@ def hzAnalyticDipoleF(r, freq, sigma, secondary=True): """ r = np.abs(r) - k = np.sqrt(-1j*2.*np.pi*freq*mu_0*sigma) + k = np.sqrt(-1j*2.*np.pi*freq*mu*sigma) m = 1 front = m / (2. * np.pi * (k**2) * (r**5) ) @@ -34,11 +35,14 @@ def hzAnalyticDipoleF(r, freq, sigma, secondary=True): if secondary: hp =-1/(4*np.pi*r**3) - return hz-hp + hz = hz-hp + + if hz.ndim == 1: + hz = Utils.mkvc(hz,2) return hz -def AnalyticMagDipoleWholeSpace(XYZ, txLoc, sig, f, m=1., orientation='X'): +def AnalyticMagDipoleWholeSpace(XYZ, srcLoc, sig, f, moment=1., orientation='X', mu = mu_0): """ Analytical solution for a dipole in a whole-space. @@ -67,15 +71,15 @@ def AnalyticMagDipoleWholeSpace(XYZ, txLoc, sig, f, m=1., orientation='X'): XYZ = Utils.asArray_N_x_Dim(XYZ, 3) - dx = XYZ[:,0]-txLoc[0] - dy = XYZ[:,1]-txLoc[1] - dz = XYZ[:,2]-txLoc[2] + dx = XYZ[:,0]-srcLoc[0] + dy = XYZ[:,1]-srcLoc[1] + dz = XYZ[:,2]-srcLoc[2] r = np.sqrt( dx**2. + dy**2. + dz**2.) - k = np.sqrt( -1j*2.*np.pi*f*mu_0*sig ) + k = np.sqrt( -1j*2.*np.pi*f*mu*sig ) kr = k*r - front = m / (4.*pi * r**3.) * np.exp(-1j*kr) + front = moment / (4.*pi * r**3.) * np.exp(-1j*kr) mid = -kr**2. + 3.*1j*kr + 3. if orientation.upper() == 'X': @@ -93,8 +97,57 @@ def AnalyticMagDipoleWholeSpace(XYZ, txLoc, sig, f, m=1., orientation='X'): Hy = front*( (dy*dz/r**2.) * mid ) Hz = front*( (dz/r)**2. * mid + (kr**2. - 1j*kr - 1.) ) - Bx = mu_0*Hx - By = mu_0*Hy - Bz = mu_0*Hz + Bx = mu*Hx + By = mu*Hy + Bz = mu*Hz + + if Bx.ndim is 1: + Bx = Utils.mkvc(Bx,2) + + if By.ndim is 1: + By = Utils.mkvc(By,2) + + if Bz.ndim is 1: + Bz = Utils.mkvc(Bz,2) + return Bx, By, Bz + +def ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): + XYZ = Utils.asArray_N_x_Dim(XYZ, 3) + + dx = XYZ[:,0]-srcLoc[0] + dy = XYZ[:,1]-srcLoc[1] + dz = XYZ[:,2]-srcLoc[2] + + r = np.sqrt( dx**2. + dy**2. + dz**2.) + k = np.sqrt( -1j*2.*np.pi*f*mu*sig ) + kr = k*r + + front = current * length / (4. * np.pi * sig * r**3) * np.exp(-1j*k*r) + mid = -k**2 * r**2 + 3*1j*k*r + 3 + + # Ex = front*((dx**2 / r**2)*mid + (k**2 * r**2 -1j*k*r)) + # Ey = front*(dx*dy / r**2)*mid + # Ez = front*(dx*dz / r**2)*mid + + if orientation.upper() == 'X': + Ex = front*((dx**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.)) + Ey = front*(dx*dy / r**2)*mid + Ez = front*(dx*dz / r**2)*mid + return Ex, Ey, Ez + + elif orientation.upper() == 'Y': + # x--> y, y--> z, z-->x + Ey = front*((dy**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.)) + Ez = front*(dy*dz / r**2)*mid + Ex = front*(dy*dx / r**2)*mid + return Ex, Ey, Ez + + elif orientation.upper() == 'Z': + # x --> z, y --> x, z --> y + Ez = front*((dz**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.)) + Ex = front*(dz*dx / r**2)*mid + Ey = front*(dz*dy / r**2)*mid + return Ex, Ey, Ez + # return Ey, Ez, Ex \ No newline at end of file diff --git a/simpegEM/Analytics/FDEMcasing.py b/simpegEM/Analytics/FDEMcasing.py new file mode 100644 index 00000000..b6c7acf1 --- /dev/null +++ b/simpegEM/Analytics/FDEMcasing.py @@ -0,0 +1,98 @@ +from SimPEG import Utils, np +from scipy.constants import mu_0, epsilon_0 +from simpegEM.Utils.EMUtils import k + +def getKc(freq,sigma,a,b,mu=mu_0,eps=epsilon_0): + a = float(a) + b = float(b) + # return 1./(2*np.pi) * np.sqrt(b / a) * np.exp(-1j*k(freq,sigma,mu,eps)*(b-a)) + return np.sqrt(b / a) * np.exp(-1j*k(freq,sigma,mu,eps)*(b-a)) + +def _r2(xyz): + return np.sum(xyz**2,1) + +def _getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.): + Kc1 = getKc(freq,sigma[1],a,b,mu[1],eps) + + nobs = obsloc.shape[0] + dxyz = obsloc - np.c_[np.ones(nobs)]*np.r_[srcloc] + + r2 = _r2(dxyz[:,:2]) + sqrtr2z2 = np.sqrt(r2 + dxyz[:,2]**2) + k2 = k(freq,sigma[2],mu[2],eps) + + return Kc1 * moment / (4.*np.pi) *np.exp(-1j*k2*sqrtr2z2) / sqrtr2z2 + + +def _getCasingHertzMagDipoleDeriv_r(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.): + HertzZ = _getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment) + + nobs = obsloc.shape[0] + dxyz = obsloc - np.c_[np.ones(nobs)]*np.r_[srcloc] + + r2 = _r2(dxyz[:,:2]) + sqrtr2z2 = np.sqrt(r2 + dxyz[:,2]**2) + k2 = k(freq,sigma[2],mu[2],eps) + + return -HertzZ * np.sqrt(r2) / sqrtr2z2 * (1j*k2 + 1./ sqrtr2z2) + + +def _getCasingHertzMagDipoleDeriv_z(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.): + HertzZ = _getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment) + + nobs = obsloc.shape[0] + dxyz = obsloc - np.c_[np.ones(nobs)]*np.r_[srcloc] + + r2z2 = _r2(dxyz) + sqrtr2z2 = np.sqrt(r2z2) + k2 = k(freq,sigma[2],mu[2],eps) + + return -HertzZ*dxyz[:,2] /sqrtr2z2 * (1j*k2 + 1./sqrtr2z2) + +def _getCasingHertzMagDipole2Deriv_z_r(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.): + HertzZ = _getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment) + dHertzZdr = _getCasingHertzMagDipoleDeriv_r(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment) + + nobs = obsloc.shape[0] + dxyz = obsloc - np.c_[np.ones(nobs)]*np.r_[srcloc] + + r2 = _r2(dxyz[:,:2]) + r = np.sqrt(r2) + z = dxyz[:,2] + sqrtr2z2 = np.sqrt(r2 + z**2) + k2 = k(freq,sigma[2],mu[2],eps) + + return dHertzZdr*(-z/sqrtr2z2)*(1j*k2+1./sqrtr2z2) + HertzZ*(z*r/sqrtr2z2**3)*(1j*k2 + 2./sqrtr2z2) + +def _getCasingHertzMagDipole2Deriv_z_z(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.): + HertzZ = _getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment) + dHertzZdz = _getCasingHertzMagDipoleDeriv_z(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment) + + nobs = obsloc.shape[0] + dxyz = obsloc - np.c_[np.ones(nobs)]*np.r_[srcloc] + + r2 = _r2(dxyz[:,:2]) + r = np.sqrt(r2) + z = dxyz[:,2] + sqrtr2z2 = np.sqrt(r2 + z**2) + k2 = k(freq,sigma[2],mu[2],eps) + + return (dHertzZdz*z + HertzZ)/sqrtr2z2*(-1j*k2 - 1./sqrtr2z2) + HertzZ*z/sqrtr2z2**3*(1j*k2*z + 2.*z/sqrtr2z2) + +def getCasingEphiMagDipole(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.): + return 1j * omega(freq) * mu * _getCasingHertzMagDipoleDeriv_r(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment) + +def getCasingHrMagDipole(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.): + return _getCasingHertzMagDipole2Deriv_z_r(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment) + +def getCasingHzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.): + d2HertzZdz2 = _getCasingHertzMagDipole2Deriv_z_z(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment) + k2 = k(freq,sigma[2],mu[2],eps) + HertzZ = _getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment) + return d2HertzZdz2 + k2**2 * HertzZ + +def getCasingBrMagDipole(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.): + return mu_0 * getCasingHrMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment) + +def getCasingBzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.): + return mu_0 * getCasingHzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu,eps,moment) \ No newline at end of file diff --git a/simpegEM/Analytics/__init__.py b/simpegEM/Analytics/__init__.py index 5828ba5c..5b7a8851 100644 --- a/simpegEM/Analytics/__init__.py +++ b/simpegEM/Analytics/__init__.py @@ -1,2 +1,3 @@ from TDEM import hzAnalyticDipoleT from FDEM import hzAnalyticDipoleF +from FDEMcasing import * diff --git a/simpegEM/Base.py b/simpegEM/Base.py index d1e6bb08..690b2463 100644 --- a/simpegEM/Base.py +++ b/simpegEM/Base.py @@ -1,16 +1,24 @@ -from SimPEG import Survey, Problem, Utils, Models, np, sp, Solver as SimpegSolver +from SimPEG import Survey, Problem, Utils, Models, Maps, PropMaps, np, sp, Solver as SimpegSolver from scipy.constants import mu_0 +class EMPropMap(Maps.PropMap): + sigma = Maps.Property("Electrical Conductivity", defaultInvProp = True, propertyLink=('rho',Maps.ReciprocalMap)) + mu = Maps.Property("Inverse Magnetic Permeability", defaultVal = mu_0, propertyLink=('mui',Maps.ReciprocalMap)) + + rho = Maps.Property("Electrical Resistivity", propertyLink=('sigma', Maps.ReciprocalMap)) + mui = Maps.Property("Inverse Magnetic Permeability", defaultVal = 1./mu_0, propertyLink=('mu', Maps.ReciprocalMap)) + + class BaseEMProblem(Problem.BaseProblem): def __init__(self, mesh, **kwargs): Problem.BaseProblem.__init__(self, mesh, **kwargs) - solType = None - storeTheseFields = ['e', 'b'] surveyPair = Survey.BaseSurvey dataPair = Survey.Data + + PropMap = EMPropMap Solver = SimpegSolver solverOpts = {} @@ -26,84 +34,112 @@ class BaseEMProblem(Problem.BaseProblem): self.__makeASymmetric = True return self.__makeASymmetric - #################################################### - # Mu Model - #################################################### - @property - def mu(self): - if getattr(self, '_mu', None) is None: - self._mu = mu_0 - return self._mu - @mu.setter - def mu(self, value): - if getattr(self, '_MfMui', None) is not None: - del self._MfMui - if getattr(self, '_MeMu', None) is not None: - del delf._MeMu - if getattr(self, '_MeMuI', None) is not None: - del self._MeMuI - self._mu = value - - #################################################### # Mass Matrices #################################################### @property - def MfMui(self): - if getattr(self, '_MfMui', None) is None: - self._MfMui = self.mesh.getFaceInnerProduct(1/self.mu) - return self._MfMui - - @property - def MeMuI(self): - # TODO: Assuming isotropic mu - if getattr(self, '_MeMuI', None) is None: - self._MeMuI = self.mesh.getEdgeInnerProduct(self.mu, invMat=True) - return self._MeMuI - - @property - def MeMu(self): - #TODO: Assuming isotropic mu - if getattr(self, '_MeMu', None) is None: - self._MeMu = self.mesh.getEdgeInnerProduct(self.mu) - return self._MeMu - + def deleteTheseOnModelUpdate(self): + toDelete = [] + if self.mapping.sigmaMap is not None or self.mapping.rhoMap is not None: + toDelete += ['_MeSigma', '_MeSigmaI','_MfRho','_MfRhoI'] + if self.mapping.muMap is not None or self.mapping.muiMap is not None: + toDelete += ['_MeMu', '_MeMuI','_MfMui','_MfMuiI'] + return toDelete + @property def Me(self): if getattr(self, '_Me', None) is None: self._Me = self.mesh.getEdgeInnerProduct() return self._Me + @property + def Mf(self): + if getattr(self, '_Mf', None) is None: + self._Mf = self.mesh.getFaceInnerProduct() + return self._Mf + + + # ----- Magnetic Permeability ----- # + @property + def MfMui(self): + if getattr(self, '_MfMui', None) is None: + self._MfMui = self.mesh.getFaceInnerProduct(self.curModel.mui) + return self._MfMui + + @property + def MfMuiI(self): + if getattr(self, '_MfMuiI', None) is None: + self._MfMuiI = self.mesh.getFaceInnerProduct(self.curModel.mui, invMat=True) + return self._MfMuiI + + @property + def MeMu(self): + if getattr(self, '_MeMu', None) is None: + self._MeMu = self.mesh.getEdgeInnerProduct(self.curModel.mu) + return self._MeMu + + @property + def MeMuI(self): + if getattr(self, '_MeMuI', None) is None: + self._MeMuI = self.mesh.getEdgeInnerProduct(self.curModel.mu, invMat=True) + return self._MeMuI + + # ----- Electrical Conductivity ----- # + #TODO: hardcoded to sigma as the model @property def MeSigma(self): - #TODO: hardcoded to sigma as the model if getattr(self, '_MeSigma', None) is None: - sigma = self.curModel.transform - self._MeSigma = self.mesh.getEdgeInnerProduct(sigma) + self._MeSigma = self.mesh.getEdgeInnerProduct(self.curModel.sigma) return self._MeSigma + # TODO: This should take a vector + def MeSigmaDeriv(self, u): + """ + Deriv of MeSigma wrt sigma + """ + return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma)(u) * self.curModel.sigmaDeriv + + @property def MeSigmaI(self): - #TODO: hardcoded to sigma as the model if getattr(self, '_MeSigmaI', None) is None: - sigma = self.curModel.transform - self._MeSigmaI = self.mesh.getEdgeInnerProduct(sigma, invMat=True) + self._MeSigmaI = self.mesh.getEdgeInnerProduct(self.curModel.sigma, invMat=True) return self._MeSigmaI + # TODO: This should take a vector + def MeSigmaIDeriv(self, u): + """ + Deriv of MeSigma wrt sigma + """ + # TODO: only works for diagonal tensors. getEdgeInnerProductDeriv, invMat=True should be implemented in SimPEG + + dMeSigmaI_dI = -self.MeSigmaI**2 + dMe_dsig = self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma)(u) + dsig_dm = self.curModel.sigmaDeriv + return dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm)) + # return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma, invMat=True)(u) + + @property - def MfSigmai(self): - #TODO: hardcoded to sigma as the model - #TODO: hardcoded to sigma diagonal - if getattr(self, '_MfSigmai', None) is None: - sigma = self.curModel.transform - self._MfSigmai = self.mesh.getFaceInnerProduct(1/sigma) - return self._MfSigmai + def MfRho(self): + if getattr(self, '_MfRho', None) is None: + self._MfRho = self.mesh.getFaceInnerProduct(self.curModel.rho) + return self._MfRho - deleteTheseOnModelUpdate = ['_MeSigma', '_MeSigmaI','_MfSigmai'] + # TODO: This should take a vector + def MfRhoDeriv(self,u): + return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) * (-Utils.sdiag(self.curModel.rho**2) * self.curModel.sigmaDeriv) + # self.curModel.rhoDeriv - def fields(self, m): - self.curModel = m - F = self.forward(m, self.getRHS, self.calcFields) - return F + @property + def MfRhoI(self): + if getattr(self, '_MfRhoI', None) is None: + self._MfRhoI = self.mesh.getFaceInnerProduct(self.curModel.rho, invMat=True) + return self._MfRhoI + + # TODO: This isn't going to work yet + # TODO: This should take a vector + def MfRhoIDeriv(self,u): + return self.mesh.getFaceInnerProductDeriv(self.curModel.rho, invMat=True)(u) * self.curModel.rhoDeriv diff --git a/simpegEM/Examples/CylInversion.py b/simpegEM/Examples/CylInversion.py index c944d6ef..4abd51fc 100644 --- a/simpegEM/Examples/CylInversion.py +++ b/simpegEM/Examples/CylInversion.py @@ -36,8 +36,8 @@ if plotIt: rxOffset=1e-3 rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 30]]), np.logspace(-5,-3, 31), 'bz') -tx = EM.TDEM.TxTDEM(np.array([0., 0., 80]), 'VMD_MVP', [rx]) -survey = EM.TDEM.SurveyTDEM([tx]) +src = EM.TDEM.SrcTDEM_VMD_MVP([rx], np.array([0., 0., 80])) +survey = EM.TDEM.SurveyTDEM([src]) prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping) prb.Solver = SolverLU diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index 1a45d139..c54137fd 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -1,108 +1,10 @@ from SimPEG import Survey, Problem, Utils, np, sp, Solver as SimpegSolver from scipy.constants import mu_0 -from SurveyFDEM import SurveyFDEM, FieldsFDEM -from simpegEM import Sources +from SurveyFDEM import SurveyFDEM +from FieldsFDEM import FieldsFDEM, FieldsFDEM_e, FieldsFDEM_b, FieldsFDEM_h, FieldsFDEM_j from simpegEM.Base import BaseEMProblem +from simpegEM.Utils.EMUtils import omega -def omega(freq): - """Change frequency to angular frequency, omega""" - return 2.*np.pi*freq - -def getSource(self,freq): - """ - :param float freq: Frequency - :rtype: numpy.ndarray (nE, nTx) - :return: RHS - """ - Txs = self.survey.getTransmitters(freq) - rhs = range(len(Txs)) - - solType = self.solType - - if solType == 'e' or solType == 'b': - gridEJx = self.mesh.gridEx - gridEJy = self.mesh.gridEy - gridEJz = self.mesh.gridEz - nEJ = self.mesh.nE - - gridBHx = self.mesh.gridFx - gridBHy = self.mesh.gridFy - gridBHz = self.mesh.gridFz - nBH = self.mesh.nF - - - C = self.mesh.edgeCurl - mui = self.MfMui - - elif solType == 'h' or solType == 'j': - gridEJx = self.mesh.gridFx - gridEJy = self.mesh.gridFy - gridEJz = self.mesh.gridFz - nEJ = self.mesh.nF - - gridBHx = self.mesh.gridEx - gridBHy = self.mesh.gridEy - gridBHz = self.mesh.gridEz - nBH = self.mesh.nE - - C = self.mesh.edgeCurl.T - mui = self.MeMuI - - else: - NotImplementedError('Only E or F sources') - - for i, tx in enumerate(Txs): - if self.mesh._meshType is 'CYL': - if self.mesh.isSymmetric: - if tx.txType == 'VMD': - SRC = Sources.MagneticDipoleVectorPotential(tx.loc, gridEJy, 'y') - elif tx.txType =='CircularLoop': - SRC = Sources.MagneticLoopVectorPotential(tx.loc, gridEJy, 'y', tx.radius) - else: - raise NotImplementedError('Only VMD and CircularLoop') - else: - raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') - - elif self.mesh._meshType is 'TENSOR': - - if tx.txType == 'VMD': - src = Sources.MagneticDipoleVectorPotential - SRCx = src(tx.loc, gridEJx, 'x') - SRCy = src(tx.loc, gridEJy, 'y') - SRCz = src(tx.loc, gridEJz, 'z') - - elif tx.txType == 'VMD_B': - src = Sources.MagneticDipoleFields - SRCx = src(tx.loc, gridBHx, 'x') - SRCy = src(tx.loc, gridBHy, 'y') - SRCz = src(tx.loc, gridBHz, 'z') - - elif tx.txType == 'CircularLoop': - src = Sources.MagneticLoopVectorPotential - SRCx = src(tx.loc, gridEJx, 'x', tx.radius) - SRCy = src(tx.loc, gridEJy, 'y', tx.radius) - SRCz = src(tx.loc, gridEJz, 'z', tx.radius) - else: - - raise NotImplemented('%s txType is not implemented' % tx.txType) - SRC = np.concatenate((SRCx, SRCy, SRCz)) - - else: - raise Exception('Unknown mesh for VMD') - - rhs[i] = SRC - - # b-forumlation - if tx.txType == 'VMD_B': - b_0 = np.concatenate(rhs).reshape((nBH, len(Txs)), order='E') - else: - a = np.concatenate(rhs).reshape((nEJ, len(Txs)), order='F') - b_0 = C*a - - if solType == 'b' or solType == 'h': - return b_0 - elif solType == 'e' or solType == 'j': - return C.T*mui*b_0 class BaseFDEMProblem(BaseEMProblem): """ @@ -110,58 +12,71 @@ class BaseFDEMProblem(BaseEMProblem): .. math:: - \\nabla \\times \\vec{E} + i \\omega \\vec{B} = 0 \\\\ - \\nabla \\times \\mu^{-1} \\vec{B} - \\sigma \\vec{E} = \\vec{J_s} + \\nabla \\times \\vec{E} + i \\omega \\vec{B} = \\vec{S_m} \\\\ + \\nabla \\times \\mu^{-1} \\vec{B} - \\sigma \\vec{E} = \\vec{S_e} """ surveyPair = SurveyFDEM + fieldsPair = FieldsFDEM - def forward(self, m, RHS, CalcFields): - - F = FieldsFDEM(self.mesh, self.survey) + def fields(self, m=None): + self.curModel = m + F = self.fieldsPair(self.mesh, self.survey) for freq in self.survey.freqs: A = self.getA(freq) - rhs = RHS(freq) + rhs = self.getRHS(freq) Ainv = self.Solver(A, **self.solverOpts) sol = Ainv * rhs - for fieldType in self.storeTheseFields: - Txs = self.survey.getTransmitters(freq) - F[Txs, fieldType] = CalcFields(sol, freq, fieldType) + Srcs = self.survey.getSrcByFreq(freq) + ftype = self._fieldType + 'Solution' + F[Srcs, ftype] = sol return F - def Jvec(self, m, v, u=None): - if u is None: - u = self.fields(m) + def Jvec(self, m, v, f=None): + if f is None: + f = self.fields(m) self.curModel = m Jv = self.dataPair(self.survey) for freq in self.survey.freqs: - A = self.getA(freq) - Ainv = self.Solver(A, **self.solverOpts) + dA_du = self.getA(freq) # + dA_duI = self.Solver(dA_du, **self.solverOpts) - for tx in self.survey.getTransmitters(freq): - u_tx = u[tx, self.solType] - w = self.getADeriv(freq, u_tx, v) - Ainvw = Ainv * w - for rx in tx.rxList: - fAinvw = self.calcFields(Ainvw, freq, rx.projField) - P = lambda v: rx.projectFieldsDeriv(tx, self.mesh, u, v) + for src in self.survey.getSrcByFreq(freq): + ftype = self._fieldType + 'Solution' + u_src = f[src, ftype] + dA_dm = self.getADeriv_m(freq, u_src, v) + dRHS_dm = self.getRHSDeriv_m(src, v) + if dRHS_dm is None: + du_dm = dA_duI * ( - dA_dm ) + else: + du_dm = dA_duI * ( - dA_dm + dRHS_dm ) + for rx in src.rxList: + # df_duFun = u.deriv_u(rx.fieldsUsed, m) + df_duFun = getattr(f, '_%sDeriv_u'%rx.projField, None) + df_du = df_duFun(src, du_dm, adjoint=False) + if df_du is not None: + du_dm = df_du - df_dm = self.calcFieldsDeriv(u_tx, freq, rx.projField, v) - if df_dm is None: - Jv[tx, rx] = - P(fAinvw) - else: - Jv[tx, rx] = - P(fAinvw) + P(df_dm) + df_dmFun = getattr(f, '_%sDeriv_m'%rx.projField, None) + df_dm = df_dmFun(src, v, adjoint=False) + if df_dm is not None: + du_dm += df_dm + + P = lambda v: rx.projectFieldsDeriv(src, self.mesh, f, v) # wrt u, also have wrt m + + + Jv[src, rx] = P(du_dm) return Utils.mkvc(Jv) - def Jtvec(self, m, v, u=None): - if u is None: - u = self.fields(m) + def Jtvec(self, m, v, f=None): + if f is None: + f = self.fields(m) self.curModel = m @@ -169,39 +84,73 @@ class BaseFDEMProblem(BaseEMProblem): if not isinstance(v, self.dataPair): v = self.dataPair(self.survey, v) - Jtv = np.zeros(self.mapping.nP) + Jtv = np.zeros(m.size) for freq in self.survey.freqs: AT = self.getA(freq).T ATinv = self.Solver(AT, **self.solverOpts) - for tx in self.survey.getTransmitters(freq): - u_tx = u[tx, self.solType] + for src in self.survey.getSrcByFreq(freq): + ftype = self._fieldType + 'Solution' + u_src = f[src, ftype] - for rx in tx.rxList: - PTv = rx.projectFieldsDeriv(tx, self.mesh, u, v[tx, rx], adjoint=True) - fPTv = self.calcFields(PTv, freq, rx.projField, adjoint=True) + for rx in src.rxList: + PTv = rx.projectFieldsDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt u, need possibility wrt m - w = ATinv * fPTv - Jtv_rx = - self.getADeriv(freq, u_tx, w, adjoint=True) + df_duTFun = getattr(f, '_%sDeriv_u'%rx.projField, None) + df_duT = df_duTFun(src, PTv, adjoint=True) + if df_duT is not None: + dA_duIT = ATinv * df_duT + else: + dA_duIT = ATinv * PTv - df_dm = self.calcFieldsDeriv(u_tx, freq, rx.projField, PTv, adjoint=True) + dA_dmT = self.getADeriv_m(freq, u_src, dA_duIT, adjoint=True) - if df_dm is not None: - Jtv_rx += df_dm + dRHS_dmT = self.getRHSDeriv_m(src, dA_duIT, adjoint=True) + + if dRHS_dmT is None: + du_dmT = - dA_dmT + else: + du_dmT = -dA_dmT + dRHS_dmT + + df_dmFun = getattr(f, '_%sDeriv_m'%rx.projField, None) + dfT_dm = df_dmFun(src, PTv, adjoint=True) + if dfT_dm is not None: + du_dmT += dfT_dm real_or_imag = rx.projComp if real_or_imag == 'real': - Jtv += Jtv_rx.real + Jtv += du_dmT.real elif real_or_imag == 'imag': - Jtv += - Jtv_rx.real + Jtv += - du_dmT.real else: raise Exception('Must be real or imag') return Jtv - def getSource(self,freq): - return self.getSource(freq) + def getSourceTerm(self, freq): + """ + :param float freq: Frequency + :rtype: numpy.ndarray (nE or nF, nSrc) + :return: RHS + """ + Srcs = self.survey.getSrcByFreq(freq) + if self._eqLocs is 'FE': + S_m = np.zeros((self.mesh.nF,len(Srcs)), dtype=complex) + S_e = np.zeros((self.mesh.nE,len(Srcs)), dtype=complex) + elif self._eqLocs is 'EF': + S_m = np.zeros((self.mesh.nE,len(Srcs)), dtype=complex) + S_e = np.zeros((self.mesh.nF,len(Srcs)), dtype=complex) + + for i, src in enumerate(Srcs): + smi, sei = src.eval(self) + if smi is not None: + S_m[:,i] = Utils.mkvc(smi) + if sei is not None: + S_e[:,i] = Utils.mkvc(sei) + + return S_m, S_e + ########################################################################################## ################################ E-B Formulation ######################################### @@ -225,10 +174,13 @@ class ProblemFDEM_e(BaseFDEMProblem): """ - solType = 'e' - def __init__(self, model, **kwargs): - BaseFDEMProblem.__init__(self, model, **kwargs) + _fieldType = 'e' + _eqLocs = 'FE' + fieldsPair = FieldsFDEM_e + + def __init__(self, mesh, **kwargs): + BaseFDEMProblem.__init__(self, mesh, **kwargs) def getA(self, freq): """ @@ -236,62 +188,77 @@ class ProblemFDEM_e(BaseFDEMProblem): :rtype: scipy.sparse.csr_matrix :return: A """ - mui = self.MfMui - sig = self.MeSigma + MfMui = self.MfMui + MeSigma = self.MeSigma C = self.mesh.edgeCurl - return C.T*mui*C + 1j*omega(freq)*sig + return C.T*MfMui*C + 1j*omega(freq)*MeSigma - def getADeriv(self, freq, u, v, adjoint=False): - sig = self.curModel.transform - dsig_dm = self.curModel.transformDeriv - dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig)(u) + def getADeriv_m(self, freq, u, v, adjoint=False): + dsig_dm = self.curModel.sigmaDeriv + dMe_dsig = self.MeSigmaDeriv(u) if adjoint: - return 1j * omega(freq) * ( dsig_dm.T * ( dMe_dsig.T * v ) ) + return 1j * omega(freq) * ( dMe_dsig.T * v ) - return 1j * omega(freq) * ( dMe_dsig * ( dsig_dm * v ) ) + return 1j * omega(freq) * ( dMe_dsig * v ) def getRHS(self, freq): """ :param float freq: Frequency - :rtype: numpy.ndarray (nE, nTx) + :rtype: numpy.ndarray (nE, nSrc) :return: RHS """ - j_s = getSource(self,freq) - return -1j*omega(freq)*j_s + S_m, S_e = self.getSourceTerm(freq) + C = self.mesh.edgeCurl + MfMui = self.MfMui - def calcFields(self, sol, freq, fieldType, adjoint=False): - e = sol - if fieldType == 'e': - return e - elif fieldType == 'b': - if not adjoint: - b = -(1./(1j*omega(freq))) * ( self.mesh.edgeCurl * e ) + RHS = C.T * (MfMui * S_m) -1j*omega(freq)*S_e + + return RHS + + def getRHSDeriv_m(self, src, v, adjoint=False): + C = self.mesh.edgeCurl + MfMui = self.MfMui + S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint) + + if adjoint: + dRHS = MfMui * (C * v) + S_mDerivv = S_mDeriv(dRHS) + S_eDerivv = S_eDeriv(v) + if S_mDerivv is not None and S_eDerivv is not None: + return S_mDerivv - 1j*omega(freq)*S_eDerivv + elif S_mDerivv is not None: + return S_mDerivv + elif S_eDerivv is not None: + return - 1j*omega(freq)*S_eDerivv else: - b = -(1./(1j*omega(freq))) * ( self.mesh.edgeCurl.T * e ) - return b - raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) + return None + else: + S_mDerivv, S_eDerivv = S_mDeriv(v), S_eDeriv(v) - def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False): - e = sol - if fieldType == 'e': - return None - elif fieldType == 'b': - return None - raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) + if S_mDerivv is not None and S_eDerivv is not None: + return C.T * (MfMui * S_mDerivv) -1j*omega(freq)*S_eDerivv + elif S_mDerivv is not None: + return C.T * (MfMui * S_mDerivv) + elif S_eDerivv is not None: + return -1j*omega(freq)*S_eDerivv + else: + return None class ProblemFDEM_b(BaseFDEMProblem): """ Solving for b! """ - solType = 'b' + _fieldType = 'b' + _eqLocs = 'FE' + fieldsPair = FieldsFDEM_b - def __init__(self, model, **kwargs): - BaseFDEMProblem.__init__(self, model, **kwargs) + def __init__(self, mesh, **kwargs): + BaseFDEMProblem.__init__(self, mesh, **kwargs) def getA(self, freq): """ @@ -299,87 +266,100 @@ class ProblemFDEM_b(BaseFDEMProblem): :rtype: scipy.sparse.csr_matrix :return: A """ - mui = self.MfMui - sigI = self.MeSigmaI + MfMui = self.MfMui + MeSigmaI = self.MeSigmaI C = self.mesh.edgeCurl iomega = 1j * omega(freq) * sp.eye(self.mesh.nF) - A = C*sigI*C.T*mui + iomega + A = C*MeSigmaI*C.T*MfMui + iomega if self._makeASymmetric is True: - return mui.T*A - return A + return MfMui.T*A + return A - def getADeriv(self, freq, u, v, adjoint=False): + def getADeriv_m(self, freq, u, v, adjoint=False): - mui = self.MfMui + MfMui = self.MfMui C = self.mesh.edgeCurl - sig = self.curModel.transform - dsig_dm = self.curModel.transformDeriv - #TODO: This only works if diagonal (no tensors)... - dMeSigmaI_dI = - self.MeSigmaI**2 + MeSigmaIDeriv = self.MeSigmaIDeriv + vec = C.T*(MfMui*u) - vec = (C.T*(mui*u)) - dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig)(vec) + MeSigmaIDeriv = MeSigmaIDeriv(vec) if adjoint: if self._makeASymmetric is True: - v = mui * v - return dsig_dm.T * ( dMe_dsig.T * ( dMeSigmaI_dI.T * ( C.T * v ) ) ) + v = MfMui * v + return MeSigmaIDeriv.T * (C.T * v) if self._makeASymmetric is True: - return mui.T * ( C * ( dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm * v ) ) ) ) - return C * ( dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm * v ) ) ) + return MfMui.T * ( C * ( MeSigmaIDeriv * v ) ) + return C * ( MeSigmaIDeriv * v ) def getRHS(self, freq): """ :param float freq: Frequency - :rtype: numpy.ndarray (nE, nTx) + :rtype: numpy.ndarray (nE, nSrc) :return: RHS """ - - b_0 = getSource(self,freq) - rhs = -1j*omega(freq)*b_0 + S_m, S_e = self.getSourceTerm(freq) + C = self.mesh.edgeCurl + MeSigmaI = self.MeSigmaI + + RHS = S_m + C * ( MeSigmaI * S_e ) + if self._makeASymmetric is True: - mui = self.MfMui - return mui.T*rhs - return rhs + MfMui = self.MfMui + return MfMui.T*RHS - def calcFields(self, sol, freq, fieldType, adjoint=False): - b = sol - if fieldType == 'e': + return RHS + + def getRHSDeriv_m(self, src, v, adjoint=False): + C = self.mesh.edgeCurl + S_m, S_e = src.eval(self) + MfMui = self.MfMui + + if self._makeASymmetric and adjoint: + v = self.MfMui * v + + if S_e is not None: + MeSigmaIDeriv = self.MeSigmaIDeriv(Utils.mkvc(S_e)) if not adjoint: - e = self.MeSigmaI * ( self.mesh.edgeCurl.T * ( self.MfMui * b ) ) - else: - e = self.MfMui.T * ( self.mesh.edgeCurl * ( self.MeSigmaI.T * b ) ) - return e - elif fieldType == 'b': - return b - raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) + RHSderiv = C * (MeSigmaIDeriv * v) + elif adjoint: + RHSderiv = MeSigmaIDeriv.T * (C.T * v) + else: + RHSderiv = None - def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False): - b = sol - if fieldType == 'e': - sig = self.curModel.transform - dsig_dm = self.curModel.transformDeriv - - C = self.mesh.edgeCurl - mui = self.MfMui - - #TODO: This only works if diagonal (no tensors)... - dMeSigmaI_dI = - self.MeSigmaI**2 - - vec = C.T * ( mui * b ) - dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig)(vec) + S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint) + S_mDeriv, S_eDeriv = S_mDeriv(v), S_eDeriv(v) + if S_mDeriv is not None and S_eDeriv is not None: if not adjoint: - return dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm * v ) ) - else: - return dsig_dm.T * ( dMe_dsig.T * ( dMeSigmaI_dI.T * v ) ) - elif fieldType == 'b': - return None - raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) + SrcDeriv = S_mDeriv + C * (self.MeSigmaI * S_eDeriv) + elif adjoint: + SrcDeriv = S_mDeriv + Self.MeSigmaI.T * ( C.T * S_eDeriv) + elif S_mDeriv is not None: + SrcDeriv = S_mDeriv + elif S_eDeriv is not None: + if not adjoint: + SrcDeriv = C * (self.MeSigmaI * S_eDeriv) + elif adjoint: + SrcDeriv = self.MeSigmaI.T * ( C.T * S_eDeriv) + else: + SrcDeriv = None + + if RHSderiv is not None and SrcDeriv is not None: + RHSderiv += SrcDeriv + elif SrcDeriv is not None: + RHSderiv = SrcDeriv + + if RHSderiv is not None: + if self._makeASymmetric is True and not adjoint: + return MfMui.T * RHSderiv + + return RHSderiv + ########################################################################################## @@ -401,7 +381,7 @@ class ProblemFDEM_j(BaseFDEMProblem): .. math:: \\nabla \\times ( \\mu^{-1} \\nabla \\times \\sigma^{-1} \\vec{J} ) + i\\omega \\vec{J} = - i\\omega\\vec{J_s} - + We discretize this to: .. math:: @@ -412,15 +392,16 @@ class ProblemFDEM_j(BaseFDEMProblem): """ - solType = 'j' - storeTheseFields = ['j','h'] + _fieldType = 'j' + _eqLocs = 'EF' + fieldsPair = FieldsFDEM_j - def __init__(self, model, **kwargs): - BaseFDEMProblem.__init__(self, model, **kwargs) + def __init__(self, mesh, **kwargs): + BaseFDEMProblem.__init__(self, mesh, **kwargs) def getA(self, freq): """ - Here, we form the operator \(\\mathbf{A}\) to solce + Here, we form the operator \(\\mathbf{A}\) to solce .. math:: \\mathbf{A} = \\mathbf{C} \\mathbf{M^e_{mu^{-1}}} \\mathbf{C^T} \\mathbf{M^f_{\\sigma^{-1}}} + i\\omega @@ -430,96 +411,98 @@ class ProblemFDEM_j(BaseFDEMProblem): """ MeMuI = self.MeMuI - MfSigi = self.MfSigmai + MfRho = self.MfRho C = self.mesh.edgeCurl iomega = 1j * omega(freq) * sp.eye(self.mesh.nF) - A = C * MeMuI * C.T * MfSigi + iomega + A = C * MeMuI * C.T * MfRho + iomega if self._makeASymmetric is True: - return MfSigi.T*A - return A + return MfRho.T*A + return A - def getADeriv(self, freq, u, v, adjoint=False): + def getADeriv_m(self, freq, u, v, adjoint=False): """ In this case, we assume that electrical conductivity, \(\\sigma\) is the physical property of interest (i.e. \(\sigma\) = model.transform). Then we want - .. math:: + .. math:: \\frac{\mathbf{A(\\sigma)} \mathbf{v}}{d \\mathbf{m}} &= \\mathbf{C} \\mathbf{M^e_{mu^{-1}}} \\mathbf{C^T} \\frac{d \\mathbf{M^f_{\\sigma^{-1}}}}{d \\mathbf{m}} &= \\mathbf{C} \\mathbf{M^e_{mu}^{-1}} \\mathbf{C^T} \\frac{d \\mathbf{M^f_{\\sigma^{-1}}}}{d \\mathbf{\\sigma^{-1}}} \\frac{d \\mathbf{\\sigma^{-1}}}{d \\mathbf{\\sigma}} \\frac{d \\mathbf{\\sigma}}{d \\mathbf{m}} """ MeMuI = self.MeMuI - MfSigi = self.MfSigmai + MfRho = self.MfRho C = self.mesh.edgeCurl - sig = self.curModel.transform - sigi = 1/sig - dsig_dm = self.curModel.transformDeriv - dsigi_dsig = -Utils.sdiag(sigi)**2 - dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(u) + MfRhoDeriv_m = self.MfRhoDeriv(u) if adjoint: if self._makeASymmetric is True: - v = MfSigi * v - return dsig_dm.T * ( dsigi_dsig.T *( dMf_dsigi.T * ( C * ( MeMuI.T * ( C.T * v ) ) ) ) ) + v = MfRho * v + return MfRhoDeriv_m.T * (C * (MeMuI.T * (C.T * v))) - if self._makeASymmetric is True: - return MfSigi.T * ( C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) ) ) - return C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) ) + if self._makeASymmetric is True: + return MfRho.T * (C * ( MeMuI * (C.T * (MfRhoDeriv_m * v) ))) + return C * (MeMuI * (C.T * (MfRhoDeriv_m * v))) def getRHS(self, freq): """ :param float freq: Frequency - :rtype: numpy.ndarray (nE, nTx) + :rtype: numpy.ndarray (nE, nSrc) :return: RHS """ - j_s = getSource(self,freq) - rhs = -1j*omega(freq)*j_s + S_m, S_e = self.getSourceTerm(freq) + C = self.mesh.edgeCurl + MeMuI = self.MeMuI + + + RHS = C * (MeMuI * S_m) - 1j * omega(freq) * S_e if self._makeASymmetric is True: - MfSigi = self.MfSigmai - return MfSigi.T*rhs - return rhs + MfRho = self.MfRho + return MfRho.T*RHS - def calcFields(self, sol, freq, fieldType, adjoint=False): - j = sol - if fieldType == 'j': - return j - elif fieldType == 'h': - MeMuI = self.MeMuI - C = self.mesh.edgeCurl - MfSigi = self.MfSigmai - if not adjoint: - h = -(1./(1j*omega(freq))) * MeMuI * ( C.T * ( MfSigi * j ) ) + return RHS + + def getRHSDeriv_m(self, src, v, adjoint=False): + C = self.mesh.edgeCurl + MeMuI = self.MeMuI + S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint) + + if adjoint: + if self._makeASymmetric: + MfRho = self.MfRho + v = MfRho*v + S_mDerivv = S_mDeriv(MeMuI.T * (C.T * v)) + S_eDerivv = S_eDeriv(v) + if S_mDerivv is not None and S_eDerivv is not None: + return S_mDerivv - 1j*omega(freq)*S_eDerivv + elif S_mDerivv is not None: + return S_mDerivv + elif S_eDerivv is not None: + return - 1j*omega(freq)*S_eDerivv else: - h = -(1./(1j*omega(freq))) * MfSigi.T * ( C * ( MeMuI.T * j ) ) - return h - raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) + return None + else: + S_mDerivv, S_eDerivv = S_mDeriv(v), S_eDeriv(v) - def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False): - j = sol - if fieldType == 'j': - return None - elif fieldType == 'h': - MeMuI = self.MeMuI - C = self.mesh.edgeCurl - sig = self.curModel.transform - sigi = 1/sig - dsig_dm = self.curModel.transformDeriv - dsigi_dsig = -Utils.sdiag(sigi)**2 - dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(j) - sigi = self.MfSigmai - if not adjoint: - return -(1./(1j*omega(freq))) * MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) + if S_mDerivv is not None and S_eDerivv is not None: + RHSDeriv = C * (MeMuI * S_mDerivv) - 1j * omega(freq) * S_eDerivv + elif S_mDerivv is not None: + RHSDeriv = C * (MeMuI * S_mDerivv) + elif S_eDerivv is not None: + RHSDeriv = - 1j * omega(freq) * S_eDerivv else: - return -(1./(1j*omega(freq))) * dsig_dm.T * ( dsigi_dsig.T * ( dMf_dsigi.T * ( C * ( MeMuI.T * v ) ) ) ) - raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) + return None + + if self._makeASymmetric: + MfRho = self.MfRho + return MfRho.T * RHSDeriv + return RHSDeriv -# Solving for h! - using primary- secondary approach class ProblemFDEM_h(BaseFDEMProblem): """ Using the H-J formulation of Maxwell's equations @@ -545,11 +528,12 @@ class ProblemFDEM_h(BaseFDEMProblem): """ - solType = 'h' - storeTheseFields = ['j','h'] + _fieldType = 'h' + _eqLocs = 'EF' + fieldsPair = FieldsFDEM_h - def __init__(self, model, **kwargs): - BaseFDEMProblem.__init__(self, model, **kwargs) + def __init__(self, mesh, **kwargs): + BaseFDEMProblem.__init__(self, mesh, **kwargs) def getA(self, freq): """ @@ -559,47 +543,64 @@ class ProblemFDEM_h(BaseFDEMProblem): """ MeMu = self.MeMu - MfSigi = self.MfSigmai + MfRho = self.MfRho C = self.mesh.edgeCurl - return C.T * MfSigi * C + 1j*omega(freq)*MeMu + return C.T * MfRho * C + 1j*omega(freq)*MeMu - def getADeriv(self, freq, u, v, adjoint=False): + def getADeriv_m(self, freq, u, v, adjoint=False): MeMu = self.MeMu C = self.mesh.edgeCurl - sig = self.curModel.transform - sigi = 1/sig - dsig_dm = self.curModel.transformDeriv - dsigi_dsig = -Utils.sdiag(sigi)**2 - - dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(C*u) + MfRhoDeriv_m = self.MfRhoDeriv(C*u) if adjoint: - return (dsig_dm.T * (dsigi_dsig.T * (dMf_dsigi.T * (C * v)))) - return (C.T * (dMf_dsigi * (dsigi_dsig * (dsig_dm * v)))) - - + return MfRhoDeriv_m.T * (C * v) + return C.T * (MfRhoDeriv_m * v) def getRHS(self, freq): """ :param float freq: Frequency - :rtype: numpy.ndarray (nE, nTx) + :rtype: numpy.ndarray (nE, nSrc) :return: RHS """ - b_0 = getSource(self,freq) - return -1j*omega(freq)*b_0 - - def calcFields(self, sol, freq, fieldType, adjoint=False): - h = sol - if fieldType == 'j': - C = self.mesh.edgeCurl - if adjoint: - return C.T*h - return C*h - elif fieldType == 'h': - return h - raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) - def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False): - return None \ No newline at end of file + S_m, S_e = self.getSourceTerm(freq) + C = self.mesh.edgeCurl + MfRho = self.MfRho + + RHS = S_m + C.T * ( MfRho * S_e ) + + return RHS + + def getRHSDeriv_m(self, src, v, adjoint=False): + _, S_e = src.eval(self) + C = self.mesh.edgeCurl + MfRho = self.MfRho + + RHSDeriv = None + + if S_e is not None: + MfRhoDeriv = self.MfRhoDeriv(S_e) + if not adjoint: + RHSDeriv = C.T * (MfRhoDeriv * v) + elif adjoint: + RHSDeriv = MfRhoDeriv.T * (C * v) + + S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint) + + S_mDeriv = S_mDeriv(v) + S_eDeriv = S_eDeriv(v) + if S_mDeriv is not None: + if RHSDeriv is not None: + RHSDeriv += S_mDeriv(v) + else: + RHSDeriv = S_mDeriv(v) + if S_eDeriv is not None: + if RHSDeriv is not None: + RHSDeriv += C.T * (MfRho * S_e) + else: + RHSDeriv = C.T * (MfRho * S_e) + + return RHSDeriv + diff --git a/simpegEM/FDEM/FieldsFDEM.py b/simpegEM/FDEM/FieldsFDEM.py new file mode 100644 index 00000000..99866bc3 --- /dev/null +++ b/simpegEM/FDEM/FieldsFDEM.py @@ -0,0 +1,373 @@ +from SimPEG import Survey, Problem, Utils, np, sp +from simpegEM.Utils.EMUtils import omega + + +class FieldsFDEM(Problem.Fields): + """Fancy Field Storage for a FDEM survey.""" + knownFields = {} + dtype = complex + +class FieldsFDEM_e(FieldsFDEM): + knownFields = {'eSolution':'E'} + aliasFields = { + 'e' : ['eSolution','E','_e'], + 'ePrimary' : ['eSolution','E','_ePrimary'], + 'eSecondary' : ['eSolution','E','_eSecondary'], + 'b' : ['eSolution','F','_b'], + 'bPrimary' : ['eSolution','F','_bPrimary'], + 'bSecondary' : ['eSolution','F','_bSecondary'] + } + + def __init__(self,mesh,survey,**kwargs): + FieldsFDEM.__init__(self,mesh,survey,**kwargs) + + def startup(self): + self.prob = self.survey.prob + self._edgeCurl = self.survey.prob.mesh.edgeCurl + + def _ePrimary(self, eSolution, srcList): + ePrimary = np.zeros_like(eSolution) + for i, src in enumerate(srcList): + ep = src.ePrimary(self.prob) + if ep is not None: + ePrimary[:,i] = ep + return ePrimary + + def _eSecondary(self, eSolution, srcList): + return eSolution + + def _e(self, eSolution, srcList): + return self._ePrimary(eSolution,srcList) + self._eSecondary(eSolution,srcList) + + def _eDeriv_u(self, src, v, adjoint = False): + return None + + def _eDeriv_m(self, src, v, adjoint = False): + # assuming primary does not depend on the model + return None + + def _bPrimary(self, eSolution, srcList): + bPrimary = np.zeros([self._edgeCurl.shape[0],eSolution.shape[1]],dtype = complex) + for i, src in enumerate(srcList): + bp = src.bPrimary(self.prob) + if bp is not None: + bPrimary[:,i] += bp + return bPrimary + + def _bSecondary(self, eSolution, srcList): + C = self._edgeCurl + b = (C * eSolution) + for i, src in enumerate(srcList): + b[:,i] *= - 1./(1j*omega(src.freq)) + S_m, _ = src.eval(self.prob) + if S_m is not None: + b[:,i] += 1./(1j*omega(src.freq)) * S_m + return b + + def _bSecondaryDeriv_u(self, src, v, adjoint = False): + C = self._edgeCurl + if adjoint: + return - 1./(1j*omega(src.freq)) * (C.T * v) + return - 1./(1j*omega(src.freq)) * (C * v) + + def _bSecondaryDeriv_m(self, src, v, adjoint = False): + S_mDeriv, _ = src.evalDeriv(self.prob, adjoint) + S_mDeriv = S_mDeriv(v) + if S_mDeriv is not None: + return 1./(1j * omega(src.freq)) * S_mDeriv + return None + + def _b(self, eSolution, srcList): + return self._bPrimary(eSolution, srcList) + self._bSecondary(eSolution, srcList) + + def _bDeriv_u(self, src, v, adjoint=False): + # Primary does not depend on u + return self._bSecondaryDeriv_u(src, v, adjoint) + + def _bDeriv_m(self, src, v, adjoint=False): + # Assuming the primary does not depend on the model + return self._bSecondaryDeriv_m(src, v, adjoint) + + +class FieldsFDEM_b(FieldsFDEM): + knownFields = {'bSolution':'F'} + aliasFields = { + 'b' : ['bSolution','F','_b'], + 'bPrimary' : ['bSolution','F','_bPrimary'], + 'bSecondary' : ['bSolution','F','_bSecondary'], + 'e' : ['bSolution','E','_e'], + 'ePrimary' : ['bSolution','E','_ePrimary'], + 'eSecondary' : ['bSolution','E','_eSecondary'], + } + + def __init__(self,mesh,survey,**kwargs): + FieldsFDEM.__init__(self,mesh,survey,**kwargs) + + def startup(self): + self.prob = self.survey.prob + self._edgeCurl = self.survey.prob.mesh.edgeCurl + self._MeSigmaI = self.survey.prob.MeSigmaI + self._MfMui = self.survey.prob.MfMui + self._MeSigmaIDeriv = self.survey.prob.MeSigmaIDeriv + + def _bPrimary(self, bSolution, srcList): + bPrimary = np.zeros_like(bSolution) + for i, src in enumerate(srcList): + bp = src.bPrimary(self.prob) + if bp is not None: + bPrimary[:,i] = bp + return bPrimary + + def _bSecondary(self, bSolution, srcList): + return bSolution + + def _b(self, bSolution, srcList): + return self._bPrimary(bSolution, srcList) + self._bSecondary(bSolution, srcList) + + def _bDeriv_u(self, src, v, adjoint=False): + return None + + def _bDeriv_m(self, src, v, adjoint=False): + # assuming primary does not depend on the model + return None + + def _ePrimary(self, bSolution, srcList): + ePrimary = np.zeros([self._edgeCurl.shape[1],bSolution.shape[1]],dtype = complex) + for i,src in enumerate(srcList): + ep = src.ePrimary(self.prob) + if ep is not None: + ePrimary[:,i] = ep + return ePrimary + + def _eSecondary(self, bSolution, srcList): + e = self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * bSolution)) + for i,src in enumerate(srcList): + _,S_e = src.eval(self.prob) + if S_e is not None: + e[:,i] += -self._MeSigmaI*S_e + return e + + def _eSecondaryDeriv_u(self, src, v, adjoint=False): + if not adjoint: + return self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * v) ) + else: + return self._MfMui.T * (self._edgeCurl * (self._MeSigmaI.T * v)) + + def _eSecondaryDeriv_m(self, src, v, adjoint=False): + bSolution = self[[src],'bSolution'] + _,S_e = src.eval(self.prob) + + w = self._edgeCurl.T * (self._MfMui * bSolution) + if S_e is not None: + w += -Utils.mkvc(S_e,2) + + if not adjoint: + de_dm = self._MeSigmaIDeriv(w) * v + elif adjoint: + de_dm = self._MeSigmaIDeriv(w).T * v + + _, S_eDeriv = src.evalDeriv(self.prob, adjoint) + Se_Deriv = S_eDeriv(v) + + if Se_Deriv is not None: + de_dm += -self._MeSigmaI * Se_Deriv + + return de_dm + + def _e(self, bSolution, srcList): + return self._ePrimary(bSolution, srcList) + self._eSecondary(bSolution, srcList) + + def _eDeriv_u(self, src, v, adjoint=False): + return self._eSecondaryDeriv_u(src, v, adjoint) + + def _eDeriv_m(self, src, v, adjoint=False): + # assuming primary doesn't depend on model + return self._eSecondaryDeriv_m(src, v, adjoint) + + +class FieldsFDEM_j(FieldsFDEM): + knownFields = {'jSolution':'F'} + aliasFields = { + 'j' : ['jSolution','F','_j'], + 'jPrimary' : ['jSolution','F','_jPrimary'], + 'jSecondary' : ['jSolution','F','_jSecondary'], + 'h' : ['jSolution','E','_h'], + 'hPrimary' : ['jSolution','E','_hPrimary'], + 'hSecondary' : ['jSolution','E','_hSecondary'], + } + + def __init__(self,mesh,survey,**kwargs): + FieldsFDEM.__init__(self,mesh,survey,**kwargs) + + def startup(self): + self.prob = self.survey.prob + self._edgeCurl = self.survey.prob.mesh.edgeCurl + self._MeMuI = self.survey.prob.MeMuI + self._MfRho = self.survey.prob.MfRho + self._MfRhoDeriv = self.survey.prob.MfRhoDeriv + + def _jPrimary(self, jSolution, srcList): + jPrimary = np.zeros_like(jSolution,dtype = complex) + for i, src in enumerate(srcList): + jp = src.jPrimary(self.prob) + if jp is not None: + jPrimary[:,i] += jp + return jPrimary + + def _jSecondary(self, jSolution, srcList): + return jSolution + + def _j(self, jSolution, srcList): + return self._jPrimary(jSolution, srcList) + self._jSecondary(jSolution, srcList) + + def _jDeriv_u(self, src, v, adjoint=False): + return None + + def _jDeriv_m(self, src, v, adjoint=False): + # assuming primary does not depend on the model + return None + + def _hPrimary(self, jSolution, srcList): + hPrimary = np.zeros([self._edgeCurl.shape[1],jSolution.shape[1]],dtype = complex) + for i, src in enumerate(srcList): + hp = src.hPrimary(self.prob) + if hp is not None: + hPrimary[:,i] = hp + return hPrimary + + def _hSecondary(self, jSolution, srcList): + MeMuI = self._MeMuI + C = self._edgeCurl + MfRho = self._MfRho + h = MeMuI * (C.T * (MfRho * jSolution) ) + for i, src in enumerate(srcList): + h[:,i] *= -1./(1j*omega(src.freq)) + S_m,_ = src.eval(self.prob) + if S_m is not None: + h[:,i] += 1./(1j*omega(src.freq)) * MeMuI * S_m + return h + + def _hSecondaryDeriv_u(self, src, v, adjoint=False): + MeMuI = self._MeMuI + C = self._edgeCurl + MfRho = self._MfRho + if not adjoint: + return -1./(1j*omega(src.freq)) * MeMuI * (C.T * (MfRho * v) ) + elif adjoint: + return -1./(1j*omega(src.freq)) * MfRho.T * (C * ( MeMuI.T * v)) + + def _hSecondaryDeriv_m(self, src, v, adjoint=False): + jSolution = self[[src],'jSolution'] + MeMuI = self._MeMuI + C = self._edgeCurl + MfRho = self._MfRho + MfRhoDeriv = self._MfRhoDeriv + + if not adjoint: + hDeriv_m = -1./(1j*omega(src.freq)) * MeMuI * (C.T * (MfRhoDeriv(jSolution)*v ) ) + elif adjoint: + hDeriv_m = -1./(1j*omega(src.freq)) * MfRhoDeriv(jSolution).T * ( C * (MeMuI.T * v ) ) + + S_mDeriv,_ = src.evalDeriv(self.prob, adjoint) + + if not adjoint: + S_mDeriv = S_mDeriv(v) + if S_mDeriv is not None: + hDeriv_m += 1./(1j*omega(src.freq)) * MeMuI * S_mDeriv + elif adjoint: + S_mDeriv = S_mDeriv(MeMuI.T * v) + if S_mDeriv is not None: + hDeriv_m += 1./(1j*omega(src.freq)) * S_mDeriv + return hDeriv_m + + + def _h(self, jSolution, srcList): + return self._hPrimary(jSolution, srcList) + self._hSecondary(jSolution, srcList) + + def _hDeriv_u(self, src, v, adjoint=False): + return self._hSecondaryDeriv_u(src, v, adjoint) + + def _hDeriv_m(self, src, v, adjoint=False): + # assuming the primary doesn't depend on the model + return self._hSecondaryDeriv_m(src, v, adjoint) + + +class FieldsFDEM_h(FieldsFDEM): + knownFields = {'hSolution':'E'} + aliasFields = { + 'h' : ['hSolution','E','_h'], + 'hPrimary' : ['hSolution','E','_hPrimary'], + 'hSecondary' : ['hSolution','E','_hSecondary'], + 'j' : ['hSolution','F','_j'], + 'jPrimary' : ['hSolution','F','_jPrimary'], + 'jSecondary' : ['hSolution','F','_jSecondary'] + } + + def __init__(self,mesh,survey,**kwargs): + FieldsFDEM.__init__(self,mesh,survey,**kwargs) + + def startup(self): + self.prob = self.survey.prob + self._edgeCurl = self.survey.prob.mesh.edgeCurl + self._MeMuI = self.survey.prob.MeMuI + self._MfRho = self.survey.prob.MfRho + + def _hPrimary(self, hSolution, srcList): + hPrimary = np.zeros_like(hSolution,dtype = complex) + for i, src in enumerate(srcList): + hp = src.hPrimary(self.prob) + if hp is not None: + hPrimary[:,i] += hp + return hPrimary + + def _hSecondary(self, hSolution, srcList): + return hSolution + + def _h(self, hSolution, srcList): + return self._hPrimary(hSolution, srcList) + self._hSecondary(hSolution, srcList) + + def _hDeriv_u(self, src, v, adjoint=False): + return None + + def _hDeriv_m(self, src, v, adjoint=False): + # assuming primary does not depend on the model + return None + + def _jPrimary(self, hSolution, srcList): + jPrimary = np.zeros([self._edgeCurl.shape[0], hSolution.shape[1]]) + for i, src in enumerate(srcList): + jp = src.jPrimary(self.prob) + if jp is not None: + jPrimary[:,i] = jp + return jPrimary + + def _jSecondary(self, hSolution, srcList): + j = self._edgeCurl*hSolution + for i, src in enumerate(srcList): + _,S_e = src.eval(self.prob) + if S_e is not None: + j[:,i] += -S_e + return j + + def _jSecondaryDeriv_u(self, src, v, adjoint=False): + if not adjoint: + return self._edgeCurl*v + elif adjoint: + return self._edgeCurl.T*v + + def _jSecondaryDeriv_m(self, src, v, adjoint=False): + _,S_eDeriv = src.evalDeriv(self.prob, adjoint) + S_eDeriv = S_eDeriv(v) + if S_eDeriv is not None: + return -S_eDeriv + return None + + def _j(self, hSolution, srcList): + return self._jPrimary(hSolution, srcList) + self._jSecondary(hSolution, srcList) + + def _jDeriv_u(self, src, v, adjoint=False): + return self._jSecondaryDeriv_u(src,v,adjoint) + + def _jDeriv_m(self, src, v, adjoint=False): + # assuming the primary does not depend on the model + return self._jSecondaryDeriv_m(src,v,adjoint) \ No newline at end of file diff --git a/simpegEM/FDEM/SurveyFDEM.py b/simpegEM/FDEM/SurveyFDEM.py index e9000b7e..b47266cf 100644 --- a/simpegEM/FDEM/SurveyFDEM.py +++ b/simpegEM/FDEM/SurveyFDEM.py @@ -1,4 +1,11 @@ from SimPEG import Survey, Problem, Utils, np, sp +from simpegEM.Utils import SrcUtils +from simpegEM.Utils.EMUtils import omega, e_from_j, j_from_e, b_from_h, h_from_b +from scipy.constants import mu_0 + +#################################################### +# Receivers +#################################################### class RxFDEM(Survey.BaseRx): @@ -51,15 +58,15 @@ class RxFDEM(Survey.BaseRx): """Component projection (real/imag)""" return self.knownRxTypes[self.rxType][2] - def projectFields(self, tx, mesh, u): + def projectFields(self, src, mesh, u): P = self.getP(mesh) - u_part_complex = u[tx, self.projField] + u_part_complex = u[src, self.projField] # get the real or imag component real_or_imag = self.projComp u_part = getattr(u_part_complex, real_or_imag) return P*u_part - def projectFieldsDeriv(self, tx, mesh, u, v, adjoint=False): + def projectFieldsDeriv(self, src, mesh, u, v, adjoint=False): P = self.getP(mesh) if not adjoint: @@ -80,45 +87,330 @@ class RxFDEM(Survey.BaseRx): return Pv -class TxFDEM(Survey.BaseTx): - - freq = None #: Frequency (float) +#################################################### +# Sources +#################################################### +class SrcFDEM(Survey.BaseSrc): + freq = None rxPair = RxFDEM - knownTxTypes = ['VMD', 'VMD_B', 'CircularLoop'] + def eval(self, prob): + S_m = self.S_m(prob) + S_e = self.S_e(prob) + return S_m, S_e - radius = None + def evalDeriv(self, prob, v, adjoint=False): + return lambda v: self.S_mDeriv(prob,v,adjoint), lambda v: self.S_eDeriv(prob,v,adjoint) - def __init__(self, loc, txType, freq, rxList): + def bPrimary(self, prob): + return None + + def hPrimary(self, prob): + return None + + def ePrimary(self, prob): + return None + + def jPrimary(self, prob): + return None + + def S_m(self, prob): + return None + + def S_e(self, prob): + return None + + def S_mDeriv(self, prob, v, adjoint = False): + return None + + def S_eDeriv(self, prob, v, adjoint = False): + return None + + +class SrcFDEM_RawVec_e(SrcFDEM): + """ + RawVec electric source. It is defined by the user provided vector S_e + + :param numpy.array S_e: electric source term + :param float freq: frequency + :param rxList: receiver list + """ + + def __init__(self, rxList, freq, S_e): + self._S_e = np.array(S_e,dtype=complex) self.freq = float(freq) - Survey.BaseTx.__init__(self, loc, txType, rxList) + SrcFDEM.__init__(self, rxList) + + def S_e(self, prob): + return self._S_e +class SrcFDEM_RawVec_m(SrcFDEM): + """ + RawVec magnetic source. It is defined by the user provided vector S_m -class FieldsFDEM(Problem.Fields): - """Fancy Field Storage for a FDEM survey.""" - knownFields = {'b': 'F', 'e': 'E', 'j': 'F', 'h': 'E'} - dtype = complex + :param numpy.array S_m: magnetic source term + :param float freq: frequency + :param rxList: receiver list + """ + def __init__(self, rxList, freq, S_m): + self._S_m = np.array(S_m,dtype=complex) + self.freq = float(freq) + SrcFDEM.__init__(self, rxList) + + def S_m(self, prob): + return self._S_m + + +class SrcFDEM_RawVec(SrcFDEM): + """ + RawVec source. It is defined by the user provided vectors S_m, S_e + + :param numpy.array S_m: magnetic source term + :param numpy.array S_e: electric source term + :param float freq: frequency + :param rxList: receiver list + """ + def __init__(self, rxList, freq, S_m, S_e): + self._S_m = np.array(S_m,dtype=complex) + self._S_e = np.array(S_e,dtype=complex) + self.freq = float(freq) + SrcFDEM.__init__(self, rxList) + + def S_m(self, prob): + return self._S_m + + def S_e(self, prob): + return self._S_e + + +class SrcFDEM_MagDipole(SrcFDEM): + + #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that + def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0): + self.freq = float(freq) + self.loc = loc + self.orientation = orientation + self.moment = moment + self.mu = mu + SrcFDEM.__init__(self, rxList) + + def bPrimary(self, prob): + eqLocs = prob._eqLocs + + if eqLocs is 'FE': + gridX = prob.mesh.gridEx + gridY = prob.mesh.gridEy + gridZ = prob.mesh.gridEz + C = prob.mesh.edgeCurl + + elif eqLocs is 'EF': + gridX = prob.mesh.gridFx + gridY = prob.mesh.gridFy + gridZ = prob.mesh.gridFz + C = prob.mesh.edgeCurl.T + + + if prob.mesh._meshType is 'CYL': + if not prob.mesh.isSymmetric: + # TODO ? + raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') + a = SrcUtils.MagneticDipoleVectorPotential(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) + + else: + srcfct = SrcUtils.MagneticDipoleVectorPotential + ax = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment) + ay = srcfct(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) + az = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment) + a = np.concatenate((ax, ay, az)) + + return C*a + + def hPrimary(self, prob): + b = self.bPrimary(prob) + return h_from_b(prob,b) + + def S_m(self, prob): + b_p = self.bPrimary(prob) + return -1j*omega(self.freq)*b_p + + def S_e(self, prob): + if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): + return None + else: + eqLocs = prob._eqLocs + + if eqLocs is 'FE': + mui_s = prob.curModel.mui - 1./self.mu + MMui_s = prob.mesh.getFaceInnerProduct(mui_s) + C = prob.mesh.edgeCurl + elif eqLocs is 'EF': + mu_s = prob.curModel.mu - self.mu + MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) + C = prob.mesh.edgeCurl.T + + return -C.T * (MMui_s * self.bPrimary(prob)) + + +class SrcFDEM_MagDipole_Bfield(SrcFDEM): + + #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that + #TODO: neither does moment + def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0): + self.freq = float(freq) + self.loc = loc + self.orientation = orientation + self.moment = moment + self.mu = mu + SrcFDEM.__init__(self, rxList) + + def bPrimary(self, prob): + eqLocs = prob._eqLocs + + if eqLocs is 'FE': + gridX = prob.mesh.gridFx + gridY = prob.mesh.gridFy + gridZ = prob.mesh.gridFz + C = prob.mesh.edgeCurl + + elif eqLocs is 'EF': + gridX = prob.mesh.gridEx + gridY = prob.mesh.gridEy + gridZ = prob.mesh.gridEz + C = prob.mesh.edgeCurl.T + + srcfct = SrcUtils.MagneticDipoleFields + if prob.mesh._meshType is 'CYL': + if not prob.mesh.isSymmetric: + # TODO ? + raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') + bx = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment) + bz = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment) + b = np.concatenate((bx,bz)) + else: + bx = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment) + by = srcfct(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) + bz = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment) + b = np.concatenate((bx,by,bz)) + + return b + + def hPrimary(self, prob): + b = self.bPrimary(prob) + return h_from_b(prob, b) + + def S_m(self, prob): + b = self.bPrimary(prob) + return -1j*omega(self.freq)*b + + def S_e(self, prob): + if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): + return None + else: + eqLocs = prob._eqLocs + + if eqLocs is 'FE': + mui_s = prob.curModel.mui - 1./self.mu + MMui_s = prob.mesh.getFaceInnerProduct(mui_s) + C = prob.mesh.edgeCurl + elif eqLocs is 'EF': + mu_s = prob.curModel.mu - self.mu + MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) + C = prob.mesh.edgeCurl.T + + return -C.T * (MMui_s * self.bPrimary(prob)) + + +class SrcFDEM_CircularLoop(SrcFDEM): + + #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that + def __init__(self, rxList, freq, loc, orientation='Z', radius = 1., mu=mu_0): + self.freq = float(freq) + self.orientation = orientation + self.radius = radius + self.mu = mu + self.loc = loc + SrcFDEM.__init__(self, rxList) + + def bPrimary(self, prob): + eqLocs = prob._eqLocs + + if eqLocs is 'FE': + gridX = prob.mesh.gridEx + gridY = prob.mesh.gridEy + gridZ = prob.mesh.gridEz + C = prob.mesh.edgeCurl + + elif eqLocs is 'EF': + gridX = prob.mesh.gridFx + gridY = prob.mesh.gridFy + gridZ = prob.mesh.gridFz + C = prob.mesh.edgeCurl.T + + if prob.mesh._meshType is 'CYL': + if not prob.mesh.isSymmetric: + # TODO ? + raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') + a = SrcUtils.MagneticDipoleVectorPotential(self.loc, gridY, 'y', moment=self.radius, mu=self.mu) + + else: + srcfct = SrcUtils.MagneticDipoleVectorPotential + ax = srcfct(self.loc, gridX, 'x', self.radius, mu=self.mu) + ay = srcfct(self.loc, gridY, 'y', self.radius, mu=self.mu) + az = srcfct(self.loc, gridZ, 'z', self.radius, mu=self.mu) + a = np.concatenate((ax, ay, az)) + + return C*a + + def hPrimary(self, prob): + b = self.bPrimary(prob) + return 1./self.mu*b + + def S_m(self, prob): + b = self.bPrimary(prob) + return -1j*omega(self.freq)*b + + def S_e(self, prob): + if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): + return None + else: + eqLocs = prob._eqLocs + + if eqLocs is 'FE': + mui_s = prob.curModel.mui - 1./self.mu + MMui_s = prob.mesh.getFaceInnerProduct(mui_s) + C = prob.mesh.edgeCurl + elif eqLocs is 'EF': + mu_s = prob.curModel.mu - self.mu + MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) + C = prob.mesh.edgeCurl.T + + return -C.T * (MMui_s * self.bPrimary(prob)) + + +#################################################### +# Survey +#################################################### class SurveyFDEM(Survey.BaseSurvey): """ docstring for SurveyFDEM """ - txPair = TxFDEM + srcPair = SrcFDEM - def __init__(self, txList, **kwargs): + def __init__(self, srcList, **kwargs): # Sort these by frequency - self.txList = txList + self.srcList = srcList Survey.BaseSurvey.__init__(self, **kwargs) _freqDict = {} - for tx in txList: - if tx.freq not in _freqDict: - _freqDict[tx.freq] = [] - _freqDict[tx.freq] += [tx] + for src in srcList: + if src.freq not in _freqDict: + _freqDict[src.freq] = [] + _freqDict[src.freq] += [src] self._freqDict = _freqDict self._freqs = sorted([f for f in self._freqDict]) @@ -134,24 +426,24 @@ class SurveyFDEM(Survey.BaseSurvey): return len(self._freqDict) @property - def nTxByFreq(self): - if getattr(self, '_nTxByFreq', None) is None: - self._nTxByFreq = {} + def nSrcByFreq(self): + if getattr(self, '_nSrcByFreq', None) is None: + self._nSrcByFreq = {} for freq in self.freqs: - self._nTxByFreq[freq] = len(self.getTransmitters(freq)) - return self._nTxByFreq + self._nSrcByFreq[freq] = len(self.getSrcByFreq(freq)) + return self._nSrcByFreq - def getTransmitters(self, freq): - """Returns the transmitters associated with a specific frequency.""" + def getSrcByFreq(self, freq): + """Returns the sources associated with a specific frequency.""" assert freq in self._freqDict, "The requested frequency is not in this survey." return self._freqDict[freq] def projectFields(self, u): data = Survey.Data(self) - for tx in self.txList: - for rx in tx.rxList: - data[tx, rx] = rx.projectFields(tx, self.mesh, u) + for src in self.srcList: + for rx in src.rxList: + data[src, rx] = rx.projectFields(src, self.mesh, u) return data def projectFieldsDeriv(self, u): - raise Exception('Use Transmitters to project fields deriv.') + raise Exception('Use Sources to project fields deriv.') diff --git a/simpegEM/FDEM/__init__.py b/simpegEM/FDEM/__init__.py index 562b9218..110b4d1e 100644 --- a/simpegEM/FDEM/__init__.py +++ b/simpegEM/FDEM/__init__.py @@ -1,2 +1,3 @@ from SurveyFDEM import * -from FDEM import ProblemFDEM_e, ProblemFDEM_b, ProblemFDEM_j, ProblemFDEM_h +from FDEM import BaseFDEMProblem, ProblemFDEM_e, ProblemFDEM_b, ProblemFDEM_j, ProblemFDEM_h +from FieldsFDEM import * \ No newline at end of file diff --git a/simpegEM/Sources/CircularLoop.py b/simpegEM/Sources/CircularLoop.py deleted file mode 100644 index 6f5f2233..00000000 --- a/simpegEM/Sources/CircularLoop.py +++ /dev/null @@ -1,101 +0,0 @@ -from SimPEG import * -from scipy.special import ellipk, ellipe - -def MagneticLoopVectorPotential(txLoc, obsLoc, component, radius): - """ - Calculate the vector potential of horizontal circular loop - at given locations - - :param numpy.ndarray txLoc: Location of the transmitter(s) (x, y, z) - :param numpy.ndarray,SimPEG.Mesh obsLoc: Where the potentials will be calculated (x, y, z) or a SimPEG Mesh - :param str,list component: The component to calculate - 'x', 'y', or 'z' if an array, or grid type if mesh, can be a list - :param numpy.ndarray I: Input current of the loop - :param numpy.ndarray radius: radius of the loop - :rtype: numpy.ndarray - :return: The vector potential each dipole at each observation location - """ - - if type(component) in [list, tuple]: - out = range(len(component)) - for i, comp in enumerate(component): - out[i] = MagneticLoopVectorPotential(txLoc, obsLoc, comp, radius) - return np.concatenate(out) - - if isinstance(obsLoc, Mesh.BaseMesh): - mesh = obsLoc - assert component in ['Ex','Ey','Ez','Fx','Fy','Fz'], "Components must be in: ['Ex','Ey','Ez','Fx','Fy','Fz']" - return MagneticLoopVectorPotential(txLoc, getattr(mesh,'grid'+component), component[1], radius) - - txLoc = np.atleast_2d(txLoc) - obsLoc = np.atleast_2d(obsLoc) - - n = obsLoc.shape[0] - nTx = txLoc.shape[0] - - if component=='z': - A = np.zeros((n, nTx)) - if nTx ==1: - return A.flatten() - return A - - else: - - A = np.zeros((n, nTx)) - for i in range (nTx): - x = obsLoc[:, 0] - txLoc[i, 0] - y = obsLoc[:, 1] - txLoc[i, 1] - z = obsLoc[:, 2] - txLoc[i, 2] - r = np.sqrt(x**2 + y**2) - m = (4 * radius * r) / ((radius + r)**2 + z**2) - m[m > 1.] = 1. - # m might be slightly larger than 1 due to rounding errors - # but ellipke requires 0 <= m <= 1 - K = ellipk(m) - E = ellipe(m) - ind = (r > 0) & (m < 1) - # % 1/r singular at r = 0 and K(m) singular at m = 1 - Aphi = np.zeros(n) - # % Common factor is (mu * I) / pi with I = 1 and mu = 4e-7 * pi. - Aphi[ind] = 4e-7 / np.sqrt(m[ind]) * np.sqrt(radius / r[ind]) *((1. - m[ind] / 2.) * K[ind] - E[ind]) - if component == 'x': - A[ind, i] = Aphi[ind] * (-y[ind] / r[ind] ) - elif component == 'y': - A[ind, i] = Aphi[ind] * ( x[ind] / r[ind] ) - else: - raise ValueError('Invalid component') - - if nTx == 1: - return A.flatten() - return A - -if __name__ == '__main__': - from SimPEG import Mesh - import matplotlib.pyplot as plt - cs = 20 - ncx, ncy, ncz = 41, 41, 40 - hx = np.ones(ncx)*cs - hy = np.ones(ncy)*cs - hz = np.ones(ncz)*cs - mesh = Mesh.TensorMesh([hx, hy, hz], 'CCC') - txLoc = np.r_[0., 0., 0.] - Ax = MagneticLoopVectorPotential(txLoc, mesh.gridEx, 'x', 200) - Ay = MagneticLoopVectorPotential(txLoc, mesh.gridEy, 'y', 200) - Az = MagneticLoopVectorPotential(txLoc, mesh.gridEz, 'z', 200) - A = np.r_[Ax, Ay, Az] - B0 = mesh.edgeCurl*A - J0 = mesh.edgeCurl.T*B0 - - # mesh.plotImage(A, vType = 'Ex') - # mesh.plotImage(A, vType = 'Ey') - - mesh.plotImage(B0, vType = 'Fx') - mesh.plotImage(B0, vType = 'Fy') - mesh.plotImage(B0, vType = 'Fz') - - # # mesh.plotImage(J0, vType = 'Ex') - # mesh.plotImage(J0, vType = 'Ey') - # mesh.plotImage(J0, vType = 'Ez') - - plt.show() - - diff --git a/simpegEM/Sources/__init__.py b/simpegEM/Sources/__init__.py deleted file mode 100644 index 1f04379e..00000000 --- a/simpegEM/Sources/__init__.py +++ /dev/null @@ -1,3 +0,0 @@ -from magneticDipole import MagneticDipoleVectorPotential -from CircularLoop import MagneticLoopVectorPotential -from magneticDipole import MagneticDipoleFields diff --git a/simpegEM/Sources/magneticDipole.py b/simpegEM/Sources/magneticDipole.py deleted file mode 100644 index 3a102e74..00000000 --- a/simpegEM/Sources/magneticDipole.py +++ /dev/null @@ -1,100 +0,0 @@ -import numpy as np -from scipy.constants import mu_0, pi -from SimPEG import Mesh - -def MagneticDipoleVectorPotential(txLoc, obsLoc, component, dipoleMoment=(0., 0., 1.)): - """ - Calculate the vector potential of a set of magnetic dipoles - at given locations 'ref. ' - - :param numpy.ndarray txLoc: Location of the transmitter(s) (x, y, z) - :param numpy.ndarray,SimPEG.Mesh obsLoc: Where the potentials will be calculated (x, y, z) or a SimPEG Mesh - :param str,list component: The component to calculate - 'x', 'y', or 'z' if an array, or grid type if mesh, can be a list - :param numpy.ndarray dipoleMoment: The vector dipole moment - :rtype: numpy.ndarray - :return: The vector potential each dipole at each observation location - """ - - if type(component) in [list, tuple]: - out = range(len(component)) - for i, comp in enumerate(component): - out[i] = MagneticDipoleVectorPotential(txLoc, obsLoc, comp, dipoleMoment=dipoleMoment) - return np.concatenate(out) - - if isinstance(obsLoc, Mesh.BaseMesh): - mesh = obsLoc - assert component in ['Ex','Ey','Ez','Fx','Fy','Fz'], "Components must be in: ['Ex','Ey','Ez','Fx','Fy','Fz']" - return MagneticDipoleVectorPotential(txLoc, getattr(mesh,'grid'+component), component[1], dipoleMoment=dipoleMoment) - - if component == 'x': - dimInd = 0 - elif component == 'y': - dimInd = 1 - elif component == 'z': - dimInd = 2 - else: - raise ValueError('Invalid component') - - txLoc = np.atleast_2d(txLoc) - obsLoc = np.atleast_2d(obsLoc) - dipoleMoment = np.atleast_2d(dipoleMoment) - - nEdges = obsLoc.shape[0] - nTx = txLoc.shape[0] - - m = np.array(dipoleMoment).repeat(nEdges, axis=0) - A = np.empty((nEdges, nTx)) - for i in range(nTx): - dR = obsLoc - txLoc[i, np.newaxis].repeat(nEdges, axis=0) - mCr = np.cross(m, dR) - r = np.sqrt((dR**2).sum(axis=1)) - A[:, i] = +(mu_0/(4*pi)) * mCr[:,dimInd]/(r**3) - if nTx == 1: - return A.flatten() - return A - -def MagneticDipoleFields(txLoc, obsLoc, component, dipoleMoment=1.): - """ - Calculate the vector potential of a set of magnetic dipoles - at given locations 'ref. ' - - :param numpy.ndarray txLoc: Location of the transmitter(s) (x, y, z) - :param numpy.ndarray obsLoc: Where the potentials will be calculated (x, y, z) - :param str component: The component to calculate - 'x', 'y', or 'z' - :param numpy.ndarray dipoleMoment: The vector dipole moment (vertical) - :rtype: numpy.ndarray - :return: The vector potential each dipole at each observation location - """ - - if component=='x': - dimInd = 0 - elif component=='y': - dimInd = 1 - elif component=='z': - dimInd = 2 - else: - raise ValueError('Invalid component') - - txLoc = np.atleast_2d(txLoc) - obsLoc = np.atleast_2d(obsLoc) - dipoleMoment = np.atleast_2d(dipoleMoment) - - nFaces = obsLoc.shape[0] - nTx = txLoc.shape[0] - - m = np.array(dipoleMoment).repeat(nFaces, axis=0) - B = np.empty((nFaces, nTx)) - for i in range(nTx): - dR = obsLoc - txLoc[i, np.newaxis].repeat(nFaces, axis=0) - r = np.sqrt((dR**2).sum(axis=1)) - if dimInd == 0: - B[:, i] = +(mu_0/(4*pi)) /(r**3) * (3*dR[:,2]*dR[:,0]/r**2) - elif dimInd == 1: - B[:, i] = +(mu_0/(4*pi)) /(r**3) * (3*dR[:,2]*dR[:,1]/r**2) - elif dimInd == 2: - B[:, i] = +(mu_0/(4*pi)) /(r**3) * (3*dR[:,2]**2/r**2-1) - else: - raise Exception("Not Implemented") - if nTx == 1: - return B.flatten() - return B diff --git a/simpegEM/TDEM/BaseTDEM.py b/simpegEM/TDEM/BaseTDEM.py index 7626d944..8f4be4cd 100644 --- a/simpegEM/TDEM/BaseTDEM.py +++ b/simpegEM/TDEM/BaseTDEM.py @@ -1,6 +1,6 @@ from SimPEG import Solver, Problem from SimPEG.Problem import BaseTimeProblem -from simpegEM import Sources +from simpegEM.Utils import SrcUtils from scipy.constants import mu_0 from SimPEG.Utils import sdiag, mkvc from SimPEG import Utils, Mesh @@ -13,21 +13,21 @@ class FieldsTDEM(Problem.TimeFields): knownFields = {'b': 'F', 'e': 'E'} def tovec(self): - nTx, nF, nE = self.survey.nTx, self.mesh.nF, self.mesh.nE - u = np.empty(0 if nTx == 1 else (0, nTx)) + nSrc, nF, nE = self.survey.nSrc, self.mesh.nF, self.mesh.nE + u = np.empty((0,nSrc)) #((0,1) if nSrc == 1 else (0, nSrc)) for i in range(self.survey.prob.nT): if 'b' in self: b = self[:,'b',i+1] else: - b = np.zeros(nF if nTx == 1 else (nF, nTx)) + b = np.zeros((nF,nSrc)) # if nSrc == 1 else (nF, nSrc)) if 'e' in self: e = self[:,'e',i+1] else: - e = np.zeros(nE if nTx == 1 else (nE, nTx)) + e = np.zeros((nE,nSrc)) # if nSrc == 1 else (nE, nSrc)) u = np.concatenate((u, b, e)) - return Utils.mkvc(u) + return Utils.mkvc(u,nSrc) class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem): @@ -42,9 +42,9 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem): self.curModel = m # Create a fields storage object F = self._FieldsForward_pair(self.mesh, self.survey) - for tx in self.survey.txList: + for src in self.survey.srcList: # Set the initial conditions - F[tx,:,0] = tx.getInitialFields(self.mesh) + 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) return F diff --git a/simpegEM/TDEM/SurveyTDEM.py b/simpegEM/TDEM/SurveyTDEM.py index 86bc2b5a..8547847b 100644 --- a/simpegEM/TDEM/SurveyTDEM.py +++ b/simpegEM/TDEM/SurveyTDEM.py @@ -1,6 +1,6 @@ from SimPEG import Utils, Survey, np from SimPEG.Survey import BaseSurvey -from simpegEM import Sources +from simpegEM.Utils import SrcUtils from BaseTDEM import FieldsTDEM @@ -51,76 +51,90 @@ class RxTDEM(Survey.BaseTimeRx): else: return timeMesh.getInterpolationMat(self.times, self.projTLoc) - def projectFields(self, tx, mesh, timeMesh, u): + def projectFields(self, src, mesh, timeMesh, u): P = self.getP(mesh, timeMesh) - u_part = Utils.mkvc(u[tx, self.projField, :]) + u_part = Utils.mkvc(u[src, self.projField, :]) return P*u_part - def projectFieldsDeriv(self, tx, mesh, timeMesh, u, v, adjoint=False): + def projectFieldsDeriv(self, src, mesh, timeMesh, u, v, adjoint=False): P = self.getP(mesh, timeMesh) if not adjoint: - return P * Utils.mkvc(v[tx, self.projField, :]) + return P * Utils.mkvc(v[src, self.projField, :]) elif adjoint: - return P.T * v[tx, self] + return P.T * v[src, self] -class TxTDEM(Survey.BaseTx): +class SrcTDEM(Survey.BaseSrc): rxPair = RxTDEM radius = None - knownTxTypes = ['VMD_MVP', 'CircularLoop_MVP'] def getInitialFields(self, mesh): - F0 = getattr(self, '_getInitialFields_' + self.txType)(mesh) + F0 = getattr(self, '_getInitialFields_' + self.srcType)(mesh) return F0 - def _getInitialFields_VMD_MVP(self, mesh): + def getJs(self, mesh, time): + return None + + +class SrcTDEM_VMD_MVP(SrcTDEM): + + def __init__(self,rxList,loc): + self.loc = loc + SrcTDEM.__init__(self,rxList) + + def getInitialFields(self, mesh): """Vertical magnetic dipole, magnetic vector potential""" if mesh._meshType is 'CYL': if mesh.isSymmetric: - MVP = Sources.MagneticDipoleVectorPotential(self.loc, mesh, 'Ey') + MVP = SrcUtils.MagneticDipoleVectorPotential(self.loc, mesh, 'Ey') else: raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') elif mesh._meshType is 'TENSOR': - MVP = Sources.MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez']) + MVP = SrcUtils.MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez']) else: raise Exception('Unknown mesh for VMD') return {"b": mesh.edgeCurl*MVP} - def _getInitialFields_CircularLoop_MVP(self, mesh): + +class SrcTDEM_CircularLoop_MVP(SrcTDEM): + + def __init__(self,rxList,loc): + self.loc = loc + SrcTDEM.__init__(self,rxList) + + def getInitialFields_(self, mesh): """Circular Loop, magnetic vector potential""" if mesh._meshType is 'CYL': if mesh.isSymmetric: - MVP = Sources.MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius) + MVP = SrcUtils.MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius) else: raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') elif mesh._meshType is 'TENSOR': - MVP = Sources.MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius) + MVP = SrcUtils.MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius) else: raise Exception('Unknown mesh for CircularLoop') return {"b": mesh.edgeCurl*MVP} - def getJs(self, mesh, time): - return None class SurveyTDEM(Survey.BaseSurvey): """ docstring for SurveyTDEM """ - txPair = TxTDEM + srcPair = SrcTDEM - def __init__(self, txList, **kwargs): + def __init__(self, srcList, **kwargs): # Sort these by frequency - self.txList = txList + self.srcList = srcList Survey.BaseSurvey.__init__(self, **kwargs) def projectFields(self, u): data = Survey.Data(self) - for tx in self.txList: - for rx in tx.rxList: - data[tx, rx] = rx.projectFields(tx, self.mesh, self.prob.timeMesh, u) + for src in self.srcList: + for rx in src.rxList: + data[src, rx] = rx.projectFields(src, self.mesh, self.prob.timeMesh, u) return data def projectFieldsDeriv(self, u, v=None, adjoint=False): @@ -128,20 +142,20 @@ class SurveyTDEM(Survey.BaseSurvey): if not adjoint: data = Survey.Data(self) - for tx in self.txList: - for rx in tx.rxList: - data[tx, rx] = rx.projectFieldsDeriv(tx, self.mesh, self.prob.timeMesh, u, v) + for src in self.srcList: + for rx in src.rxList: + data[src, rx] = rx.projectFieldsDeriv(src, self.mesh, self.prob.timeMesh, u, v) return data else: f = FieldsTDEM(self.mesh, self) - for tx in self.txList: - for rx in tx.rxList: - Ptv = rx.projectFieldsDeriv(tx, self.mesh, self.prob.timeMesh, u, v, adjoint=True) + for src in self.srcList: + for rx in src.rxList: + Ptv = rx.projectFieldsDeriv(src, self.mesh, self.prob.timeMesh, u, v, adjoint=True) Ptv = Ptv.reshape((-1, self.prob.timeMesh.nN), order='F') if rx.projField not in f: # first time we are projecting - f[tx, rx.projField, :] = Ptv + f[src, rx.projField, :] = Ptv else: # there are already fields, so let's add to them! - f[tx, rx.projField, :] += Ptv + f[src, rx.projField, :] += Ptv return f diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py index ccaff98c..cd38660c 100644 --- a/simpegEM/TDEM/TDEM_b.py +++ b/simpegEM/TDEM/TDEM_b.py @@ -14,7 +14,7 @@ class FieldsTDEM_e_from_b(FieldsTDEM): self.edgeCurlT = self.survey.prob.mesh.edgeCurl.T self.MfMui = self.survey.prob.MfMui - def e_from_b(self, b, txInd, timeInd): + def e_from_b(self, b, srcInd, timeInd): # TODO: implement non-zero js return self.MeSigmaI*(self.edgeCurlT*(self.MfMui*b)) @@ -32,10 +32,10 @@ class FieldsTDEM_e_from_b_Ah(FieldsTDEM): self.edgeCurlT = self.survey.prob.mesh.edgeCurl.T self.MfMui = self.survey.prob.MfMui - def e_from_b(self, y_b, txInd, tInd): + def e_from_b(self, y_b, srcInd, tInd): y_e = self.MeSigmaI*(self.edgeCurlT*(self.MfMui*y_b)) if 'e' in self.p: - y_e = y_e - self.MeSigmaI*self.p[txInd,'e',tInd] + y_e = y_e - self.MeSigmaI*self.p[srcInd,'e',tInd] return y_e class ProblemTDEM_b(BaseTDEMProblem): @@ -73,8 +73,10 @@ class ProblemTDEM_b(BaseTDEMProblem): def getRHS(self, tInd, F): dt = self.timeSteps[tInd] - B_n = np.c_[[F[tx,'b',tInd] for tx in self.survey.txList]].T - RHS = (1.0/dt)*self.MfMui*B_n + B_n = np.c_[[F[src,'b',tInd] for src in self.survey.srcList]].T + if B_n.shape[0] is not 1: + raise NotImplementedError('getRHS not implemented for this shape of B_n') + RHS = (1.0/dt)*self.MfMui*B_n[0,:,:] #TODO: This is a hack return RHS #################################################### @@ -95,7 +97,7 @@ class ProblemTDEM_b(BaseTDEMProblem): u = self.fields(m) self.curModel = m - # Note: Fields has shape (nF/E, nTx, nT+1) + # Note: Fields has shape (nF/E, nSrc, nT+1) # However, p will only really fill (:,:,1:nT+1) # meaning the 'initial fields' are zero (:,:,0) p = FieldsTDEM(self.mesh, self.survey) @@ -106,15 +108,17 @@ class ProblemTDEM_b(BaseTDEMProblem): # fake initial 'e' fields p[:, 'e', 0] = 0.0 - dMdsig = self.mesh.getEdgeInnerProductDeriv(self.curModel.transform) - dsigdm_x_v = self.curModel.transformDeriv*vec + dMdsig = self.MeSigmaDeriv + # self.mesh.getEdgeInnerProductDeriv(self.curModel.transform) + # dsigdm_x_v = self.curModel.sigmaDeriv*vec + # dsigdm_x_v = self.curModel.transformDeriv*vec for i in range(1,self.nT+1): # TODO: G[1] may be dependent on the model # for a galvanic source (deriv of the dc problem) # - # Do multiplication for all tx in self.survey.txList - for tx in self.survey.txList: - p[tx, 'e', i] = - dMdsig(u[tx,'e',i]) * dsigdm_x_v + # Do multiplication for all src in self.survey.srcList + for src in self.survey.srcList: + p[src, 'e', i] = - dMdsig(u[src,'e',i]) * vec return p def Gtvec(self, m, vec, u=None): @@ -130,20 +134,21 @@ class ProblemTDEM_b(BaseTDEMProblem): if u is None: u = self.fields(m) self.curModel = m - dMdsig = self.mesh.getEdgeInnerProductDeriv(self.curModel.transform) - dsigdm = self.curModel.transformDeriv + # dMdsig = self.mesh.getEdgeInnerProductDeriv(self.curModel.transform) + # dsigdm = self.curModel.transformDeriv + MeSigmaDeriv = self.MeSigmaDeriv - nTx = self.survey.nTx + nSrc = self.survey.nSrc VUs = None # Here we can do internal multiplications of Gt*v and then multiply by MsigDeriv.T in one go. for i in range(1,self.nT+1): vu = None - for tx in self.survey.txList: - vutx = dMdsig(u[tx,'e',i]).T * vec[tx,'e',i] - vu = vutx if vu is None else vu + vutx + for src in self.survey.srcList: + vusrc = MeSigmaDeriv(u[src,'e',i]).T * vec[src,'e',i] + vu = vusrc if vu is None else vu + vusrc VUs = vu if VUs is None else VUs + vu - p = -dsigdm.T*VUs - return p + # p = -dsigdm.T*VUs + return -VUs def solveAh(self, m, p): """ @@ -241,8 +246,8 @@ class ProblemTDEM_b(BaseTDEMProblem): # 1 (tInd=1 uses fields 2 and 3) def AhtRHS(tInd, y): - nTx, nF = self.survey.nTx, self.mesh.nF - rhs = np.zeros(nF if nTx == 1 else (nF, nTx)) + nSrc, nF = self.survey.nSrc, self.mesh.nF + rhs = np.zeros((nF,1) if nSrc == 1 else (nF, nSrc)) if 'e' in p: rhs += self.MfMui*(self.mesh.edgeCurl*(self.MeSigmaI*p[:,'e',tInd+1])) diff --git a/simpegEM/TDEM/__init__.py b/simpegEM/TDEM/__init__.py index 15ff3f43..dd5a8bce 100644 --- a/simpegEM/TDEM/__init__.py +++ b/simpegEM/TDEM/__init__.py @@ -1,3 +1,3 @@ -from SurveyTDEM import SurveyTDEM, RxTDEM, TxTDEM +from SurveyTDEM import * #SurveyTDEM, RxTDEM, SrcTDEM from BaseTDEM import BaseTDEMProblem, FieldsTDEM from TDEM_b import ProblemTDEM_b diff --git a/simpegEM/Tests/test_FDEM.py b/simpegEM/Tests/test_FDEM.py index 3a1c4c92..24c9205f 100644 --- a/simpegEM/Tests/test_FDEM.py +++ b/simpegEM/Tests/test_FDEM.py @@ -1,9 +1,17 @@ import unittest from SimPEG import * import simpegEM as EM -import sys +import sys from scipy.constants import mu_0 +testDerivs = True +testCrossCheck = True +testAdjoint = True +testEB = True +testHJ = True + +verbose = False + TOL = 1e-4 FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order CONDUCTIVITY = 1e1 @@ -11,6 +19,9 @@ MU = mu_0 freq = 1e-1 addrandoms = True +SrcType = 'MagDipole' #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec' + + def getProblem(fdemType, comp): cs = 5. ncx, ncy, ncz = 6, 6, 6 @@ -22,22 +33,64 @@ def getProblem(fdemType, comp): mapping = Maps.ExpMap(mesh) - x = np.array([np.linspace(-30,-15,3),np.linspace(15,30,3)]) #don't sample right by the transmitter + x = np.array([np.linspace(-30,-15,3),np.linspace(15,30,3)]) #don't sample right by the source XYZ = Utils.ndgrid(x,x,np.r_[0.]) Rx0 = EM.FDEM.RxFDEM(XYZ, comp) - Tx0 = EM.FDEM.TxFDEM(np.r_[0.,0.,0.], 'VMD', freq, [Rx0]) - survey = EM.FDEM.SurveyFDEM([Tx0]) + if SrcType is 'MagDipole': + Src = EM.FDEM.SrcFDEM_MagDipole([Rx0], freq=freq, loc=np.r_[0.,0.,0.]) + elif SrcType is 'MagDipole_Bfield': + Src = EM.FDEM.SrcFDEM_MagDipole_Bfield([Rx0], freq=freq, loc=np.r_[0.,0.,0.]) + elif SrcType is 'CircularLoop': + Src2 = EM.FDEM.SrcFDEM_CircularLoop([Rx0], freq=freq, loc=np.r_[0.,0.,0.]) + if verbose: + print ' Fetching %s problem' % (fdemType) if fdemType == 'e': + if SrcType is 'RawVec': + S_m = np.zeros(mesh.nF) + S_e = np.zeros(mesh.nE) + S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1. + S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1. + Src = EM.FDEM.SrcFDEM_RawVec([Rx0], freq, S_m, S_e) + + survey = EM.FDEM.SurveyFDEM([Src]) prb = EM.FDEM.ProblemFDEM_e(mesh, mapping=mapping) + elif fdemType == 'b': + if SrcType is 'RawVec': + S_m = np.zeros(mesh.nF) + S_e = np.zeros(mesh.nE) + S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1. + S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1. + Src = EM.FDEM.SrcFDEM_RawVec([Rx0], freq, S_m, S_e) + + survey = EM.FDEM.SurveyFDEM([Src]) prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping) + elif fdemType == 'j': + if SrcType is 'RawVec': + S_m = np.zeros(mesh.nE) + S_e = np.zeros(mesh.nF) + S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1. + S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1. + Src = EM.FDEM.SrcFDEM_RawVec([Rx0], freq, S_m, S_e) + + survey = EM.FDEM.SurveyFDEM([Src]) prb = EM.FDEM.ProblemFDEM_j(mesh, mapping=mapping) + elif fdemType == 'h': + if SrcType is 'RawVec': + S_m = np.zeros(mesh.nE) + S_e = np.zeros(mesh.nF) + S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1. + S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1. + Src = EM.FDEM.SrcFDEM_RawVec([Rx0], freq, S_m, S_e) + + survey = EM.FDEM.SurveyFDEM([Src]) prb = EM.FDEM.ProblemFDEM_h(mesh, mapping=mapping) + else: raise NotImplementedError() prb.pair(survey) @@ -54,38 +107,42 @@ def adjointTest(fdemType, comp): prb = getProblem(fdemType, comp) print 'Adjoint %s formulation - %s' % (fdemType, comp) - m = np.log(np.ones(prb.mesh.nC)*CONDUCTIVITY) - mu = np.log(np.ones(prb.mesh.nC)*MU) + m = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY) + mu = np.ones(prb.mesh.nC)*MU if addrandoms is True: - m = m + np.random.randn(prb.mesh.nC)*CONDUCTIVITY*1e-1 + m = m + np.random.randn(prb.mapping.nP)*np.log(CONDUCTIVITY)*1e-1 mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-1 - prb.mu = mu survey = prb.survey - + # prb.PropMap.PropModel.mu = mu + # prb.PropMap.PropModel.mui = 1./mu u = prb.fields(m) v = np.random.rand(survey.nD) - w = np.random.rand(prb.mapping.nP) + w = np.random.rand(prb.mesh.nC) - vJw = v.dot(prb.Jvec(m, w, u=u)) - wJtv = w.dot(prb.Jtvec(m, v, u=u)) + vJw = v.dot(prb.Jvec(m, w, u)) + wJtv = w.dot(prb.Jtvec(m, v, u)) tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR]) print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol return np.abs(vJw - wJtv) < tol + def derivTest(fdemType, comp): + prb = getProblem(fdemType, comp) print '%s formulation - %s' % (fdemType, comp) - x0 = np.log(np.ones(prb.mesh.nC)*CONDUCTIVITY) + x0 = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY) mu = np.log(np.ones(prb.mesh.nC)*MU) if addrandoms is True: - x0 = x0 + np.random.randn(prb.mesh.nC)*CONDUCTIVITY*1e-1 - mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-1 + x0 = x0 + np.random.randn(prb.mapping.nP)*np.log(CONDUCTIVITY)*1e-1 + mu = mu + np.random.randn(prb.mapping.nP)*MU*1e-1 + + # prb.PropMap.PropModel.mu = mu + # prb.PropMap.PropModel.mui = 1./mu - prb.mu = mu survey = prb.survey def fun(x): return survey.dpred(x), lambda x: prb.Jvec(x0, x) @@ -103,16 +160,16 @@ def crossCheckTest(fdemType, comp): mu = np.log(np.ones(mesh.nC)*MU) if addrandoms is True: - m = m + np.random.randn(mesh.nC)*CONDUCTIVITY*1e-1 + m = m + np.random.randn(mesh.nC)*np.log(CONDUCTIVITY)*1e-1 mu = mu + np.random.randn(mesh.nC)*MU*1e-1 - prb1.mu = mu + # prb1.PropMap.PropModel.mu = mu + # prb1.PropMap.PropModel.mui = 1./mu survey1 = prb1.survey + d1 = survey1.dpred(m) - u1 = prb1.fields(m) - d1 = Utils.mkvc(survey1.projectFields(u1)) - - prb1.unpair + if verbose: + print ' Problem 1 solved' if fdemType == 'e': prb2 = getProblem('b', comp) @@ -125,11 +182,12 @@ def crossCheckTest(fdemType, comp): else: raise NotImplementedError() - prb2.mu = mu + # prb2.mu = mu survey2 = prb2.survey + d2 = survey2.dpred(m) - u2 = prb2.fields(m) - d2 = Utils.mkvc(survey2.projectFields(u2)) + if verbose: + print ' Problem 2 solved' r = d2-d1 l2r = l2norm(r) @@ -144,265 +202,277 @@ def crossCheckTest(fdemType, comp): class FDEM_DerivTests(unittest.TestCase): - def test_Jvec_exr_Eform(self): - self.assertTrue(derivTest('e', 'exr')) - def test_Jvec_exr_Bform(self): - self.assertTrue(derivTest('b', 'exr')) - def test_Jvec_eyr_Eform(self): - self.assertTrue(derivTest('e', 'eyr')) - def test_Jvec_eyr_Bform(self): - self.assertTrue(derivTest('b', 'eyr')) - def test_Jvec_ezr_Eform(self): - self.assertTrue(derivTest('e', 'ezr')) - def test_Jvec_ezr_Bform(self): - self.assertTrue(derivTest('b', 'ezr')) - def test_Jvec_exi_Eform(self): - self.assertTrue(derivTest('e', 'exi')) - def test_Jvec_exi_Bform(self): - self.assertTrue(derivTest('b', 'exi')) - def test_Jvec_eyi_Eform(self): - self.assertTrue(derivTest('e', 'eyi')) - def test_Jvec_eyi_Bform(self): - self.assertTrue(derivTest('b', 'eyi')) - def test_Jvec_ezi_Eform(self): - self.assertTrue(derivTest('e', 'ezi')) - def test_Jvec_ezi_Bform(self): - self.assertTrue(derivTest('b', 'ezi')) + if testDerivs: + if testEB: + def test_Jvec_exr_Eform(self): + self.assertTrue(derivTest('e', 'exr')) + def test_Jvec_eyr_Eform(self): + self.assertTrue(derivTest('e', 'eyr')) + def test_Jvec_ezr_Eform(self): + self.assertTrue(derivTest('e', 'ezr')) + def test_Jvec_exi_Eform(self): + self.assertTrue(derivTest('e', 'exi')) + def test_Jvec_eyi_Eform(self): + self.assertTrue(derivTest('e', 'eyi')) + def test_Jvec_ezi_Eform(self): + self.assertTrue(derivTest('e', 'ezi')) - def test_Jvec_bxr_Eform(self): - self.assertTrue(derivTest('e', 'bxr')) - def test_Jvec_bxr_Bform(self): - self.assertTrue(derivTest('b', 'bxr')) - def test_Jvec_byr_Eform(self): - self.assertTrue(derivTest('e', 'byr')) - def test_Jvec_byr_Bform(self): - self.assertTrue(derivTest('b', 'byr')) - def test_Jvec_bzr_Eform(self): - self.assertTrue(derivTest('e', 'bzr')) - def test_Jvec_bzr_Bform(self): - self.assertTrue(derivTest('b', 'bzr')) - def test_Jvec_bxi_Eform(self): - self.assertTrue(derivTest('e', 'bxi')) - def test_Jvec_bxi_Bform(self): - self.assertTrue(derivTest('b', 'bxi')) - def test_Jvec_byi_Eform(self): - self.assertTrue(derivTest('e', 'byi')) - def test_Jvec_byi_Bform(self): - self.assertTrue(derivTest('b', 'byi')) - def test_Jvec_bzi_Eform(self): - self.assertTrue(derivTest('e', 'bzi')) - def test_Jvec_bzi_Bform(self): - self.assertTrue(derivTest('b', 'bzi')) + def test_Jvec_bxr_Eform(self): + self.assertTrue(derivTest('e', 'bxr')) + def test_Jvec_byr_Eform(self): + self.assertTrue(derivTest('e', 'byr')) + def test_Jvec_bzr_Eform(self): + self.assertTrue(derivTest('e', 'bzr')) + def test_Jvec_bxi_Eform(self): + self.assertTrue(derivTest('e', 'bxi')) + def test_Jvec_byi_Eform(self): + self.assertTrue(derivTest('e', 'byi')) + def test_Jvec_bzi_Eform(self): + self.assertTrue(derivTest('e', 'bzi')) + + def test_Jvec_exr_Bform(self): + self.assertTrue(derivTest('b', 'exr')) + def test_Jvec_eyr_Bform(self): + self.assertTrue(derivTest('b', 'eyr')) + def test_Jvec_ezr_Bform(self): + self.assertTrue(derivTest('b', 'ezr')) + def test_Jvec_exi_Bform(self): + self.assertTrue(derivTest('b', 'exi')) + def test_Jvec_eyi_Bform(self): + self.assertTrue(derivTest('b', 'eyi')) + def test_Jvec_ezi_Bform(self): + self.assertTrue(derivTest('b', 'ezi')) + + def test_Jvec_bxr_Bform(self): + self.assertTrue(derivTest('b', 'bxr')) + def test_Jvec_byr_Bform(self): + self.assertTrue(derivTest('b', 'byr')) + def test_Jvec_bzr_Bform(self): + self.assertTrue(derivTest('b', 'bzr')) + def test_Jvec_bxi_Bform(self): + self.assertTrue(derivTest('b', 'bxi')) + def test_Jvec_byi_Bform(self): + self.assertTrue(derivTest('b', 'byi')) + def test_Jvec_bzi_Bform(self): + self.assertTrue(derivTest('b', 'bzi')) + + if testHJ: + def test_Jvec_jxr_Jform(self): + self.assertTrue(derivTest('j', 'jxr')) + def test_Jvec_jyr_Jform(self): + self.assertTrue(derivTest('j', 'jyr')) + def test_Jvec_jzr_Jform(self): + self.assertTrue(derivTest('j', 'jzr')) + def test_Jvec_jxi_Jform(self): + self.assertTrue(derivTest('j', 'jxi')) + def test_Jvec_jyi_Jform(self): + self.assertTrue(derivTest('j', 'jyi')) + def test_Jvec_jzi_Jform(self): + self.assertTrue(derivTest('j', 'jzi')) + + def test_Jvec_hxr_Jform(self): + self.assertTrue(derivTest('j', 'hxr')) + def test_Jvec_hyr_Jform(self): + self.assertTrue(derivTest('j', 'hyr')) + def test_Jvec_hzr_Jform(self): + self.assertTrue(derivTest('j', 'hzr')) + def test_Jvec_hxi_Jform(self): + self.assertTrue(derivTest('j', 'hxi')) + def test_Jvec_hyi_Jform(self): + self.assertTrue(derivTest('j', 'hyi')) + def test_Jvec_hzi_Jform(self): + self.assertTrue(derivTest('j', 'hzi')) + + def test_Jvec_hxr_Hform(self): + self.assertTrue(derivTest('h', 'hxr')) + def test_Jvec_hyr_Hform(self): + self.assertTrue(derivTest('h', 'hyr')) + def test_Jvec_hzr_Hform(self): + self.assertTrue(derivTest('h', 'hzr')) + def test_Jvec_hxi_Hform(self): + self.assertTrue(derivTest('h', 'hxi')) + def test_Jvec_hyi_Hform(self): + self.assertTrue(derivTest('h', 'hyi')) + def test_Jvec_hzi_Hform(self): + self.assertTrue(derivTest('h', 'hzi')) + + def test_Jvec_hxr_Hform(self): + self.assertTrue(derivTest('h', 'jxr')) + def test_Jvec_hyr_Hform(self): + self.assertTrue(derivTest('h', 'jyr')) + def test_Jvec_hzr_Hform(self): + self.assertTrue(derivTest('h', 'jzr')) + def test_Jvec_hxi_Hform(self): + self.assertTrue(derivTest('h', 'jxi')) + def test_Jvec_hyi_Hform(self): + self.assertTrue(derivTest('h', 'jyi')) + def test_Jvec_hzi_Hform(self): + self.assertTrue(derivTest('h', 'jzi')) + if testAdjoint: + if testEB: + def test_Jtvec_adjointTest_exr_Eform(self): + self.assertTrue(adjointTest('e', 'exr')) + def test_Jtvec_adjointTest_eyr_Eform(self): + self.assertTrue(adjointTest('e', 'eyr')) + def test_Jtvec_adjointTest_ezr_Eform(self): + self.assertTrue(adjointTest('e', 'ezr')) + def test_Jtvec_adjointTest_exi_Eform(self): + self.assertTrue(adjointTest('e', 'exi')) + def test_Jtvec_adjointTest_eyi_Eform(self): + self.assertTrue(adjointTest('e', 'eyi')) + def test_Jtvec_adjointTest_ezi_Eform(self): + self.assertTrue(adjointTest('e', 'ezi')) - def test_Jtvec_adjointTest_exr_Eform(self): - self.assertTrue(adjointTest('e', 'exr')) - def test_Jtvec_adjointTest_exr_Bform(self): - self.assertTrue(adjointTest('b', 'exr')) - def test_Jtvec_adjointTest_eyr_Eform(self): - self.assertTrue(adjointTest('e', 'eyr')) - def test_Jtvec_adjointTest_eyr_Bform(self): - self.assertTrue(adjointTest('b', 'eyr')) - def test_Jtvec_adjointTest_ezr_Eform(self): - self.assertTrue(adjointTest('e', 'ezr')) - def test_Jtvec_adjointTest_ezr_Bform(self): - self.assertTrue(adjointTest('b', 'ezr')) - def test_Jtvec_adjointTest_exi_Eform(self): - self.assertTrue(adjointTest('e', 'exi')) - def test_Jtvec_adjointTest_exi_Bform(self): - self.assertTrue(adjointTest('b', 'exi')) - def test_Jtvec_adjointTest_eyi_Eform(self): - self.assertTrue(adjointTest('e', 'eyi')) - def test_Jtvec_adjointTest_eyi_Bform(self): - self.assertTrue(adjointTest('b', 'eyi')) - def test_Jtvec_adjointTest_ezi_Eform(self): - self.assertTrue(adjointTest('e', 'ezi')) - def test_Jtvec_adjointTest_ezi_Bform(self): - self.assertTrue(adjointTest('b', 'ezi')) + def test_Jtvec_adjointTest_bxr_Eform(self): + self.assertTrue(adjointTest('e', 'bxr')) + def test_Jtvec_adjointTest_byr_Eform(self): + self.assertTrue(adjointTest('e', 'byr')) + def test_Jtvec_adjointTest_bzr_Eform(self): + self.assertTrue(adjointTest('e', 'bzr')) + def test_Jtvec_adjointTest_bxi_Eform(self): + self.assertTrue(adjointTest('e', 'bxi')) + def test_Jtvec_adjointTest_byi_Eform(self): + self.assertTrue(adjointTest('e', 'byi')) + def test_Jtvec_adjointTest_bzi_Eform(self): + self.assertTrue(adjointTest('e', 'bzi')) - def test_Jtvec_adjointTest_bxr_Eform(self): - self.assertTrue(adjointTest('e', 'bxr')) - def test_Jtvec_adjointTest_bxr_Bform(self): - self.assertTrue(adjointTest('b', 'bxr')) - def test_Jtvec_adjointTest_byr_Eform(self): - self.assertTrue(adjointTest('e', 'byr')) - def test_Jtvec_adjointTest_byr_Bform(self): - self.assertTrue(adjointTest('b', 'byr')) - def test_Jtvec_adjointTest_bzr_Eform(self): - self.assertTrue(adjointTest('e', 'bzr')) - def test_Jtvec_adjointTest_bzr_Bform(self): - self.assertTrue(adjointTest('b', 'bzr')) - def test_Jtvec_adjointTest_bxi_Eform(self): - self.assertTrue(adjointTest('e', 'bxi')) - def test_Jtvec_adjointTest_bxi_Bform(self): - self.assertTrue(adjointTest('b', 'bxi')) - def test_Jtvec_adjointTest_byi_Eform(self): - self.assertTrue(adjointTest('e', 'byi')) - def test_Jtvec_adjointTest_byi_Bform(self): - self.assertTrue(adjointTest('b', 'byi')) - def test_Jtvec_adjointTest_bzi_Eform(self): - self.assertTrue(adjointTest('e', 'bzi')) - def test_Jtvec_adjointTest_bzi_Bform(self): - self.assertTrue(adjointTest('b', 'bzi')) + def test_Jtvec_adjointTest_exr_Bform(self): + self.assertTrue(adjointTest('b', 'exr')) + def test_Jtvec_adjointTest_eyr_Bform(self): + self.assertTrue(adjointTest('b', 'eyr')) + def test_Jtvec_adjointTest_ezr_Bform(self): + self.assertTrue(adjointTest('b', 'ezr')) + def test_Jtvec_adjointTest_exi_Bform(self): + self.assertTrue(adjointTest('b', 'exi')) + def test_Jtvec_adjointTest_eyi_Bform(self): + self.assertTrue(adjointTest('b', 'eyi')) + def test_Jtvec_adjointTest_ezi_Bform(self): + self.assertTrue(adjointTest('b', 'ezi')) + def test_Jtvec_adjointTest_bxr_Bform(self): + self.assertTrue(adjointTest('b', 'bxr')) + def test_Jtvec_adjointTest_byr_Bform(self): + self.assertTrue(adjointTest('b', 'byr')) + def test_Jtvec_adjointTest_bzr_Bform(self): + self.assertTrue(adjointTest('b', 'bzr')) + def test_Jtvec_adjointTest_bxi_Bform(self): + self.assertTrue(adjointTest('b', 'bxi')) + def test_Jtvec_adjointTest_byi_Bform(self): + self.assertTrue(adjointTest('b', 'byi')) + def test_Jtvec_adjointTest_bzi_Bform(self): + self.assertTrue(adjointTest('b', 'bzi')) - def test_Jvec_jxr_Jform(self): - self.assertTrue(derivTest('j', 'jxr')) - def test_Jvec_jyr_Jform(self): - self.assertTrue(derivTest('j', 'jyr')) - def test_Jvec_jzr_Jform(self): - self.assertTrue(derivTest('j', 'jzr')) - def test_Jvec_jxi_Jform(self): - self.assertTrue(derivTest('j', 'jxi')) - def test_Jvec_jyi_Jform(self): - self.assertTrue(derivTest('j', 'jyi')) - def test_Jvec_jzi_Jform(self): - self.assertTrue(derivTest('j', 'jzi')) + if testHJ: + def test_Jtvec_adjointTest_jxr_Jform(self): + self.assertTrue(adjointTest('j', 'jxr')) + def test_Jtvec_adjointTest_jyr_Jform(self): + self.assertTrue(adjointTest('j', 'jyr')) + def test_Jtvec_adjointTest_jzr_Jform(self): + self.assertTrue(adjointTest('j', 'jzr')) + def test_Jtvec_adjointTest_jxi_Jform(self): + self.assertTrue(adjointTest('j', 'jxi')) + def test_Jtvec_adjointTest_jyi_Jform(self): + self.assertTrue(adjointTest('j', 'jyi')) + def test_Jtvec_adjointTest_jzi_Jform(self): + self.assertTrue(adjointTest('j', 'jzi')) - def test_Jvec_hxr_Jform(self): - self.assertTrue(derivTest('j', 'hxr')) - def test_Jvec_hyr_Jform(self): - self.assertTrue(derivTest('j', 'hyr')) - def test_Jvec_hzr_Jform(self): - self.assertTrue(derivTest('j', 'hzr')) - def test_Jvec_hxi_Jform(self): - self.assertTrue(derivTest('j', 'hxi')) - def test_Jvec_hyi_Jform(self): - self.assertTrue(derivTest('j', 'hyi')) - def test_Jvec_hzi_Jform(self): - self.assertTrue(derivTest('j', 'hzi')) + def test_Jtvec_adjointTest_hxr_Jform(self): + self.assertTrue(adjointTest('j', 'hxr')) + def test_Jtvec_adjointTest_hyr_Jform(self): + self.assertTrue(adjointTest('j', 'hyr')) + def test_Jtvec_adjointTest_hzr_Jform(self): + self.assertTrue(adjointTest('j', 'hzr')) + def test_Jtvec_adjointTest_hxi_Jform(self): + self.assertTrue(adjointTest('j', 'hxi')) + def test_Jtvec_adjointTest_hyi_Jform(self): + self.assertTrue(adjointTest('j', 'hyi')) + def test_Jtvec_adjointTest_hzi_Jform(self): + self.assertTrue(adjointTest('j', 'hzi')) - def test_Jvec_hxr_Hform(self): - self.assertTrue(derivTest('h', 'hxr')) - def test_Jvec_hyr_Hform(self): - self.assertTrue(derivTest('h', 'hyr')) - def test_Jvec_hzr_Hform(self): - self.assertTrue(derivTest('h', 'hzr')) - def test_Jvec_hxi_Hform(self): - self.assertTrue(derivTest('h', 'hxi')) - def test_Jvec_hyi_Hform(self): - self.assertTrue(derivTest('h', 'hyi')) - def test_Jvec_hzi_Hform(self): - self.assertTrue(derivTest('h', 'hzi')) + def test_Jtvec_adjointTest_hxr_Hform(self): + self.assertTrue(adjointTest('h', 'hxr')) + def test_Jtvec_adjointTest_hyr_Hform(self): + self.assertTrue(adjointTest('h', 'hyr')) + def test_Jtvec_adjointTest_hzr_Hform(self): + self.assertTrue(adjointTest('h', 'hzr')) + def test_Jtvec_adjointTest_hxi_Hform(self): + self.assertTrue(adjointTest('h', 'hxi')) + def test_Jtvec_adjointTest_hyi_Hform(self): + self.assertTrue(adjointTest('h', 'hyi')) + def test_Jtvec_adjointTest_hzi_Hform(self): + self.assertTrue(adjointTest('h', 'hzi')) - def test_Jvec_hxr_Hform(self): - self.assertTrue(derivTest('h', 'jxr')) - def test_Jvec_hyr_Hform(self): - self.assertTrue(derivTest('h', 'jyr')) - def test_Jvec_hzr_Hform(self): - self.assertTrue(derivTest('h', 'jzr')) - def test_Jvec_hxi_Hform(self): - self.assertTrue(derivTest('h', 'jxi')) - def test_Jvec_hyi_Hform(self): - self.assertTrue(derivTest('h', 'jyi')) - def test_Jvec_hzi_Hform(self): - self.assertTrue(derivTest('h', 'jzi')) - - def test_Jtvec_adjointTest_jxr_Jform(self): - self.assertTrue(adjointTest('j', 'jxr')) - def test_Jtvec_adjointTest_jyr_Jform(self): - self.assertTrue(adjointTest('j', 'jyr')) - def test_Jtvec_adjointTest_jzr_Jform(self): - self.assertTrue(adjointTest('j', 'jzr')) - def test_Jtvec_adjointTest_jxi_Jform(self): - self.assertTrue(adjointTest('j', 'jxi')) - def test_Jtvec_adjointTest_jyi_Jform(self): - self.assertTrue(adjointTest('j', 'jyi')) - def test_Jtvec_adjointTest_jzi_Jform(self): - self.assertTrue(adjointTest('j', 'jzi')) - - def test_Jtvec_adjointTest_hxr_Jform(self): - self.assertTrue(adjointTest('j', 'hxr')) - def test_Jtvec_adjointTest_hyr_Jform(self): - self.assertTrue(adjointTest('j', 'hyr')) - def test_Jtvec_adjointTest_hzr_Jform(self): - self.assertTrue(adjointTest('j', 'hzr')) - def test_Jtvec_adjointTest_hxi_Jform(self): - self.assertTrue(adjointTest('j', 'hxi')) - def test_Jtvec_adjointTest_hyi_Jform(self): - self.assertTrue(adjointTest('j', 'hyi')) - def test_Jtvec_adjointTest_hzi_Jform(self): - self.assertTrue(adjointTest('j', 'hzi')) - - def test_Jtvec_adjointTest_hxr_Hform(self): - self.assertTrue(adjointTest('h', 'hxr')) - def test_Jtvec_adjointTest_hyr_Hform(self): - self.assertTrue(adjointTest('h', 'hyr')) - def test_Jtvec_adjointTest_hzr_Hform(self): - self.assertTrue(adjointTest('h', 'hzr')) - def test_Jtvec_adjointTest_hxi_Hform(self): - self.assertTrue(adjointTest('h', 'hxi')) - def test_Jtvec_adjointTest_hyi_Hform(self): - self.assertTrue(adjointTest('h', 'hyi')) - def test_Jtvec_adjointTest_hzi_Hform(self): - self.assertTrue(adjointTest('h', 'hzi')) - - def test_Jtvec_adjointTest_hxr_Hform(self): - self.assertTrue(adjointTest('h', 'jxr')) - def test_Jtvec_adjointTest_hyr_Hform(self): - self.assertTrue(adjointTest('h', 'jyr')) - def test_Jtvec_adjointTest_hzr_Hform(self): - self.assertTrue(adjointTest('h', 'jzr')) - def test_Jtvec_adjointTest_hxi_Hform(self): - self.assertTrue(adjointTest('h', 'jxi')) - def test_Jtvec_adjointTest_hyi_Hform(self): - self.assertTrue(adjointTest('h', 'jyi')) - def test_Jtvec_adjointTest_hzi_Hform(self): - self.assertTrue(adjointTest('h', 'jzi')) + def test_Jtvec_adjointTest_hxr_Hform(self): + self.assertTrue(adjointTest('h', 'jxr')) + def test_Jtvec_adjointTest_hyr_Hform(self): + self.assertTrue(adjointTest('h', 'jyr')) + def test_Jtvec_adjointTest_hzr_Hform(self): + self.assertTrue(adjointTest('h', 'jzr')) + def test_Jtvec_adjointTest_hxi_Hform(self): + self.assertTrue(adjointTest('h', 'jxi')) + def test_Jtvec_adjointTest_hyi_Hform(self): + self.assertTrue(adjointTest('h', 'jyi')) + def test_Jtvec_adjointTest_hzi_Hform(self): + self.assertTrue(adjointTest('h', 'jzi')) - def test_EB_CrossCheck_exr_Eform(self): - self.assertTrue(crossCheckTest('e', 'exr')) - def test_EB_CrossCheck_eyr_Eform(self): - self.assertTrue(crossCheckTest('e', 'eyr')) - def test_EB_CrossCheck_ezr_Eform(self): - self.assertTrue(crossCheckTest('e', 'ezr')) - def test_EB_CrossCheck_exi_Eform(self): - self.assertTrue(crossCheckTest('e', 'exi')) - def test_EB_CrossCheck_eyi_Eform(self): - self.assertTrue(crossCheckTest('e', 'eyi')) - def test_EB_CrossCheck_ezi_Eform(self): - self.assertTrue(crossCheckTest('e', 'ezi')) + if testCrossCheck: + if testEB: + def test_EB_CrossCheck_exr_Eform(self): + self.assertTrue(crossCheckTest('e', 'exr')) + def test_EB_CrossCheck_eyr_Eform(self): + self.assertTrue(crossCheckTest('e', 'eyr')) + def test_EB_CrossCheck_ezr_Eform(self): + self.assertTrue(crossCheckTest('e', 'ezr')) + def test_EB_CrossCheck_exi_Eform(self): + self.assertTrue(crossCheckTest('e', 'exi')) + def test_EB_CrossCheck_eyi_Eform(self): + self.assertTrue(crossCheckTest('e', 'eyi')) + def test_EB_CrossCheck_ezi_Eform(self): + self.assertTrue(crossCheckTest('e', 'ezi')) - def test_EB_CrossCheck_bxr_Eform(self): - self.assertTrue(crossCheckTest('e', 'bxr')) - def test_EB_CrossCheck_byr_Eform(self): - self.assertTrue(crossCheckTest('e', 'byr')) - # def test_EB_CrossCheck_bzr_Eform(self): - # self.assertTrue(crossCheckTest('e', 'bzr')) # Doesn't make sense to test this for p-s approach - def test_EB_CrossCheck_bxi_Eform(self): - self.assertTrue(crossCheckTest('e', 'bxi')) - def test_EB_CrossCheck_byi_Eform(self): - self.assertTrue(crossCheckTest('e', 'byi')) - def test_EB_CrossCheck_bzi_Eform(self): - self.assertTrue(crossCheckTest('e', 'bzi')) + def test_EB_CrossCheck_bxr_Eform(self): + self.assertTrue(crossCheckTest('e', 'bxr')) + def test_EB_CrossCheck_byr_Eform(self): + self.assertTrue(crossCheckTest('e', 'byr')) + def test_EB_CrossCheck_bzr_Eform(self): + self.assertTrue(crossCheckTest('e', 'bzr')) + def test_EB_CrossCheck_bxi_Eform(self): + self.assertTrue(crossCheckTest('e', 'bxi')) + def test_EB_CrossCheck_byi_Eform(self): + self.assertTrue(crossCheckTest('e', 'byi')) + def test_EB_CrossCheck_bzi_Eform(self): + self.assertTrue(crossCheckTest('e', 'bzi')) - def test_HJ_CrossCheck_jxr_Jform(self): - self.assertTrue(crossCheckTest('j', 'jxr')) - def test_HJ_CrossCheck_jyr_Jform(self): - self.assertTrue(crossCheckTest('j', 'jyr')) - def test_HJ_CrossCheck_jzr_Jform(self): - self.assertTrue(crossCheckTest('j', 'jzr')) - def test_HJ_CrossCheck_jxi_Jform(self): - self.assertTrue(crossCheckTest('j', 'jxi')) - def test_HJ_CrossCheck_jyi_Jform(self): - self.assertTrue(crossCheckTest('j', 'jyi')) - def test_HJ_CrossCheck_jzi_Jform(self): - self.assertTrue(crossCheckTest('j', 'jzi')) + if testHJ: + def test_HJ_CrossCheck_jxr_Jform(self): + self.assertTrue(crossCheckTest('j', 'jxr')) + def test_HJ_CrossCheck_jyr_Jform(self): + self.assertTrue(crossCheckTest('j', 'jyr')) + def test_HJ_CrossCheck_jzr_Jform(self): + self.assertTrue(crossCheckTest('j', 'jzr')) + def test_HJ_CrossCheck_jxi_Jform(self): + self.assertTrue(crossCheckTest('j', 'jxi')) + def test_HJ_CrossCheck_jyi_Jform(self): + self.assertTrue(crossCheckTest('j', 'jyi')) + def test_HJ_CrossCheck_jzi_Jform(self): + self.assertTrue(crossCheckTest('j', 'jzi')) + + def test_HJ_CrossCheck_hxr_Jform(self): + self.assertTrue(crossCheckTest('j', 'hxr')) + def test_HJ_CrossCheck_hyr_Jform(self): + self.assertTrue(crossCheckTest('j', 'hyr')) + def test_HJ_CrossCheck_hzr_Jform(self): + self.assertTrue(crossCheckTest('j', 'hzr')) + def test_HJ_CrossCheck_hxi_Jform(self): + self.assertTrue(crossCheckTest('j', 'hxi')) + def test_HJ_CrossCheck_hyi_Jform(self): + self.assertTrue(crossCheckTest('j', 'hyi')) + def test_HJ_CrossCheck_hzi_Jform(self): + self.assertTrue(crossCheckTest('j', 'hzi')) - def test_HJ_CrossCheck_hxr_Jform(self): - self.assertTrue(crossCheckTest('j', 'hxr')) - def test_HJ_CrossCheck_hyr_Jform(self): - self.assertTrue(crossCheckTest('j', 'hyr')) - # def test_HJ_CrossCheck_hzr_Jform(self): - # self.assertTrue(crossCheckTest('j', 'hzr')) # Doesn't make sense to test this for p-s approach - def test_HJ_CrossCheck_hxi_Jform(self): - self.assertTrue(crossCheckTest('j', 'hxi')) - def test_HJ_CrossCheck_hyi_Jform(self): - self.assertTrue(crossCheckTest('j', 'hyi')) - def test_HJ_CrossCheck_hzi_Jform(self): - self.assertTrue(crossCheckTest('j', 'hzi')) if __name__ == '__main__': unittest.main() diff --git a/simpegEM/Tests/test_FDEMCasing.py b/simpegEM/Tests/test_FDEMCasing.py new file mode 100644 index 00000000..12eedb14 --- /dev/null +++ b/simpegEM/Tests/test_FDEMCasing.py @@ -0,0 +1,61 @@ +from SimPEG import Tests, Utils, np +import simpegEM.Analytics.FDEMcasing as Casing +import unittest +from scipy.constants import mu_0 + + +n = 50. +freq = 1. +a = 5e-2 +b = a + 1e-2 +sigma = np.r_[10., 5.5e6, 1e-1] +mu = mu_0*np.r_[1.,100.,1.] +srcloc = np.r_[0., 0., 0.] +xobs = np.random.rand(n)+10. +yobs = np.zeros(n) +zobs = np.random.randn(n) +plotit = False + +def CasingMagDipoleDeriv_r(x): + obsloc = np.vstack([x, yobs, zobs]).T + + f = Casing._getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu) + g = Utils.sdiag(Casing._getCasingHertzMagDipoleDeriv_r(srcloc,obsloc,freq,sigma,a,b,mu)) + + return f,g + +def CasingMagDipoleDeriv_z(z): + obsloc = np.vstack([xobs, yobs, z]).T + + f = Casing._getCasingHertzMagDipole(srcloc,obsloc,freq,sigma,a,b,mu) + g = Utils.sdiag(Casing._getCasingHertzMagDipoleDeriv_z(srcloc,obsloc,freq,sigma,a,b,mu)) + + return f,g + +def CasingMagDipole2Deriv_z_r(x): + obsloc = np.vstack([x, yobs, zobs]).T + + f = Casing._getCasingHertzMagDipoleDeriv_z(srcloc,obsloc,freq,sigma,a,b,mu) + g = Utils.sdiag(Casing._getCasingHertzMagDipole2Deriv_z_r(srcloc,obsloc,freq,sigma,a,b,mu)) + + return f,g + +def CasingMagDipole2Deriv_z_z(z): + obsloc = np.vstack([xobs, yobs, z]).T + + f = Casing._getCasingHertzMagDipoleDeriv_z(srcloc,obsloc,freq,sigma,a,b,mu) + g = Utils.sdiag(Casing._getCasingHertzMagDipole2Deriv_z_z(srcloc,obsloc,freq,sigma,a,b,mu)) + + return f,g + + + +class Casing_DerivTest(unittest.TestCase): + Tests.checkDerivative(CasingMagDipoleDeriv_r,np.ones(n)*10+np.random.randn(n),plotIt=False) + Tests.checkDerivative(CasingMagDipoleDeriv_z,np.random.randn(n),plotIt=False) + Tests.checkDerivative(CasingMagDipole2Deriv_z_r,np.ones(n)*10+np.random.randn(n),plotIt=False) + Tests.checkDerivative(CasingMagDipole2Deriv_z_z,np.random.randn(n),plotIt=False) + + +if __name__ == '__main__': + unittest.main() \ No newline at end of file diff --git a/simpegEM/Tests/test_FDEM_analytics.py b/simpegEM/Tests/test_FDEM_analytics.py index cae6d006..c0fc4210 100644 --- a/simpegEM/Tests/test_FDEM_analytics.py +++ b/simpegEM/Tests/test_FDEM_analytics.py @@ -4,6 +4,7 @@ import simpegEM as EM from scipy.constants import mu_0 plotIt = False +freq = 1e2 class FDEM_analyticTests(unittest.TestCase): @@ -22,9 +23,9 @@ class FDEM_analyticTests(unittest.TestCase): x = np.linspace(-10,10,5) XYZ = Utils.ndgrid(x,np.r_[0],np.r_[0]) rxList = EM.FDEM.RxFDEM(XYZ, 'exi') - Tx0 = EM.FDEM.TxFDEM(np.r_[0.,0.,0.], 'VMD', 1e2, [rxList]) + Src0 = EM.FDEM.SrcFDEM_MagDipole([rxList],loc=np.r_[0.,0.,0.], freq=freq) - survey = EM.FDEM.SurveyFDEM([Tx0]) + survey = EM.FDEM.SurveyFDEM([Src0]) prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping) prb.pair(survey) @@ -43,7 +44,7 @@ class FDEM_analyticTests(unittest.TestCase): self.prb = prb self.mesh = mesh self.m = m - self.Tx0 = Tx0 + self.Src0 = Src0 self.sig = sig def test_Transect(self): @@ -51,20 +52,19 @@ class FDEM_analyticTests(unittest.TestCase): u = self.prb.fields(self.m) - bfz = self.mesh.r(u[self.Tx0, 'b'],'F','Fz','M') - + bfz = self.mesh.r(u[self.Src0, 'b'],'F','Fz','M') x = np.linspace(-55,55,12) XYZ = Utils.ndgrid(x,np.r_[0],np.r_[0]) P = self.mesh.getInterpolationMat(XYZ, 'Fz') - an = EM.Analytics.FDEM.hzAnalyticDipoleF(x, self.Tx0.freq, self.sig) + an = EM.Analytics.FDEM.hzAnalyticDipoleF(x, self.Src0.freq, self.sig) - diff = np.log10(np.abs(P*np.imag(u[self.Tx0, 'b']) - mu_0*np.imag(an))) + diff = np.log10(np.abs(P*np.imag(u[self.Src0, 'b']) - mu_0*np.imag(an))) if plotIt: import matplotlib.pyplot as plt - plt.plot(x,np.log10(np.abs(P*np.imag(u[self.Tx0, 'b'])))) + plt.plot(x,np.log10(np.abs(P*np.imag(u[self.Src0, 'b'])))) plt.plot(x,np.log10(np.abs(mu_0*np.imag(an))), 'r') plt.plot(x,diff,'g') plt.show() diff --git a/simpegEM/Tests/test_FieldsObject.py b/simpegEM/Tests/test_FieldsObject.py deleted file mode 100644 index f28d1ad5..00000000 --- a/simpegEM/Tests/test_FieldsObject.py +++ /dev/null @@ -1,103 +0,0 @@ -import unittest -from SimPEG import * -import simpegEM as EM - -class FieldsTest(unittest.TestCase): - - def setUp(self): - mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10,11,12]],[0,0,-30]) - x = np.linspace(5,10,3) - XYZ = Utils.ndgrid(x,x,np.r_[0.]) - txLoc = np.r_[0,0,0.] - rxList0 = EM.FDEM.RxFDEM(XYZ, 'exi') - Tx0 = EM.FDEM.TxFDEM(txLoc, 'VMD', 3., [rxList0]) - rxList1 = EM.FDEM.RxFDEM(XYZ, 'bxi') - Tx1 = EM.FDEM.TxFDEM(txLoc, 'VMD', 3., [rxList1]) - rxList2 = EM.FDEM.RxFDEM(XYZ, 'bxi') - Tx2 = EM.FDEM.TxFDEM(txLoc, 'VMD', 2., [rxList2]) - rxList3 = EM.FDEM.RxFDEM(XYZ, 'bxi') - Tx3 = EM.FDEM.TxFDEM(txLoc, 'VMD', 2., [rxList3]) - Tx4 = EM.FDEM.TxFDEM(txLoc, 'VMD', 1., [rxList0, rxList1, rxList2, rxList3]) - txList = [Tx0,Tx1,Tx2,Tx3,Tx4] - survey = EM.FDEM.SurveyFDEM(txList) - self.F = EM.FDEM.FieldsFDEM(mesh, survey) - self.Tx0 = Tx0 - self.Tx1 = Tx1 - self.mesh = mesh - self.XYZ = XYZ - - def test_SetGet(self): - F = self.F - for freq in F.survey.freqs: - nFreq = F.survey.nTxByFreq[freq] - Txs = F.survey.getTransmitters(freq) - e = np.random.rand(F.mesh.nE, nFreq) - F[Txs, 'e'] = e - b = np.random.rand(F.mesh.nF, nFreq) - F[Txs, 'b'] = b - if nFreq == 1: - F[Txs, 'b'] = Utils.mkvc(b) - if e.shape[1] == 1: - e, b = Utils.mkvc(e), Utils.mkvc(b) - self.assertTrue(np.all(F[Txs, 'e'] == e)) - self.assertTrue(np.all(F[Txs, 'b'] == b)) - F[Txs] = {'b':b,'e':e} - self.assertTrue(np.all(F[Txs, 'e'] == e)) - self.assertTrue(np.all(F[Txs, 'b'] == b)) - - lastFreq = F[Txs] - self.assertTrue(type(lastFreq) is dict) - self.assertTrue(sorted([k for k in lastFreq]) == ['b','e']) - self.assertTrue(np.all(lastFreq['b'] == b)) - self.assertTrue(np.all(lastFreq['e'] == e)) - - Tx_f3 = F.survey.getTransmitters(3.) - self.assertTrue(F[Tx_f3,'b'].shape == (F.mesh.nF, 2)) - - b = np.random.rand(F.mesh.nF, 2) - Tx_f0 = F.survey.getTransmitters(self.Tx0.freq) - F[Tx_f0,'b'] = b - self.assertTrue(F[self.Tx0]['b'].shape == (F.mesh.nF,)) - self.assertTrue(F[self.Tx0,'b'].shape == (F.mesh.nF,)) - self.assertTrue(np.all(F[self.Tx0,'b'] == b[:,0])) - self.assertTrue(np.all(F[self.Tx1,'b'] == b[:,1])) - - def test_assertions(self): - freq = self.F.survey.freqs[0] - Txs = self.F.survey.getTransmitters(freq) - bWrongSize = np.random.rand(self.F.mesh.nE, self.F.survey.nTxByFreq[freq]) - def fun(): self.F[Txs, 'b'] = bWrongSize - self.assertRaises(ValueError, fun) - def fun(): self.F[-999.] - self.assertRaises(KeyError, fun) - def fun(): self.F['notRight'] - self.assertRaises(KeyError, fun) - def fun(): self.F[Txs,'notThere'] - self.assertRaises(KeyError, fun) - - def test_FieldProjections(self): - F = self.F - for freq in F.survey.freqs: - nFreq = F.survey.nTxByFreq[freq] - Txs = F.survey.getTransmitters(freq) - e = np.random.rand(F.mesh.nE, nFreq) - b = np.random.rand(F.mesh.nF, nFreq) - F[Txs] = {'b':b,'e':e} - - Txs = F.survey.getTransmitters(freq) - for ii, tx in enumerate(Txs): - for jj, rx in enumerate(tx.rxList): - dat = rx.projectFields(tx, self.mesh, F) - self.assertTrue(dat.dtype == float) - fieldType = rx.projField - u = {'b':b[:,ii], 'e': e[:,ii]}[fieldType] - real_or_imag = rx.projComp - u = getattr(u, real_or_imag) - gloc = rx.projGLoc - d = self.mesh.getInterpolationMat(self.XYZ, gloc)*u - self.assertTrue(np.all(dat == d)) - - - -if __name__ == '__main__': - unittest.main() diff --git a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py index 9f6050d5..0efcce2d 100644 --- a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py +++ b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py @@ -3,6 +3,7 @@ from SimPEG import * import simpegEM as EM plotIt = False +tol = 1e-6 class TDEM_bDerivTests(unittest.TestCase): @@ -22,9 +23,9 @@ class TDEM_bDerivTests(unittest.TestCase): rxOffset = 40. rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz') - tx = EM.TDEM.TxTDEM(np.array([0., 0., 0.]), 'VMD_MVP', [rx]) + src = EM.TDEM.SrcTDEM_VMD_MVP([rx], loc=np.array([0., 0., 0.])) - survey = EM.TDEM.SurveyTDEM([tx]) + survey = EM.TDEM.SurveyTDEM([src]) self.prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping) # self.prb.timeSteps = [1e-5] @@ -60,7 +61,7 @@ class TDEM_bDerivTests(unittest.TestCase): self.assertLess(np.linalg.norm(V1-V2)/np.linalg.norm(V2), 1.e-6) V1 = Ahu[:,'e',1] - self.assertLess(np.linalg.norm(V1), 1.e-6) + return np.linalg.norm(V1) < 1.e-6 for i in range(2,prb.nT): @@ -74,7 +75,7 @@ class TDEM_bDerivTests(unittest.TestCase): V1 = Ahu[:,'e',i] V2 = prb.MeSigma*u[:,'e',i] # print np.linalg.norm(V1), np.linalg.norm(V2) - self.assertLess(np.linalg.norm(V1)/np.linalg.norm(V2), 1.e-6) + return np.linalg.norm(V1)/np.linalg.norm(V2), 1.e-6 def test_AhVecVSMat_OneTS(self): @@ -120,7 +121,7 @@ class TDEM_bDerivTests(unittest.TestCase): u1 = prb.solveAh(sigma,f).tovec().flatten() u2 = sp.linalg.spsolve(A.tocsr(),f.tovec()) - self.assertLess(np.linalg.norm(u1-u2),1e-8) + self.assertTrue(np.linalg.norm(u1-u2)<1e-8) def test_solveAhVsAhVec(self): @@ -156,7 +157,7 @@ class TDEM_bDerivTests(unittest.TestCase): derChk = lambda m: [self.prb._AhVec(m, f).tovec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).tovec()] print '\ntest_DerivG' passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20) - self.assertTrue(passed) + return passed def test_Deriv_dUdM(self): @@ -171,8 +172,7 @@ class TDEM_bDerivTests(unittest.TestCase): derChk = lambda m: [self.prb.fields(m).tovec(), lambda mx: -prb.solveAh(sigma, prb.Gvec(sigma, mx, u=f)).tovec()] print '\n' print 'test_Deriv_dUdM' - passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20) - self.assertTrue(passed) + Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20) def test_Deriv_J(self): @@ -188,8 +188,7 @@ class TDEM_bDerivTests(unittest.TestCase): derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(sigma, mx)] print '\n' print 'test_Deriv_J' - passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20) - self.assertTrue(passed) + Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20) def test_projectAdjoint(self): prb = self.prb @@ -208,7 +207,7 @@ class TDEM_bDerivTests(unittest.TestCase): V1 = d_vec.dot(survey.projectFieldsDeriv(None, v=f).tovec()) V2 = f.tovec().dot(survey.projectFieldsDeriv(None, v=d, adjoint=True).tovec()) - self.assertLess((V1-V2)/np.abs(V1), 1e-6) + self.assertTrue((V1-V2)/np.abs(V1) < tol) def test_adjointAhVsAht(self): prb = self.prb @@ -227,7 +226,7 @@ class TDEM_bDerivTests(unittest.TestCase): V1 = f2.tovec().dot(prb._AhVec(sigma, f1).tovec()) V2 = f1.tovec().dot(prb._AhtVec(sigma, f2).tovec()) - self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6) + self.assertTrue(np.abs(V1-V2)/np.abs(V1) < tol) # def test_solveAhtVsAhtVec(self): # prb = self.prb @@ -292,7 +291,7 @@ class TDEM_bDerivTests(unittest.TestCase): V1 = m.dot(prb.Gtvec(sigma, v, u)) V2 = v.tovec().dot(prb.Gvec(sigma, m, u).tovec()) - self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6) + self.assertTrue(np.abs(V1-V2)/np.abs(V1) < tol) def test_adjointJvecVsJtvec(self): mesh = self.mesh @@ -304,8 +303,10 @@ class TDEM_bDerivTests(unittest.TestCase): V1 = d.dot(prb.Jvec(sigma, m)) V2 = m.dot(prb.Jtvec(sigma, d)) - print 'AdjointTest', V1, V2 - self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6) + passed = np.abs(V1-V2)/np.abs(V1) < tol + print 'AdjointTest', V1, V2, passed + self.assertTrue(passed) + diff --git a/simpegEM/Tests/test_TDEM_b_MultiTx_DerivAdjoint.py b/simpegEM/Tests/test_TDEM_b_MultiSrc_DerivAdjoint.py similarity index 80% rename from simpegEM/Tests/test_TDEM_b_MultiTx_DerivAdjoint.py rename to simpegEM/Tests/test_TDEM_b_MultiSrc_DerivAdjoint.py index e17343e0..dc613a2a 100644 --- a/simpegEM/Tests/test_TDEM_b_MultiTx_DerivAdjoint.py +++ b/simpegEM/Tests/test_TDEM_b_MultiSrc_DerivAdjoint.py @@ -22,11 +22,11 @@ class TDEM_bDerivTests(unittest.TestCase): rxOffset = 40. rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz') - tx = EM.TDEM.TxTDEM(np.array([0., 0., 0.]), 'VMD_MVP', [rx]) + src = EM.TDEM.SrcTDEM_VMD_MVP( [rx], loc=np.array([0., 0., 0.])) rx2 = EM.TDEM.RxTDEM(np.array([[rxOffset-10, 0., 0.]]), np.logspace(-5,-4, 25), 'bz') - tx2 = EM.TDEM.TxTDEM(np.array([0., 0., 0.]), 'VMD_MVP', [rx2]) + src2 = EM.TDEM.SrcTDEM_VMD_MVP( [rx2], loc=np.array([0., 0., 0.])) - survey = EM.TDEM.SurveyTDEM([tx,tx2]) + survey = EM.TDEM.SurveyTDEM([src,src2]) self.prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping) # self.prb.timeSteps = [1e-5] @@ -60,8 +60,7 @@ class TDEM_bDerivTests(unittest.TestCase): derChk = lambda m: [self.prb._AhVec(m, f).tovec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).tovec()] print '\ntest_DerivG' - passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20) - self.assertTrue(passed) + Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20) def test_Deriv_dUdM(self): @@ -76,8 +75,7 @@ class TDEM_bDerivTests(unittest.TestCase): derChk = lambda m: [self.prb.fields(m).tovec(), lambda mx: -prb.solveAh(sigma, prb.Gvec(sigma, mx, u=f)).tovec()] print '\n' print 'test_Deriv_dUdM' - passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20) - self.assertTrue(passed) + Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20) def test_Deriv_J(self): @@ -93,28 +91,27 @@ class TDEM_bDerivTests(unittest.TestCase): derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(sigma, mx)] print '\n' print 'test_Deriv_J' - passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20) - self.assertTrue(passed) + Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20) def test_projectAdjoint(self): prb = self.prb survey = prb.survey - nTx = survey.nTx + nSrc = survey.nSrc mesh = self.mesh # Generate random fields and data f = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) for i in range(prb.nT): - f[:,'b',i] = np.random.rand(mesh.nF, nTx) - f[:,'e',i] = np.random.rand(mesh.nE, nTx) + f[:,'b',i] = np.random.rand(mesh.nF, nSrc) + f[:,'e',i] = np.random.rand(mesh.nE, nSrc) d_vec = np.random.rand(survey.nD) d = Survey.Data(survey,v=d_vec) # Check that d.T*Q*f = f.T*Q.T*d V1 = d_vec.dot(survey.projectFieldsDeriv(None, v=f).tovec()) - V2 = f.tovec().dot(survey.projectFieldsDeriv(None, v=d, adjoint=True).tovec()) + V2 = np.sum((f.tovec())*(survey.projectFieldsDeriv(None, v=d, adjoint=True).tovec())) - self.assertLess((V1-V2)/np.abs(V1), 1e-6) + self.assertTrue((V1-V2)/np.abs(V1) < 1e-6) def test_adjointGvecVsGtvec(self): mesh = self.mesh @@ -134,8 +131,8 @@ class TDEM_bDerivTests(unittest.TestCase): v[:,'e',i] = np.random.rand(mesh.nE, 2) V1 = m.dot(prb.Gtvec(sigma, v, u)) - V2 = v.tovec().dot(prb.Gvec(sigma, m, u).tovec()) - self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6) + V2 = np.sum(v.tovec()*prb.Gvec(sigma, m, u).tovec()) + self.assertTrue(np.abs(V1-V2)/np.abs(V1) <1e-6) def test_adjointJvecVsJtvec(self): mesh = self.mesh @@ -148,7 +145,7 @@ class TDEM_bDerivTests(unittest.TestCase): V1 = d.dot(prb.Jvec(sigma, m)) V2 = m.dot(prb.Jtvec(sigma, d)) print 'AdjointTest', V1, V2 - self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6) + self.assertTrue(np.abs(V1-V2)/np.abs(V1) < 1e-6) diff --git a/simpegEM/Tests/test_TDEM_combos.py b/simpegEM/Tests/test_TDEM_combos.py index 69f9e397..bca3c07e 100644 --- a/simpegEM/Tests/test_TDEM_combos.py +++ b/simpegEM/Tests/test_TDEM_combos.py @@ -4,7 +4,7 @@ import simpegEM as EM plotIt = False -def getProb(meshType='CYL',rxTypes='bx,bz',nTx=1): +def getProb(meshType='CYL',rxTypes='bx,bz',nSrc=1): cs = 5. ncx = 20 ncy = 6 @@ -19,12 +19,12 @@ def getProb(meshType='CYL',rxTypes='bx,bz',nTx=1): rxOffset = 40. - txs = [] - for ii in range(nTx): + srcs = [] + for ii in range(nSrc): rxs = [EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20 + ii), rxType) for rxType in rxTypes.split(',')] - txs += [EM.TDEM.TxTDEM(np.array([0., 0., 0.]), 'VMD_MVP', rxs)] + srcs += [EM.TDEM.SrcTDEM_VMD_MVP(rxs,np.array([0., 0., 0.]))] - survey = EM.TDEM.SurveyTDEM(txs) + survey = EM.TDEM.SurveyTDEM(srcs) prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping) # prb.timeSteps = [1e-5] @@ -68,8 +68,8 @@ class TDEM_bDerivTests(unittest.TestCase): def test_Jvec_bxbz(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx,bz'))) def test_Adjoint_bxbz(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx,bz'))) - def test_Jvec_bxbz_2tx(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx,bz',nTx=2))) - def test_Adjoint_bxbz_2tx(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx,bz',nTx=2))) + def test_Jvec_bxbz_2src(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx,bz',nSrc=2))) + def test_Adjoint_bxbz_2src(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx,bz',nSrc=2))) def test_Jvec_bxbzbz(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx,bz,bz'))) def test_Adjoint_bxbzbz(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx,bz,bz'))) diff --git a/simpegEM/Tests/test_TDEM_forward_Analytic.py b/simpegEM/Tests/test_TDEM_forward_Analytic.py index 4557c5cf..11db6ab8 100644 --- a/simpegEM/Tests/test_TDEM_forward_Analytic.py +++ b/simpegEM/Tests/test_TDEM_forward_Analytic.py @@ -28,9 +28,10 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=[1e-5, mapping = Maps.ExpMap(mesh) * Maps.Vertical1DMap(mesh) * actMap rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-5,-4, 21), 'bz') - tx = EM.TDEM.TxTDEM(np.array([0., 0., 0.]), 'VMD_MVP', [rx]) + src = EM.TDEM.SrcTDEM_VMD_MVP([rx], loc=np.array([0., 0., 0.])) + # src = EM.TDEM.SrcTDEM([rx], loc=np.array([0., 0., 0.])) - survey = EM.TDEM.SurveyTDEM([tx]) + survey = EM.TDEM.SurveyTDEM([src]) prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping) prb.Solver = MumpsSolver @@ -60,26 +61,26 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=[1e-5, class TDEM_bTests(unittest.TestCase): - def test_analitic_p2_CYL_50m(self): + def test_analytic_p2_CYL_50m(self): self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e+2) < 0.01) - def test_analitic_p1_CYL_50m(self): + def test_analytic_p1_CYL_50m(self): self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e+1) < 0.01) - def test_analitic_p0_CYL_50m(self): + def test_analytic_p0_CYL_50m(self): self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e+0) < 0.01) - def test_analitic_m1_CYL_50m(self): + def test_analytic_m1_CYL_50m(self): self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e-1) < 0.01) - def test_analitic_m2_CYL_50m(self): + def test_analytic_m2_CYL_50m(self): self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e-2) < 0.01) - def test_analitic_m3_CYL_50m(self): + def test_analytic_m3_CYL_50m(self): self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=50., sig_half=1e-3) < 0.02) - def test_analitic_p0_CYL_1m(self): + def test_analytic_p0_CYL_1m(self): self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=1.0, sig_half=1e+0) < 0.01) - def test_analitic_m1_CYL_1m(self): + def test_analytic_m1_CYL_1m(self): self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=1.0, sig_half=1e-1) < 0.01) - def test_analitic_m2_CYL_1m(self): + def test_analytic_m2_CYL_1m(self): self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=1.0, sig_half=1e-2) < 0.01) - def test_analitic_m3_CYL_1m(self): + def test_analytic_m3_CYL_1m(self): self.assertTrue(halfSpaceProblemAnaDiff('CYL', rxOffset=1.0, sig_half=1e-3) < 0.02) diff --git a/simpegEM/Utils/EMUtils.py b/simpegEM/Utils/EMUtils.py new file mode 100644 index 00000000..4a342acb --- /dev/null +++ b/simpegEM/Utils/EMUtils.py @@ -0,0 +1,49 @@ +import numpy as np +from scipy.constants import mu_0, epsilon_0 + +# useful params +def omega(freq): + """Angular frequency, omega""" + return 2.*np.pi*freq + +def k(freq, sigma, mu=mu_0, eps=epsilon_0): + """ Eq 1.47 - 1.49 in Ward and Hohmann """ + w = omega(freq) + alp = w * np.sqrt( mu*eps/2 * ( np.sqrt(1. + (sigma / (eps*w))**2 ) + 1) ) + beta = w * np.sqrt( mu*eps/2 * ( np.sqrt(1. + (sigma / (eps*w))**2 ) - 1) ) + return alp - 1j*beta + +# Constitutive relations +def e_from_j(prob,j): + eqLocs = prob._eqLocs + if eqLocs is 'FE': + MSigmaI = prob.MeSigmaI + elif eqLocs is 'EF': + MSigmaI = prob.MfRho + return MSigmaI*j + +def j_from_e(prob,e): + eqLocs = prob._eqLocs + if eqLocs is 'FE': + MSigma = prob.MeSigma + elif eqLocs is 'EF': + MSigma = prob.MfRhoI + return MSigma*e + +def b_from_h(prob,h): + eqLocs = prob._eqLocs + if eqLocs is 'FE': + MMu = prob.MfMuiI + elif eqLocs is 'EF': + MMu = prob.MeMu + return MMu*h + +def h_from_b(prob,b): + eqLocs = prob._eqLocs + if eqLocs is 'FE': + MMuI = prob.MfMui + elif eqLocs is 'EF': + MMuI = prob.MeMuI + return MMuI*b + + diff --git a/simpegEM/Utils/SrcUtils.py b/simpegEM/Utils/SrcUtils.py new file mode 100644 index 00000000..1827f6b2 --- /dev/null +++ b/simpegEM/Utils/SrcUtils.py @@ -0,0 +1,203 @@ +from SimPEG import * +from scipy.special import ellipk, ellipe +from scipy.constants import mu_0, pi + +def MagneticDipoleVectorPotential(srcLoc, obsLoc, component, moment=1., dipoleMoment=(0., 0., 1.), mu = mu_0): + """ + Calculate the vector potential of a set of magnetic dipoles + at given locations 'ref. ' + + :param numpy.ndarray srcLoc: Location of the source(s) (x, y, z) + :param numpy.ndarray,SimPEG.Mesh obsLoc: Where the potentials will be calculated (x, y, z) or a SimPEG Mesh + :param str,list component: The component to calculate - 'x', 'y', or 'z' if an array, or grid type if mesh, can be a list + :param numpy.ndarray dipoleMoment: The vector dipole moment + :rtype: numpy.ndarray + :return: The vector potential each dipole at each observation location + """ + #TODO: break this out! + + if type(component) in [list, tuple]: + out = range(len(component)) + for i, comp in enumerate(component): + out[i] = MagneticDipoleVectorPotential(srcLoc, obsLoc, comp, dipoleMoment=dipoleMoment) + return np.concatenate(out) + + if isinstance(obsLoc, Mesh.BaseMesh): + mesh = obsLoc + assert component in ['Ex','Ey','Ez','Fx','Fy','Fz'], "Components must be in: ['Ex','Ey','Ez','Fx','Fy','Fz']" + return MagneticDipoleVectorPotential(srcLoc, getattr(mesh,'grid'+component), component[1], dipoleMoment=dipoleMoment) + + if component == 'x': + dimInd = 0 + elif component == 'y': + dimInd = 1 + elif component == 'z': + dimInd = 2 + else: + raise ValueError('Invalid component') + + srcLoc = np.atleast_2d(srcLoc) + obsLoc = np.atleast_2d(obsLoc) + dipoleMoment = np.atleast_2d(dipoleMoment) + + nEdges = obsLoc.shape[0] + nSrc = srcLoc.shape[0] + + m = np.array(dipoleMoment).repeat(nEdges, axis=0) + A = np.empty((nEdges, nSrc)) + for i in range(nSrc): + dR = obsLoc - srcLoc[i, np.newaxis].repeat(nEdges, axis=0) + mCr = np.cross(m, dR) + r = np.sqrt((dR**2).sum(axis=1)) + A[:, i] = +(mu/(4*pi)) * mCr[:,dimInd]/(r**3) + if nSrc == 1: + return A.flatten() + return A + + +def MagneticDipoleFields(srcLoc, obsLoc, component, moment=1., mu = mu_0): + """ + Calculate the vector potential of a set of magnetic dipoles + at given locations 'ref. ' + + :param numpy.ndarray srcLoc: Location of the source(s) (x, y, z) + :param numpy.ndarray obsLoc: Where the potentials will be calculated (x, y, z) + :param str component: The component to calculate - 'x', 'y', or 'z' + :param numpy.ndarray moment: The vector dipole moment (vertical) + :rtype: numpy.ndarray + :return: The vector potential each dipole at each observation location + """ + + if component=='x': + dimInd = 0 + elif component=='y': + dimInd = 1 + elif component=='z': + dimInd = 2 + else: + raise ValueError('Invalid component') + + srcLoc = np.atleast_2d(srcLoc) + obsLoc = np.atleast_2d(obsLoc) + moment = np.atleast_2d(moment) + + nFaces = obsLoc.shape[0] + nSrc = srcLoc.shape[0] + + m = np.array(moment).repeat(nFaces, axis=0) + B = np.empty((nFaces, nSrc)) + for i in range(nSrc): + dR = obsLoc - srcLoc[i, np.newaxis].repeat(nFaces, axis=0) + r = np.sqrt((dR**2).sum(axis=1)) + if dimInd == 0: + B[:, i] = +(mu/(4*pi)) /(r**3) * (3*dR[:,2]*dR[:,0]/r**2) + elif dimInd == 1: + B[:, i] = +(mu/(4*pi)) /(r**3) * (3*dR[:,2]*dR[:,1]/r**2) + elif dimInd == 2: + B[:, i] = +(mu/(4*pi)) /(r**3) * (3*dR[:,2]**2/r**2-1) + else: + raise Exception("Not Implemented") + if nSrc == 1: + return B.flatten() + return B + + + +def MagneticLoopVectorPotential(srcLoc, obsLoc, component, radius, mu=mu_0): + """ + Calculate the vector potential of horizontal circular loop + at given locations + + :param numpy.ndarray srcLoc: Location of the source(s) (x, y, z) + :param numpy.ndarray,SimPEG.Mesh obsLoc: Where the potentials will be calculated (x, y, z) or a SimPEG Mesh + :param str,list component: The component to calculate - 'x', 'y', or 'z' if an array, or grid type if mesh, can be a list + :param numpy.ndarray I: Input current of the loop + :param numpy.ndarray radius: radius of the loop + :rtype: numpy.ndarray + :return: The vector potential each dipole at each observation location + """ + + if type(component) in [list, tuple]: + out = range(len(component)) + for i, comp in enumerate(component): + out[i] = MagneticLoopVectorPotential(srcLoc, obsLoc, comp, radius, mu) + return np.concatenate(out) + + if isinstance(obsLoc, Mesh.BaseMesh): + mesh = obsLoc + assert component in ['Ex','Ey','Ez','Fx','Fy','Fz'], "Components must be in: ['Ex','Ey','Ez','Fx','Fy','Fz']" + return MagneticLoopVectorPotential(srcLoc, getattr(mesh,'grid'+component), component[1], radius, mu) + + srcLoc = np.atleast_2d(srcLoc) + obsLoc = np.atleast_2d(obsLoc) + + n = obsLoc.shape[0] + nSrc = srcLoc.shape[0] + + if component=='z': + A = np.zeros((n, nSrc)) + if nSrc ==1: + return A.flatten() + return A + + else: + + A = np.zeros((n, nSrc)) + for i in range (nSrc): + x = obsLoc[:, 0] - srcLoc[i, 0] + y = obsLoc[:, 1] - srcLoc[i, 1] + z = obsLoc[:, 2] - srcLoc[i, 2] + r = np.sqrt(x**2 + y**2) + m = (4 * radius * r) / ((radius + r)**2 + z**2) + m[m > 1.] = 1. + # m might be slightly larger than 1 due to rounding errors + # but ellipke requires 0 <= m <= 1 + K = ellipk(m) + E = ellipe(m) + ind = (r > 0) & (m < 1) + # % 1/r singular at r = 0 and K(m) singular at m = 1 + Aphi = np.zeros(n) + # % Common factor is (mu * I) / pi with I = 1 and mu = 4e-7 * pi. + Aphi[ind] = 4e-7 / np.sqrt(m[ind]) * np.sqrt(radius / r[ind]) *((1. - m[ind] / 2.) * K[ind] - E[ind]) + if component == 'x': + A[ind, i] = Aphi[ind] * (-y[ind] / r[ind] ) + elif component == 'y': + A[ind, i] = Aphi[ind] * ( x[ind] / r[ind] ) + else: + raise ValueError('Invalid component') + + if nSrc == 1: + return A.flatten() + return A + +if __name__ == '__main__': + from SimPEG import Mesh + import matplotlib.pyplot as plt + cs = 20 + ncx, ncy, ncz = 41, 41, 40 + hx = np.ones(ncx)*cs + hy = np.ones(ncy)*cs + hz = np.ones(ncz)*cs + mesh = Mesh.TensorMesh([hx, hy, hz], 'CCC') + srcLoc = np.r_[0., 0., 0.] + Ax = MagneticLoopVectorPotential(srcLoc, mesh.gridEx, 'x', 200) + Ay = MagneticLoopVectorPotential(srcLoc, mesh.gridEy, 'y', 200) + Az = MagneticLoopVectorPotential(srcLoc, mesh.gridEz, 'z', 200) + A = np.r_[Ax, Ay, Az] + B0 = mesh.edgeCurl*A + J0 = mesh.edgeCurl.T*B0 + + # mesh.plotImage(A, vType = 'Ex') + # mesh.plotImage(A, vType = 'Ey') + + mesh.plotImage(B0, vType = 'Fx') + mesh.plotImage(B0, vType = 'Fy') + mesh.plotImage(B0, vType = 'Fz') + + # # mesh.plotImage(J0, vType = 'Ex') + # mesh.plotImage(J0, vType = 'Ey') + # mesh.plotImage(J0, vType = 'Ez') + + plt.show() + + diff --git a/simpegEM/Utils/__init__.py b/simpegEM/Utils/__init__.py index 2e736a8c..6e430cf9 100644 --- a/simpegEM/Utils/__init__.py +++ b/simpegEM/Utils/__init__.py @@ -1,3 +1,5 @@ # import Sources # import Ana # import Solver +import EMUtils +import SrcUtils \ No newline at end of file diff --git a/simpegEM/__init__.py b/simpegEM/__init__.py index 90a1dcf7..6a1ca774 100644 --- a/simpegEM/__init__.py +++ b/simpegEM/__init__.py @@ -2,5 +2,5 @@ import TDEM import FDEM import Base -import Sources import Analytics +import Utils