diff --git a/SimPEG/EM/Analytics/DC.py b/SimPEG/EM/Analytics/DC.py index 69d17090..4fb03fd3 100644 --- a/SimPEG/EM/Analytics/DC.py +++ b/SimPEG/EM/Analytics/DC.py @@ -1,7 +1,7 @@ import numpy as np from scipy.constants import mu_0, pi -def DCAnalytic(txloc, rxlocs, sigma, flag="wholespace"): +def DCAnalyticHalf(txloc, rxlocs, sigma, flag="wholespace"): """ Analytic solution for electric potential from a postive pole @@ -32,3 +32,82 @@ def DCAnalytic(txloc, rxlocs, sigma, flag="wholespace"): return phi +deg2rad = lambda deg: deg/180.*np.pi +rad2deg = lambda rad: rad*180./np.pi + +def DCAnalyticSphere(txloc, rxloc, xc, radius, sigma, sigma1, \ + flag = "sec", order=12): +# def DCSpherePointCurrent(txloc, rxloc, xc, radius, rho, rho1, \ +# flag = "sec", order=12): + """ + + Parameters: + + txloc (array) : current electrode location (x,y,z) + xc (float) : x center of depressed sphere + rxloc (array) : electrode locations + (Nx3 array, # of electrodes) + radius (float): radius of the sphere (m) + rho (float) : resistivity of the background (ohm-m) + rho1 (float) : resistivity of the sphere + flag (string) : "sec", "total", "prim" + (default="sec") + "sec": secondary potential only due to sphere + "prim": primary potential from the point source + "total": "sec"+"prim" + order (float) : maximum order of Legendre polynomial + (default=12) + + Written by Seogi Kang (skang@eos.ubc.ca) + Ph.D. Candidate of University of British Columbia, Canada + + """ + + Pleg = [] + # Compute Legendre Polynomial + for i in range(order): + Pleg.append(special.legendre(i, monic=0)) + + + rho = 1./sigma + rho1 = 1./sigma1 + + # Center of the sphere should be aligned in txloc in y-direction + yc = txloc[1] + xyz = np.c_[rxloc[:,0]-xc, rxloc[:,1]-yc, rxloc[:,2]] + r = np.sqrt( (xyz**2).sum(axis=1) ) + + x0 = abs(txloc[0]-xc) + + costheta = xyz[:,0]/r * (txloc[0]-xc)/x0 + phi = np.zeros_like(r) + R = (r**2+x0**2.-2.*r*x0*costheta)**0.5 + # primary potential in a whole space + prim = rho*1./(4*np.pi*R) + + if flag =="prim": + return prim + + sphind = r < radius + out = np.zeros_like(r) + for n in range(order): + An, Bn = AnBnfun(n, radius, x0, rho, rho1) + dumout = An*r[~sphind]**(-n-1.)*Pleg[n](costheta[~sphind]) + out[~sphind] += dumout + dumin = Bn*r[sphind]**(n)*Pleg[n](costheta[sphind]) + out[sphind] += dumin + + out[~sphind] += prim[~sphind] + + if flag == "sec": + return out-prim + elif flag == "total": + return out + +def AnBnfun(n, radius, x0, rho, rho1, I=1.): + const = I*rho/(4*np.pi) + bunmo = n*rho + (n+1)*rho1 + An = const * radius**(2*n+1) / x0 ** (n+1.) * n * \ + (rho1-rho) / bunmo + Bn = const * 1. / x0 ** (n+1.) * (2*n+1) * (rho1) / bunmo + return An, Bn diff --git a/SimPEG/EM/Analytics/__init__.py b/SimPEG/EM/Analytics/__init__.py index d251f205..9df2aef7 100644 --- a/SimPEG/EM/Analytics/__init__.py +++ b/SimPEG/EM/Analytics/__init__.py @@ -1,4 +1,4 @@ from TDEM import hzAnalyticDipoleT from FDEM import hzAnalyticDipoleF from FDEMcasing import * -from DC import DCAnalytic +from DC import DCAnalyticHalf, DCAnalyticSphere diff --git a/SimPEG/EM/Static/DC/SrcDC.py b/SimPEG/EM/Static/DC/SrcDC.py index 1e64835e..4eda7fe2 100644 --- a/SimPEG/EM/Static/DC/SrcDC.py +++ b/SimPEG/EM/Static/DC/SrcDC.py @@ -1,6 +1,6 @@ import SimPEG # from SimPEG.EM.Base import BaseEMSurvey -from SimPEG.Utils import Zero, closestPoints +from SimPEG.Utils import Zero, closestPoints, mkvc import numpy as np class BaseSrc(SimPEG.Survey.BaseSrc): @@ -27,13 +27,13 @@ class Dipole(BaseSrc): def eval(self, prob): if prob._formulation == 'HJ': - inds = closestPoints(prob.mesh, self.loc) + inds = closestPoints(prob.mesh, self.loc, gridLoc='CC') q = np.zeros(prob.mesh.nC) q[inds] = self.current * np.r_[1., -1.] elif prob._formulation == 'EB': - inds = closestPoints(prob.mesh, self.loc) - q = np.zeros(prob.mesh.nN) - q[inds] = self.current * np.r_[1., -1.] + qa = prob.mesh.getInterpolationMat(self.loc[0], locType='N').todense() + qb = -prob.mesh.getInterpolationMat(self.loc[1], locType='N').todense() + q = mkvc(qa+qb) return q # def bc_contribution diff --git a/SimPEG/Mesh/View.py b/SimPEG/Mesh/View.py index 089d7d9a..6f009dc9 100644 --- a/SimPEG/Mesh/View.py +++ b/SimPEG/Mesh/View.py @@ -206,7 +206,7 @@ class TensorView(object): return out viewOpts = ['real','imag','abs','vec'] normalOpts = ['X', 'Y', 'Z'] - vTypeOpts = ['CC', 'CCv','F','E','Fx','Fy','Fz','E','Ex','Ey','Ez'] + vTypeOpts = ['CC', 'CCv','N','F','E','Fx','Fy','Fz','E','Ex','Ey','Ez'] # Some user error checking assert vType in vTypeOpts, "vType must be in ['%s']" % "','".join(vTypeOpts) diff --git a/tests/em/static/test_DC_analytic.py b/tests/em/static/test_DC_analytic.py new file mode 100644 index 00000000..3755c6ba --- /dev/null +++ b/tests/em/static/test_DC_analytic.py @@ -0,0 +1,68 @@ +import unittest +from SimPEG import Mesh, Utils, EM, Maps, np +import SimPEG.EM.Static.DC as DC + +class DCProblemAnalyticTests(unittest.TestCase): + + def setUp(self): + + cs = 25. + hx = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)] + hy = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)] + hz = [(cs,7, -1.3),(cs,20)] + mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCN") + sigma = np.ones(mesh.nC)*1e-2 + + x = mesh.vectorCCx[(mesh.vectorCCx>-155.)&(mesh.vectorCCx<155.)] + y = mesh.vectorCCx[(mesh.vectorCCy>-155.)&(mesh.vectorCCy<155.)] + Aloc = np.r_[-200., 0., 0.] + Bloc = np.r_[200., 0., 0.] + M = Utils.ndgrid(x-25.,y, np.r_[0.]) + N = Utils.ndgrid(x+25.,y, np.r_[0.]) + phiA = EM.Analytics.DCAnalyticHalf(Aloc, [M,N], 1e-2, flag="halfspace") + phiB = EM.Analytics.DCAnalyticHalf(Bloc, [M,N], 1e-2, flag="halfspace") + data_anal = phiA-phiB + + rx = DC.Rx.Dipole(M, N) + src = DC.Src.Dipole([rx], Aloc, Bloc) + survey = DC.Survey([src]) + + self.survey = survey + self.mesh = mesh + self.sigma = sigma + self.data_anal = data_anal + + try: + from pymatsolver import MumpsSolver + self.Solver = MumpsSolver + except ImportError, e: + self.Solver = SolverLU + + def test_N(self): + problem = DC.Problem3D_N(self.mesh) + problem.Solver = self.Solver + problem.pair(self.survey) + data = self.survey.dpred(self.sigma) + err= np.linalg.norm(data-self.data_anal)/np.linalg.norm(self.data_anal) + if err < 0.2: + passed = True + else: + passed = False + self.assertTrue(passed) + + def test_CC(self): + problem = DC.Problem3D_N(self.mesh) + problem.Solver = self.Solver + problem.pair(self.survey) + data = self.survey.dpred(self.sigma) + err= np.linalg.norm(data-self.data_anal)/np.linalg.norm(self.data_anal) + if err < 0.2: + passed = True + print ">> DC analytic test for Problem3D_CC is pased" + else: + passed = False + self.assertTrue(passed) + +if __name__ == '__main__': + unittest.main() + diff --git a/tests/em/static/test_DC_deriv.py b/tests/em/static/test_DC_deriv.py deleted file mode 100644 index e69de29b..00000000 diff --git a/tests/em/static/test_DC.py b/tests/em/static/test_DC_jvecjtvecadj.py similarity index 100% rename from tests/em/static/test_DC.py rename to tests/em/static/test_DC_jvecjtvecadj.py