mirror of
https://github.com/wassname/simpeg.git
synced 2026-09-13 13:03:14 +08:00
Compare commits
6
Commits
dev
...
feat/Python3
| Author | SHA1 | Date | |
|---|---|---|---|
|
|
714afa860f | ||
|
|
f7a70aa6a7 | ||
|
|
7189ec5b2f | ||
|
|
4c2a61bf54 | ||
|
|
ca8d8f8c2d | ||
|
|
362975d2bd |
+2
-1
@@ -1,6 +1,7 @@
|
||||
language: python
|
||||
python:
|
||||
- 2.7
|
||||
- 3.4
|
||||
|
||||
sudo: false
|
||||
|
||||
@@ -38,7 +39,7 @@ before_install:
|
||||
-O miniconda.sh; fi
|
||||
- chmod +x miniconda.sh
|
||||
- ./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
|
||||
|
||||
install:
|
||||
|
||||
@@ -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 *
|
||||
|
||||
class FieldsDC_CC(Problem.Fields):
|
||||
@@ -61,7 +69,7 @@ class SrcDipole(Survey.BaseSrc):
|
||||
pts = [self.loc[0], self.loc[1]]
|
||||
inds = Utils.closestPoints(prob.mesh, pts)
|
||||
q = np.zeros(prob.mesh.nC)
|
||||
q[inds] = - self.current * ( np.r_[1., -1.] / prob.mesh.vol[inds] )
|
||||
q[inds] = - self.current * (np.r_[1., -1.] / prob.mesh.vol[inds])
|
||||
# self._rhsDict[mesh] = q
|
||||
# return self._rhsDict[mesh]
|
||||
return q
|
||||
|
||||
@@ -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 BaseDC import SurveyDC, FieldsDC_CC
|
||||
from .BaseDC import SurveyDC, FieldsDC_CC
|
||||
|
||||
class SurveyIP(SurveyDC):
|
||||
"""
|
||||
@@ -52,7 +59,7 @@ class ProblemIP(Problem.BaseProblem):
|
||||
# sigma = self.curModel.transform
|
||||
sigma = self.sigma
|
||||
Av = self.mesh.aveF2CC
|
||||
self._Msig = Utils.sdiag(1/(self.mesh.dim * Av.T * (1/sigma)))
|
||||
self._Msig = Utils.sdiag(1//(self.mesh.dim * Av.T * (1/sigma)))
|
||||
return self._Msig
|
||||
|
||||
@property
|
||||
|
||||
+31
-21
@@ -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
|
||||
import BaseDC as DC
|
||||
import BaseDC as IP
|
||||
from . import BaseDC as DC
|
||||
from . import BaseDC as IP
|
||||
import warnings
|
||||
|
||||
def getActiveindfromTopo(mesh, topo):
|
||||
@@ -67,7 +77,7 @@ def readUBC_DC3Dobstopo(filename,mesh,topo,probType="CC"):
|
||||
if "!" in line.split(): 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
|
||||
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
|
||||
@@ -228,10 +238,10 @@ def plot_pseudoSection(DCsurvey, axs, surveyType='dipole-dipole', unitType='volt
|
||||
|
||||
elif surveyType == 'dipole-dipole':
|
||||
|
||||
leg = data * 2*np.pi / ( 1/MA - 1/MB - 1/NB + 1/NA )
|
||||
leg = data * 2*np.pi / (1/MA - 1/MB - 1/NB + 1/NA)
|
||||
|
||||
else:
|
||||
print """unitType must be 'pole-dipole' | 'dipole-dipole' """
|
||||
print("""unitType must be 'pole-dipole' | 'dipole-dipole' """)
|
||||
break
|
||||
|
||||
|
||||
@@ -246,11 +256,11 @@ def plot_pseudoSection(DCsurvey, axs, surveyType='dipole-dipole', unitType='volt
|
||||
rho = np.hstack([rho,leg])
|
||||
|
||||
else:
|
||||
print """unitType must be 'appResistivity' | 'appConductivity' | 'volt' """
|
||||
print("""unitType must be 'appResistivity' | 'appConductivity' | 'volt' """)
|
||||
break
|
||||
|
||||
midx = np.hstack([midx, ( Cmid + Pmid )/2 ])
|
||||
midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + (Tx[0][2] + Tx[1][2])/2 ])
|
||||
midx = np.hstack([midx, (Cmid + Pmid)/2])
|
||||
midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + (Tx[0][2] + Tx[1][2])/2])
|
||||
|
||||
# Grid points
|
||||
grid_x, grid_z = np.mgrid[np.min(midx):np.max(midx), np.min(midz):np.max(midz)]
|
||||
@@ -340,11 +350,11 @@ def gen_DCIPsurvey(endl, mesh, surveyType, AM_sep, MN_sep, nrx):
|
||||
dl_x = ( endl[1,0] - endl[0,0] ) / dl_len
|
||||
dl_y = ( endl[1,1] - endl[0,1] ) / dl_len
|
||||
|
||||
nstn = np.floor( dl_len / AM_sep )
|
||||
nstn = np.floor(dl_len / AM_sep)
|
||||
|
||||
# Compute discrete pole location along line
|
||||
stn_x = endl[0,0] + np.array(range(int(nstn)))*dl_x*AM_sep
|
||||
stn_y = endl[0,1] + np.array(range(int(nstn)))*dl_y*AM_sep
|
||||
stn_x = endl[0,0] + np.array(list(range(int(nstn))))*dl_x*AM_sep
|
||||
stn_y = endl[0,1] + np.array(list(range(int(nstn))))*dl_y*AM_sep
|
||||
|
||||
# Create line of P1 locations
|
||||
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])
|
||||
|
||||
# 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
|
||||
if nstn <= 0:
|
||||
continue
|
||||
|
||||
# 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_y = N[ii,1] + dl_y*MN_sep + np.array(range(int(nstn)))*dl_y*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(list(range(int(nstn))))*dl_y*AM_sep
|
||||
|
||||
# Create receiver poles
|
||||
# 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_w = box_l/2.
|
||||
|
||||
nstn = np.floor( box_l / AM_sep )
|
||||
nstn = box_l // AM_sep
|
||||
|
||||
# Compute discrete pole location along line
|
||||
stn_x = min_x + np.array(range(int(nstn)))*dl_x*AM_sep
|
||||
stn_y = min_y + np.array(range(int(nstn)))*dl_y*AM_sep
|
||||
stn_x = min_x + np.array(list(range(int(nstn))))*dl_x*AM_sep
|
||||
stn_y = min_y + np.array(list(range(int(nstn))))*dl_y*AM_sep
|
||||
|
||||
# Define number of cross lines
|
||||
nlin = int(np.floor( box_w / AM_sep ))
|
||||
lind = range(-nlin,nlin+1)
|
||||
nlin = int(box_w // AM_sep)
|
||||
lind = list(range(-nlin,nlin+1))
|
||||
|
||||
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,:])
|
||||
SrcList.append(srcClass)
|
||||
else:
|
||||
print """surveyType must be either 'pole-dipole', 'dipole-dipole' or 'gradient'. """
|
||||
print("""surveyType must be either 'pole-dipole', 'dipole-dipole' or 'gradient'. """)
|
||||
|
||||
survey = DC.SurveyDC(SrcList)
|
||||
return survey, Tx, Rx
|
||||
@@ -668,7 +678,7 @@ def readUBC_DC3Dobs(fileName, rtype = 'DC'):
|
||||
obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!')
|
||||
|
||||
else:
|
||||
print "rtype must be 'DC'(default) | 'IP'"
|
||||
print("rtype must be 'DC'(default) | 'IP'")
|
||||
|
||||
# Pre-allocate
|
||||
srcLists = []
|
||||
|
||||
@@ -1,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
|
||||
|
||||
def WennerSrcList(nElecs, aSpacing, in2D=False, plotIt=False):
|
||||
@@ -5,7 +12,7 @@ def WennerSrcList(nElecs, aSpacing, in2D=False, plotIt=False):
|
||||
import SimPEG.DCIP as DC
|
||||
|
||||
elocs = np.arange(0,aSpacing*nElecs,aSpacing)
|
||||
elocs -= (nElecs*aSpacing - aSpacing)/2
|
||||
elocs -= (nElecs*aSpacing - aSpacing) / 2
|
||||
space = 1
|
||||
WENNER = np.zeros((0,),dtype=int)
|
||||
for ii in range(nElecs):
|
||||
|
||||
+10
-4
@@ -1,4 +1,10 @@
|
||||
from BaseDC import *
|
||||
from BaseIP import *
|
||||
from DCIPUtils import *
|
||||
import Utils
|
||||
from __future__ import absolute_import
|
||||
from __future__ import unicode_literals
|
||||
from __future__ import print_function
|
||||
from __future__ import division
|
||||
from future import standard_library
|
||||
standard_library.install_aliases()
|
||||
from .BaseDC import *
|
||||
from .BaseIP import *
|
||||
from .DCIPUtils import *
|
||||
from . import Utils
|
||||
|
||||
+13
-6
@@ -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
|
||||
|
||||
.. note::
|
||||
@@ -9,8 +18,6 @@ class BaseDataMisfit(object):
|
||||
You should inherit from this class to create your own data misfit term.
|
||||
"""
|
||||
|
||||
__metaclass__ = Utils.SimPEGMetaClass
|
||||
|
||||
debug = False #: Print debugging information
|
||||
counter = None #: Set this to a SimPEG.Utils.Counter() if you want to count things
|
||||
|
||||
@@ -93,11 +100,11 @@ class l2_DataMisfit(BaseDataMisfit):
|
||||
survey = self.survey
|
||||
|
||||
if getattr(survey,'std', None) is None:
|
||||
print 'SimPEG.DataMisfit.l2_DataMisfit assigning default std of 5%'
|
||||
print('SimPEG.DataMisfit.l2_DataMisfit assigning default std of 5%')
|
||||
survey.std = 0.05
|
||||
|
||||
if getattr(survey, 'eps', None) is None:
|
||||
print 'SimPEG.DataMisfit.l2_DataMisfit assigning default eps of 1e-5 * ||dobs||'
|
||||
print('SimPEG.DataMisfit.l2_DataMisfit assigning default eps of 1e-5 * ||dobs||')
|
||||
survey.eps = np.linalg.norm(Utils.mkvc(survey.dobs),2)*1e-5
|
||||
|
||||
self._Wd = Utils.sdiag(1/(abs(survey.dobs)*survey.std+survey.eps))
|
||||
|
||||
+31
-20
@@ -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):
|
||||
"""InversionDirective"""
|
||||
@@ -15,7 +26,7 @@ class InversionDirective(object):
|
||||
@inversion.setter
|
||||
def inversion(self, i):
|
||||
if getattr(self,'_inversion',None) is not None:
|
||||
print 'Warning: InversionDirective %s has switched to a new inversion.' % self.__name__
|
||||
print('Warning: InversionDirective %s has switched to a new inversion.' % self.__name__)
|
||||
self._inversion = i
|
||||
|
||||
@property
|
||||
@@ -68,7 +79,7 @@ class DirectiveList(object):
|
||||
def inversion(self, i):
|
||||
if self.inversion is i: return
|
||||
if getattr(self,'_inversion',None) is not None:
|
||||
print 'Warning: %s has switched to a new inversion.' % self.__name__
|
||||
print('Warning: %s has switched to a new inversion.' % self.__name__)
|
||||
for d in self.dList:
|
||||
d.inversion = i
|
||||
self._inversion = i
|
||||
@@ -120,7 +131,7 @@ class BetaEstimate_ByEig(InversionDirective):
|
||||
:return: beta0
|
||||
"""
|
||||
|
||||
if self.debug: print 'Calculating the beta0 parameter.'
|
||||
if self.debug: print('Calculating the beta0 parameter.')
|
||||
|
||||
m = self.invProb.curModel
|
||||
f = self.invProb.getFields(m, store=True, deleteWarmstart=False)
|
||||
@@ -141,7 +152,7 @@ class BetaSchedule(InversionDirective):
|
||||
|
||||
def endIter(self):
|
||||
if self.opt.iter > 0 and self.opt.iter % self.coolingRate == 0:
|
||||
if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter
|
||||
if self.debug: print('BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter)
|
||||
self.invProb.beta /= self.coolingFactor
|
||||
|
||||
|
||||
@@ -192,7 +203,7 @@ class SaveModelEveryIteration(SaveEveryIteration):
|
||||
"""SaveModelEveryIteration"""
|
||||
|
||||
def initialize(self):
|
||||
print "SimPEG.SaveModelEveryIteration will save your models as: '###-%s.npy'"%self.fileName
|
||||
print("SimPEG.SaveModelEveryIteration will save your models as: '###-%s.npy'"%self.fileName)
|
||||
|
||||
def endIter(self):
|
||||
np.save('%03d-%s' % (self.opt.iter, self.fileName), self.opt.xc)
|
||||
@@ -202,7 +213,7 @@ class SaveOutputEveryIteration(SaveEveryIteration):
|
||||
"""SaveModelEveryIteration"""
|
||||
|
||||
def initialize(self):
|
||||
print "SimPEG.SaveOutputEveryIteration will save your inversion progress as: '###-%s.txt'"%self.fileName
|
||||
print("SimPEG.SaveOutputEveryIteration will save your inversion progress as: '###-%s.txt'"%self.fileName)
|
||||
f = open(self.fileName+'.txt', 'w')
|
||||
f.write(" # beta phi_d phi_m f\n")
|
||||
f.close()
|
||||
@@ -216,7 +227,7 @@ class SaveOutputDictEveryIteration(SaveEveryIteration):
|
||||
"""SaveOutputDictEveryIteration"""
|
||||
|
||||
def initialize(self):
|
||||
print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '###-%s.npz'"%self.fileName
|
||||
print("SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '###-%s.npz'"%self.fileName)
|
||||
|
||||
def endIter(self):
|
||||
# Save the data.
|
||||
@@ -294,7 +305,7 @@ class Update_IRLS(InversionDirective):
|
||||
|
||||
# After reaching target misfit with l2-norm, switch to IRLS (mode:2)
|
||||
if self.invProb.phi_d < self.target and self.mode == 1:
|
||||
print "Convergence with smooth l2-norm regularization: Start IRLS steps..."
|
||||
print("Convergence with smooth l2-norm regularization: Start IRLS steps...")
|
||||
|
||||
self.mode = 2
|
||||
|
||||
@@ -302,17 +313,17 @@ class Update_IRLS(InversionDirective):
|
||||
# model values
|
||||
if getattr(self, 'reg.eps', None) is None:
|
||||
self.reg.eps_p = np.percentile(np.abs(self.invProb.curModel),self.prctile)
|
||||
else:
|
||||
else:
|
||||
self.reg.eps_p = self.eps[0]
|
||||
|
||||
|
||||
if getattr(self, 'reg.eps', None) is None:
|
||||
self.reg.eps_q = np.percentile(np.abs(self.reg.regmesh.cellDiffxStencil*(self.reg.mapping * self.invProb.curModel)),self.prctile)
|
||||
else:
|
||||
else:
|
||||
self.reg.eps_q = self.eps[1]
|
||||
|
||||
print "L[p qx qy qz]-norm : " + str(self.reg.norms)
|
||||
print "eps_p: " + str(self.reg.eps_p) + " eps_q: " + str(self.reg.eps_q)
|
||||
|
||||
|
||||
print("L[p qx qy qz]-norm : " + str(self.reg.norms))
|
||||
print("eps_p: " + str(self.reg.eps_p) + " eps_q: " + str(self.reg.eps_q))
|
||||
|
||||
self.reg.norms = self.norms
|
||||
self.coolingFactor = 1.
|
||||
self.coolingRate = 1
|
||||
@@ -328,7 +339,7 @@ class Update_IRLS(InversionDirective):
|
||||
|
||||
# Beta Schedule
|
||||
if self.opt.iter > 0 and self.opt.iter % self.coolingRate == 0:
|
||||
if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter
|
||||
if self.debug: print('BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter)
|
||||
self.invProb.beta /= self.coolingFactor
|
||||
|
||||
|
||||
@@ -340,17 +351,17 @@ class Update_IRLS(InversionDirective):
|
||||
phim_new = self.reg.eval(self.invProb.curModel)
|
||||
self.f_change = np.abs(self.f_old - phim_new) / self.f_old
|
||||
|
||||
print "Regularization decrease: %6.3e" % (self.f_change)
|
||||
print("Regularization decrease: %6.3e" % (self.f_change))
|
||||
|
||||
# Check for maximum number of IRLS cycles
|
||||
if self.IRLSiter == self.maxIRLSiter:
|
||||
print "Reach maximum number of IRLS cycles: %i" % self.maxIRLSiter
|
||||
print("Reach maximum number of IRLS cycles: %i" % self.maxIRLSiter)
|
||||
self.opt.stopNextIteration = True
|
||||
return
|
||||
|
||||
# Check if the function has changed enough
|
||||
if self.f_change < self.f_min_change and self.IRLSiter > 1:
|
||||
print "Minimum decrease in regularization. End of IRLS"
|
||||
print("Minimum decrease in regularization. End of IRLS")
|
||||
self.opt.stopNextIteration = True
|
||||
return
|
||||
else:
|
||||
|
||||
@@ -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
|
||||
from scipy.constants import mu_0, pi
|
||||
from scipy import special
|
||||
|
||||
@@ -1,4 +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
|
||||
from scipy.constants import mu_0, pi
|
||||
from scipy.special import erf
|
||||
|
||||
@@ -1,4 +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
|
||||
from scipy.constants import mu_0, pi, epsilon_0
|
||||
from scipy.special import erf
|
||||
|
||||
@@ -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 scipy.constants import mu_0, epsilon_0
|
||||
from SimPEG.EM.Utils.EMUtils import k
|
||||
@@ -34,7 +40,7 @@ def _getCasingHertzMagDipoleDeriv_r(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones
|
||||
sqrtr2z2 = np.sqrt(r2 + dxyz[:,2]**2)
|
||||
k2 = k(freq,sigma[2],mu[2],eps)
|
||||
|
||||
return -HertzZ * np.sqrt(r2) / sqrtr2z2 * (1j*k2 + 1./ sqrtr2z2)
|
||||
return -HertzZ * np.sqrt(r2) / sqrtr2z2 * (1j*k2 + 1. / sqrtr2z2)
|
||||
|
||||
|
||||
def _getCasingHertzMagDipoleDeriv_z(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.):
|
||||
|
||||
@@ -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
|
||||
from scipy.constants import mu_0, pi
|
||||
from scipy.special import erf
|
||||
|
||||
@@ -1,5 +1,11 @@
|
||||
from TDEM import hzAnalyticDipoleT
|
||||
from FDEM import hzAnalyticDipoleF
|
||||
from FDEMcasing import *
|
||||
from DC import DCAnalyticHalf, DCAnalyticSphere
|
||||
from FDEMDipolarfields 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 .TDEM import hzAnalyticDipoleT
|
||||
from .FDEM import hzAnalyticDipoleF
|
||||
from .FDEMcasing import *
|
||||
from .DC import DCAnalyticHalf, DCAnalyticSphere
|
||||
from .FDEMDipolarfields 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()
|
||||
from SimPEG import Survey, Problem, Utils, Models, Maps, PropMaps, np, sp, Solver as SimpegSolver
|
||||
from scipy.constants import mu_0
|
||||
|
||||
|
||||
@@ -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 scipy.sparse as sp
|
||||
import SimPEG
|
||||
@@ -42,7 +49,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
|
||||
:return: total electric field
|
||||
"""
|
||||
if getattr(self, '_ePrimary', None) is None or getattr(self, '_eSecondary', None) is None:
|
||||
raise NotImplementedError ('Getting e from %s is not implemented' %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)
|
||||
|
||||
@@ -56,7 +63,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
|
||||
:return: total magnetic flux density
|
||||
"""
|
||||
if getattr(self, '_bPrimary', None) is None or getattr(self, '_bSecondary', None) is None:
|
||||
raise NotImplementedError ('Getting b from %s is not implemented' %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)
|
||||
|
||||
@@ -70,7 +77,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
|
||||
:return: total magnetic field
|
||||
"""
|
||||
if getattr(self, '_hPrimary', None) is None or getattr(self, '_hSecondary', None) is None:
|
||||
raise NotImplementedError ('Getting h from %s is not implemented' %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)
|
||||
|
||||
@@ -84,7 +91,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
|
||||
:return: total current density
|
||||
"""
|
||||
if getattr(self, '_jPrimary', None) is None or getattr(self, '_jSecondary', None) is None:
|
||||
raise NotImplementedError ('Getting j from %s is not implemented' %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)
|
||||
|
||||
@@ -100,7 +107,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
|
||||
:return: derivative times a vector (or tuple for adjoint)
|
||||
"""
|
||||
if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None:
|
||||
raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0])
|
||||
raise NotImplementedError ('Getting eDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
|
||||
|
||||
if 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)
|
||||
"""
|
||||
if getattr(self, '_bDeriv_u', None) is None or getattr(self, '_bDeriv_m', None) is None:
|
||||
raise NotImplementedError ('Getting bDerivs from %s is not implemented' %self.knownFields.keys()[0])
|
||||
raise NotImplementedError ('Getting bDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
|
||||
|
||||
if 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)
|
||||
"""
|
||||
if getattr(self, '_hDeriv_u', None) is None or getattr(self, '_hDeriv_m', None) is None:
|
||||
raise NotImplementedError ('Getting hDerivs from %s is not implemented' %self.knownFields.keys()[0])
|
||||
raise NotImplementedError ('Getting hDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
|
||||
|
||||
if 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)
|
||||
"""
|
||||
if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None:
|
||||
raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0])
|
||||
raise NotImplementedError ('Getting jDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
|
||||
|
||||
if adjoint:
|
||||
return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint)
|
||||
@@ -285,7 +292,7 @@ class Fields3D_e(FieldsFDEM):
|
||||
C = self._edgeCurl
|
||||
b = (C * eSolution)
|
||||
for i, src in enumerate(srcList):
|
||||
b[:,i] *= - 1./(1j*omega(src.freq))
|
||||
b[:,i] *= -1./(1j*omega(src.freq))
|
||||
s_m, _ = src.eval(self.prob)
|
||||
b[:,i] = b[:,i]+ 1./(1j*omega(src.freq)) * s_m
|
||||
return b
|
||||
@@ -345,7 +352,7 @@ class Fields3D_e(FieldsFDEM):
|
||||
:rtype: numpy.ndarray
|
||||
:return: product of the derivative of the current density with respect to the field we solved for with a vector
|
||||
"""
|
||||
n = int(self._aveE2CCV.shape[0] / self._nC) # number of components (instead of checking if cyl or not)
|
||||
n = int(self._aveE2CCV.shape[0] // self._nC) # number of components (instead of checking if cyl or not)
|
||||
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
|
||||
|
||||
if adjoint:
|
||||
@@ -382,8 +389,8 @@ class Fields3D_e(FieldsFDEM):
|
||||
:rtype: numpy.ndarray
|
||||
:return: magnetic field
|
||||
"""
|
||||
n = int(self._aveF2CCV.shape[0] / self._nC) # Number of Components
|
||||
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
|
||||
n = int(self._aveF2CCV.shape[0] // self._nC) # Number of Components
|
||||
VI = sdiag(np.kron(np.ones(n), 1. // self.prob.mesh.vol))
|
||||
|
||||
return VI * (self._aveF2CCV * (self._MfMui * self._b(eSolution, srcList)))
|
||||
|
||||
@@ -397,7 +404,7 @@ class Fields3D_e(FieldsFDEM):
|
||||
:rtype: numpy.ndarray
|
||||
:return: product of the derivative of the magnetic field with respect to the field we solved for with a vector
|
||||
"""
|
||||
n = int(self._aveF2CCV.shape[0] / self._nC) # Number of Components
|
||||
n = int(self._aveF2CCV.shape[0] // self._nC) # Number of Components
|
||||
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
|
||||
if adjoint:
|
||||
v = self._MfMui.T * (self._aveF2CCV.T * (VI.T * du_dm_v))
|
||||
@@ -414,7 +421,7 @@ class Fields3D_e(FieldsFDEM):
|
||||
:rtype: numpy.ndarray
|
||||
:return: product of the magnetic field derivative with respect to the inversion model with a vector
|
||||
"""
|
||||
n = int(self._aveF2CCV.shape[0] / self._nC) # Number of Components
|
||||
n = int(self._aveF2CCV.shape[0] // self._nC) # Number of Components
|
||||
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
|
||||
if adjoint:
|
||||
v = self._MfMui.T * (self._aveF2CCV.T * (VI.T * v))
|
||||
@@ -607,7 +614,7 @@ class Fields3D_b(FieldsFDEM):
|
||||
:return: primary current density
|
||||
"""
|
||||
|
||||
n = int(self._aveE2CCV.shape[0] / self._nC) # number of components
|
||||
n = int(self._aveE2CCV.shape[0] // self._nC) # number of components
|
||||
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
|
||||
|
||||
return VI * (self._aveE2CCV * ( self._MeSigma * self._e(bSolution,srcList ) ) )
|
||||
@@ -624,7 +631,7 @@ class Fields3D_b(FieldsFDEM):
|
||||
:rtype: numpy.ndarray
|
||||
:return: product of the derivative of the current density with respect to the field we solved for with a vector
|
||||
"""
|
||||
n = int(self._aveE2CCV.shape[0] / self._nC) # number of components
|
||||
n = int(self._aveE2CCV.shape[0] // self._nC) # number of components
|
||||
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
|
||||
if adjoint:
|
||||
return self._MfMui.T * ( self._edgeCurl * ( self._aveE2CCV.T * (VI.T * du_dm_v) ) )
|
||||
@@ -652,7 +659,7 @@ class Fields3D_b(FieldsFDEM):
|
||||
:rtype: numpy.ndarray
|
||||
:return: magnetic field
|
||||
"""
|
||||
n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
|
||||
n = int(self._aveF2CCV.shape[0] // self._nC) #number of components
|
||||
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
|
||||
return VI * (self._aveF2CCV * (self._MfMui * self._b(bSolution, srcList)))
|
||||
|
||||
@@ -667,7 +674,7 @@ class Fields3D_b(FieldsFDEM):
|
||||
:rtype: numpy.ndarray
|
||||
:return: product of the derivative of the magnetic field with respect to the field we solved for with a vector
|
||||
"""
|
||||
n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
|
||||
n = int(self._aveF2CCV.shape[0] // self._nC) #number of components
|
||||
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
|
||||
|
||||
if adjoint:
|
||||
@@ -890,7 +897,7 @@ class Fields3D_j(FieldsFDEM):
|
||||
:rtype: numpy.ndarray
|
||||
:return: electric field
|
||||
"""
|
||||
n = int(self._aveF2CCV.shape[0] / self._nC) # number of components
|
||||
n = int(self._aveF2CCV.shape[0] // self._nC) # number of components
|
||||
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
|
||||
return VI * (self._aveF2CCV * (self._MfRho * self._j(jSolution, srcList)))
|
||||
|
||||
@@ -904,7 +911,7 @@ class Fields3D_j(FieldsFDEM):
|
||||
:rtype: numpy.ndarray
|
||||
:return: product of the derivative of the electric field with respect to the field we solved for with a vector
|
||||
"""
|
||||
n = int(self._aveF2CCV.shape[0] / self._nC) # number of components
|
||||
n = int(self._aveF2CCV.shape[0] // self._nC) # number of components
|
||||
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
|
||||
if adjoint:
|
||||
return self._MfRho.T * ( self._aveF2CCV.T * ( VI.T * du_dm_v ) )
|
||||
@@ -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
|
||||
"""
|
||||
jSolution = Utils.mkvc(self[src,'jSolution'])
|
||||
n = int(self._aveF2CCV.shape[0] / self._nC) # number of components
|
||||
n = int(self._aveF2CCV.shape[0] // self._nC) # number of components
|
||||
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
|
||||
if adjoint:
|
||||
return self._MfRhoDeriv(jSolution).T * ( self._aveF2CCV.T * ( VI.T * v ) )
|
||||
@@ -936,7 +943,7 @@ class Fields3D_j(FieldsFDEM):
|
||||
:rtype: numpy.ndarray
|
||||
:return: secondary magnetic flux density
|
||||
"""
|
||||
n = int(self._aveE2CCV.shape[0] / self._nC) # number of components
|
||||
n = int(self._aveE2CCV.shape[0] // self._nC) # number of components
|
||||
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
|
||||
|
||||
return VI * (self._aveE2CCV * ( self._MeMu * self._h(jSolution,srcList)) )
|
||||
@@ -951,7 +958,7 @@ class Fields3D_j(FieldsFDEM):
|
||||
:rtype: numpy.ndarray
|
||||
:return: product of the derivative of the magnetic flux density with respect to the field we solved for with a vector
|
||||
"""
|
||||
n = int(self._aveF2CCV.shape[0] / self._nC) # number of components
|
||||
n = int(self._aveF2CCV.shape[0] // self._nC) # number of components
|
||||
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
|
||||
|
||||
if adjoint:
|
||||
@@ -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
|
||||
"""
|
||||
jSolution = self[src,'jSolution']
|
||||
n = int(self._aveE2CCV.shape[0] / self._nC) # number of components
|
||||
n = int(self._aveE2CCV.shape[0] // self._nC) # number of components
|
||||
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
|
||||
s_mDeriv,_ = src.evalDeriv(self.prob, adjoint = adjoint)
|
||||
|
||||
@@ -1151,7 +1158,7 @@ class Fields3D_h(FieldsFDEM):
|
||||
:rtype: numpy.ndarray
|
||||
:return: electric field
|
||||
"""
|
||||
n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
|
||||
n = int(self._aveF2CCV.shape[0] // self._nC) #number of components
|
||||
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
|
||||
return VI * (self._aveF2CCV * (self._MfRho * self._j(hSolution, srcList)))
|
||||
|
||||
@@ -1165,7 +1172,7 @@ class Fields3D_h(FieldsFDEM):
|
||||
:rtype: numpy.ndarray
|
||||
:return: product of the derivative of the electric field with respect to the field we solved for with a vector
|
||||
"""
|
||||
n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
|
||||
n = int(self._aveF2CCV.shape[0] // self._nC) #number of components
|
||||
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
|
||||
if adjoint:
|
||||
return self._edgeCurl.T * ( self._MfRho.T * ( self._aveF2CCV.T * ( VI.T * du_dm_v ) ) )
|
||||
@@ -1182,7 +1189,7 @@ class Fields3D_h(FieldsFDEM):
|
||||
:return: product of the electric field derivative with respect to the inversion model with a vector
|
||||
"""
|
||||
hSolution = Utils.mkvc(self[src,'hSolution'])
|
||||
n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
|
||||
n = int(self._aveF2CCV.shape[0] // self._nC) #number of components
|
||||
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
|
||||
if adjoint:
|
||||
return ( self._MfRhoDeriv(self._edgeCurl * hSolution).T * ( self._aveF2CCV.T * (VI.T * v) ) )
|
||||
@@ -1198,7 +1205,7 @@ class Fields3D_h(FieldsFDEM):
|
||||
:return: magnetic flux density
|
||||
"""
|
||||
h = self._h(hSolution, srcList)
|
||||
n = int(self._aveE2CCV.shape[0] / self._nC) #number of components
|
||||
n = int(self._aveE2CCV.shape[0] // self._nC) #number of components
|
||||
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
|
||||
|
||||
return VI * (self._aveE2CCV * (self._MeMu * h))
|
||||
@@ -1213,7 +1220,7 @@ class Fields3D_h(FieldsFDEM):
|
||||
:rtype: numpy.ndarray
|
||||
:return: product of the derivative of the magnetic flux density with respect to the field we solved for with a vector
|
||||
"""
|
||||
n = int(self._aveE2CCV.shape[0] / self._nC) #number of components
|
||||
n = int(self._aveE2CCV.shape[0] // self._nC) #number of components
|
||||
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
|
||||
if adjoint:
|
||||
return self._MeMu.T * (self._aveE2CCV.T * ( VI.T * du_dm_v ))
|
||||
|
||||
@@ -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 scipy.constants import mu_0
|
||||
from SurveyFDEM import Survey as SurveyFDEM
|
||||
from FieldsFDEM import FieldsFDEM, Fields3D_e, Fields3D_b, Fields3D_h, Fields3D_j
|
||||
from .SurveyFDEM import Survey as SurveyFDEM
|
||||
from .FieldsFDEM import FieldsFDEM, Fields3D_e, Fields3D_b, Fields3D_h, Fields3D_j
|
||||
from SimPEG.EM.Base import BaseEMProblem
|
||||
from SimPEG.EM.Utils import omega
|
||||
|
||||
|
||||
@@ -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
|
||||
from SimPEG import sp
|
||||
|
||||
|
||||
@@ -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 scipy.constants import mu_0
|
||||
from SimPEG.EM.Utils import *
|
||||
|
||||
@@ -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
|
||||
from SimPEG.EM.Utils import *
|
||||
from SimPEG.EM.Base import BaseEMSurvey
|
||||
from scipy.constants import mu_0
|
||||
from SimPEG.Utils import Zero, Identity
|
||||
import SrcFDEM as Src
|
||||
import RxFDEM as Rx
|
||||
from . import SrcFDEM as Src
|
||||
from . import RxFDEM as Rx
|
||||
from SimPEG import sp
|
||||
|
||||
class Survey(BaseEMSurvey):
|
||||
|
||||
@@ -1,5 +1,11 @@
|
||||
from SurveyFDEM import Survey
|
||||
import SrcFDEM as Src
|
||||
import RxFDEM as Rx
|
||||
from ProblemFDEM import Problem3D_e, Problem3D_b, Problem3D_j, Problem3D_h
|
||||
from FieldsFDEM import Fields3D_e, Fields3D_b, Fields3D_j, Fields3D_h
|
||||
from __future__ import absolute_import
|
||||
from __future__ import unicode_literals
|
||||
from __future__ import print_function
|
||||
from __future__ import division
|
||||
from future import standard_library
|
||||
standard_library.install_aliases()
|
||||
from .SurveyFDEM import Survey
|
||||
from . import SrcFDEM as Src
|
||||
from . import RxFDEM as Rx
|
||||
from .ProblemFDEM import Problem3D_e, Problem3D_b, Problem3D_j, Problem3D_h
|
||||
from .FieldsFDEM import Fields3D_e, Fields3D_b, Fields3D_j, Fields3D_h
|
||||
|
||||
@@ -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
|
||||
|
||||
def getxBCyBC_CC(mesh, alpha, beta, gamma):
|
||||
|
||||
@@ -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
|
||||
from SimPEG.Utils import Identity, Zero
|
||||
import numpy as np
|
||||
@@ -9,7 +15,7 @@ class Fields(SimPEG.Problem.Fields):
|
||||
|
||||
def _phiDeriv(self, src, du_dm_v, v, adjoint=False):
|
||||
if getattr(self, '_phiDeriv_u', None) is None or getattr(self, '_phiDeriv_m', None) is None:
|
||||
raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %self.knownFields.keys()[0])
|
||||
raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
|
||||
|
||||
if 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):
|
||||
if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None:
|
||||
raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0])
|
||||
raise NotImplementedError ('Getting eDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
|
||||
|
||||
if 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):
|
||||
if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None:
|
||||
raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0])
|
||||
raise NotImplementedError ('Getting jDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
|
||||
|
||||
if adjoint:
|
||||
return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint)
|
||||
|
||||
@@ -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
|
||||
from SimPEG.Utils import Identity, Zero
|
||||
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):
|
||||
if getattr(self, '_phiDeriv_u', None) is None or getattr(self, '_phiDeriv_m', None) is None:
|
||||
raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %self.knownFields.keys()[0])
|
||||
raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
|
||||
|
||||
if 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):
|
||||
if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None:
|
||||
raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0])
|
||||
raise NotImplementedError ('Getting eDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
|
||||
|
||||
if 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):
|
||||
if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None:
|
||||
raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0])
|
||||
raise NotImplementedError ('Getting jDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
|
||||
|
||||
if adjoint:
|
||||
return self._jDeriv_u(kyInd, src, v, adjoint), self._jDeriv_m(kyInd, src, v, adjoint)
|
||||
|
||||
@@ -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.EM.Base import BaseEMProblem
|
||||
from SurveyDC import Survey
|
||||
from FieldsDC import Fields, Fields_CC, Fields_N
|
||||
from .SurveyDC import Survey
|
||||
from .FieldsDC import Fields, Fields_CC, Fields_N
|
||||
from SimPEG.Utils import sdiag
|
||||
import numpy as np
|
||||
from SimPEG.Utils import Zero
|
||||
from BoundaryUtils import getxBCyBC_CC
|
||||
from .BoundaryUtils import getxBCyBC_CC
|
||||
|
||||
class BaseDCProblem(BaseEMProblem):
|
||||
|
||||
|
||||
@@ -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.EM.Base import BaseEMProblem
|
||||
from SurveyDC import Survey, Survey_ky
|
||||
from FieldsDC_2D import Fields_ky, Fields_ky_CC, Fields_ky_N
|
||||
from .SurveyDC import Survey, Survey_ky
|
||||
from .FieldsDC_2D import Fields_ky, Fields_ky_CC, Fields_ky_N
|
||||
from SimPEG.Utils import sdiag
|
||||
import numpy as np
|
||||
from SimPEG.Utils import Zero
|
||||
from BoundaryUtils import getxBCyBC_CC
|
||||
from .BoundaryUtils import getxBCyBC_CC
|
||||
|
||||
class BaseDCProblem_2D(BaseEMProblem):
|
||||
|
||||
|
||||
@@ -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 numpy as np
|
||||
from SimPEG.Utils import Zero, closestPoints
|
||||
|
||||
@@ -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
|
||||
# from SimPEG.EM.Base import BaseEMSurvey
|
||||
from SimPEG.Utils import Zero, closestPoints, mkvc
|
||||
|
||||
@@ -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
|
||||
from SimPEG.EM.Base import BaseEMSurvey
|
||||
from SimPEG import sp, Survey
|
||||
from SimPEG.Utils import Zero, Identity
|
||||
from RxDC import BaseRx
|
||||
from SrcDC import BaseSrc
|
||||
from .RxDC import BaseRx
|
||||
from .SrcDC import BaseSrc
|
||||
|
||||
class Survey(BaseEMSurvey):
|
||||
rxPair = BaseRx
|
||||
|
||||
@@ -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
|
||||
|
||||
def WennerSrcList(nElecs, aSpacing, in2D=False, plotIt=False):
|
||||
|
||||
@@ -1,8 +1,14 @@
|
||||
from ProblemDC import Problem3D_CC, Problem3D_N
|
||||
from ProblemDC_2D import Problem2D_CC, Problem2D_N
|
||||
from SurveyDC import Survey, Survey_ky
|
||||
import SrcDC as Src #Pole
|
||||
import RxDC as Rx
|
||||
from FieldsDC import Fields_CC
|
||||
from BoundaryUtils import getxBCyBC_CC
|
||||
import Utils
|
||||
from __future__ import absolute_import
|
||||
from __future__ import unicode_literals
|
||||
from __future__ import print_function
|
||||
from __future__ import division
|
||||
from future import standard_library
|
||||
standard_library.install_aliases()
|
||||
from .ProblemDC import Problem3D_CC, Problem3D_N
|
||||
from .ProblemDC_2D import Problem2D_CC, Problem2D_N
|
||||
from .SurveyDC import Survey, Survey_ky
|
||||
from . import SrcDC as Src #Pole
|
||||
from . import RxDC as Rx
|
||||
from .FieldsDC import Fields_CC
|
||||
from .BoundaryUtils import getxBCyBC_CC
|
||||
from . import Utils
|
||||
|
||||
@@ -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.EM.Base import BaseEMProblem
|
||||
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
|
||||
from SimPEG.Utils import Zero
|
||||
from SimPEG.EM.Static.DC import getxBCyBC_CC
|
||||
from SurveyIP import Survey
|
||||
from .SurveyIP import Survey
|
||||
|
||||
class IPPropMap(Maps.PropMap):
|
||||
"""
|
||||
|
||||
@@ -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
|
||||
from SimPEG.EM.Base import BaseEMSurvey
|
||||
from SimPEG import sp, Survey
|
||||
|
||||
@@ -1,2 +1,8 @@
|
||||
from ProblemIP import Problem3D_CC, Problem3D_N
|
||||
from SurveyIP import Survey
|
||||
from __future__ import absolute_import
|
||||
from __future__ import unicode_literals
|
||||
from __future__ import print_function
|
||||
from __future__ import division
|
||||
from future import standard_library
|
||||
standard_library.install_aliases()
|
||||
from .ProblemIP import Problem3D_CC, Problem3D_N
|
||||
from .SurveyIP import Survey
|
||||
|
||||
@@ -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.EM.Base import BaseEMProblem
|
||||
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
|
||||
from SimPEG.Utils import Zero
|
||||
from SimPEG.EM.Static.DC import getxBCyBC_CC
|
||||
from SurveySIP import Survey, Data
|
||||
from .SurveySIP import Survey, Data
|
||||
|
||||
class ColeColePropMap(Maps.PropMap):
|
||||
"""
|
||||
@@ -105,7 +113,7 @@ class BaseSIPProblem(BaseEMProblem):
|
||||
JvAll = []
|
||||
#Assume only eta and tau (eta first then tau)
|
||||
# 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)):
|
||||
t = self.survey.times[tind]
|
||||
|
||||
@@ -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.Regularization import BaseRegularization, Simple
|
||||
|
||||
|
||||
@@ -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 numpy as np
|
||||
from SimPEG.Utils import Zero, closestPoints
|
||||
|
||||
@@ -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
|
||||
# from SimPEG.EM.Base import BaseEMSurvey
|
||||
from SimPEG.Utils import Zero, closestPoints, mkvc
|
||||
|
||||
@@ -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
|
||||
from SimPEG.EM.Base import BaseEMSurvey
|
||||
from SimPEG import np, sp, Survey, Utils
|
||||
|
||||
@@ -1,5 +1,11 @@
|
||||
from ProblemSIP import Problem3D_CC, Problem3D_N
|
||||
from SurveySIP import Survey, Data
|
||||
import SrcSIP as Src #Pole
|
||||
import RxSIP as Rx
|
||||
from Regularization import MultiRegularization
|
||||
from __future__ import absolute_import
|
||||
from __future__ import unicode_literals
|
||||
from __future__ import print_function
|
||||
from __future__ import division
|
||||
from future import standard_library
|
||||
standard_library.install_aliases()
|
||||
from .ProblemSIP import Problem3D_CC, Problem3D_N
|
||||
from .SurveySIP import Survey, Data
|
||||
from . import SrcSIP as Src #Pole
|
||||
from . import RxSIP as Rx
|
||||
from .Regularization import MultiRegularization
|
||||
|
||||
@@ -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.EM.Static import DC, IP
|
||||
|
||||
@@ -85,10 +93,10 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None):
|
||||
|
||||
elif stype == 'dpdp':
|
||||
|
||||
leg = data * 2*np.pi / ( 1/MA - 1/MB + 1/NB - 1/NA )
|
||||
LEG.append(1./(2*np.pi) *( 1/MA - 1/MB + 1/NB - 1/NA ))
|
||||
leg = data * 2*np.pi / (1/MA - 1/MB + 1/NB - 1/NA)
|
||||
LEG.append(1./(2*np.pi) * (1/MA - 1/MB + 1/NB - 1/NA))
|
||||
else:
|
||||
print """dtype must be 'pdp'(pole-dipole) | 'dpdp' (dipole-dipole) """
|
||||
print("""dtype must be 'pdp'(pole-dipole) | 'dpdp' (dipole-dipole) """)
|
||||
break
|
||||
|
||||
|
||||
@@ -103,7 +111,7 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None):
|
||||
rho = np.hstack([rho,leg])
|
||||
|
||||
else:
|
||||
print """dtype must be 'appr' | 'appc' | 'volt' """
|
||||
print("""dtype must be 'appr' | 'appc' | 'volt' """)
|
||||
break
|
||||
|
||||
|
||||
@@ -184,14 +192,14 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
|
||||
# Mesure survey length and direction
|
||||
dl_len = xy_2_r(endl[0,0],endl[1,0],endl[0,1],endl[1,1])
|
||||
|
||||
dl_x = ( endl[1,0] - endl[0,0] ) / dl_len
|
||||
dl_y = ( endl[1,1] - endl[0,1] ) / dl_len
|
||||
dl_x = (endl[1,0] - endl[0,0]) / dl_len
|
||||
dl_y = (endl[1,1] - endl[0,1]) / dl_len
|
||||
|
||||
nstn = np.floor( dl_len / a )
|
||||
nstn = np.floor(dl_len / a)
|
||||
|
||||
# Compute discrete pole location along line
|
||||
stn_x = endl[0,0] + np.array(range(int(nstn)))*dl_x*a
|
||||
stn_y = endl[0,1] + np.array(range(int(nstn)))*dl_y*a
|
||||
stn_x = endl[0,0] + np.array(list(range(int(nstn))))*dl_x*a
|
||||
stn_y = endl[0,1] + np.array(list(range(int(nstn))))*dl_y*a
|
||||
|
||||
if mesh.dim==2:
|
||||
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])
|
||||
|
||||
# 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
|
||||
if nstn <= 0:
|
||||
continue
|
||||
|
||||
# Compute discrete pole location along line
|
||||
stn_x = N[ii,0] + dl_x*b + np.array(range(int(nstn)))*dl_x*a
|
||||
stn_y = N[ii,1] + dl_y*b + np.array(range(int(nstn)))*dl_y*a
|
||||
stn_x = N[ii,0] + dl_x*b + np.array(list(range(int(nstn))))*dl_x*a
|
||||
stn_y = N[ii,1] + dl_y*b + np.array(list(range(int(nstn))))*dl_y*a
|
||||
|
||||
# Create receiver poles
|
||||
|
||||
@@ -275,17 +283,17 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
|
||||
max_y = endl[1,1] - dl_y * b
|
||||
|
||||
box_l = np.sqrt( (min_x - max_x)**2 + (min_y - max_y)**2 )
|
||||
box_w = box_l/2.
|
||||
box_w = box_l / 2.
|
||||
|
||||
nstn = np.floor( box_l / a )
|
||||
nstn = np.floor(box_l / a)
|
||||
|
||||
# Compute discrete pole location along line
|
||||
stn_x = min_x + np.array(range(int(nstn)))*dl_x*a
|
||||
stn_y = min_y + np.array(range(int(nstn)))*dl_y*a
|
||||
stn_x = min_x + np.array(list(range(int(nstn))))*dl_x*a
|
||||
stn_y = min_y + np.array(list(range(int(nstn))))*dl_y*a
|
||||
|
||||
# Define number of cross lines
|
||||
nlin = int(np.floor( box_w / a ))
|
||||
lind = range(-nlin,nlin+1)
|
||||
nlin = int(box_w // a)
|
||||
lind = list(range(-nlin,nlin+1))
|
||||
|
||||
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,:])
|
||||
SrcList.append(srcClass)
|
||||
else:
|
||||
print """stype must be either 'pdp', 'dpdp' or 'gradient'. """
|
||||
print("""stype must be either 'pdp', 'dpdp' or 'gradient'. """)
|
||||
|
||||
|
||||
return SrcList
|
||||
|
||||
@@ -1 +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 *
|
||||
|
||||
@@ -1,3 +1,9 @@
|
||||
import DC
|
||||
import IP
|
||||
import SIP
|
||||
from __future__ import absolute_import
|
||||
from __future__ import unicode_literals
|
||||
from __future__ import print_function
|
||||
from __future__ import division
|
||||
from future import standard_library
|
||||
standard_library.install_aliases()
|
||||
from . import DC
|
||||
from . import IP
|
||||
from . import SIP
|
||||
|
||||
+21
-14
@@ -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.Problem import BaseTimeProblem
|
||||
from SimPEG.EM.Utils import *
|
||||
@@ -47,7 +54,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
|
||||
self.waveformType = "GENERAL"
|
||||
|
||||
def fields(self, m):
|
||||
if self.verbose: print '%s\nCalculating fields(m)\n%s'%('*'*50,'*'*50)
|
||||
if self.verbose: print('%s\nCalculating fields(m)\n%s'%('*'*50,'*'*50))
|
||||
self.curModel = m
|
||||
# Create a fields storage object
|
||||
F = self._FieldsForward_pair(self.mesh, self.survey)
|
||||
@@ -55,7 +62,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
|
||||
# Set the initial conditions
|
||||
F[src,:,0] = src.getInitialFields(self.mesh)
|
||||
F = self.forward(m, self.getRHS, F=F)
|
||||
if self.verbose: print '%s\nDone calculating fields(m)\n%s'%('*'*50,'*'*50)
|
||||
if self.verbose: print('%s\nDone calculating fields(m)\n%s'%('*'*50,'*'*50))
|
||||
return F
|
||||
|
||||
def forward(self, m, RHS, F=None):
|
||||
@@ -70,13 +77,13 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
|
||||
if Ainv is not None:
|
||||
Ainv.clean()
|
||||
A = self.getA(tInd)
|
||||
if self.verbose: print 'Factoring... (dt = %e)'%dt
|
||||
if self.verbose: print('Factoring... (dt = %e)'%dt)
|
||||
Ainv = self.Solver(A, **self.solverOpts)
|
||||
if self.verbose: print 'Done'
|
||||
if self.verbose: print('Done')
|
||||
rhs = RHS(tInd, F)
|
||||
if self.verbose: print ' Solving... (tInd = %d)'%tInd
|
||||
if self.verbose: print(' Solving... (tInd = %d)'%tInd)
|
||||
sol = Ainv * rhs
|
||||
if self.verbose: print ' Done...'
|
||||
if self.verbose: print(' Done...')
|
||||
if sol.ndim == 1:
|
||||
sol.shape = (sol.size,1)
|
||||
F[:,self.solType,tInd+1] = sol
|
||||
@@ -95,13 +102,13 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
|
||||
if Ainv is not None:
|
||||
Ainv.clean()
|
||||
A = self.getA(tInd)
|
||||
if self.verbose: print 'Factoring (Adjoint)... (dt = %e)'%dt
|
||||
if self.verbose: print('Factoring (Adjoint)... (dt = %e)'%dt)
|
||||
Ainv = self.Solver(A, **self.solverOpts)
|
||||
if self.verbose: print 'Done'
|
||||
if self.verbose: print('Done')
|
||||
rhs = RHS(tInd, F)
|
||||
if self.verbose: print ' Solving (Adjoint)... (tInd = %d)'%tInd
|
||||
if self.verbose: print(' Solving (Adjoint)... (tInd = %d)'%tInd)
|
||||
sol = Ainv * rhs
|
||||
if self.verbose: print ' Done...'
|
||||
if self.verbose: print(' Done...')
|
||||
if sol.ndim == 1:
|
||||
sol.shape = (sol.size,1)
|
||||
F[:,self.solType,tInd+1] = sol
|
||||
@@ -123,14 +130,14 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
|
||||
* Compute \\\(\\\\vec{w} = -\\\mathbf{Q} \\\\vec{y}\\\)
|
||||
|
||||
"""
|
||||
if self.verbose: print '%s\nCalculating J(v)\n%s'%('*'*50,'*'*50)
|
||||
if self.verbose: print('%s\nCalculating J(v)\n%s'%('*'*50,'*'*50))
|
||||
self.curModel = m
|
||||
if f is None:
|
||||
f = self.fields(m)
|
||||
p = self.Gvec(m, v, f)
|
||||
y = self.solveAh(m, p)
|
||||
Jv = self.survey.evalDeriv(f, v=y)
|
||||
if self.verbose: print '%s\nDone calculating J(v)\n%s'%('*'*50,'*'*50)
|
||||
if self.verbose: print('%s\nDone calculating J(v)\n%s'%('*'*50,'*'*50))
|
||||
return - mkvc(Jv)
|
||||
|
||||
def Jtvec(self, m, v, f=None):
|
||||
@@ -148,7 +155,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
|
||||
* Compute \\\(\\\\vec{w} = -\\\mathbf{G}^\\\\top y\\\)
|
||||
|
||||
"""
|
||||
if self.verbose: print '%s\nCalculating J^T(v)\n%s'%('*'*50,'*'*50)
|
||||
if self.verbose: print('%s\nCalculating J^T(v)\n%s'%('*'*50,'*'*50))
|
||||
self.curModel = m
|
||||
if f is None:
|
||||
f = self.fields(m)
|
||||
@@ -159,6 +166,6 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
|
||||
p = self.survey.evalDeriv(f, v=v, adjoint=True)
|
||||
y = self.solveAht(m, p)
|
||||
w = self.Gtvec(m, y, f)
|
||||
if self.verbose: print '%s\nDone calculating J^T(v)\n%s'%('*'*50,'*'*50)
|
||||
if self.verbose: print('%s\nDone calculating J^T(v)\n%s'%('*'*50,'*'*50))
|
||||
return - mkvc(w)
|
||||
|
||||
|
||||
@@ -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.Survey import BaseSurvey
|
||||
from SimPEG.EM.Utils import *
|
||||
from BaseTDEM import FieldsTDEM
|
||||
from .BaseTDEM import FieldsTDEM
|
||||
|
||||
|
||||
class RxTDEM(Survey.BaseTimeRx):
|
||||
@@ -87,7 +93,7 @@ class SrcTDEM_VMD_MVP(SrcTDEM):
|
||||
def getInitialFields(self, mesh):
|
||||
"""Vertical magnetic dipole, magnetic vector potential"""
|
||||
if self.waveformType == "STEPOFF":
|
||||
print ">> Step waveform: Non-zero initial condition"
|
||||
print(">> Step waveform: Non-zero initial condition")
|
||||
if mesh._meshType is 'CYL':
|
||||
if mesh.isSymmetric:
|
||||
MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
|
||||
@@ -99,7 +105,7 @@ class SrcTDEM_VMD_MVP(SrcTDEM):
|
||||
raise Exception('Unknown mesh for VMD')
|
||||
return {"b": mesh.edgeCurl*MVP}
|
||||
elif self.waveformType == "GENERAL":
|
||||
print ">> General waveform: Zero initial condition"
|
||||
print(">> General waveform: Zero initial condition")
|
||||
return {"b": np.zeros(mesh.nF)}
|
||||
else:
|
||||
raise NotImplementedError("Only use STEPOFF or GENERAL")
|
||||
@@ -127,7 +133,7 @@ class SrcTDEM_CircularLoop_MVP(SrcTDEM):
|
||||
def getInitialFields(self, mesh):
|
||||
"""Circular Loop, magnetic vector potential"""
|
||||
if self.waveformType == "STEPOFF":
|
||||
print ">> Step waveform: Non-zero initial condition"
|
||||
print(">> Step waveform: Non-zero initial condition")
|
||||
if mesh._meshType is 'CYL':
|
||||
if mesh.isSymmetric:
|
||||
MVP = MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
|
||||
@@ -139,7 +145,7 @@ class SrcTDEM_CircularLoop_MVP(SrcTDEM):
|
||||
raise Exception('Unknown mesh for CircularLoop')
|
||||
return {"b": mesh.edgeCurl*MVP}
|
||||
elif self.waveformType == "GENERAL":
|
||||
print ">> General waveform: Zero initial condition"
|
||||
print(">> General waveform: Zero initial condition")
|
||||
return {"b": np.zeros(mesh.nF)}
|
||||
else:
|
||||
raise NotImplementedError("Only use STEPOFF or GENERAL")
|
||||
|
||||
@@ -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
|
||||
import numpy as np
|
||||
from SurveyTDEM import SurveyTDEM
|
||||
from .SurveyTDEM import SurveyTDEM
|
||||
|
||||
|
||||
class FieldsTDEM_e_from_b(FieldsTDEM):
|
||||
|
||||
@@ -1,3 +1,9 @@
|
||||
from SurveyTDEM import * #SurveyTDEM, RxTDEM, SrcTDEM
|
||||
from BaseTDEM import BaseTDEMProblem, FieldsTDEM
|
||||
from TDEM_b import ProblemTDEM_b
|
||||
from __future__ import absolute_import
|
||||
from __future__ import unicode_literals
|
||||
from __future__ import print_function
|
||||
from __future__ import division
|
||||
from future import standard_library
|
||||
standard_library.install_aliases()
|
||||
from .SurveyTDEM import * #SurveyTDEM, RxTDEM, SrcTDEM
|
||||
from .BaseTDEM import BaseTDEMProblem, FieldsTDEM
|
||||
from .TDEM_b import ProblemTDEM_b
|
||||
|
||||
@@ -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 scipy.special import ellipk, ellipe
|
||||
from scipy.constants import mu_0, pi
|
||||
@@ -17,7 +24,7 @@ def MagneticDipoleVectorPotential(srcLoc, obsLoc, component, moment=1., dipoleMo
|
||||
#TODO: break this out!
|
||||
|
||||
if type(component) in [list, tuple]:
|
||||
out = range(len(component))
|
||||
out = list(range(len(component)))
|
||||
for i, comp in enumerate(component):
|
||||
out[i] = MagneticDipoleVectorPotential(srcLoc, obsLoc, comp, dipoleMoment=dipoleMoment)
|
||||
return np.concatenate(out)
|
||||
@@ -118,7 +125,7 @@ def MagneticLoopVectorPotential(srcLoc, obsLoc, component, radius, mu=mu_0):
|
||||
"""
|
||||
|
||||
if type(component) in [list, tuple]:
|
||||
out = range(len(component))
|
||||
out = list(range(len(component)))
|
||||
for i, comp in enumerate(component):
|
||||
out[i] = MagneticLoopVectorPotential(srcLoc, obsLoc, comp, radius, mu)
|
||||
return np.concatenate(out)
|
||||
@@ -158,11 +165,11 @@ def MagneticLoopVectorPotential(srcLoc, obsLoc, component, radius, mu=mu_0):
|
||||
# % 1/r singular at r = 0 and K(m) singular at m = 1
|
||||
Aphi = np.zeros(n)
|
||||
# % Common factor is (mu * I) / pi with I = 1 and mu = 4e-7 * pi.
|
||||
Aphi[ind] = 4e-7 / np.sqrt(m[ind]) * np.sqrt(radius / r[ind]) *((1. - m[ind] / 2.) * K[ind] - E[ind])
|
||||
Aphi[ind] = 4e-7 / np.sqrt(m[ind]) * np.sqrt(radius/ r[ind]) *((1. - m[ind] / 2.) * K[ind] - E[ind])
|
||||
if component == 'x':
|
||||
A[ind, i] = Aphi[ind] * (-y[ind] / r[ind] )
|
||||
A[ind, i] = Aphi[ind] * (-y[ind] / r[ind])
|
||||
elif component == 'y':
|
||||
A[ind, i] = Aphi[ind] * ( x[ind] / r[ind] )
|
||||
A[ind, i] = Aphi[ind] * (x[ind] / r[ind])
|
||||
else:
|
||||
raise ValueError('Invalid component')
|
||||
|
||||
|
||||
@@ -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
|
||||
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):
|
||||
""" Eq 1.47 - 1.49 in Ward and Hohmann """
|
||||
w = omega(freq)
|
||||
alp = w * np.sqrt( mu*eps/2 * ( np.sqrt(1. + (sigma / (eps*w))**2 ) + 1) )
|
||||
beta = w * np.sqrt( mu*eps/2 * ( np.sqrt(1. + (sigma / (eps*w))**2 ) - 1) )
|
||||
alp = w * np.sqrt( mu*eps/2 * ( np.sqrt(1. + (sigma / (eps*w)))**2 ) + 1)
|
||||
beta = w * np.sqrt( mu*eps/2 * ( np.sqrt(1. + (sigma / (eps*w)))**2 ) - 1)
|
||||
return alp - 1j*beta
|
||||
|
||||
|
||||
|
||||
@@ -1,2 +1,8 @@
|
||||
from EMUtils import omega, k
|
||||
from AnalyticUtils import MagneticDipoleFields, MagneticDipoleVectorPotential, MagneticLoopVectorPotential
|
||||
from __future__ import absolute_import
|
||||
from __future__ import unicode_literals
|
||||
from __future__ import print_function
|
||||
from __future__ import division
|
||||
from future import standard_library
|
||||
standard_library.install_aliases()
|
||||
from .EMUtils import omega, k
|
||||
from .AnalyticUtils import MagneticDipoleFields, MagneticDipoleVectorPotential, MagneticLoopVectorPotential
|
||||
@@ -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
|
||||
from SimPEG import *
|
||||
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))
|
||||
|
||||
if verbose:
|
||||
print ' Fetching %s problem' % (fdemType)
|
||||
print(' Fetching %s problem' % (fdemType))
|
||||
|
||||
if fdemType == 'e':
|
||||
survey = EM.FDEM.Survey(Src)
|
||||
@@ -83,7 +90,7 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False):
|
||||
try:
|
||||
from pymatsolver import MumpsSolver
|
||||
prb.Solver = MumpsSolver
|
||||
except ImportError, e:
|
||||
except ImportError as e:
|
||||
prb.Solver = SolverLU
|
||||
|
||||
return prb
|
||||
@@ -94,7 +101,7 @@ def crossCheckTest(SrcList, fdemType1, fdemType2, comp, addrandoms = False, useM
|
||||
|
||||
prb1 = getFDEMProblem(fdemType1, comp, SrcList, freq, useMu, verbose)
|
||||
mesh = prb1.mesh
|
||||
print 'Cross Checking Forward: %s, %s formulations - %s' % (fdemType1, fdemType2, comp)
|
||||
print('Cross Checking Forward: %s, %s formulations - %s' % (fdemType1, fdemType2, comp))
|
||||
|
||||
logsig = np.log(np.ones(mesh.nC)*CONDUCTIVITY)
|
||||
mu = np.ones(mesh.nC)*MU
|
||||
@@ -112,7 +119,7 @@ def crossCheckTest(SrcList, fdemType1, fdemType2, comp, addrandoms = False, useM
|
||||
d1 = survey1.dpred(m)
|
||||
|
||||
if verbose:
|
||||
print ' Problem 1 solved'
|
||||
print(' Problem 1 solved')
|
||||
|
||||
|
||||
prb2 = getFDEMProblem(fdemType2, comp, SrcList, freq, useMu, verbose)
|
||||
@@ -121,11 +128,11 @@ def crossCheckTest(SrcList, fdemType1, fdemType2, comp, addrandoms = False, useM
|
||||
d2 = survey2.dpred(m)
|
||||
|
||||
if verbose:
|
||||
print ' Problem 2 solved'
|
||||
print(' Problem 2 solved')
|
||||
|
||||
r = d2-d1
|
||||
l2r = l2norm(r)
|
||||
|
||||
tol = np.max([TOL*(10**int(np.log10(0.5* (l2norm(d1) + l2norm(d2)) ))),FLR])
|
||||
print l2norm(d1), l2norm(d2), l2r , tol, l2r < tol
|
||||
print(l2norm(d1), l2norm(d2), l2r , tol, l2r < tol)
|
||||
return l2r < tol
|
||||
|
||||
+12
-6
@@ -1,7 +1,13 @@
|
||||
import TDEM
|
||||
import FDEM
|
||||
import Static
|
||||
import Base
|
||||
import Analytics
|
||||
import Utils
|
||||
from __future__ import absolute_import
|
||||
from __future__ import unicode_literals
|
||||
from __future__ import print_function
|
||||
from __future__ import division
|
||||
from future import standard_library
|
||||
standard_library.install_aliases()
|
||||
from . import TDEM
|
||||
from . import FDEM
|
||||
from . import Static
|
||||
from . import Base
|
||||
from . import Analytics
|
||||
from . import Utils
|
||||
from scipy.constants import mu_0, epsilon_0
|
||||
|
||||
@@ -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 *
|
||||
import SimPEG.EM.Static.DC as DC
|
||||
|
||||
@@ -29,7 +35,7 @@ def run(plotIt=True):
|
||||
try:
|
||||
from pymatsolver import MumpsSolver
|
||||
problem.Solver = MumpsSolver
|
||||
except Exception, e:
|
||||
except Exception as e:
|
||||
pass
|
||||
data = survey.dpred(sigma)
|
||||
|
||||
@@ -61,8 +67,8 @@ def run(plotIt=True):
|
||||
ax[0].set_title('Computed')
|
||||
plt.show()
|
||||
|
||||
return np.linalg.norm(data-data_ana)/np.linalg.norm(data_ana)
|
||||
return np.linalg.norm(data-data_ana) / np.linalg.norm(data_ana)
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
print run()
|
||||
print(run())
|
||||
|
||||
@@ -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
|
||||
import SimPEG.DCIP as DC
|
||||
import time
|
||||
@@ -57,7 +65,7 @@ def run(loc=None, sig=None, radi=None, param=None, surveyType='dipole-dipole', u
|
||||
model[ind] = sig[2]
|
||||
|
||||
# Get index of the center
|
||||
indy = int(mesh.nCy/2)
|
||||
indy = int(mesh.nCy // 2)
|
||||
|
||||
# Plot the model for reference
|
||||
# Define core mesh extent
|
||||
@@ -124,10 +132,10 @@ def run(loc=None, sig=None, radi=None, param=None, surveyType='dipole-dipole', u
|
||||
tx = np.squeeze(Tx[ii][:,0:1])
|
||||
tinf = tx + np.array([dl_x,dl_y,0])*dl_len*2
|
||||
inds = Utils.closestPoints(mesh, np.c_[tx,tinf].T)
|
||||
RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*( [-1] / mesh.vol[inds] )
|
||||
RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*([-1] / mesh.vol[inds])
|
||||
else:
|
||||
inds = Utils.closestPoints(mesh, np.asarray(Tx[ii]).T )
|
||||
RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*( [-1,1] / mesh.vol[inds] )
|
||||
RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*([-1,1] / mesh.vol[inds])
|
||||
|
||||
# Iterative Solve
|
||||
Ainvb = sp.linalg.bicgstab(P*A,P*RHS, tol=1e-5)
|
||||
@@ -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
|
||||
|
||||
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 'Forward completed'
|
||||
print('Transmitter {0} of {1}'.format(ii,len(Tx)))
|
||||
print('Forward completed')
|
||||
|
||||
# Let's just convert the 3D format into 2D (distance along line) and plot
|
||||
survey2D = DC.convertObs_DC3D_to_2D(survey, np.ones(survey.nSrc) , 'Xloc')
|
||||
|
||||
@@ -1,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 *
|
||||
import SimPEG.EM as EM
|
||||
from SimPEG.EM import mu_0
|
||||
@@ -56,7 +62,7 @@ def run(plotIt=True):
|
||||
try:
|
||||
from pymatsolver import MumpsSolver
|
||||
prb.Solver = MumpsSolver
|
||||
except ImportError, e:
|
||||
except ImportError as e:
|
||||
prb.Solver = SolverLU
|
||||
|
||||
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 *
|
||||
import SimPEG.EM as EM
|
||||
|
||||
|
||||
@@ -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.EM import FDEM, Analytics, mu_0
|
||||
import time
|
||||
@@ -67,8 +74,8 @@ def run(plotIt=True):
|
||||
casing_l = 300 # length of the casing
|
||||
|
||||
casing_r = 0.1
|
||||
casing_a = casing_r - casing_t/2. # inner radius
|
||||
casing_b = casing_r + casing_t/2. # outer radius
|
||||
casing_a = casing_r - casing_t / 2. # inner radius
|
||||
casing_b = casing_r + casing_t / 2. # outer radius
|
||||
casing_z = np.r_[-casing_l,0.]
|
||||
|
||||
|
||||
@@ -78,25 +85,25 @@ def run(plotIt=True):
|
||||
src_loc = np.r_[0.,0.,dsz]
|
||||
inf_loc = np.r_[0.,0.,1e4]
|
||||
|
||||
print 'Skin Depth: ', [(500./np.sqrt(sigmaback*_)) for _ in freqs]
|
||||
print('Skin Depth: ', [(500. / np.sqrt(sigmaback*_)) for _ in freqs])
|
||||
|
||||
|
||||
# ------------------ MESH ------------------
|
||||
# fine cells near well bore
|
||||
csx1, csx2 = 2e-3, 60.
|
||||
pfx1, pfx2 = 1.3, 1.3
|
||||
ncx1 = np.ceil(casing_b/csx1+2)
|
||||
ncx1 = np.ceil(casing_b/csx1)+2
|
||||
|
||||
# pad nicely to second cell size
|
||||
npadx1 = np.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)])
|
||||
dx1 = sum(hx1a)+sum(hx1b)
|
||||
dx1 = np.floor(dx1/csx2)
|
||||
hx1b *= (dx1*csx2 - sum(hx1a))/sum(hx1b)
|
||||
dx1 = dx1 // csx2
|
||||
hx1b *= (dx1*csx2 - sum(hx1a)) / sum(hx1b)
|
||||
|
||||
# second chunk of mesh
|
||||
dx2 = 300. # uniform mesh out to here
|
||||
ncx2 = np.ceil((dx2 - dx1)/csx2)
|
||||
ncx2 = np.ceil((dx2 - dx1) / csx2)
|
||||
npadx2 = 45
|
||||
hx2a, hx2b = Utils.meshTensor([(csx2,ncx2)]), Utils.meshTensor([(csx2,npadx2,pfx2)])
|
||||
hx = np.hstack([hx1a,hx1b,hx2a,hx2b])
|
||||
@@ -110,8 +117,8 @@ def run(plotIt=True):
|
||||
# Mesh
|
||||
mesh = Mesh.CylMesh([hx,1.,hz], [0.,0.,-np.sum(hz[:npadzu+ncz-nza])])
|
||||
|
||||
print 'Mesh Extent xmax: %f,: zmin: %f, zmax: %f'%(mesh.vectorCCx.max(), mesh.vectorCCz.min(), mesh.vectorCCz.max())
|
||||
print 'Number of cells', mesh.nC
|
||||
print('Mesh Extent xmax: %f,: zmin: %f, zmax: %f'%(mesh.vectorCCx.max(), mesh.vectorCCz.min(), mesh.vectorCCz.max()))
|
||||
print('Number of cells', mesh.nC)
|
||||
|
||||
if plotIt is True:
|
||||
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
|
||||
@@ -224,7 +231,7 @@ def run(plotIt=True):
|
||||
# ------------- Solve ---------------------------
|
||||
t0 = time.time()
|
||||
fieldsCasing = problem.fields(sigCasing)
|
||||
print 'Time to solve 2 sources', time.time() - t0
|
||||
print('Time to solve 2 sources', time.time() - t0)
|
||||
|
||||
# Plot current
|
||||
|
||||
@@ -251,9 +258,9 @@ def run(plotIt=True):
|
||||
in1_in = in1[np.r_[inds]]
|
||||
z_in = mesh.gridFz[inds_fz,2]
|
||||
|
||||
in0_in = in0_in.reshape([in0_in.shape[0]/3,3])
|
||||
in1_in = in1_in.reshape([in1_in.shape[0]/3,3])
|
||||
z_in = z_in.reshape([z_in.shape[0]/3,3])
|
||||
in0_in = in0_in.reshape([in0_in.shape[0]//3,3])
|
||||
in1_in = in1_in.reshape([in1_in.shape[0]//3,3])
|
||||
z_in = z_in.reshape([z_in.shape[0]//3,3])
|
||||
|
||||
I0 = in0_in.sum(1).real
|
||||
I1 = in1_in.sum(1).real
|
||||
|
||||
@@ -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 *
|
||||
import SimPEG.EM as EM
|
||||
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.FLOW import Richards
|
||||
|
||||
|
||||
@@ -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 *
|
||||
|
||||
|
||||
@@ -47,7 +55,7 @@ def run(N=100, plotIt=True):
|
||||
|
||||
# Distance weighting
|
||||
wr = np.sum(prob.G**2.,axis=0)**0.5
|
||||
wr = ( wr/np.max(wr) )
|
||||
wr = ( wr/np.max(wr))
|
||||
|
||||
dmis = DataMisfit.l2_DataMisfit(survey)
|
||||
dmis.Wd = 1./wd
|
||||
@@ -59,15 +67,15 @@ def run(N=100, plotIt=True):
|
||||
reg.cell_weights = wr
|
||||
|
||||
reg.mref = np.zeros(mesh.nC)
|
||||
|
||||
|
||||
|
||||
opt = Optimization.ProjectedGNCG(maxIter=100 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 10, tolCG = 1e-3)
|
||||
invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
|
||||
update_Jacobi = Directives.Update_lin_PreCond()
|
||||
|
||||
|
||||
# Set the IRLS directive, penalize the lowest 25 percentile of model values
|
||||
# Start with an l2-l2, then switch to lp-norms
|
||||
norms = [0., 0., 2., 2.]
|
||||
norms = [0., 0., 2., 2.]
|
||||
IRLS = Directives.Update_IRLS( norms=norms, prctile = 25, maxIRLSiter = 15, minGNiter=3)
|
||||
|
||||
inv = Inversion.BaseInversion(invProb, directiveList=[IRLS,betaest,update_Jacobi])
|
||||
@@ -75,7 +83,7 @@ def run(N=100, plotIt=True):
|
||||
# Run inversion
|
||||
mrec = inv.run(m0)
|
||||
|
||||
print "Final misfit:" + str(invProb.dmisfit.eval(mrec))
|
||||
print("Final misfit:" + str(invProb.dmisfit.eval(mrec)))
|
||||
|
||||
|
||||
if plotIt:
|
||||
|
||||
@@ -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 *
|
||||
|
||||
|
||||
|
||||
@@ -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 numpy as np
|
||||
import SimPEG.MT as MT
|
||||
|
||||
@@ -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.
|
||||
|
||||
# Import
|
||||
|
||||
@@ -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
|
||||
|
||||
def run(plotIt=True):
|
||||
|
||||
@@ -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
|
||||
|
||||
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 Mesh, Utils, np, SolverLU
|
||||
|
||||
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 *
|
||||
|
||||
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 *
|
||||
|
||||
def run(plotIt=True):
|
||||
|
||||
@@ -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 *
|
||||
|
||||
def run(plotIt=True, n=60):
|
||||
@@ -87,7 +95,7 @@ def run(plotIt=True, n=60):
|
||||
if elapsed > capture[jj]:
|
||||
PHIS += [(elapsed, phi.copy())]
|
||||
jj += 1
|
||||
if ii % 10 == 0: print ii, elapsed
|
||||
if ii % 10 == 0: print(ii, elapsed)
|
||||
ii += 1
|
||||
|
||||
if plotIt:
|
||||
|
||||
@@ -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 *
|
||||
|
||||
def run(plotIt=True):
|
||||
|
||||
@@ -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 *
|
||||
|
||||
def run(plotIt=True, n=60):
|
||||
@@ -28,15 +36,15 @@ def run(plotIt=True, n=60):
|
||||
axes[0].set_xlim([-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],'%d'%ii, color='r')
|
||||
|
||||
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],'%d'%ii, color='g')
|
||||
|
||||
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,'%d'%(ii+M.nFx), color='m')
|
||||
|
||||
axes[1].spy(M.faceDiv)
|
||||
|
||||
@@ -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 *
|
||||
|
||||
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 *
|
||||
|
||||
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.Utils import surface2ind_topo
|
||||
|
||||
|
||||
+30
-23
@@ -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.
|
||||
|
||||
##### AUTOIMPORTS #####
|
||||
import DC_Analytic_Dipole
|
||||
import DC_Forward_PseudoSection
|
||||
import EM_FDEM_1D_Inversion
|
||||
import EM_FDEM_Analytic_MagDipoleWholespace
|
||||
import EM_Schenkel_Morrison_Casing
|
||||
import EM_TDEM_1D_Inversion
|
||||
import FLOW_Richards_1D_Celia1990
|
||||
import Inversion_IRLS
|
||||
import Inversion_Linear
|
||||
import Maps_ComboMaps
|
||||
import Maps_Mesh2Mesh
|
||||
import Mesh_Basic_ForwardDC
|
||||
import Mesh_Basic_PlotImage
|
||||
import Mesh_Basic_Types
|
||||
import Mesh_Operators_CahnHilliard
|
||||
import Mesh_QuadTree_Creation
|
||||
import Mesh_QuadTree_FaceDiv
|
||||
import Mesh_QuadTree_HangingNodes
|
||||
import Mesh_Tensor_Creation
|
||||
import MT_1D_ForwardAndInversion
|
||||
import MT_3D_Foward
|
||||
import Utils_surface2ind_topo
|
||||
from . import DC_Analytic_Dipole
|
||||
from . import DC_Forward_PseudoSection
|
||||
from . import EM_FDEM_1D_Inversion
|
||||
from . import EM_FDEM_Analytic_MagDipoleWholespace
|
||||
from . import EM_Schenkel_Morrison_Casing
|
||||
from . import EM_TDEM_1D_Inversion
|
||||
from . import FLOW_Richards_1D_Celia1990
|
||||
from . import Inversion_IRLS
|
||||
from . import Inversion_Linear
|
||||
from . import Maps_ComboMaps
|
||||
from . import Maps_Mesh2Mesh
|
||||
from . import Mesh_Basic_ForwardDC
|
||||
from . import Mesh_Basic_PlotImage
|
||||
from . import Mesh_Basic_Types
|
||||
from . import Mesh_Operators_CahnHilliard
|
||||
from . import Mesh_QuadTree_Creation
|
||||
from . import Mesh_QuadTree_FaceDiv
|
||||
from . import Mesh_QuadTree_HangingNodes
|
||||
from . import Mesh_Tensor_Creation
|
||||
from . import MT_1D_ForwardAndInversion
|
||||
from . import MT_3D_Foward
|
||||
from . import Utils_surface2ind_topo
|
||||
|
||||
__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"]
|
||||
|
||||
@@ -104,7 +111,7 @@ if __name__ == '__main__':
|
||||
|
||||
rst = os.path.sep.join((filePath.split(os.path.sep)[:-3] + ['docs', 'content', 'examples', name + '.rst']))
|
||||
|
||||
print 'Creating: %s.rst'%name
|
||||
print('Creating: %s.rst'%name)
|
||||
f = open(rst, 'w')
|
||||
f.write(out)
|
||||
f.close()
|
||||
|
||||
@@ -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 future.utils import with_metaclass
|
||||
|
||||
|
||||
class NonLinearMap(object):
|
||||
class NonLinearMap(with_metaclass(Utils.SimPEGMetaClass, object)):
|
||||
"""
|
||||
SimPEG NonLinearMap
|
||||
|
||||
"""
|
||||
|
||||
__metaclass__ = Utils.SimPEGMetaClass
|
||||
|
||||
counter = None #: A SimPEG.Utils.Counter object
|
||||
mesh = None #: A SimPEG Mesh
|
||||
|
||||
@@ -186,7 +192,7 @@ class _haverkamp_theta(NonLinearMap):
|
||||
|
||||
def transformDerivU(self, u, m):
|
||||
self.setModel(m)
|
||||
g = (self.alpha*((self.theta_s - self.theta_r)/
|
||||
g = (self.alpha*((self.theta_s - self.theta_r) /
|
||||
(self.alpha + abs(u)**self.beta)**2)
|
||||
*(-self.beta*abs(u)**(self.beta-1)*np.sign(u)))
|
||||
g[u >= 0] = 0
|
||||
@@ -273,7 +279,7 @@ class _vangenuchten_theta(NonLinearMap):
|
||||
def transform(self, u, m):
|
||||
self.setModel(m)
|
||||
m = 1 - 1.0/self.n
|
||||
f = (( self.theta_s - self.theta_r )/
|
||||
f = (( self.theta_s - self.theta_r ) /
|
||||
((1+abs(self.alpha*u)**self.n)**m) + self.theta_r)
|
||||
if Utils.isScalar(self.theta_s):
|
||||
f[u >= 0] = self.theta_s
|
||||
@@ -343,7 +349,7 @@ class _vangenuchten_k(NonLinearMap):
|
||||
Ks = self.Ks
|
||||
m = 1.0 - 1.0/n
|
||||
|
||||
g = I*alpha*n*np.exp(Ks)*abs(alpha*u)**(n - 1.0)*np.sign(alpha*u)*(1.0/n - 1.0)*((abs(alpha*u)**n + 1)**(1.0/n - 1))**(I - 1)*((1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)**2*(abs(alpha*u)**n + 1)**(1.0/n - 2) - (2*alpha*n*np.exp(Ks)*abs(alpha*u)**(n - 1)*np.sign(alpha*u)*(1.0/n - 1)*((abs(alpha*u)**n + 1)**(1.0/n - 1))**I*((1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)*(abs(alpha*u)**n + 1)**(1.0/n - 2))/(((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1) + 1)*(1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1.0/n))
|
||||
g = I*alpha*n*np.exp(Ks)*abs(alpha*u)**(n - 1.0)*np.sign(alpha*u)*(1.0/n - 1.0)*((abs(alpha*u)**n + 1)**(1.0/n - 1))**(I - 1)*((1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)**2*(abs(alpha*u)**n + 1)**(1.0/n - 2) - (2*alpha*n*np.exp(Ks)*abs(alpha*u)**(n - 1)*np.sign(alpha*u)*(1.0/n - 1)*((abs(alpha*u)**n + 1)**(1.0/n - 1))**I*((1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)*(abs(alpha*u)**n + 1)**(1.0/n - 2))/(((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1))+ 1)*(1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1.0/n)
|
||||
g[u >= 0] = 0
|
||||
g = Utils.sdiag(g)
|
||||
return g
|
||||
|
||||
@@ -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 Empirical import RichardsMap
|
||||
from .Empirical import RichardsMap
|
||||
import time
|
||||
|
||||
|
||||
@@ -61,7 +68,7 @@ class RichardsSurvey(Survey.BaseSurvey):
|
||||
|
||||
@Utils.requires('prob')
|
||||
def eval(self, U, m):
|
||||
Ds = range(len(self.rxList))
|
||||
Ds = list(range(len(self.rxList)))
|
||||
for ii, rx in enumerate(self.rxList):
|
||||
Ds[ii] = rx.eval(U, m,
|
||||
self.prob.mapping,
|
||||
@@ -73,7 +80,7 @@ class RichardsSurvey(Survey.BaseSurvey):
|
||||
@Utils.requires('prob')
|
||||
def evalDeriv(self, U, m):
|
||||
"""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):
|
||||
Ds[ii] = rx.evalDeriv(U, m,
|
||||
self.prob.mapping,
|
||||
@@ -135,12 +142,12 @@ class RichardsProblem(Problem.BaseTimeProblem):
|
||||
@Utils.timeIt
|
||||
def fields(self, m):
|
||||
tic = time.time()
|
||||
u = range(self.nT+1)
|
||||
u = list(range(self.nT+1))
|
||||
u[0] = self.initialConditions
|
||||
for ii, dt in enumerate(self.timeSteps):
|
||||
bc = self.getBoundaryConditions(ii, u[ii])
|
||||
u[ii+1] = self.rootFinder.root(lambda hn1m, return_g=True: self.getResidual(m, u[ii], hn1m, dt, bc, return_g=return_g), u[ii])
|
||||
if self.debug: print "Solving Fields (%4d/%d - %3.1f%% Done) %d Iterations, %4.2f seconds"%(ii+1, self.nT, 100.0*(ii+1)/self.nT, self.rootFinder.iter, time.time() - tic)
|
||||
if self.debug: print("Solving Fields (%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
|
||||
|
||||
@Utils.timeIt
|
||||
@@ -238,7 +245,7 @@ class RichardsProblem(Problem.BaseTimeProblem):
|
||||
f = self.fields(m)
|
||||
|
||||
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):
|
||||
dt = self.timeSteps[ii]
|
||||
bc = self.getBoundaryConditions(ii, f[ii])
|
||||
@@ -263,7 +270,7 @@ class RichardsProblem(Problem.BaseTimeProblem):
|
||||
if f is None:
|
||||
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.
|
||||
bc = self.getBoundaryConditions(0, f[0])
|
||||
@@ -295,7 +302,7 @@ class RichardsProblem(Problem.BaseTimeProblem):
|
||||
bc = self.getBoundaryConditions(ii-1, f[ii-1])
|
||||
Asub, Adiag, B = self.diagsJacobian(m, f[ii-1], f[ii], self.timeSteps[ii-1], bc)
|
||||
#select the correct part of v
|
||||
vpart = 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)
|
||||
JTvC = AdiaginvT * (PTv[vpart] - minus)
|
||||
minus = Asub.T*JTvC # this is now the super diagonal.
|
||||
|
||||
@@ -1,2 +1,8 @@
|
||||
import Empirical
|
||||
from RichardsProblem 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 . import Empirical
|
||||
from .RichardsProblem import *
|
||||
|
||||
@@ -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
|
||||
|
||||
+11
-2
@@ -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):
|
||||
"""Fancy Field Storage
|
||||
@@ -244,7 +253,7 @@ class TimeFields(Fields):
|
||||
out = func(pointerFields, srcII, timeII)
|
||||
else: #loop over the time steps
|
||||
nT = pointerShape[2]
|
||||
out = range(nT)
|
||||
out = list(range(nT))
|
||||
for i, TIND_i in enumerate(timeII):
|
||||
fieldI = pointerFields[:,:,i]
|
||||
if fieldI.shape[0] == fieldI.size:
|
||||
|
||||
+19
-12
@@ -1,14 +1,21 @@
|
||||
import Utils, Survey, Problem, numpy as np, scipy.sparse as sp, gc
|
||||
from Utils.SolverUtils import *
|
||||
import DataMisfit
|
||||
import Regularization
|
||||
from __future__ import print_function
|
||||
from __future__ import absolute_import
|
||||
from __future__ import unicode_literals
|
||||
from __future__ import division
|
||||
from future import standard_library
|
||||
standard_library.install_aliases()
|
||||
from builtins import object
|
||||
from . import Utils, Survey, Problem
|
||||
import numpy as np, scipy.sparse as sp, gc
|
||||
from .Utils.SolverUtils import *
|
||||
from . import DataMisfit
|
||||
from . import Regularization
|
||||
from future.utils import with_metaclass
|
||||
|
||||
|
||||
class BaseInvProblem(object):
|
||||
class BaseInvProblem(with_metaclass(Utils.SimPEGMetaClass, object)):
|
||||
"""BaseInvProblem(dmisfit, reg, opt)"""
|
||||
|
||||
__metaclass__ = Utils.SimPEGMetaClass
|
||||
|
||||
beta = 1.0 #: Trade-off parameter
|
||||
|
||||
debug = False #: Print debugging information
|
||||
@@ -54,10 +61,10 @@ class BaseInvProblem(object):
|
||||
|
||||
Called when inversion is first starting.
|
||||
"""
|
||||
if self.debug: print 'Calling InvProblem.startup'
|
||||
if self.debug: print('Calling InvProblem.startup')
|
||||
|
||||
if self.reg.mref is None:
|
||||
print 'SimPEG.InvProblem will set Regularization.mref to m0.'
|
||||
print('SimPEG.InvProblem will set Regularization.mref to m0.')
|
||||
self.reg.mref = m0
|
||||
|
||||
self.phi_d = np.nan
|
||||
@@ -65,8 +72,8 @@ class BaseInvProblem(object):
|
||||
|
||||
self.curModel = m0
|
||||
|
||||
print """SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
|
||||
***Done using same Solver and solverOpts as the problem***"""
|
||||
print("""SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
|
||||
***Done using same Solver and solverOpts as the problem***""")
|
||||
self.opt.bfgsH0 = self.prob.Solver(self.reg.eval2Deriv(self.curModel), **self.prob.solverOpts)
|
||||
|
||||
@property
|
||||
@@ -87,7 +94,7 @@ class BaseInvProblem(object):
|
||||
for mtest, u_ofmtest in self.warmstart:
|
||||
if m is mtest:
|
||||
f = u_ofmtest
|
||||
if self.debug: print 'InvProb is Warm Starting!'
|
||||
if self.debug: print('InvProb is Warm Starting!')
|
||||
break
|
||||
|
||||
if f is None:
|
||||
|
||||
+11
-5
@@ -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
|
||||
from SimPEG import Utils, sp, np
|
||||
from Optimization import Remember, IterationPrinters, StoppingCriteria
|
||||
import Directives
|
||||
from .Optimization import Remember, IterationPrinters, StoppingCriteria
|
||||
from . import Directives
|
||||
from future.utils import with_metaclass
|
||||
|
||||
|
||||
class BaseInversion(object):
|
||||
class BaseInversion(with_metaclass(Utils.SimPEGMetaClass, object)):
|
||||
"""
|
||||
|
||||
Inversion Class.
|
||||
|
||||
"""
|
||||
|
||||
__metaclass__ = Utils.SimPEGMetaClass
|
||||
|
||||
name = 'BaseInversion'
|
||||
|
||||
debug = False #: Print debugging information
|
||||
|
||||
+8
-2
@@ -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.EM.FDEM.ProblemFDEM import BaseFDEMProblem
|
||||
from SurveyMT import Survey, Data
|
||||
from FieldsMT import BaseMTFields
|
||||
from .SurveyMT import Survey, Data
|
||||
from .FieldsMT import BaseMTFields
|
||||
|
||||
|
||||
class BaseMTProblem(BaseFDEMProblem):
|
||||
|
||||
+12
-6
@@ -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 scipy.constants import mu_0
|
||||
import sys
|
||||
@@ -63,7 +69,7 @@ class Fields1D_e(BaseMTFields):
|
||||
C = self.mesh.nodalGrad
|
||||
b = (C * eSolution)
|
||||
for i, src in enumerate(srcList):
|
||||
b[:,i] *= - 1./(1j*omega(src.freq))
|
||||
b[:,i] *= -1./(1j*omega(src.freq))
|
||||
# There is no magnetic source in the MT problem
|
||||
# S_m, _ = src.eval(self.survey.prob)
|
||||
# if S_m is not None:
|
||||
@@ -188,7 +194,7 @@ class Fields3D_e(BaseMTFields):
|
||||
# adjoint: returns a 2*nE long vector with zero's for py
|
||||
return np.vstack((v,np.zeros_like(v)))
|
||||
# Not adjoint: return only the px part of the vector
|
||||
return v[:len(v)/2]
|
||||
return v[:len(v)//2]
|
||||
|
||||
def _e_pyDeriv_u(self, src, v, adjoint = False):
|
||||
'''
|
||||
@@ -198,7 +204,7 @@ class Fields3D_e(BaseMTFields):
|
||||
# adjoint: returns a 2*nE long vector with zero's for px
|
||||
return np.vstack((np.zeros_like(v),v))
|
||||
# Not adjoint: return only the px part of the vector
|
||||
return v[len(v)/2::]
|
||||
return v[len(v)//2::]
|
||||
|
||||
def _e_pxDeriv_m(self, src, v, adjoint = False):
|
||||
# assuming primary does not depend on the model
|
||||
@@ -227,7 +233,7 @@ class Fields3D_e(BaseMTFields):
|
||||
C = self.mesh.edgeCurl
|
||||
b = (C * e_pxSolution)
|
||||
for i, src in enumerate(srcList):
|
||||
b[:,i] *= - 1./(1j*omega(src.freq))
|
||||
b[:,i] *= -1./(1j*omega(src.freq))
|
||||
# There is no magnetic source in the MT problem
|
||||
# S_m, _ = src.eval(self.survey.prob)
|
||||
# if S_m is not None:
|
||||
@@ -238,7 +244,7 @@ class Fields3D_e(BaseMTFields):
|
||||
C = self.mesh.edgeCurl
|
||||
b = (C * e_pySolution)
|
||||
for i, src in enumerate(srcList):
|
||||
b[:,i] *= - 1./(1j*omega(src.freq))
|
||||
b[:,i] *= -1./(1j*omega(src.freq))
|
||||
# There is no magnetic source in the MT problem
|
||||
# S_m, _ = src.eval(self.survey.prob)
|
||||
# if S_m is not None:
|
||||
@@ -348,4 +354,4 @@ class Fields3D_e(BaseMTFields):
|
||||
This function stacks the fields derivatives appropriately
|
||||
"""
|
||||
# The fields have no dependance to the model.
|
||||
return None
|
||||
return None
|
||||
|
||||
@@ -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 import mkvc
|
||||
from scipy.constants import mu_0
|
||||
@@ -46,7 +52,7 @@ class eForm_psField(BaseMTProblem):
|
||||
Edge inner product matrix
|
||||
"""
|
||||
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
|
||||
|
||||
@property
|
||||
@@ -142,7 +148,7 @@ class eForm_psField(BaseMTProblem):
|
||||
for freq in self.survey.freqs:
|
||||
if self.verbose:
|
||||
startTime = time.time()
|
||||
print 'Starting work for {:.3e}'.format(freq)
|
||||
print('Starting work for {:.3e}'.format(freq))
|
||||
sys.stdout.flush()
|
||||
A = self.getA(freq)
|
||||
rhs = self.getRHS(freq)
|
||||
@@ -158,7 +164,7 @@ class eForm_psField(BaseMTProblem):
|
||||
# b = -( self.mesh.nodalGrad * e )/( 1j*omega(freq) )
|
||||
# F[Src, 'b_1d'] = b[:,1]
|
||||
if self.verbose:
|
||||
print 'Ran for {:f} seconds'.format(time.time()-startTime)
|
||||
print('Ran for {:f} seconds'.format(time.time()-startTime))
|
||||
sys.stdout.flush()
|
||||
return F
|
||||
|
||||
@@ -191,7 +197,7 @@ class eForm_TotalField(BaseMTProblem):
|
||||
Edge inner product matrix
|
||||
"""
|
||||
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
|
||||
|
||||
@property
|
||||
@@ -249,7 +255,7 @@ class eForm_TotalField(BaseMTProblem):
|
||||
Ed, Eu, Hd, Hu = getEHfields(self.mesh,self.curModel.sigma,freq,self.mesh.vectorNx)
|
||||
Etot = (Ed + Eu)
|
||||
sourceAmp = 1.0
|
||||
Etot = ((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
|
||||
eBC = np.r_[Etot[0],Etot[-1]]
|
||||
# The right hand side
|
||||
@@ -274,7 +280,7 @@ class eForm_TotalField(BaseMTProblem):
|
||||
for freq in self.survey.freqs:
|
||||
if self.verbose:
|
||||
startTime = time.time()
|
||||
print 'Starting work for {:.3e}'.format(freq)
|
||||
print('Starting work for {:.3e}'.format(freq))
|
||||
sys.stdout.flush()
|
||||
A = self.getA(freq)
|
||||
rhs, e_o = self.getRHS(freq)
|
||||
@@ -286,6 +292,6 @@ class eForm_TotalField(BaseMTProblem):
|
||||
# NOTE: only store e fields
|
||||
F[Src, 'e_1dSolution'] = e[:,0]
|
||||
if self.verbose:
|
||||
print 'Ran for {:f} seconds'.format(time.time()-startTime)
|
||||
print('Ran for {:f} seconds'.format(time.time()-startTime))
|
||||
sys.stdout.flush()
|
||||
return F
|
||||
|
||||
@@ -1 +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
|
||||
@@ -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
|
||||
@@ -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.EM.Utils import omega
|
||||
from scipy.constants import mu_0
|
||||
@@ -115,7 +121,7 @@ class eForm_ps(BaseMTProblem):
|
||||
for freq in self.survey.freqs:
|
||||
if self.verbose:
|
||||
startTime = time.time()
|
||||
print 'Starting work for {:.3e}'.format(freq)
|
||||
print('Starting work for {:.3e}'.format(freq))
|
||||
sys.stdout.flush()
|
||||
A = self.getA(freq)
|
||||
rhs = self.getRHS(freq)
|
||||
@@ -131,7 +137,7 @@ class eForm_ps(BaseMTProblem):
|
||||
# Note curl e = -iwb so b = -curl/iw
|
||||
|
||||
if self.verbose:
|
||||
print 'Ran for {:f} seconds'.format(time.time()-startTime)
|
||||
print('Ran for {:f} seconds'.format(time.time()-startTime))
|
||||
sys.stdout.flush()
|
||||
Ainv.clean()
|
||||
return F
|
||||
|
||||
@@ -1 +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
|
||||
+10
-4
@@ -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.EM.FDEM.SrcFDEM import BaseSrc as FDEMBaseSrc
|
||||
from SimPEG.EM.Utils import omega
|
||||
from scipy.constants import mu_0
|
||||
from numpy.lib import recfunctions as recFunc
|
||||
from Utils.sourceUtils import homo1DModelSource
|
||||
from Utils import rec2ndarr
|
||||
from .Utils.sourceUtils import homo1DModelSource
|
||||
from .Utils import rec2ndarr
|
||||
import sys
|
||||
|
||||
#################
|
||||
@@ -78,7 +84,7 @@ class polxy_1Dprimary(BaseMTSrc):
|
||||
C = problem.mesh.nodalGrad
|
||||
elif problem.mesh.dim == 3:
|
||||
C = problem.mesh.edgeCurl
|
||||
bBG_bp = (- C * self.ePrimary(problem) )*(1/( 1j*omega(self.freq) ))
|
||||
bBG_bp = (- C * self.ePrimary(problem) )*(1/(1j*omega(self.freq)))
|
||||
return bBG_bp
|
||||
|
||||
def S_e(self,problem):
|
||||
@@ -155,7 +161,7 @@ class polxy_3Dprimary(BaseMTSrc):
|
||||
C = problem.mesh.nodalGrad
|
||||
elif problem.mesh.dim == 3:
|
||||
C = problem.mesh.edgeCurl
|
||||
bBG_bp = (- C * self.ePrimary(problem) )*(1/( 1j*omega(self.freq) ))
|
||||
bBG_bp = (- C * self.ePrimary(problem) )*(1/(1j*omega(self.freq)))
|
||||
return bBG_bp
|
||||
|
||||
def S_e(self,problem):
|
||||
|
||||
+27
-21
@@ -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.EM.FDEM.SrcFDEM import BaseSrc as FDEMBaseSrc
|
||||
from SimPEG.EM.Utils import omega
|
||||
from scipy.constants import mu_0
|
||||
from numpy.lib import recfunctions as recFunc
|
||||
from Utils import rec2ndarr
|
||||
import SrcMT
|
||||
from .Utils import rec2ndarr
|
||||
from . import SrcMT
|
||||
import sys
|
||||
|
||||
#################
|
||||
@@ -76,7 +82,7 @@ class Rx(SimPEGsurvey.BaseRx):
|
||||
bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0
|
||||
# Note: Has a minus sign in front, to comply with quadrant calculations.
|
||||
# Can be derived from zyx case for the 3D case.
|
||||
f_part_complex = -ex/bx
|
||||
f_part_complex = old_div(-ex,bx)
|
||||
# elif self.projType is 'Z2D':
|
||||
elif self.projType is 'Z3D':
|
||||
## NOTE: Assumes that e is on edges and b on the faces. Need to generalize that or use a prop of fields to determine that.
|
||||
@@ -103,13 +109,13 @@ class Rx(SimPEGsurvey.BaseRx):
|
||||
hy_py = Pby*f[src,'b_py']/mu_0
|
||||
# Make the complex data
|
||||
if 'zxx' in self.rxType:
|
||||
f_part_complex = ( 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:
|
||||
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:
|
||||
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:
|
||||
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':
|
||||
if self.locs.ndim == 3:
|
||||
horLoc = self.locs[:,:,0]
|
||||
@@ -127,9 +133,9 @@ class Rx(SimPEGsurvey.BaseRx):
|
||||
by_py = Pby*f[src,'b_py']
|
||||
bz_py = Pbz*f[src,'b_py']
|
||||
if 'tzx' in self.rxType:
|
||||
f_part_complex = (- 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:
|
||||
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:
|
||||
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')
|
||||
# ex = Pex*mkvc(f[src,'e_1d'],2)
|
||||
# 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_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_de = -mkvc(Utils.sdiag(old_div(1.,(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)))*(Pex*v),2)
|
||||
dP_db = mkvc( Utils.sdiag(Pex*mkvc(f[src,'e_1d'],2))*(Utils.sdiag(old_div(1.,(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))).T*Utils.sdiag(old_div(1.,(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))))*(Pbx*f._bDeriv_u(src,v)/mu_0),2)
|
||||
PDeriv_complex = np.sum(np.hstack((dP_de,dP_db)),1)
|
||||
elif self.projType is 'Z2D':
|
||||
raise NotImplementedError('Has not been implement for 2D impedance tensor')
|
||||
@@ -198,7 +204,7 @@ class Rx(SimPEGsurvey.BaseRx):
|
||||
# Update the input vector
|
||||
sDiag = lambda t: Utils.sdiag(mkvc(t,2))
|
||||
# Define the components of the derivative
|
||||
Hd = sDiag(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)
|
||||
# Calculate components
|
||||
if 'zxx' in self.rxType:
|
||||
@@ -247,7 +253,7 @@ class Rx(SimPEGsurvey.BaseRx):
|
||||
# Update the input vector
|
||||
sDiag = lambda t: Utils.sdiag(mkvc(t,2))
|
||||
# Define the components of the derivative
|
||||
Hd = sDiag(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)
|
||||
if 'tzx' in self.rxType:
|
||||
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')
|
||||
# ex = Pex*mkvc(f[src,'e_1d'],2)
|
||||
# bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0
|
||||
dP_deTv = -mkvc(Pex.T*Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)).T*v,2)
|
||||
db_duv = Pbx.T/mu_0*Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))*(Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))).T*Utils.sdiag(Pex*mkvc(f[src,'e_1d'],2)).T*v
|
||||
dP_deTv = -mkvc(Pex.T*Utils.sdiag(old_div(1.,(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))).T*v,2)
|
||||
db_duv = Pbx.T/mu_0*Utils.sdiag(old_div(1.,(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)))*(Utils.sdiag(old_div(1.,(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)))).T*Utils.sdiag(Pex*mkvc(f[src,'e_1d'],2)).T*v
|
||||
dP_dbTv = mkvc(f._bDeriv_u(src,db_duv,adjoint=True),2)
|
||||
PDeriv_real = np.sum(np.hstack((dP_deTv,dP_dbTv)),1)
|
||||
elif self.projType is 'Z2D':
|
||||
@@ -300,17 +306,17 @@ class Rx(SimPEGsurvey.BaseRx):
|
||||
aey_px_u = lambda vec: f._e_pxDeriv_u(src,Pey.T*vec,adjoint=True)
|
||||
aex_py_u = lambda vec: f._e_pyDeriv_u(src,Pex.T*vec,adjoint=True)
|
||||
aey_py_u = lambda vec: f._e_pyDeriv_u(src,Pey.T*vec,adjoint=True)
|
||||
ahx_px_u = lambda vec: f._b_pxDeriv_u(src,Pbx.T*vec,adjoint=True)/mu_0
|
||||
ahy_px_u = lambda vec: f._b_pxDeriv_u(src,Pby.T*vec,adjoint=True)/mu_0
|
||||
ahx_py_u = lambda vec: f._b_pyDeriv_u(src,Pbx.T*vec,adjoint=True)/mu_0
|
||||
ahy_py_u = lambda vec: f._b_pyDeriv_u(src,Pby.T*vec,adjoint=True)/mu_0
|
||||
ahx_px_u = lambda vec: old_div(f._b_pxDeriv_u(src,Pbx.T*vec,adjoint=True),mu_0)
|
||||
ahy_px_u = lambda vec: old_div(f._b_pxDeriv_u(src,Pby.T*vec,adjoint=True),mu_0)
|
||||
ahx_py_u = lambda vec: old_div(f._b_pyDeriv_u(src,Pbx.T*vec,adjoint=True),mu_0)
|
||||
ahy_py_u = lambda vec: old_div(f._b_pyDeriv_u(src,Pby.T*vec,adjoint=True),mu_0)
|
||||
|
||||
# Update the input vector
|
||||
# Define shortcuts
|
||||
sDiag = lambda t: Utils.sdiag(mkvc(t,2))
|
||||
sVec = lambda t: Utils.sp.csr_matrix(mkvc(t,2))
|
||||
# Define the components of the derivative
|
||||
aHd = sDiag(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)
|
||||
# Need to fix this to reflect the adjoint
|
||||
if 'zxx' in self.rxType:
|
||||
@@ -362,7 +368,7 @@ class Rx(SimPEGsurvey.BaseRx):
|
||||
sDiag = lambda t: Utils.sdiag(mkvc(t,2))
|
||||
sVec = lambda t: Utils.sp.csr_matrix(mkvc(t,2))
|
||||
# Define the components of the derivative
|
||||
aHd = sDiag(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)
|
||||
# Need to fix this to reflect the adjoint
|
||||
if 'tzx' in self.rxType:
|
||||
|
||||
@@ -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
|
||||
|
||||
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
|
||||
for lnr, h in enumerate(m1d.hx): # lnr-number of layer, h-thickness of the layer
|
||||
# Calculate
|
||||
yp1 = k[lnr]/(w*mu[lnr]) # Admittance of the layer below the current layer
|
||||
zp = (w*mu[lnr+1])/k[lnr+1] # Impedance in the current layer
|
||||
yp1 = old_div(k[lnr],(w*mu[lnr])) # Admittance of the layer below the current layer
|
||||
zp = old_div((w*mu[lnr+1]),k[lnr+1]) # Impedance in the current layer
|
||||
# Build the propagation matrix
|
||||
|
||||
# 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])
|
||||
|
||||
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
|
||||
Ed = np.empty((zd.size,),dtype=complex)
|
||||
@@ -62,14 +69,14 @@ def getEHfields(m1d,sigma,freq,zd,scaleUD=True):
|
||||
dind = dup >= zd
|
||||
Ed[dind] = UDp[1,0]*np.exp(-1j*k[0]*(dup-zd[dind]))
|
||||
Eu[dind] = UDp[0,0]*np.exp(1j*k[0]*(dup-zd[dind]))
|
||||
Hd[dind] = (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]))
|
||||
Hd[dind] = (old_div(k[0],(w*mu[0])))*UDp[1,0]*np.exp(-1j*k[0]*(dup-zd[dind]))
|
||||
Hu[dind] = -(old_div(k[0],(w*mu[0])))*UDp[0,0]*np.exp(1j*k[0]*(dup-zd[dind]))
|
||||
for ki,mui,epsi,dlow,dup,Up,Dp in zip(k[1::],mu[1::],eps[1::],m1d.vectorNx[:-1],m1d.vectorNx[1::],UDp[0,1::],UDp[1,1::]):
|
||||
dind = np.logical_and(dup >= zd, zd > dlow)
|
||||
Ed[dind] = Dp*np.exp(-1j*ki*(dup-zd[dind]))
|
||||
Eu[dind] = Up*np.exp(1j*ki*(dup-zd[dind]))
|
||||
Hd[dind] = (ki/(w*mui))*Dp*np.exp(-1j*ki*(dup-zd[dind]))
|
||||
Hu[dind] = -(ki/(w*mui))*Up*np.exp(1j*ki*(dup-zd[dind]))
|
||||
Hd[dind] = (old_div(ki,(w*mui)))*Dp*np.exp(-1j*ki*(dup-zd[dind]))
|
||||
Hu[dind] = -(old_div(ki,(w*mui)))*Up*np.exp(1j*ki*(dup-zd[dind]))
|
||||
|
||||
# Return return the fields
|
||||
return Ed, Eu, Hd, Hu
|
||||
@@ -92,15 +99,15 @@ def getImpedance(m1d,sigma,freq):
|
||||
om = 2*np.pi*fr
|
||||
Zall = np.empty(len(h)+1,dtype='complex')
|
||||
# Calculate the impedance for the bottom layer
|
||||
Zall[0] = (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):
|
||||
# Calculate the wave number
|
||||
# print nr,sigma[nr]
|
||||
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()
|
||||
Z1d[nrFr] = Zall[-1]
|
||||
|
||||
@@ -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
|
||||
from MT1Danalytic import getEHfields
|
||||
from .MT1Danalytic import getEHfields
|
||||
from scipy.constants import mu_0
|
||||
|
||||
def get1DEfields(m1d,sigma,freq,sourceAmp=1.0):
|
||||
@@ -9,7 +15,7 @@ def get1DEfields(m1d,sigma,freq,sourceAmp=1.0):
|
||||
G = m1d.nodalGrad
|
||||
# Mass matrices
|
||||
# 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
|
||||
Msig = m1d.getFaceInnerProduct(sigma)
|
||||
# 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)
|
||||
Etot = (Ed + Eu)
|
||||
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
|
||||
bc = np.r_[Etot[0],Etot[-1]]
|
||||
# The right hand side
|
||||
|
||||
@@ -1,4 +1,10 @@
|
||||
from MT1Dsolutions import * # Add the names of the functions
|
||||
from MT1Danalytic import *
|
||||
from dataUtils import *
|
||||
from ediFilesUtils 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 .MT1Dsolutions import * # Add the names of the functions
|
||||
from .MT1Danalytic import *
|
||||
from .dataUtils import *
|
||||
from .ediFilesUtils import *
|
||||
|
||||
@@ -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,
|
||||
import numpy as np, matplotlib.pyplot as plt, sys
|
||||
import SimPEG as simpeg
|
||||
@@ -45,13 +51,13 @@ def rotateData(MTdata, rotAngle):
|
||||
|
||||
|
||||
def appResPhs(freq, z):
|
||||
app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
|
||||
app_phs = np.arctan2(z.imag,z.real)*(180/np.pi)
|
||||
app_res = (old_div((old_div(1.,(8e-7*np.pi**2))),freq))*np.abs(z)**2
|
||||
app_phs = np.arctan2(z.imag,z.real)*(old_div(180,np.pi))
|
||||
return app_res, app_phs
|
||||
|
||||
def skindepth(rho, freq):
|
||||
''' 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):
|
||||
return x.view((dt, len(x.dtype.names)))
|
||||
@@ -64,7 +70,7 @@ def makeAnalyticSolution(mesh, model, elev, freqs):
|
||||
anaE = anaEd+anaEu
|
||||
anaH = anaHd+anaHu
|
||||
|
||||
anaZ = anaE/anaH
|
||||
anaZ = old_div(anaE,anaH)
|
||||
# Add to the list
|
||||
data1D.append((freq,0,0,elev,anaZ[0]))
|
||||
dataRec = np.array(data1D,dtype=[('freq',float),('x',float),('y',float),('z',float),('zyx',complex)])
|
||||
@@ -97,7 +103,7 @@ def plotMT1DModelData(problem, models, symList=None):
|
||||
|
||||
# if not symList:
|
||||
# symList = ['x']*len(models)
|
||||
import plotDataTypes as pDt
|
||||
from . import plotDataTypes as pDt
|
||||
# Loop through the models.
|
||||
modelList = [problem.survey.mtrue]
|
||||
modelList.extend(models)
|
||||
@@ -110,14 +116,14 @@ def plotMT1DModelData(problem, models, symList=None):
|
||||
else:
|
||||
data1D = problem.dataPair(problem.survey,problem.survey.dpred(model)).toRecArray('Complex')
|
||||
# Plot the data and the model
|
||||
colRat = nr/((len(modelList)-1.999)*1.)
|
||||
colRat = old_div(nr,((len(modelList)-1.999)*1.))
|
||||
if colRat > 1.:
|
||||
col = 'k'
|
||||
else:
|
||||
col = plt.cm.seismic(1-colRat)
|
||||
# The model - make the pts to plot
|
||||
meshPts = np.concatenate((problem.mesh.gridN[0:1],np.kron(problem.mesh.gridN[1::],np.ones(2))[:-1]))
|
||||
modelPts = np.kron(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)
|
||||
|
||||
## Data
|
||||
@@ -144,7 +150,7 @@ def plotMT1DModelData(problem, models, symList=None):
|
||||
|
||||
# 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)
|
||||
[ l.set_rotation(90) for l in axM.get_yticklabels()]
|
||||
[ l.set_rotation(90) for l in axR.get_yticklabels()]
|
||||
@@ -157,7 +163,7 @@ def plotMT1DModelData(problem, models, symList=None):
|
||||
|
||||
def printTime():
|
||||
import time
|
||||
print time.strftime("%a, %d %b %Y %H:%M:%S +0000", time.localtime())
|
||||
print(time.strftime("%a, %d %b %Y %H:%M:%S +0000", time.localtime()))
|
||||
|
||||
def convert3Dto1Dobject(MTdata,rxType3D='zyx'):
|
||||
from SimPEG import MT
|
||||
|
||||
@@ -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.
|
||||
from SimPEG import mkvc
|
||||
from scipy.constants import mu_0
|
||||
@@ -9,7 +18,7 @@ import numpy as np
|
||||
import os, sys, re
|
||||
|
||||
|
||||
class EDIimporter:
|
||||
class EDIimporter(object):
|
||||
"""
|
||||
A class to import EDIfiles.
|
||||
|
||||
@@ -18,7 +27,7 @@ class EDIimporter:
|
||||
|
||||
# Define data converters
|
||||
_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
|
||||
filesList = None
|
||||
@@ -116,7 +125,7 @@ class EDIimporter:
|
||||
try:
|
||||
import osr
|
||||
except ImportError as e:
|
||||
print 'Could not import osr, missing the gdal package\nCan not project coordinates'
|
||||
print('Could not import osr, missing the gdal package\nCan not project coordinates')
|
||||
raise e
|
||||
# Coordinates convertor
|
||||
if self._2out is None:
|
||||
@@ -126,7 +135,7 @@ class EDIimporter:
|
||||
if self._outEPSG is None:
|
||||
# Find the UTM EPSG number
|
||||
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
|
||||
out.ImportFromEPSG(self._outEPSG)
|
||||
self._2out = osr.CoordinateTransformation(src,out)
|
||||
|
||||
@@ -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
|
||||
|
||||
|
||||
@@ -35,7 +41,7 @@ def plotIsoFreqNSimpedance(ax,freq,array,flag,par='abs',colorbar=True,colorNorm=
|
||||
cmap = plt.get_cmap('RdYlBu')
|
||||
level = np.concatenate((-np.logspace(0,-10,31),np.logspace(-10,0,31)))
|
||||
clevel = np.concatenate((-np.logspace(0,-8,5),np.logspace(-8,0,5)))
|
||||
plotNorm = colors.SymLogNorm(1e-10,linscale=2)
|
||||
plotNorm = colors.SymLogNorm(1e-10,linscale=2)
|
||||
if cLevel:
|
||||
level = np.concatenate((-np.logspace(0,-10,31),np.logspace(-10,0,31)))
|
||||
clevel = np.concatenate((-np.logspace(0,-8,5),np.logspace(-8,0,5)))
|
||||
@@ -65,45 +71,45 @@ def plotIsoFreqNSDiff(ax,freq,arrayList,flag,par='abs',colorbar=True,cLevel=True
|
||||
x, y = arrayList[0]['x'][indUniFreq0],arrayList[0]['y'][indUniFreq0]
|
||||
if par == 'abs':
|
||||
if useLog:
|
||||
zPlot = (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:
|
||||
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:
|
||||
maskInd = np.logical_or(np.abs(arrayList[0][flag][indUniFreq0])< 1e-3,np.abs(arrayList[1][flag][indUniFreq1]) < 1e-3)
|
||||
zPlot = np.ma.array(zPlot)
|
||||
zPlot[maskInd] = mask
|
||||
if cLevel:
|
||||
level = np.arange(-200,201,10)
|
||||
level = np.arange(-200,201,10)
|
||||
clevel = np.arange(-200,201,25)
|
||||
else:
|
||||
level = np.linspace(zPlot.min(),zPlot.max(),100)
|
||||
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
|
||||
elif par == 'real':
|
||||
if useLog:
|
||||
zPlot = (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:
|
||||
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:
|
||||
maskInd = np.logical_or(np.abs(np.real(arrayList[0][flag][indUniFreq0])) < 1e-3,np.abs(np.real(arrayList[1][flag][indUniFreq1])) < 1e-3)
|
||||
zPlot = np.ma.array(zPlot)
|
||||
zPlot[maskInd] = mask
|
||||
if cLevel:
|
||||
level = np.arange(-200,201,10)
|
||||
level = np.arange(-200,201,10)
|
||||
clevel = np.arange(-200,201,25)
|
||||
else:
|
||||
level = np.linspace(zPlot.min(),zPlot.max(),100)
|
||||
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
|
||||
elif par == 'imag':
|
||||
if useLog:
|
||||
zPlot = (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:
|
||||
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:
|
||||
maskInd = np.logical_or(np.abs(np.imag(arrayList[0][flag][indUniFreq0])) < 1e-3,np.abs(np.imag(arrayList[1][flag][indUniFreq1])) < 1e-3)
|
||||
zPlot = np.ma.array(zPlot)
|
||||
zPlot[maskInd] = mask
|
||||
if cLevel:
|
||||
level = np.arange(-200,201,10)
|
||||
level = np.arange(-200,201,10)
|
||||
clevel = np.arange(-200,201,25)
|
||||
else:
|
||||
level = np.linspace(zPlot.min(),zPlot.max(),100)
|
||||
@@ -174,13 +180,13 @@ def plotIsoFreqNStipper(ax,freq,array,flag,par='abs',colorbar=True,colorNorm='Sy
|
||||
plt.colorbar(cs,cax=ax.cax,ticks=clevel,format='%1.2e')
|
||||
ax.set_title(flag+' '+par,fontsize=8)
|
||||
|
||||
def plotIsoStaImpedance(ax,loc,array,flag,par='abs',pSym='s',pColor=None):
|
||||
|
||||
appResFact = 1/(8*np.pi**2*10**(-7))
|
||||
def plotIsoStaImpedance(ax,loc,array,flag,par='abs',pSym='s',pColor=None):
|
||||
|
||||
appResFact = old_div(1,(8*np.pi**2*10**(-7)))
|
||||
treshold = 1.0 # 1 meter
|
||||
indUniSta = np.sqrt(np.sum((rec2nd(array[['x','y']])-loc)**2,axis=1)) < treshold
|
||||
freq = array['freq'][indUniSta]
|
||||
|
||||
|
||||
if par == 'abs':
|
||||
zPlot = np.abs(array[flag][indUniSta])
|
||||
elif par == 'real':
|
||||
@@ -188,9 +194,9 @@ def plotIsoStaImpedance(ax,loc,array,flag,par='abs',pSym='s',pColor=None):
|
||||
elif par == 'imag':
|
||||
zPlot = np.imag(array[flag][indUniSta])
|
||||
elif par == 'res':
|
||||
zPlot = (appResFact/freq)*np.abs(array[flag][indUniSta])**2
|
||||
zPlot = (old_div(appResFact,freq))*np.abs(array[flag][indUniSta])**2
|
||||
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 'xx' in flag:
|
||||
@@ -198,7 +204,7 @@ def plotIsoStaImpedance(ax,loc,array,flag,par='abs',pSym='s',pColor=None):
|
||||
pColor = 'g'
|
||||
elif 'xy' in flag:
|
||||
lab = 'XY'
|
||||
pColor = 'r'
|
||||
pColor = 'r'
|
||||
elif 'yx' in flag:
|
||||
lab = 'YX'
|
||||
pColor = 'b'
|
||||
@@ -207,14 +213,14 @@ def plotIsoStaImpedance(ax,loc,array,flag,par='abs',pSym='s',pColor=None):
|
||||
pColor = 'y'
|
||||
|
||||
ax.plot(freq,zPlot,color=pColor,marker=pSym,label=flag)
|
||||
|
||||
|
||||
|
||||
def plotPsudoSectNSimpedance(ax,sectDict,array,flag,par='abs',colorbar=True,colorNorm='None',cLevel=None,contour=True):
|
||||
|
||||
indSect = np.where(sectDict.values()[0]==array[sectDict.keys()[0]])
|
||||
indSect = np.where(list(sectDict.values())[0]==array[list(sectDict.keys())[0]])
|
||||
|
||||
# Define the plot axes
|
||||
if 'x' in sectDict.keys()[0]:
|
||||
if 'x' in list(sectDict.keys())[0]:
|
||||
x = array['y'][indSect]
|
||||
else:
|
||||
x = array['x'][indSect]
|
||||
@@ -229,9 +235,9 @@ def plotPsudoSectNSimpedance(ax,sectDict,array,flag,par='abs',colorbar=True,colo
|
||||
else:
|
||||
level = np.linspace(zPlot.min(),zPlot.max(),100,endpoint=True)
|
||||
clevel = np.linspace(zPlot.min(),zPlot.max(),10,endpoint=True)
|
||||
|
||||
elif par == 'ares':
|
||||
zPlot = np.abs(array[flag][indSect])**2/(8*np.pi**2*10**(-7)*array['freq'][indSect])
|
||||
|
||||
elif par == 'ares':
|
||||
zPlot = old_div(np.abs(array[flag][indSect])**2,(8*np.pi**2*10**(-7)*array['freq'][indSect]))
|
||||
cmap = plt.get_cmap('RdYlBu')#seismic)
|
||||
if cLevel:
|
||||
zMax = np.log10(cLevel[1])
|
||||
@@ -244,7 +250,7 @@ def plotPsudoSectNSimpedance(ax,sectDict,array,flag,par='abs',colorbar=True,colo
|
||||
plotNorm = colors.LogNorm()
|
||||
|
||||
elif par == 'aphs':
|
||||
zPlot = np.arctan2(array[flag][indSect].imag,array[flag][indSect].real)*(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)
|
||||
if cLevel:
|
||||
zMax = cLevel[1]
|
||||
@@ -271,7 +277,7 @@ def plotPsudoSectNSimpedance(ax,sectDict,array,flag,par='abs',colorbar=True,colo
|
||||
elif par == 'imag':
|
||||
zPlot = np.imag(array[flag][indSect])
|
||||
cmap = plt.get_cmap('Spectral') #('RdYlBu')
|
||||
|
||||
|
||||
if cLevel:
|
||||
zMax = np.log10(cLevel[1])
|
||||
zMin = np.log10(cLevel[0])
|
||||
@@ -281,7 +287,7 @@ def plotPsudoSectNSimpedance(ax,sectDict,array,flag,par='abs',colorbar=True,colo
|
||||
level = np.concatenate((-np.logspace(zMax,zMin-.125,(zMax-zMin)*8+1,endpoint=True),np.logspace(zMin-.125,zMax,(zMax-zMin)*8+1,endpoint=True)))
|
||||
clevel = np.concatenate((-np.logspace(zMax,zMin,(zMax-zMin)*1+1,endpoint=True),np.logspace(zMin,zMax,(zMax-zMin)*1+1,endpoint=True)))
|
||||
plotNorm = colors.SymLogNorm(np.abs(level).min(),linscale=0.1)
|
||||
|
||||
|
||||
if colorNorm=='SymLog':
|
||||
plotNorm = colors.SymLogNorm(np.abs(level).min(),linscale=0.1)
|
||||
elif colorNorm=='Lin':
|
||||
@@ -307,14 +313,14 @@ def plotPsudoSectNSDiff(ax,sectDict,arrayList,flag,par='abs',colorbar=True,color
|
||||
def sortInArr(arr):
|
||||
return np.sort(arr,order=['freq','x','y','z'])
|
||||
# Find the index for the slice
|
||||
indSect0 = np.where(sectDict.values()[0]==arrayList[0][sectDict.keys()[0]])
|
||||
indSect1 = np.where(sectDict.values()[0]==arrayList[1][sectDict.keys()[0]])
|
||||
indSect0 = np.where(list(sectDict.values())[0]==arrayList[0][list(sectDict.keys())[0]])
|
||||
indSect1 = np.where(list(sectDict.values())[0]==arrayList[1][list(sectDict.keys())[0]])
|
||||
# Extract and sort the mats
|
||||
arr0 = sortInArr(arrayList[0][indSect0])
|
||||
arr1 = sortInArr(arrayList[1][indSect1])
|
||||
|
||||
# Define the plot axes
|
||||
if 'x' in sectDict.keys()[0]:
|
||||
if 'x' in list(sectDict.keys())[0]:
|
||||
x0 = arr0['y']
|
||||
x1 = arr1['y']
|
||||
else:
|
||||
@@ -322,35 +328,35 @@ def plotPsudoSectNSDiff(ax,sectDict,arrayList,flag,par='abs',colorbar=True,color
|
||||
x1 = arr1['x']
|
||||
y0 = arr0['freq']
|
||||
y1 = arr1['freq']
|
||||
|
||||
|
||||
|
||||
if par == 'abs':
|
||||
if useLog:
|
||||
zPlot = (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:
|
||||
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:
|
||||
maskInd = np.logical_or(np.abs(arr0[flag])< 1e-3,np.abs(arr1[flag]) < 1e-3)
|
||||
zPlot = np.ma.array(zPlot)
|
||||
zPlot[maskInd] = mask
|
||||
cmap = plt.get_cmap('RdYlBu')#seismic)
|
||||
elif par == 'ares':
|
||||
arF = 1/(8*np.pi**2*10**(-7))
|
||||
arF = old_div(1,(8*np.pi**2*10**(-7)))
|
||||
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:
|
||||
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:
|
||||
maskInd = np.logical_or(np.abs(arr0[flag])< 1e-3,np.abs(arr1[flag]) < 1e-3)
|
||||
zPlot = np.ma.array(zPlot)
|
||||
zPlot[maskInd] = mask
|
||||
cmap = plt.get_cmap('Spectral')#seismic)
|
||||
|
||||
|
||||
elif par == 'aphs':
|
||||
if useLog:
|
||||
zPlot = (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:
|
||||
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:
|
||||
maskInd = np.logical_or(np.abs(arr0[flag])< 1e-3,np.abs(arr1[flag]) < 1e-3)
|
||||
zPlot = np.ma.array(zPlot)
|
||||
@@ -358,26 +364,26 @@ def plotPsudoSectNSDiff(ax,sectDict,arrayList,flag,par='abs',colorbar=True,color
|
||||
cmap = plt.get_cmap('Spectral')#seismic)
|
||||
elif par == 'real':
|
||||
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:
|
||||
zPlot = (arr0[flag].real - arr1[flag].real)/arr1[flag].real
|
||||
zPlot = old_div((arr0[flag].real - arr1[flag].real),arr1[flag].real)
|
||||
if mask:
|
||||
maskInd = np.logical_or(arr0[flag].real< 1e-3,arr1[flag].real < 1e-3)
|
||||
zPlot = np.ma.array(zPlot)
|
||||
zPlot[maskInd] = mask
|
||||
cmap = plt.get_cmap('Spectral') #('Spectral')
|
||||
|
||||
|
||||
elif par == 'imag':
|
||||
if useLog:
|
||||
zPlot = (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:
|
||||
zPlot = (arr0[flag].imag - arr1[flag].imag)/arr1[flag].imag
|
||||
zPlot = old_div((arr0[flag].imag - arr1[flag].imag),arr1[flag].imag)
|
||||
if mask:
|
||||
maskInd = np.logical_or(arr0[flag].imag< 1e-3,arr1[flag].imag < 1e-3)
|
||||
zPlot = np.ma.array(zPlot)
|
||||
zPlot[maskInd] = mask
|
||||
cmap = plt.get_cmap('Spectral') #('RdYlBu')
|
||||
|
||||
|
||||
if cLevel:
|
||||
zMax = np.log10(cLevel[1])
|
||||
zMin = np.log10(cLevel[0])
|
||||
@@ -385,18 +391,18 @@ def plotPsudoSectNSDiff(ax,sectDict,arrayList,flag,par='abs',colorbar=True,color
|
||||
zMax = (np.ceil(np.log10(np.abs(zPlot).max())))
|
||||
zMin = (np.floor(np.log10(np.abs(zPlot).min())))
|
||||
|
||||
|
||||
|
||||
if colorNorm=='SymLog':
|
||||
level = np.concatenate((-np.logspace(zMax,zMin-.125,(zMax-zMin)*8+1,endpoint=True),np.logspace(zMin-.125,zMax,(zMax-zMin)*8+1,endpoint=True)))
|
||||
clevel = np.concatenate((-np.logspace(zMax,zMin,(zMax-zMin)*1+1,endpoint=True),np.logspace(zMin,zMax,(zMax-zMin)*1+1,endpoint=True)))
|
||||
plotNorm = colors.SymLogNorm(np.abs(level).min(),linscale=0.1)
|
||||
elif colorNorm=='Lin':
|
||||
if cLevel:
|
||||
level = np.arange(cLevel[0],cLevel[1]+.1,(cLevel[1] - cLevel[0])/50.)
|
||||
clevel = np.arange(cLevel[0],cLevel[1]+.1,(cLevel[1] - cLevel[0])/10.)
|
||||
level = np.arange(cLevel[0],cLevel[1]+.1,old_div((cLevel[1] - cLevel[0]),50.))
|
||||
clevel = np.arange(cLevel[0],cLevel[1]+.1,old_div((cLevel[1] - cLevel[0]),10.))
|
||||
else:
|
||||
level = np.arange(zPlot.min(),zPlot.max(),(zPlot.max() - zPlot.min())/50.)
|
||||
clevel = np.arange(zPlot.min(),zPlot.max(),(zPlot.max() - zPlot.min())/10.)
|
||||
level = np.arange(zPlot.min(),zPlot.max(),old_div((zPlot.max() - zPlot.min()),50.))
|
||||
clevel = np.arange(zPlot.min(),zPlot.max(),old_div((zPlot.max() - zPlot.min()),10.))
|
||||
plotNorm = colors.Normalize()
|
||||
elif colorNorm=='Log':
|
||||
level = np.logspace(zMin-.125,zMax,(zMax-zMin)*8+1,endpoint=True)
|
||||
|
||||
@@ -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
|
||||
|
||||
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)
|
||||
# Make the fields into a dictionary of location and the fields
|
||||
e0_1d = Eu+Ed
|
||||
E1dFieldDict = dict(zip(mesh.vectorNz,e0_1d))
|
||||
E1dFieldDict = dict(list(zip(mesh.vectorNz,e0_1d)))
|
||||
if mesh.dim == 1:
|
||||
eBG_px = simpeg.mkvc(e0_1d,2)
|
||||
eBG_py = -simpeg.mkvc(e0_1d,2) # added a minus to make the results in the correct quadrents.
|
||||
|
||||
@@ -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
|
||||
|
||||
def homo1DModelSource(mesh,freq,m_back):
|
||||
|
||||
+11
-5
@@ -1,5 +1,11 @@
|
||||
import Utils
|
||||
from SurveyMT import Rx, Survey, Data
|
||||
from FieldsMT import Fields1D_e, Fields3D_e
|
||||
import Problem1D, Problem2D, Problem3D
|
||||
import SrcMT
|
||||
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 Utils
|
||||
from .SurveyMT import Rx, Survey, Data
|
||||
from .FieldsMT import Fields1D_e, Fields3D_e
|
||||
from . import Problem1D, Problem2D, Problem3D
|
||||
from . import SrcMT
|
||||
Some files were not shown because too many files have changed in this diff Show More
Reference in New Issue
Block a user