Compare commits

..
Author SHA1 Message Date
Lindsey Heagy 8f1df12252 copy docs badges from readme to docs index 2016-07-21 22:30:49 -07:00
Lindsey Heagy 002554209b listing and remove wildcard inputs from tes_Mixed_boundaryPoisson 2016-07-21 22:28:20 -07:00
Lindsey Heagy 76f140f89b pep8 linting 2016-07-21 21:58:42 -07:00
Lindsey Heagy 5ce26a97ae remove wildcard imports from test_Fields 2016-07-21 21:39:11 -07:00
Lindsey Heagy a087c6e5b1 fix broken link 2016-07-21 21:34:17 -07:00
Lindsey Heagy 44ec0623d6 fixed formatting changes from quantified code, test links in docs 2016-07-21 14:38:05 -07:00
Lindsey Heagy 57d0c84b41 use @property for nodal and cell centered curvi grids 2016-07-21 13:19:23 -07:00
Lindsey Heagy a7befd6783 fix typo in cell grad 2016-07-21 12:38:21 -07:00
Lindsey Heagy 7acf46ddb8 fix typo in faceDivz 2016-07-21 12:32:49 -07:00
Lindsey Heagy b54d2494f5 add codacy to readme 2016-07-21 12:20:58 -07:00
Cody 45f5906554 Migrated % string formating 2016-07-21 12:08:19 -07:00
Lindsey Heagy 8093288391 fix double negative 2016-07-21 11:59:59 -07:00
Lindsey Heagy 7457912d84 avoid double import of numpy 2016-07-21 11:58:52 -07:00
Lindsey Heagy 910d4214b3 remove extra import of unittest 2016-07-21 11:56:49 -07:00
Lindsey Heagy 20e80ed983 use @property decorator in DiffOperators.py 2016-07-21 11:55:14 -07:00
Lindsey Heagy 3e8580f28d update readme 2016-07-21 09:59:43 -07:00
Lindsey Heagy 3ad2c98d43 use @property decorator in curve mesh (see: https://www.quantifiedcode.com/app/project/933aa3decf444538aa432c8817169b6d?groups=code_patterns%3A3bECxdfc%3Af0&tab=basics) 2016-07-21 09:58:14 -07:00
206 changed files with 2219 additions and 3273 deletions
+1 -2
View File
@@ -1,7 +1,6 @@
language: python
python:
- 2.7
- 3.4
sudo: false
@@ -39,7 +38,7 @@ before_install:
-O miniconda.sh; fi
- chmod +x miniconda.sh
- ./miniconda.sh -b
- export PATH=/home/travis/anaconda/bin:/home/travis/anaconda3/bin:/home/travis/miniconda/bin:/home/travis/miniconda3/bin:$PATH
- export PATH=/home/travis/anaconda/bin:/home/travis/miniconda/bin:$PATH
- conda update --yes conda
install:
+8 -5
View File
@@ -21,17 +21,20 @@ SimPEG
:target: https://travis-ci.org/simpeg/simpeg
:alt: Travis CI build status
.. image:: https://img.shields.io/coveralls/simpeg/simpeg.svg
:target: https://coveralls.io/r/simpeg/simpeg?branch=master
:alt: Coverage status
.. image:: http://img.shields.io/badge/GITTER-JOIN_CHAT-brightgreen.svg?style=flat-square
:alt: gitter chat room at https://gitter.im/simpeg/simpeg
:target: https://gitter.im/simpeg/simpeg
.. image:: https://codecov.io/gh/simpeg/simpeg/branch/master/graph/badge.svg
  :target: https://codecov.io/gh/simpeg/simpeg
   :target: https://codecov.io/gh/simpeg/simpeg
.. image:: https://www.quantifiedcode.com/api/v1/project/933aa3decf444538aa432c8817169b6d/badge.svg
:target: https://www.quantifiedcode.com/app/project/933aa3decf444538aa432c8817169b6d
:alt: Code issues
.. image:: https://api.codacy.com/project/badge/Grade/4fc959a5294a418fa21fc7bc3b3aa078
:target: https://www.codacy.com/app/lindseyheagy/simpeg?utm_source=github.com&utm_medium=referral&utm_content=simpeg/simpeg&utm_campaign=Badge_Grade
:alt: codacy
Simulation and Parameter Estimation in Geophysics - A python package for simulation and gradient based parameter estimation in the context of geophysical applications.
+1 -9
View File
@@ -1,11 +1,3 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from builtins import super
from future import standard_library
standard_library.install_aliases()
from builtins import range
from SimPEG import *
class FieldsDC_CC(Problem.Fields):
@@ -69,7 +61,7 @@ class SrcDipole(Survey.BaseSrc):
pts = [self.loc[0], self.loc[1]]
inds = Utils.closestPoints(prob.mesh, pts)
q = np.zeros(prob.mesh.nC)
q[inds] = - self.current * (np.r_[1., -1.] / prob.mesh.vol[inds])
q[inds] = - self.current * ( np.r_[1., -1.] / prob.mesh.vol[inds] )
# self._rhsDict[mesh] = q
# return self._rhsDict[mesh]
return q
+2 -9
View File
@@ -1,12 +1,5 @@
from __future__ import absolute_import
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from future import standard_library
standard_library.install_aliases()
from builtins import range
from SimPEG import *
from .BaseDC import SurveyDC, FieldsDC_CC
from BaseDC import SurveyDC, FieldsDC_CC
class SurveyIP(SurveyDC):
"""
@@ -59,7 +52,7 @@ class ProblemIP(Problem.BaseProblem):
# sigma = self.curModel.transform
sigma = self.sigma
Av = self.mesh.aveF2CC
self._Msig = Utils.sdiag(1//(self.mesh.dim * Av.T * (1/sigma)))
self._Msig = Utils.sdiag(1/(self.mesh.dim * Av.T * (1/sigma)))
return self._Msig
@property
+28 -38
View File
@@ -1,16 +1,6 @@
from __future__ import print_function
from __future__ import absolute_import
from __future__ import division
from __future__ import unicode_literals
from builtins import open
from builtins import int
from future import standard_library
standard_library.install_aliases()
from builtins import map
from builtins import range
from SimPEG import np, Utils
from . import BaseDC as DC
from . import BaseDC as IP
import BaseDC as DC
import BaseDC as IP
import warnings
def getActiveindfromTopo(mesh, topo):
@@ -77,7 +67,7 @@ def readUBC_DC3Dobstopo(filename,mesh,topo,probType="CC"):
if "!" in line.split(): continue
elif line == '\n': continue
elif line == ' \n': continue
temp = list(map(float, line.split()))
temp = map(float, line.split())
# Read a line for the current electrode
if len(temp) == 5: # SRC: Only X and Y are provided (assume no topography)
#TODO consider topography and assign the closest cell center in the earth
@@ -238,10 +228,10 @@ def plot_pseudoSection(DCsurvey, axs, surveyType='dipole-dipole', unitType='volt
elif surveyType == 'dipole-dipole':
leg = data * 2*np.pi / (1/MA - 1/MB - 1/NB + 1/NA)
leg = data * 2*np.pi / ( 1/MA - 1/MB - 1/NB + 1/NA )
else:
print("""unitType must be 'pole-dipole' | 'dipole-dipole' """)
print """unitType must be 'pole-dipole' | 'dipole-dipole' """
break
@@ -256,11 +246,11 @@ def plot_pseudoSection(DCsurvey, axs, surveyType='dipole-dipole', unitType='volt
rho = np.hstack([rho,leg])
else:
print("""unitType must be 'appResistivity' | 'appConductivity' | 'volt' """)
print """unitType must be 'appResistivity' | 'appConductivity' | 'volt' """
break
midx = np.hstack([midx, (Cmid + Pmid)/2])
midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + (Tx[0][2] + Tx[1][2])/2])
midx = np.hstack([midx, ( Cmid + Pmid )/2 ])
midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + (Tx[0][2] + Tx[1][2])/2 ])
# Grid points
grid_x, grid_z = np.mgrid[np.min(midx):np.max(midx), np.min(midz):np.max(midz)]
@@ -350,11 +340,11 @@ def gen_DCIPsurvey(endl, mesh, surveyType, AM_sep, MN_sep, nrx):
dl_x = ( endl[1,0] - endl[0,0] ) / dl_len
dl_y = ( endl[1,1] - endl[0,1] ) / dl_len
nstn = np.floor(dl_len / AM_sep)
nstn = np.floor( dl_len / AM_sep )
# Compute discrete pole location along line
stn_x = endl[0,0] + np.array(list(range(int(nstn))))*dl_x*AM_sep
stn_y = endl[0,1] + np.array(list(range(int(nstn))))*dl_y*AM_sep
stn_x = endl[0,0] + np.array(range(int(nstn)))*dl_x*AM_sep
stn_y = endl[0,1] + np.array(range(int(nstn)))*dl_y*AM_sep
# Create line of P1 locations
M = np.c_[stn_x, stn_y, np.ones(nstn).T*mesh.vectorNz[-1]]
@@ -386,15 +376,15 @@ def gen_DCIPsurvey(endl, mesh, surveyType, AM_sep, MN_sep, nrx):
AB = xy_2_r(tx[0,1],endl[1,0],tx[1,1],endl[1,1])
# Number of receivers to fit
nstn = np.min([(AB - MN_sep) // AM_sep, nrx])
nstn = np.min([np.floor( (AB - MN_sep) / AM_sep ) , nrx])
# Check if there is enough space, else break the loop
if nstn <= 0:
continue
# Compute discrete pole location along line
stn_x = N[ii,0] + dl_x*MN_sep + np.array(list(range(int(nstn))))*dl_x*AM_sep
stn_y = N[ii,1] + dl_y*MN_sep + np.array(list(range(int(nstn))))*dl_y*AM_sep
stn_x = N[ii,0] + dl_x*MN_sep + np.array(range(int(nstn)))*dl_x*AM_sep
stn_y = N[ii,1] + dl_y*MN_sep + np.array(range(int(nstn)))*dl_y*AM_sep
# Create receiver poles
# Create line of P1 locations
@@ -429,15 +419,15 @@ def gen_DCIPsurvey(endl, mesh, surveyType, AM_sep, MN_sep, nrx):
box_l = np.sqrt( (min_x - max_x)**2 + (min_y - max_y)**2 )
box_w = box_l/2.
nstn = box_l // AM_sep
nstn = np.floor( box_l / AM_sep )
# Compute discrete pole location along line
stn_x = min_x + np.array(list(range(int(nstn))))*dl_x*AM_sep
stn_y = min_y + np.array(list(range(int(nstn))))*dl_y*AM_sep
stn_x = min_x + np.array(range(int(nstn)))*dl_x*AM_sep
stn_y = min_y + np.array(range(int(nstn)))*dl_y*AM_sep
# Define number of cross lines
nlin = int(box_w // AM_sep)
lind = list(range(-nlin,nlin+1))
nlin = int(np.floor( box_w / AM_sep ))
lind = range(-nlin,nlin+1)
ngrad = nstn * len(lind)
@@ -459,7 +449,7 @@ def gen_DCIPsurvey(endl, mesh, surveyType, AM_sep, MN_sep, nrx):
srcClass = DC.SrcDipole([rxClass], M[0,:], N[-1,:])
SrcList.append(srcClass)
else:
print("""surveyType must be either 'pole-dipole', 'dipole-dipole' or 'gradient'. """)
print """surveyType must be either 'pole-dipole', 'dipole-dipole' or 'gradient'. """
survey = DC.SurveyDC(SrcList)
return survey, Tx, Rx
@@ -486,7 +476,7 @@ def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType, iptype = 0):
fid.write('! ' + surveyType + ' FORMAT\n')
if iptype!=0:
fid.write('IPTYPE=%i\n'%iptype)
fid.write('IPTYPE={0:d}\n'.format(iptype))
else:
fid.write('! ' + stype + ' FORMAT\n')
@@ -522,7 +512,7 @@ def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType, iptype = 0):
if surveyType == 'SURFACE':
fid.writelines("%f " % ii for ii in mkvc(tx[0,:]))
fid.writelines("{0:f} ".format(ii) for ii in mkvc(tx[0,:]))
M = M[:,0]
N = N[:,0]
@@ -531,7 +521,7 @@ def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType, iptype = 0):
# Flip sign for z-elevation to depth
tx[2::2,:] = -tx[2::2,:]
fid.writelines("%e " % ii for ii in mkvc(tx[::2,:]))
fid.writelines("{0:e} ".format(ii) for ii in mkvc(tx[::2,:]))
M = M[:,0::2]
N = N[:,0::2]
@@ -539,22 +529,22 @@ def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType, iptype = 0):
M[:,1::2] = -M[:,1::2]
N[:,1::2] = -N[:,1::2]
fid.write('%i\n'% nD)
fid.write('{0:d}\n'.format(nD))
np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%f',delimiter=' ',newline='\n')
if dim=='3D':
if surveyType == 'SURFACE':
fid.writelines("%e " % ii for ii in mkvc(tx[0:2,:]))
fid.writelines("{0:e} ".format(ii) for ii in mkvc(tx[0:2,:]))
M = M[:,0:2]
N = N[:,0:2]
if surveyType == 'GENERAL':
fid.writelines("%e " % ii for ii in mkvc(tx[0:3,:]))
fid.writelines("{0:e} ".format(ii) for ii in mkvc(tx[0:3,:]))
fid.write('%i\n'% nD)
fid.write('{0:d}\n'.format(nD))
np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n')
fid.write('\n')
@@ -678,7 +668,7 @@ def readUBC_DC3Dobs(fileName, rtype = 'DC'):
obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!')
else:
print("rtype must be 'DC'(default) | 'IP'")
print "rtype must be 'DC'(default) | 'IP'"
# Pre-allocate
srcLists = []
+1 -8
View File
@@ -1,10 +1,3 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
import numpy as np
def WennerSrcList(nElecs, aSpacing, in2D=False, plotIt=False):
@@ -12,7 +5,7 @@ def WennerSrcList(nElecs, aSpacing, in2D=False, plotIt=False):
import SimPEG.DCIP as DC
elocs = np.arange(0,aSpacing*nElecs,aSpacing)
elocs -= (nElecs*aSpacing - aSpacing) / 2
elocs -= (nElecs*aSpacing - aSpacing)/2
space = 1
WENNER = np.zeros((0,),dtype=int)
for ii in range(nElecs):
+4 -10
View File
@@ -1,10 +1,4 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from .BaseDC import *
from .BaseIP import *
from .DCIPUtils import *
from . import Utils
from BaseDC import *
from BaseIP import *
from DCIPUtils import *
import Utils
+6 -13
View File
@@ -1,16 +1,7 @@
from __future__ import print_function
from __future__ import division
from __future__ import unicode_literals
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import object
from . import Utils, Survey, Problem
import numpy as np, scipy.sparse as sp, gc
from future.utils import with_metaclass
import Utils, Survey, Problem, numpy as np, scipy.sparse as sp, gc
class BaseDataMisfit(with_metaclass(Utils.SimPEGMetaClass, object)):
class BaseDataMisfit(object):
"""BaseDataMisfit
.. note::
@@ -18,6 +9,8 @@ class BaseDataMisfit(with_metaclass(Utils.SimPEGMetaClass, object)):
You should inherit from this class to create your own data misfit term.
"""
__metaclass__ = Utils.SimPEGMetaClass
debug = False #: Print debugging information
counter = None #: Set this to a SimPEG.Utils.Counter() if you want to count things
@@ -100,11 +93,11 @@ class l2_DataMisfit(BaseDataMisfit):
survey = self.survey
if getattr(survey,'std', None) is None:
print('SimPEG.DataMisfit.l2_DataMisfit assigning default std of 5%')
print 'SimPEG.DataMisfit.l2_DataMisfit assigning default std of 5%'
survey.std = 0.05
if getattr(survey, 'eps', None) is None:
print('SimPEG.DataMisfit.l2_DataMisfit assigning default eps of 1e-5 * ||dobs||')
print 'SimPEG.DataMisfit.l2_DataMisfit assigning default eps of 1e-5 * ||dobs||'
survey.eps = np.linalg.norm(Utils.mkvc(survey.dobs),2)*1e-5
self._Wd = Utils.sdiag(1/(abs(survey.dobs)*survey.std+survey.eps))
+25 -36
View File
@@ -1,15 +1,4 @@
from __future__ import print_function
from __future__ import division
from __future__ import unicode_literals
from __future__ import absolute_import
from builtins import open
from builtins import int
from future import standard_library
standard_library.install_aliases()
from builtins import str
from builtins import object
from . import Utils
import numpy as np
import Utils, numpy as np
class InversionDirective(object):
"""InversionDirective"""
@@ -26,7 +15,7 @@ class InversionDirective(object):
@inversion.setter
def inversion(self, i):
if getattr(self,'_inversion',None) is not None:
print('Warning: InversionDirective %s has switched to a new inversion.' % self.__name__)
print 'Warning: InversionDirective {0!s} has switched to a new inversion.'.format(self.__name__)
self._inversion = i
@property
@@ -58,7 +47,7 @@ class DirectiveList(object):
def __init__(self, *directives, **kwargs):
self.dList = []
for d in directives:
assert isinstance(d, InversionDirective), 'All directives must be InversionDirectives not %s' % d.__name__
assert isinstance(d, InversionDirective), 'All directives must be InversionDirectives not {0!s}'.format(d.__name__)
self.dList.append(d)
Utils.setKwargs(self, **kwargs)
@@ -79,7 +68,7 @@ class DirectiveList(object):
def inversion(self, i):
if self.inversion is i: return
if getattr(self,'_inversion',None) is not None:
print('Warning: %s has switched to a new inversion.' % self.__name__)
print 'Warning: {0!s} has switched to a new inversion.'.format(self.__name__)
for d in self.dList:
d.inversion = i
self._inversion = i
@@ -90,7 +79,7 @@ class DirectiveList(object):
return
directives = ['initialize', 'endIter', 'finish']
assert ruleType in directives, 'Directive type must be in ["%s"]' % '", "'.join(directives)
assert ruleType in directives, 'Directive type must be in ["{0!s}"]'.format('", "'.join(directives))
for r in self.dList:
getattr(r, ruleType)()
@@ -131,7 +120,7 @@ class BetaEstimate_ByEig(InversionDirective):
:return: beta0
"""
if self.debug: print('Calculating the beta0 parameter.')
if self.debug: print 'Calculating the beta0 parameter.'
m = self.invProb.curModel
f = self.invProb.getFields(m, store=True, deleteWarmstart=False)
@@ -152,7 +141,7 @@ class BetaSchedule(InversionDirective):
def endIter(self):
if self.opt.iter > 0 and self.opt.iter % self.coolingRate == 0:
if self.debug: print('BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter)
if self.debug: print 'BetaSchedule is cooling Beta. Iteration: {0:d}'.format(self.opt.iter)
self.invProb.beta /= self.coolingFactor
@@ -192,7 +181,7 @@ class SaveEveryIteration(InversionDirective):
def fileName(self):
if getattr(self, '_fileName', None) is None:
from datetime import datetime
self._fileName = '%s-%s'%(self.name, datetime.now().strftime('%Y-%m-%d-%H-%M'))
self._fileName = '{0!s}-{1!s}'.format(self.name, datetime.now().strftime('%Y-%m-%d-%H-%M'))
return self._fileName
@fileName.setter
def fileName(self, value):
@@ -203,31 +192,31 @@ class SaveModelEveryIteration(SaveEveryIteration):
"""SaveModelEveryIteration"""
def initialize(self):
print("SimPEG.SaveModelEveryIteration will save your models as: '###-%s.npy'"%self.fileName)
print "SimPEG.SaveModelEveryIteration will save your models as: '###-{0!s}.npy'".format(self.fileName)
def endIter(self):
np.save('%03d-%s' % (self.opt.iter, self.fileName), self.opt.xc)
np.save('{0:03d}-{1!s}'.format(self.opt.iter, self.fileName), self.opt.xc)
class SaveOutputEveryIteration(SaveEveryIteration):
"""SaveModelEveryIteration"""
def initialize(self):
print("SimPEG.SaveOutputEveryIteration will save your inversion progress as: '###-%s.txt'"%self.fileName)
print "SimPEG.SaveOutputEveryIteration will save your inversion progress as: '###-{0!s}.txt'".format(self.fileName)
f = open(self.fileName+'.txt', 'w')
f.write(" # beta phi_d phi_m f\n")
f.close()
def endIter(self):
f = open(self.fileName+'.txt', 'a')
f.write(' %3d %1.4e %1.4e %1.4e %1.4e\n'%(self.opt.iter, self.invProb.beta, self.invProb.phi_d, self.invProb.phi_m, self.opt.f))
f.write(' {0:3d} {1:1.4e} {2:1.4e} {3:1.4e} {4:1.4e}\n'.format(self.opt.iter, self.invProb.beta, self.invProb.phi_d, self.invProb.phi_m, self.opt.f))
f.close()
class SaveOutputDictEveryIteration(SaveEveryIteration):
"""SaveOutputDictEveryIteration"""
def initialize(self):
print("SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '###-%s.npz'"%self.fileName)
print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '###-{0!s}.npz'".format(self.fileName)
def endIter(self):
# Save the data.
@@ -305,7 +294,7 @@ class Update_IRLS(InversionDirective):
# After reaching target misfit with l2-norm, switch to IRLS (mode:2)
if self.invProb.phi_d < self.target and self.mode == 1:
print("Convergence with smooth l2-norm regularization: Start IRLS steps...")
print "Convergence with smooth l2-norm regularization: Start IRLS steps..."
self.mode = 2
@@ -313,17 +302,17 @@ class Update_IRLS(InversionDirective):
# model values
if getattr(self, 'reg.eps', None) is None:
self.reg.eps_p = np.percentile(np.abs(self.invProb.curModel),self.prctile)
else:
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:
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))
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
@@ -339,7 +328,7 @@ class Update_IRLS(InversionDirective):
# Beta Schedule
if self.opt.iter > 0 and self.opt.iter % self.coolingRate == 0:
if self.debug: print('BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter)
if self.debug: print 'BetaSchedule is cooling Beta. Iteration: {0:d}'.format(self.opt.iter)
self.invProb.beta /= self.coolingFactor
@@ -351,17 +340,17 @@ class Update_IRLS(InversionDirective):
phim_new = self.reg.eval(self.invProb.curModel)
self.f_change = np.abs(self.f_old - phim_new) / self.f_old
print("Regularization decrease: %6.3e" % (self.f_change))
print "Regularization decrease: {0:6.3e}".format((self.f_change))
# Check for maximum number of IRLS cycles
if self.IRLSiter == self.maxIRLSiter:
print("Reach maximum number of IRLS cycles: %i" % self.maxIRLSiter)
print "Reach maximum number of IRLS cycles: {0:d}".format(self.maxIRLSiter)
self.opt.stopNextIteration = True
return
# Check if the function has changed enough
if self.f_change < self.f_min_change and self.IRLSiter > 1:
print("Minimum decrease in regularization. End of IRLS")
print "Minimum decrease in regularization. End of IRLS"
self.opt.stopNextIteration = True
return
else:
-7
View File
@@ -1,10 +1,3 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
import numpy as np
from scipy.constants import mu_0, pi
from scipy import special
-5
View File
@@ -1,9 +1,4 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
import numpy as np
from scipy.constants import mu_0, pi
from scipy.special import erf
-5
View File
@@ -1,9 +1,4 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
import numpy as np
from scipy.constants import mu_0, pi, epsilon_0
from scipy.special import erf
+1 -7
View File
@@ -1,9 +1,3 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import Utils, np
from scipy.constants import mu_0, epsilon_0
from SimPEG.EM.Utils.EMUtils import k
@@ -40,7 +34,7 @@ def _getCasingHertzMagDipoleDeriv_r(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones
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)
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.):
-6
View File
@@ -1,9 +1,3 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
import numpy as np
from scipy.constants import mu_0, pi
from scipy.special import erf
+5 -11
View File
@@ -1,11 +1,5 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from .TDEM import hzAnalyticDipoleT
from .FDEM import hzAnalyticDipoleF
from .FDEMcasing import *
from .DC import DCAnalyticHalf, DCAnalyticSphere
from .FDEMDipolarfields import *
from TDEM import hzAnalyticDipoleT
from FDEM import hzAnalyticDipoleF
from FDEMcasing import *
from DC import DCAnalyticHalf, DCAnalyticSphere
from FDEMDipolarfields import *
-6
View File
@@ -1,9 +1,3 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import Survey, Problem, Utils, Models, Maps, PropMaps, np, sp, Solver as SimpegSolver
from scipy.constants import mu_0
+29 -36
View File
@@ -1,10 +1,3 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from builtins import int
from future import standard_library
standard_library.install_aliases()
import numpy as np
import scipy.sparse as sp
import SimPEG
@@ -49,7 +42,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
:return: total electric field
"""
if getattr(self, '_ePrimary', None) is None or getattr(self, '_eSecondary', None) is None:
raise NotImplementedError ('Getting e from %s is not implemented' %list(self.knownFields.keys())[0])
raise NotImplementedError ('Getting e from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
return self._ePrimary(solution,srcList) + self._eSecondary(solution,srcList)
@@ -63,7 +56,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
:return: total magnetic flux density
"""
if getattr(self, '_bPrimary', None) is None or getattr(self, '_bSecondary', None) is None:
raise NotImplementedError ('Getting b from %s is not implemented' %list(self.knownFields.keys())[0])
raise NotImplementedError ('Getting b from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
return self._bPrimary(solution, srcList) + self._bSecondary(solution, srcList)
@@ -77,7 +70,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
:return: total magnetic field
"""
if getattr(self, '_hPrimary', None) is None or getattr(self, '_hSecondary', None) is None:
raise NotImplementedError ('Getting h from %s is not implemented' %list(self.knownFields.keys())[0])
raise NotImplementedError ('Getting h from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
return self._hPrimary(solution, srcList) + self._hSecondary(solution, srcList)
@@ -91,7 +84,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
:return: total current density
"""
if getattr(self, '_jPrimary', None) is None or getattr(self, '_jSecondary', None) is None:
raise NotImplementedError ('Getting j from %s is not implemented' %list(self.knownFields.keys())[0])
raise NotImplementedError ('Getting j from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
return self._jPrimary(solution, srcList) + self._jSecondary(solution, srcList)
@@ -107,7 +100,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
:return: derivative times a vector (or tuple for adjoint)
"""
if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None:
raise NotImplementedError ('Getting eDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
raise NotImplementedError ('Getting eDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
if adjoint:
return self._eDeriv_u(src, v, adjoint), self._eDeriv_m(src, v, adjoint)
@@ -125,7 +118,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
:return: derivative times a vector (or tuple for adjoint)
"""
if getattr(self, '_bDeriv_u', None) is None or getattr(self, '_bDeriv_m', None) is None:
raise NotImplementedError ('Getting bDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
raise NotImplementedError ('Getting bDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
if adjoint:
return self._bDeriv_u(src, v, adjoint), self._bDeriv_m(src, v, adjoint)
@@ -143,7 +136,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
:return: derivative times a vector (or tuple for adjoint)
"""
if getattr(self, '_hDeriv_u', None) is None or getattr(self, '_hDeriv_m', None) is None:
raise NotImplementedError ('Getting hDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
raise NotImplementedError ('Getting hDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
if adjoint:
return self._hDeriv_u(src, v, adjoint), self._hDeriv_m(src, v, adjoint)
@@ -161,7 +154,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
:return: derivative times a vector (or tuple for adjoint)
"""
if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None:
raise NotImplementedError ('Getting jDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
raise NotImplementedError ('Getting jDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
if adjoint:
return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint)
@@ -292,7 +285,7 @@ class Fields3D_e(FieldsFDEM):
C = self._edgeCurl
b = (C * eSolution)
for i, src in enumerate(srcList):
b[:,i] *= -1./(1j*omega(src.freq))
b[:,i] *= - 1./(1j*omega(src.freq))
s_m, _ = src.eval(self.prob)
b[:,i] = b[:,i]+ 1./(1j*omega(src.freq)) * s_m
return b
@@ -352,7 +345,7 @@ class Fields3D_e(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the derivative of the current density with respect to the field we solved for with a vector
"""
n = int(self._aveE2CCV.shape[0] // self._nC) # number of components (instead of checking if cyl or not)
n = int(self._aveE2CCV.shape[0] / self._nC) # number of components (instead of checking if cyl or not)
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
@@ -389,8 +382,8 @@ class Fields3D_e(FieldsFDEM):
:rtype: numpy.ndarray
:return: magnetic field
"""
n = int(self._aveF2CCV.shape[0] // self._nC) # Number of Components
VI = sdiag(np.kron(np.ones(n), 1. // self.prob.mesh.vol))
n = int(self._aveF2CCV.shape[0] / self._nC) # Number of Components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (self._aveF2CCV * (self._MfMui * self._b(eSolution, srcList)))
@@ -404,7 +397,7 @@ class Fields3D_e(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the derivative of the magnetic field with respect to the field we solved for with a vector
"""
n = int(self._aveF2CCV.shape[0] // self._nC) # Number of Components
n = int(self._aveF2CCV.shape[0] / self._nC) # Number of Components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
v = self._MfMui.T * (self._aveF2CCV.T * (VI.T * du_dm_v))
@@ -421,7 +414,7 @@ class Fields3D_e(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the magnetic field derivative with respect to the inversion model with a vector
"""
n = int(self._aveF2CCV.shape[0] // self._nC) # Number of Components
n = int(self._aveF2CCV.shape[0] / self._nC) # Number of Components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
v = self._MfMui.T * (self._aveF2CCV.T * (VI.T * v))
@@ -614,7 +607,7 @@ class Fields3D_b(FieldsFDEM):
:return: primary current density
"""
n = int(self._aveE2CCV.shape[0] // self._nC) # number of components
n = int(self._aveE2CCV.shape[0] / self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (self._aveE2CCV * ( self._MeSigma * self._e(bSolution,srcList ) ) )
@@ -631,7 +624,7 @@ class Fields3D_b(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the derivative of the current density with respect to the field we solved for with a vector
"""
n = int(self._aveE2CCV.shape[0] // self._nC) # number of components
n = int(self._aveE2CCV.shape[0] / self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
return self._MfMui.T * ( self._edgeCurl * ( self._aveE2CCV.T * (VI.T * du_dm_v) ) )
@@ -659,7 +652,7 @@ class Fields3D_b(FieldsFDEM):
:rtype: numpy.ndarray
:return: magnetic field
"""
n = int(self._aveF2CCV.shape[0] // self._nC) #number of components
n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (self._aveF2CCV * (self._MfMui * self._b(bSolution, srcList)))
@@ -674,7 +667,7 @@ class Fields3D_b(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the derivative of the magnetic field with respect to the field we solved for with a vector
"""
n = int(self._aveF2CCV.shape[0] // self._nC) #number of components
n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
@@ -897,7 +890,7 @@ class Fields3D_j(FieldsFDEM):
:rtype: numpy.ndarray
:return: electric field
"""
n = int(self._aveF2CCV.shape[0] // self._nC) # number of components
n = int(self._aveF2CCV.shape[0] / self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (self._aveF2CCV * (self._MfRho * self._j(jSolution, srcList)))
@@ -911,7 +904,7 @@ class Fields3D_j(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the derivative of the electric field with respect to the field we solved for with a vector
"""
n = int(self._aveF2CCV.shape[0] // self._nC) # number of components
n = int(self._aveF2CCV.shape[0] / self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
return self._MfRho.T * ( self._aveF2CCV.T * ( VI.T * du_dm_v ) )
@@ -928,7 +921,7 @@ class Fields3D_j(FieldsFDEM):
:return: product of the derivative of the electric field with respect to the model with a vector
"""
jSolution = Utils.mkvc(self[src,'jSolution'])
n = int(self._aveF2CCV.shape[0] // self._nC) # number of components
n = int(self._aveF2CCV.shape[0] / self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
return self._MfRhoDeriv(jSolution).T * ( self._aveF2CCV.T * ( VI.T * v ) )
@@ -943,7 +936,7 @@ class Fields3D_j(FieldsFDEM):
:rtype: numpy.ndarray
:return: secondary magnetic flux density
"""
n = int(self._aveE2CCV.shape[0] // self._nC) # number of components
n = int(self._aveE2CCV.shape[0] / self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (self._aveE2CCV * ( self._MeMu * self._h(jSolution,srcList)) )
@@ -958,7 +951,7 @@ class Fields3D_j(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the derivative of the magnetic flux density with respect to the field we solved for with a vector
"""
n = int(self._aveF2CCV.shape[0] // self._nC) # number of components
n = int(self._aveF2CCV.shape[0] / self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
@@ -976,7 +969,7 @@ class Fields3D_j(FieldsFDEM):
:return: product of the derivative of the magnetic flux density with respect to the model with a vector
"""
jSolution = self[src,'jSolution']
n = int(self._aveE2CCV.shape[0] // self._nC) # number of components
n = int(self._aveE2CCV.shape[0] / self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
s_mDeriv,_ = src.evalDeriv(self.prob, adjoint = adjoint)
@@ -1158,7 +1151,7 @@ class Fields3D_h(FieldsFDEM):
:rtype: numpy.ndarray
:return: electric field
"""
n = int(self._aveF2CCV.shape[0] // self._nC) #number of components
n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (self._aveF2CCV * (self._MfRho * self._j(hSolution, srcList)))
@@ -1172,7 +1165,7 @@ class Fields3D_h(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the derivative of the electric field with respect to the field we solved for with a vector
"""
n = int(self._aveF2CCV.shape[0] // self._nC) #number of components
n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
return self._edgeCurl.T * ( self._MfRho.T * ( self._aveF2CCV.T * ( VI.T * du_dm_v ) ) )
@@ -1189,7 +1182,7 @@ class Fields3D_h(FieldsFDEM):
:return: product of the electric field derivative with respect to the inversion model with a vector
"""
hSolution = Utils.mkvc(self[src,'hSolution'])
n = int(self._aveF2CCV.shape[0] // self._nC) #number of components
n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
return ( self._MfRhoDeriv(self._edgeCurl * hSolution).T * ( self._aveF2CCV.T * (VI.T * v) ) )
@@ -1205,7 +1198,7 @@ class Fields3D_h(FieldsFDEM):
:return: magnetic flux density
"""
h = self._h(hSolution, srcList)
n = int(self._aveE2CCV.shape[0] // self._nC) #number of components
n = int(self._aveE2CCV.shape[0] / self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (self._aveE2CCV * (self._MeMu * h))
@@ -1220,7 +1213,7 @@ class Fields3D_h(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the derivative of the magnetic flux density with respect to the field we solved for with a vector
"""
n = int(self._aveE2CCV.shape[0] // self._nC) #number of components
n = int(self._aveE2CCV.shape[0] / self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
return self._MeMu.T * (self._aveE2CCV.T * ( VI.T * du_dm_v ))
+2 -8
View File
@@ -1,13 +1,7 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from SimPEG import Problem, Utils, np, sp, Solver as SimpegSolver
from scipy.constants import mu_0
from .SurveyFDEM import Survey as SurveyFDEM
from .FieldsFDEM import FieldsFDEM, Fields3D_e, Fields3D_b, Fields3D_h, Fields3D_j
from SurveyFDEM import Survey as SurveyFDEM
from FieldsFDEM import FieldsFDEM, Fields3D_e, Fields3D_b, Fields3D_h, Fields3D_j
from SimPEG.EM.Base import BaseEMProblem
from SimPEG.EM.Utils import omega
+2 -9
View File
@@ -1,10 +1,3 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from builtins import super
from future import standard_library
standard_library.install_aliases()
import SimPEG
from SimPEG import sp
@@ -18,8 +11,8 @@ class BaseRx(SimPEG.Survey.BaseRx):
"""
def __init__(self, locs, orientation=None, component=None):
assert(orientation in ['x','y','z']), "Orientation %s not known. Orientation must be in 'x', 'y', 'z'. Arbitrary orientations have not yet been implemented."%orientation
assert(component in ['real', 'imag']), "'component' must be 'real' or 'imag', not %s"%component
assert(orientation in ['x','y','z']), "Orientation {0!s} not known. Orientation must be in 'x', 'y', 'z'. Arbitrary orientations have not yet been implemented.".format(orientation)
assert(component in ['real', 'imag']), "'component' must be 'real' or 'imag', not {0!s}".format(component)
self.projComp = orientation
self.component = component
-6
View File
@@ -1,9 +1,3 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import Survey, Problem, Utils, np, sp
from scipy.constants import mu_0
from SimPEG.EM.Utils import *
+2 -8
View File
@@ -1,16 +1,10 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
import SimPEG
from SimPEG.EM.Utils import *
from SimPEG.EM.Base import BaseEMSurvey
from scipy.constants import mu_0
from SimPEG.Utils import Zero, Identity
from . import SrcFDEM as Src
from . import RxFDEM as Rx
import SrcFDEM as Src
import RxFDEM as Rx
from SimPEG import sp
class Survey(BaseEMSurvey):
+5 -11
View File
@@ -1,11 +1,5 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from .SurveyFDEM import Survey
from . import SrcFDEM as Src
from . import RxFDEM as Rx
from .ProblemFDEM import Problem3D_e, Problem3D_b, Problem3D_j, Problem3D_h
from .FieldsFDEM import Fields3D_e, Fields3D_b, Fields3D_j, Fields3D_h
from SurveyFDEM import Survey
import SrcFDEM as Src
import RxFDEM as Rx
from ProblemFDEM import Problem3D_e, Problem3D_b, Problem3D_j, Problem3D_h
from FieldsFDEM import Fields3D_e, Fields3D_b, Fields3D_j, Fields3D_h
-6
View File
@@ -1,9 +1,3 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
import numpy as np
def getxBCyBC_CC(mesh, alpha, beta, gamma):
+3 -9
View File
@@ -1,9 +1,3 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
import SimPEG
from SimPEG.Utils import Identity, Zero
import numpy as np
@@ -15,7 +9,7 @@ class Fields(SimPEG.Problem.Fields):
def _phiDeriv(self, src, du_dm_v, v, adjoint=False):
if getattr(self, '_phiDeriv_u', None) is None or getattr(self, '_phiDeriv_m', None) is None:
raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
raise NotImplementedError ('Getting phiDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
if adjoint:
return self._phiDeriv_u(src, v, adjoint=adjoint), self._phiDeriv_m(src, v, adjoint=adjoint)
@@ -24,7 +18,7 @@ class Fields(SimPEG.Problem.Fields):
def _eDeriv(self, src, du_dm_v, v, adjoint=False):
if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None:
raise NotImplementedError ('Getting eDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
raise NotImplementedError ('Getting eDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
if adjoint:
return self._eDeriv_u(src, v, adjoint), self._eDeriv_m(src, v, adjoint)
@@ -32,7 +26,7 @@ class Fields(SimPEG.Problem.Fields):
def _jDeriv(self, src, du_dm_v, v, adjoint=False):
if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None:
raise NotImplementedError ('Getting jDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
raise NotImplementedError ('Getting jDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
if adjoint:
return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint)
+3 -9
View File
@@ -1,9 +1,3 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
import SimPEG
from SimPEG.Utils import Identity, Zero
import numpy as np
@@ -38,7 +32,7 @@ class Fields_ky(SimPEG.Problem.TimeFields):
def _phiDeriv(self,kyInd, src, du_dm_v, v, adjoint=False):
if getattr(self, '_phiDeriv_u', None) is None or getattr(self, '_phiDeriv_m', None) is None:
raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
raise NotImplementedError ('Getting phiDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
if adjoint:
return self._phiDeriv_u(kyInd, src, v, adjoint=adjoint), self._phiDeriv_m(kyInd, src, v, adjoint=adjoint)
@@ -47,7 +41,7 @@ class Fields_ky(SimPEG.Problem.TimeFields):
def _eDeriv(self,kyInd, src, du_dm_v, v, adjoint=False):
if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None:
raise NotImplementedError ('Getting eDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
raise NotImplementedError ('Getting eDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
if adjoint:
return self._eDeriv_u(kyInd, src, v, adjoint), self._eDeriv_m(kyInd, src, v, adjoint)
@@ -55,7 +49,7 @@ class Fields_ky(SimPEG.Problem.TimeFields):
def _jDeriv(self,kyInd, src, du_dm_v, v, adjoint=False):
if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None:
raise NotImplementedError ('Getting jDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
raise NotImplementedError ('Getting jDerivs from {0!s} is not implemented'.format(self.knownFields.keys()[0]))
if adjoint:
return self._jDeriv_u(kyInd, src, v, adjoint), self._jDeriv_m(kyInd, src, v, adjoint)
+5 -11
View File
@@ -1,17 +1,11 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from SimPEG import Problem, Utils
from SimPEG.EM.Base import BaseEMProblem
from .SurveyDC import Survey
from .FieldsDC import Fields, Fields_CC, Fields_N
from SurveyDC import Survey
from FieldsDC import Fields, Fields_CC, Fields_N
from SimPEG.Utils import sdiag
import numpy as np
from SimPEG.Utils import Zero
from .BoundaryUtils import getxBCyBC_CC
from BoundaryUtils import getxBCyBC_CC
class BaseDCProblem(BaseEMProblem):
@@ -52,7 +46,7 @@ class BaseDCProblem(BaseEMProblem):
du_dm_v = self.Ainv * ( - dA_dm_v + dRHS_dm_v )
for rx in src.rxList:
df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
df_dmFun = getattr(f, '_{0!s}Deriv'.format(rx.projField), None)
df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
return Utils.mkvc(Jv)
@@ -75,7 +69,7 @@ class BaseDCProblem(BaseEMProblem):
u_src = f[src, self._solutionType]
for rx in src.rxList:
PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m
df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
df_duTFun = getattr(f, '_{0!s}Deriv'.format(rx.projField), None)
df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True)
ATinvdf_duT = self.Ainv * df_duT
+5 -12
View File
@@ -1,18 +1,11 @@
from __future__ import absolute_import
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from future import standard_library
standard_library.install_aliases()
from builtins import range
from SimPEG import Problem, Utils
from SimPEG.EM.Base import BaseEMProblem
from .SurveyDC import Survey, Survey_ky
from .FieldsDC_2D import Fields_ky, Fields_ky_CC, Fields_ky_N
from SurveyDC import Survey, Survey_ky
from FieldsDC_2D import Fields_ky, Fields_ky_CC, Fields_ky_N
from SimPEG.Utils import sdiag
import numpy as np
from SimPEG.Utils import Zero
from .BoundaryUtils import getxBCyBC_CC
from BoundaryUtils import getxBCyBC_CC
class BaseDCProblem_2D(BaseEMProblem):
@@ -67,7 +60,7 @@ class BaseDCProblem_2D(BaseEMProblem):
dRHS_dm_v = self.getRHSDeriv(ky, src, v)
du_dm_v = self.Ainv[iky] * ( - dA_dm_v + dRHS_dm_v )
for rx in src.rxList:
df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
df_dmFun = getattr(f, '_{0!s}Deriv'.format(rx.projField), None)
df_dm_v = df_dmFun(iky, src, du_dm_v, v, adjoint=False)
# Trapezoidal intergration
Jv1_temp = 1./np.pi*rx.evalDeriv(ky, src, self.mesh, f, df_dm_v)
@@ -108,7 +101,7 @@ class BaseDCProblem_2D(BaseEMProblem):
ky = self.kys[iky]
AT = self.getA(ky)
PTv = rx.evalDeriv(ky, src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m
df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
df_duTFun = getattr(f, '_{0!s}Deriv'.format(rx.projField), None)
df_duT, df_dmT = df_duTFun(iky, src, None, PTv, adjoint=True)
ATinvdf_duT = self.Ainv[iky] * df_duT
-7
View File
@@ -1,10 +1,3 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
import SimPEG
import numpy as np
from SimPEG.Utils import Zero, closestPoints
-6
View File
@@ -1,9 +1,3 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
import SimPEG
# from SimPEG.EM.Base import BaseEMSurvey
from SimPEG.Utils import Zero, closestPoints, mkvc
+2 -8
View File
@@ -1,15 +1,9 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
import SimPEG
from SimPEG.EM.Base import BaseEMSurvey
from SimPEG import sp, Survey
from SimPEG.Utils import Zero, Identity
from .RxDC import BaseRx
from .SrcDC import BaseSrc
from RxDC import BaseRx
from SrcDC import BaseSrc
class Survey(BaseEMSurvey):
rxPair = BaseRx
-7
View File
@@ -1,10 +1,3 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
import numpy as np
def WennerSrcList(nElecs, aSpacing, in2D=False, plotIt=False):
+8 -14
View File
@@ -1,14 +1,8 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from .ProblemDC import Problem3D_CC, Problem3D_N
from .ProblemDC_2D import Problem2D_CC, Problem2D_N
from .SurveyDC import Survey, Survey_ky
from . import SrcDC as Src #Pole
from . import RxDC as Rx
from .FieldsDC import Fields_CC
from .BoundaryUtils import getxBCyBC_CC
from . import Utils
from ProblemDC import Problem3D_CC, Problem3D_N
from ProblemDC_2D import Problem2D_CC, Problem2D_N
from SurveyDC import Survey, Survey_ky
import SrcDC as Src #Pole
import RxDC as Rx
from FieldsDC import Fields_CC
from BoundaryUtils import getxBCyBC_CC
import Utils
+3 -9
View File
@@ -1,9 +1,3 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from SimPEG import Problem, Utils, Maps, Mesh
from SimPEG.EM.Base import BaseEMProblem
from SimPEG.EM.Static.DC.FieldsDC import Fields, Fields_CC, Fields_N
@@ -11,7 +5,7 @@ from SimPEG.Utils import sdiag
import numpy as np
from SimPEG.Utils import Zero
from SimPEG.EM.Static.DC import getxBCyBC_CC
from .SurveyIP import Survey
from SurveyIP import Survey
class IPPropMap(Maps.PropMap):
"""
@@ -62,7 +56,7 @@ class BaseIPProblem(BaseEMProblem):
du_dm_v = self.Ainv * ( - dA_dm_v + dRHS_dm_v )
for rx in src.rxList:
df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
df_dmFun = getattr(f, '_{0!s}Deriv'.format(rx.projField), None)
df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
# Conductivity (d u / d log sigma)
@@ -89,7 +83,7 @@ class BaseIPProblem(BaseEMProblem):
u_src = f[src, self._solutionType]
for rx in src.rxList:
PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m
df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
df_duTFun = getattr(f, '_{0!s}Deriv'.format(rx.projField), None)
df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True)
ATinvdf_duT = self.Ainv * df_duT
dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True)
-6
View File
@@ -1,9 +1,3 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
import SimPEG
from SimPEG.EM.Base import BaseEMSurvey
from SimPEG import sp, Survey
+2 -8
View File
@@ -1,8 +1,2 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from .ProblemIP import Problem3D_CC, Problem3D_N
from .SurveyIP import Survey
from ProblemIP import Problem3D_CC, Problem3D_N
from SurveyIP import Survey
+5 -13
View File
@@ -1,11 +1,3 @@
from __future__ import absolute_import
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from builtins import int
from future import standard_library
standard_library.install_aliases()
from builtins import range
from SimPEG import Problem, Utils, Maps, Mesh
from SimPEG.EM.Base import BaseEMProblem
from SimPEG.EM.Static.DC.FieldsDC import Fields, Fields_CC, Fields_N
@@ -13,7 +5,7 @@ from SimPEG.Utils import sdiag
import numpy as np
from SimPEG.Utils import Zero
from SimPEG.EM.Static.DC import getxBCyBC_CC
from .SurveySIP import Survey, Data
from SurveySIP import Survey, Data
class ColeColePropMap(Maps.PropMap):
"""
@@ -91,7 +83,7 @@ class BaseSIPProblem(BaseEMProblem):
for rx in src.rxList:
timeindex = rx.getTimeP(self.survey.times)
if timeindex[tind]:
df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
df_dmFun = getattr(f, '_{0!s}Deriv'.format(rx.projField), None)
df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
Jv[src, rx, t] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
@@ -113,7 +105,7 @@ class BaseSIPProblem(BaseEMProblem):
JvAll = []
#Assume only eta and tau (eta first then tau)
# v = [2*Mx1]
v = v.reshape((v.size//2), 2), order='F')
v = v.reshape((int(v.size/2), 2), order='F')
for tind in range(len(self.survey.times)):
t = self.survey.times[tind]
@@ -130,7 +122,7 @@ class BaseSIPProblem(BaseEMProblem):
for rx in src.rxList:
timeindex = rx.getTimeP(self.survey.times)
if timeindex[tind]:
df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
df_dmFun = getattr(f, '_{0!s}Deriv'.format(rx.projField), None)
df_dm_v0 = df_dmFun(src, du_dm_v0, v0, adjoint=False)
df_dm_v1 = df_dmFun(src, du_dm_v1, v1, adjoint=False)
Jv[src, rx, t] = rx.evalDeriv(src, self.mesh, f, df_dm_v0)
@@ -161,7 +153,7 @@ class BaseSIPProblem(BaseEMProblem):
timeindex = rx.getTimeP(self.survey.times)
if timeindex[tind]:
PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx, t], adjoint=True) # wrt f, need possibility wrt m
df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
df_duTFun = getattr(f, '_{0!s}Deriv'.format(rx.projField), None)
df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True)
ATinvdf_duT = self.Ainv * df_duT
dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True)
-7
View File
@@ -1,10 +1,3 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
from SimPEG import Utils, Maps, Mesh, sp, np
from SimPEG.Regularization import BaseRegularization, Simple
-6
View File
@@ -1,9 +1,3 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
import SimPEG
import numpy as np
from SimPEG.Utils import Zero, closestPoints
-6
View File
@@ -1,9 +1,3 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
import SimPEG
# from SimPEG.EM.Base import BaseEMSurvey
from SimPEG.Utils import Zero, closestPoints, mkvc
-7
View File
@@ -1,10 +1,3 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import str
import SimPEG
from SimPEG.EM.Base import BaseEMSurvey
from SimPEG import np, sp, Survey, Utils
+5 -11
View File
@@ -1,11 +1,5 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from .ProblemSIP import Problem3D_CC, Problem3D_N
from .SurveySIP import Survey, Data
from . import SrcSIP as Src #Pole
from . import RxSIP as Rx
from .Regularization import MultiRegularization
from ProblemSIP import Problem3D_CC, Problem3D_N
from SurveySIP import Survey, Data
import SrcSIP as Src #Pole
import RxSIP as Rx
from Regularization import MultiRegularization
+19 -27
View File
@@ -1,11 +1,3 @@
from __future__ import print_function
from __future__ import division
from __future__ import unicode_literals
from __future__ import absolute_import
from builtins import int
from future import standard_library
standard_library.install_aliases()
from builtins import range
from SimPEG import np
from SimPEG.EM.Static import DC, IP
@@ -93,10 +85,10 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None):
elif stype == 'dpdp':
leg = data * 2*np.pi / (1/MA - 1/MB + 1/NB - 1/NA)
LEG.append(1./(2*np.pi) * (1/MA - 1/MB + 1/NB - 1/NA))
leg = data * 2*np.pi / ( 1/MA - 1/MB + 1/NB - 1/NA )
LEG.append(1./(2*np.pi) *( 1/MA - 1/MB + 1/NB - 1/NA ))
else:
print("""dtype must be 'pdp'(pole-dipole) | 'dpdp' (dipole-dipole) """)
print """dtype must be 'pdp'(pole-dipole) | 'dpdp' (dipole-dipole) """
break
@@ -111,7 +103,7 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None):
rho = np.hstack([rho,leg])
else:
print("""dtype must be 'appr' | 'appc' | 'volt' """)
print """dtype must be 'appr' | 'appc' | 'volt' """
break
@@ -192,14 +184,14 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
# Mesure survey length and direction
dl_len = xy_2_r(endl[0,0],endl[1,0],endl[0,1],endl[1,1])
dl_x = (endl[1,0] - endl[0,0]) / dl_len
dl_y = (endl[1,1] - endl[0,1]) / dl_len
dl_x = ( endl[1,0] - endl[0,0] ) / dl_len
dl_y = ( endl[1,1] - endl[0,1] ) / dl_len
nstn = np.floor(dl_len / a)
nstn = np.floor( dl_len / a )
# Compute discrete pole location along line
stn_x = endl[0,0] + np.array(list(range(int(nstn))))*dl_x*a
stn_y = endl[0,1] + np.array(list(range(int(nstn))))*dl_y*a
stn_x = endl[0,0] + np.array(range(int(nstn)))*dl_x*a
stn_y = endl[0,1] + np.array(range(int(nstn)))*dl_y*a
if mesh.dim==2:
ztop = mesh.vectorNy[-1]
@@ -238,15 +230,15 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
AB = xy_2_r(tx[0,1],endl[1,0],tx[1,1],endl[1,1])
# Number of receivers to fit
nstn = np.min([(AB - b) // a, n])
nstn = np.min([np.floor( (AB - b) / a ) , n])
# Check if there is enough space, else break the loop
if nstn <= 0:
continue
# Compute discrete pole location along line
stn_x = N[ii,0] + dl_x*b + np.array(list(range(int(nstn))))*dl_x*a
stn_y = N[ii,1] + dl_y*b + np.array(list(range(int(nstn))))*dl_y*a
stn_x = N[ii,0] + dl_x*b + np.array(range(int(nstn)))*dl_x*a
stn_y = N[ii,1] + dl_y*b + np.array(range(int(nstn)))*dl_y*a
# Create receiver poles
@@ -283,17 +275,17 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
max_y = endl[1,1] - dl_y * b
box_l = np.sqrt( (min_x - max_x)**2 + (min_y - max_y)**2 )
box_w = box_l / 2.
box_w = box_l/2.
nstn = np.floor(box_l / a)
nstn = np.floor( box_l / a )
# Compute discrete pole location along line
stn_x = min_x + np.array(list(range(int(nstn))))*dl_x*a
stn_y = min_y + np.array(list(range(int(nstn))))*dl_y*a
stn_x = min_x + np.array(range(int(nstn)))*dl_x*a
stn_y = min_y + np.array(range(int(nstn)))*dl_y*a
# Define number of cross lines
nlin = int(box_w // a)
lind = list(range(-nlin,nlin+1))
nlin = int(np.floor( box_w / a ))
lind = range(-nlin,nlin+1)
ngrad = nstn * len(lind)
@@ -318,7 +310,7 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
srcClass = DC.Src.Dipole([rxClass], M[0,:], N[-1,:])
SrcList.append(srcClass)
else:
print("""stype must be either 'pdp', 'dpdp' or 'gradient'. """)
print """stype must be either 'pdp', 'dpdp' or 'gradient'. """
return SrcList
+1 -7
View File
@@ -1,7 +1 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from .StaticUtils import *
from StaticUtils import *
+3 -9
View File
@@ -1,9 +1,3 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from . import DC
from . import IP
from . import SIP
import DC
import IP
import SIP
+14 -21
View File
@@ -1,10 +1,3 @@
from __future__ import print_function
from __future__ import unicode_literals
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
from SimPEG import Solver, Problem
from SimPEG.Problem import BaseTimeProblem
from SimPEG.EM.Utils import *
@@ -54,7 +47,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
self.waveformType = "GENERAL"
def fields(self, m):
if self.verbose: print('%s\nCalculating fields(m)\n%s'%('*'*50,'*'*50))
if self.verbose: print '{0!s}\nCalculating fields(m)\n{1!s}'.format('*'*50, '*'*50)
self.curModel = m
# Create a fields storage object
F = self._FieldsForward_pair(self.mesh, self.survey)
@@ -62,7 +55,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
# Set the initial conditions
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))
if self.verbose: print '{0!s}\nDone calculating fields(m)\n{1!s}'.format('*'*50, '*'*50)
return F
def forward(self, m, RHS, F=None):
@@ -77,13 +70,13 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
if Ainv is not None:
Ainv.clean()
A = self.getA(tInd)
if self.verbose: print('Factoring... (dt = %e)'%dt)
if self.verbose: print 'Factoring... (dt = {0:e})'.format(dt)
Ainv = self.Solver(A, **self.solverOpts)
if self.verbose: print('Done')
if self.verbose: print 'Done'
rhs = RHS(tInd, F)
if self.verbose: print(' Solving... (tInd = %d)'%tInd)
if self.verbose: print ' Solving... (tInd = {0:d})'.format(tInd)
sol = Ainv * rhs
if self.verbose: print(' Done...')
if self.verbose: print ' Done...'
if sol.ndim == 1:
sol.shape = (sol.size,1)
F[:,self.solType,tInd+1] = sol
@@ -102,13 +95,13 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
if Ainv is not None:
Ainv.clean()
A = self.getA(tInd)
if self.verbose: print('Factoring (Adjoint)... (dt = %e)'%dt)
if self.verbose: print 'Factoring (Adjoint)... (dt = {0:e})'.format(dt)
Ainv = self.Solver(A, **self.solverOpts)
if self.verbose: print('Done')
if self.verbose: print 'Done'
rhs = RHS(tInd, F)
if self.verbose: print(' Solving (Adjoint)... (tInd = %d)'%tInd)
if self.verbose: print ' Solving (Adjoint)... (tInd = {0:d})'.format(tInd)
sol = Ainv * rhs
if self.verbose: print(' Done...')
if self.verbose: print ' Done...'
if sol.ndim == 1:
sol.shape = (sol.size,1)
F[:,self.solType,tInd+1] = sol
@@ -130,14 +123,14 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
* Compute \\\(\\\\vec{w} = -\\\mathbf{Q} \\\\vec{y}\\\)
"""
if self.verbose: print('%s\nCalculating J(v)\n%s'%('*'*50,'*'*50))
if self.verbose: print '{0!s}\nCalculating J(v)\n{1!s}'.format('*'*50, '*'*50)
self.curModel = m
if f is None:
f = self.fields(m)
p = self.Gvec(m, v, f)
y = self.solveAh(m, p)
Jv = self.survey.evalDeriv(f, v=y)
if self.verbose: print('%s\nDone calculating J(v)\n%s'%('*'*50,'*'*50))
if self.verbose: print '{0!s}\nDone calculating J(v)\n{1!s}'.format('*'*50, '*'*50)
return - mkvc(Jv)
def Jtvec(self, m, v, f=None):
@@ -155,7 +148,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
* Compute \\\(\\\\vec{w} = -\\\mathbf{G}^\\\\top y\\\)
"""
if self.verbose: print('%s\nCalculating J^T(v)\n%s'%('*'*50,'*'*50))
if self.verbose: print '{0!s}\nCalculating J^T(v)\n{1!s}'.format('*'*50, '*'*50)
self.curModel = m
if f is None:
f = self.fields(m)
@@ -166,6 +159,6 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
p = self.survey.evalDeriv(f, v=v, adjoint=True)
y = self.solveAht(m, p)
w = self.Gtvec(m, y, f)
if self.verbose: print('%s\nDone calculating J^T(v)\n%s'%('*'*50,'*'*50))
if self.verbose: print '{0!s}\nDone calculating J^T(v)\n{1!s}'.format('*'*50, '*'*50)
return - mkvc(w)
+5 -11
View File
@@ -1,13 +1,7 @@
from __future__ import print_function
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from SimPEG import Utils, Survey, np
from SimPEG.Survey import BaseSurvey
from SimPEG.EM.Utils import *
from .BaseTDEM import FieldsTDEM
from BaseTDEM import FieldsTDEM
class RxTDEM(Survey.BaseTimeRx):
@@ -93,7 +87,7 @@ class SrcTDEM_VMD_MVP(SrcTDEM):
def getInitialFields(self, mesh):
"""Vertical magnetic dipole, magnetic vector potential"""
if self.waveformType == "STEPOFF":
print(">> Step waveform: Non-zero initial condition")
print ">> Step waveform: Non-zero initial condition"
if mesh._meshType is 'CYL':
if mesh.isSymmetric:
MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
@@ -105,7 +99,7 @@ class SrcTDEM_VMD_MVP(SrcTDEM):
raise Exception('Unknown mesh for VMD')
return {"b": mesh.edgeCurl*MVP}
elif self.waveformType == "GENERAL":
print(">> General waveform: Zero initial condition")
print ">> General waveform: Zero initial condition"
return {"b": np.zeros(mesh.nF)}
else:
raise NotImplementedError("Only use STEPOFF or GENERAL")
@@ -133,7 +127,7 @@ class SrcTDEM_CircularLoop_MVP(SrcTDEM):
def getInitialFields(self, mesh):
"""Circular Loop, magnetic vector potential"""
if self.waveformType == "STEPOFF":
print(">> Step waveform: Non-zero initial condition")
print ">> Step waveform: Non-zero initial condition"
if mesh._meshType is 'CYL':
if mesh.isSymmetric:
MVP = MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
@@ -145,7 +139,7 @@ class SrcTDEM_CircularLoop_MVP(SrcTDEM):
raise Exception('Unknown mesh for CircularLoop')
return {"b": mesh.edgeCurl*MVP}
elif self.waveformType == "GENERAL":
print(">> General waveform: Zero initial condition")
print ">> General waveform: Zero initial condition"
return {"b": np.zeros(mesh.nF)}
else:
raise NotImplementedError("Only use STEPOFF or GENERAL")
+2 -9
View File
@@ -1,14 +1,7 @@
from __future__ import absolute_import
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from future import standard_library
standard_library.install_aliases()
from builtins import range
from .BaseTDEM import BaseTDEMProblem, FieldsTDEM
from BaseTDEM import BaseTDEMProblem, FieldsTDEM
from SimPEG.Utils import mkvc, sdiag
import numpy as np
from .SurveyTDEM import SurveyTDEM
from SurveyTDEM import SurveyTDEM
class FieldsTDEM_e_from_b(FieldsTDEM):
+3 -9
View File
@@ -1,9 +1,3 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from .SurveyTDEM import * #SurveyTDEM, RxTDEM, SrcTDEM
from .BaseTDEM import BaseTDEMProblem, FieldsTDEM
from .TDEM_b import ProblemTDEM_b
from SurveyTDEM import * #SurveyTDEM, RxTDEM, SrcTDEM
from BaseTDEM import BaseTDEMProblem, FieldsTDEM
from TDEM_b import ProblemTDEM_b
+5 -12
View File
@@ -1,10 +1,3 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
from SimPEG import *
from scipy.special import ellipk, ellipe
from scipy.constants import mu_0, pi
@@ -24,7 +17,7 @@ def MagneticDipoleVectorPotential(srcLoc, obsLoc, component, moment=1., dipoleMo
#TODO: break this out!
if type(component) in [list, tuple]:
out = list(range(len(component)))
out = range(len(component))
for i, comp in enumerate(component):
out[i] = MagneticDipoleVectorPotential(srcLoc, obsLoc, comp, dipoleMoment=dipoleMoment)
return np.concatenate(out)
@@ -125,7 +118,7 @@ def MagneticLoopVectorPotential(srcLoc, obsLoc, component, radius, mu=mu_0):
"""
if type(component) in [list, tuple]:
out = list(range(len(component)))
out = range(len(component))
for i, comp in enumerate(component):
out[i] = MagneticLoopVectorPotential(srcLoc, obsLoc, comp, radius, mu)
return np.concatenate(out)
@@ -165,11 +158,11 @@ def MagneticLoopVectorPotential(srcLoc, obsLoc, component, radius, mu=mu_0):
# % 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])
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])
A[ind, i] = Aphi[ind] * (-y[ind] / r[ind] )
elif component == 'y':
A[ind, i] = Aphi[ind] * (x[ind] / r[ind])
A[ind, i] = Aphi[ind] * ( x[ind] / r[ind] )
else:
raise ValueError('Invalid component')
+2 -8
View File
@@ -1,9 +1,3 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
import numpy as np
from scipy.constants import mu_0, epsilon_0
@@ -15,8 +9,8 @@ def omega(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)
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
+2 -8
View File
@@ -1,8 +1,2 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from .EMUtils import omega, k
from .AnalyticUtils import MagneticDipoleFields, MagneticDipoleVectorPotential, MagneticLoopVectorPotential
from EMUtils import omega, k
from AnalyticUtils import MagneticDipoleFields, MagneticDipoleVectorPotential, MagneticLoopVectorPotential
+6 -13
View File
@@ -1,10 +1,3 @@
from __future__ import print_function
from __future__ import division
from __future__ import unicode_literals
from __future__ import absolute_import
from builtins import int
from future import standard_library
standard_library.install_aliases()
import unittest
from SimPEG import *
from SimPEG import EM
@@ -65,7 +58,7 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False):
Src.append(EM.FDEM.Src.RawVec([rx0], freq, mesh.getEdgeInnerProduct()*S_m, S_e))
if verbose:
print(' Fetching %s problem' % (fdemType))
print ' Fetching {0!s} problem'.format((fdemType))
if fdemType == 'e':
survey = EM.FDEM.Survey(Src)
@@ -90,7 +83,7 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False):
try:
from pymatsolver import MumpsSolver
prb.Solver = MumpsSolver
except ImportError as e:
except ImportError, e:
prb.Solver = SolverLU
return prb
@@ -101,7 +94,7 @@ def crossCheckTest(SrcList, fdemType1, fdemType2, comp, addrandoms = False, useM
prb1 = getFDEMProblem(fdemType1, comp, SrcList, freq, useMu, verbose)
mesh = prb1.mesh
print('Cross Checking Forward: %s, %s formulations - %s' % (fdemType1, fdemType2, comp))
print 'Cross Checking Forward: {0!s}, {1!s} formulations - {2!s}'.format(fdemType1, fdemType2, comp)
logsig = np.log(np.ones(mesh.nC)*CONDUCTIVITY)
mu = np.ones(mesh.nC)*MU
@@ -119,7 +112,7 @@ def crossCheckTest(SrcList, fdemType1, fdemType2, comp, addrandoms = False, useM
d1 = survey1.dpred(m)
if verbose:
print(' Problem 1 solved')
print ' Problem 1 solved'
prb2 = getFDEMProblem(fdemType2, comp, SrcList, freq, useMu, verbose)
@@ -128,11 +121,11 @@ def crossCheckTest(SrcList, fdemType1, fdemType2, comp, addrandoms = False, useM
d2 = survey2.dpred(m)
if verbose:
print(' Problem 2 solved')
print ' Problem 2 solved'
r = d2-d1
l2r = l2norm(r)
tol = np.max([TOL*(10**int(np.log10(0.5* (l2norm(d1) + l2norm(d2)) ))),FLR])
print(l2norm(d1), l2norm(d2), l2r , tol, l2r < tol)
print l2norm(d1), l2norm(d2), l2r , tol, l2r < tol
return l2r < tol
+6 -12
View File
@@ -1,13 +1,7 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from . import TDEM
from . import FDEM
from . import Static
from . import Base
from . import Analytics
from . import Utils
import TDEM
import FDEM
import Static
import Base
import Analytics
import Utils
from scipy.constants import mu_0, epsilon_0
+3 -9
View File
@@ -1,9 +1,3 @@
from __future__ import print_function
from __future__ import division
from __future__ import unicode_literals
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import *
import SimPEG.EM.Static.DC as DC
@@ -35,7 +29,7 @@ def run(plotIt=True):
try:
from pymatsolver import MumpsSolver
problem.Solver = MumpsSolver
except Exception as e:
except Exception, e:
pass
data = survey.dpred(sigma)
@@ -67,8 +61,8 @@ def run(plotIt=True):
ax[0].set_title('Computed')
plt.show()
return np.linalg.norm(data-data_ana) / np.linalg.norm(data_ana)
return np.linalg.norm(data-data_ana)/np.linalg.norm(data_ana)
if __name__ == '__main__':
print(run())
print run()
+6 -14
View File
@@ -1,11 +1,3 @@
from __future__ import print_function
from __future__ import division
from __future__ import unicode_literals
from __future__ import absolute_import
from builtins import int
from future import standard_library
standard_library.install_aliases()
from builtins import range
from SimPEG import Mesh, Utils, np, sp
import SimPEG.DCIP as DC
import time
@@ -65,7 +57,7 @@ def run(loc=None, sig=None, radi=None, param=None, surveyType='dipole-dipole', u
model[ind] = sig[2]
# Get index of the center
indy = int(mesh.nCy // 2)
indy = int(mesh.nCy/2)
# Plot the model for reference
# Define core mesh extent
@@ -132,10 +124,10 @@ def run(loc=None, sig=None, radi=None, param=None, surveyType='dipole-dipole', u
tx = np.squeeze(Tx[ii][:,0:1])
tinf = tx + np.array([dl_x,dl_y,0])*dl_len*2
inds = Utils.closestPoints(mesh, np.c_[tx,tinf].T)
RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*([-1] / mesh.vol[inds])
RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*( [-1] / mesh.vol[inds] )
else:
inds = Utils.closestPoints(mesh, np.asarray(Tx[ii]).T )
RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*([-1,1] / mesh.vol[inds])
RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*( [-1,1] / mesh.vol[inds] )
# Iterative Solve
Ainvb = sp.linalg.bicgstab(P*A,P*RHS, tol=1e-5)
@@ -151,10 +143,10 @@ def run(loc=None, sig=None, radi=None, param=None, surveyType='dipole-dipole', u
dtemp = (P1*phi - P2*phi)*np.pi
data.append( dtemp )
print('\rTransmitter {0} of {1} -> Time:{2} sec'.format(ii,len(Tx),time.time()- start_time), end=' ')
print '\rTransmitter {0} of {1} -> Time:{2} sec'.format(ii,len(Tx),time.time()- start_time),
print('Transmitter {0} of {1}'.format(ii,len(Tx)))
print('Forward completed')
print 'Transmitter {0} of {1}'.format(ii,len(Tx))
print 'Forward completed'
# Let's just convert the 3D format into 2D (distance along line) and plot
survey2D = DC.convertObs_DC3D_to_2D(survey, np.ones(survey.nSrc) , 'Xloc')
+1 -7
View File
@@ -1,9 +1,3 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import *
import SimPEG.EM as EM
from SimPEG.EM import mu_0
@@ -62,7 +56,7 @@ def run(plotIt=True):
try:
from pymatsolver import MumpsSolver
prb.Solver = MumpsSolver
except ImportError as e:
except ImportError, e:
prb.Solver = SolverLU
prb.pair(survey)
@@ -1,9 +1,3 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import *
import SimPEG.EM as EM
+14 -21
View File
@@ -1,10 +1,3 @@
from __future__ import print_function
from __future__ import division
from __future__ import unicode_literals
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
from SimPEG import *
from SimPEG.EM import FDEM, Analytics, mu_0
import time
@@ -74,8 +67,8 @@ def run(plotIt=True):
casing_l = 300 # length of the casing
casing_r = 0.1
casing_a = casing_r - casing_t / 2. # inner radius
casing_b = casing_r + casing_t / 2. # outer radius
casing_a = casing_r - casing_t/2. # inner radius
casing_b = casing_r + casing_t/2. # outer radius
casing_z = np.r_[-casing_l,0.]
@@ -85,25 +78,25 @@ def run(plotIt=True):
src_loc = np.r_[0.,0.,dsz]
inf_loc = np.r_[0.,0.,1e4]
print('Skin Depth: ', [(500. / np.sqrt(sigmaback*_)) for _ in freqs])
print 'Skin Depth: ', [(500./np.sqrt(sigmaback*_)) for _ in freqs]
# ------------------ MESH ------------------
# fine cells near well bore
csx1, csx2 = 2e-3, 60.
pfx1, pfx2 = 1.3, 1.3
ncx1 = np.ceil(casing_b/csx1)+2
ncx1 = np.ceil(casing_b/csx1+2)
# pad nicely to second cell size
npadx1 = np.log(csx2/csx1) // np.log(pfx1)
npadx1 = np.floor(np.log(csx2/csx1) / np.log(pfx1))
hx1a,hx1b = Utils.meshTensor([(csx1,ncx1)]),Utils.meshTensor([(csx1,npadx1,pfx1)])
dx1 = sum(hx1a)+sum(hx1b)
dx1 = dx1 // csx2
hx1b *= (dx1*csx2 - sum(hx1a)) / sum(hx1b)
dx1 = np.floor(dx1/csx2)
hx1b *= (dx1*csx2 - sum(hx1a))/sum(hx1b)
# second chunk of mesh
dx2 = 300. # uniform mesh out to here
ncx2 = np.ceil((dx2 - dx1) / csx2)
ncx2 = np.ceil((dx2 - dx1)/csx2)
npadx2 = 45
hx2a, hx2b = Utils.meshTensor([(csx2,ncx2)]), Utils.meshTensor([(csx2,npadx2,pfx2)])
hx = np.hstack([hx1a,hx1b,hx2a,hx2b])
@@ -117,8 +110,8 @@ def run(plotIt=True):
# Mesh
mesh = Mesh.CylMesh([hx,1.,hz], [0.,0.,-np.sum(hz[:npadzu+ncz-nza])])
print('Mesh Extent xmax: %f,: zmin: %f, zmax: %f'%(mesh.vectorCCx.max(), mesh.vectorCCz.min(), mesh.vectorCCz.max()))
print('Number of cells', mesh.nC)
print 'Mesh Extent xmax: {0:f},: zmin: {1:f}, zmax: {2:f}'.format(mesh.vectorCCx.max(), mesh.vectorCCz.min(), mesh.vectorCCz.max())
print 'Number of cells', mesh.nC
if plotIt is True:
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
@@ -231,7 +224,7 @@ def run(plotIt=True):
# ------------- Solve ---------------------------
t0 = time.time()
fieldsCasing = problem.fields(sigCasing)
print('Time to solve 2 sources', time.time() - t0)
print 'Time to solve 2 sources', time.time() - t0
# Plot current
@@ -258,9 +251,9 @@ def run(plotIt=True):
in1_in = in1[np.r_[inds]]
z_in = mesh.gridFz[inds_fz,2]
in0_in = in0_in.reshape([in0_in.shape[0]//3,3])
in1_in = in1_in.reshape([in1_in.shape[0]//3,3])
z_in = z_in.reshape([z_in.shape[0]//3,3])
in0_in = in0_in.reshape([in0_in.shape[0]/3,3])
in1_in = in1_in.reshape([in1_in.shape[0]/3,3])
z_in = z_in.reshape([z_in.shape[0]/3,3])
I0 = in0_in.sum(1).real
I1 = in1_in.sum(1).real
-6
View File
@@ -1,9 +1,3 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import *
import SimPEG.EM as EM
from SimPEG.EM import mu_0
@@ -1,9 +1,3 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import *
from SimPEG.FLOW import Richards
+5 -13
View File
@@ -1,11 +1,3 @@
from __future__ import print_function
from __future__ import division
from __future__ import unicode_literals
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import str
from builtins import range
from SimPEG import *
@@ -55,7 +47,7 @@ def run(N=100, plotIt=True):
# Distance weighting
wr = np.sum(prob.G**2.,axis=0)**0.5
wr = ( wr/np.max(wr))
wr = ( wr/np.max(wr) )
dmis = DataMisfit.l2_DataMisfit(survey)
dmis.Wd = 1./wd
@@ -67,15 +59,15 @@ def run(N=100, plotIt=True):
reg.cell_weights = wr
reg.mref = np.zeros(mesh.nC)
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()
# 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.]
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])
@@ -83,7 +75,7 @@ def run(N=100, plotIt=True):
# Run inversion
mrec = inv.run(m0)
print("Final misfit:" + str(invProb.dmisfit.eval(mrec)))
print "Final misfit:" + str(invProb.dmisfit.eval(mrec))
if plotIt:
-7
View File
@@ -1,10 +1,3 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
from SimPEG import *
@@ -1,9 +1,3 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
import SimPEG as simpeg
import numpy as np
import SimPEG.MT as MT
-6
View File
@@ -1,9 +1,3 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
# Test script to use SimPEG.MT platform to forward model synthetic data.
# Import
-7
View File
@@ -1,10 +1,3 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from builtins import dict
from future import standard_library
standard_library.install_aliases()
from SimPEG import Mesh, Maps, np
def run(plotIt=True):
-6
View File
@@ -1,9 +1,3 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import Mesh, Maps, Utils
def run(plotIt=True):
-6
View File
@@ -1,9 +1,3 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import Mesh, Utils, np, SolverLU
def run(plotIt=True):
-6
View File
@@ -1,9 +1,3 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import *
def run(plotIt=True):
-6
View File
@@ -1,9 +1,3 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import *
def run(plotIt=True):
+2 -10
View File
@@ -1,11 +1,3 @@
from __future__ import print_function
from __future__ import unicode_literals
from __future__ import division
from __future__ import absolute_import
from builtins import int
from future import standard_library
standard_library.install_aliases()
from builtins import zip
from SimPEG import *
def run(plotIt=True, n=60):
@@ -95,7 +87,7 @@ def run(plotIt=True, n=60):
if elapsed > capture[jj]:
PHIS += [(elapsed, phi.copy())]
jj += 1
if ii % 10 == 0: print(ii, elapsed)
if ii % 10 == 0: print ii, elapsed
ii += 1
if plotIt:
@@ -106,7 +98,7 @@ def run(plotIt=True, n=60):
ii = int(ii)
out = M.plotImage(PHIS[ii][1],ax=ax)
ax.axis('off')
ax.set_title('Elapsed Time: %4.1f'%PHIS[ii][0])
ax.set_title('Elapsed Time: {0:4.1f}'.format(PHIS[ii][0]))
plt.show()
if __name__ == '__main__':
@@ -1,10 +1,3 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
from SimPEG import *
def run(plotIt=True):
+6 -14
View File
@@ -1,11 +1,3 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import zip
from builtins import range
from SimPEG import *
def run(plotIt=True, n=60):
@@ -36,16 +28,16 @@ def run(plotIt=True, n=60):
axes[0].set_xlim([-1,17])
axes[0].set_ylim([-1,17])
for ii, loc in zip(list(range(M.nC)),M.gridCC):
axes[0].text(loc[0]+0.2,loc[1],'%d'%ii, color='r')
for ii, loc in zip(range(M.nC),M.gridCC):
axes[0].text(loc[0]+0.2,loc[1],'{0:d}'.format(ii), color='r')
axes[0].plot(M.gridFx[:,0],M.gridFx[:,1], 'g>')
for ii, loc in zip(list(range(M.nFx)),M.gridFx):
axes[0].text(loc[0]+0.2,loc[1],'%d'%ii, color='g')
for ii, loc in zip(range(M.nFx),M.gridFx):
axes[0].text(loc[0]+0.2,loc[1],'{0:d}'.format(ii), color='g')
axes[0].plot(M.gridFy[:,0],M.gridFy[:,1], 'm^')
for ii, loc in zip(list(range(M.nFy)),M.gridFy):
axes[0].text(loc[0]+0.2,loc[1]+0.2,'%d'%(ii+M.nFx), color='m')
for ii, loc in zip(range(M.nFy),M.gridFy):
axes[0].text(loc[0]+0.2,loc[1]+0.2,'{0:d}'.format((ii+M.nFx)), color='m')
axes[1].spy(M.faceDiv)
axes[1].set_title('Face Divergence')
@@ -1,9 +1,3 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import *
def run(plotIt=True):
-6
View File
@@ -1,9 +1,3 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import *
def run(plotIt=True):
@@ -1,9 +1,3 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import *
from SimPEG.Utils import surface2ind_topo
+30 -37
View File
@@ -1,35 +1,28 @@
from __future__ import print_function
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import division
from builtins import open
from future import standard_library
standard_library.install_aliases()
# Run this file to add imports.
##### AUTOIMPORTS #####
from . import DC_Analytic_Dipole
from . import DC_Forward_PseudoSection
from . import EM_FDEM_1D_Inversion
from . import EM_FDEM_Analytic_MagDipoleWholespace
from . import EM_Schenkel_Morrison_Casing
from . import EM_TDEM_1D_Inversion
from . import FLOW_Richards_1D_Celia1990
from . import Inversion_IRLS
from . import Inversion_Linear
from . import Maps_ComboMaps
from . import Maps_Mesh2Mesh
from . import Mesh_Basic_ForwardDC
from . import Mesh_Basic_PlotImage
from . import Mesh_Basic_Types
from . import Mesh_Operators_CahnHilliard
from . import Mesh_QuadTree_Creation
from . import Mesh_QuadTree_FaceDiv
from . import Mesh_QuadTree_HangingNodes
from . import Mesh_Tensor_Creation
from . import MT_1D_ForwardAndInversion
from . import MT_3D_Foward
from . import Utils_surface2ind_topo
import DC_Analytic_Dipole
import DC_Forward_PseudoSection
import EM_FDEM_1D_Inversion
import EM_FDEM_Analytic_MagDipoleWholespace
import EM_Schenkel_Morrison_Casing
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
import Mesh_Operators_CahnHilliard
import Mesh_QuadTree_Creation
import Mesh_QuadTree_FaceDiv
import Mesh_QuadTree_HangingNodes
import Mesh_Tensor_Creation
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", "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"]
@@ -66,7 +59,7 @@ if __name__ == '__main__':
if line == "##### AUTOIMPORTS #####\n":
inimports = not inimports
if inimports:
out += '\n'.join(["import %s"%_ for _ in exfiles])
out += '\n'.join(["import {0!s}".format(_) for _ in exfiles])
out += '\n\n__examples__ = ["' + '", "'.join(exfiles)+ '"]\n'
out += '\n##### AUTOIMPORTS #####\n'
f.close()
@@ -83,11 +76,11 @@ if __name__ == '__main__':
docstr = runFunction.__doc__
if docstr is None:
doc = '%s\n%s'%(name.replace('_',' '),'='*len(name))
doc = '{0!s}\n{1!s}'.format(name.replace('_',' '), '='*len(name))
else:
doc = '\n'.join([_[8:].rstrip() for _ in docstr.split('\n')])
out = """.. _examples_%s:
out = """.. _examples_{0!s}:
.. --------------------------------- ..
.. ..
@@ -97,21 +90,21 @@ if __name__ == '__main__':
.. ..
.. --------------------------------- ..
%s
{1!s}
.. plot::
from SimPEG import Examples
Examples.%s.run()
Examples.{2!s}.run()
.. literalinclude:: ../../../SimPEG/Examples/%s.py
.. literalinclude:: ../../../SimPEG/Examples/{3!s}.py
:language: python
:linenos:
"""%(name,doc,name,name)
""".format(name, doc, name, name)
rst = os.path.sep.join((filePath.split(os.path.sep)[:-3] + ['docs', 'content', 'examples', name + '.rst']))
print('Creating: %s.rst'%name)
print 'Creating: {0!s}.rst'.format(name)
f = open(rst, 'w')
f.write(out)
f.close()
+7 -13
View File
@@ -1,20 +1,14 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import object
from SimPEG import Mesh, Maps, Utils, np
from future.utils import with_metaclass
class NonLinearMap(with_metaclass(Utils.SimPEGMetaClass, object)):
class NonLinearMap(object):
"""
SimPEG NonLinearMap
"""
__metaclass__ = Utils.SimPEGMetaClass
counter = None #: A SimPEG.Utils.Counter object
mesh = None #: A SimPEG Mesh
@@ -122,7 +116,7 @@ class RichardsMap(object):
ax.semilogx(self.k(h, m), h)
def _assertMatchesPair(self, pair):
assert isinstance(self, pair), "Mapping object must be an instance of a %s class."%(pair.__name__)
assert isinstance(self, pair), "Mapping object must be an instance of a {0!s} class.".format((pair.__name__))
@@ -192,7 +186,7 @@ class _haverkamp_theta(NonLinearMap):
def transformDerivU(self, u, m):
self.setModel(m)
g = (self.alpha*((self.theta_s - self.theta_r) /
g = (self.alpha*((self.theta_s - self.theta_r)/
(self.alpha + abs(u)**self.beta)**2)
*(-self.beta*abs(u)**(self.beta-1)*np.sign(u)))
g[u >= 0] = 0
@@ -279,7 +273,7 @@ class _vangenuchten_theta(NonLinearMap):
def transform(self, u, m):
self.setModel(m)
m = 1 - 1.0/self.n
f = (( self.theta_s - self.theta_r ) /
f = (( self.theta_s - self.theta_r )/
((1+abs(self.alpha*u)**self.n)**m) + self.theta_r)
if Utils.isScalar(self.theta_s):
f[u >= 0] = self.theta_s
@@ -349,7 +343,7 @@ class _vangenuchten_k(NonLinearMap):
Ks = self.Ks
m = 1.0 - 1.0/n
g = I*alpha*n*np.exp(Ks)*abs(alpha*u)**(n - 1.0)*np.sign(alpha*u)*(1.0/n - 1.0)*((abs(alpha*u)**n + 1)**(1.0/n - 1))**(I - 1)*((1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)**2*(abs(alpha*u)**n + 1)**(1.0/n - 2) - (2*alpha*n*np.exp(Ks)*abs(alpha*u)**(n - 1)*np.sign(alpha*u)*(1.0/n - 1)*((abs(alpha*u)**n + 1)**(1.0/n - 1))**I*((1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)*(abs(alpha*u)**n + 1)**(1.0/n - 2))/(((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1))+ 1)*(1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1.0/n)
g = I*alpha*n*np.exp(Ks)*abs(alpha*u)**(n - 1.0)*np.sign(alpha*u)*(1.0/n - 1.0)*((abs(alpha*u)**n + 1)**(1.0/n - 1))**(I - 1)*((1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)**2*(abs(alpha*u)**n + 1)**(1.0/n - 2) - (2*alpha*n*np.exp(Ks)*abs(alpha*u)**(n - 1)*np.sign(alpha*u)*(1.0/n - 1)*((abs(alpha*u)**n + 1)**(1.0/n - 1))**I*((1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)*(abs(alpha*u)**n + 1)**(1.0/n - 2))/(((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1) + 1)*(1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1.0/n))
g[u >= 0] = 0
g = Utils.sdiag(g)
return g
+8 -15
View File
@@ -1,12 +1,5 @@
from __future__ import print_function
from __future__ import absolute_import
from __future__ import division
from __future__ import unicode_literals
from future import standard_library
standard_library.install_aliases()
from builtins import range
from SimPEG import *
from .Empirical import RichardsMap
from Empirical import RichardsMap
import time
@@ -68,7 +61,7 @@ class RichardsSurvey(Survey.BaseSurvey):
@Utils.requires('prob')
def eval(self, U, m):
Ds = list(range(len(self.rxList)))
Ds = range(len(self.rxList))
for ii, rx in enumerate(self.rxList):
Ds[ii] = rx.eval(U, m,
self.prob.mapping,
@@ -80,7 +73,7 @@ class RichardsSurvey(Survey.BaseSurvey):
@Utils.requires('prob')
def evalDeriv(self, U, m):
"""The Derivative with respect to the fields."""
Ds = list(range(len(self.rxList)))
Ds = range(len(self.rxList))
for ii, rx in enumerate(self.rxList):
Ds[ii] = rx.evalDeriv(U, m,
self.prob.mapping,
@@ -142,12 +135,12 @@ class RichardsProblem(Problem.BaseTimeProblem):
@Utils.timeIt
def fields(self, m):
tic = time.time()
u = list(range(self.nT+1))
u = range(self.nT+1)
u[0] = self.initialConditions
for ii, dt in enumerate(self.timeSteps):
bc = self.getBoundaryConditions(ii, u[ii])
u[ii+1] = self.rootFinder.root(lambda hn1m, return_g=True: self.getResidual(m, u[ii], hn1m, dt, bc, return_g=return_g), u[ii])
if self.debug: print("Solving Fields (%4d/%d - %3.1f%% Done) %d Iterations, %4.2f seconds"%(ii+1, self.nT, 100.0*(ii+1)/self.nT, self.rootFinder.iter, time.time() - tic))
if self.debug: print "Solving Fields ({0:4d}/{1:d} - {2:3.1f}% Done) {3:d} Iterations, {4:4.2f} seconds".format(ii+1, self.nT, 100.0*(ii+1)/self.nT, self.rootFinder.iter, time.time() - tic)
return u
@Utils.timeIt
@@ -245,7 +238,7 @@ class RichardsProblem(Problem.BaseTimeProblem):
f = self.fields(m)
nn = len(f)-1
Asubs, Adiags, Bs = list(range(nn)), list(range(nn)), list(range(nn))
Asubs, Adiags, Bs = range(nn), range(nn), range(nn)
for ii in range(nn):
dt = self.timeSteps[ii]
bc = self.getBoundaryConditions(ii, f[ii])
@@ -270,7 +263,7 @@ class RichardsProblem(Problem.BaseTimeProblem):
if f is None:
f = self.fields(m)
JvC = list(range(len(f)-1)) # Cell to hold each row of the long vector.
JvC = range(len(f)-1) # Cell to hold each row of the long vector.
# This is done via forward substitution.
bc = self.getBoundaryConditions(0, f[0])
@@ -302,7 +295,7 @@ class RichardsProblem(Problem.BaseTimeProblem):
bc = self.getBoundaryConditions(ii-1, f[ii-1])
Asub, Adiag, B = self.diagsJacobian(m, f[ii-1], f[ii], self.timeSteps[ii-1], bc)
#select the correct part of v
vpart = list(range((ii)*Adiag.shape[0], (ii+1)*Adiag.shape[0]))
vpart = range((ii)*Adiag.shape[0], (ii+1)*Adiag.shape[0])
AdiaginvT = self.Solver(Adiag.T, **self.solverOpts)
JTvC = AdiaginvT * (PTv[vpart] - minus)
minus = Asub.T*JTvC # this is now the super diagonal.
+2 -8
View File
@@ -1,8 +1,2 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from . import Empirical
from .RichardsProblem import *
import Empirical
from RichardsProblem import *
+1 -7
View File
@@ -1,7 +1 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from . import Richards
import Richards
+6 -15
View File
@@ -1,13 +1,4 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
from builtins import object
from . import Utils
import numpy as np, scipy.sparse as sp
import Utils, numpy as np, scipy.sparse as sp
class Fields(object):
"""Fancy Field Storage
@@ -46,7 +37,7 @@ class Fields(object):
for f in self.knownFields:
loc =self.knownFields[f]
sz += np.array(self._storageShape(loc)).prod()*8.0/(1024**2)
return "%e MB"%sz
return "{0:e} MB".format(sz)
def _storageShape(self, loc):
nSrc = self.survey.nSrc
@@ -93,12 +84,12 @@ class Fields(object):
return
if accessType=='set' and name not in self.knownFields:
if name in self.aliasFields:
raise KeyError("Invalid field name (%s) for setter, you can't set an aliased property"%name)
raise KeyError("Invalid field name ({0!s}) for setter, you can't set an aliased property".format(name))
else:
raise KeyError('Invalid field name (%s) for setter'%name)
raise KeyError('Invalid field name ({0!s}) for setter'.format(name))
elif accessType=='get' and (name not in self.knownFields and name not in self.aliasFields):
raise KeyError('Invalid field name (%s) for getter'%name)
raise KeyError('Invalid field name ({0!s}) for getter'.format(name))
return name
def _indexAndNameFromKey(self, key, accessType):
@@ -253,7 +244,7 @@ class TimeFields(Fields):
out = func(pointerFields, srcII, timeII)
else: #loop over the time steps
nT = pointerShape[2]
out = list(range(nT))
out = range(nT)
for i, TIND_i in enumerate(timeII):
fieldI = pointerFields[:,:,i]
if fieldI.shape[0] == fieldI.size:
+12 -19
View File
@@ -1,21 +1,14 @@
from __future__ import print_function
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from builtins import object
from . import Utils, Survey, Problem
import numpy as np, scipy.sparse as sp, gc
from .Utils.SolverUtils import *
from . import DataMisfit
from . import Regularization
from future.utils import with_metaclass
import Utils, Survey, Problem, numpy as np, scipy.sparse as sp, gc
from Utils.SolverUtils import *
import DataMisfit
import Regularization
class BaseInvProblem(with_metaclass(Utils.SimPEGMetaClass, object)):
class BaseInvProblem(object):
"""BaseInvProblem(dmisfit, reg, opt)"""
__metaclass__ = Utils.SimPEGMetaClass
beta = 1.0 #: Trade-off parameter
debug = False #: Print debugging information
@@ -61,10 +54,10 @@ class BaseInvProblem(with_metaclass(Utils.SimPEGMetaClass, object)):
Called when inversion is first starting.
"""
if self.debug: print('Calling InvProblem.startup')
if self.debug: print 'Calling InvProblem.startup'
if self.reg.mref is None:
print('SimPEG.InvProblem will set Regularization.mref to m0.')
print 'SimPEG.InvProblem will set Regularization.mref to m0.'
self.reg.mref = m0
self.phi_d = np.nan
@@ -72,8 +65,8 @@ class BaseInvProblem(with_metaclass(Utils.SimPEGMetaClass, object)):
self.curModel = m0
print("""SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
***Done using same Solver and solverOpts as the problem***""")
print """SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
***Done using same Solver and solverOpts as the problem***"""
self.opt.bfgsH0 = self.prob.Solver(self.reg.eval2Deriv(self.curModel), **self.prob.solverOpts)
@property
@@ -94,7 +87,7 @@ class BaseInvProblem(with_metaclass(Utils.SimPEGMetaClass, object)):
for mtest, u_ofmtest in self.warmstart:
if m is mtest:
f = u_ofmtest
if self.debug: print('InvProb is Warm Starting!')
if self.debug: print 'InvProb is Warm Starting!'
break
if f is None:
+5 -11
View File
@@ -1,24 +1,18 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from builtins import object
import SimPEG
from SimPEG import Utils, sp, np
from .Optimization import Remember, IterationPrinters, StoppingCriteria
from . import Directives
from future.utils import with_metaclass
from Optimization import Remember, IterationPrinters, StoppingCriteria
import Directives
class BaseInversion(with_metaclass(Utils.SimPEGMetaClass, object)):
class BaseInversion(object):
"""
Inversion Class.
"""
__metaclass__ = Utils.SimPEGMetaClass
name = 'BaseInversion'
debug = False #: Print debugging information
+2 -8
View File
@@ -1,13 +1,7 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from SimPEG import SolverLU as SimpegSolver, PropMaps, Utils, mkvc, sp, np
from SimPEG.EM.FDEM.ProblemFDEM import BaseFDEMProblem
from .SurveyMT import Survey, Data
from .FieldsMT import BaseMTFields
from SurveyMT import Survey, Data
from FieldsMT import BaseMTFields
class BaseMTProblem(BaseFDEMProblem):
+6 -12
View File
@@ -1,9 +1,3 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import Survey, Utils, Problem, np, sp, mkvc
from scipy.constants import mu_0
import sys
@@ -69,7 +63,7 @@ class Fields1D_e(BaseMTFields):
C = self.mesh.nodalGrad
b = (C * eSolution)
for i, src in enumerate(srcList):
b[:,i] *= -1./(1j*omega(src.freq))
b[:,i] *= - 1./(1j*omega(src.freq))
# There is no magnetic source in the MT problem
# S_m, _ = src.eval(self.survey.prob)
# if S_m is not None:
@@ -194,7 +188,7 @@ class Fields3D_e(BaseMTFields):
# adjoint: returns a 2*nE long vector with zero's for py
return np.vstack((v,np.zeros_like(v)))
# Not adjoint: return only the px part of the vector
return v[:len(v)//2]
return v[:len(v)/2]
def _e_pyDeriv_u(self, src, v, adjoint = False):
'''
@@ -204,7 +198,7 @@ class Fields3D_e(BaseMTFields):
# adjoint: returns a 2*nE long vector with zero's for px
return np.vstack((np.zeros_like(v),v))
# Not adjoint: return only the px part of the vector
return v[len(v)//2::]
return v[len(v)/2::]
def _e_pxDeriv_m(self, src, v, adjoint = False):
# assuming primary does not depend on the model
@@ -233,7 +227,7 @@ class Fields3D_e(BaseMTFields):
C = self.mesh.edgeCurl
b = (C * e_pxSolution)
for i, src in enumerate(srcList):
b[:,i] *= -1./(1j*omega(src.freq))
b[:,i] *= - 1./(1j*omega(src.freq))
# There is no magnetic source in the MT problem
# S_m, _ = src.eval(self.survey.prob)
# if S_m is not None:
@@ -244,7 +238,7 @@ class Fields3D_e(BaseMTFields):
C = self.mesh.edgeCurl
b = (C * e_pySolution)
for i, src in enumerate(srcList):
b[:,i] *= -1./(1j*omega(src.freq))
b[:,i] *= - 1./(1j*omega(src.freq))
# There is no magnetic source in the MT problem
# S_m, _ = src.eval(self.survey.prob)
# if S_m is not None:
@@ -354,4 +348,4 @@ class Fields3D_e(BaseMTFields):
This function stacks the fields derivatives appropriately
"""
# The fields have no dependance to the model.
return None
return None
+7 -13
View File
@@ -1,9 +1,3 @@
from __future__ import print_function
from __future__ import division
from __future__ import unicode_literals
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG.EM.Utils import omega
from SimPEG import mkvc
from scipy.constants import mu_0
@@ -52,7 +46,7 @@ class eForm_psField(BaseMTProblem):
Edge inner product matrix
"""
if getattr(self, '_MeMui', None) is None:
self._MeMui = self.mesh.getEdgeInnerProduct(old_div(1.0,mu_0))
self._MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
return self._MeMui
@property
@@ -148,7 +142,7 @@ class eForm_psField(BaseMTProblem):
for freq in self.survey.freqs:
if self.verbose:
startTime = time.time()
print('Starting work for {:.3e}'.format(freq))
print 'Starting work for {:.3e}'.format(freq)
sys.stdout.flush()
A = self.getA(freq)
rhs = self.getRHS(freq)
@@ -164,7 +158,7 @@ class eForm_psField(BaseMTProblem):
# b = -( self.mesh.nodalGrad * e )/( 1j*omega(freq) )
# F[Src, 'b_1d'] = b[:,1]
if self.verbose:
print('Ran for {:f} seconds'.format(time.time()-startTime))
print 'Ran for {:f} seconds'.format(time.time()-startTime)
sys.stdout.flush()
return F
@@ -197,7 +191,7 @@ class eForm_TotalField(BaseMTProblem):
Edge inner product matrix
"""
if getattr(self, '_MeMui', None) is None:
self._MeMui = self.mesh.getEdgeInnerProduct(old_div(1.0,mu_0))
self._MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
return self._MeMui
@property
@@ -255,7 +249,7 @@ class eForm_TotalField(BaseMTProblem):
Ed, Eu, Hd, Hu = getEHfields(self.mesh,self.curModel.sigma,freq,self.mesh.vectorNx)
Etot = (Ed + Eu)
sourceAmp = 1.0
Etot = ((old_div(Etot,Etot[-1]))*sourceAmp) # Scale the fields to be equal to sourceAmp at the top
Etot = ((Etot/Etot[-1])*sourceAmp) # Scale the fields to be equal to sourceAmp at the top
## Note: The analytic solution is derived with e^iwt
eBC = np.r_[Etot[0],Etot[-1]]
# The right hand side
@@ -280,7 +274,7 @@ class eForm_TotalField(BaseMTProblem):
for freq in self.survey.freqs:
if self.verbose:
startTime = time.time()
print('Starting work for {:.3e}'.format(freq))
print 'Starting work for {:.3e}'.format(freq)
sys.stdout.flush()
A = self.getA(freq)
rhs, e_o = self.getRHS(freq)
@@ -292,6 +286,6 @@ class eForm_TotalField(BaseMTProblem):
# NOTE: only store e fields
F[Src, 'e_1dSolution'] = e[:,0]
if self.verbose:
print('Ran for {:f} seconds'.format(time.time()-startTime))
print 'Ran for {:f} seconds'.format(time.time()-startTime)
sys.stdout.flush()
return F
+1 -7
View File
@@ -1,7 +1 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from .Probs import eForm_TotalField, eForm_psField
from Probs import eForm_TotalField, eForm_psField
-6
View File
@@ -1,7 +1 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
pass
+2 -8
View File
@@ -1,9 +1,3 @@
from __future__ import print_function
from __future__ import unicode_literals
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import Survey, Problem, Utils, Models, np, sp, mkvc, SolverLU as SimpegSolver
from SimPEG.EM.Utils import omega
from scipy.constants import mu_0
@@ -121,7 +115,7 @@ class eForm_ps(BaseMTProblem):
for freq in self.survey.freqs:
if self.verbose:
startTime = time.time()
print('Starting work for {:.3e}'.format(freq))
print 'Starting work for {:.3e}'.format(freq)
sys.stdout.flush()
A = self.getA(freq)
rhs = self.getRHS(freq)
@@ -137,7 +131,7 @@ class eForm_ps(BaseMTProblem):
# Note curl e = -iwb so b = -curl/iw
if self.verbose:
print('Ran for {:f} seconds'.format(time.time()-startTime))
print 'Ran for {:f} seconds'.format(time.time()-startTime)
sys.stdout.flush()
Ainv.clean()
return F
+1 -7
View File
@@ -1,7 +1 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from .Probs import eForm_ps
from Probs import eForm_ps
+4 -10
View File
@@ -1,16 +1,10 @@
from __future__ import absolute_import
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from future import standard_library
standard_library.install_aliases()
from SimPEG import Utils, Problem, Maps, np, sp, mkvc
from SimPEG.EM.FDEM.SrcFDEM import BaseSrc as FDEMBaseSrc
from SimPEG.EM.Utils import omega
from scipy.constants import mu_0
from numpy.lib import recfunctions as recFunc
from .Utils.sourceUtils import homo1DModelSource
from .Utils import rec2ndarr
from Utils.sourceUtils import homo1DModelSource
from Utils import rec2ndarr
import sys
#################
@@ -84,7 +78,7 @@ class polxy_1Dprimary(BaseMTSrc):
C = problem.mesh.nodalGrad
elif problem.mesh.dim == 3:
C = problem.mesh.edgeCurl
bBG_bp = (- C * self.ePrimary(problem) )*(1/(1j*omega(self.freq)))
bBG_bp = (- C * self.ePrimary(problem) )*(1/( 1j*omega(self.freq) ))
return bBG_bp
def S_e(self,problem):
@@ -161,7 +155,7 @@ class polxy_3Dprimary(BaseMTSrc):
C = problem.mesh.nodalGrad
elif problem.mesh.dim == 3:
C = problem.mesh.edgeCurl
bBG_bp = (- C * self.ePrimary(problem) )*(1/(1j*omega(self.freq)))
bBG_bp = (- C * self.ePrimary(problem) )*(1/( 1j*omega(self.freq) ))
return bBG_bp
def S_e(self,problem):
+21 -27
View File
@@ -1,16 +1,10 @@
from __future__ import absolute_import
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from future import standard_library
standard_library.install_aliases()
from SimPEG import Survey as SimPEGsurvey, Utils, Problem, Maps, np, sp, mkvc
from SimPEG.EM.FDEM.SrcFDEM import BaseSrc as FDEMBaseSrc
from SimPEG.EM.Utils import omega
from scipy.constants import mu_0
from numpy.lib import recfunctions as recFunc
from .Utils import rec2ndarr
from . import SrcMT
from Utils import rec2ndarr
import SrcMT
import sys
#################
@@ -82,7 +76,7 @@ class Rx(SimPEGsurvey.BaseRx):
bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0
# Note: Has a minus sign in front, to comply with quadrant calculations.
# Can be derived from zyx case for the 3D case.
f_part_complex = old_div(-ex,bx)
f_part_complex = -ex/bx
# elif self.projType is 'Z2D':
elif self.projType is 'Z3D':
## NOTE: Assumes that e is on edges and b on the faces. Need to generalize that or use a prop of fields to determine that.
@@ -109,13 +103,13 @@ class Rx(SimPEGsurvey.BaseRx):
hy_py = Pby*f[src,'b_py']/mu_0
# Make the complex data
if 'zxx' in self.rxType:
f_part_complex = old_div(( ex_px*hy_py - ex_py*hy_px),(hx_px*hy_py - hx_py*hy_px))
f_part_complex = ( ex_px*hy_py - ex_py*hy_px)/(hx_px*hy_py - hx_py*hy_px)
elif 'zxy' in self.rxType:
f_part_complex = old_div((-ex_px*hx_py + ex_py*hx_px),(hx_px*hy_py - hx_py*hy_px))
f_part_complex = (-ex_px*hx_py + ex_py*hx_px)/(hx_px*hy_py - hx_py*hy_px)
elif 'zyx' in self.rxType:
f_part_complex = old_div(( ey_px*hy_py - ey_py*hy_px),(hx_px*hy_py - hx_py*hy_px))
f_part_complex = ( ey_px*hy_py - ey_py*hy_px)/(hx_px*hy_py - hx_py*hy_px)
elif 'zyy' in self.rxType:
f_part_complex = old_div((-ey_px*hx_py + ey_py*hx_px),(hx_px*hy_py - hx_py*hy_px))
f_part_complex = (-ey_px*hx_py + ey_py*hx_px)/(hx_px*hy_py - hx_py*hy_px)
elif self.projType is 'T3D':
if self.locs.ndim == 3:
horLoc = self.locs[:,:,0]
@@ -133,9 +127,9 @@ class Rx(SimPEGsurvey.BaseRx):
by_py = Pby*f[src,'b_py']
bz_py = Pbz*f[src,'b_py']
if 'tzx' in self.rxType:
f_part_complex = old_div((- by_px*bz_py + by_py*bz_px),(bx_px*by_py - bx_py*by_px))
f_part_complex = (- by_px*bz_py + by_py*bz_px)/(bx_px*by_py - bx_py*by_px)
if 'tzy' in self.rxType:
f_part_complex = old_div(( bx_px*bz_py - bx_py*bz_px),(bx_px*by_py - bx_py*by_px))
f_part_complex = ( bx_px*bz_py - bx_py*bz_px)/(bx_px*by_py - bx_py*by_px)
else:
NotImplementedError('Projection of {:s} receiver type is not implemented.'.format(self.rxType))
@@ -163,8 +157,8 @@ class Rx(SimPEGsurvey.BaseRx):
Pbx = mesh.getInterpolationMat(self.locs[:,-1],'Ex')
# ex = Pex*mkvc(f[src,'e_1d'],2)
# bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0
dP_de = -mkvc(Utils.sdiag(old_div(1.,(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)))*(Pex*v),2)
dP_db = mkvc( Utils.sdiag(Pex*mkvc(f[src,'e_1d'],2))*(Utils.sdiag(old_div(1.,(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))).T*Utils.sdiag(old_div(1.,(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))))*(Pbx*f._bDeriv_u(src,v)/mu_0),2)
dP_de = -mkvc(Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))*(Pex*v),2)
dP_db = mkvc( Utils.sdiag(Pex*mkvc(f[src,'e_1d'],2))*(Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)).T*Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)))*(Pbx*f._bDeriv_u(src,v)/mu_0),2)
PDeriv_complex = np.sum(np.hstack((dP_de,dP_db)),1)
elif self.projType is 'Z2D':
raise NotImplementedError('Has not been implement for 2D impedance tensor')
@@ -204,7 +198,7 @@ class Rx(SimPEGsurvey.BaseRx):
# Update the input vector
sDiag = lambda t: Utils.sdiag(mkvc(t,2))
# Define the components of the derivative
Hd = sDiag(old_div(1.,(sDiag(hx_px)*hy_py - sDiag(hx_py)*hy_px)))
Hd = sDiag(1./(sDiag(hx_px)*hy_py - sDiag(hx_py)*hy_px))
Hd_uV = sDiag(hy_py)*hx_px_u(v) + sDiag(hx_px)*hy_py_u(v) - sDiag(hx_py)*hy_px_u(v) - sDiag(hy_px)*hx_py_u(v)
# Calculate components
if 'zxx' in self.rxType:
@@ -253,7 +247,7 @@ class Rx(SimPEGsurvey.BaseRx):
# Update the input vector
sDiag = lambda t: Utils.sdiag(mkvc(t,2))
# Define the components of the derivative
Hd = sDiag(old_div(1.,(sDiag(bx_px)*by_py - sDiag(bx_py)*by_px)))
Hd = sDiag(1./(sDiag(bx_px)*by_py - sDiag(bx_py)*by_px))
Hd_uV = sDiag(by_py)*bx_px_u(v) + sDiag(bx_px)*by_py_u(v) - sDiag(bx_py)*by_px_u(v) - sDiag(by_px)*bx_py_u(v)
if 'tzx' in self.rxType:
Tij = sDiag(Hd*( - sDiag(by_px)*bz_py + sDiag(by_py)*bz_px ))
@@ -273,8 +267,8 @@ class Rx(SimPEGsurvey.BaseRx):
Pbx = mesh.getInterpolationMat(self.locs[:,-1],'Ex')
# ex = Pex*mkvc(f[src,'e_1d'],2)
# bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0
dP_deTv = -mkvc(Pex.T*Utils.sdiag(old_div(1.,(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))).T*v,2)
db_duv = Pbx.T/mu_0*Utils.sdiag(old_div(1.,(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)))*(Utils.sdiag(old_div(1.,(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)))).T*Utils.sdiag(Pex*mkvc(f[src,'e_1d'],2)).T*v
dP_deTv = -mkvc(Pex.T*Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)).T*v,2)
db_duv = Pbx.T/mu_0*Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))*(Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))).T*Utils.sdiag(Pex*mkvc(f[src,'e_1d'],2)).T*v
dP_dbTv = mkvc(f._bDeriv_u(src,db_duv,adjoint=True),2)
PDeriv_real = np.sum(np.hstack((dP_deTv,dP_dbTv)),1)
elif self.projType is 'Z2D':
@@ -306,17 +300,17 @@ class Rx(SimPEGsurvey.BaseRx):
aey_px_u = lambda vec: f._e_pxDeriv_u(src,Pey.T*vec,adjoint=True)
aex_py_u = lambda vec: f._e_pyDeriv_u(src,Pex.T*vec,adjoint=True)
aey_py_u = lambda vec: f._e_pyDeriv_u(src,Pey.T*vec,adjoint=True)
ahx_px_u = lambda vec: old_div(f._b_pxDeriv_u(src,Pbx.T*vec,adjoint=True),mu_0)
ahy_px_u = lambda vec: old_div(f._b_pxDeriv_u(src,Pby.T*vec,adjoint=True),mu_0)
ahx_py_u = lambda vec: old_div(f._b_pyDeriv_u(src,Pbx.T*vec,adjoint=True),mu_0)
ahy_py_u = lambda vec: old_div(f._b_pyDeriv_u(src,Pby.T*vec,adjoint=True),mu_0)
ahx_px_u = lambda vec: f._b_pxDeriv_u(src,Pbx.T*vec,adjoint=True)/mu_0
ahy_px_u = lambda vec: f._b_pxDeriv_u(src,Pby.T*vec,adjoint=True)/mu_0
ahx_py_u = lambda vec: f._b_pyDeriv_u(src,Pbx.T*vec,adjoint=True)/mu_0
ahy_py_u = lambda vec: f._b_pyDeriv_u(src,Pby.T*vec,adjoint=True)/mu_0
# Update the input vector
# Define shortcuts
sDiag = lambda t: Utils.sdiag(mkvc(t,2))
sVec = lambda t: Utils.sp.csr_matrix(mkvc(t,2))
# Define the components of the derivative
aHd = sDiag(old_div(1.,(sDiag(ahx_px)*ahy_py - sDiag(ahx_py)*ahy_px)))
aHd = sDiag(1./(sDiag(ahx_px)*ahy_py - sDiag(ahx_py)*ahy_px))
aHd_uV = lambda x: ahx_px_u(sDiag(ahy_py)*x) + ahx_px_u(sDiag(ahy_py)*x) - ahy_px_u(sDiag(ahx_py)*x) - ahx_py_u(sDiag(ahy_px)*x)
# Need to fix this to reflect the adjoint
if 'zxx' in self.rxType:
@@ -368,7 +362,7 @@ class Rx(SimPEGsurvey.BaseRx):
sDiag = lambda t: Utils.sdiag(mkvc(t,2))
sVec = lambda t: Utils.sp.csr_matrix(mkvc(t,2))
# Define the components of the derivative
aHd = sDiag(old_div(1.,(sDiag(abx_px)*aby_py - sDiag(abx_py)*aby_px)))
aHd = sDiag(1./(sDiag(abx_px)*aby_py - sDiag(abx_py)*aby_px))
aHd_uV = lambda x: abx_px_u(sDiag(aby_py)*x) + abx_px_u(sDiag(aby_py)*x) - aby_px_u(sDiag(abx_py)*x) - abx_py_u(sDiag(aby_px)*x)
# Need to fix this to reflect the adjoint
if 'tzx' in self.rxType:
+10 -17
View File
@@ -1,10 +1,3 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import zip
# Analytic solution of EM fields due to a plane wave
import numpy as np, SimPEG as simpeg
@@ -40,8 +33,8 @@ def getEHfields(m1d,sigma,freq,zd,scaleUD=True):
# Loop over all the layers, starting at the bottom layer
for lnr, h in enumerate(m1d.hx): # lnr-number of layer, h-thickness of the layer
# Calculate
yp1 = old_div(k[lnr],(w*mu[lnr])) # Admittance of the layer below the current layer
zp = old_div((w*mu[lnr+1]),k[lnr+1]) # Impedance in the current layer
yp1 = k[lnr]/(w*mu[lnr]) # Admittance of the layer below the current layer
zp = (w*mu[lnr+1])/k[lnr+1] # Impedance in the current layer
# Build the propagation matrix
# Convert fields to down/up going components in layer below current layer
@@ -55,7 +48,7 @@ def getEHfields(m1d,sigma,freq,zd,scaleUD=True):
UDp[:,lnr+1] = elamh.dot(Pjinv.dot(Pj1)).dot(UDp[:,lnr])
if scaleUD:
UDp[:,lnr+1::-1] = old_div(UDp[:,lnr+1::-1],UDp[1,lnr+1])
UDp[:,lnr+1::-1] = UDp[:,lnr+1::-1]/UDp[1,lnr+1]
# Calculate the fields
Ed = np.empty((zd.size,),dtype=complex)
@@ -69,14 +62,14 @@ def getEHfields(m1d,sigma,freq,zd,scaleUD=True):
dind = dup >= zd
Ed[dind] = UDp[1,0]*np.exp(-1j*k[0]*(dup-zd[dind]))
Eu[dind] = UDp[0,0]*np.exp(1j*k[0]*(dup-zd[dind]))
Hd[dind] = (old_div(k[0],(w*mu[0])))*UDp[1,0]*np.exp(-1j*k[0]*(dup-zd[dind]))
Hu[dind] = -(old_div(k[0],(w*mu[0])))*UDp[0,0]*np.exp(1j*k[0]*(dup-zd[dind]))
Hd[dind] = (k[0]/(w*mu[0]))*UDp[1,0]*np.exp(-1j*k[0]*(dup-zd[dind]))
Hu[dind] = -(k[0]/(w*mu[0]))*UDp[0,0]*np.exp(1j*k[0]*(dup-zd[dind]))
for ki,mui,epsi,dlow,dup,Up,Dp in zip(k[1::],mu[1::],eps[1::],m1d.vectorNx[:-1],m1d.vectorNx[1::],UDp[0,1::],UDp[1,1::]):
dind = np.logical_and(dup >= zd, zd > dlow)
Ed[dind] = Dp*np.exp(-1j*ki*(dup-zd[dind]))
Eu[dind] = Up*np.exp(1j*ki*(dup-zd[dind]))
Hd[dind] = (old_div(ki,(w*mui)))*Dp*np.exp(-1j*ki*(dup-zd[dind]))
Hu[dind] = -(old_div(ki,(w*mui)))*Up*np.exp(1j*ki*(dup-zd[dind]))
Hd[dind] = (ki/(w*mui))*Dp*np.exp(-1j*ki*(dup-zd[dind]))
Hu[dind] = -(ki/(w*mui))*Up*np.exp(1j*ki*(dup-zd[dind]))
# Return return the fields
return Ed, Eu, Hd, Hu
@@ -99,15 +92,15 @@ def getImpedance(m1d,sigma,freq):
om = 2*np.pi*fr
Zall = np.empty(len(h)+1,dtype='complex')
# Calculate the impedance for the bottom layer
Zall[0] = old_div((mu_0*om),np.sqrt(mu_0*eps_0*(om)**2 - 1j*mu_0*sigma[0]*om))
Zall[0] = (mu_0*om)/np.sqrt(mu_0*eps_0*(om)**2 - 1j*mu_0*sigma[0]*om)
for nr,hi in enumerate(h):
# Calculate the wave number
# print nr,sigma[nr]
k = np.sqrt(mu_0*eps_0*om**2 - 1j*mu_0*sigma[nr]*om)
Z = old_div((mu_0*om),k)
Z = (mu_0*om)/k
Zall[nr+1] = Z *(old_div((Zall[nr] + Z*np.tanh(1j*k*hi)),(Z + Zall[nr]*np.tanh(1j*k*hi))))
Zall[nr+1] = Z *((Zall[nr] + Z*np.tanh(1j*k*hi))/(Z + Zall[nr]*np.tanh(1j*k*hi)))
#pdb.set_trace()
Z1d[nrFr] = Zall[-1]
+3 -9
View File
@@ -1,11 +1,5 @@
from __future__ import absolute_import
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from future import standard_library
standard_library.install_aliases()
import numpy as np, SimPEG as simpeg
from .MT1Danalytic import getEHfields
from MT1Danalytic import getEHfields
from scipy.constants import mu_0
def get1DEfields(m1d,sigma,freq,sourceAmp=1.0):
@@ -15,7 +9,7 @@ def get1DEfields(m1d,sigma,freq,sourceAmp=1.0):
G = m1d.nodalGrad
# Mass matrices
# Magnetic permeability
Mmu = simpeg.Utils.sdiag(m1d.vol*(old_div(1.0,mu_0)))
Mmu = simpeg.Utils.sdiag(m1d.vol*(1.0/mu_0))
# Conductivity
Msig = m1d.getFaceInnerProduct(sigma)
# Set up the solution matrix
@@ -29,7 +23,7 @@ def get1DEfields(m1d,sigma,freq,sourceAmp=1.0):
Ed, Eu, Hd, Hu = getEHfields(m1d,sigma,freq,m1d.vectorNx)
Etot = (Ed + Eu)
if sourceAmp is not None:
Etot = ((old_div(Etot,Etot[-1]))*sourceAmp) # Scale the fields to be equal to sourceAmp at the top
Etot = ((Etot/Etot[-1])*sourceAmp) # Scale the fields to be equal to sourceAmp at the top
## Note: The analytic solution is derived with e^iwt
bc = np.r_[Etot[0],Etot[-1]]
# The right hand side
+4 -10
View File
@@ -1,10 +1,4 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from .MT1Dsolutions import * # Add the names of the functions
from .MT1Danalytic import *
from .dataUtils import *
from .ediFilesUtils import *
from MT1Dsolutions import * # Add the names of the functions
from MT1Danalytic import *
from dataUtils import *
from ediFilesUtils import *
+9 -15
View File
@@ -1,9 +1,3 @@
from __future__ import print_function
from __future__ import absolute_import
from __future__ import division
from __future__ import unicode_literals
from future import standard_library
standard_library.install_aliases()
# Utils used for the data,
import numpy as np, matplotlib.pyplot as plt, sys
import SimPEG as simpeg
@@ -51,13 +45,13 @@ def rotateData(MTdata, rotAngle):
def appResPhs(freq, z):
app_res = (old_div((old_div(1.,(8e-7*np.pi**2))),freq))*np.abs(z)**2
app_phs = np.arctan2(z.imag,z.real)*(old_div(180,np.pi))
app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
app_phs = np.arctan2(z.imag,z.real)*(180/np.pi)
return app_res, app_phs
def skindepth(rho, freq):
''' Function to calculate the skindepth of EM waves'''
return np.sqrt( (rho*((old_div(1,(freq * mu_0 * np.pi ))))))
return np.sqrt( (rho*((1/(freq * mu_0 * np.pi )))))
def rec2ndarr(x, dt=float):
return x.view((dt, len(x.dtype.names)))
@@ -70,7 +64,7 @@ def makeAnalyticSolution(mesh, model, elev, freqs):
anaE = anaEd+anaEu
anaH = anaHd+anaHu
anaZ = old_div(anaE,anaH)
anaZ = anaE/anaH
# Add to the list
data1D.append((freq,0,0,elev,anaZ[0]))
dataRec = np.array(data1D,dtype=[('freq',float),('x',float),('y',float),('z',float),('zyx',complex)])
@@ -103,7 +97,7 @@ def plotMT1DModelData(problem, models, symList=None):
# if not symList:
# symList = ['x']*len(models)
from . import plotDataTypes as pDt
import plotDataTypes as pDt
# Loop through the models.
modelList = [problem.survey.mtrue]
modelList.extend(models)
@@ -116,14 +110,14 @@ def plotMT1DModelData(problem, models, symList=None):
else:
data1D = problem.dataPair(problem.survey,problem.survey.dpred(model)).toRecArray('Complex')
# Plot the data and the model
colRat = old_div(nr,((len(modelList)-1.999)*1.))
colRat = nr/((len(modelList)-1.999)*1.)
if colRat > 1.:
col = 'k'
else:
col = plt.cm.seismic(1-colRat)
# The model - make the pts to plot
meshPts = np.concatenate((problem.mesh.gridN[0:1],np.kron(problem.mesh.gridN[1::],np.ones(2))[:-1]))
modelPts = np.kron(old_div(1.,(problem.mapping.sigmaMap*model)),np.ones(2,))
modelPts = np.kron(1./(problem.mapping.sigmaMap*model),np.ones(2,))
axM.semilogx(modelPts,meshPts,color=col)
## Data
@@ -150,7 +144,7 @@ def plotMT1DModelData(problem, models, symList=None):
# Fix labels and ticks
yMtick = [old_div(l,1000) for l in axM.get_yticks().tolist()]
yMtick = [l/1000 for l in axM.get_yticks().tolist()]
axM.set_yticklabels(yMtick)
[ l.set_rotation(90) for l in axM.get_yticklabels()]
[ l.set_rotation(90) for l in axR.get_yticklabels()]
@@ -163,7 +157,7 @@ def plotMT1DModelData(problem, models, symList=None):
def printTime():
import time
print(time.strftime("%a, %d %b %Y %H:%M:%S +0000", time.localtime()))
print time.strftime("%a, %d %b %Y %H:%M:%S +0000", time.localtime())
def convert3Dto1Dobject(MTdata,rxType3D='zyx'):
from SimPEG import MT
+4 -13
View File
@@ -1,12 +1,3 @@
from __future__ import print_function
from __future__ import division
from __future__ import unicode_literals
from __future__ import absolute_import
from builtins import open
from builtins import int
from future import standard_library
standard_library.install_aliases()
from builtins import object
# Functions to import and export MT EDI files.
from SimPEG import mkvc
from scipy.constants import mu_0
@@ -18,7 +9,7 @@ import numpy as np
import os, sys, re
class EDIimporter(object):
class EDIimporter:
"""
A class to import EDIfiles.
@@ -27,7 +18,7 @@ class EDIimporter(object):
# Define data converters
_impUnitEDI2SI = 4*np.pi*1e-4 # Convert Z[mV/km/nT] (as in EDI)to Z[V/A] SI unit
_impUnitSI2EDI = old_div(1.,_impUnitEDI2SI) # ConvertZ[V/A] SI unit to Z[mV/km/nT] (as in EDI)
_impUnitSI2EDI = 1./_impUnitEDI2SI # ConvertZ[V/A] SI unit to Z[mV/km/nT] (as in EDI)
# Properties
filesList = None
@@ -125,7 +116,7 @@ class EDIimporter(object):
try:
import osr
except ImportError as e:
print('Could not import osr, missing the gdal package\nCan not project coordinates')
print 'Could not import osr, missing the gdal package\nCan not project coordinates'
raise e
# Coordinates convertor
if self._2out is None:
@@ -135,7 +126,7 @@ class EDIimporter(object):
if self._outEPSG is None:
# Find the UTM EPSG number
Nnr = 700 if latD < 0.0 else 600
utmZ = int(1+old_div((longD+180.0),6.0))
utmZ = int(1+(longD+180.0)/6.0)
self._outEPSG = 32000 + Nnr + utmZ
out.ImportFromEPSG(self._outEPSG)
self._2out = osr.CoordinateTransformation(src,out)
+49 -55
View File
@@ -1,9 +1,3 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from matplotlib import pyplot as plt, colors, numpy as np
@@ -41,7 +35,7 @@ def plotIsoFreqNSimpedance(ax,freq,array,flag,par='abs',colorbar=True,colorNorm=
cmap = plt.get_cmap('RdYlBu')
level = np.concatenate((-np.logspace(0,-10,31),np.logspace(-10,0,31)))
clevel = np.concatenate((-np.logspace(0,-8,5),np.logspace(-8,0,5)))
plotNorm = colors.SymLogNorm(1e-10,linscale=2)
plotNorm = colors.SymLogNorm(1e-10,linscale=2)
if cLevel:
level = np.concatenate((-np.logspace(0,-10,31),np.logspace(-10,0,31)))
clevel = np.concatenate((-np.logspace(0,-8,5),np.logspace(-8,0,5)))
@@ -71,45 +65,45 @@ def plotIsoFreqNSDiff(ax,freq,arrayList,flag,par='abs',colorbar=True,cLevel=True
x, y = arrayList[0]['x'][indUniFreq0],arrayList[0]['y'][indUniFreq0]
if par == 'abs':
if useLog:
zPlot = old_div((np.log10(np.abs(arrayList[0][flag][indUniFreq0])) - np.log10(np.abs(arrayList[1][flag][indUniFreq1]))),np.log10(np.abs(arrayList[1][flag][indUniFreq1])))
zPlot = (np.log10(np.abs(arrayList[0][flag][indUniFreq0])) - np.log10(np.abs(arrayList[1][flag][indUniFreq1])))/np.log10(np.abs(arrayList[1][flag][indUniFreq1]))
else:
zPlot = old_div((np.abs(arrayList[0][flag][indUniFreq0]) - np.abs(arrayList[1][flag][indUniFreq1])),np.abs(arrayList[1][flag][indUniFreq1]))
zPlot = (np.abs(arrayList[0][flag][indUniFreq0]) - np.abs(arrayList[1][flag][indUniFreq1]))/np.abs(arrayList[1][flag][indUniFreq1])
if mask:
maskInd = np.logical_or(np.abs(arrayList[0][flag][indUniFreq0])< 1e-3,np.abs(arrayList[1][flag][indUniFreq1]) < 1e-3)
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
if cLevel:
level = np.arange(-200,201,10)
level = np.arange(-200,201,10)
clevel = np.arange(-200,201,25)
else:
level = np.linspace(zPlot.min(),zPlot.max(),100)
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
elif par == 'real':
if useLog:
zPlot = old_div((np.log10(np.real(arrayList[0][flag][indUniFreq0])) -np.log10(np.real(arrayList[1][flag][indUniFreq1]))),np.log10(np.abs((np.real(arrayList[1][flag][indUniFreq1])))))
zPlot = (np.log10(np.real(arrayList[0][flag][indUniFreq0])) -np.log10(np.real(arrayList[1][flag][indUniFreq1])))/np.log10(np.abs((np.real(arrayList[1][flag][indUniFreq1]))))
else:
zPlot = old_div((np.real(arrayList[0][flag][indUniFreq0]) -np.real(arrayList[1][flag][indUniFreq1])),np.abs((np.real(arrayList[1][flag][indUniFreq1]))))
zPlot = (np.real(arrayList[0][flag][indUniFreq0]) -np.real(arrayList[1][flag][indUniFreq1]))/np.abs((np.real(arrayList[1][flag][indUniFreq1])))
if mask:
maskInd = np.logical_or(np.abs(np.real(arrayList[0][flag][indUniFreq0])) < 1e-3,np.abs(np.real(arrayList[1][flag][indUniFreq1])) < 1e-3)
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
if cLevel:
level = np.arange(-200,201,10)
level = np.arange(-200,201,10)
clevel = np.arange(-200,201,25)
else:
level = np.linspace(zPlot.min(),zPlot.max(),100)
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
elif par == 'imag':
if useLog:
zPlot = old_div((np.log10(np.imag(arrayList[0][flag][indUniFreq0])) -np.log10(np.imag(arrayList[1][flag][indUniFreq1]))),np.log10(np.abs((np.imag(arrayList[1][flag][indUniFreq1])))))
zPlot = (np.log10(np.imag(arrayList[0][flag][indUniFreq0])) -np.log10(np.imag(arrayList[1][flag][indUniFreq1])))/np.log10(np.abs((np.imag(arrayList[1][flag][indUniFreq1]))))
else:
zPlot = old_div((np.imag(arrayList[0][flag][indUniFreq0]) -np.imag(arrayList[1][flag][indUniFreq1])),np.abs((np.imag(arrayList[1][flag][indUniFreq1]))))
zPlot = (np.imag(arrayList[0][flag][indUniFreq0]) -np.imag(arrayList[1][flag][indUniFreq1]))/np.abs((np.imag(arrayList[1][flag][indUniFreq1])))
if mask:
maskInd = np.logical_or(np.abs(np.imag(arrayList[0][flag][indUniFreq0])) < 1e-3,np.abs(np.imag(arrayList[1][flag][indUniFreq1])) < 1e-3)
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
if cLevel:
level = np.arange(-200,201,10)
level = np.arange(-200,201,10)
clevel = np.arange(-200,201,25)
else:
level = np.linspace(zPlot.min(),zPlot.max(),100)
@@ -180,13 +174,13 @@ def plotIsoFreqNStipper(ax,freq,array,flag,par='abs',colorbar=True,colorNorm='Sy
plt.colorbar(cs,cax=ax.cax,ticks=clevel,format='%1.2e')
ax.set_title(flag+' '+par,fontsize=8)
def plotIsoStaImpedance(ax,loc,array,flag,par='abs',pSym='s',pColor=None):
appResFact = old_div(1,(8*np.pi**2*10**(-7)))
def plotIsoStaImpedance(ax,loc,array,flag,par='abs',pSym='s',pColor=None):
appResFact = 1/(8*np.pi**2*10**(-7))
treshold = 1.0 # 1 meter
indUniSta = np.sqrt(np.sum((rec2nd(array[['x','y']])-loc)**2,axis=1)) < treshold
freq = array['freq'][indUniSta]
if par == 'abs':
zPlot = np.abs(array[flag][indUniSta])
elif par == 'real':
@@ -194,9 +188,9 @@ def plotIsoStaImpedance(ax,loc,array,flag,par='abs',pSym='s',pColor=None):
elif par == 'imag':
zPlot = np.imag(array[flag][indUniSta])
elif par == 'res':
zPlot = (old_div(appResFact,freq))*np.abs(array[flag][indUniSta])**2
zPlot = (appResFact/freq)*np.abs(array[flag][indUniSta])**2
elif par == 'phs':
zPlot = np.arctan2(array[flag][indUniSta].imag,array[flag][indUniSta].real)*(old_div(180,np.pi))
zPlot = np.arctan2(array[flag][indUniSta].imag,array[flag][indUniSta].real)*(180/np.pi)
if not pColor:
if 'xx' in flag:
@@ -204,7 +198,7 @@ def plotIsoStaImpedance(ax,loc,array,flag,par='abs',pSym='s',pColor=None):
pColor = 'g'
elif 'xy' in flag:
lab = 'XY'
pColor = 'r'
pColor = 'r'
elif 'yx' in flag:
lab = 'YX'
pColor = 'b'
@@ -213,14 +207,14 @@ def plotIsoStaImpedance(ax,loc,array,flag,par='abs',pSym='s',pColor=None):
pColor = 'y'
ax.plot(freq,zPlot,color=pColor,marker=pSym,label=flag)
def plotPsudoSectNSimpedance(ax,sectDict,array,flag,par='abs',colorbar=True,colorNorm='None',cLevel=None,contour=True):
indSect = np.where(list(sectDict.values())[0]==array[list(sectDict.keys())[0]])
indSect = np.where(sectDict.values()[0]==array[sectDict.keys()[0]])
# Define the plot axes
if 'x' in list(sectDict.keys())[0]:
if 'x' in sectDict.keys()[0]:
x = array['y'][indSect]
else:
x = array['x'][indSect]
@@ -235,9 +229,9 @@ def plotPsudoSectNSimpedance(ax,sectDict,array,flag,par='abs',colorbar=True,colo
else:
level = np.linspace(zPlot.min(),zPlot.max(),100,endpoint=True)
clevel = np.linspace(zPlot.min(),zPlot.max(),10,endpoint=True)
elif par == 'ares':
zPlot = old_div(np.abs(array[flag][indSect])**2,(8*np.pi**2*10**(-7)*array['freq'][indSect]))
elif par == 'ares':
zPlot = np.abs(array[flag][indSect])**2/(8*np.pi**2*10**(-7)*array['freq'][indSect])
cmap = plt.get_cmap('RdYlBu')#seismic)
if cLevel:
zMax = np.log10(cLevel[1])
@@ -250,7 +244,7 @@ def plotPsudoSectNSimpedance(ax,sectDict,array,flag,par='abs',colorbar=True,colo
plotNorm = colors.LogNorm()
elif par == 'aphs':
zPlot = np.arctan2(array[flag][indSect].imag,array[flag][indSect].real)*(old_div(180,np.pi))
zPlot = np.arctan2(array[flag][indSect].imag,array[flag][indSect].real)*(180/np.pi)
cmap = plt.get_cmap('RdYlBu')#seismic)
if cLevel:
zMax = cLevel[1]
@@ -277,7 +271,7 @@ def plotPsudoSectNSimpedance(ax,sectDict,array,flag,par='abs',colorbar=True,colo
elif par == 'imag':
zPlot = np.imag(array[flag][indSect])
cmap = plt.get_cmap('Spectral') #('RdYlBu')
if cLevel:
zMax = np.log10(cLevel[1])
zMin = np.log10(cLevel[0])
@@ -287,7 +281,7 @@ def plotPsudoSectNSimpedance(ax,sectDict,array,flag,par='abs',colorbar=True,colo
level = np.concatenate((-np.logspace(zMax,zMin-.125,(zMax-zMin)*8+1,endpoint=True),np.logspace(zMin-.125,zMax,(zMax-zMin)*8+1,endpoint=True)))
clevel = np.concatenate((-np.logspace(zMax,zMin,(zMax-zMin)*1+1,endpoint=True),np.logspace(zMin,zMax,(zMax-zMin)*1+1,endpoint=True)))
plotNorm = colors.SymLogNorm(np.abs(level).min(),linscale=0.1)
if colorNorm=='SymLog':
plotNorm = colors.SymLogNorm(np.abs(level).min(),linscale=0.1)
elif colorNorm=='Lin':
@@ -313,14 +307,14 @@ def plotPsudoSectNSDiff(ax,sectDict,arrayList,flag,par='abs',colorbar=True,color
def sortInArr(arr):
return np.sort(arr,order=['freq','x','y','z'])
# Find the index for the slice
indSect0 = np.where(list(sectDict.values())[0]==arrayList[0][list(sectDict.keys())[0]])
indSect1 = np.where(list(sectDict.values())[0]==arrayList[1][list(sectDict.keys())[0]])
indSect0 = np.where(sectDict.values()[0]==arrayList[0][sectDict.keys()[0]])
indSect1 = np.where(sectDict.values()[0]==arrayList[1][sectDict.keys()[0]])
# Extract and sort the mats
arr0 = sortInArr(arrayList[0][indSect0])
arr1 = sortInArr(arrayList[1][indSect1])
# Define the plot axes
if 'x' in list(sectDict.keys())[0]:
if 'x' in sectDict.keys()[0]:
x0 = arr0['y']
x1 = arr1['y']
else:
@@ -328,35 +322,35 @@ def plotPsudoSectNSDiff(ax,sectDict,arrayList,flag,par='abs',colorbar=True,color
x1 = arr1['x']
y0 = arr0['freq']
y1 = arr1['freq']
if par == 'abs':
if useLog:
zPlot = old_div((np.log10(np.abs(arr0[flag])) - np.log10(np.abs(arr1[flag]))),np.log10(np.abs(arr1[flag])))
zPlot = (np.log10(np.abs(arr0[flag])) - np.log10(np.abs(arr1[flag])))/np.log10(np.abs(arr1[flag]))
else:
zPlot = old_div((np.abs(arr0[flag]) - np.abs(arr1[flag])),np.abs(arr1[flag]))
zPlot = (np.abs(arr0[flag]) - np.abs(arr1[flag]))/np.abs(arr1[flag])
if mask:
maskInd = np.logical_or(np.abs(arr0[flag])< 1e-3,np.abs(arr1[flag]) < 1e-3)
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
cmap = plt.get_cmap('RdYlBu')#seismic)
elif par == 'ares':
arF = old_div(1,(8*np.pi**2*10**(-7)))
arF = 1/(8*np.pi**2*10**(-7))
if useLog:
zPlot = old_div((np.log10((old_div(arF,arr0['freq']))*np.abs(arr0[flag])**2) - np.log10((old_div(arF,arr1['freq']))*np.abs(arr1[flag])**2)),np.log10((old_div(arF,arr1['freq']))*np.abs(arr1[flag])**2))
zPlot = (np.log10((arF/arr0['freq'])*np.abs(arr0[flag])**2) - np.log10((arF/arr1['freq'])*np.abs(arr1[flag])**2))/np.log10((arF/arr1['freq'])*np.abs(arr1[flag])**2)
else:
zPlot = old_div(((old_div(arF,arr0['freq']))*np.abs(arr0[flag])**2 - (old_div(arF,arr1['freq']))*np.abs(arr1[flag])**2),((old_div(arF,arr1['freq']))*np.abs(arr1[flag])**2))
zPlot = ((arF/arr0['freq'])*np.abs(arr0[flag])**2 - (arF/arr1['freq'])*np.abs(arr1[flag])**2)/((arF/arr1['freq'])*np.abs(arr1[flag])**2)
if mask:
maskInd = np.logical_or(np.abs(arr0[flag])< 1e-3,np.abs(arr1[flag]) < 1e-3)
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
cmap = plt.get_cmap('Spectral')#seismic)
elif par == 'aphs':
if useLog:
zPlot = old_div((np.log10(np.arctan2(arr0[flag].imag,arr0[flag].real)*(old_div(180,np.pi))) - np.log10(np.arctan2(arr1[flag].imag,arr1[flag].real)*(old_div(180,np.pi))) ),np.log10(np.arctan2(arr1[flag].imag,arr1[flag].real)*(old_div(180,np.pi))))
zPlot = (np.log10(np.arctan2(arr0[flag].imag,arr0[flag].real)*(180/np.pi)) - np.log10(np.arctan2(arr1[flag].imag,arr1[flag].real)*(180/np.pi)) )/np.log10(np.arctan2(arr1[flag].imag,arr1[flag].real)*(180/np.pi))
else:
zPlot = old_div(( np.arctan2(arr0[flag].imag,arr0[flag].real)*(old_div(180,np.pi)) - np.arctan2(arr1[flag].imag,arr1[flag].real)*(old_div(180,np.pi)) ),(np.arctan2(arr1[flag].imag,arr1[flag].real)*(old_div(180,np.pi))))
zPlot = ( np.arctan2(arr0[flag].imag,arr0[flag].real)*(180/np.pi) - np.arctan2(arr1[flag].imag,arr1[flag].real)*(180/np.pi) )/(np.arctan2(arr1[flag].imag,arr1[flag].real)*(180/np.pi))
if mask:
maskInd = np.logical_or(np.abs(arr0[flag])< 1e-3,np.abs(arr1[flag]) < 1e-3)
zPlot = np.ma.array(zPlot)
@@ -364,26 +358,26 @@ def plotPsudoSectNSDiff(ax,sectDict,arrayList,flag,par='abs',colorbar=True,color
cmap = plt.get_cmap('Spectral')#seismic)
elif par == 'real':
if useLog:
zPlot = old_div((np.log10(arr0[flag].real) - np.log10(arr1[flag].real)),np.log10(arr1[flag].real))
zPlot = (np.log10(arr0[flag].real) - np.log10(arr1[flag].real))/np.log10(arr1[flag].real)
else:
zPlot = old_div((arr0[flag].real - arr1[flag].real),arr1[flag].real)
zPlot = (arr0[flag].real - arr1[flag].real)/arr1[flag].real
if mask:
maskInd = np.logical_or(arr0[flag].real< 1e-3,arr1[flag].real < 1e-3)
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
cmap = plt.get_cmap('Spectral') #('Spectral')
elif par == 'imag':
if useLog:
zPlot = old_div((np.log10(arr0[flag].imag) - np.log10(arr1[flag].imag)),np.log10(arr1[flag].imag))
zPlot = (np.log10(arr0[flag].imag) - np.log10(arr1[flag].imag))/np.log10(arr1[flag].imag)
else:
zPlot = old_div((arr0[flag].imag - arr1[flag].imag),arr1[flag].imag)
zPlot = (arr0[flag].imag - arr1[flag].imag)/arr1[flag].imag
if mask:
maskInd = np.logical_or(arr0[flag].imag< 1e-3,arr1[flag].imag < 1e-3)
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
cmap = plt.get_cmap('Spectral') #('RdYlBu')
if cLevel:
zMax = np.log10(cLevel[1])
zMin = np.log10(cLevel[0])
@@ -391,18 +385,18 @@ def plotPsudoSectNSDiff(ax,sectDict,arrayList,flag,par='abs',colorbar=True,color
zMax = (np.ceil(np.log10(np.abs(zPlot).max())))
zMin = (np.floor(np.log10(np.abs(zPlot).min())))
if colorNorm=='SymLog':
level = np.concatenate((-np.logspace(zMax,zMin-.125,(zMax-zMin)*8+1,endpoint=True),np.logspace(zMin-.125,zMax,(zMax-zMin)*8+1,endpoint=True)))
clevel = np.concatenate((-np.logspace(zMax,zMin,(zMax-zMin)*1+1,endpoint=True),np.logspace(zMin,zMax,(zMax-zMin)*1+1,endpoint=True)))
plotNorm = colors.SymLogNorm(np.abs(level).min(),linscale=0.1)
elif colorNorm=='Lin':
if cLevel:
level = np.arange(cLevel[0],cLevel[1]+.1,old_div((cLevel[1] - cLevel[0]),50.))
clevel = np.arange(cLevel[0],cLevel[1]+.1,old_div((cLevel[1] - cLevel[0]),10.))
level = np.arange(cLevel[0],cLevel[1]+.1,(cLevel[1] - cLevel[0])/50.)
clevel = np.arange(cLevel[0],cLevel[1]+.1,(cLevel[1] - cLevel[0])/10.)
else:
level = np.arange(zPlot.min(),zPlot.max(),old_div((zPlot.max() - zPlot.min()),50.))
clevel = np.arange(zPlot.min(),zPlot.max(),old_div((zPlot.max() - zPlot.min()),10.))
level = np.arange(zPlot.min(),zPlot.max(),(zPlot.max() - zPlot.min())/50.)
clevel = np.arange(zPlot.min(),zPlot.max(),(zPlot.max() - zPlot.min())/10.)
plotNorm = colors.Normalize()
elif colorNorm=='Log':
level = np.logspace(zMin-.125,zMax,(zMax-zMin)*8+1,endpoint=True)
+1 -9
View File
@@ -1,11 +1,3 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from builtins import dict
from future import standard_library
standard_library.install_aliases()
from builtins import zip
import SimPEG as simpeg, numpy as np
def homo1DModelSource(mesh,freq,sigma_1d):
@@ -98,7 +90,7 @@ def analytic1DModelSource(mesh,freq,sigma_1d):
Eu, Ed, _, _ = getEHfields(mesh1d,sigma_1d,freq,mesh.vectorNz)
# Make the fields into a dictionary of location and the fields
e0_1d = Eu+Ed
E1dFieldDict = dict(list(zip(mesh.vectorNz,e0_1d)))
E1dFieldDict = dict(zip(mesh.vectorNz,e0_1d))
if mesh.dim == 1:
eBG_px = simpeg.mkvc(e0_1d,2)
eBG_py = -simpeg.mkvc(e0_1d,2) # added a minus to make the results in the correct quadrents.
-6
View File
@@ -1,9 +1,3 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
import SimPEG as simpeg, numpy as np
def homo1DModelSource(mesh,freq,m_back):

Some files were not shown because too many files have changed in this diff Show More