From 28d67e311205ed1ba12b24f06a891df8180eb514 Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Fri, 27 May 2016 11:26:47 -0700 Subject: [PATCH 01/19] Start of ED !! --- SimPEG/EM/Analytics/FDEM_fields.py | 51 ++++++++++++++++++++++++++++++ 1 file changed, 51 insertions(+) create mode 100644 SimPEG/EM/Analytics/FDEM_fields.py diff --git a/SimPEG/EM/Analytics/FDEM_fields.py b/SimPEG/EM/Analytics/FDEM_fields.py new file mode 100644 index 00000000..adf10e4e --- /dev/null +++ b/SimPEG/EM/Analytics/FDEM_fields.py @@ -0,0 +1,51 @@ +from __future__ import division +import numpy as np +from scipy.constants import mu_0, pi +from scipy.special import erf +from SimPEG import Utils + + +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 + +def A_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): +def J_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): +def H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): +def B_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): + mu*H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0) From 1960b52dfd58a9fc81c225b5ad005d98ee9ca74a Mon Sep 17 00:00:00 2001 From: micmitch Date: Fri, 27 May 2016 14:56:48 -0700 Subject: [PATCH 02/19] First stab at analytic functions for the fields from a harmonic electric dipole source. Not sure about the exception that I try to throw if multiple frequencies and multiple evaluation locations are both specified. --- SimPEG/EM/Analytics/FDEM_fields.py | 119 +++++++++++++++++++++++++---- 1 file changed, 106 insertions(+), 13 deletions(-) diff --git a/SimPEG/EM/Analytics/FDEM_fields.py b/SimPEG/EM/Analytics/FDEM_fields.py index adf10e4e..abf956e8 100644 --- a/SimPEG/EM/Analytics/FDEM_fields.py +++ b/SimPEG/EM/Analytics/FDEM_fields.py @@ -5,24 +5,27 @@ from scipy.special import erf from SimPEG import Utils -def ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): +def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): + epsilon = 8.854187817*(10.**-12) + omega = 2.*np.pi*f + XYZ = Utils.asArray_N_x_Dim(XYZ, 3) + # Check + if XYZ.shape[0] > 1 & f.length > 1 + except: + print "I/O type error: For multiple field locations only a single frequency can be specified." 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 + # k = np.sqrt( -1j*2.*np.pi*f*mu*sig ) + k = np.sqrt( omega**2. *mu*epsilon -1j*omega*mu*sig ) - front = current * length / (4. * np.pi * sig * r**3) * np.exp(-1j*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 @@ -42,10 +45,100 @@ def ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orienta Ex = front*(dz*dx / r**2)*mid Ey = front*(dz*dy / r**2)*mid return Ex, Ey, Ez - # return Ey, Ez, Ex + +def J_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): + Ex, Ey, Ez = E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0) + Jx = sig*Ex + Jy = sig*Ey + Jz = sig*Ez + return Jx, Jy, Jz + + +def H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): + epsilon = 8.854187817*(10.**-12) + omega = 2.*np.pi*f + + XYZ = Utils.asArray_N_x_Dim(XYZ, 3) + # Check + if XYZ.shape[0] > 1 & f.length > 1 + except: + print "I/O type error: For multiple field locations only a single frequency can be specified." + + 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 ) + k = np.sqrt( omega**2. *mu*epsilon -1j*omega*mu*sig ) + + front = current * length / (4.*np.pi* r**2) * (-1j*k*r + 1) * np.exp(-1j*k*r) + + if orientation.upper() == 'X': + Hy = front*(-dz / r) + Hz = front*(dy / r) + Hx = np.zeros_like(Hy) + return Hx, Hy, Hz + + elif orientation.upper() == 'Y': + Hx = front*(dz / r) + Hz = front*(-dx / r) + Hy = np.zeros_like(Hx) + return Hx, Hy, Hz + + elif orientation.upper() == 'Z': + Hx = front*(-dy / r) + Hy = front*(dx / r) + Hz = np.zeros_like(Hx) + return Hx, Hy, Hz + + +def B_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): + Hx, Hy, Hz = H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0) + Bx = mu*Hx + By = mu*Hy + Bz = mu*Hz + return Bx, By, Bz + def A_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): -def J_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): -def H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): -def B_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): - mu*H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0) + epsilon = 8.854187817*(10.**-12) + omega = 2.*np.pi*f + + XYZ = Utils.asArray_N_x_Dim(XYZ, 3) + # Check + if XYZ.shape[0] > 1 & f.length > 1 + except: + print "I/O type error: For multiple field locations only a single frequency can be specified." + + 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( omega**2. *mu*epsilon -1j*omega*mu*sig ) + + front = current * length / (4.*np.pi*r) + + if orientation.upper() == 'X': + Ax = front*np.exp(-1j*k*r) + Ay = np.zeros_like(Ax) + Az = np.zeros_like(Ax) + return Ax, Ay, Az + + elif orientation.upper() == 'Y': + Ay = front*np.exp(-1j*k*r) + Ax = np.zeros_like(Ay) + Az = np.zeros_like(Ay) + return Ax, Ay, Az + + elif orientation.upper() == 'Z': + Az = front*np.exp(-1j*k*r) + Ax = np.zeros_like(Ay) + Ay = np.zeros_like(Ay) + return Ax, Ay, Az + + + + + From ef12a3674a3394c9c9f05adbf2f1c2f374bfc99f Mon Sep 17 00:00:00 2001 From: D Fournier Date: Mon, 6 Jun 2016 12:28:01 -0700 Subject: [PATCH 03/19] Automate the epsilon picking based on percentile of model values for DEFAULT mode. Fix example. Fix bug with Maps using array of values --- SimPEG/Directives.py | 39 ++++++++++++++++++++----------- SimPEG/Examples/Inversion_IRLS.py | 36 ++++++---------------------- SimPEG/Maps.py | 4 +++- 3 files changed, 36 insertions(+), 43 deletions(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index 932b8fe2..a1e4f8c1 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -253,8 +253,7 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration): class Update_IRLS(InversionDirective): eps_min = None - eps_p = None - eps_q = None + eps = None norms = [2.,2.,2.,2.] factor = None gamma = None @@ -263,6 +262,7 @@ class Update_IRLS(InversionDirective): f_old = None f_min_change = 1e-2 beta_tol = 5e-2 + prctile = 95 # Solving parameter for IRLS (mode:2) IRLSiter = 0 @@ -297,9 +297,22 @@ class Update_IRLS(InversionDirective): print "Convergence with smooth l2-norm regularization: Start IRLS steps..." self.mode = 2 - print self.eps_p, self.eps_q, self.norms - self.reg.eps_p = self.eps_p - self.reg.eps_q = self.eps_q + + # Either use the supplied epsilon, or fix base on distribution of + # model values + if getattr(self, 'reg.eps', None) is None: + self.reg.eps_p = np.percentile(np.abs(self.invProb.curModel),self.prctile) + else: + self.reg.eps_p = self.eps[0] + + if getattr(self, 'reg.eps', None) is None: + self.reg.eps_q = np.percentile(np.abs(self.reg.regmesh.cellDiffxStencil*(self.reg.mapping * self.invProb.curModel)),self.prctile) + else: + self.reg.eps_q = self.eps[1] + + print "L[p qx qy qz]-norm : " + str(self.reg.norms) + print "eps_p: " + str(self.reg.eps_p) + " eps_q: " + str(self.reg.eps_q) + self.reg.norms = self.norms self.coolingFactor = 1. self.coolingRate = 1 @@ -343,14 +356,14 @@ class Update_IRLS(InversionDirective): else: self.f_old = phim_new - # Cool the threshold parameter if required - if getattr(self, 'factor', None) is not None: - eps = self.reg.eps / self.factor - - if getattr(self, 'eps_min', None) is not None: - self.reg.eps = np.max([self.eps_min,eps]) - else: - self.reg.eps = eps +# # Cool the threshold parameter if required +# if getattr(self, 'factor', None) is not None: +# eps = self.reg.eps / self.factor +# +# if getattr(self, 'eps_min', None) is not None: +# self.reg.eps = np.max([self.eps_min,eps]) +# else: +# self.reg.eps = eps # Get phi_m at the end of current iteration self.phi_m_last = self.invProb.phi_m_last diff --git a/SimPEG/Examples/Inversion_IRLS.py b/SimPEG/Examples/Inversion_IRLS.py index afd90525..0a0757cf 100644 --- a/SimPEG/Examples/Inversion_IRLS.py +++ b/SimPEG/Examples/Inversion_IRLS.py @@ -42,55 +42,33 @@ def run(N=100, plotIt=True): survey = Survey.LinearSurvey() survey.pair(prob) survey.dobs = prob.fields(mtrue) + std_noise * np.random.randn(nk) - #survey.makeSyntheticData(mtrue, std=std_noise) wd = np.ones(nk) * std_noise - #print survey.std[0] - #M = prob.mesh # Distance weighting wr = np.sum(prob.G**2.,axis=0)**0.5 wr = ( wr/np.max(wr) ) -# reg = Regularization.Simple(mesh) -# reg.mref = mref -# reg.cell_weights = wr -# dmis = DataMisfit.l2_DataMisfit(survey) dmis.Wd = 1./wd -# -# opt = Optimization.ProjectedGNCG(maxIter=20,lower=-2.,upper=2., maxIterCG= 10, tolCG = 1e-4) -# invProb = InvProblem.BaseInvProblem(dmis, reg, opt) -# invProb.curModel = m0 -# -# beta = Directives.BetaSchedule(coolingFactor=2, coolingRate=1) -# target = Directives.TargetMisfit() -# + betaest = Directives.BetaEstimate_ByEig() -# inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest, target]) -# -# -# mrec = inv.run(m0) -# ml2 = mrec -# print "Final misfit:" + str(invProb.dmisfit.eval(mrec)) -# -# # Switch regularization to sparse -# phim = invProb.phi_m_last -# phid = invProb.phi_d reg = Regularization.Sparse(mesh) reg.mref = mref reg.cell_weights = wr reg.mref = np.zeros(mesh.nC) - eps_p = 5e-2 - eps_q = 5e-2 - norms = [0., 0., 2., 2.] + opt = Optimization.ProjectedGNCG(maxIter=100 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 10, tolCG = 1e-3) invProb = InvProblem.BaseInvProblem(dmis, reg, opt) update_Jacobi = Directives.Update_lin_PreCond() - IRLS = Directives.Update_IRLS( norms=norms, eps_p=eps_p, eps_q=eps_q) + + # Set the IRLS directive, penalize the lowest 25 percentile of model values + # Start with an l2-l2, then switch to lp-norms + norms = [0., 0., 2., 2.] + IRLS = Directives.Update_IRLS( norms=norms, prctile = 25, maxIRLSiter = 15, minGNiter=3) inv = Inversion.BaseInversion(invProb, directiveList=[IRLS,betaest,update_Jacobi]) diff --git a/SimPEG/Maps.py b/SimPEG/Maps.py index 3e97499f..982ed0ef 100644 --- a/SimPEG/Maps.py +++ b/SimPEG/Maps.py @@ -502,7 +502,9 @@ class InjectActiveCells(IdentityMap): if Utils.isScalar(valInactive): self.valInactive = np.ones(self.nC)*float(valInactive) else: - self.valInactive = valInactive.copy() + self.valInactive = np.ones(self.nC) + self.valInactive[self.indInactive] = valInactive.copy() + self.valInactive[self.indActive] = 0 inds = np.nonzero(self.indActive)[0] From a3a5c860088bfb58d388c47b6d2d26c933573e49 Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Sat, 11 Jun 2016 21:58:03 +0200 Subject: [PATCH 04/19] Fix couple bugs in FDEM analytics --- SimPEG/EM/Analytics/FDEM_fields.py | 21 +++++++++------------ SimPEG/EM/Analytics/__init__.py | 1 + 2 files changed, 10 insertions(+), 12 deletions(-) diff --git a/SimPEG/EM/Analytics/FDEM_fields.py b/SimPEG/EM/Analytics/FDEM_fields.py index abf956e8..e0ac1b34 100644 --- a/SimPEG/EM/Analytics/FDEM_fields.py +++ b/SimPEG/EM/Analytics/FDEM_fields.py @@ -11,9 +11,8 @@ def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., XYZ = Utils.asArray_N_x_Dim(XYZ, 3) # Check - if XYZ.shape[0] > 1 & f.length > 1 - except: - print "I/O type error: For multiple field locations only a single frequency can be specified." + if XYZ.shape[0] > 1 & f.shape[0] > 1: + raise Exception("I/O type error: For multiple field locations only a single frequency can be specified.") dx = XYZ[:,0]-srcLoc[0] dy = XYZ[:,1]-srcLoc[1] @@ -51,7 +50,7 @@ def J_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., Jx = sig*Ex Jy = sig*Ey Jz = sig*Ez - return Jx, Jy, Jz + return Jx, Jy, Jz def H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): @@ -60,10 +59,9 @@ def H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., XYZ = Utils.asArray_N_x_Dim(XYZ, 3) # Check - if XYZ.shape[0] > 1 & f.length > 1 - except: - print "I/O type error: For multiple field locations only a single frequency can be specified." - + if XYZ.shape[0] > 1 & f.shape[0] > 1: + raise Exception("I/O type error: For multiple field locations only a single frequency can be specified.") + dx = XYZ[:,0]-srcLoc[0] dy = XYZ[:,1]-srcLoc[1] dz = XYZ[:,2]-srcLoc[2] @@ -99,7 +97,7 @@ def B_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., By = mu*Hy Bz = mu*Hz return Bx, By, Bz - + def A_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): epsilon = 8.854187817*(10.**-12) @@ -107,9 +105,8 @@ def A_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., XYZ = Utils.asArray_N_x_Dim(XYZ, 3) # Check - if XYZ.shape[0] > 1 & f.length > 1 - except: - print "I/O type error: For multiple field locations only a single frequency can be specified." + if XYZ.shape[0] > 1 & f.shape[0] > 1: + raise Exception("I/O type error: For multiple field locations only a single frequency can be specified.") dx = XYZ[:,0]-srcLoc[0] dy = XYZ[:,1]-srcLoc[1] diff --git a/SimPEG/EM/Analytics/__init__.py b/SimPEG/EM/Analytics/__init__.py index 9df2aef7..9c5e5e28 100644 --- a/SimPEG/EM/Analytics/__init__.py +++ b/SimPEG/EM/Analytics/__init__.py @@ -2,3 +2,4 @@ from TDEM import hzAnalyticDipoleT from FDEM import hzAnalyticDipoleF from FDEMcasing import * from DC import DCAnalyticHalf, DCAnalyticSphere +from FDEM_fields import * From 5e3c1da8e2626db11b14578b489faa17e0bad90a Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Mon, 13 Jun 2016 16:41:59 -0700 Subject: [PATCH 05/19] add place holder for galvanic and inductive electric fields... --- SimPEG/EM/Analytics/FDEM_fields.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/SimPEG/EM/Analytics/FDEM_fields.py b/SimPEG/EM/Analytics/FDEM_fields.py index e0ac1b34..4a521ba1 100644 --- a/SimPEG/EM/Analytics/FDEM_fields.py +++ b/SimPEG/EM/Analytics/FDEM_fields.py @@ -5,6 +5,9 @@ from scipy.special import erf from SimPEG import Utils +# def E_galvanic_from_ElectricDipoleWholeSpaced +# def E_inductive_from_ElectricDipoleWholeSpaced + def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): epsilon = 8.854187817*(10.**-12) omega = 2.*np.pi*f @@ -45,6 +48,10 @@ def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., Ey = front*(dz*dy / r**2)*mid return Ex, Ey, Ez + +# def J_galvanic_from_ElectricDipoleWholeSpaced +# def J_inductive_from_ElectricDipoleWholeSpaced + def J_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): Ex, Ey, Ez = E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0) Jx = sig*Ex From 093f441331c034e8dfb27ce91d4495c13e690efb Mon Sep 17 00:00:00 2001 From: micmitch Date: Mon, 13 Jun 2016 17:35:52 -0700 Subject: [PATCH 06/19] Added functions to split electric field into "galvanic" and "inductive" portions. --- SimPEG/EM/Analytics/FDEM_fields.py | 102 +++++++++++++++++++++++++++-- 1 file changed, 96 insertions(+), 6 deletions(-) diff --git a/SimPEG/EM/Analytics/FDEM_fields.py b/SimPEG/EM/Analytics/FDEM_fields.py index 4a521ba1..ac18da7f 100644 --- a/SimPEG/EM/Analytics/FDEM_fields.py +++ b/SimPEG/EM/Analytics/FDEM_fields.py @@ -4,10 +4,6 @@ from scipy.constants import mu_0, pi from scipy.special import erf from SimPEG import Utils - -# def E_galvanic_from_ElectricDipoleWholeSpaced -# def E_inductive_from_ElectricDipoleWholeSpaced - def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): epsilon = 8.854187817*(10.**-12) omega = 2.*np.pi*f @@ -49,8 +45,86 @@ def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., return Ex, Ey, Ez -# def J_galvanic_from_ElectricDipoleWholeSpaced -# def J_inductive_from_ElectricDipoleWholeSpaced +def E_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): + epsilon = 8.854187817*(10.**-12) + omega = 2.*np.pi*f + + XYZ = Utils.asArray_N_x_Dim(XYZ, 3) + # Check + if XYZ.shape[0] > 1 & f.shape[0] > 1: + raise Exception("I/O type error: For multiple field locations only a single frequency can be specified.") + + 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 ) + k = np.sqrt( omega**2. *mu*epsilon -1j*omega*mu*sig ) + + front = current * length / (4.*np.pi*sig* r**3) * np.exp(-1j*k*r) + mid = -k**2 * r**2 + 3*1j*k*r + 3 + + if orientation.upper() == 'X': + Ex_galvanic = front*((dx**2 / r**2)*mid + (-1j*k*r-1.)) + Ey_galvanic = front*(dx*dy / r**2)*mid + Ez_galvanic = front*(dx*dz / r**2)*mid + return Ex_galvanic, Ey_galvanic, Ez_galvanic + + elif orientation.upper() == 'Y': + # x--> y, y--> z, z-->x + Ey_galvanic = front*((dy**2 / r**2)*mid + (-1j*k*r-1.)) + Ez_galvanic = front*(dy*dz / r**2)*mid + Ex_galvanic = front*(dy*dx / r**2)*mid + return Ex_galvanic, Ey_galvanic, Ez_galvanic + + elif orientation.upper() == 'Z': + # x --> z, y --> x, z --> y + Ez_galvanic = front*((dz**2 / r**2)*mid + (-1j*k*r-1.)) + Ex_galvanic = front*(dz*dx / r**2)*mid + Ey_galvanic = front*(dz*dy / r**2)*mid + return Ex_galvanic, Ey_galvanic, Ez_galvanic + + +def E_inductive_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): + epsilon = 8.854187817*(10.**-12) + omega = 2.*np.pi*f + + XYZ = Utils.asArray_N_x_Dim(XYZ, 3) + # Check + if XYZ.shape[0] > 1 & f.shape[0] > 1: + raise Exception("I/O type error: For multiple field locations only a single frequency can be specified.") + + 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 ) + k = np.sqrt( omega**2. *mu*epsilon -1j*omega*mu*sig ) + + front = current * length / (4.*np.pi*sig* r**3) * np.exp(-1j*k*r) + + if orientation.upper() == 'X': + Ex_inductive = front*(k**2 * r**2) + Ey_inductive = 0 + Ez_inductive = 0 + return Ex_inductive, Ey_inductive, Ez_inductive + + elif orientation.upper() == 'Y': + # x--> y, y--> z, z-->x + Ey_inductive = front*(k**2 * r**2) + Ez_inductive = 0 + Ex_inductive = 0 + return Ex_inductive, Ey_inductive, Ez_inductive + + elif orientation.upper() == 'Z': + # x --> z, y --> x, z --> y + Ez_inductive = front*(k**2 * r**2) + Ex_inductive = 0 + Ey_inductive = 0 + return Ex_inductive, Ey_inductive, Ez_inductive + def J_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): Ex, Ey, Ez = E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0) @@ -60,6 +134,22 @@ def J_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., return Jx, Jy, Jz +def J_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): + Ex_galvanic, Ey_galvanic, Ez_galvanic = E_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0) + Jx_galvanic = sig*Ex_galvanic + Jy_galvanic = sig*Ey_galvanic + Jz_galvanic = sig*Ez_galvanic + return Jx_galvanic, Jy_galvanic, Jz_galvanic + + +def J_inductive_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): + Ex_inductive, Ey_inductive, Ez_inductive = E_inductive_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0) + Jx_inductive = sig*Ex_inductive + Jy_inductive = sig*Ey_inductive + Jz_inductive = sig*Ez_inductive + return Jx_inductive, Jy_inductive, Jz_inductive + + def H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): epsilon = 8.854187817*(10.**-12) omega = 2.*np.pi*f From b9d30af4a88903da8941e4afa6355cece6c939eb Mon Sep 17 00:00:00 2001 From: micmitch Date: Mon, 13 Jun 2016 17:47:46 -0700 Subject: [PATCH 07/19] Changed \sigma to \hat{\sigma} = \sigma + i \omega \epsilon in the E field calculations. --- SimPEG/EM/Analytics/FDEM_fields.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/SimPEG/EM/Analytics/FDEM_fields.py b/SimPEG/EM/Analytics/FDEM_fields.py index ac18da7f..467407af 100644 --- a/SimPEG/EM/Analytics/FDEM_fields.py +++ b/SimPEG/EM/Analytics/FDEM_fields.py @@ -7,6 +7,7 @@ from SimPEG import Utils def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): epsilon = 8.854187817*(10.**-12) omega = 2.*np.pi*f + sig_hat = sig + 1j*omega*epsilon XYZ = Utils.asArray_N_x_Dim(XYZ, 3) # Check @@ -21,7 +22,7 @@ def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., # k = np.sqrt( -1j*2.*np.pi*f*mu*sig ) k = np.sqrt( omega**2. *mu*epsilon -1j*omega*mu*sig ) - front = current * length / (4.*np.pi*sig* r**3) * np.exp(-1j*k*r) + front = current * length / (4.*np.pi*sig_hat* r**3) * np.exp(-1j*k*r) mid = -k**2 * r**2 + 3*1j*k*r + 3 if orientation.upper() == 'X': @@ -48,6 +49,7 @@ def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., def E_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): epsilon = 8.854187817*(10.**-12) omega = 2.*np.pi*f + sig_hat = sig + 1j*omega*epsilon XYZ = Utils.asArray_N_x_Dim(XYZ, 3) # Check @@ -62,7 +64,7 @@ def E_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., l # k = np.sqrt( -1j*2.*np.pi*f*mu*sig ) k = np.sqrt( omega**2. *mu*epsilon -1j*omega*mu*sig ) - front = current * length / (4.*np.pi*sig* r**3) * np.exp(-1j*k*r) + front = current * length / (4.*np.pi*sig_hat* r**3) * np.exp(-1j*k*r) mid = -k**2 * r**2 + 3*1j*k*r + 3 if orientation.upper() == 'X': @@ -89,6 +91,7 @@ def E_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., l def E_inductive_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): epsilon = 8.854187817*(10.**-12) omega = 2.*np.pi*f + sig_hat = sig + 1j*omega*epsilon XYZ = Utils.asArray_N_x_Dim(XYZ, 3) # Check @@ -103,7 +106,7 @@ def E_inductive_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., # k = np.sqrt( -1j*2.*np.pi*f*mu*sig ) k = np.sqrt( omega**2. *mu*epsilon -1j*omega*mu*sig ) - front = current * length / (4.*np.pi*sig* r**3) * np.exp(-1j*k*r) + front = current * length / (4.*np.pi*sig_hat* r**3) * np.exp(-1j*k*r) if orientation.upper() == 'X': Ex_inductive = front*(k**2 * r**2) From e815ddaec79d59e142a9a3899d0fc1bac7434340 Mon Sep 17 00:00:00 2001 From: micmitch Date: Tue, 14 Jun 2016 15:15:53 -0700 Subject: [PATCH 08/19] Removed d typos from the end of function names. --- SimPEG/EM/Analytics/FDEM_fields.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/SimPEG/EM/Analytics/FDEM_fields.py b/SimPEG/EM/Analytics/FDEM_fields.py index 467407af..bd1ee368 100644 --- a/SimPEG/EM/Analytics/FDEM_fields.py +++ b/SimPEG/EM/Analytics/FDEM_fields.py @@ -46,7 +46,7 @@ def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., return Ex, Ey, Ez -def E_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): +def E_galvanic_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): epsilon = 8.854187817*(10.**-12) omega = 2.*np.pi*f sig_hat = sig + 1j*omega*epsilon @@ -88,7 +88,7 @@ def E_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., l return Ex_galvanic, Ey_galvanic, Ez_galvanic -def E_inductive_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): +def E_inductive_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): epsilon = 8.854187817*(10.**-12) omega = 2.*np.pi*f sig_hat = sig + 1j*omega*epsilon @@ -137,7 +137,7 @@ def J_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., return Jx, Jy, Jz -def J_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): +def J_galvanic_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): Ex_galvanic, Ey_galvanic, Ez_galvanic = E_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0) Jx_galvanic = sig*Ex_galvanic Jy_galvanic = sig*Ey_galvanic @@ -145,7 +145,7 @@ def J_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., l return Jx_galvanic, Jy_galvanic, Jz_galvanic -def J_inductive_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): +def J_inductive_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): Ex_inductive, Ey_inductive, Ez_inductive = E_inductive_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0) Jx_inductive = sig*Ex_inductive Jy_inductive = sig*Ey_inductive From f0944362c8afded6dc263629681775705fb73256 Mon Sep 17 00:00:00 2001 From: micmitch Date: Tue, 14 Jun 2016 15:22:36 -0700 Subject: [PATCH 09/19] Silly mistake... needed zero arrays instead of scalars. --- SimPEG/EM/Analytics/FDEM_fields.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/SimPEG/EM/Analytics/FDEM_fields.py b/SimPEG/EM/Analytics/FDEM_fields.py index bd1ee368..ca250ad6 100644 --- a/SimPEG/EM/Analytics/FDEM_fields.py +++ b/SimPEG/EM/Analytics/FDEM_fields.py @@ -110,22 +110,22 @@ def E_inductive_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., l if orientation.upper() == 'X': Ex_inductive = front*(k**2 * r**2) - Ey_inductive = 0 - Ez_inductive = 0 + Ey_inductive = np.zeros_like(Ex_inductive) + Ez_inductive = np.zeros_like(Ex_inductive) return Ex_inductive, Ey_inductive, Ez_inductive elif orientation.upper() == 'Y': # x--> y, y--> z, z-->x Ey_inductive = front*(k**2 * r**2) - Ez_inductive = 0 - Ex_inductive = 0 + Ez_inductive = np.zeros_like(Ey_inductive) + Ex_inductive = np.zeros_like(Ey_inductive) return Ex_inductive, Ey_inductive, Ez_inductive elif orientation.upper() == 'Z': # x --> z, y --> x, z --> y Ez_inductive = front*(k**2 * r**2) - Ex_inductive = 0 - Ey_inductive = 0 + Ex_inductive = np.zeros_like(Ez_inductive) + Ey_inductive = np.zeros_like(Ez_inductive) return Ex_inductive, Ey_inductive, Ez_inductive From a54713f54688c268ecb7d5e3ed261c3ff2ae3881 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 22 Jun 2016 11:32:54 -0600 Subject: [PATCH 10/19] Change to NotImplementedError. --- SimPEG/Survey.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/SimPEG/Survey.py b/SimPEG/Survey.py index fbc88276..65c46972 100644 --- a/SimPEG/Survey.py +++ b/SimPEG/Survey.py @@ -311,7 +311,6 @@ class BaseSurvey(object): if f is None: f = self.prob.fields(m) return Utils.mkvc(self.eval(f)) - @Utils.count def eval(self, f): """eval(f) @@ -322,7 +321,7 @@ class BaseSurvey(object): d_\\text{pred} = \mathbf{P} f(m) """ - raise NotImplemented('eval is not yet implemented.') + raise NotImplementedError('eval is not yet implemented.') @Utils.count def evalDeriv(self, f): @@ -334,7 +333,7 @@ class BaseSurvey(object): \\frac{\partial d_\\text{pred}}{\partial u} = \mathbf{P} """ - raise NotImplemented('eval is not yet implemented.') + raise NotImplementedError('eval is not yet implemented.') @Utils.count def residual(self, m, f=None): From e1ba80883d1dcaaac8f7804812ec216b8b974fea Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Thu, 23 Jun 2016 09:10:50 -0700 Subject: [PATCH 11/19] Incorporate Lindsey's suggestoins --- .../{FDEM_fields.py => FDEMDipolarfields.py} | 125 +++++++++++++----- 1 file changed, 93 insertions(+), 32 deletions(-) rename SimPEG/EM/Analytics/{FDEM_fields.py => FDEMDipolarfields.py} (69%) diff --git a/SimPEG/EM/Analytics/FDEM_fields.py b/SimPEG/EM/Analytics/FDEMDipolarfields.py similarity index 69% rename from SimPEG/EM/Analytics/FDEM_fields.py rename to SimPEG/EM/Analytics/FDEMDipolarfields.py index ca250ad6..ef724542 100644 --- a/SimPEG/EM/Analytics/FDEM_fields.py +++ b/SimPEG/EM/Analytics/FDEMDipolarfields.py @@ -1,13 +1,24 @@ from __future__ import division import numpy as np -from scipy.constants import mu_0, pi +from scipy.constants import mu_0, pi, epsilon_0 from scipy.special import erf from SimPEG import Utils -def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): - epsilon = 8.854187817*(10.**-12) - omega = 2.*np.pi*f - sig_hat = sig + 1j*omega*epsilon +omega = lambda f: 2.*np.pi*f +# TODO: +# r = lambda dx, dy, dz: np.sqrt( dx**2. + dy**2. + dz**2.) +# k = lambda f, mu, epsilon, sig: np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig ) + +def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=0., epsr=1.): + + """ + Computing Analytic Electric fields from Electrical Dipole in a Wholespace + TODO: + Add description of parameters + """ + mu = mu_0*(1+kappa) + epsilon = epsilon_0*epsr + sig_hat = sig + 1j*omeg*epsilon XYZ = Utils.asArray_N_x_Dim(XYZ, 3) # Check @@ -20,7 +31,7 @@ def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., r = np.sqrt( dx**2. + dy**2. + dz**2.) # k = np.sqrt( -1j*2.*np.pi*f*mu*sig ) - k = np.sqrt( omega**2. *mu*epsilon -1j*omega*mu*sig ) + k = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig ) front = current * length / (4.*np.pi*sig_hat* r**3) * np.exp(-1j*k*r) mid = -k**2 * r**2 + 3*1j*k*r + 3 @@ -46,10 +57,16 @@ def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., return Ex, Ey, Ez -def E_galvanic_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): - epsilon = 8.854187817*(10.**-12) - omega = 2.*np.pi*f - sig_hat = sig + 1j*omega*epsilon +def E_galvanic_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.): + + """ + Computing Galvanic portion of Electric fields from Electrical Dipole in a Wholespace + TODO: + Add description of parameters + """ + mu = mu_0*(1+kappa) + epsilon = epsilon_0*epsr + sig_hat = sig + 1j*omeg*epsilon XYZ = Utils.asArray_N_x_Dim(XYZ, 3) # Check @@ -62,7 +79,7 @@ def E_galvanic_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., le r = np.sqrt( dx**2. + dy**2. + dz**2.) # k = np.sqrt( -1j*2.*np.pi*f*mu*sig ) - k = np.sqrt( omega**2. *mu*epsilon -1j*omega*mu*sig ) + k = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig ) front = current * length / (4.*np.pi*sig_hat* r**3) * np.exp(-1j*k*r) mid = -k**2 * r**2 + 3*1j*k*r + 3 @@ -88,10 +105,16 @@ def E_galvanic_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., le return Ex_galvanic, Ey_galvanic, Ez_galvanic -def E_inductive_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): - epsilon = 8.854187817*(10.**-12) - omega = 2.*np.pi*f - sig_hat = sig + 1j*omega*epsilon +def E_inductive_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.): + + """ + Computing Inductive portion of Electric fields from Electrical Dipole in a Wholespace + TODO: + Add description of parameters + """ + mu = mu_0*(1+kappa) + epsilon = epsilon_0*epsr + sig_hat = sig + 1j*omeg*epsilon XYZ = Utils.asArray_N_x_Dim(XYZ, 3) # Check @@ -104,7 +127,7 @@ def E_inductive_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., l r = np.sqrt( dx**2. + dy**2. + dz**2.) # k = np.sqrt( -1j*2.*np.pi*f*mu*sig ) - k = np.sqrt( omega**2. *mu*epsilon -1j*omega*mu*sig ) + k = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig ) front = current * length / (4.*np.pi*sig_hat* r**3) * np.exp(-1j*k*r) @@ -129,34 +152,60 @@ def E_inductive_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., l return Ex_inductive, Ey_inductive, Ez_inductive -def J_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): - Ex, Ey, Ez = E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0) +def J_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.): + + """ + Computing Current densities from Electrical Dipole in a Wholespace + TODO: + Add description of parameters + """ + + Ex, Ey, Ez = E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.) Jx = sig*Ex Jy = sig*Ey Jz = sig*Ez return Jx, Jy, Jz -def J_galvanic_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): - Ex_galvanic, Ey_galvanic, Ez_galvanic = E_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0) +def J_galvanic_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.): + + """ + Computing Galvanic portion of Current densities from Electrical Dipole in a Wholespace + TODO: + Add description of parameters + """ + + Ex_galvanic, Ey_galvanic, Ez_galvanic = E_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.) Jx_galvanic = sig*Ex_galvanic Jy_galvanic = sig*Ey_galvanic Jz_galvanic = sig*Ez_galvanic return Jx_galvanic, Jy_galvanic, Jz_galvanic -def J_inductive_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): - Ex_inductive, Ey_inductive, Ez_inductive = E_inductive_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0) +def J_inductive_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.): + + """ + Computing Inductive portion of Current densities from Electrical Dipole in a Wholespace + TODO: + Add description of parameters + """ + + Ex_inductive, Ey_inductive, Ez_inductive = E_inductive_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.) Jx_inductive = sig*Ex_inductive Jy_inductive = sig*Ey_inductive Jz_inductive = sig*Ez_inductive return Jx_inductive, Jy_inductive, Jz_inductive -def H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): - epsilon = 8.854187817*(10.**-12) - omega = 2.*np.pi*f +def H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.): + """ + Computing Magnetic fields from Electrical Dipole in a Wholespace + TODO: + Add description of parameters + """ + mu = mu_0*(1+kappa) + epsilon = epsilon_0*epsr XYZ = Utils.asArray_N_x_Dim(XYZ, 3) # Check if XYZ.shape[0] > 1 & f.shape[0] > 1: @@ -168,7 +217,7 @@ def H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., r = np.sqrt( dx**2. + dy**2. + dz**2.) # k = np.sqrt( -1j*2.*np.pi*f*mu*sig ) - k = np.sqrt( omega**2. *mu*epsilon -1j*omega*mu*sig ) + k = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig ) front = current * length / (4.*np.pi* r**2) * (-1j*k*r + 1) * np.exp(-1j*k*r) @@ -191,18 +240,30 @@ def H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., return Hx, Hy, Hz -def B_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): - Hx, Hy, Hz = H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0) +def B_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.): + + """ + Computing Magnetic flux densites from Electrical Dipole in a Wholespace + TODO: + Add description of parameters + """ + + Hx, Hy, Hz = H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.) Bx = mu*Hx By = mu*Hy Bz = mu*Hz return Bx, By, Bz -def A_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0): - epsilon = 8.854187817*(10.**-12) - omega = 2.*np.pi*f +def A_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.): + """ + Computing Electric vector potentials from Electrical Dipole in a Wholespace + TODO: + Add description of parameters + """ + mu = mu_0*(1+kappa) + epsilon = epsilon_0*epsr XYZ = Utils.asArray_N_x_Dim(XYZ, 3) # Check if XYZ.shape[0] > 1 & f.shape[0] > 1: @@ -213,7 +274,7 @@ def A_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., dz = XYZ[:,2]-srcLoc[2] r = np.sqrt( dx**2. + dy**2. + dz**2.) - k = np.sqrt( omega**2. *mu*epsilon -1j*omega*mu*sig ) + k = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig ) front = current * length / (4.*np.pi*r) From 8b44f8d96b91dec91a1f0d5dc1cad1c96e59cf43 Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Thu, 23 Jun 2016 10:14:13 -0700 Subject: [PATCH 12/19] change FDEM_fields.py to FDEMDipolarfield.py --- SimPEG/EM/Analytics/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/EM/Analytics/__init__.py b/SimPEG/EM/Analytics/__init__.py index 9c5e5e28..8331501d 100644 --- a/SimPEG/EM/Analytics/__init__.py +++ b/SimPEG/EM/Analytics/__init__.py @@ -2,4 +2,4 @@ from TDEM import hzAnalyticDipoleT from FDEM import hzAnalyticDipoleF from FDEMcasing import * from DC import DCAnalyticHalf, DCAnalyticSphere -from FDEM_fields import * +from FDEMDipolarfields import * From 0763925743e4ff26ce2e4831824870b78ab5e0c7 Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Thu, 23 Jun 2016 14:12:49 -0700 Subject: [PATCH 13/19] fix minor bugs in analytics --- SimPEG/EM/Analytics/FDEMDipolarfields.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/SimPEG/EM/Analytics/FDEMDipolarfields.py b/SimPEG/EM/Analytics/FDEMDipolarfields.py index ef724542..7b79a680 100644 --- a/SimPEG/EM/Analytics/FDEMDipolarfields.py +++ b/SimPEG/EM/Analytics/FDEMDipolarfields.py @@ -18,7 +18,7 @@ def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., """ mu = mu_0*(1+kappa) epsilon = epsilon_0*epsr - sig_hat = sig + 1j*omeg*epsilon + sig_hat = sig + 1j*omega(f)*epsilon XYZ = Utils.asArray_N_x_Dim(XYZ, 3) # Check @@ -66,7 +66,7 @@ def E_galvanic_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., le """ mu = mu_0*(1+kappa) epsilon = epsilon_0*epsr - sig_hat = sig + 1j*omeg*epsilon + sig_hat = sig + 1j*omega(f)*epsilon XYZ = Utils.asArray_N_x_Dim(XYZ, 3) # Check From 3960cfc3134af4461ab5240054110f93266e72c5 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Sun, 26 Jun 2016 16:17:30 -0600 Subject: [PATCH 14/19] Update index.rst --- docs/index.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/index.rst b/docs/index.rst index fb6fb4f3..fd6959bc 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -1,4 +1,4 @@ -.. image:: https://raw.github.com/simpeg/simpeg/master/docs/simpeg-logo.png +.. image:: https://raw.github.com/simpeg/simpeg/master/docs/images/simpeg-logo.png :alt: SimPEG Logo SimPEG Documentation From 8ed3ec18fab625c1f8a5128a257ae5bd71c54982 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Sun, 26 Jun 2016 16:29:20 -0600 Subject: [PATCH 15/19] Update README.rst --- README.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.rst b/README.rst index c9b30a78..85777c4c 100644 --- a/README.rst +++ b/README.rst @@ -1,4 +1,4 @@ -.. image:: https://raw.github.com/simpeg/simpeg/master/docs/simpeg-logo.png +.. image:: https://raw.github.com/simpeg/simpeg/master/docs/images/simpeg-logo.png :alt: SimPEG Logo ====== From ba173674ecf513742b82a7530e1dcce4982504e2 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Sun, 26 Jun 2016 17:07:07 -0600 Subject: [PATCH 16/19] Mesh2Mesh and Combo Map examples. Also fixed plotting codes to show the plots by default. --- .gitignore | 1 + SimPEG/Examples/DC_Analytic_Dipole.py | 4 +- SimPEG/Examples/Maps_ComboMaps.py | 62 +++++++++++++++++++++++ SimPEG/Examples/Maps_Mesh2Mesh.py | 41 +++++++++++++++ SimPEG/Examples/Utils_surface2ind_topo.py | 2 +- SimPEG/Examples/__init__.py | 4 +- docs/content/api_core/api_Maps.rst | 47 ++--------------- docs/content/examples/Inversion_IRLS.rst | 26 ++++++++++ docs/content/examples/Maps_ComboMaps.rst | 48 ++++++++++++++++++ docs/content/examples/Maps_Mesh2Mesh.rst | 27 ++++++++++ 10 files changed, 215 insertions(+), 47 deletions(-) create mode 100644 SimPEG/Examples/Maps_ComboMaps.py create mode 100644 SimPEG/Examples/Maps_Mesh2Mesh.py create mode 100644 docs/content/examples/Inversion_IRLS.rst create mode 100644 docs/content/examples/Maps_ComboMaps.rst create mode 100644 docs/content/examples/Maps_Mesh2Mesh.rst diff --git a/.gitignore b/.gitignore index e1264b7d..5fd423b1 100644 --- a/.gitignore +++ b/.gitignore @@ -40,3 +40,4 @@ nosetests.xml docs/_build/ Makefile docs/warnings.txt +.DS_Store diff --git a/SimPEG/Examples/DC_Analytic_Dipole.py b/SimPEG/Examples/DC_Analytic_Dipole.py index ba999fa3..347446e3 100644 --- a/SimPEG/Examples/DC_Analytic_Dipole.py +++ b/SimPEG/Examples/DC_Analytic_Dipole.py @@ -1,7 +1,7 @@ from SimPEG import * import SimPEG.EM.Static.DC as DC -def run(plotIt=False): +def run(plotIt=True): cs = 25. hx = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)] hy = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)] @@ -65,4 +65,4 @@ def run(plotIt=False): if __name__ == '__main__': - print run(plotIt=True) + print run() diff --git a/SimPEG/Examples/Maps_ComboMaps.py b/SimPEG/Examples/Maps_ComboMaps.py new file mode 100644 index 00000000..3e495023 --- /dev/null +++ b/SimPEG/Examples/Maps_ComboMaps.py @@ -0,0 +1,62 @@ +from SimPEG import Mesh, Maps, np + +def run(plotIt=True): + """ + + Maps: ComboMaps + =============== + + We will use an example where we want a 1D layered earth as + our model, but we want to map this to a 2D discretization to do our forward + modeling. We will also assume that we are working in log conductivity still, + so after the transformation we want to map to conductivity space. + To do this we will introduce the vertical 1D map (:class:`SimPEG.Maps.SurjectVertical1D`), + which does the first part of what we just described. The second part will be + done by the :class:`SimPEG.Maps.ExpMap` described above. + + .. code-block:: python + :linenos: + + M = Mesh.TensorMesh([7,5]) + v1dMap = Maps.SurjectVertical1D(M) + expMap = Maps.ExpMap(M) + myMap = expMap * v1dMap + m = np.r_[0.2,1,0.1,2,2.9] # only 5 model parameters! + sig = myMap * m + + If you noticed, it was pretty easy to combine maps. What is even cooler is + that the derivatives also are made for you (if everything goes right). + Just to be sure that the derivative is correct, you should always run the test + on the mapping that you create. + + """ + + + M = Mesh.TensorMesh([7,5]) + v1dMap = Maps.SurjectVertical1D(M) + expMap = Maps.ExpMap(M) + myMap = expMap * v1dMap + m = np.r_[0.2,1,0.1,2,2.9] # only 5 model parameters! + sig = myMap * m + + if not plotIt: return + + import matplotlib.pyplot as plt + figs, axs = plt.subplots(1,2) + axs[0].plot(m, M.vectorCCy, 'b-o') + axs[0].set_title('Model') + axs[0].set_ylabel('Depth, y') + axs[0].set_xlabel('Value, $m_i$') + axs[0].set_xlim(0,3) + axs[0].set_ylim(0,1) + clbar = plt.colorbar(M.plotImage(sig,ax=axs[1],grid=True,gridOpts=dict(color='grey'))[0]) + axs[1].set_title('Physical Property') + axs[1].set_ylabel('Depth, y') + clbar.set_label('$\sigma = \exp(\mathbf{P}m)$') + plt.tight_layout() + plt.show() + + +if __name__ == '__main__': + run() + diff --git a/SimPEG/Examples/Maps_Mesh2Mesh.py b/SimPEG/Examples/Maps_Mesh2Mesh.py new file mode 100644 index 00000000..254dc036 --- /dev/null +++ b/SimPEG/Examples/Maps_Mesh2Mesh.py @@ -0,0 +1,41 @@ +from SimPEG import Mesh, Maps, Utils + +def run(plotIt=True): + """ + + Maps: Mesh2Mesh + =============== + + This mapping allows you to go from one mesh to another. + + """ + + M = Mesh.TensorMesh([100,100]) + h1 = Utils.meshTensor([(6,7,-1.5),(6,10),(6,7,1.5)]) + h1 = h1/h1.sum() + M2 = Mesh.TensorMesh([h1,h1]) + V = Utils.ModelBuilder.randomModel(M.vnC, seed=79, its=50) + v = Utils.mkvc(V) + modh = Maps.Mesh2Mesh([M,M2]) + modH = Maps.Mesh2Mesh([M2,M]) + H = modH * v + h = modh * H + + if not plotIt: return + + import matplotlib.pyplot as plt + ax = plt.subplot(131) + M.plotImage(v, ax=ax) + ax.set_title('Fine Mesh (Original)') + ax = plt.subplot(132) + M2.plotImage(H,clim=[0,1],ax=ax) + ax.set_title('Course Mesh') + ax = plt.subplot(133) + M.plotImage(h,clim=[0,1],ax=ax) + ax.set_title('Fine Mesh (Interpolated)') + plt.show() + + +if __name__ == '__main__': + run() + diff --git a/SimPEG/Examples/Utils_surface2ind_topo.py b/SimPEG/Examples/Utils_surface2ind_topo.py index e4f748cb..07e3e4aa 100644 --- a/SimPEG/Examples/Utils_surface2ind_topo.py +++ b/SimPEG/Examples/Utils_surface2ind_topo.py @@ -2,7 +2,7 @@ from SimPEG import * from SimPEG.Utils import surface2ind_topo -def run(plotIt=False, nx=5, ny=5): +def run(plotIt=True, nx=5, ny=5): """ Utils: surface2ind_topo diff --git a/SimPEG/Examples/__init__.py b/SimPEG/Examples/__init__.py index 2ee4eb6d..3c9cbca9 100644 --- a/SimPEG/Examples/__init__.py +++ b/SimPEG/Examples/__init__.py @@ -10,6 +10,8 @@ import EM_TDEM_1D_Inversion import FLOW_Richards_1D_Celia1990 import Inversion_IRLS import Inversion_Linear +import Maps_ComboMaps +import Maps_Mesh2Mesh import Mesh_Basic_ForwardDC import Mesh_Basic_PlotImage import Mesh_Basic_Types @@ -22,7 +24,7 @@ import MT_1D_ForwardAndInversion import MT_3D_Foward import Utils_surface2ind_topo -__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Inversion_IRLS", "Inversion_Linear", "Mesh_Basic_ForwardDC", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward", "Utils_surface2ind_topo"] +__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Inversion_IRLS", "Inversion_Linear", "Maps_ComboMaps", "Maps_Mesh2Mesh", "Mesh_Basic_ForwardDC", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward", "Utils_surface2ind_topo"] ##### AUTOIMPORTS ##### diff --git a/docs/content/api_core/api_Maps.rst b/docs/content/api_core/api_Maps.rst index a5de8a2f..e1e744c6 100644 --- a/docs/content/api_core/api_Maps.rst +++ b/docs/content/api_core/api_Maps.rst @@ -63,26 +63,8 @@ done by the :class:`SimPEG.Maps.ExpMap` described above. .. plot:: - from SimPEG import * - import matplotlib.pyplot as plt - M = Mesh.TensorMesh([7,5]) - v1dMap = Maps.SurjectVertical1D(M) - expMap = Maps.ExpMap(M) - myMap = expMap * v1dMap - m = np.r_[0.2,1,0.1,2,2.9] # only 5 model parameters! - sig = myMap * m - figs, axs = plt.subplots(1,2) - axs[0].plot(m, M.vectorCCy, 'b-o') - axs[0].set_title('Model') - axs[0].set_ylabel('Depth, y') - axs[0].set_xlabel('Value, $m_i$') - axs[0].set_xlim(0,3) - axs[0].set_ylim(0,1) - clbar = plt.colorbar(M.plotImage(sig,ax=axs[1],grid=True,gridOpts=dict(color='grey'))[0]) - axs[1].set_title('Physical Property') - axs[1].set_ylabel('Depth, y') - clbar.set_label('$\sigma = \exp(\mathbf{P}m)$') - plt.tight_layout() + from SimPEG import Examples + Examples.Maps_ComboMaps.run() If you noticed, it was pretty easy to combine maps. What is even cooler is that the derivatives also are made for you (if everything goes right). @@ -167,31 +149,10 @@ Map 2D Cross-Section to 3D Model Mesh to Mesh Map ---------------- - .. plot:: - from SimPEG import * - import matplotlib.pyplot as plt - M = Mesh.TensorMesh([100,100]) - h1 = Utils.meshTensor([(6,7,-1.5),(6,10),(6,7,1.5)]) - h1 = h1/h1.sum() - M2 = Mesh.TensorMesh([h1,h1]) - V = Utils.ModelBuilder.randomModel(M.vnC, seed=79, its=50) - v = Utils.mkvc(V) - modh = Maps.Mesh2Mesh([M,M2]) - modH = Maps.Mesh2Mesh([M2,M]) - H = modH * v - h = modh * H - ax = plt.subplot(131) - M.plotImage(v, ax=ax) - ax.set_title('Fine Mesh (Original)') - ax = plt.subplot(132) - M2.plotImage(H,clim=[0,1],ax=ax) - ax.set_title('Course Mesh') - ax = plt.subplot(133) - M.plotImage(h,clim=[0,1],ax=ax) - ax.set_title('Fine Mesh (Interpolated)') - plt.show() + from SimPEG import Examples + Examples.Maps_Mesh2Mesh.run() .. autoclass:: SimPEG.Maps.Mesh2Mesh diff --git a/docs/content/examples/Inversion_IRLS.rst b/docs/content/examples/Inversion_IRLS.rst new file mode 100644 index 00000000..192ada50 --- /dev/null +++ b/docs/content/examples/Inversion_IRLS.rst @@ -0,0 +1,26 @@ +.. _examples_Inversion_IRLS: + +.. --------------------------------- .. +.. .. +.. THIS FILE IS AUTO GENEREATED .. +.. .. +.. SimPEG/Examples/__init__.py .. +.. .. +.. --------------------------------- .. + + +Inversion: Linear Problem +========================= + +Here we go over the basics of creating a linear problem and inversion. + + + +.. plot:: + + from SimPEG import Examples + Examples.Inversion_IRLS.run() + +.. literalinclude:: ../../../SimPEG/Examples/Inversion_IRLS.py + :language: python + :linenos: diff --git a/docs/content/examples/Maps_ComboMaps.rst b/docs/content/examples/Maps_ComboMaps.rst new file mode 100644 index 00000000..69ebbf39 --- /dev/null +++ b/docs/content/examples/Maps_ComboMaps.rst @@ -0,0 +1,48 @@ +.. _examples_Maps_ComboMaps: + +.. --------------------------------- .. +.. .. +.. THIS FILE IS AUTO GENEREATED .. +.. .. +.. SimPEG/Examples/__init__.py .. +.. .. +.. --------------------------------- .. + + + +Maps: ComboMaps +=============== + +We will use an example where we want a 1D layered earth as +our model, but we want to map this to a 2D discretization to do our forward +modeling. We will also assume that we are working in log conductivity still, +so after the transformation we want to map to conductivity space. +To do this we will introduce the vertical 1D map (:class:`SimPEG.Maps.SurjectVertical1D`), +which does the first part of what we just described. The second part will be +done by the :class:`SimPEG.Maps.ExpMap` described above. + +.. code-block:: python + :linenos: + + M = Mesh.TensorMesh([7,5]) + v1dMap = Maps.SurjectVertical1D(M) + expMap = Maps.ExpMap(M) + myMap = expMap * v1dMap + m = np.r_[0.2,1,0.1,2,2.9] # only 5 model parameters! + sig = myMap * m + +If you noticed, it was pretty easy to combine maps. What is even cooler is +that the derivatives also are made for you (if everything goes right). +Just to be sure that the derivative is correct, you should always run the test +on the mapping that you create. + + + +.. plot:: + + from SimPEG import Examples + Examples.Maps_ComboMaps.run() + +.. literalinclude:: ../../../SimPEG/Examples/Maps_ComboMaps.py + :language: python + :linenos: diff --git a/docs/content/examples/Maps_Mesh2Mesh.rst b/docs/content/examples/Maps_Mesh2Mesh.rst new file mode 100644 index 00000000..60233b57 --- /dev/null +++ b/docs/content/examples/Maps_Mesh2Mesh.rst @@ -0,0 +1,27 @@ +.. _examples_Maps_Mesh2Mesh: + +.. --------------------------------- .. +.. .. +.. THIS FILE IS AUTO GENEREATED .. +.. .. +.. SimPEG/Examples/__init__.py .. +.. .. +.. --------------------------------- .. + + + +Maps: Mesh2Mesh +=============== + +This mapping allows you to go from one mesh to another. + + + +.. plot:: + + from SimPEG import Examples + Examples.Maps_Mesh2Mesh.run() + +.. literalinclude:: ../../../SimPEG/Examples/Maps_Mesh2Mesh.py + :language: python + :linenos: From eda23944113a6525e7a332c9c003f075ececd4cd Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Mon, 27 Jun 2016 13:04:30 -0700 Subject: [PATCH 17/19] fix bug for omega. --- SimPEG/EM/Analytics/FDEMDipolarfields.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/EM/Analytics/FDEMDipolarfields.py b/SimPEG/EM/Analytics/FDEMDipolarfields.py index 7b79a680..b643e649 100644 --- a/SimPEG/EM/Analytics/FDEMDipolarfields.py +++ b/SimPEG/EM/Analytics/FDEMDipolarfields.py @@ -114,7 +114,7 @@ def E_inductive_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., l """ mu = mu_0*(1+kappa) epsilon = epsilon_0*epsr - sig_hat = sig + 1j*omeg*epsilon + sig_hat = sig + 1j*omega(f)*epsilon XYZ = Utils.asArray_N_x_Dim(XYZ, 3) # Check From 334cd8e454bcb8a7b24b0c6dbbcc75f464364baf Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Wed, 29 Jun 2016 09:49:29 -0700 Subject: [PATCH 18/19] =?UTF-8?q?Bump=20version:=200.1.11=20=E2=86=92=200.?= =?UTF-8?q?1.12?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .bumpversion.cfg | 2 +- SimPEG/__init__.py | 2 +- docs/conf.py | 4 ++-- setup.py | 2 +- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/.bumpversion.cfg b/.bumpversion.cfg index 3ca3f780..31842d79 100644 --- a/.bumpversion.cfg +++ b/.bumpversion.cfg @@ -1,4 +1,4 @@ [bumpversion] -current_version = 0.1.11 +current_version = 0.1.12 files = setup.py SimPEG/__init__.py docs/conf.py diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index 23d4de4a..536a8103 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -15,7 +15,7 @@ import Directives import Inversion import Tests -__version__ = '0.1.11' +__version__ = '0.1.12' __author__ = 'Rowan Cockett' __license__ = 'MIT' __copyright__ = 'Copyright 2014 Rowan Cockett' diff --git a/docs/conf.py b/docs/conf.py index 143ea545..d696464d 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -51,9 +51,9 @@ copyright = u'2013 - 2016, SimPEG Developers' # built documents. # # The short X.Y version. -version = '0.1.11' +version = '0.1.12' # The full version, including alpha/beta/rc tags. -release = '0.1.11' +release = '0.1.12' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/setup.py b/setup.py index 06cd9b27..a52bfe2c 100644 --- a/setup.py +++ b/setup.py @@ -83,7 +83,7 @@ with open("README.rst") as f: setup( name = "SimPEG", - version = "0.1.11", + version = "0.1.12", packages = find_packages(), install_requires = ['numpy>=1.7', 'scipy>=0.13', From a289b656cd50d0602bd8ada1e2a004cd3181c5dc Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Wed, 29 Jun 2016 13:09:11 -0700 Subject: [PATCH 19/19] Fixes for kwargs variables in FDEMDipolarfields.py --- SimPEG/EM/Analytics/FDEMDipolarfields.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/SimPEG/EM/Analytics/FDEMDipolarfields.py b/SimPEG/EM/Analytics/FDEMDipolarfields.py index b643e649..e65bbdd1 100644 --- a/SimPEG/EM/Analytics/FDEMDipolarfields.py +++ b/SimPEG/EM/Analytics/FDEMDipolarfields.py @@ -160,7 +160,7 @@ def J_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., Add description of parameters """ - Ex, Ey, Ez = E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.) + Ex, Ey, Ez = E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=current, length=length, orientation=orientation, kappa=kappa, epsr=epsr) Jx = sig*Ex Jy = sig*Ey Jz = sig*Ez @@ -175,7 +175,7 @@ def J_galvanic_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., le Add description of parameters """ - Ex_galvanic, Ey_galvanic, Ez_galvanic = E_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.) + Ex_galvanic, Ey_galvanic, Ez_galvanic = E_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=current, length=length, orientation=orientation, kappa=kappa, epsr=epsr) Jx_galvanic = sig*Ex_galvanic Jy_galvanic = sig*Ey_galvanic Jz_galvanic = sig*Ez_galvanic @@ -190,7 +190,7 @@ def J_inductive_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., l Add description of parameters """ - Ex_inductive, Ey_inductive, Ez_inductive = E_inductive_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.) + Ex_inductive, Ey_inductive, Ez_inductive = E_inductive_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=current, length=length, orientation=orientation, kappa=kappa, epsr=epsr) Jx_inductive = sig*Ex_inductive Jy_inductive = sig*Ey_inductive Jz_inductive = sig*Ez_inductive @@ -248,7 +248,7 @@ def B_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., Add description of parameters """ - Hx, Hy, Hz = H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.) + Hx, Hy, Hz = H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=current, length=length, orientation=orientation, kappa=kappa, epsr=epsr) Bx = mu*Hx By = mu*Hy Bz = mu*Hz