Compare commits

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

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