Compare commits

..
Author SHA1 Message Date
GudniRos 2fcdabf3d5 Fixed directive to save iteration dictionary. 2016-06-09 12:39:25 -07:00
GudniRos b965c96242 Updating examples by running python __init__.py 2016-06-07 13:48:29 -07:00
GudniRos f06ba238f8 Fixing Examples/__init__.py after merge conflict 2016-06-07 13:14:36 -07:00
GudniRos 11408a3788 Fixing an FDEM import to be more explicit 2016-06-07 11:52:29 -07:00
GudniRos 358e2f96af Pull request #328 on github.
Fixing MT namespace to NSEM.
2016-06-07 10:54:16 -07:00
GudniRos cf5181016f Made MT_1d_analytic_nLayer_Earth.run have default n layers of 3.
Trying to fix example testing.
2016-06-07 09:45:02 -07:00
GudniRos d1f7d0c37a Adding ipywidgets to the conda install list for travis 2016-06-01 21:54:19 -07:00
GudniRos 50c4a6caf8 Refactor and updated tests 2016-06-01 20:53:12 -07:00
GudniRos b3029b697b Adjust settings for the MT1D inversion example 2016-06-01 12:00:52 -07:00
GudniRos 8d4581e92f Cleaned MT_1D inverison example 2016-06-01 11:55:12 -07:00
Gudni Karl Rosenkjaer 0a44796d4c Merge pull request #329 from simpeg/em/dev
Em/dev into mt/dev
2016-06-01 11:17:19 -07:00
GudniRos 8117ee3b06 Fixing MT_1D example, runs but inversion results should be better. 2016-06-01 10:22:18 -07:00
GudniRos fb7f5a53d4 Working on examples 2016-06-01 01:05:21 -07:00
GudniRos f34314bba2 Fixed tests and Jvec 2016-05-31 22:33:57 -07:00
GudniRos 7087632701 Fixing name space, updating utils 2016-05-31 22:06:47 -07:00
GudniRos 6165e619ed Moving Problem1D/3D folders to NSEM file
Looking at making Jvec more general.
2016-05-10 16:07:45 -07:00
GudniRos 1ab67b5790 Updated MT-->NSEM for all the classes.
Changed Vertical1DMap --> SurjectVertical1D
2016-05-10 14:41:20 -07:00
Gudni Karl Rosenkjaer a96e9e08d7 Merge pull request #316 from simpeg/em/dev
Em/dev into mt/NSEMrefact
2016-05-10 14:09:13 -07:00
GudniRos b1569e5734 Changes to MT1D analytic example 2016-05-10 14:00:48 -07:00
GudniRos 22f0a742e7 Fixed UBC mesh read in function 2016-05-02 05:16:56 -07:00
GudniRos 18357a11da Fixing analytic function 2016-04-15 14:48:58 -07:00
GudniRos abc5d72725 Merge branch 'em/dev' into mt/dev 2016-04-15 13:29:39 -07:00
GudniRos 6482b94cf1 Merge branch 'master' into mt/dev 2016-04-05 10:26:03 -07:00
GudniRos 009806f2ba Merge remote-tracking branch 'origin/dev' into mt/dev 2016-03-29 16:52:01 -07:00
GudniRos 127c51974f Fixed errors in analytic solution 2016-03-10 08:10:21 -08:00
GudniRos e96945b991 Indexing 2016-03-09 14:30:39 -08:00
GudniRos fc444a345f Fixing dtypes in the MT1Danalytic 2016-03-09 14:24:15 -08:00
GudniRos 9a739d8380 Merge branch 'mt/dev' of https://github.com/simpeg/simpeg into mt/dev 2016-03-09 13:27:50 -08:00
GudniRos 664a7bb484 Updated plotting functions for MT 2016-03-09 13:26:50 -08:00
151 changed files with 2691 additions and 3227 deletions
+1 -1
View File
@@ -1,4 +1,4 @@
[bumpversion]
current_version = 0.1.12
current_version = 0.1.10
files = setup.py SimPEG/__init__.py docs/conf.py
-2
View File
@@ -39,5 +39,3 @@ nosetests.xml
*.sublime-workspace
docs/_build/
Makefile
docs/warnings.txt
.DS_Store
+3 -26
View File
@@ -24,25 +24,18 @@ env:
- TEST_DIR=tests/examples
- TEST_DIR=tests/em/fdem/inverse/adjoint
- TEST_DIR=tests/em/fdem/forward
- TEST_DIR=tests/docs;
GAE_PYTHONPATH=${HOME}/.cache/google_appengine;
PATH=$PATH:${HOME}/google-cloud-sdk/bin;
PYTHONPATH=${PYTHONPATH}:${GAE_PYTHONPATH};
CLOUDSDK_CORE_DISABLE_PROMPTS=1
# Setup anaconda
before_install:
# Install packages
- if [ ${TRAVIS_PYTHON_VERSION:0:1} == "2" ]; then wget http://repo.continuum.io/miniconda/Miniconda-3.8.3-Linux-x86_64.sh
-O miniconda.sh; else wget http://repo.continuum.io/miniconda/Miniconda3-3.8.3-Linux-x86_64.sh
-O miniconda.sh; fi
- if [ ${TRAVIS_PYTHON_VERSION:0:1} == "2" ]; then wget http://repo.continuum.io/miniconda/Miniconda-3.8.3-Linux-x86_64.sh -O miniconda.sh; else wget http://repo.continuum.io/miniconda/Miniconda3-3.8.3-Linux-x86_64.sh -O miniconda.sh; fi
- chmod +x miniconda.sh
- ./miniconda.sh -b
- export PATH=/home/travis/anaconda/bin:/home/travis/miniconda/bin:$PATH
- conda update --yes conda
# Install packages
install:
- conda install --yes pip python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib cython ipython nose vtk sphinx
- conda install --yes pip python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib cython ipython ipywidgets nose vtk
- pip install nose-cov python-coveralls
- git clone https://github.com/rowanc1/pymatsolver.git
@@ -53,28 +46,12 @@ install:
# Run test
script:
# test docs
- nosetests $TEST_DIR --with-cov --cov SimPEG --cov-config .coveragerc -v -s
# Calculate coverage
after_success:
- coveralls --config_file .coveragerc
- if [ "$TRAVIS_BRANCH" = "master" -a "$TRAVIS_PULL_REQUEST" = "false" ]; then
if [ ${TEST_DIR} == "tests/docs" ]; then
python scripts/fetch_gae_sdk.py $(dirname "${GAE_PYTHONPATH}");
openssl aes-256-cbc -K $encrypted_93066031461c_key -iv $encrypted_93066031461c_iv
-in docs/credentials.tar.gz.enc -out credentials.tar.gz -d ;
if [ ! -d ${HOME}/google-cloud-sdk ]; then curl https://sdk.cloud.google.com | bash; fi ;
tar -xzf credentials.tar.gz ;
gcloud auth activate-service-account --key-file client-secret.json ;
gcloud config set project simpegdocs;
gcloud -q components update gae-python;
gcloud -q preview app deploy ./docs/app.yaml --version ${TRAVIS_COMMIT} --promote;
fi;
fi
notifications:
email:
- rowanc1@gmail.com
+1 -1
View File
@@ -1,4 +1,4 @@
.. image:: https://raw.github.com/simpeg/simpeg/master/docs/images/simpeg-logo.png
.. image:: https://raw.github.com/simpeg/simpeg/master/docs/simpeg-logo.png
:alt: SimPEG Logo
======
+2 -2
View File
@@ -162,8 +162,8 @@ class ProblemDC_CC(Problem.BaseProblem):
"""
Makes the matrix A(m) for the DC resistivity problem.
:param numpy.ndarray m: model
:rtype: scipy.sparse.csc_matrix
:param numpy.array m: model
:rtype: scipy.csc_matrix
:return: A(m)
.. math::
+1 -1
View File
@@ -71,7 +71,7 @@ class ProblemIP(Problem.BaseProblem):
Makes the matrix A(m) for the DC resistivity problem.
:param numpy.array m: model
:rtype: scipy.sparse.csc_matrix
:rtype: scipy.csc_matrix
:return: A(m)
.. math::
+36 -58
View File
@@ -167,7 +167,7 @@ class TargetMisfit(InversionDirective):
class SaveEveryIteration(InversionDirective):
class _SaveEveryIteration(InversionDirective):
@property
def name(self):
if getattr(self, '_name', None) is None:
@@ -188,7 +188,7 @@ class SaveEveryIteration(InversionDirective):
self._fileName = value
class SaveModelEveryIteration(SaveEveryIteration):
class SaveModelEveryIteration(_SaveEveryIteration):
"""SaveModelEveryIteration"""
def initialize(self):
@@ -198,7 +198,7 @@ class SaveModelEveryIteration(SaveEveryIteration):
np.save('%03d-%s' % (self.opt.iter, self.fileName), self.opt.xc)
class SaveOutputEveryIteration(SaveEveryIteration):
class SaveOutputEveryIteration(_SaveEveryIteration):
"""SaveModelEveryIteration"""
def initialize(self):
@@ -212,48 +212,40 @@ class SaveOutputEveryIteration(SaveEveryIteration):
f.write(' %3d %1.4e %1.4e %1.4e %1.4e\n'%(self.opt.iter, self.invProb.beta, self.invProb.phi_d, self.invProb.phi_m, self.opt.f))
f.close()
class SaveOutputDictEveryIteration(SaveEveryIteration):
"""SaveOutputDictEveryIteration"""
class SaveOutputDictEveryIteration(_SaveEveryIteration):
"""
Saves inversion parameters at every iteraion.
"""
def initialize(self):
print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '###-%s.npz'"%self.fileName
def endIter(self):
# Save the data.
ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.reg.mref) )
phi_ms = 0.5*ms.dot(ms)
if self.reg.mrefInSmooth == True:
mref = self.reg.mref
else:
mref = 0
mx = self.reg.Wx * ( self.reg.mapping * (self.invProb.curModel - mref) )
phi_mx = 0.5 * mx.dot(mx)
if self.prob.mesh.dim >= 2:
my = self.reg.Wy * ( self.reg.mapping * (self.invProb.curModel - mref) )
phi_my = 0.5 * my.dot(my)
else:
phi_my = 'NaN'
if self.prob.mesh.dim==3:
mz = self.reg.Wz * ( self.reg.mapping * (self.invProb.curModel - mref) )
phi_mz = 0.5 * mz.dot(mz)
else:
phi_mz = 'NaN'
# Initialize the output dict
outDict = {}
# Save the data.
outDict['iter'] = self.opt.iter
outDict['beta'] = self.invProb.beta
outDict['phi_d'] = self.invProb.phi_d
outDict['phi_ms'] = self.reg._evalSmall(self.invProb.curModel)
outDict['phi_mx'] = self.reg._evalSmoothx(self.invProb.curModel)
outDict['phi_my'] = self.reg._evalSmoothy(self.invProb.curModel) if self.prob.mesh.dim >= 2 else 'NaN'
outDict['phi_mz'] = self.reg._evalSmoothz(self.invProb.curModel) if self.prob.mesh.dim==3 else 'NaN'
outDict['f'] = self.opt.f
outDict['m'] = self.invProb.curModel
outDict['dpred'] = self.invProb.dpred
# Save the file as a npz
np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred)
# mref = getattr(self, 'm_prev', None)
# if mref is None:
# if self.debug: print 'UpdateReferenceModel is using mref0'
# mref = self.mref0
# self.m_prev = self.invProb.m_current
# return mref
np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), outDict)
class Update_IRLS(InversionDirective):
eps_min = None
eps = None
eps_p = None
eps_q = None
norms = [2.,2.,2.,2.]
factor = None
gamma = None
@@ -262,7 +254,6 @@ class Update_IRLS(InversionDirective):
f_old = None
f_min_change = 1e-2
beta_tol = 5e-2
prctile = 95
# Solving parameter for IRLS (mode:2)
IRLSiter = 0
@@ -297,22 +288,9 @@ class Update_IRLS(InversionDirective):
print "Convergence with smooth l2-norm regularization: Start IRLS steps..."
self.mode = 2
# Either use the supplied epsilon, or fix base on distribution of
# model values
if getattr(self, 'reg.eps', None) is None:
self.reg.eps_p = np.percentile(np.abs(self.invProb.curModel),self.prctile)
else:
self.reg.eps_p = self.eps[0]
if getattr(self, 'reg.eps', None) is None:
self.reg.eps_q = np.percentile(np.abs(self.reg.regmesh.cellDiffxStencil*(self.reg.mapping * self.invProb.curModel)),self.prctile)
else:
self.reg.eps_q = self.eps[1]
print "L[p qx qy qz]-norm : " + str(self.reg.norms)
print "eps_p: " + str(self.reg.eps_p) + " eps_q: " + str(self.reg.eps_q)
print self.eps_p, self.eps_q, self.norms
self.reg.eps_p = self.eps_p
self.reg.eps_q = self.eps_q
self.reg.norms = self.norms
self.coolingFactor = 1.
self.coolingRate = 1
@@ -356,14 +334,14 @@ class Update_IRLS(InversionDirective):
else:
self.f_old = phim_new
# # Cool the threshold parameter if required
# if getattr(self, 'factor', None) is not None:
# eps = self.reg.eps / self.factor
#
# if getattr(self, 'eps_min', None) is not None:
# self.reg.eps = np.max([self.eps_min,eps])
# else:
# self.reg.eps = eps
# Cool the threshold parameter if required
if getattr(self, 'factor', None) is not None:
eps = self.reg.eps / self.factor
if getattr(self, 'eps_min', None) is not None:
self.reg.eps = np.max([self.eps_min,eps])
else:
self.reg.eps = eps
# Get phi_m at the end of current iteration
self.phi_m_last = self.invProb.phi_m_last
-302
View File
@@ -1,302 +0,0 @@
from __future__ import division
import numpy as np
from scipy.constants import mu_0, pi, epsilon_0
from scipy.special import erf
from SimPEG import Utils
omega = lambda f: 2.*np.pi*f
# TODO:
# r = lambda dx, dy, dz: np.sqrt( dx**2. + dy**2. + dz**2.)
# k = lambda f, mu, epsilon, sig: np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig )
def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=0., epsr=1.):
"""
Computing Analytic Electric fields from Electrical Dipole in a Wholespace
TODO:
Add description of parameters
"""
mu = mu_0*(1+kappa)
epsilon = epsilon_0*epsr
sig_hat = sig + 1j*omega(f)*epsilon
XYZ = Utils.asArray_N_x_Dim(XYZ, 3)
# Check
if XYZ.shape[0] > 1 & f.shape[0] > 1:
raise Exception("I/O type error: For multiple field locations only a single frequency can be specified.")
dx = XYZ[:,0]-srcLoc[0]
dy = XYZ[:,1]-srcLoc[1]
dz = XYZ[:,2]-srcLoc[2]
r = np.sqrt( dx**2. + dy**2. + dz**2.)
# k = np.sqrt( -1j*2.*np.pi*f*mu*sig )
k = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig )
front = current * length / (4.*np.pi*sig_hat* r**3) * np.exp(-1j*k*r)
mid = -k**2 * r**2 + 3*1j*k*r + 3
if orientation.upper() == 'X':
Ex = front*((dx**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.))
Ey = front*(dx*dy / r**2)*mid
Ez = front*(dx*dz / r**2)*mid
return Ex, Ey, Ez
elif orientation.upper() == 'Y':
# x--> y, y--> z, z-->x
Ey = front*((dy**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.))
Ez = front*(dy*dz / r**2)*mid
Ex = front*(dy*dx / r**2)*mid
return Ex, Ey, Ez
elif orientation.upper() == 'Z':
# x --> z, y --> x, z --> y
Ez = front*((dz**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.))
Ex = front*(dz*dx / r**2)*mid
Ey = front*(dz*dy / r**2)*mid
return Ex, Ey, Ez
def E_galvanic_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.):
"""
Computing Galvanic portion of Electric fields from Electrical Dipole in a Wholespace
TODO:
Add description of parameters
"""
mu = mu_0*(1+kappa)
epsilon = epsilon_0*epsr
sig_hat = sig + 1j*omega(f)*epsilon
XYZ = Utils.asArray_N_x_Dim(XYZ, 3)
# Check
if XYZ.shape[0] > 1 & f.shape[0] > 1:
raise Exception("I/O type error: For multiple field locations only a single frequency can be specified.")
dx = XYZ[:,0]-srcLoc[0]
dy = XYZ[:,1]-srcLoc[1]
dz = XYZ[:,2]-srcLoc[2]
r = np.sqrt( dx**2. + dy**2. + dz**2.)
# k = np.sqrt( -1j*2.*np.pi*f*mu*sig )
k = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig )
front = current * length / (4.*np.pi*sig_hat* r**3) * np.exp(-1j*k*r)
mid = -k**2 * r**2 + 3*1j*k*r + 3
if orientation.upper() == 'X':
Ex_galvanic = front*((dx**2 / r**2)*mid + (-1j*k*r-1.))
Ey_galvanic = front*(dx*dy / r**2)*mid
Ez_galvanic = front*(dx*dz / r**2)*mid
return Ex_galvanic, Ey_galvanic, Ez_galvanic
elif orientation.upper() == 'Y':
# x--> y, y--> z, z-->x
Ey_galvanic = front*((dy**2 / r**2)*mid + (-1j*k*r-1.))
Ez_galvanic = front*(dy*dz / r**2)*mid
Ex_galvanic = front*(dy*dx / r**2)*mid
return Ex_galvanic, Ey_galvanic, Ez_galvanic
elif orientation.upper() == 'Z':
# x --> z, y --> x, z --> y
Ez_galvanic = front*((dz**2 / r**2)*mid + (-1j*k*r-1.))
Ex_galvanic = front*(dz*dx / r**2)*mid
Ey_galvanic = front*(dz*dy / r**2)*mid
return Ex_galvanic, Ey_galvanic, Ez_galvanic
def E_inductive_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.):
"""
Computing Inductive portion of Electric fields from Electrical Dipole in a Wholespace
TODO:
Add description of parameters
"""
mu = mu_0*(1+kappa)
epsilon = epsilon_0*epsr
sig_hat = sig + 1j*omega(f)*epsilon
XYZ = Utils.asArray_N_x_Dim(XYZ, 3)
# Check
if XYZ.shape[0] > 1 & f.shape[0] > 1:
raise Exception("I/O type error: For multiple field locations only a single frequency can be specified.")
dx = XYZ[:,0]-srcLoc[0]
dy = XYZ[:,1]-srcLoc[1]
dz = XYZ[:,2]-srcLoc[2]
r = np.sqrt( dx**2. + dy**2. + dz**2.)
# k = np.sqrt( -1j*2.*np.pi*f*mu*sig )
k = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig )
front = current * length / (4.*np.pi*sig_hat* r**3) * np.exp(-1j*k*r)
if orientation.upper() == 'X':
Ex_inductive = front*(k**2 * r**2)
Ey_inductive = np.zeros_like(Ex_inductive)
Ez_inductive = np.zeros_like(Ex_inductive)
return Ex_inductive, Ey_inductive, Ez_inductive
elif orientation.upper() == 'Y':
# x--> y, y--> z, z-->x
Ey_inductive = front*(k**2 * r**2)
Ez_inductive = np.zeros_like(Ey_inductive)
Ex_inductive = np.zeros_like(Ey_inductive)
return Ex_inductive, Ey_inductive, Ez_inductive
elif orientation.upper() == 'Z':
# x --> z, y --> x, z --> y
Ez_inductive = front*(k**2 * r**2)
Ex_inductive = np.zeros_like(Ez_inductive)
Ey_inductive = np.zeros_like(Ez_inductive)
return Ex_inductive, Ey_inductive, Ez_inductive
def J_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.):
"""
Computing Current densities from Electrical Dipole in a Wholespace
TODO:
Add description of parameters
"""
Ex, Ey, Ez = E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=current, length=length, orientation=orientation, kappa=kappa, epsr=epsr)
Jx = sig*Ex
Jy = sig*Ey
Jz = sig*Ez
return Jx, Jy, Jz
def J_galvanic_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.):
"""
Computing Galvanic portion of Current densities from Electrical Dipole in a Wholespace
TODO:
Add description of parameters
"""
Ex_galvanic, Ey_galvanic, Ez_galvanic = E_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=current, length=length, orientation=orientation, kappa=kappa, epsr=epsr)
Jx_galvanic = sig*Ex_galvanic
Jy_galvanic = sig*Ey_galvanic
Jz_galvanic = sig*Ez_galvanic
return Jx_galvanic, Jy_galvanic, Jz_galvanic
def J_inductive_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.):
"""
Computing Inductive portion of Current densities from Electrical Dipole in a Wholespace
TODO:
Add description of parameters
"""
Ex_inductive, Ey_inductive, Ez_inductive = E_inductive_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=current, length=length, orientation=orientation, kappa=kappa, epsr=epsr)
Jx_inductive = sig*Ex_inductive
Jy_inductive = sig*Ey_inductive
Jz_inductive = sig*Ez_inductive
return Jx_inductive, Jy_inductive, Jz_inductive
def H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.):
"""
Computing Magnetic fields from Electrical Dipole in a Wholespace
TODO:
Add description of parameters
"""
mu = mu_0*(1+kappa)
epsilon = epsilon_0*epsr
XYZ = Utils.asArray_N_x_Dim(XYZ, 3)
# Check
if XYZ.shape[0] > 1 & f.shape[0] > 1:
raise Exception("I/O type error: For multiple field locations only a single frequency can be specified.")
dx = XYZ[:,0]-srcLoc[0]
dy = XYZ[:,1]-srcLoc[1]
dz = XYZ[:,2]-srcLoc[2]
r = np.sqrt( dx**2. + dy**2. + dz**2.)
# k = np.sqrt( -1j*2.*np.pi*f*mu*sig )
k = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig )
front = current * length / (4.*np.pi* r**2) * (-1j*k*r + 1) * np.exp(-1j*k*r)
if orientation.upper() == 'X':
Hy = front*(-dz / r)
Hz = front*(dy / r)
Hx = np.zeros_like(Hy)
return Hx, Hy, Hz
elif orientation.upper() == 'Y':
Hx = front*(dz / r)
Hz = front*(-dx / r)
Hy = np.zeros_like(Hx)
return Hx, Hy, Hz
elif orientation.upper() == 'Z':
Hx = front*(-dy / r)
Hy = front*(dx / r)
Hz = np.zeros_like(Hx)
return Hx, Hy, Hz
def B_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.):
"""
Computing Magnetic flux densites from Electrical Dipole in a Wholespace
TODO:
Add description of parameters
"""
Hx, Hy, Hz = H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=current, length=length, orientation=orientation, kappa=kappa, epsr=epsr)
Bx = mu*Hx
By = mu*Hy
Bz = mu*Hz
return Bx, By, Bz
def A_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.):
"""
Computing Electric vector potentials from Electrical Dipole in a Wholespace
TODO:
Add description of parameters
"""
mu = mu_0*(1+kappa)
epsilon = epsilon_0*epsr
XYZ = Utils.asArray_N_x_Dim(XYZ, 3)
# Check
if XYZ.shape[0] > 1 & f.shape[0] > 1:
raise Exception("I/O type error: For multiple field locations only a single frequency can be specified.")
dx = XYZ[:,0]-srcLoc[0]
dy = XYZ[:,1]-srcLoc[1]
dz = XYZ[:,2]-srcLoc[2]
r = np.sqrt( dx**2. + dy**2. + dz**2.)
k = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig )
front = current * length / (4.*np.pi*r)
if orientation.upper() == 'X':
Ax = front*np.exp(-1j*k*r)
Ay = np.zeros_like(Ax)
Az = np.zeros_like(Ax)
return Ax, Ay, Az
elif orientation.upper() == 'Y':
Ay = front*np.exp(-1j*k*r)
Ax = np.zeros_like(Ay)
Az = np.zeros_like(Ay)
return Ax, Ay, Az
elif orientation.upper() == 'Z':
Az = front*np.exp(-1j*k*r)
Ax = np.zeros_like(Ay)
Ay = np.zeros_like(Ay)
return Ax, Ay, Az
-1
View File
@@ -2,4 +2,3 @@ from TDEM import hzAnalyticDipoleT
from FDEM import hzAnalyticDipoleF
from FDEMcasing import *
from DC import DCAnalyticHalf, DCAnalyticSphere
from FDEMDipolarfields import *
+3 -4
View File
@@ -20,10 +20,10 @@ class BaseEMProblem(Problem.BaseProblem):
Problem.BaseProblem.__init__(self, mesh, **kwargs)
surveyPair = Survey.BaseSurvey #: The survey to pair with.
dataPair = Survey.Data #: The data to pair with.
surveyPair = Survey.BaseSurvey
dataPair = Survey.Data
PropMap = EMPropMap #: The property mapping
PropMap = EMPropMap
Solver = SimpegSolver
solverOpts = {}
@@ -217,7 +217,6 @@ class BaseEMSurvey(Survey.BaseSurvey):
def eval(self, f):
"""
Project fields to receiver locations
:param Fields u: fields object
:rtype: numpy.ndarray
:return: data
+30 -18
View File
@@ -6,11 +6,11 @@ from SimPEG.EM.Utils import omega
from SimPEG.Utils import Zero, Identity, sdiag
class FieldsFDEM(SimPEG.Problem.Fields):
class Fields(SimPEG.Problem.Fields):
"""
Fancy Field Storage for a FDEM survey. Only one field type is stored for
each problem, the rest are computed. The fields object acts like an array and is indexed by
each problem, the rest are computed. The fields obejct acts like an array and is indexed by
.. code-block:: python
@@ -92,7 +92,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
"""
Total derivative of e with respect to the inversion model. Returns :math:`d\mathbf{e}/d\mathbf{m}` for forward and (:math:`d\mathbf{e}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint
:param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: source
:param Src src: sorce
:param numpy.ndarray du_dm_v: derivative of the solution vector with respect to the model times a vector (is None for adjoint)
:param numpy.ndarray v: vector to take sensitivity product with
:param bool adjoint: adjoint?
@@ -110,7 +110,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
"""
Total derivative of b with respect to the inversion model. Returns :math:`d\mathbf{b}/d\mathbf{m}` for forward and (:math:`d\mathbf{b}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint
:param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: source
:param Src src: sorce
:param numpy.ndarray du_dm_v: derivative of the solution vector with respect to the model times a vector (is None for adjoint)
:param numpy.ndarray v: vector to take sensitivity product with
:param bool adjoint: adjoint?
@@ -128,7 +128,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
"""
Total derivative of h with respect to the inversion model. Returns :math:`d\mathbf{h}/d\mathbf{m}` for forward and (:math:`d\mathbf{h}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint
:param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: source
:param Src src: sorce
:param numpy.ndarray du_dm_v: derivative of the solution vector with respect to the model times a vector (is None for adjoint)
:param numpy.ndarray v: vector to take sensitivity product with
:param bool adjoint: adjoint?
@@ -146,7 +146,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
"""
Total derivative of j with respect to the inversion model. Returns :math:`d\mathbf{j}/d\mathbf{m}` for forward and (:math:`d\mathbf{j}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint
:param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: source
:param Src src: sorce
:param numpy.ndarray du_dm_v: derivative of the solution vector with respect to the model times a vector (is None for adjoint)
:param numpy.ndarray v: vector to take sensitivity product with
:param bool adjoint: adjoint?
@@ -160,12 +160,12 @@ class FieldsFDEM(SimPEG.Problem.Fields):
return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint)
return np.array(self._jDeriv_u(src, du_dm_v, adjoint) + self._jDeriv_m(src, v, adjoint), dtype = complex)
class Fields3D_e(FieldsFDEM):
class Fields3D_e(Fields):
"""
Fields object for Problem3D_e.
:param BaseMesh mesh: mesh
:param SimPEG.EM.FDEM.SurveyFDEM.Survey survey: survey
:param Mesh mesh: mesh
:param Survey survey: survey
"""
knownFields = {'eSolution':'E'}
@@ -180,6 +180,9 @@ class Fields3D_e(FieldsFDEM):
'h' : ['eSolution','CCV','_h'],
}
def __init__(self, mesh, survey, **kwargs):
Fields.__init__(self, mesh, survey, **kwargs)
def startup(self):
self.prob = self.survey.prob
self._edgeCurl = self.survey.prob.mesh.edgeCurl
@@ -423,12 +426,12 @@ class Fields3D_e(FieldsFDEM):
class Fields3D_b(FieldsFDEM):
class Fields3D_b(Fields):
"""
Fields object for Problem3D_b.
:param BaseMesh mesh: mesh
:param SimPEG.EM.FDEM.SurveyFDEM.Survey survey: survey
:param Mesh mesh: mesh
:param Survey survey: survey
"""
knownFields = {'bSolution':'F'}
@@ -443,6 +446,9 @@ class Fields3D_b(FieldsFDEM):
'h' : ['bSolution','CCV','_h'],
}
def __init__(self,mesh,survey,**kwargs):
Fields.__init__(self,mesh,survey,**kwargs)
def startup(self):
self.prob = self.survey.prob
self._edgeCurl = self.survey.prob.mesh.edgeCurl
@@ -687,12 +693,12 @@ class Fields3D_b(FieldsFDEM):
return Zero()
class Fields3D_j(FieldsFDEM):
class Fields3D_j(Fields):
"""
Fields object for Problem3D_j.
:param BaseMesh mesh: mesh
:param SimPEG.EM.FDEM.SurveyFDEM.Survey survey: survey
:param Mesh mesh: mesh
:param Survey survey: survey
"""
knownFields = {'jSolution':'F'}
@@ -707,6 +713,9 @@ class Fields3D_j(FieldsFDEM):
'b' : ['jSolution','CCV','_b'],
}
def __init__(self,mesh,survey,**kwargs):
Fields.__init__(self,mesh,survey,**kwargs)
def startup(self):
self.prob = self.survey.prob
self._edgeCurl = self.survey.prob.mesh.edgeCurl
@@ -979,12 +988,12 @@ class Fields3D_j(FieldsFDEM):
return 1./(1j * omega(src.freq)) * VI * (self._aveE2CCV * ( s_mDeriv(v) - self._edgeCurl.T * ( self._MfRhoDeriv(jSolution) * v ) ) )
class Fields3D_h(FieldsFDEM):
class Fields3D_h(Fields):
"""
Fields object for Problem3D_h.
:param BaseMesh mesh: mesh
:param SimPEG.EM.FDEM.SurveyFDEM.Survey survey: survey
:param Mesh mesh: mesh
:param Survey survey: survey
"""
knownFields = {'hSolution':'E'}
@@ -999,6 +1008,9 @@ class Fields3D_h(FieldsFDEM):
'b' : ['hSolution','CCV','_b'],
}
def __init__(self,mesh,survey,**kwargs):
Fields.__init__(self,mesh,survey,**kwargs)
def startup(self):
self.prob = self.survey.prob
self._edgeCurl = self.survey.prob.mesh.edgeCurl
+16 -21
View File
@@ -1,7 +1,7 @@
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 FieldsFDEM import Fields, Fields3D_e, Fields3D_b, Fields3D_h, Fields3D_j
from SimPEG.EM.Base import BaseEMProblem
from SimPEG.EM.Utils import omega
@@ -31,11 +31,10 @@ class BaseFDEMProblem(BaseEMProblem):
if using the H-J formulation (:code:`Problem3D_j` or :code:`Problem3D_h`). Note that here, :math:`\mathbf{s_m}` is an integrated quantity.
The problem performs the elimination so that we are solving the system for \\\(\\\mathbf{e},\\\mathbf{b},\\\mathbf{j} \\\) or \\\(\\\mathbf{h}\\\)
"""
surveyPair = SurveyFDEM
fieldsPair = FieldsFDEM
fieldsPair = Fields
def fields(self, m):
"""
@@ -65,7 +64,7 @@ class BaseFDEMProblem(BaseEMProblem):
:param numpy.array m: inversion model (nP,)
:param numpy.array v: vector which we take sensitivity product with (nP,)
:param SimPEG.EM.FDEM.FieldsFDEM.FieldsFDEM u: fields object
:param SimPEG.EM.FDEM.Fields u: fields object
:rtype numpy.array:
:return: Jv (ndata,)
"""
@@ -100,7 +99,7 @@ class BaseFDEMProblem(BaseEMProblem):
:param numpy.array m: inversion model (nP,)
:param numpy.array v: vector which we take adjoint product with (nP,)
:param SimPEG.EM.FDEM.FieldsFDEM.FieldsFDEM u: fields object
:param SimPEG.EM.FDEM.Fields u: fields object
:rtype numpy.array:
:return: Jv (ndata,)
"""
@@ -154,8 +153,8 @@ class BaseFDEMProblem(BaseEMProblem):
Evaluates the sources for a given frequency and puts them in matrix form
:param float freq: Frequency
:rtype: tuple
:return: (s_m, s_e) (nE or nF, nSrc)
:rtype: (numpy.ndarray, numpy.ndarray)
:return: s_m, s_e (nE or nF, nSrc)
"""
Srcs = self.survey.getSrcByFreq(freq)
if self._formulation is 'EB':
@@ -195,7 +194,7 @@ class Problem3D_e(BaseFDEMProblem):
which we solve for :math:`\mathbf{e}`.
:param SimPEG.Mesh.BaseMesh.BaseMesh mesh: mesh
:param SimPEG.Mesh mesh: mesh
"""
_solutionType = 'eSolution'
@@ -270,7 +269,7 @@ class Problem3D_e(BaseFDEMProblem):
Derivative of the right hand side with respect to the model
:param float freq: frequency
:param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
:param SimPEG.EM.FDEM.Src src: FDEM source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
@@ -306,7 +305,7 @@ class Problem3D_b(BaseFDEMProblem):
.. note ::
The inverse problem will not work with full anisotropy
:param SimPEG.Mesh.BaseMesh.BaseMesh mesh: mesh
:param SimPEG.Mesh mesh: mesh
"""
_solutionType = 'bSolution'
@@ -401,7 +400,7 @@ class Problem3D_b(BaseFDEMProblem):
Derivative of the right hand side with respect to the model
:param float freq: frequency
:param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
:param SimPEG.EM.FDEM.Src src: FDEM source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
@@ -445,7 +444,6 @@ class Problem3D_j(BaseFDEMProblem):
\mathbf{h} = \\frac{1}{i \omega} \mathbf{M_{\mu}^e}^{-1} \\left(-\mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{j} + \mathbf{M^e} \mathbf{s_m} \\right)
and solve for \\\(\\\mathbf{j}\\\) using
.. math ::
@@ -455,7 +453,7 @@ class Problem3D_j(BaseFDEMProblem):
.. note::
This implementation does not yet work with full anisotropy!!
:param SimPEG.Mesh.BaseMesh.BaseMesh mesh: mesh
:param SimPEG.Mesh mesh: mesh
"""
_solutionType = 'jSolution'
@@ -531,8 +529,8 @@ class Problem3D_j(BaseFDEMProblem):
\mathbf{RHS} = \mathbf{C} \mathbf{M_{\mu}^e}^{-1}\mathbf{s_m} -i\omega \mathbf{s_e}
:param float freq: Frequency
:rtype: numpy.ndarray
:return: RHS (nE, nSrc)
:rtype: numpy.ndarray (nE, nSrc)
:return: RHS
"""
s_m, s_e = self.getSourceTerm(freq)
@@ -551,7 +549,7 @@ class Problem3D_j(BaseFDEMProblem):
Derivative of the right hand side with respect to the model
:param float freq: frequency
:param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
:param SimPEG.EM.FDEM.Src src: FDEM source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
@@ -593,7 +591,7 @@ class Problem3D_h(BaseFDEMProblem):
\\left(\mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{C} + i \omega \mathbf{M_{\mu}^e}\\right) \mathbf{h} = \mathbf{M^e} \mathbf{s_m} + \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{s_e}
:param SimPEG.Mesh.BaseMesh.BaseMesh mesh: mesh
:param SimPEG.Mesh mesh: mesh
"""
_solutionType = 'hSolution'
@@ -610,11 +608,9 @@ class Problem3D_h(BaseFDEMProblem):
.. math::
\mathbf{A} = \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{C} + i \omega \mathbf{M_{\mu}^e}
:param float freq: Frequency
:rtype: scipy.sparse.csr_matrix
:return: A
"""
MeMu = self.MeMu
@@ -657,7 +653,6 @@ class Problem3D_h(BaseFDEMProblem):
:param float freq: Frequency
:rtype: numpy.ndarray
:return: RHS (nE, nSrc)
"""
s_m, s_e = self.getSourceTerm(freq)
@@ -671,7 +666,7 @@ class Problem3D_h(BaseFDEMProblem):
Derivative of the right hand side with respect to the model
:param float freq: frequency
:param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
:param SimPEG.EM.FDEM.Src src: FDEM source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
+5 -5
View File
@@ -25,10 +25,10 @@ class BaseRx(SimPEG.Survey.BaseRx):
def eval(self, src, mesh, f):
"""
Project fields to receivers to get data.
Project fields to recievers to get data.
:param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
:param BaseMesh mesh: mesh used
:param Source src: FDEM source
:param Mesh mesh: mesh used
:param Fields f: fields object
:rtype: numpy.ndarray
:return: fields projected to recievers
@@ -44,8 +44,8 @@ class BaseRx(SimPEG.Survey.BaseRx):
"""
Derivative of projected fields with respect to the inversion model times a vector.
:param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
:param BaseMesh mesh: mesh used
:param Source src: FDEM source
:param Mesh mesh: mesh used
:param Fields f: fields object
:param numpy.ndarray v: vector to multiply
:rtype: numpy.ndarray
+28 -28
View File
@@ -23,8 +23,8 @@ class BaseSrc(Survey.BaseSrc):
- :math:`s_m` : magnetic source term
- :math:`s_e` : electric source term
:param BaseFDEMProblem prob: FDEM Problem
:rtype: tuple
:param Problem prob: FDEM Problem
:rtype: (numpy.ndarray, numpy.ndarray)
:return: tuple with magnetic source term and electric source term
"""
s_m = self.s_m(prob)
@@ -37,10 +37,10 @@ class BaseSrc(Survey.BaseSrc):
- :code:`s_mDeriv` : derivative of the magnetic source term
- :code:`s_eDeriv` : derivative of the electric source term
:param BaseFDEMProblem prob: FDEM Problem
:param Problem prob: FDEM Problem
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: tuple
:rtype: (numpy.ndarray, numpy.ndarray)
:return: tuple with magnetic source term and electric source term derivatives times a vector
"""
if v is not None:
@@ -52,7 +52,7 @@ class BaseSrc(Survey.BaseSrc):
"""
Primary magnetic flux density
:param BaseFDEMProblem prob: FDEM Problem
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
:return: primary magnetic flux density
"""
@@ -64,7 +64,7 @@ class BaseSrc(Survey.BaseSrc):
"""
Primary magnetic field
:param BaseFDEMProblem prob: FDEM Problem
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
@@ -76,7 +76,7 @@ class BaseSrc(Survey.BaseSrc):
"""
Primary electric field
:param BaseFDEMProblem prob: FDEM Problem
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
:return: primary electric field
"""
@@ -88,7 +88,7 @@ class BaseSrc(Survey.BaseSrc):
"""
Primary current density
:param BaseFDEMProblem prob: FDEM Problem
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
:return: primary current density
"""
@@ -100,7 +100,7 @@ class BaseSrc(Survey.BaseSrc):
"""
Magnetic source term
:param BaseFDEMProblem prob: FDEM Problem
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
:return: magnetic source term on mesh
"""
@@ -110,7 +110,7 @@ class BaseSrc(Survey.BaseSrc):
"""
Electric source term
:param BaseFDEMProblem prob: FDEM Problem
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
:return: electric source term on mesh
"""
@@ -120,7 +120,7 @@ class BaseSrc(Survey.BaseSrc):
"""
Derivative of magnetic source term with respect to the inversion model
:param BaseFDEMProblem prob: FDEM Problem
:param Problem prob: FDEM Problem
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
@@ -133,7 +133,7 @@ class BaseSrc(Survey.BaseSrc):
"""
Derivative of electric source term with respect to the inversion model
:param BaseFDEMProblem prob: FDEM Problem
:param Problem prob: FDEM Problem
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
@@ -162,7 +162,7 @@ class RawVec_e(BaseSrc):
"""
Electric source term
:param BaseFDEMProblem prob: FDEM Problem
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
:return: electric source term on mesh
"""
@@ -191,7 +191,7 @@ class RawVec_m(BaseSrc):
"""
Magnetic source term
:param BaseFDEMProblem prob: FDEM Problem
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
:return: magnetic source term on mesh
"""
@@ -220,7 +220,7 @@ class RawVec(BaseSrc):
"""
Magnetic source term
:param BaseFDEMProblem prob: FDEM Problem
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
:return: magnetic source term on mesh
"""
@@ -232,7 +232,7 @@ class RawVec(BaseSrc):
"""
Electric source term
:param BaseFDEMProblem prob: FDEM Problem
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
:return: electric source term on mesh
"""
@@ -301,7 +301,7 @@ class MagDipole(BaseSrc):
"""
The primary magnetic flux density from a magnetic vector potential
:param BaseFDEMProblem prob: FDEM problem
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
@@ -339,7 +339,7 @@ class MagDipole(BaseSrc):
"""
The primary magnetic field from a magnetic vector potential
:param BaseFDEMProblem prob: FDEM problem
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
@@ -350,7 +350,7 @@ class MagDipole(BaseSrc):
"""
The magnetic source term
:param BaseFDEMProblem prob: FDEM problem
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
@@ -364,7 +364,7 @@ class MagDipole(BaseSrc):
"""
The electric source term
:param BaseFDEMProblem prob: FDEM problem
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
@@ -416,7 +416,7 @@ class MagDipole_Bfield(BaseSrc):
"""
The primary magnetic flux density from the analytic solution for magnetic fields from a dipole
:param BaseFDEMProblem prob: FDEM problem
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
@@ -455,7 +455,7 @@ class MagDipole_Bfield(BaseSrc):
"""
The primary magnetic field from a magnetic vector potential
:param BaseFDEMProblem prob: FDEM problem
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
@@ -466,7 +466,7 @@ class MagDipole_Bfield(BaseSrc):
"""
The magnetic source term
:param BaseFDEMProblem prob: FDEM problem
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
@@ -479,7 +479,7 @@ class MagDipole_Bfield(BaseSrc):
"""
The electric source term
:param BaseFDEMProblem prob: FDEM problem
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
@@ -530,7 +530,7 @@ class CircularLoop(BaseSrc):
"""
The primary magnetic flux density from a magnetic vector potential
:param BaseFDEMProblem prob: FDEM problem
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
@@ -567,7 +567,7 @@ class CircularLoop(BaseSrc):
"""
The primary magnetic field from a magnetic vector potential
:param BaseFDEMProblem prob: FDEM problem
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
@@ -578,7 +578,7 @@ class CircularLoop(BaseSrc):
"""
The magnetic source term
:param BaseFDEMProblem prob: FDEM problem
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
@@ -591,7 +591,7 @@ class CircularLoop(BaseSrc):
"""
The electric source term
:param BaseFDEMProblem prob: FDEM problem
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
+3 -3
View File
@@ -112,7 +112,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
"""
:param numpy.array m: Conductivity model
:param numpy.ndarray v: vector (model object)
:param FieldsTDEM f: Fields resulting from m
:param simpegEM.TDEM.FieldsTDEM f: Fields resulting from m
:rtype: numpy.ndarray
:return: w (data object)
@@ -136,8 +136,8 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
def Jtvec(self, m, v, f=None):
"""
:param numpy.array m: Conductivity model
:param numpy.ndarray v: vector (or a :class:`SimPEG.Survey.Data` object)
:param FieldsTDEM u: Fields resulting from m
:param numpy.ndarray,SimPEG.Survey.Data v: vector (data object)
:param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m
:rtype: numpy.ndarray
:return: w (model object)
+13 -13
View File
@@ -87,8 +87,8 @@ class ProblemTDEM_b(BaseTDEMProblem):
"""
:param numpy.array m: Conductivity model
:param numpy.array vec: vector (like a model)
:param FieldsTDEM u: Fields resulting from m
:rtype: FieldsTDEM
:param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m
:rtype: simpegEM.TDEM.FieldsTDEM
:return: f
Multiply G by a vector
@@ -125,9 +125,9 @@ class ProblemTDEM_b(BaseTDEMProblem):
"""
:param numpy.array m: Conductivity model
:param numpy.array vec: vector (like a fields)
:param FieldsTDEM u: Fields resulting from m
:rtype: numpy.ndarray
:return: p (like a model)
:param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m
:rtype: np.ndarray (like a model)
:return: p
Multiply G.T by a vector
"""
@@ -153,8 +153,8 @@ class ProblemTDEM_b(BaseTDEMProblem):
def solveAh(self, m, p):
"""
:param numpy.array m: Conductivity model
:param FieldsTDEM p: Fields object
:rtype: FieldsTDEM
:param simpegEM.TDEM.FieldsTDEM p: Fields object
:rtype: simpegEM.TDEM.FieldsTDEM
:return: y
Solve the block-matrix system \\\(\\\hat{A} \\\hat{y} = \\\hat{p}\\\):
@@ -200,8 +200,8 @@ class ProblemTDEM_b(BaseTDEMProblem):
def solveAht(self, m, p):
"""
:param numpy.array m: Conductivity model
:param FieldsTDEM p: Fields object
:rtype: FieldsTDEM
:param simpegEM.TDEM.FieldsTDEM p: Fields object
:rtype: simpegEM.TDEM.FieldsTDEM
:return: y
Solve the block-matrix system \\\(\\\hat{A}^\\\\top \\\hat{y} = \\\hat{p}\\\):
@@ -270,8 +270,8 @@ class ProblemTDEM_b(BaseTDEMProblem):
def _AhVec(self, m, vec):
"""
:param numpy.array m: Conductivity model
:param FieldsTDEM vec: Fields object
:rtype: FieldsTDEM
:param simpegEM.TDEM.FieldsTDEM vec: Fields object
:rtype: simpegEM.TDEM.FieldsTDEM
:return: f
Multiply the matrix \\\(\\\hat{A}\\\) by a fields vector where
@@ -315,8 +315,8 @@ class ProblemTDEM_b(BaseTDEMProblem):
def _AhtVec(self, m, vec):
"""
:param numpy.array m: Conductivity model
:param FieldsTDEM vec: Fields object
:rtype: FieldsTDEM
:param simpegEM.TDEM.FieldsTDEM vec: Fields object
:rtype: simpegEM.TDEM.FieldsTDEM
:return: f
Multiply the matrix \\\(\\\hat{A}\\\) by a fields vector where
+7 -7
View File
@@ -1,7 +1,7 @@
from SimPEG import *
import SimPEG.EM.Static.DC as DC
import SimPEG.DCIP as DC
def run(plotIt=True):
def run(plotIt=False):
cs = 25.
hx = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
hy = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
@@ -21,10 +21,10 @@ def run(plotIt=True):
# ax.plot(xyz_rxP[:,0],xyz_rxP[:,1], 'w.')
# ax.plot(xyz_rxN[:,0],xyz_rxN[:,1], 'r.', ms = 3)
rx = DC.Rx.Dipole(xyz_rxP, xyz_rxN)
src = DC.Src.Dipole([rx], np.r_[-200, 0, -12.5], np.r_[+200, 0, -12.5])
survey = DC.Survey([src])
problem = DC.Problem3D_CC(mesh)
rx = DC.RxDipole(xyz_rxP, xyz_rxN)
src = DC.SrcDipole([rx], [-200, 0, -12.5], [+200, 0, -12.5])
survey = DC.SurveyDC([src])
problem = DC.ProblemDC_CC(mesh)
problem.pair(survey)
try:
from pymatsolver import MumpsSolver
@@ -65,4 +65,4 @@ def run(plotIt=True):
if __name__ == '__main__':
print run()
print run(plotIt=True)
@@ -19,13 +19,10 @@ def run(plotIt=True):
Morrison Casing Model, and the results are used in a 2016 SEG abstract by
Yang et al.
.. code-block:: text
Schenkel, C.J., and H.F. Morrison, 1990, Effects of well casing on potential field measurements using downhole current sources: Geophysical prospecting, 38, 663-686.
- Schenkel, C.J., and H.F. Morrison, 1990, Effects of well casing on potential field measurements using downhole current sources: Geophysical prospecting, 38, 663-686.
The model consists of:
- Air: Conductivity 1e-8 S/m, above z = 0
- Background: conductivity 1e-2 S/m, below z = 0
- Casing: conductivity 1e6 S/m
@@ -218,7 +215,7 @@ def run(plotIt=True):
# ------------ Problem and Survey ---------------
survey = FDEM.Survey(sg_p + dg_p)
mapping = [('sigma', Maps.IdentityMap(mesh))]
problem = FDEM.Problem3D_h(mesh, mapping=mapping, Solver=solver)
problem = FDEM.Problem3D_h(mesh, mapping=mapping)
problem.pair(survey)
# ------------- Solve ---------------------------
+29 -7
View File
@@ -42,33 +42,55 @@ def run(N=100, plotIt=True):
survey = Survey.LinearSurvey()
survey.pair(prob)
survey.dobs = prob.fields(mtrue) + std_noise * np.random.randn(nk)
#survey.makeSyntheticData(mtrue, std=std_noise)
wd = np.ones(nk) * std_noise
#print survey.std[0]
#M = prob.mesh
# Distance weighting
wr = np.sum(prob.G**2.,axis=0)**0.5
wr = ( wr/np.max(wr) )
# reg = Regularization.Simple(mesh)
# reg.mref = mref
# reg.cell_weights = wr
#
dmis = DataMisfit.l2_DataMisfit(survey)
dmis.Wd = 1./wd
#
# opt = Optimization.ProjectedGNCG(maxIter=20,lower=-2.,upper=2., maxIterCG= 10, tolCG = 1e-4)
# invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
# invProb.curModel = m0
#
# beta = Directives.BetaSchedule(coolingFactor=2, coolingRate=1)
# target = Directives.TargetMisfit()
#
betaest = Directives.BetaEstimate_ByEig()
# inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest, target])
#
#
# mrec = inv.run(m0)
# ml2 = mrec
# print "Final misfit:" + str(invProb.dmisfit.eval(mrec))
#
# # Switch regularization to sparse
# phim = invProb.phi_m_last
# phid = invProb.phi_d
reg = Regularization.Sparse(mesh)
reg.mref = mref
reg.cell_weights = wr
reg.mref = np.zeros(mesh.nC)
eps_p = 5e-2
eps_q = 5e-2
norms = [0., 0., 2., 2.]
opt = Optimization.ProjectedGNCG(maxIter=100 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 10, tolCG = 1e-3)
invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
update_Jacobi = Directives.Update_lin_PreCond()
# Set the IRLS directive, penalize the lowest 25 percentile of model values
# Start with an l2-l2, then switch to lp-norms
norms = [0., 0., 2., 2.]
IRLS = Directives.Update_IRLS( norms=norms, prctile = 25, maxIRLSiter = 15, minGNiter=3)
IRLS = Directives.Update_IRLS( norms=norms, eps_p=eps_p, eps_q=eps_q)
inv = Inversion.BaseInversion(invProb, directiveList=[IRLS,betaest,update_Jacobi])
+30 -24
View File
@@ -1,13 +1,15 @@
import SimPEG as simpeg
import numpy as np
import SimPEG.MT as MT
from SimPEG import NSEM
from scipy.constants import mu_0
import matplotlib.pyplot as plt
np.random.seed(1983)
def run(plotIt=True):
"""
MT: 1D: Inversion
=================
=======================
Forward model 1D MT data.
Setup and run a MT 1D inversion.
@@ -17,13 +19,13 @@ def run(plotIt=True):
## Setup the forward modeling
# Setting up 1D mesh and conductivity models to forward model data.
# Frequency
nFreq = 31
freqs = np.logspace(3,-3,nFreq)
nFreq = 26
freqs = np.logspace(2,-3,nFreq)
# Set mesh parameters
ct = 20
air = simpeg.Utils.meshTensor([(ct,16,1.4)])
ct = 10
air = simpeg.Utils.meshTensor([(ct,25,1.4)])
core = np.concatenate( ( np.kron(simpeg.Utils.meshTensor([(ct,10,-1.3)]),np.ones((5,))) , simpeg.Utils.meshTensor([(ct,5)]) ) )
bot = simpeg.Utils.meshTensor([(core[0],10,-1.4)])
bot = simpeg.Utils.meshTensor([(core[0],25,-1.4)])
x0 = -np.array([np.sum(np.concatenate((core,bot)))])
# Make the model
m1d = simpeg.Mesh.TensorMesh([np.concatenate((bot,core,air))], x0=x0)
@@ -33,7 +35,7 @@ def run(plotIt=True):
layer1 = (m1d.vectorCCx<-500.) & (m1d.vectorCCx>=-800.)
layer2 = (m1d.vectorCCx<-3500.) & (m1d.vectorCCx>=-5000.)
# Set the conductivity values
sig_half = 2e-3
sig_half = 1e-2
sig_air = 1e-8
sig_layer1 = .2
sig_layer2 = .2
@@ -57,31 +59,31 @@ def run(plotIt=True):
# Receivers
rxList = []
for rxType in ['z1dr','z1di']:
rxList.append(MT.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType))
rxList.append(NSEM.Rx(simpeg.mkvc(np.array([-0.5]),2).T,rxType))
# Source list
srcList =[]
for freq in freqs:
srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq))
# Make the survey
survey = MT.Survey(srcList)
survey = NSEM.Survey(srcList)
survey.mtrue = m_true
## Set the problem
problem = MT.Problem1D.eForm_psField(m1d,sigmaPrimary=sigma_0,mapping=mappingExpAct)
problem = NSEM.Problem1D_ePrimSec(m1d,sigmaPrimary=sigma_0,mapping=mappingExpAct)
problem.pair(survey)
## Forward model data
# Project the data
survey.dtrue = survey.dpred(m_true)
survey.dobs = survey.dtrue + 0.025*abs(survey.dtrue)*np.random.randn(*survey.dtrue.shape)
survey.dobs = survey.dtrue + 0.01*abs(survey.dtrue)*np.random.randn(*survey.dtrue.shape)
if plotIt:
fig = MT.Utils.dataUtils.plotMT1DModelData(problem, [m_0])
fig = NSEM.Utils.dataUtils.plotMT1DModelData(problem,[])
fig.suptitle('Target - smooth true')
# Assign uncertainties
std = 0.05 # 5% std
std = 0.025 # 5% std
survey.std = np.abs(survey.dobs*std)
# Assign the data weight
Wd = 1./survey.std
@@ -90,30 +92,33 @@ def run(plotIt=True):
# Define a counter
C = simpeg.Utils.Counter()
# Set the optimization
opt = simpeg.Optimization.InexactGaussNewton(maxIter = 30)
opt = simpeg.Optimization.ProjectedGNCG(maxIter = 25)
opt.counter = C
opt.LSshorten = 0.5
opt.lower = np.log(1e-4)
opt.upper = np.log(5)
opt.LSshorten = 0.1
opt.remember('xc')
# Data misfit
dmis = simpeg.DataMisfit.l2_DataMisfit(survey)
dmis.Wd = Wd
# Regularization - with a regularization mesh
regMesh = simpeg.Mesh.TensorMesh([m1d.hx[problem.mapping.sigmaMap.maps[-1].indActive]],m1d.x0)
regMesh = simpeg.Mesh.TensorMesh([m1d.hx[active]],m1d.x0)
reg = simpeg.Regularization.Tikhonov(regMesh)
reg.mrefInSmooth = True
reg.alpha_s = 1e-7
reg.alpha_s = 1e-1
reg.alpha_x = 1.
# Inversion problem
invProb = simpeg.InvProblem.BaseInvProblem(dmis, reg, opt)
invProb.counter = C
# Beta cooling
beta = simpeg.Directives.BetaSchedule()
beta.coolingRate = 4
betaest = simpeg.Directives.BetaEstimate_ByEig(beta0_ratio=0.75)
beta.coolingRate = 4.
beta.coolingFactor = 4.
betaest = simpeg.Directives.BetaEstimate_ByEig(beta0_ratio=1.)
betaest.beta0 = 1.
targmis = simpeg.Directives.TargetMisfit()
targmis.target = survey.nD
saveModel = simpeg.Directives.SaveModelEveryIteration()
saveModel.fileName = 'Inversion_TargMisEqnD_smoothTrue'
# Create an inversion object
inv = simpeg.Inversion.BaseInversion(invProb, directiveList=[beta,betaest,targmis])
@@ -121,8 +126,9 @@ def run(plotIt=True):
mopt = inv.run(m_0)
if plotIt:
fig = MT.Utils.dataUtils.plotMT1DModelData(problem,[mopt])
fig = NSEM.Utils.dataUtils.plotMT1DModelData(problem,[mopt])
fig.suptitle('Target - smooth true')
fig.axes[0].set_ylim([-10000,500])
plt.show()
if __name__ == '__main__':
@@ -0,0 +1,428 @@
from scipy.constants import epsilon_0, mu_0
import matplotlib.pyplot as plt
import numpy as np
from ipywidgets import *
from SimPEG.EM.Utils import k, omega
"""
MT1D: n layered earth problem
*****************************
Author: Thibaut Astic
Contact: thast@eos.ubc.ca
Date: January 2016
This code compute the analytic response of a n-layered Earth to a plane wave (Magneto-Tellurics).
We start by looking at Maxwell's equations in the electric
field \\\(\\\mathbf{E}\\) and the magnetic flux
\\\(\\\mathbf{H}\\) to write the wave equations
\\(\\ \nabla ^2 \mathbf{E_x} + k^2 \mathbf{E_x} = 0 \\) &
\\(\\ \nabla ^2 \mathbf{H_y} + k^2 \mathbf{H_y} = 0 \\)
Then solving the equations in each layer "j" between z_{j-1} and z_j in the form of
\\(\\ E_{x,j} (z) = U_j e^{i k (z-z_{j-1})} + D_j e^{-i k (z-z_{j-1})} \\)
\\(\\ H_{y,j} (z) = \frac{1}{Z_j} (D_j e^{-i k (z-z_{j-1})} - U_j e^{i k (z-z_{j-1})}) \\)
With U and D the Up and Down components of the E-field.
The iteration from one layer to another is ensure by:
\\(\\ \left(\begin{matrix} E_{x,j} \\ H_{y,j} \end{matrix} \right) =
P_j T_j P^{-1}_J \left(\begin{matrix} E_{x,j+1} \\ H_{y,j+1} \end{matrix} \right) \\)
And the Boundary Condition is set for the E-field in the last layer, with no Up component (=0)
and only a down component (=1 then normalized by the highest amplitude to ensure numeric stability)
The layer 0 is assumed to be the air layer.
"""
#Define a frquency range for a survey
frange = lambda minfreq, maxfreq, step: np.logspace(minfreq,maxfreq,num = step, base = 10.)
#Functions to create random physical Perties for a n-layered earth
thick = lambda minthick, maxthick, nlayer: np.append(np.array([1.2*10.**5]),
np.ndarray.round(minthick + (maxthick-minthick)* np.random.rand(nlayer-1,1)
,decimals =1))
sig = lambda minsig, maxsig, nlayer: np.append(np.array([0.]),
np.ndarray.round(10.**minsig + (10.**maxsig-10.**minsig)* np.random.rand(nlayer,1)
,decimals=3))
mu = lambda minmu, maxmu, nlayer: np.append(np.array([1.]),
np.ndarray.round(minmu + (maxmu-minmu)* np.random.rand(nlayer,1)
,decimals=1))
eps = lambda mineps, maxeps, nlayer: np.append(np.array([1.]),
np.ndarray.round(mineps + (maxeps-mineps)* np.random.rand(nlayer,1)
,decimals=1))
#Evaluate Impedance Z of a layer
ImpZ = lambda f, mu, k: omega(f)*mu*mu_0/k
#Complex Cole-Cole Conductivity - EM utils
PCC= lambda siginf,m,t,c,f: siginf*(1.-(m/(1.+(1j*omega(f)*t)**c)))
#Converted thickness array into top of layer array
top = lambda thick: np.cumsum(thick)
#Propagation Matrix and theirs inverses
#matrix T for transition of Up and Down components accross a layer
T = lambda h,k: np.matrix([[np.exp(1j*k*h),0.],[0.,np.exp(-1j*k*h)]],dtype='complex_')
Tinv = lambda h,k: np.matrix([[np.exp(-1j*k*h),0.],[0.,np.exp(1j*k*h)]],dtype='complex_')
#transition of Up and Down components accross a layer
UD_Z = lambda UD,z,zj,k : T((z-zj),k)*UD
#matrix P relating Up and Down components with E and H fields
P = lambda z: np.matrix([[1.,1,],[-1./z,1./z]],dtype='complex_')
Pinv = lambda z: np.matrix([[1.,-z],[1.,z]],dtype='complex_')/2.
#Time Variation of E and H
E_ZT = lambda U,D,f,t : np.exp(1j*omega(f)*t)*(U+D)
H_ZT = lambda U,D,Z,f,t : (1./Z)*np.exp(1j*omega(f)*t)*(D-U)
#Plot the configuration of the problem
def PlotConfiguration(thick,sig,eps,mu,ax,widthg,z):
topn = top(thick)
widthn = np.arange(-widthg,widthg+widthg/10.,widthg/10.)
ax.set_ylim([z.min(),z.max()])
ax.set_xlim([-widthg,widthg])
ax.set_ylabel("Depth (m)", fontsize=16.)
ax.yaxis.tick_right()
ax.yaxis.set_label_position("right")
#define filling for the different layers
hatches=['/' , '+', 'x', '|' , '\\', '-' , 'o' , 'O' , '.' , '*' ]
#Write the physical properties of air
ax.annotate(("Air, $\sigma$ =%1.0f mS/m")%(sig[0]*10**(3)),
xy=(-widthg/2., -np.abs(z.max())/2.), xycoords='data',
xytext=(-widthg/2., -np.abs(z.max())/2.), textcoords='data',
fontsize=14.)
ax.annotate(("$\epsilon_r$= %1i")%(eps[0]),
xy=(-widthg/2., -np.abs(z.max())/3.), xycoords='data',
xytext=(-widthg/2., -np.abs(z.max())/3.), textcoords='data',
fontsize=14.)
ax.annotate(("$\mu_r$= %1i")%(mu[0]),
xy=(-widthg/2., -np.abs(z.max())/3.), xycoords='data',
xytext=(0, -np.abs(z.max())/3.), textcoords='data',
fontsize=14.)
#Write the physical properties of the differents layers up to the (n-1)-th and fill it with pattern
for i in range(1,len(topn)-1,1):
if topn[i] == topn[i+1]:
pass
else:
ax.annotate(("$\sigma$ =%3.3f mS/m")%(sig[i]*10**(3)),
xy=(0., (2.*topn[i]+topn[i+1])/3), xycoords='data',
xytext=(0., (2.*topn[i]+topn[i+1])/3), textcoords='data',
fontsize=14.)
ax.annotate(("$\epsilon_r$= %1i")%(eps[i]),
xy=(-widthg/1.1, (2.*topn[i]+topn[i+1])/3), xycoords='data',
xytext=(-widthg/1.1, (2.*topn[i]+topn[i+1])/3), textcoords='data',
fontsize=14.)
ax.annotate(("$\mu_r$= %1.2f")%(mu[i]),
xy=(-widthg/2., (2.*topn[i]+topn[i+1])/3), xycoords='data',
xytext=(-widthg/2., (2.*topn[i]+topn[i+1])/3), textcoords='data',
fontsize=14.)
ax.plot(widthn,topn[i]*np.ones_like(widthn),color='black')
ax.fill_between(widthn,topn[i],topn[i+1],alpha=0.3,color="none",edgecolor='black', hatch=hatches[(i-1)%10])
#Write the physical properties of the n-th layer and fill it with pattern
ax.plot(widthn,topn[-1]*np.ones_like(widthn),color='black')
ax.fill_between(widthn,topn[-1],z.max(),alpha=0.3,color="none",edgecolor='black', hatch=hatches[(len(topn)-2)%10])
ax.annotate(("$\sigma$ =%3.3f mS/m")%(sig[-1]*10**(3)),
xy=(0., (2.*topn[-1]+z.max())/3), xycoords='data',
xytext=(0., (2.*topn[-1]+z.max())/3), textcoords='data',
fontsize=14.)
ax.annotate(("$\epsilon_r$= %1i")%(eps[-1]),
xy=(-widthg/1.1, (2.*topn[-1]+z.max())/3), xycoords='data',
xytext=(-widthg/1.1, (2.*topn[-1]+z.max())/3), textcoords='data',
fontsize=14.)
ax.annotate(("$\mu_r$= %1.2f")%(mu[-1]),
xy=(-widthg/2., (2.*topn[-1]+z.max())/3), xycoords='data',
xytext=(-widthg/2., (2.*topn[-1]+z.max())/3), textcoords='data',
fontsize=14.)
#plot Trees!
ax.annotate("",
xy=(widthg/2., -1.*z.max()/5.), xycoords='data',
xytext=(widthg/2., 0.), textcoords='data',
arrowprops=dict(arrowstyle='->, head_width=1.2,head_length=1.2',color='green',linewidth=2.)
)
ax.annotate("",
xy=(widthg/2., -3./4.*z.max()/5.), xycoords='data',
xytext=(widthg/2., 0.), textcoords='data',
arrowprops=dict(arrowstyle='->, head_width=1.4,head_length=1.4',color='green',linewidth=2.)
)
ax.annotate("",
xy=(widthg/2., -1./2.*z.max()/5.), xycoords='data',
xytext=(widthg/2., 0.), textcoords='data',
arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.)
)
ax.annotate("",
xy=(1.2*widthg/2., -1.*z.max()/5.), xycoords='data',
xytext=(1.2*widthg/2., 0.), textcoords='data',
arrowprops=dict(arrowstyle='->, head_width=1.2,head_length=1.2',color='green',linewidth=2.)
)
ax.annotate("",
xy=(1.2*widthg/2., -3./4.*z.max()/5.), xycoords='data',
xytext=(1.2*widthg/2., 0.), textcoords='data',
arrowprops=dict(arrowstyle='->, head_width=1.4,head_length=1.4',color='green',linewidth=2.)
)
ax.annotate("",
xy=(1.2*widthg/2., -1./2.*z.max()/5.), xycoords='data',
xytext=(1.2*widthg/2., 0.), textcoords='data',
arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.)
)
ax.annotate("",
xy=(1.5*widthg/2., -1.*z.max()/5.), xycoords='data',
xytext=(1.5*widthg/2., 0.), textcoords='data',
arrowprops=dict(arrowstyle='->, head_width=1.2,head_length=1.2',color='green',linewidth=2.)
)
ax.annotate("",
xy=(1.5*widthg/2., -3./4.*z.max()/5.), xycoords='data',
xytext=(1.5*widthg/2., 0.), textcoords='data',
arrowprops=dict(arrowstyle='->, head_width=1.4,head_length=1.4',color='green',linewidth=2.)
)
ax.annotate("",
xy=(1.5*widthg/2., -1./2.*z.max()/5.), xycoords='data',
xytext=(1.5*widthg/2., 0.), textcoords='data',
arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.)
)
ax.invert_yaxis()
return ax
#Propagate Up and Down component for a certain frequency & evaluate E and H field
def Propagate(f,H,sig,chg,taux,c,mu,eps,n):
sigcm = np.zeros_like(sig,dtype='complex_')
for j in range(1,len(sig)):
sigcm[j]=PCC(sig[j],chg[j],taux[j],c[j],f)
K = k(f, sigcm, mu, eps)
Z = ImpZ(f,mu,K)
EH = np.matrix(np.zeros((2,n+1),dtype = 'complex_'),dtype = 'complex_')
UD = np.matrix(np.zeros((2,n+1),dtype = 'complex_'),dtype = 'complex_')
UD[1,-1] = 1.
for i in range(-2,-(n+2),-1):
UD[:,i] = Tinv(H[i+1],K[i])*Pinv(Z[i])*P(Z[i+1])*UD[:,i+1]
UD = UD/((np.abs(UD[0,:]+UD[1,:])).max())
for j in range(0,n+1):
EH[:,j] = np.matrix([[1.,1,],[-1./Z[j],1./Z[j]]])*UD[:,j]
return UD, EH, Z ,K
#Evaluate the apparent resistivity and phase for a frequency range
def appres(F,H,sig,chg,taux,c,mu,eps,n):
Res = np.zeros_like(F)
Phase = np.zeros_like(F)
App_ImpZ= np.zeros_like(F,dtype='complex_')
for i in range(0,len(F)):
UD,EH,Z ,K = Propagate(F[i],H,sig,chg,taux,c,mu,eps,n)
App_ImpZ[i] = EH[0,1]/EH[1,1]
Res[i] = np.abs(App_ImpZ[i])**2./(mu_0*omega(F[i]))
Phase[i] = np.angle(App_ImpZ[i], deg = True)
return Res,Phase
#Evaluate Up, Down components, E and H field, for a frequency range,
#a discretized depth range and a time range (use to calculate envelope)
def calculateEHzt(F,H,sig,chg,taux,c,mu,eps,n,zsample,tsample):
topc = top(H)
layer = np.zeros(len(zsample),dtype=np.int)-1
Exzt = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
Hyzt = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
Uz = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
Dz = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
UDaux = np.matrix(np.zeros((2,len(zsample)),dtype = 'complex_'),dtype = 'complex_')
for i in range(0,n+1,1):
layer = layer+(zsample>=topc[i])*1
for j in range(0,len(F)):
UD,EH,Z ,K = Propagate(F[j],H,sig,chg,taux,c,mu,eps,n)
for p in range(0,len(zsample)):
UDaux[:,p] = UD_Z(UD[:,layer[p]],zsample[p],topc[layer[p]],K[layer[p]])
for q in range(0,len(tsample)):
Exzt[p,q] = Exzt[p,q] + E_ZT(UDaux[0,p],UDaux[1,p],F[j],tsample[q])/len(F)
Hyzt[p,q] = Hyzt[p,q] + H_ZT(UDaux[0,p],UDaux[1,p],Z[layer[p]],F[j],tsample[q])/len(F)
Uz[p,q] = Uz[p,q] + UDaux[0,p]*np.exp(1j*omega(F[j])*tsample[q])/len(F)
Dz[p,q] = Dz[p,q] + UDaux[1,p]*np.exp(1j*omega(F[j])*tsample[q])/len(F)
return Exzt,Hyzt,Uz,Dz,UDaux,layer
#Function to Plot Apparent Resistivity and Phase
def PlotAppRes(F,H,sig,chg,taux,c,mu,eps,n,fenvelope,PlotEnvelope):
Res, Phase = appres(F,H,sig,chg,taux,c,mu,eps,n)
fig,ax = plt.subplots(1,2,figsize=(16,10))
ax[0].scatter(Res,F,color='black')
ax[0].set_xscale('Log')
ax[0].set_yscale('Log')
ax[0].set_xlim([10.**(np.log10(Res.min())-1.),10.**(np.log10(Res.max())+1.)])
ax[0].set_ylim([F.min(),F.max()])
ax[0].set_xlabel('Apparent Resistivity (Ohm*m)',fontsize=16.,color="black")
ax[0].set_ylabel('Frequency (Hz)',fontsize=16.)
ax[0].grid(which='major')
ax0 = ax[0].twiny()
ax0.set_xlim([0.,90.])
ax0.set_ylim([F.min(),F.max()])
ax0.scatter(Phase,F,color='purple')
ax0.set_xlabel('Phase (Degrees)',fontsize=16.,color="purple")
zc=np.arange(-(H[1:].max()+10)*n,(H[1:].max()+10)*n,10.)
ax[0].tick_params(labelsize=16)
ax[1].tick_params(labelsize=16)
ax0.tick_params(labelsize=16)
if PlotEnvelope:
widthn=np.logspace(np.log10(Res.min())-1., np.log10(Res.max())+1., num=100, endpoint=True, base=10.0)
fenvelope1n=np.ones(100)*fenvelope
ax[0].plot(widthn,fenvelope1n,linestyle='dashed',color='black')
tc=np.arange(0.,1./fenvelope,0.01/(fenvelope))
Exzt,Hyzt,Uz,Dz,UDaux,layer = calculateEHzt(np.array([fenvelope]),H,sig,chg,taux,c,mu,eps,n,zc,tc)
ax1=ax[1].twiny()
ax[1].tick_params(labelsize=16)
ax1.tick_params(labelsize=16)
ax[1].set_xlabel('Amplitude Electric Field E (V/m)',color='blue',fontsize=16)
ax1.set_xlabel('Amplitude Magnetic Field H (A/m)',color='red',fontsize=16)
ax[1].fill_betweenx(zc,np.squeeze(np.asarray(np.real(Exzt.min(axis=1)))),
np.squeeze(np.asarray(np.real(Exzt.max(axis=1)))),
color='blue', alpha=0.1)
ax1.fill_betweenx(zc,np.squeeze(np.asarray(np.real(Hyzt.min(axis=1)))),
np.squeeze(np.asarray(np.real(Hyzt.max(axis=1)))),
color='red', alpha=0.1)
ax[1] = PlotConfiguration(H,sig,eps,mu,ax[1],(1.5*np.abs(Exzt).max()),zc)
ax1.set_xlim([-1.5*np.abs(Hyzt).max(),1.5*np.abs(Hyzt).max()])
ax1.set_xlim([-1.5*np.abs(Hyzt).max(),1.5*np.abs(Hyzt).max()])
else:
print 'No envelop (if True, might be slow)'
ax[1] = PlotConfiguration(H,sig,eps,mu,ax[1],1.,zc)
ax[1].get_xaxis().set_ticks([])
plt.show()
#Interactive MT for Notebook
def PlotAppRes3LayersInteract(h1,h2,sigl1,sigl2,sigl3,mul1,mul2,mul3,epsl1,epsl2,epsl3,PlotEnvelope,F_Envelope):
frangn=frange(-5,5,100.)
sig3= np.array([0.,0.001,0.1, 0.001])
thick3 = np.array([120000.,50.,50.])
eps3=np.array([1.,1.,1.,1])
mu3=np.array([1.,1.,1.,1])
chg3=np.array([0.,0.1,0.,0.2])
chg3_0=np.array([0.,0.1,0.,0.])
taux3=np.array([0.,0.1,0.,0.1])
c3=np.array([1.,1.,1.,1.])
sig3[1]=sigl1
sig3[1]=10.**sig3[1]
sig3[2]=sigl2
sig3[2]=10.**sig3[2]
sig3[3]=sigl3
sig3[3]=10.**sig3[3]
mu3[1]=mul1
mu3[2]=mul2
mu3[3]=mul3
eps3[1]=epsl1
eps3[2]=epsl2
eps3[3]=epsl3
thick3[1]=h1
thick3[2]=h2
PlotAppRes(frangn,thick3,sig3,chg3_0,taux3,c3,mu3,eps3,3,F_Envelope,PlotEnvelope)
def run(n=3,plotIt=True):
# something to make a plot
F = frange(-5.,5.,20)
H = thick(50.,100.,n)
sign = sig(-5.,0.,n)
mun = mu(1.,2.,n)
epsn = eps(1.,9.,n)
chg = np.zeros_like(sign)
taux = np.zeros_like(sign)
c = np.zeros_like(sign)
Res, Phase = appres(F,H,sign,chg,taux,c,mun,epsn,n)
if plotIt:
PlotAppRes(F, H, sign, chg, taux, c, mun, epsn, n, fenvelope=1000., PlotEnvelope=True)
return Res, Phase
if __name__ == '__main__':
run()
+8 -7
View File
@@ -2,7 +2,7 @@
# Import
import SimPEG as simpeg
from SimPEG import MT
from SimPEG import NSEM
import numpy as np
try:
from pymatsolver import MumpsSolver as Solver
@@ -12,7 +12,7 @@ except:
def run(plotIt=True, nFreq=1):
"""
MT: 3D: Forward
===============
=======================
Forward model 3D MT data.
@@ -37,24 +37,25 @@ def run(plotIt=True, nFreq=1):
for loc in rx_loc:
# NOTE: loc has to be a (1,3) np.ndarray otherwise errors accure
for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi','tzxr','tzxi','tzyr','tzyi']:
rxList.append(MT.Rx(simpeg.mkvc(loc,2).T,rxType))
rxList.append(NSEM.Rx(simpeg.mkvc(loc,2).T,rxType))
# Source list
srcList =[]
for freq in np.logspace(3,-3,nFreq):
srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq))
# Survey MT
survey = MT.Survey(srcList)
survey = NSEM.Survey(srcList)
## Setup the problem object
problem = MT.Problem3D.eForm_ps(M, sigmaPrimary=sigBG, Solver=Solver)
problem = NSEM.Problem3D_ePrimSec(M, sigmaPrimary=sigBG)
problem.pair(survey)
problem.Solver = Solver
# Calculate the data
fields = problem.fields(sig)
dataVec = survey.eval(fields)
# Make the data
mtData = MT.Data(survey, dataVec)
mtData = NSEM.Data(survey,dataVec)
# Add plots
if plotIt:
pass
-62
View File
@@ -1,62 +0,0 @@
from SimPEG import Mesh, Maps, np
def run(plotIt=True):
"""
Maps: ComboMaps
===============
We will use an example where we want a 1D layered earth as
our model, but we want to map this to a 2D discretization to do our forward
modeling. We will also assume that we are working in log conductivity still,
so after the transformation we want to map to conductivity space.
To do this we will introduce the vertical 1D map (:class:`SimPEG.Maps.SurjectVertical1D`),
which does the first part of what we just described. The second part will be
done by the :class:`SimPEG.Maps.ExpMap` described above.
.. code-block:: python
:linenos:
M = Mesh.TensorMesh([7,5])
v1dMap = Maps.SurjectVertical1D(M)
expMap = Maps.ExpMap(M)
myMap = expMap * v1dMap
m = np.r_[0.2,1,0.1,2,2.9] # only 5 model parameters!
sig = myMap * m
If you noticed, it was pretty easy to combine maps. What is even cooler is
that the derivatives also are made for you (if everything goes right).
Just to be sure that the derivative is correct, you should always run the test
on the mapping that you create.
"""
M = Mesh.TensorMesh([7,5])
v1dMap = Maps.SurjectVertical1D(M)
expMap = Maps.ExpMap(M)
myMap = expMap * v1dMap
m = np.r_[0.2,1,0.1,2,2.9] # only 5 model parameters!
sig = myMap * m
if not plotIt: return
import matplotlib.pyplot as plt
figs, axs = plt.subplots(1,2)
axs[0].plot(m, M.vectorCCy, 'b-o')
axs[0].set_title('Model')
axs[0].set_ylabel('Depth, y')
axs[0].set_xlabel('Value, $m_i$')
axs[0].set_xlim(0,3)
axs[0].set_ylim(0,1)
clbar = plt.colorbar(M.plotImage(sig,ax=axs[1],grid=True,gridOpts=dict(color='grey'))[0])
axs[1].set_title('Physical Property')
axs[1].set_ylabel('Depth, y')
clbar.set_label('$\sigma = \exp(\mathbf{P}m)$')
plt.tight_layout()
plt.show()
if __name__ == '__main__':
run()
-41
View File
@@ -1,41 +0,0 @@
from SimPEG import Mesh, Maps, Utils
def run(plotIt=True):
"""
Maps: Mesh2Mesh
===============
This mapping allows you to go from one mesh to another.
"""
M = Mesh.TensorMesh([100,100])
h1 = Utils.meshTensor([(6,7,-1.5),(6,10),(6,7,1.5)])
h1 = h1/h1.sum()
M2 = Mesh.TensorMesh([h1,h1])
V = Utils.ModelBuilder.randomModel(M.vnC, seed=79, its=50)
v = Utils.mkvc(V)
modh = Maps.Mesh2Mesh([M,M2])
modH = Maps.Mesh2Mesh([M2,M])
H = modH * v
h = modh * H
if not plotIt: return
import matplotlib.pyplot as plt
ax = plt.subplot(131)
M.plotImage(v, ax=ax)
ax.set_title('Fine Mesh (Original)')
ax = plt.subplot(132)
M2.plotImage(H,clim=[0,1],ax=ax)
ax.set_title('Course Mesh')
ax = plt.subplot(133)
M.plotImage(h,clim=[0,1],ax=ax)
ax.set_title('Fine Mesh (Interpolated)')
plt.show()
if __name__ == '__main__':
run()
+7 -9
View File
@@ -2,12 +2,8 @@ from SimPEG import *
from SimPEG.Utils import surface2ind_topo
def run(plotIt=True, nx=5, ny=5):
def run(plotIt=False, nx = 5, ny = 5):
"""
Utils: surface2ind_topo
=======================
Here we show how to use :code:`Utils.surface2ind_topo` to identify cells below
a topographic surface.
@@ -17,25 +13,27 @@ def run(plotIt=True, nx=5, ny=5):
xtopo = np.linspace(mesh.gridN[:,0].min(), mesh.gridN[:,0].max())
topo = 0.4*np.sin(xtopo*5) # define a topographic surface
Topo = np.hstack([Utils.mkvc(xtopo,2), Utils.mkvc(topo,2)]) #make it an array
Topo = np.hstack([Utils.mkvc(xtopo,2),Utils.mkvc(topo,2)]) #make it an array
indcc = surface2ind_topo(mesh, Topo, 'CC')
indcc = surface2ind_topo(mesh, Topo,'CC')
if plotIt:
from matplotlib.pylab import plt
from scipy.interpolate import interp1d
fig, ax = plt.subplots(1,1, figsize=(6,6))
fig, ax = plt.subplots(1,1,figsize=(6,6))
mesh.plotGrid(ax=ax, nodes=True, centers=True)
ax.plot(xtopo,topo,'k',linewidth=1)
# ax.plot(mesh.vectorNx, interp1d(xtopo,topo)(mesh.vectorNx),'--k',linewidth=3)
ax.plot(mesh.vectorCCx, interp1d(xtopo,topo)(mesh.vectorCCx),'--k',linewidth=3)
aveN2CC = Utils.sdiag(mesh.aveN2CC.T.sum(1))*mesh.aveN2CC.T
a = aveN2CC * indcc
a[a > 0] = 1.
a[a < 0.25] = np.nan
a = a.reshape(mesh.vnN, order='F')
masked_array = np.ma.array(a, mask=np.isnan(a))
ax.pcolor(mesh.vectorNx,mesh.vectorNy,masked_array.T, cmap=plt.cm.gray, alpha=0.2)
ax.pcolor(mesh.vectorNx,mesh.vectorNy,masked_array.T, cmap = plt.cm.gray,alpha=0.2)
plt.show()
+18 -19
View File
@@ -1,30 +1,29 @@
# 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 Mesh_QuadTree_Creation
import EM_TDEM_1D_Inversion
import Mesh_QuadTree_FaceDiv
import Mesh_Tensor_Creation
import FLOW_Richards_1D_Celia1990
import DC_Forward_PseudoSection
import Mesh_Operators_CahnHilliard
import Mesh_Basic_Types
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 EM_Schenkel_Morrison_Casing
import MT_3D_Foward
import Mesh_Basic_ForwardDC
import MT_1D_ForwardAndInversion
import Utils_surface2ind_topo
import MT_1D_analytic_nlayer_Earth
import EM_FDEM_Analytic_MagDipoleWholespace
import Mesh_Basic_PlotImage
import DC_Analytic_Dipole
import Mesh_QuadTree_HangingNodes
__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Inversion_IRLS", "Inversion_Linear", "Maps_ComboMaps", "Maps_Mesh2Mesh", "Mesh_Basic_ForwardDC", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward", "Utils_surface2ind_topo"]
__examples__ = ["EM_FDEM_1D_Inversion", "Mesh_QuadTree_Creation", "EM_TDEM_1D_Inversion", "Mesh_QuadTree_FaceDiv", "Mesh_Tensor_Creation", "FLOW_Richards_1D_Celia1990", "DC_Forward_PseudoSection", "Mesh_Operators_CahnHilliard", "Mesh_Basic_Types", "Inversion_IRLS", "Inversion_Linear", "EM_Schenkel_Morrison_Casing", "MT_3D_Foward", "Mesh_Basic_ForwardDC", "MT_1D_ForwardAndInversion", "Utils_surface2ind_topo", "MT_1D_analytic_nlayer_Earth", "EM_FDEM_Analytic_MagDipoleWholespace", "Mesh_Basic_PlotImage", "DC_Analytic_Dipole", "Mesh_QuadTree_HangingNodes"]
##### AUTOIMPORTS #####
@@ -40,7 +39,7 @@ if __name__ == '__main__':
# Create the examples dir in the docs folder.
fName = os.path.realpath(__file__)
docExamplesDir = os.path.sep.join(fName.split(os.path.sep)[:-3] + ['docs', 'content', 'examples'])
docExamplesDir = os.path.sep.join(fName.split(os.path.sep)[:-3] + ['docs', 'examples'])
shutil.rmtree(docExamplesDir)
os.makedirs(docExamplesDir)
@@ -97,12 +96,12 @@ if __name__ == '__main__':
from SimPEG import Examples
Examples.%s.run()
.. literalinclude:: ../../../SimPEG/Examples/%s.py
.. literalinclude:: ../../SimPEG/Examples/%s.py
:language: python
:linenos:
"""%(name,doc,name,name)
rst = os.path.sep.join((filePath.split(os.path.sep)[:-3] + ['docs', 'content', 'examples', name + '.rst']))
rst = os.path.sep.join((filePath.split(os.path.sep)[:-3] + ['docs', 'examples', name + '.rst']))
print 'Creating: %s.rst'%name
f = open(rst, 'w')
+2 -2
View File
@@ -31,7 +31,7 @@ class NonLinearMap(object):
"""
:param numpy.array u: fields
:param numpy.array m: model
:rtype: scipy.sparse.csr_matrix
:rtype: scipy.csr_matrix
:return: derivative of transformed model
The *transform* changes the model into the physical property.
@@ -44,7 +44,7 @@ class NonLinearMap(object):
"""
:param numpy.array u: fields
:param numpy.array m: model
:rtype: scipy.sparse.csr_matrix
:rtype: scipy.csr_matrix
:return: derivative of transformed model
The *transform* changes the model into the physical property.
-132
View File
@@ -1,132 +0,0 @@
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
class BaseMTProblem(BaseFDEMProblem):
"""
Base class for all Natural source problems.
"""
def __init__(self, mesh, **kwargs):
BaseFDEMProblem.__init__(self, mesh, **kwargs)
Utils.setKwargs(self, **kwargs)
# Set the default pairs of the problem
surveyPair = Survey
dataPair = Data
fieldsPair = BaseMTFields
# Set the solver
Solver = SimpegSolver
solverOpts = {}
verbose = False
# Notes:
# Use the forward and devs from BaseFDEMProblem
# Might need to add more stuff here.
## NEED to clean up the Jvec and Jtvec to use Zero and Identities for None components.
def Jvec(self, m, v, f=None):
"""
Function to calculate the data sensitivities dD/dm times a vector.
:param numpy.ndarray m (nC, 1) - conductive model
:param numpy.ndarray v (nC, 1) - random vector
:param MTfields object (optional) - MT fields object, if not given it is calculated
:rtype: MTdata object
:return: Data sensitivities wrt m
"""
# Calculate the fields
if f is None:
f= self.fields(m)
# Set current model
self.curModel = m
# Initiate the Jv object
Jv = self.dataPair(self.survey)
# Loop all the frequenies
for freq in self.survey.freqs:
dA_du = self.getA(freq) #
dA_duI = self.Solver(dA_du, **self.solverOpts)
for src in self.survey.getSrcByFreq(freq):
# We need fDeriv_m = df/du*du/dm + df/dm
# Construct du/dm, it requires a solve
# NOTE: need to account for the 2 polarizations in the derivatives.
f_src = f[src,:]
# dA_dm and dRHS_dm should be of size nE,2, so that we can multiply by dA_duI. The 2 columns are each of the polarizations.
dA_dm = self.getADeriv_m(freq, f_src, v) # Size: nE,2 (u_px,u_py) in the columns.
dRHS_dm = self.getRHSDeriv_m(freq, v) # Size: nE,2 (u_px,u_py) in the columns.
if dRHS_dm is None:
du_dm = dA_duI * ( -dA_dm )
else:
du_dm = dA_duI * ( -dA_dm + dRHS_dm )
# Calculate the projection derivatives
for rx in src.rxList:
# Get the projection derivative
# v should be of size 2*nE (for 2 polarizations)
PDeriv_u = lambda t: rx.evalDeriv(src, self.mesh, f, t) # wrt u, we don't have have PDeriv wrt m
Jv[src, rx] = PDeriv_u(mkvc(du_dm))
dA_duI.clean()
# Return the vectorized sensitivities
return mkvc(Jv)
def Jtvec(self, m, v, f=None):
"""
Function to calculate the transpose of the data sensitivities (dD/dm)^T times a vector.
:param numpy.ndarray m (nC, 1) - conductive model
:param numpy.ndarray v (nD, 1) - vector
:param MTfields object u (optional) - MT fields object, if not given it is calculated
:rtype: MTdata object
:return: Data sensitivities wrt m
"""
if f is None:
f = self.fields(m)
self.curModel = m
# Ensure v is a data object.
if not isinstance(v, self.dataPair):
v = self.dataPair(self.survey, v)
Jtv = np.zeros(m.size)
for freq in self.survey.freqs:
AT = self.getA(freq).T
ATinv = self.Solver(AT, **self.solverOpts)
for src in self.survey.getSrcByFreq(freq):
ftype = self._fieldType + 'Solution'
f_src = f[src, :]
for rx in src.rxList:
# Get the adjoint evalDeriv
# PTv needs to be nE,
PTv = rx.evalDeriv(src, self.mesh, f, mkvc(v[src, rx],2), adjoint=True) # wrt u, need possibility wrt m
# Get the
dA_duIT = ATinv * PTv
dA_dmT = self.getADeriv_m(freq, f_src, mkvc(dA_duIT), adjoint=True)
dRHS_dmT = self.getRHSDeriv_m(freq, mkvc(dA_duIT), adjoint=True)
# Make du_dmT
if dRHS_dmT is None:
du_dmT = -dA_dmT
else:
du_dmT = -dA_dmT + dRHS_dmT
# Select the correct component
# du_dmT needs to be of size nC,
real_or_imag = rx.projComp
if real_or_imag == 'real':
Jtv += du_dmT.real
elif real_or_imag == 'imag':
Jtv += -du_dmT.real
else:
raise Exception('Must be real or imag')
# Clean the factorization, clear memory.
ATinv.clean()
return Jtv
-291
View File
@@ -1,291 +0,0 @@
from SimPEG.EM.Utils import omega
from SimPEG import mkvc
from scipy.constants import mu_0
from SimPEG.MT.BaseMT import BaseMTProblem
from SimPEG.MT.SurveyMT import Survey, Data
from SimPEG.MT.FieldsMT import Fields1D_e
from SimPEG.MT.Utils.MT1Danalytic import getEHfields
import numpy as np
import multiprocessing, sys, time
class eForm_psField(BaseMTProblem):
"""
A MT problem soving a e formulation and primary/secondary fields decomposion.
By eliminating the magnetic flux density using
.. math ::
\mathbf{b} = \\frac{1}{i \omega}\\left(-\mathbf{C} \mathbf{e} \\right)
we can write Maxwell's equations as a second order system in \\\(\\\mathbf{e}\\\) only:
.. math ::
\\left(\mathbf{C}^T \mathbf{M^e_{\mu^{-1}}} \mathbf{C} + i \omega \mathbf{M^f_\sigma}] \mathbf{e}_{s} =& i \omega \mathbf{M^f_{\delta \sigma}} \mathbf{e}_{p}
which we solve for \\\(\\\mathbf{e_s}\\\). The total field \\\mathbf{e}\\ = \\\mathbf{e_p}\\ + \\\mathbf{e_s}\\.
The primary field is estimated from a background model (commonly half space ).
"""
# From FDEMproblem: Used to project the fields. Currently not used for MTproblem.
_fieldType = 'e_1d'
_eqLocs = 'EF'
_sigmaPrimary = None
def __init__(self, mesh, **kwargs):
BaseMTProblem.__init__(self, mesh, **kwargs)
self.fieldsPair = Fields1D_e
# self._sigmaPrimary = sigmaPrimary
@property
def MeMui(self):
"""
Edge inner product matrix
"""
if getattr(self, '_MeMui', None) is None:
self._MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
return self._MeMui
@property
def MfSigma(self):
"""
Edge inner product matrix
"""
if getattr(self, '_MfSigma', None) is None:
self._MfSigma = self.mesh.getFaceInnerProduct(self.curModel.sigma)
return self._MfSigma
@property
def sigmaPrimary(self):
"""
A background model, use for the calculation of the primary fields.
"""
return self._sigmaPrimary
@sigmaPrimary.setter
def sigmaPrimary(self, val):
# Note: TODO add logic for val, make sure it is the correct size.
self._sigmaPrimary = val
def getA(self, freq):
"""
Function to get the A matrix.
:param float freq: Frequency
:rtype: scipy.sparse.csr_matrix
:return: A
"""
# Note: need to use the code above since in the 1D problem I want
# e to live on Faces(nodes) and h on edges(cells). Might need to rethink this
# Possible that _fieldType and _eqLocs can fix this
MeMui = self.MeMui
MfSigma = self.MfSigma
C = self.mesh.nodalGrad
# Make A
A = C.T*MeMui*C + 1j*omega(freq)*MfSigma
# Either return full or only the inner part of A
return A
def getADeriv_m(self, freq, u, v, adjoint=False):
"""
The derivative of A wrt sigma
"""
dsig_dm = self.curModel.sigmaDeriv
MeMui = self.MeMui
#
u_src = u['e_1dSolution']
dMfSigma_dm = self.mesh.getFaceInnerProductDeriv(self.curModel.sigma)(u_src) * self.curModel.sigmaDeriv
if adjoint:
return 1j * omega(freq) * ( dMfSigma_dm.T * v )
# Note: output has to be nN/nF, not nC/nE.
# v should be nC
return 1j * omega(freq) * ( dMfSigma_dm * v )
def getRHS(self, freq):
"""
Function to return the right hand side for the system.
:param float freq: Frequency
:rtype: numpy.ndarray (nF, 1), numpy.ndarray (nF, 1)
:return: RHS for 1 polarizations, primary fields
"""
# Get sources for the frequncy(polarizations)
Src = self.survey.getSrcByFreq(freq)[0]
S_e = Src.S_e(self)
return -1j * omega(freq) * S_e
def getRHSDeriv_m(self, freq, v, adjoint=False):
"""
The derivative of the RHS wrt sigma
"""
Src = self.survey.getSrcByFreq(freq)[0]
S_eDeriv = Src.S_eDeriv_m(self, v, adjoint)
return -1j * omega(freq) * S_eDeriv
def fields(self, m):
'''
Function to calculate all the fields for the model m.
:param np.ndarray (nC,) m: Conductivity model
'''
# Set the current model
self.curModel = m
F = Fields1D_e(self.mesh, self.survey)
for freq in self.survey.freqs:
if self.verbose:
startTime = time.time()
print 'Starting work for {:.3e}'.format(freq)
sys.stdout.flush()
A = self.getA(freq)
rhs = self.getRHS(freq)
Ainv = self.Solver(A, **self.solverOpts)
e_s = Ainv * rhs
# Store the fields
Src = self.survey.getSrcByFreq(freq)[0]
# NOTE: only store the e_solution(secondary), all other components calculated in the fields object
F[Src, 'e_1dSolution'] = e_s[:,-1] # Only storing the yx polarization as 1d
# Note curl e = -iwb so b = -curl e /iw
# 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)
sys.stdout.flush()
return F
# Note this is not fully functional.
# Missing:
# Fields class corresponding to the fields
# Update Jvec and Jtvec to include all the derivatives components
# Other things ...
class eForm_TotalField(BaseMTProblem):
"""
A MT problem solving a e formulation and a Total bondary domain decompostion.
Solves the equation:
Math:
"""
# From FDEMproblem: Used to project the fields. Currently not used for MTproblem.
_fieldType = 'e'
_eqLocs = 'EF'
def __init__(self, mesh, **kwargs):
BaseMTProblem.__init__(self, mesh, **kwargs)
@property
def MeMui(self):
"""
Edge inner product matrix
"""
if getattr(self, '_MeMui', None) is None:
self._MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
return self._MeMui
@property
def MfSigma(self):
"""
Edge inner product matrix
"""
if getattr(self, '_MfSigma', None) is None:
self._MfSigma = self.mesh.getFaceInnerProduct(self.curModel.sigma)
return self._MfSigma
def getA(self, freq, full=False):
"""
Function to get the A matrix.
:param float freq: Frequency
:param logic full: Return full A or the inner part
:rtype: scipy.sparse.csr_matrix
:return: A
"""
MeMui = self.MeMui
MfSigma = self.MfSigma
# Note: need to use the code above since in the 1D problem I want
# e to live on Faces(nodes) and h on edges(cells). Might need to rethink this
# Possible that _fieldType and _eqLocs can fix this
# MeMui = self.MfMui
# MfSigma = self.MfSigma
C = self.mesh.nodalGrad
# Make A
A = C.T*MeMui*C + 1j*omega(freq)*MfSigma
# Either return full or only the inner part of A
if full:
return A
else:
return A[1:-1,1:-1]
def getADeriv_m(self, freq, u, v, adjoint=False):
raise NotImplementedError('getADeriv is not implemented')
def getRHS(self, freq):
"""
Function to return the right hand side for the system.
:param float freq: Frequency
:rtype: numpy.ndarray (nE, 2), numpy.ndarray (nE, 2)
:return: RHS for both polarizations, primary fields
"""
# Get sources for the frequency
# NOTE: Need to use the source information, doesn't really apply in 1D
src = self.survey.getSrcByFreq(freq)
# Get the full A
A = self.getA(freq,full=True)
# Define the outer part of the solution matrix
Aio = A[1:-1,[0,-1]]
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
## Note: The analytic solution is derived with e^iwt
eBC = np.r_[Etot[0],Etot[-1]]
# The right hand side
return -Aio*eBC, eBC
def getRHSderiv_m(self, freq, backSigma, u, v, adjoint=False):
raise NotImplementedError('getRHSDeriv not implemented yet')
return None
def fields(self, m):
'''
Function to calculate all the fields for the model m.
:param np.ndarray (nC,) m: Conductivity model
:param np.ndarray (nC,) m_back: Background conductivity model
'''
self.curModel = m
# RHS, CalcFields = self.getRHS(freq,m_back), self.calcFields
F = Fields1D_e(self.mesh, self.survey)
for freq in self.survey.freqs:
if self.verbose:
startTime = time.time()
print 'Starting work for {:.3e}'.format(freq)
sys.stdout.flush()
A = self.getA(freq)
rhs, e_o = self.getRHS(freq)
Ainv = self.Solver(A, **self.solverOpts)
e_i = Ainv * rhs
e = mkvc(np.r_[e_o[0], e_i, e_o[1]],2)
# Store the fields
Src = self.survey.getSrcByFreq(freq)
# NOTE: only store e fields
F[Src, 'e_1dSolution'] = e[:,0]
if self.verbose:
print 'Ran for {:f} seconds'.format(time.time()-startTime)
sys.stdout.flush()
return F
-1
View File
@@ -1 +0,0 @@
from Probs import eForm_TotalField, eForm_psField
View File
-1
View File
@@ -1 +0,0 @@
pass
-138
View File
@@ -1,138 +0,0 @@
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
from SimPEG.MT.BaseMT import BaseMTProblem
from SimPEG.MT.SurveyMT import Survey, Data
from SimPEG.MT.FieldsMT import Fields3D_e
import multiprocessing, sys, time
class eForm_ps(BaseMTProblem):
"""
A MT problem solving a e formulation and a primary/secondary fields decompostion.
By eliminating the magnetic flux density using
.. math ::
\mathbf{b} = \\frac{1}{i \omega}\\left(-\mathbf{C} \mathbf{e} \\right)
we can write Maxwell's equations as a second order system in \\\(\\\mathbf{e}\\\) only:
.. math ::
\\left(\mathbf{C}^T \mathbf{M^f_{\mu^{-1}}} \mathbf{C} + i \omega \mathbf{M^e_\sigma}] \mathbf{e}_{s} =& i \omega \mathbf{M^e_{\delta \sigma}} \mathbf{e}_{p}
which we solve for \\\(\\\mathbf{e_s}\\\). The total field \\\mathbf{e}\\ = \\\mathbf{e_p}\\ + \\\mathbf{e_s}\\.
The primary field is estimated from a background model (commonly as a 1D model).
"""
# From FDEMproblem: Used to project the fields. Currently not used for MTproblem.
_fieldType = 'e'
_eqLocs = 'FE'
fieldsPair = Fields3D_e
_sigmaPrimary = None
def __init__(self, mesh, **kwargs):
BaseMTProblem.__init__(self, mesh, **kwargs)
@property
def sigmaPrimary(self):
"""
A background model, use for the calculation of the primary fields.
"""
return self._sigmaPrimary
@sigmaPrimary.setter
def sigmaPrimary(self, val):
# Note: TODO add logic for val, make sure it is the correct size.
self._sigmaPrimary = val
def getA(self, freq):
"""
Function to get the A system.
:param float freq: Frequency
:rtype: scipy.sparse.csr_matrix
:return: A
"""
Mmui = self.MfMui
Msig = self.MeSigma
C = self.mesh.edgeCurl
return C.T*Mmui*C + 1j*omega(freq)*Msig
def getADeriv_m(self, freq, u, v, adjoint=False):
"""
Calculate the derivative of A wrt m.
"""
# This considers both polarizations and returns a nE,2 matrix for each polarization
if adjoint:
dMe_dsigV = sp.hstack(( self.MeSigmaDeriv( u['e_pxSolution'] ).T, self.MeSigmaDeriv(u['e_pySolution'] ).T ))*v
else:
# Need a nE,2 matrix to be returned
dMe_dsigV = np.hstack(( mkvc(self.MeSigmaDeriv( u['e_pxSolution'] )*v,2), mkvc( self.MeSigmaDeriv(u['e_pySolution'] )*v,2) ))
return 1j * omega(freq) * dMe_dsigV
def getRHS(self, freq):
"""
Function to return the right hand side for the system.
:param float freq: Frequency
:rtype: numpy.ndarray (nE, 2), numpy.ndarray (nE, 2)
:return: RHS for both polarizations, primary fields
"""
# Get sources for the frequncy(polarizations)
Src = self.survey.getSrcByFreq(freq)[0]
S_e = Src.S_e(self)
return -1j * omega(freq) * S_e
def getRHSDeriv_m(self, freq, v, adjoint=False):
"""
The derivative of the RHS with respect to sigma
"""
Src = self.survey.getSrcByFreq(freq)[0]
S_eDeriv = Src.S_eDeriv_m(self, v, adjoint)
return -1j * omega(freq) * S_eDeriv
def fields(self, m):
'''
Function to calculate all the fields for the model m.
:param np.ndarray (nC,) m: Conductivity model
'''
# Set the current model
self.curModel = m
F = Fields3D_e(self.mesh, self.survey)
for freq in self.survey.freqs:
if self.verbose:
startTime = time.time()
print 'Starting work for {:.3e}'.format(freq)
sys.stdout.flush()
A = self.getA(freq)
rhs = self.getRHS(freq)
# Solve the system
Ainv = self.Solver(A, **self.solverOpts)
e_s = Ainv * rhs
# Store the fields
Src = self.survey.getSrcByFreq(freq)[0]
# Store the fieldss
F[Src, 'e_pxSolution'] = e_s[:,0]
F[Src, 'e_pySolution'] = e_s[:,1]
# Note curl e = -iwb so b = -curl/iw
if self.verbose:
print 'Ran for {:f} seconds'.format(time.time()-startTime)
sys.stdout.flush()
Ainv.clean()
return F
-1
View File
@@ -1 +0,0 @@
from Probs import eForm_ps
-4
View File
@@ -1,4 +0,0 @@
from MT1Dsolutions import * # Add the names of the functions
from MT1Danalytic import *
from dataUtils import *
from ediFilesUtils import *
-46
View File
@@ -1,46 +0,0 @@
import SimPEG as simpeg, numpy as np
def homo1DModelSource(mesh,freq,m_back):
'''
Function that calculates and return background fields for a 3D mesh and model.
The calculuations use 1D field solution for a vertical slice throught model (south-western most column),
which is assigned at the fields everywhere for the respective polarizations.2
:param Simpeg mesh object mesh: Holds information on the discretization
:param float freq: The frequency to solve at
:param np.array m_back: Background model of conductivity to base the calculations on.
:rtype: numpy.ndarray (mesh.nE,2)
:return: eBG_bp, E fields for the background model at both polarizations.
'''
# import
from SimPEG.MT.Utils import get1DEfields
# Get a 1d solution for a halfspace background
mesh1d = simpeg.Mesh.TensorMesh([mesh.hz],np.array([mesh.x0[2]]))
# Note: Everything is using e^iwt
e0_1d = get1DEfields(mesh1d,mesh.r(m_back,'CC','CC','M')[0,0,:],freq)
# Setup x (east) polarization (_x)
ex_px = np.zeros(mesh.vnEx,dtype=complex)
ey_px = np.zeros((mesh.nEy,1),dtype=complex)
ez_px = np.zeros((mesh.nEz,1),dtype=complex)
# Assign the source to ex_x
for i in np.arange(mesh.vnEx[0]):
for j in np.arange(mesh.vnEx[1]):
ex_px[i,j,:] = -e0_1d
eBG_px = np.vstack((simpeg.Utils.mkvc(ex_px,2),ey_px,ez_px))
# Setup y (north) polarization (_py)
ex_py = np.zeros((mesh.nEx,1), dtype='complex128')
ey_py = np.zeros(mesh.vnEy, dtype='complex128')
ez_py = np.zeros((mesh.nEz,1), dtype='complex128')
# Assign the source to ey_py
for i in np.arange(mesh.vnEy[0]):
for j in np.arange(mesh.vnEy[1]):
ey_py[i,j,:] = e0_1d
# ey_py[1:-1,1:-1,1:-1] = 0
eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py))
# Return the electric fields
eBG_bp = np.hstack((eBG_px,eBG_py))
return eBG_bp
-5
View File
@@ -1,5 +0,0 @@
import Utils
from SurveyMT import Rx, Survey, Data
from FieldsMT import Fields1D_e, Fields3D_e
import Problem1D, Problem2D, Problem3D
import SrcMT
+167 -313
View File
@@ -1,6 +1,4 @@
import Utils
import numpy as np
import scipy.sparse as sp
import Utils, numpy as np, scipy.sparse as sp
from scipy.sparse.linalg import LinearOperator
from Tests import checkDerivative
from PropMaps import PropMap, Property
@@ -8,7 +6,6 @@ from numpy.polynomial import polynomial
from scipy.interpolate import UnivariateSpline
import warnings
class IdentityMap(object):
"""
SimPEG Map
@@ -20,11 +17,10 @@ class IdentityMap(object):
Utils.setKwargs(self, **kwargs)
if nP is not None:
assert type(nP) in [int, long], 'Number of parameters '
'must be an integer.'
assert type(nP) in [int, long], ' Number of parameters must be an integer.'
self.mesh = mesh
self._nP = nP
self._nP = nP
@property
def nP(self):
@@ -45,8 +41,8 @@ class IdentityMap(object):
If this is a meshless mapping (i.e. nP is defined independently)
the shape will be the the shape (nP,nP).
:rtype: tuple
:return: shape of the operator as a tuple (int,int)
:rtype: (int,int)
:return: shape of the operator as a tuple
"""
if self._nP is not None:
return (self.nP, self.nP)
@@ -54,14 +50,14 @@ class IdentityMap(object):
return ('*', self.nP)
return (self.mesh.nC, self.nP)
def _transform(self, m):
"""
Changes the model into the physical property.
.. note::
This can be called by the __mul__ property against a
:meth:numpy.ndarray.
This can be called by the __mul__ property against a numpy.ndarray.
:param numpy.array m: model
:rtype: numpy.array
@@ -85,25 +81,22 @@ class IdentityMap(object):
"""
raise NotImplementedError('The transformInverse is not implemented.')
def deriv(self, m, v=None):
def deriv(self, m):
"""
The derivative of the transformation.
:param numpy.array m: model
:rtype: scipy.sparse.csr_matrix
:rtype: scipy.csr_matrix
:return: derivative of transformed model
"""
if v is not None:
return v
return sp.identity(self.nP)
def test(self, m=None, **kwargs):
"""Test the derivative of the mapping.
:param numpy.array m: model
:param kwargs: key word arguments of
:meth:`SimPEG.Tests.checkDerivative`
:param kwargs: key word arguments of :meth:`SimPEG.Tests.checkDerivative`
:rtype: bool
:return: passed the test?
@@ -113,52 +106,26 @@ class IdentityMap(object):
m = abs(np.random.rand(self.nP))
if 'plotIt' not in kwargs:
kwargs['plotIt'] = False
return checkDerivative(lambda m : [self * m, self.deriv(m)], m, num=4,
**kwargs)
def testVec(self, m=None, **kwargs):
"""Test the derivative of the mapping times a vector.
:param numpy.array m: model
:param kwargs: key word arguments of
:meth:`SimPEG.Tests.checkDerivative`
:rtype: bool
:return: passed the test?
"""
print 'Testing %s' % str(self)
if m is None:
m = abs(np.random.rand(self.nP))
if 'plotIt' not in kwargs:
kwargs['plotIt'] = False
return checkDerivative(lambda m: [self*m, lambda x: self.deriv(m, x)],
m, num=4, **kwargs)
return checkDerivative(lambda m : [self * m, self.deriv(m)], m, num=4, **kwargs)
def _assertMatchesPair(self, pair):
assert (isinstance(self, pair) or
isinstance(self, ComboMap) and isinstance(self.maps[0], pair)
), ("Mapping object must be an instance of a %s"
" class."% (pair.__name__))
isinstance(self, ComboMap) and isinstance(self.maps[0], pair)
), "Mapping object must be an instance of a %s class."%(pair.__name__)
def __mul__(self, val):
if isinstance(val, IdentityMap):
if (not (self.shape[1] == '*' or val.shape[0] == '*') and not
self.shape[1] == val.shape[0]):
raise ValueError('Dimension mismatch in %s and %s.'
% (str(self), str(val)))
if not (self.shape[1] == '*' or val.shape[0] == '*') and not self.shape[1] == val.shape[0]:
raise ValueError('Dimension mismatch in %s and %s.' % (str(self), str(val)))
return ComboMap([self, val])
elif isinstance(val, np.ndarray):
if not self.shape[1] == '*' and not self.shape[1] == val.shape[0]:
raise ValueError('Dimension mismatch in %s and np.ndarray%s.'
% (str(self), str(val.shape)))
raise ValueError('Dimension mismatch in %s and np.ndarray%s.' % (str(self), str(val.shape)))
return self._transform(val)
raise Exception('Unrecognized data type to multiply. '
'Try a map or a numpy.ndarray!')
raise Exception('Unrecognized data type to multiply. Try a map or a numpy.ndarray!')
def __str__(self):
return "%s(%s,%s)" % (self.__class__.__name__, self.shape[0],
self.shape[1])
return "%s(%s,%s)" % (self.__class__.__name__, self.shape[0], self.shape[1])
class ComboMap(IdentityMap):
@@ -169,17 +136,11 @@ class ComboMap(IdentityMap):
self.maps = []
for ii, m in enumerate(maps):
assert isinstance(m, IdentityMap), "Unrecognized data type, "
"inherit from an IdentityMap or ComboMap!"
if (ii > 0 and not (self.shape[1] == '*' or m.shape[0] == '*') and
not self.shape[1] == m.shape[0]):
assert isinstance(m, IdentityMap), 'Unrecognized data type, inherit from an IdentityMap or ComboMap!'
if ii > 0 and not (self.shape[1] == '*' or m.shape[0] == '*') and not self.shape[1] == m.shape[0]:
prev = self.maps[-1]
errArgs = (prev.__class__.__name__, prev.shape[0],
prev.shape[1], m.__class__.__name__, m.shape[0],
m.shape[1])
raise ValueError('Dimension mismatch in map[%s] (%s, %s) '
'and map[%s] (%s, %s).' % errArgs)
errArgs = (prev.__class__.__name__, prev.shape[0], prev.shape[1], m.__class__.__name__, m.shape[0], m.shape[1])
raise ValueError('Dimension mismatch in map[%s] (%s, %s) and map[%s] (%s, %s).' % errArgs)
if isinstance(m, ComboMap):
self.maps += m.maps
@@ -203,13 +164,8 @@ class ComboMap(IdentityMap):
m = map_i * m
return m
def deriv(self, m, v=None):
if v is not None:
deriv = v
else:
deriv = 1
def deriv(self, m):
deriv = 1
mi = m
for map_i in reversed(self.maps):
deriv = map_i.deriv(mi) * deriv
@@ -256,10 +212,11 @@ class ExpMap(IdentityMap):
"""
return np.log(Utils.mkvc(D))
def deriv(self, m, v=None):
def deriv(self, m):
"""
:param numpy.array m: model
:rtype: scipy.sparse.csr_matrix
:rtype: scipy.csr_matrix
:return: derivative of transformed model
The *transform* changes the model into the physical property.
@@ -279,11 +236,7 @@ class ExpMap(IdentityMap):
\\frac{\partial \exp{m}}{\partial m} = \\text{sdiag}(\exp{m})
"""
deriv = Utils.sdiag(np.exp(Utils.mkvc(m)))
if v is not None:
return deriv * v
return deriv
return Utils.sdiag(np.exp(Utils.mkvc(m)))
class ReciprocalMap(IdentityMap):
"""
@@ -300,12 +253,10 @@ class ReciprocalMap(IdentityMap):
def inverse(self, D):
return 1.0 / Utils.mkvc(m)
def deriv(self, m, v=None):
def deriv(self, m):
# TODO: if this is a tensor, you might have a problem.
deriv = Utils.sdiag( - Utils.mkvc(m)**(-2) )
if v is not None:
return deriv * v
return deriv
return Utils.sdiag( - Utils.mkvc(m)**(-2) )
class LogMap(IdentityMap):
@@ -314,13 +265,13 @@ class LogMap(IdentityMap):
If \\(p\\) is the physical property and \\(m\\) is the model, then
.. math::
..math::
p = \\log(m)
and
.. math::
..math::
m = \\exp(p)
@@ -335,20 +286,17 @@ class LogMap(IdentityMap):
def _transform(self, m):
return np.log(Utils.mkvc(m))
def deriv(self, m, v=None):
def deriv(self, m):
mod = Utils.mkvc(m)
deriv = np.zeros(mod.shape)
tol = 1e-16 # zero
ind = np.greater_equal(np.abs(mod), tol)
ind = np.greater_equal(np.abs(mod),tol)
deriv[ind] = 1.0/mod[ind]
if v is not None:
return Utils.sdiag(deriv)*v
return Utils.sdiag(deriv)
def inverse(self, m):
return np.exp(Utils.mkvc(m))
class SurjectFull(IdentityMap):
"""
SurjectFull
@@ -357,8 +305,8 @@ class SurjectFull(IdentityMap):
full model space.
"""
def __init__(self, mesh, **kwargs):
IdentityMap.__init__(self, mesh, **kwargs)
def __init__(self,mesh,**kwargs):
IdentityMap.__init__(self, mesh,**kwargs)
@property
def nP(self):
@@ -370,28 +318,22 @@ class SurjectFull(IdentityMap):
:rtype: numpy.array
:return: transformed model
"""
return np.ones(self.mesh.nC) * m
return np.ones(self.mesh.nC)*m
def deriv(self, m, v=None):
def deriv(self, m):
"""
:param numpy.array m: model
:rtype: numpy.array
:return: derivative of transformed model
"""
deriv = np.ones([self.mesh.nC,1])
if v is not None:
return deriv * v
return deriv
return np.ones([self.mesh.nC,1])
class FullMap(SurjectFull):
"""FullMap is depreciated. Use SurjectVertical1DMap instead.
"""
def __init__(self, mesh, **kwargs):
def __init__(self,mesh,**kwargs):
warnings.warn(
"`FullMap` is deprecated and will be removed in future versions. Use `SurjectFull` instead",
FutureWarning)
SurjectFull.__init__(self, mesh, **kwargs)
SurjectFull.__init__(self,mesh,**kwargs)
class SurjectVertical1D(IdentityMap):
"""SurjectVertical1DMap
@@ -421,10 +363,10 @@ class SurjectVertical1D(IdentityMap):
repNum = self.mesh.vnC[:self.mesh.dim-1].prod()
return Utils.mkvc(m).repeat(repNum)
def deriv(self, m, v=None):
def deriv(self, m):
"""
:param numpy.array m: model
:rtype: scipy.sparse.csr_matrix
:rtype: scipy.csr_matrix
:return: derivative of transformed model
"""
repNum = self.mesh.vnC[:self.mesh.dim-1].prod()
@@ -432,22 +374,14 @@ class SurjectVertical1D(IdentityMap):
(np.ones(repNum),
(range(repNum), np.zeros(repNum))
), shape=(repNum, 1))
deriv = sp.kron(sp.identity(self.nP), repVec)
if v is not None:
return deriv * v
return deriv
return sp.kron(sp.identity(self.nP), repVec)
class Vertical1DMap(SurjectVertical1D):
"""
Vertical1DMap is depreciated. Use SurjectVertical1D instead.
"""
def __init__(self, mesh, **kwargs):
def __init__(self,mesh,**kwargs):
warnings.warn(
"`Vertical1DMap` is deprecated and will be removed in future versions. Use `SurjectVertical1D` instead",
FutureWarning)
SurjectVertical1D.__init__(self, mesh, **kwargs)
SurjectVertical1D.__init__(self,mesh,**kwargs)
class Surject2Dto3D(IdentityMap):
"""Map2Dto3D
@@ -456,12 +390,12 @@ class Surject2Dto3D(IdentityMap):
3D model space.
"""
normal = 'Y' #: The normal
normal = 'Y' #: The normal
def __init__(self, mesh, **kwargs):
assert mesh.dim == 3, 'Only works for a 3D Mesh'
IdentityMap.__init__(self, mesh, **kwargs)
assert self.normal in ['X', 'Y', 'Z'], 'For now, only "Y" normal is supported'
assert self.normal in ['X','Y','Z'], 'For now, only "Y" normal is supported'
@property
def nP(self):
@@ -490,10 +424,10 @@ class Surject2Dto3D(IdentityMap):
elif self.normal == 'X':
return Utils.mkvc(m.reshape(self.mesh.vnC[[1,2]], order='F')[np.newaxis,:,:].repeat(self.mesh.nCx,axis=0))
def deriv(self, m, v=None):
def deriv(self, m):
"""
:param numpy.array m: model
:rtype: scipy.sparse.csr_matrix
:rtype: scipy.csr_matrix
:return: derivative of transformed model
"""
inds = self * np.arange(self.nP)
@@ -502,25 +436,19 @@ class Surject2Dto3D(IdentityMap):
(np.ones(nC),
(range(nC), inds)
), shape=(nC, nP))
if v is not None:
return P * v
return P
class Map2Dto3D(Surject2Dto3D):
"""Map2Dto3D is depreciated. Use Surject2Dto3D instead
"""
def __init__(self, mesh, **kwargs):
def __init__(self,mesh,**kwargs):
warnings.warn(
"`Map2Dto3D` is deprecated and will be removed in future versions. Use `Surject2Dto3D` instead",
FutureWarning)
Surject2Dto3D.__init__(self, mesh, **kwargs)
Surject2Dto3D.__init__(self,mesh,**kwargs)
class Mesh2Mesh(IdentityMap):
"""
Takes a model on one mesh are translates it to another mesh.
"""
def __init__(self, meshes, **kwargs):
@@ -533,7 +461,7 @@ class Mesh2Mesh(IdentityMap):
self.mesh = meshes[0]
self.mesh2 = meshes[1]
self.P = self.mesh2.getInterpolationMat(self.mesh.gridCC, 'CC', zerosOutside=True)
self.P = self.mesh2.getInterpolationMat(self.mesh.gridCC,'CC',zerosOutside=True)
@property
def shape(self):
@@ -544,13 +472,9 @@ class Mesh2Mesh(IdentityMap):
def nP(self):
"""Number of parameters in the model."""
return self.mesh2.nC
def _transform(self, m):
return self.P * m
def deriv(self, m, v=None):
if v is not None:
return self.P * v
return self.P*m
def deriv(self, m):
return self.P
@@ -560,9 +484,9 @@ class InjectActiveCells(IdentityMap):
"""
indActive = None #: Active Cells
valInactive = None #: Values of inactive Cells
nC = None #: Number of cells in the full model
indActive = None #: Active Cells
valInactive = None #: Values of inactive Cells
nC = None #: Number of cells in the full model
def __init__(self, mesh, indActive, valInactive, nC=None):
self.mesh = mesh
@@ -570,7 +494,7 @@ class InjectActiveCells(IdentityMap):
self.nC = nC or mesh.nC
if indActive.dtype is not bool:
z = np.zeros(self.nC, dtype=bool)
z = np.zeros(self.nC,dtype=bool)
z[indActive] = True
indActive = z
self.indActive = indActive
@@ -578,15 +502,11 @@ class InjectActiveCells(IdentityMap):
if Utils.isScalar(valInactive):
self.valInactive = np.ones(self.nC)*float(valInactive)
else:
self.valInactive = np.ones(self.nC)
self.valInactive[self.indInactive] = valInactive.copy()
self.valInactive = valInactive.copy()
self.valInactive[self.indActive] = 0
inds = np.nonzero(self.indActive)[0]
self.P = sp.csr_matrix((np.ones(inds.size), (inds, range(inds.size))),
shape=(self.nC, self.nP)
)
self.P = sp.csr_matrix((np.ones(inds.size),(inds, range(inds.size))), shape=(self.nC, self.nP))
@property
def shape(self):
@@ -598,25 +518,18 @@ class InjectActiveCells(IdentityMap):
return self.indActive.sum()
def _transform(self, m):
return self.P * m + self.valInactive
return self.P*m + self.valInactive
def inverse(self, D):
return self.P.T*D
def deriv(self, m, v=None):
if v is not None:
return self.P * v
def deriv(self, m):
return self.P
class ActiveCells(InjectActiveCells):
"""ActiveCells is depreciated. Use InjectActiveCells instead.
"""
def __init__(self, mesh, indActive, valInactive, nC=None):
warnings.warn(
"`ActiveCells` is deprecated and will be removed in future "
"versions. Use `InjectActiveCells` instead",
"`ActiveCells` is deprecated and will be removed in future versions. Use `InjectActiveCells` instead",
FutureWarning)
InjectActiveCells.__init__(self, mesh, indActive, valInactive, nC)
@@ -624,10 +537,11 @@ class ActiveCells(InjectActiveCells):
class Weighting(IdentityMap):
"""
Model weight parameters.
"""
weights = None #: Active Cells
nC = None #: Number of cells in the full model
weights = None #: Active Cells
nC = None #: Number of cells in the full model
def __init__(self, mesh, weights=None, nC=None):
self.mesh = mesh
@@ -657,9 +571,7 @@ class Weighting(IdentityMap):
Pinv = Utils.sdiag(self.weights**(-1.))
return Pinv*D
def deriv(self, m, v=None):
if v is not None:
return self.P * v
def deriv(self, m):
return self.P
@@ -687,32 +599,26 @@ class ComplexMap(IdentityMap):
nC = self.mesh.nC
return m[:nC] + m[nC:]*1j
def deriv(self, m, v=None):
def deriv(self, m):
nC = self.nP/2
shp = (nC, nC*2)
def fwd(v):
return v[:nC] + v[nC:]*1j
def adj(v):
return np.r_[v.real, v.imag]
if v is not None:
return LinearOperator(shp, matvec=fwd, rmatvec=adj) * v
return LinearOperator(shp, matvec=fwd, rmatvec=adj)
return np.r_[v.real,v.imag]
return LinearOperator(shp,matvec=fwd,rmatvec=adj)
inverse = deriv
class ParametricCircleMap(IdentityMap):
"""ParametricCircleMap
class CircleMap(IdentityMap):
"""CircleMap
Parameterize the model space using a circle in a wholespace.
..math::
\sigma(m) = \sigma_1 + (\sigma_2 - \sigma_1)\left(
\\arctan\left(100*\sqrt{(\\vec{x}-x_0)^2 + (\\vec{y}-y_0)}-r
\\right) \pi^{-1} + 0.5\\right)
\sigma(m) = \sigma_1 + (\sigma_2 - \sigma_1)\left(\\arctan\left(100*\sqrt{(\\vec{x}-x_0)^2 + (\\vec{y}-y_0)}-r\\right) \pi^{-1} + 0.5\\right)
Define the model as:
@@ -721,69 +627,46 @@ class ParametricCircleMap(IdentityMap):
m = [\sigma_1, \sigma_2, x_0, y_0, r]
"""
def __init__(self, mesh, logSigma=True):
assert mesh.dim == 2, "Working for a 2D mesh only right now. But it isn't that hard to change.. :)"
IdentityMap.__init__(self, mesh)
self.logSigma = logSigma
slope = 1e-1
def __init__(self, mesh, logSigma=True):
assert mesh.dim == 2, "Working for a 2D mesh only right now. "
"But it isn't that hard to change.. :)"
IdentityMap.__init__(self, mesh)
# TODO: this should be done through a composition with and ExpMap
self.logSigma = logSigma
@property
def nP(self):
return 5
def _transform(self, m):
a = self.slope
sig1, sig2, x, y, r = m[0], m[1], m[2], m[3], m[4]
sig1,sig2,x,y,r = m[0],m[1],m[2],m[3],m[4]
if self.logSigma:
sig1, sig2 = np.exp(sig1), np.exp(sig2)
X = self.mesh.gridCC[:, 0]
Y = self.mesh.gridCC[:, 1]
return sig1 + (sig2 - sig1)*(np.arctan(a*(np.sqrt((X-x)**2 +
(Y-y)**2) - r))/np.pi + 0.5)
X = self.mesh.gridCC[:,0]
Y = self.mesh.gridCC[:,1]
return sig1 + (sig2 - sig1)*(np.arctan(a*(np.sqrt((X-x)**2 + (Y-y)**2) - r))/np.pi + 0.5)
def deriv(self, m, v=None):
def deriv(self, m):
a = self.slope
sig1, sig2, x, y, r = m[0], m[1], m[2], m[3], m[4]
sig1,sig2,x,y,r = m[0],m[1],m[2],m[3],m[4]
if self.logSigma:
sig1, sig2 = np.exp(sig1), np.exp(sig2)
X = self.mesh.gridCC[:, 0]
Y = self.mesh.gridCC[:, 1]
X = self.mesh.gridCC[:,0]
Y = self.mesh.gridCC[:,1]
if self.logSigma:
g1 = -(np.arctan(a*(-r + np.sqrt((X - x)**2 + (Y - y)**2)))/np.pi +
0.5)*sig1 + sig1
g2 = (np.arctan(a*(-r + np.sqrt((X - x)**2 + (Y - y)**2)))/np.pi +
0.5)*sig2
g1 = -(np.arctan(a*(-r + np.sqrt((X - x)**2 + (Y - y)**2)))/np.pi + 0.5)*sig1 + sig1
g2 = (np.arctan(a*(-r + np.sqrt((X - x)**2 + (Y - y)**2)))/np.pi + 0.5)*sig2
else:
g1 = -(np.arctan(a*(-r + np.sqrt((X - x)**2 + (Y - y)**2)))/np.pi +
0.5) + 1.0
g2 = (np.arctan(a*(-r + np.sqrt((X - x)**2 + (Y - y)**2)))/np.pi +
0.5)
g1 = -(np.arctan(a*(-r + np.sqrt((X - x)**2 + (Y - y)**2)))/np.pi + 0.5) + 1.0
g2 = (np.arctan(a*(-r + np.sqrt((X - x)**2 + (Y - y)**2)))/np.pi + 0.5)
g3 = a*(-X + x)*(-sig1 + sig2)/(np.pi*(a**2*(-r + np.sqrt((X - x)**2 + (Y - y)**2))**2 + 1)*np.sqrt((X - x)**2 + (Y - y)**2))
g4 = a*(-Y + y)*(-sig1 + sig2)/(np.pi*(a**2*(-r + np.sqrt((X - x)**2 + (Y - y)**2))**2 + 1)*np.sqrt((X - x)**2 + (Y - y)**2))
g5 = -a*(-sig1 + sig2)/(np.pi*(a**2*(-r + np.sqrt((X - x)**2 + (Y - y)**2))**2 + 1))
if v is not None:
return sp.csr_matrix(np.c_[g1, g2, g3, g4, g5]) * v
return sp.csr_matrix(np.c_[g1, g2, g3, g4, g5])
return sp.csr_matrix(np.c_[g1,g2,g3,g4,g5])
class CircleMap(ParametricCircleMap):
"""CircleMap is depreciated. Use ParametricCircleMap instead.
"""
def __init__(self, mesh, logSigma=True):
warnings.warn(
"`CircleMap` is deprecated and will be removed in future "
"versions. Use `ParametricCircleMap` instead",
FutureWarning)
ParametricCircleMap.__init__(self, mesh, logSigma)
class ParametricPolyMap(IdentityMap):
class PolyMap(IdentityMap):
"""PolyMap
@@ -802,8 +685,7 @@ class ParametricPolyMap(IdentityMap):
Can take in an actInd vector to account for topography.
"""
def __init__(self, mesh, order, logSigma=True, normal='X', actInd=None):
def __init__(self, mesh, order, logSigma=True, normal='X', actInd = None):
IdentityMap.__init__(self, mesh)
self.logSigma = logSigma
self.order = order
@@ -828,88 +710,78 @@ class ParametricPolyMap(IdentityMap):
if np.isscalar(self.order):
nP = self.order+3
else:
nP = (self.order[0]+1)*(self.order[1]+1)+2
nP =(self.order[0]+1)*(self.order[1]+1)+2
return nP
def _transform(self, m):
# Set model parameters
alpha = self.slope
sig1, sig2 = m[0], m[1]
sig1,sig2 = m[0],m[1]
c = m[2:]
if self.logSigma:
sig1, sig2 = np.exp(sig1), np.exp(sig2)
# 2D
#2D
if self.mesh.dim == 2:
X = self.mesh.gridCC[self.actInd, 0]
Y = self.mesh.gridCC[self.actInd, 1]
X = self.mesh.gridCC[self.actInd,0]
Y = self.mesh.gridCC[self.actInd,1]
if self.normal =='X':
f = polynomial.polyval(Y, c) - X
elif self.normal =='Y':
f = polynomial.polyval(X, c) - Y
else:
raise(Exception("Input for normal = X or Y or Z"))
# 3D
#3D
elif self.mesh.dim == 3:
X = self.mesh.gridCC[self.actInd, 0]
Y = self.mesh.gridCC[self.actInd, 1]
Z = self.mesh.gridCC[self.actInd, 2]
if self.normal == 'X':
f = (polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,
self.order[1]+1))) - X)
elif self.normal == 'Y':
f = (polynomial.polyval2d(X, Z, c.reshape((self.order[0]+1,
self.order[1]+1))) - Y)
elif self.normal == 'Z':
f = (polynomial.polyval2d(X, Y, c.reshape((self.order[0]+1,
self.order[1]+1))) - Z)
X = self.mesh.gridCC[self.actInd,0]
Y = self.mesh.gridCC[self.actInd,1]
Z = self.mesh.gridCC[self.actInd,2]
if self.normal =='X':
f = polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - X
elif self.normal =='Y':
f = polynomial.polyval2d(X, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - Y
elif self.normal =='Z':
f = polynomial.polyval2d(X, Y, c.reshape((self.order[0]+1,self.order[1]+1))) - Z
else:
raise(Exception("Input for normal = X or Y or Z"))
else:
raise(Exception("Only supports 2D"))
return sig1+(sig2-sig1)*(np.arctan(alpha*f)/np.pi+0.5)
def deriv(self, m, v=None):
def deriv(self, m):
alpha = self.slope
sig1, sig2, c = m[0], m[1], m[2:]
sig1,sig2, c = m[0],m[1],m[2:]
if self.logSigma:
sig1, sig2 = np.exp(sig1), np.exp(sig2)
# 2D
#2D
if self.mesh.dim == 2:
X = self.mesh.gridCC[self.actInd, 0]
Y = self.mesh.gridCC[self.actInd, 1]
X = self.mesh.gridCC[self.actInd,0]
Y = self.mesh.gridCC[self.actInd,1]
if self.normal == 'X':
if self.normal =='X':
f = polynomial.polyval(Y, c) - X
V = polynomial.polyvander(Y, len(c)-1)
elif self.normal == 'Y':
elif self.normal =='Y':
f = polynomial.polyval(X, c) - Y
V = polynomial.polyvander(X, len(c)-1)
else:
raise(Exception("Input for normal = X or Y or Z"))
# 3D
#3D
elif self.mesh.dim == 3:
X = self.mesh.gridCC[self.actInd, 0]
Y = self.mesh.gridCC[self.actInd, 1]
Z = self.mesh.gridCC[self.actInd, 2]
X = self.mesh.gridCC[self.actInd,0]
Y = self.mesh.gridCC[self.actInd,1]
Z = self.mesh.gridCC[self.actInd,2]
if self.normal == 'X':
f = (polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,
self.order[1]+1))) - X)
if self.normal =='X':
f = polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - X
V = polynomial.polyvander2d(Y, Z, self.order)
elif self.normal == 'Y':
f = (polynomial.polyval2d(X, Z, c.reshape((self.order[0]+1,
self.order[1]+1))) - Y)
elif self.normal =='Y':
f = polynomial.polyval2d(X, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - Y
V = polynomial.polyvander2d(X, Z, self.order)
elif self.normal == 'Z':
f = (polynomial.polyval2d(X, Y, c.reshape((self.order[0]+1,
self.order[1]+1))) - Z)
elif self.normal =='Z':
f = polynomial.polyval2d(X, Y, c.reshape((self.order[0]+1,self.order[1]+1))) - Z
V = polynomial.polyvander2d(X, Y, self.order)
else:
raise(Exception("Input for normal = X or Y or Z"))
@@ -923,17 +795,13 @@ class ParametricPolyMap(IdentityMap):
g3 = Utils.sdiag(alpha*(sig2-sig1)/(1.+(alpha*f)**2)/np.pi)*V
if v is not None:
return sp.csr_matrix(np.c_[g1, g2, g3]) * v
return sp.csr_matrix(np.c_[g1, g2, g3])
return sp.csr_matrix(np.c_[g1,g2,g3])
class ParametricSplineMap(IdentityMap):
class SplineMap(IdentityMap):
"""SplineMap
Parameterize the boundary of two geological units using
a spline interpolation
Parameterize the boundary of two geological units using a spline interpolation
..math::
@@ -946,10 +814,7 @@ class ParametricSplineMap(IdentityMap):
m = [\sigma_1, \sigma_2, y]
"""
slope = 1e4
def __init__(self, mesh, pts, ptsv=None, order=3, logSigma=True, normal='X'):
def __init__(self, mesh, pts, ptsv=None,order=3, logSigma=True, normal='X'):
IdentityMap.__init__(self, mesh)
self.logSigma = logSigma
self.order = order
@@ -959,6 +824,7 @@ class ParametricSplineMap(IdentityMap):
self.ptsv = ptsv
self.spl = None
slope = 1e4
@property
def nP(self):
if self.mesh.dim == 2:
@@ -971,18 +837,18 @@ class ParametricSplineMap(IdentityMap):
def _transform(self, m):
# Set model parameters
alpha = self.slope
sig1, sig2 = m[0], m[1]
sig1,sig2 = m[0],m[1]
c = m[2:]
if self.logSigma:
sig1, sig2 = np.exp(sig1), np.exp(sig2)
# 2D
#2D
if self.mesh.dim == 2:
X = self.mesh.gridCC[:, 0]
Y = self.mesh.gridCC[:, 1]
X = self.mesh.gridCC[:,0]
Y = self.mesh.gridCC[:,1]
self.spl = UnivariateSpline(self.pts, c, k=self.order, s=0)
if self.normal == 'X':
if self.normal =='X':
f = self.spl(Y) - X
elif self.normal == 'Y':
elif self.normal =='Y':
f = self.spl(X) - Y
else:
raise(Exception("Input for normal = X or Y or Z"))
@@ -994,18 +860,18 @@ class ParametricSplineMap(IdentityMap):
# Using 2D interpolation is possible
elif self.mesh.dim == 3:
X = self.mesh.gridCC[:, 0]
Y = self.mesh.gridCC[:, 1]
Z = self.mesh.gridCC[:, 2]
X = self.mesh.gridCC[:,0]
Y = self.mesh.gridCC[:,1]
Z = self.mesh.gridCC[:,2]
npts = np.size(self.pts)
if np.mod(c.size, 2):
raise(Exception("Put even points!"))
self.spl = {"splb": UnivariateSpline(self.pts, c[:npts], k=self.order, s=0),
"splt": UnivariateSpline(self.pts, c[npts:], k=self.order, s=0)}
self.spl = {"splb":UnivariateSpline(self.pts, c[:npts], k=self.order, s=0),
"splt":UnivariateSpline(self.pts, c[npts:], k=self.order, s=0)}
if self.normal == 'X':
if self.normal =='X':
zb = self.ptsv[0]
zt = self.ptsv[1]
flines = (self.spl["splt"](Y)-self.spl["splb"](Y))*(Z-zb)/(zt-zb) + self.spl["splb"](Y)
@@ -1017,29 +883,30 @@ class ParametricSplineMap(IdentityMap):
else:
raise(Exception("Only supports 2D and 3D"))
return sig1+(sig2-sig1)*(np.arctan(alpha*f)/np.pi+0.5)
def deriv(self, m, v=None):
def deriv(self, m):
alpha = self.slope
sig1, sig2, c = m[0], m[1], m[2:]
sig1,sig2, c = m[0],m[1],m[2:]
if self.logSigma:
sig1, sig2 = np.exp(sig1), np.exp(sig2)
#2D
if self.mesh.dim == 2:
X = self.mesh.gridCC[:, 0]
Y = self.mesh.gridCC[:, 1]
X = self.mesh.gridCC[:,0]
Y = self.mesh.gridCC[:,1]
if self.normal == 'X':
if self.normal =='X':
f = self.spl(Y) - X
elif self.normal == 'Y':
elif self.normal =='Y':
f = self.spl(X) - Y
else:
raise(Exception("Input for normal = X or Y or Z"))
#3D
elif self.mesh.dim == 3:
X = self.mesh.gridCC[:, 0]
Y = self.mesh.gridCC[:, 1]
Z = self.mesh.gridCC[:, 2]
X = self.mesh.gridCC[:,0]
Y = self.mesh.gridCC[:,1]
Z = self.mesh.gridCC[:,2]
if self.normal =='X':
zb = self.ptsv[0]
zt = self.ptsv[1]
@@ -1057,9 +924,10 @@ class ParametricSplineMap(IdentityMap):
g1 = -(np.arctan(alpha*f)/np.pi + 0.5) + 1.0
g2 = (np.arctan(alpha*f)/np.pi + 0.5)
if self.mesh.dim == 2:
if self.mesh.dim ==2:
g3 = np.zeros((self.mesh.nC, self.npts))
if self.normal == 'Y':
if self.normal =='Y':
# Here we use perturbation to compute sensitivity
# TODO: bit more generalization of this ...
# Modfications for X and Z directions ...
@@ -1074,11 +942,11 @@ class ParametricSplineMap(IdentityMap):
spla = UnivariateSpline(self.pts, ca, k=self.order, s=0)
splb = UnivariateSpline(self.pts, cb, k=self.order, s=0)
fderiv = (spla(X)-splb(X))/(2*dy)
g3[:, i] = Utils.sdiag(alpha*(sig2-sig1)/(1.+(alpha*f)**2)/np.pi)*fderiv
g3[:,i] = Utils.sdiag(alpha*(sig2-sig1)/(1.+(alpha*f)**2)/np.pi)*fderiv
elif self.mesh.dim == 3:
elif self.mesh.dim==3:
g3 = np.zeros((self.mesh.nC, self.npts*2))
if self.normal == 'X':
if self.normal =='X':
# Here we use perturbation to compute sensitivity
for i in range(self.npts*2):
ctemp = c[i]
@@ -1088,40 +956,26 @@ class ParametricSplineMap(IdentityMap):
dy = self.mesh.hy[ind]*1.5
ca[i] = ctemp+dy
cb[i] = ctemp-dy
# treat bottom boundary
if i < self.npts:
#treat bottom boundary
if i< self.npts:
splba = UnivariateSpline(self.pts, ca[:self.npts], k=self.order, s=0)
splbb = UnivariateSpline(self.pts, cb[:self.npts], k=self.order, s=0)
flinesa = (self.spl["splt"](Y)-splba(Y))*(Z-zb)/(zt-zb) + splba(Y) - X
flinesb = (self.spl["splt"](Y)-splbb(Y))*(Z-zb)/(zt-zb) + splbb(Y) - X
# treat top boundary
#treat top boundary
else:
splta = UnivariateSpline(self.pts, ca[self.npts:], k=self.order, s=0)
spltb = UnivariateSpline(self.pts, ca[self.npts:], k=self.order, s=0)
flinesa = (self.spl["splt"](Y)-splta(Y))*(Z-zb)/(zt-zb) + splta(Y) - X
flinesb = (self.spl["splt"](Y)-spltb(Y))*(Z-zb)/(zt-zb) + spltb(Y) - X
fderiv = (flinesa-flinesb)/(2*dy)
g3[:, i] = Utils.sdiag(alpha*(sig2-sig1)/(1.+(alpha*f)**2)/np.pi)*fderiv
g3[:,i] = Utils.sdiag(alpha*(sig2-sig1)/(1.+(alpha*f)**2)/np.pi)*fderiv
else :
raise(Exception("Not Implemented for Y and Z, your turn :)"))
return sp.csr_matrix(np.c_[g1,g2,g3])
if v is not None:
return sp.csr_matrix(np.c_[g1, g2, g3]) * v
return sp.csr_matrix(np.c_[g1, g2, g3])
class SplineMap(ParametricSplineMap):
"""SplineMap is depreciated. Use ParametricSplineMap instead.
"""
def __init__(self, mesh, pts, ptsv=None, order=3, logSigma=True,
normal='X'):
warnings.warn(
"`SplineMap` is deprecated and will be removed in future "
"versions. Use `ParametricSplineMap` instead",
FutureWarning)
ParametricSplineMap.__init__(self, mesh, pts, ptsv, order, logSigma,
normal)
+24 -26
View File
@@ -7,8 +7,8 @@ class BaseMesh(object):
BaseMesh does all the counting you don't want to do.
BaseMesh should be inherited by meshes with a regular structure.
:param numpy.array n: (or list) number of cells in each direction (dim, )
:param numpy.array x0: (or list) Origin of the mesh (dim, )
:param numpy.array,list n: number of cells in each direction (dim, )
:param numpy.array,list x0: Origin of the mesh (dim, )
"""
@@ -34,8 +34,8 @@ class BaseMesh(object):
"""
Origin of the mesh
:rtype: numpy.array
:return: x0, (dim, )
:rtype: numpy.array (dim, )
:return: x0
"""
return self._x0
@@ -116,8 +116,8 @@ class BaseMesh(object):
"""
Total number of edges in each direction
:rtype: numpy.array
:return: [nEx, nEy, nEz], (dim, )
:rtype: numpy.array (dim, )
:return: [nEx, nEy, nEz]
.. plot::
:include-source:
@@ -173,8 +173,8 @@ class BaseMesh(object):
"""
Total number of faces in each direction
:rtype: numpy.array
:return: [nFx, nFy, nFz], (dim, )
:rtype: numpy.array (dim, )
:return: [nFx, nFy, nFz]
.. plot::
:include-source:
@@ -200,8 +200,8 @@ class BaseMesh(object):
"""
Face Normals
:rtype: numpy.array
:return: normals, (sum(nF), dim)
:rtype: numpy.array (sum(nF), dim)
:return: normals
"""
if self.dim == 2:
nX = np.c_[np.ones(self.nFx), np.zeros(self.nFx)]
@@ -218,8 +218,8 @@ class BaseMesh(object):
"""
Edge Tangents
:rtype: numpy.array
:return: normals, (sum(nE), dim)
:rtype: numpy.array (sum(nE), dim)
:return: normals
"""
if self.dim == 2:
tX = np.c_[np.ones(self.nEx), np.zeros(self.nEx)]
@@ -236,9 +236,8 @@ class BaseMesh(object):
Given a vector, fV, in cartesian coordinates, this will project it onto the mesh using the normals
:param numpy.array fV: face vector with shape (nF, dim)
:rtype: numpy.array
:return: projected face vector, (nF, )
:rtype: numpy.array with shape (nF, )
:return: projected face vector
"""
assert isinstance(fV, np.ndarray), 'fV must be an ndarray'
assert len(fV.shape) == 2 and fV.shape[0] == self.nF and fV.shape[1] == self.dim, 'fV must be an ndarray of shape (nF x dim)'
@@ -249,9 +248,8 @@ class BaseMesh(object):
Given a vector, eV, in cartesian coordinates, this will project it onto the mesh using the tangents
:param numpy.array eV: edge vector with shape (nE, dim)
:rtype: numpy.array
:return: projected edge vector, (nE, )
:rtype: numpy.array with shape (nE, )
:return: projected edge vector
"""
assert isinstance(eV, np.ndarray), 'eV must be an ndarray'
assert len(eV.shape) == 2 and eV.shape[0] == self.nE and eV.shape[1] == self.dim, 'eV must be an ndarray of shape (nE x dim)'
@@ -297,7 +295,7 @@ class BaseRectangularMesh(BaseMesh):
"""
Total number of cells in each direction
:rtype: numpy.array
:rtype: numpy.array (dim, )
:return: [nCx, nCy, nCz]
"""
return np.array([x for x in [self.nCx, self.nCy, self.nCz] if not x is None])
@@ -337,7 +335,7 @@ class BaseRectangularMesh(BaseMesh):
"""
Total number of nodes in each direction
:rtype: numpy.array
:rtype: numpy.array (dim, )
:return: [nNx, nNy, nNz]
"""
return np.array([x for x in [self.nNx, self.nNy, self.nNz] if not x is None])
@@ -347,7 +345,7 @@ class BaseRectangularMesh(BaseMesh):
"""
Number of x-edges in each direction
:rtype: numpy.array
:rtype: numpy.array (dim, )
:return: vnEx
"""
return np.array([x for x in [self.nCx, self.nNy, self.nNz] if not x is None])
@@ -357,7 +355,7 @@ class BaseRectangularMesh(BaseMesh):
"""
Number of y-edges in each direction
:rtype: numpy.array
:rtype: numpy.array (dim, )
:return: vnEy or None if dim < 2
"""
return None if self.dim < 2 else np.array([x for x in [self.nNx, self.nCy, self.nNz] if not x is None])
@@ -367,7 +365,7 @@ class BaseRectangularMesh(BaseMesh):
"""
Number of z-edges in each direction
:rtype: numpy.array
:rtype: numpy.array (dim, )
:return: vnEz or None if dim < 3
"""
return None if self.dim < 3 else np.array([x for x in [self.nNx, self.nNy, self.nCz] if not x is None])
@@ -377,7 +375,7 @@ class BaseRectangularMesh(BaseMesh):
"""
Number of x-faces in each direction
:rtype: numpy.array
:rtype: numpy.array (dim, )
:return: vnFx
"""
return np.array([x for x in [self.nNx, self.nCy, self.nCz] if not x is None])
@@ -387,7 +385,7 @@ class BaseRectangularMesh(BaseMesh):
"""
Number of y-faces in each direction
:rtype: numpy.array
:rtype: numpy.array (dim, )
:return: vnFy or None if dim < 2
"""
return None if self.dim < 2 else np.array([x for x in [self.nCx, self.nNy, self.nCz] if not x is None])
@@ -397,7 +395,7 @@ class BaseRectangularMesh(BaseMesh):
"""
Number of z-faces in each direction
:rtype: numpy.array
:rtype: numpy.array (dim, )
:return: vnFz or None if dim < 3
"""
return None if self.dim < 3 else np.array([x for x in [self.nCx, self.nCy, self.nNz] if not x is None])
+6 -6
View File
@@ -68,8 +68,8 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
"""
Number of x-faces in each direction
:rtype: numpy.array
:return: vnFx, (dim, )
:rtype: numpy.array (dim, )
:return: vnFx
"""
return self.vnC
@@ -78,8 +78,8 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
"""
Number of y-edges in each direction
:rtype: numpy.array
:return: vnEy or None if dim < 2, (dim, )
:rtype: numpy.array (dim, )
:return: vnEy or None if dim < 2
"""
nNx = self.nNx if self.isSymmetric else self.nNx - 1
return np.r_[nNx, self.nCy, self.nNz]
@@ -89,8 +89,8 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
"""
Number of z-edges in each direction
:rtype: numpy.array
:return: vnEz or None if nCy > 1, (dim, )
:rtype: numpy.array (dim, )
:return: vnEz or None if nCy > 1
"""
if self.isSymmetric:
return np.r_[self.nNx, self.nNy, self.nCz]
+9 -8
View File
@@ -16,7 +16,7 @@ class InnerProducts(object):
:param bool invProp: inverts the material property
:param bool invMat: inverts the matrix
:param bool doFast: do a faster implementation if available.
:rtype: scipy.sparse.csr_matrix
:rtype: scipy.csr_matrix
:return: M, the inner product matrix (nF, nF)
"""
return self._getInnerProduct('F', prop=prop, invProp=invProp, invMat=invMat, doFast=doFast)
@@ -27,7 +27,7 @@ class InnerProducts(object):
:param bool invProp: inverts the material property
:param bool invMat: inverts the matrix
:param bool doFast: do a faster implementation if available.
:rtype: scipy.sparse.csr_matrix
:rtype: scipy.csr_matrix
:return: M, the inner product matrix (nE, nE)
"""
return self._getInnerProduct('E', prop=prop, invProp=invProp, invMat=invMat, doFast=doFast)
@@ -39,7 +39,7 @@ class InnerProducts(object):
:param bool invProp: inverts the material property
:param bool invMat: inverts the matrix
:param bool doFast: do a faster implementation if available.
:rtype: scipy.sparse.csr_matrix
:rtype: scipy.csr_matrix
:return: M, the inner product matrix (nE, nE)
"""
assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges"
@@ -115,12 +115,13 @@ class InnerProducts(object):
:param bool doFast: do a faster implementation if available.
:param bool invProp: inverts the material property
:param bool invMat: inverts the matrix
:rtype: function
:return: dMdmu(u), the derivative of the inner product matrix (u)
Given u, dMdmu returns (nF, nC*nA)
:param numpy.ndarray u: vector that multiplies dMdmu
:rtype: scipy.sparse.csr_matrix
:param np.ndarray u: vector that multiplies dMdmu
:rtype: scipy.csr_matrix
:return: dMdmu, the derivative of the inner product matrix for a certain u
"""
return self._getInnerProductDeriv(prop, 'F', doFast=doFast, invProp=invProp, invMat=invMat)
@@ -132,7 +133,7 @@ class InnerProducts(object):
:param bool doFast: do a faster implementation if available.
:param bool invProp: inverts the material property
:param bool invMat: inverts the matrix
:rtype: scipy.sparse.csr_matrix
:rtype: scipy.csr_matrix
:return: dMdm, the derivative of the inner product matrix (nE, nC*nA)
"""
return self._getInnerProductDeriv(prop, 'E', doFast=doFast, invProp=invProp, invMat=invMat)
@@ -144,7 +145,7 @@ class InnerProducts(object):
:param bool doFast: do a faster implementation if available.
:param bool invProp: inverts the material property
:param bool invMat: inverts the matrix
:rtype: scipy.sparse.csr_matrix
:rtype: scipy.csr_matrix
:return: dMdm, the derivative of the inner product matrix (nE, nC*nA)
"""
fast = None
@@ -168,7 +169,7 @@ class InnerProducts(object):
:param numpy.array v: vector to multiply (required in the general implementation)
:param list P: list of projection matrices
:param str projType: 'F' for faces 'E' for edges
:rtype: scipy.sparse.csr_matrix
:rtype: scipy.csr_matrix
:return: dMdm, the derivative of the inner product matrix (n, nC*nA)
"""
assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges"
+79 -23
View File
@@ -6,11 +6,13 @@ class TensorMeshIO(object):
@classmethod
def readUBC(TensorMesh, fileName):
"""
Read UBC GIF 3D tensor mesh and generate 3D TensorMesh in SimPEG.
Read UBC GIF 3DTensor mesh and generate 3D Tensor mesh in simpegTD
:param string fileName: path to the UBC GIF mesh file
:rtype: TensorMesh
:return: The tensor mesh for the fileName.
Input:
:param fileName, path to the UBC GIF mesh file
Output:
:param SimPEG TensorMesh object
"""
# Interal function to read cell size lines for the UBC mesh files.
@@ -46,9 +48,11 @@ class TensorMeshIO(object):
Read VTK Rectilinear (vtr xml file) and return SimPEG Tensor mesh and model
Input:
:param string fileName: path to the vtr model file to read
:rtype: tuple
:return: (TensorMesh, modelDictionary)
:param vtrFileName, path to the vtr model file to write to
Output:
:return SimPEG TensorMesh object
:return SimPEG model dictionary
"""
# Import
@@ -98,8 +102,9 @@ class TensorMeshIO(object):
Makes and saves a VTK rectilinear file (vtr) for a simpeg Tensor mesh and model.
Input:
:param string fileName: path to the output vtk file
:param dict models: dictionary of numpy.array - Name('s) and array('s). Match number of cells
:param str, path to the output vtk file
:param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK
:param models, dictionary of numpy.array - Name('s) and array('s). Match number of cells
"""
# Import
@@ -135,7 +140,6 @@ class TensorMeshIO(object):
vtkObj.GetCellData().AddArray(vtkDoubleArr)
# Set the active scalar
vtkObj.GetCellData().SetActiveScalars(models.keys()[0])
# vtkObj.Update()
# Check the extension of the fileName
ext = os.path.splitext(fileName)[1]
@@ -152,14 +156,61 @@ class TensorMeshIO(object):
vtrWriteFilter.SetFileName(fileName)
vtrWriteFilter.Update()
def _toVTRObj(mesh,models=None):
"""
Makes and saves a VTK rectilinear file (vtr) for a simpeg Tensor mesh and model.
Input:
:param str, path to the output vtk file
:param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK
:param models, dictionary of numpy.array - Name('s) and array('s). Match number of cells
"""
# Import
from vtk import vtkRectilinearGrid as rectGrid, VTK_VERSION
from vtk.util.numpy_support import numpy_to_vtk
# Deal with dimensionalities
if mesh.dim >= 1:
vX = mesh.vectorNx
xD = mesh.nNx
yD,zD = 1,1
vY, vZ = np.array([0,0])
if mesh.dim >= 2:
vY = mesh.vectorNy
yD = mesh.nNy
if mesh.dim == 3:
vZ = mesh.vectorNz
zD = mesh.nNz
# Use rectilinear VTK grid.
# Assign the spatial information.
vtkObj = rectGrid()
vtkObj.SetDimensions(xD,yD,zD)
vtkObj.SetXCoordinates(numpy_to_vtk(vX,deep=1))
vtkObj.SetYCoordinates(numpy_to_vtk(vY,deep=1))
vtkObj.SetZCoordinates(numpy_to_vtk(vZ,deep=1))
# Assign the model('s) to the object
if models is not None:
for item in models.iteritems():
# Convert numpy array
vtkDoubleArr = numpy_to_vtk(item[1],deep=1)
vtkDoubleArr.SetName(item[0])
vtkObj.GetCellData().AddArray(vtkDoubleArr)
# Set the active scalar
vtkObj.GetCellData().SetActiveScalars(models.keys()[0])
return vtkObj
def readModelUBC(mesh, fileName):
"""
Read UBC 3DTensor mesh model and generate 3D Tensor mesh model in simpeg
:param string fileName: path to the UBC GIF mesh file to read
:rtype: numpy.ndarray
:return: model with TensorMesh ordered
Input:
:param fileName, path to the UBC GIF mesh file to read
:param mesh, TensorMesh object, mesh that coresponds to the model
Output:
:return numpy array, model with TensorMesh ordered
"""
f = open(fileName, 'r')
model = np.array(map(float, f.readlines()))
@@ -175,7 +226,8 @@ class TensorMeshIO(object):
Writes a model associated with a SimPEG TensorMesh
to a UBC-GIF format model file.
:param string fileName: File to write to
:param str fileName: File to write to
:param simpeg.Mesh.TensorMesh mesh: The mesh
:param numpy.ndarray model: The model
"""
@@ -192,8 +244,8 @@ class TensorMeshIO(object):
"""
Writes a SimPEG TensorMesh to a UBC-GIF format mesh file.
:param string fileName: File to write to
:param dict models: A dictionary of the models
:param str fileName: File to write to
:param simpeg.Mesh.TensorMesh mesh: The mesh
"""
assert mesh.dim == 3
@@ -222,8 +274,9 @@ class TreeMeshIO(object):
"""
Write UBC ocTree mesh and model files from a simpeg ocTree mesh and model.
:param string fileName: File to write to
:param dict models: The models in a dictionary, where the keys is the name of the of the model file
:param str fileName: File to write to
:param simpeg.Mesh.TreeMesh mesh: The mesh
:param dictionary models: The models in a dictionary, where the keys is the name of the of the model file
"""
# Calculate information to write in the file.
@@ -276,9 +329,10 @@ class TreeMeshIO(object):
Input:
:param str meshFile: path to the UBC GIF OcTree mesh file to read
:rtype: SimPEG.Mesh.TreeMesh
:return: The octree mesh
Output:
:return SimPEG.Mesh.TreeMesh mesh: The octree mesh
:return list of ndarray's: models as a list of numpy array's
"""
## Read the file lines
@@ -324,9 +378,11 @@ class TreeMeshIO(object):
"""
Read UBC OcTree model and get vector
:param string fileName: path to the UBC GIF model file to read
:rtype: numpy.ndarray
:return: OcTree model
Input:
:param fileName, path to the UBC GIF model file to read
Output:
:return numpy array, OcTree model
"""
if type(fileName) is list:
+4 -4
View File
@@ -198,8 +198,8 @@ class BaseTensorMesh(BaseMesh):
Determines if a set of points are inside a mesh.
:param numpy.ndarray pts: Location of points to test
:rtype numpy.ndarray:
:return: inside, numpy array of booleans
:rtype numpy.ndarray
:return inside, numpy array of booleans
"""
pts = Utils.asArray_N_x_Dim(pts, self.dim)
@@ -221,7 +221,7 @@ class BaseTensorMesh(BaseMesh):
:param numpy.ndarray loc: Location of points to interpolate to
:param str locType: What to interpolate (see below)
:rtype: scipy.sparse.csr_matrix
:rtype: scipy.sparse.csr.csr_matrix
:return: M, the interpolation matrix
locType can be::
@@ -289,7 +289,7 @@ class BaseTensorMesh(BaseMesh):
:param bool returnP: returns the projection matrices
:param bool invProp: inverts the material property
:param bool invMat: inverts the matrix
:rtype: scipy.sparse.csr_matrix
:rtype: scipy.csr_matrix
:return: M, the inner product matrix (nF, nF)
"""
assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges"
+1 -1
View File
@@ -1875,7 +1875,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
:param numpy.ndarray locs: Location of points to interpolate to
:param str locType: What to interpolate (see below)
:rtype: scipy.sparse.csr_matrix
:rtype: scipy.sparse.csr.csr_matrix
:return: M, the interpolation matrix
locType can be::
@@ -4,18 +4,21 @@ import sys
from numpy.lib import recfunctions as recFunc
from SimPEG.EM.Utils import omega
##############
### Fields ###
##############
class BaseMTFields(Problem.Fields):
"""Field Storage for a MT survey."""
class BaseNSEMFields(Problem.Fields):
"""Field Storage for a NSEM survey."""
knownFields = {}
dtype = complex
class Fields1D_e(BaseMTFields):
###########
# 1D Fields
###########
class Fields1D_ePrimSec(BaseNSEMFields):
"""
Fields storage for the 1D MT solution.
Fields storage for the 1D NSEM solution.
"""
knownFields = {'e_1dSolution':'F'}
aliasFields = {
@@ -28,7 +31,119 @@ class Fields1D_e(BaseMTFields):
}
def __init__(self,mesh,survey,**kwargs):
BaseMTFields.__init__(self,mesh,survey,**kwargs)
BaseNSEMFields.__init__(self,mesh,survey,**kwargs)
def _ePrimary(self, eSolution, srcList):
ePrimary = np.zeros_like(eSolution)
for i, src in enumerate(srcList):
ep = src.ePrimary(self.survey.prob)
if ep is not None:
ePrimary[:,i] = ep[:,-1]
return ePrimary
def _eSecondary(self, eSolution, srcList):
return eSolution
def _e(self, eSolution, srcList):
return self._ePrimary(eSolution,srcList) + self._eSecondary(eSolution,srcList)
def _eDeriv_u(self, src, du_dm_v, adjoint = False):
return Utils.Identity()*du_dm_v
def _eDeriv_m(self, src, v, adjoint = False):
# assuming primary does not depend on the model
return Utils.Zero()
def _bPrimary(self, eSolution, srcList):
bPrimary = np.zeros([self.survey.mesh.nE,eSolution.shape[1]], dtype = complex)
for i, src in enumerate(srcList):
bp = src.bPrimary(self.survey.prob)
if bp is not None:
bPrimary[:,i] += bp[:,-1]
return bPrimary
def _bSecondary(self, eSolution, srcList):
C = self.mesh.nodalGrad
b = (C * eSolution)
for i, src in enumerate(srcList):
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:
# b[:,i] += 1./(1j*omega(src.freq)) * S_m
return b
def _b(self, eSolution, srcList):
return self._bPrimary(eSolution, srcList) + self._bSecondary(eSolution, srcList)
def _bSecondaryDeriv_u(self, src, v, adjoint = False):
C = self.mesh.nodalGrad
if adjoint:
return - 1./(1j*omega(src.freq)) * (C.T * v)
return - 1./(1j*omega(src.freq)) * (C * v)
def _bSecondaryDeriv_m(self, src, v, adjoint = False):
# Doesn't depend on m
# _, S_eDeriv = src.evalDeriv(self.survey.prob, adjoint)
# S_eDeriv = S_eDeriv(v)
# if S_eDeriv is not None:
# return 1./(1j * omega(src.freq)) * S_eDeriv
return None
def _bDeriv_u(self, src, v, adjoint=False):
# Primary does not depend on u
return self._bSecondaryDeriv_u(src, v, adjoint)
def _bDeriv_m(self, src, v, adjoint=False):
# Assuming the primary does not depend on the model
return self._bSecondaryDeriv_m(src, v, adjoint)
def _fDeriv_u(self, src, v, adjoint=False):
"""
Derivative of the fields object wrt u.
:param NSEMsrc src: NSEM source
:param numpy.ndarray v: random vector of f_sol.size
This function stacks the fields derivatives appropriately
return a vector of size (nreEle+nrbEle)
"""
de_du = v #Utils.spdiag(np.ones((self.nF,)))
db_du = self._bDeriv_u(src, v, adjoint)
# Return the stack
# This doesn't work...
return np.vstack((de_du,db_du))
def _fDeriv_m(self, src, v, adjoint=False):
"""
Derivative of the fields object wrt m.
This function stacks the fields derivatives appropriately
"""
return None
class Fields1D_eTotal(BaseNSEMFields):
"""
Fields storage for the 1D NSEM solution solved with for a total domain formulation.
Used in conjuction with Problem1D_eTotal.
"""
knownFields = {'e_1dSolution':'F'}
aliasFields = {
'e_1d' : ['e_1dSolution','F','_e'],
'e_1dPrimary' : ['e_1dSolution','F','_ePrimary'],
'e_1dSecondary' : ['e_1dSolution','F','_eSecondary'],
'b_1d' : ['e_1dSolution','E','_b'],
'b_1dPrimary' : ['e_1dSolution','E','_bPrimary'],
'b_1dSecondary' : ['e_1dSolution','E','_bSecondary']
}
def __init__(self,mesh,survey,**kwargs):
BaseNSEMFields.__init__(self,mesh,survey,**kwargs)
def _ePrimary(self, eSolution, srcList):
ePrimary = np.zeros_like(eSolution)
@@ -99,7 +214,7 @@ class Fields1D_e(BaseMTFields):
"""
Derivative of the fields object wrt u.
:param MTsrc src: MT source
:param NSEMsrc src: NSEM source
:param numpy.ndarray v: random vector of f_sol.size
This function stacks the fields derivatives appropriately
@@ -120,9 +235,18 @@ class Fields1D_e(BaseMTFields):
"""
return None
class Fields3D_e(BaseMTFields):
###########
# 2D Fields
###########
###########
# 3D Fields
###########
class Fields3D_ePrimSec(BaseNSEMFields):
"""
Fields storage for the 3D MT solution. Labels polarizations by px and py.
Fields storage for the 3D NSEM solution. Labels polarizations by px and py.
:param SimPEG object mesh: The solution mesh
:param SimPEG object survey: A survey object
@@ -147,7 +271,7 @@ class Fields3D_e(BaseMTFields):
}
def __init__(self,mesh,survey,**kwargs):
BaseMTFields.__init__(self,mesh,survey,**kwargs)
BaseNSEMFields.__init__(self,mesh,survey,**kwargs)
def _e_pxPrimary(self, e_pxSolution, srcList):
e_pxPrimary = np.zeros_like(e_pxSolution)
@@ -228,7 +352,7 @@ class Fields3D_e(BaseMTFields):
b = (C * e_pxSolution)
for i, src in enumerate(srcList):
b[:,i] *= - 1./(1j*omega(src.freq))
# There is no magnetic source in the MT problem
# There is no magnetic source in the NSEM problem
# S_m, _ = src.eval(self.survey.prob)
# if S_m is not None:
# b[:,i] += 1./(1j*omega(src.freq)) * S_m
@@ -239,7 +363,7 @@ class Fields3D_e(BaseMTFields):
b = (C * e_pySolution)
for i, src in enumerate(srcList):
b[:,i] *= - 1./(1j*omega(src.freq))
# There is no magnetic source in the MT problem
# There is no magnetic source in the NSEM problem
# S_m, _ = src.eval(self.survey.prob)
# if S_m is not None:
# b[:,i] += 1./(1j*omega(src.freq)) * S_m
@@ -302,7 +426,7 @@ class Fields3D_e(BaseMTFields):
"""
Derivative of the fields object wrt u.
:param MTsrc src: MT source
:param NSEMsrc src: NSEM source
:param numpy.ndarray v: random vector of f_sol.size
This function stacks the fields derivatives appropriately
@@ -319,7 +443,7 @@ class Fields3D_e(BaseMTFields):
"""
Derivative of the fields object wrt u.
:param MTsrc src: MT source
:param NSEMsrc src: NSEM source
:param numpy.ndarray v: random vector of f_sol.size
This function stacks the fields derivatives appropriately
+560
View File
@@ -0,0 +1,560 @@
from SimPEG.EM.Utils.EMUtils import omega, mu_0
from SimPEG import SolverLU as SimpegSolver, PropMaps, Utils, mkvc, sp, np
from SimPEG.EM.FDEM.ProblemFDEM import BaseFDEMProblem
from SurveyNSEM import Survey, Data
from FieldsNSEM import BaseNSEMFields, Fields1D_ePrimSec, Fields3D_ePrimSec
from SimPEG.NSEM.Utils.MT1Danalytic import getEHfields
import time, sys
class BaseNSEMProblem(BaseFDEMProblem):
"""
Base class for all Natural source problems.
"""
def __init__(self, mesh, **kwargs):
BaseFDEMProblem.__init__(self, mesh, **kwargs)
Utils.setKwargs(self, **kwargs)
# Set the default pairs of the problem
surveyPair = Survey
dataPair = Data
fieldsPair = BaseNSEMFields
# Set the solver
Solver = SimpegSolver
solverOpts = {}
verbose = False
# Notes:
# Use the forward and devs from BaseFDEMProblem
# Might need to add more stuff here.
## NEED to clean up the Jvec and Jtvec to use Zero and Identities for None components.
def Jvec(self, m, v, f=None):
"""
Function to calculate the data sensitivities dD/dm times a vector.
:param numpy.ndarray m (nC, 1) - conductive model
:param numpy.ndarray v (nC, 1) - random vector
:param NSEMfields object (optional) - NSEM fields object, if not given it is calculated
:rtype: NSEMdata object
:return: Data sensitivities wrt m
"""
# Calculate the fields
if f is None:
f= self.fields(m)
# Set current model
self.curModel = m
# Initiate the Jv object
Jv = self.dataPair(self.survey)
# Loop all the frequenies
for freq in self.survey.freqs:
dA_du = self.getA(freq) #
dA_duI = self.Solver(dA_du, **self.solverOpts)
for src in self.survey.getSrcByFreq(freq):
# We need fDeriv_m = df/du*du/dm + df/dm
# Construct du/dm, it requires a solve
# NOTE: need to account for the 2 polarizations in the derivatives.
u_src = f[src,:] # u should be a vector by definition. Need to fix this...
# dA_dm and dRHS_dm should be of size nE,2, so that we can multiply by dA_duI. The 2 columns are each of the polarizations.
dA_dm = self.getADeriv_m(freq, u_src, v) # Size: nE,2 (u_px,u_py) in the columns.
dRHS_dm = self.getRHSDeriv_m(freq, v) # Size: nE,2 (u_px,u_py) in the columns.
if dRHS_dm is None:
du_dm = dA_duI * ( -dA_dm )
else:
du_dm = dA_duI * ( -dA_dm + dRHS_dm )
# Calculate the projection derivatives
for rx in src.rxList:
# Get the projection derivative
# v should be of size 2*nE (for 2 polarizations)
PDeriv_u = lambda t: rx.evalDeriv(src, self.mesh, f, t) # wrt u, we don't have have PDeriv wrt m
Jv[src, rx] = PDeriv_u(mkvc(du_dm))
dA_duI.clean()
# Return the vectorized sensitivities
return mkvc(Jv)
def Jtvec(self, m, v, f=None):
"""
Function to calculate the transpose of the data sensitivities (dD/dm)^T times a vector.
:param numpy.ndarray m (nC, 1) - conductive model
:param numpy.ndarray v (nD, 1) - vector
:param NSEMfields object f (optional) - NSEM fields object, if not given it is calculated
:rtype: NSEMdata object
:return: Data sensitivities wrt m
"""
if f is None:
f = self.fields(m)
self.curModel = m
# Ensure v is a data object.
if not isinstance(v, self.dataPair):
v = self.dataPair(self.survey, v)
Jtv = np.zeros(m.size)
for freq in self.survey.freqs:
AT = self.getA(freq).T
ATinv = self.Solver(AT, **self.solverOpts)
for src in self.survey.getSrcByFreq(freq):
ftype = self._solutionType
f_src = f[src, :] # Need to fix this...
for rx in src.rxList:
# Get the adjoint evalDeriv
# PTv needs to be nE,
PTv = rx.evalDeriv(src, self.mesh, f, mkvc(v[src, rx],2), adjoint=True) # wrt u, need possibility wrt m
# Get the
dA_duIT = ATinv * PTv
dA_dmT = self.getADeriv_m(freq, f_src, mkvc(dA_duIT), adjoint=True)
dRHS_dmT = self.getRHSDeriv_m(freq, mkvc(dA_duIT), adjoint=True)
# Make du_dmT
if dRHS_dmT is None:
du_dmT = -dA_dmT
else:
du_dmT = -dA_dmT + dRHS_dmT
# Select the correct component
# du_dmT needs to be of size nC,
real_or_imag = rx.projComp
if real_or_imag == 'real':
Jtv += du_dmT.real
elif real_or_imag == 'imag':
Jtv += -du_dmT.real
else:
raise Exception('Must be real or imag')
# Clean the factorization, clear memory.
ATinv.clean()
return Jtv
###################################
## 1D problems
###################################
class Problem1D_ePrimSec(BaseNSEMProblem):
"""
A NSEM problem soving a e formulation and primary/secondary fields decomposion.
By eliminating the magnetic flux density using
.. math ::
\mathbf{b} = \\frac{1}{i \omega}\\left(-\mathbf{C} \mathbf{e} \\right)
we can write Maxwell's equations as a second order system in \\\(\\\mathbf{e}\\\) only:
.. math ::
\\left(\mathbf{C}^T \mathbf{M^e_{\mu^{-1}}} \mathbf{C} + i \omega \mathbf{M^f_\sigma}] \mathbf{e}_{s} =& i \omega \mathbf{M^f_{\delta \sigma}} \mathbf{e}_{p}
which we solve for \\\(\\\mathbf{e_s}\\\). The total field \\\mathbf{e}\\ = \\\mathbf{e_p}\\ + \\\mathbf{e_s}\\.
The primary field is estimated from a background model (commonly half space ).
"""
# From FDEMproblem: Used to project the fields. Currently not used for NSEMproblem.
_solutionType = 'e_1dSolution'
_formulation = 'EF'
fieldsPair = Fields1D_ePrimSec
# Initiate properties
_sigmaPrimary = None
def __init__(self, mesh, **kwargs):
BaseNSEMProblem.__init__(self, mesh, **kwargs)
# self._sigmaPrimary = sigmaPrimary
@property
def MeMui(self):
"""
Edge inner product matrix
"""
if getattr(self, '_MeMui', None) is None:
self._MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
return self._MeMui
@property
def MfSigma(self):
"""
Edge inner product matrix
"""
if getattr(self, '_MfSigma', None) is None:
self._MfSigma = self.mesh.getFaceInnerProduct(self.curModel.sigma)
return self._MfSigma
@property
def sigmaPrimary(self):
"""
A background model, use for the calculation of the primary fields.
"""
return self._sigmaPrimary
@sigmaPrimary.setter
def sigmaPrimary(self, val):
# Note: TODO add logic for val, make sure it is the correct size.
self._sigmaPrimary = val
def getA(self, freq):
"""
Function to get the A matrix.
:param float freq: Frequency
:rtype: scipy.sparse.csr_matrix
:return: A
"""
# Note: need to use the code above since in the 1D problem I want
# e to live on Faces(nodes) and h on edges(cells). Might need to rethink this
# Possible that _fieldType and _eqLocs can fix this
MeMui = self.MeMui
MfSigma = self.MfSigma
C = self.mesh.nodalGrad
# Make A
A = C.T*MeMui*C + 1j*omega(freq)*MfSigma
# Either return full or only the inner part of A
return A
def getADeriv_m(self, freq, u, v, adjoint=False):
"""
The derivative of A wrt sigma
"""
dsig_dm = self.curModel.sigmaDeriv
MeMui = self.MeMui
#
u_src = u['e_1dSolution']
dMfSigma_dm = self.mesh.getFaceInnerProductDeriv(self.curModel.sigma)(u_src) * self.curModel.sigmaDeriv
if adjoint:
return 1j * omega(freq) * ( dMfSigma_dm.T * v )
# Note: output has to be nN/nF, not nC/nE.
# v should be nC
return 1j * omega(freq) * ( dMfSigma_dm * v )
def getRHS(self, freq):
"""
Function to return the right hand side for the system.
:param float freq: Frequency
:rtype: numpy.ndarray (nF, 1), numpy.ndarray (nF, 1)
:return: RHS for 1 polarizations, primary fields
"""
# Get sources for the frequncy(polarizations)
Src = self.survey.getSrcByFreq(freq)[0]
S_e = Src.S_e(self)
return -1j * omega(freq) * S_e
def getRHSDeriv_m(self, freq, v, adjoint=False):
"""
The derivative of the RHS wrt sigma
"""
Src = self.survey.getSrcByFreq(freq)[0]
S_eDeriv = Src.S_eDeriv_m(self, v, adjoint)
return -1j * omega(freq) * S_eDeriv
def fields(self, m):
'''
Function to calculate all the fields for the model m.
:param np.ndarray (nC,) m: Conductivity model
'''
# Set the current model
self.curModel = m
# Make the fields object
F = self.fieldsPair(self.mesh, self.survey)
# Loop over the frequencies
for freq in self.survey.freqs:
if self.verbose:
startTime = time.time()
print 'Starting work for {:.3e}'.format(freq)
sys.stdout.flush()
A = self.getA(freq)
rhs = self.getRHS(freq)
Ainv = self.Solver(A, **self.solverOpts)
e_s = Ainv * rhs
# Store the fields
Src = self.survey.getSrcByFreq(freq)[0]
# NOTE: only store the e_solution(secondary), all other components calculated in the fields object
F[Src, 'e_1dSolution'] = e_s[:,-1] # Only storing the yx polarization as 1d
# Note curl e = -iwb so b = -curl e /iw
# 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)
sys.stdout.flush()
return F
# Note this is not fully functional.
# Missing:
# Fields class corresponding to the fields
# Update Jvec and Jtvec to include all the derivatives components
# Other things ...
class Problem1D_eTotal(BaseNSEMProblem):
"""
A NSEM problem solving a e formulation and a Total bondary domain decompostion.
Solves the equation:
Math:
Have to do this...
Not implement correctly.......
"""
# From FDEMproblem: Used to project the fields. Currently not used for NSEMproblem.
_solutionType = 'e_1dSolution'
_formulation = 'EF'
# fieldsPair = Fields1D_eTotal
def __init__(self, mesh, **kwargs):
BaseNSEMProblem.__init__(self, mesh, **kwargs)
@property
def MeMui(self):
"""
Edge inner product matrix
"""
if getattr(self, '_MeMui', None) is None:
self._MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
return self._MeMui
@property
def MfSigma(self):
"""
Edge inner product matrix
"""
if getattr(self, '_MfSigma', None) is None:
self._MfSigma = self.mesh.getFaceInnerProduct(self.curModel.sigma)
return self._MfSigma
def getA(self, freq, full=False):
"""
Function to get the A matrix.
:param float freq: Frequency
:param logic full: Return full A or the inner part
:rtype: scipy.sparse.csr_matrix
:return: A
"""
MeMui = self.MeMui
MfSigma = self.MfSigma
# Note: need to use the code above since in the 1D problem I want
# e to live on Faces(nodes) and h on edges(cells). Might need to rethink this
# Possible that _fieldType and _eqLocs can fix this
# MeMui = self.MfMui
# MfSigma = self.MfSigma
C = self.mesh.nodalGrad
# Make A
A = C.T*MeMui*C + 1j*omega(freq)*MfSigma
# Either return full or only the inner part of A
if full:
return A
else:
return A[1:-1,1:-1]
def getADeriv_m(self, freq, u, v, adjoint=False):
raise NotImplementedError('getADeriv is not implemented')
def getRHS(self, freq):
"""
Function to return the right hand side for the system.
:param float freq: Frequency
:rtype: numpy.ndarray (nE, 2), numpy.ndarray (nE, 2)
:return: RHS for both polarizations, primary fields
"""
# Get sources for the frequency
# NOTE: Need to use the source information, doesn't really apply in 1D
src = self.survey.getSrcByFreq(freq)
# Get the full A
A = self.getA(freq,full=True)
# Define the outer part of the solution matrix
Aio = A[1:-1,[0,-1]]
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
## Note: The analytic solution is derived with e^iwt
eBC = np.r_[Etot[0],Etot[-1]]
# The right hand side
return -Aio*eBC, eBC
def getRHSderiv_m(self, freq, backSigma, u, v, adjoint=False):
raise NotImplementedError('getRHSDeriv not implemented yet')
return None
def fields(self, m):
'''
Function to calculate all the fields for the model m.
:param np.ndarray (nC,) m: Conductivity model
:param np.ndarray (nC,) m_back: Background conductivity model
'''
self.curModel = m
# RHS, CalcFields = self.getRHS(freq,m_back), self.calcFields
F = Fields1D_eTotal(self.mesh, self.survey)
for freq in self.survey.freqs:
if self.verbose:
startTime = time.time()
print 'Starting work for {:.3e}'.format(freq)
sys.stdout.flush()
A = self.getA(freq)
rhs, e_o = self.getRHS(freq)
Ainv = self.Solver(A, **self.solverOpts)
e_i = Ainv * rhs
e = mkvc(np.r_[e_o[0], e_i, e_o[1]],2)
# Store the fields
Src = self.survey.getSrcByFreq(freq)
# NOTE: only store e fields
F[Src, 'e_1dSolution'] = e[:,0]
if self.verbose:
print 'Ran for {:f} seconds'.format(time.time()-startTime)
sys.stdout.flush()
return F
###################################
## 3D problems
###################################
class Problem3D_ePrimSec(BaseNSEMProblem):
"""
A NSEM problem solving a e formulation and a primary/secondary fields decompostion.
By eliminating the magnetic flux density using
.. math ::
\mathbf{b} = \\frac{1}{i \omega}\\left(-\mathbf{C} \mathbf{e} \\right)
we can write Maxwell's equations as a second order system in \\\(\\\mathbf{e}\\\) only:
.. math ::
\\left(\mathbf{C}^T \mathbf{M^f_{\mu^{-1}}} \mathbf{C} + i \omega \mathbf{M^e_\sigma}] \mathbf{e}_{s} =& i \omega \mathbf{M^e_{\delta \sigma}} \mathbf{e}_{p}
which we solve for \\\(\\\mathbf{e_s}\\\). The total field \\\mathbf{e}\\ = \\\mathbf{e_p}\\ + \\\mathbf{e_s}\\.
The primary field is estimated from a background model (commonly as a 1D model).
"""
# From FDEMproblem: Used to project the fields. Currently not used for NSEMproblem.
_solutionType = [ 'e_pxSolution', 'e_pySolution'] # Forces order on the object
_formulation = 'EB'
fieldsPair = Fields3D_ePrimSec
# Initiate properties
_sigmaPrimary = None
def __init__(self, mesh, **kwargs):
BaseNSEMProblem.__init__(self, mesh, **kwargs)
@property
def sigmaPrimary(self):
"""
A background model, use for the calculation of the primary fields.
"""
return self._sigmaPrimary
@sigmaPrimary.setter
def sigmaPrimary(self, val):
# Note: TODO add logic for val, make sure it is the correct size.
self._sigmaPrimary = val
def getA(self, freq):
"""
Function to get the A system.
:param float freq: Frequency
:rtype: scipy.sparse.csr_matrix
:return: A
"""
Mmui = self.MfMui
Msig = self.MeSigma
C = self.mesh.edgeCurl
return C.T*Mmui*C + 1j*omega(freq)*Msig
def getADeriv_m(self, freq, u, v, adjoint=False):
"""
Calculate the derivative of A wrt m.
"""
# Fix u to be a matrix nE,2
# This considers both polarizations and returns a nE,2 matrix for each polarization
if adjoint:
dMe_dsigV = sp.hstack(( self.MeSigmaDeriv( u['e_pxSolution'] ).T, self.MeSigmaDeriv(u['e_pySolution'] ).T ))*v
else:
# Need a nE,2 matrix to be returned
dMe_dsigV = np.hstack(( mkvc(self.MeSigmaDeriv( u['e_pxSolution'] )*v,2), mkvc( self.MeSigmaDeriv(u['e_pySolution'] )*v,2) ))
return 1j * omega(freq) * dMe_dsigV
def getRHS(self, freq):
"""
Function to return the right hand side for the system.
:param float freq: Frequency
:rtype: numpy.ndarray (nE, 2), numpy.ndarray (nE, 2)
:return: RHS for both polarizations, primary fields
"""
# Get sources for the frequncy(polarizations)
Src = self.survey.getSrcByFreq(freq)[0]
S_e = Src.S_e(self)
return -1j * omega(freq) * S_e
def getRHSDeriv_m(self, freq, v, adjoint=False):
"""
The derivative of the RHS with respect to sigma
"""
Src = self.survey.getSrcByFreq(freq)[0]
S_eDeriv = Src.S_eDeriv_m(self, v, adjoint)
return -1j * omega(freq) * S_eDeriv
def fields(self, m):
'''
Function to calculate all the fields for the model m.
:param np.ndarray (nC,) m: Conductivity model
'''
# Set the current model
self.curModel = m
F = self.fieldsPair(self.mesh, self.survey)
for freq in self.survey.freqs:
if self.verbose:
startTime = time.time()
print 'Starting work for {:.3e}'.format(freq)
sys.stdout.flush()
A = self.getA(freq)
rhs = self.getRHS(freq)
# Solve the system
Ainv = self.Solver(A, **self.solverOpts)
e_s = Ainv * rhs
# Store the fields
Src = self.survey.getSrcByFreq(freq)[0]
# Store the fields
# Use self._solutionType
F[Src, 'e_pxSolution'] = e_s[:,0]
F[Src, 'e_pySolution'] = e_s[:,1]
# Note curl e = -iwb so b = -curl/iw
if self.verbose:
print 'Ran for {:f} seconds'.format(time.time()-startTime)
sys.stdout.flush()
Ainv.clean()
return F
+11 -11
View File
@@ -11,9 +11,9 @@ import sys
### Sources ###
#################
class BaseMTSrc(FDEMBaseSrc):
class BaseNSEMSrc(FDEMBaseSrc):
'''
Sources for the MT problem.
Sources for the NSEM problem.
Use the SimPEG BaseSrc, since the source fields share properties with the transmitters.
:param float freq: The frequency of the source
@@ -29,28 +29,28 @@ class BaseMTSrc(FDEMBaseSrc):
FDEMBaseSrc.__init__(self, rxList)
# 1D sources
class polxy_1DhomotD(BaseMTSrc):
class polxy_1DhomotD(BaseNSEMSrc):
"""
MT source for both polarizations (x and y) for the total Domain.
NSEM source for both polarizations (x and y) for the total Domain.
It calculates fields calculated based on conditions on the boundary of the domain.
"""
def __init__(self, rxList, freq):
BaseMTSrc.__init__(self, rxList, freq)
BaseNSEMSrc.__init__(self, rxList, freq)
# TODO: need to add the primary fields calc and source terms into the problem.
# Need to implement such that it works for all dims.
class polxy_1Dprimary(BaseMTSrc):
class polxy_1Dprimary(BaseNSEMSrc):
"""
MT source for both polarizations (x and y) given a 1D primary models.
NSEM source for both polarizations (x and y) given a 1D primary models.
It assigns fields calculated from the 1D model as fields in the full space of the problem.
"""
def __init__(self, rxList, freq):
# assert mkvc(self.mesh.hz.shape,1) == mkvc(sigma1d.shape,1),'The number of values in the 1D background model does not match the number of vertical cells (hz).'
self.sigma1d = None
BaseMTSrc.__init__(self, rxList, freq)
BaseNSEMSrc.__init__(self, rxList, freq)
# Hidden property of the ePrimary
self._ePrimary = None
@@ -128,15 +128,15 @@ class polxy_1Dprimary(BaseMTSrc):
# v should be nC size
return MsigmaDeriv * v
class polxy_3Dprimary(BaseMTSrc):
class polxy_3Dprimary(BaseNSEMSrc):
"""
MT source for both polarizations (x and y) given a 3D primary model. It assigns fields calculated from the 1D model
NSEM source for both polarizations (x and y) given a 3D primary model. It assigns fields calculated from the 1D model
as fields in the full space of the problem.
"""
def __init__(self, rxList, freq):
# assert mkvc(self.mesh.hz.shape,1) == mkvc(sigma1d.shape,1),'The number of values in the 1D background model does not match the number of vertical cells (hz).'
self.sigmaPrimary = None
BaseMTSrc.__init__(self, rxList, freq)
BaseNSEMSrc.__init__(self, rxList, freq)
# Hidden property of the ePrimary
self._ePrimary = None
@@ -4,7 +4,7 @@ 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
import SrcNSEM
import sys
#################
@@ -63,9 +63,9 @@ class Rx(SimPEGsurvey.BaseRx):
'''
Project the fields to natural source data.
:param SrcMT src: The source of the fields to project
:param SrcNSEM src: The source of the fields to project
:param SimPEG.Mesh mesh:
:param FieldsMT f: Natural source fields object to project
:param FieldsNSEM f: Natural source fields object to project
'''
## NOTE: Assumes that e is on t
@@ -143,9 +143,9 @@ class Rx(SimPEGsurvey.BaseRx):
"""
The derivative of the projection wrt u
:param MTsrc src: MT source
:param NSEMsrc src: NSEM source
:param TensorMesh mesh: Mesh defining the topology of the problem
:param MTfields f: MT fields object of the source
:param NSEMfields f: NSEM fields object of the source
:param numpy.ndarray v: Random vector of size
"""
@@ -390,12 +390,12 @@ class Rx(SimPEGsurvey.BaseRx):
#################
class Survey(SimPEGsurvey.BaseSurvey):
"""
Survey class for MT. Contains all the sources associated with the survey.
Survey class for NSEM. Contains all the sources associated with the survey.
:param list srcList: List of sources associated with the survey
"""
srcPair = SrcMT.BaseMTSrc
srcPair = SrcNSEM.BaseNSEMSrc
def __init__(self, srcList, **kwargs):
# Sort these by frequency
@@ -443,7 +443,7 @@ class Survey(SimPEGsurvey.BaseSurvey):
#################
class Data(SimPEGsurvey.Data):
'''
Data class for MTdata. Stores the data vector indexed by the survey.
Data class for NSEMdata. Stores the data vector indexed by the survey.
:param SimPEG survey object survey:
:param v vector of the data in order matching of the survey
@@ -461,7 +461,7 @@ class Data(SimPEGsurvey.Data):
def toRecArray(self,returnType='RealImag'):
'''
Function that returns a numpy.recarray for a SimpegMT impedance data object.
Function that returns a numpy.recarray for a SimpegNSEM impedance data object.
:param str returnType: Switches between returning a rec array where the impedance is split to real and imaginary ('RealImag') or is a complex ('Complex')
@@ -483,7 +483,7 @@ class Data(SimPEGsurvey.Data):
locs = np.hstack((np.array([[0.0]]),locs))
tArrRec = np.concatenate((src.freq*np.ones((locs.shape[0],1)),locs,np.nan*np.ones((locs.shape[0],12))),axis=1).view(dtRI)
# np.array([(src.freq,rx.locs[0,0],rx.locs[0,1],rx.locs[0,2],np.nan ,np.nan ,np.nan ,np.nan ,np.nan ,np.nan ,np.nan ,np.nan ) for rx in src.rxList],dtype=dtRI)
# Get the type and the value for the DataMT object as a list
# Get the type and the value for the DataNSEM object as a list
typeList = [[rx.rxType.replace('z1d','zyx'),self[src,rx]] for rx in src.rxList]
# Insert the values to the temp array
for nr,(key,val) in enumerate(typeList):
@@ -517,17 +517,17 @@ class Data(SimPEGsurvey.Data):
@classmethod
def fromRecArray(cls, recArray, srcType='primary'):
"""
Class method that reads in a numpy record array to MTdata object.
Class method that reads in a numpy record array to NSEMdata object.
Only imports the impedance data.
"""
if srcType=='primary':
src = SrcMT.polxy_1Dprimary
src = SrcNSEM.polxy_1Dprimary
elif srcType=='total':
src = SrcMT.polxy_1DhomotD
src = SrcNSEM.polxy_1DhomotD
else:
raise NotImplementedError('{:s} is not a valid source type for MTdata')
raise NotImplementedError('{:s} is not a valid source type for NSEMdata')
# Find all the frequencies in recArray
uniFreq = np.unique(recArray['freq'])
@@ -3,7 +3,7 @@
import numpy as np, SimPEG as simpeg
from scipy.constants import mu_0, epsilon_0 as eps_0
def getEHfields(m1d,sigma,freq,zd,scaleUD=True):
def getEHfields(m1d,sigma,freq,zd,scaleUD=True,scaleValue=1):
'''Analytic solution for MT 1D layered earth. Returns E and H fields.
:param SimPEG.mesh, object m1d: Mesh object with the 1D spatial information.
@@ -12,7 +12,7 @@ def getEHfields(m1d,sigma,freq,zd,scaleUD=True):
:param numpy array, vector zd: location to calculate EH fields at
:param bollean, scaleUD: scales the output to be 1 at the top, increases numeracal stability.
Assumes a halfspace with the same conductive as the last cell below.
Assumes a halfspace with the same conductive as the deepest cell.
'''
# Note add an error check for the mesh and sigma are the same size.
@@ -29,7 +29,7 @@ def getEHfields(m1d,sigma,freq,zd,scaleUD=True):
# Initiate the propagation matrix, in the order down up.
UDp = np.zeros((2,m1d.nC+1),dtype=complex)
UDp[1,0] = 1. # Set the wave amplitude as 1 into the half-space at the bottom of the mesh
UDp[1,0] = scaleValue # Set the wave amplitude as 1 into the half-space at the bottom of the mesh
# 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
@@ -38,9 +38,9 @@ def getEHfields(m1d,sigma,freq,zd,scaleUD=True):
# Build the propagation matrix
# Convert fields to down/up going components in layer below current layer
Pj1 = np.array([[1,1],[yp1,-yp1]])
Pj1 = np.array([[1,1],[yp1,-yp1]],dtype=complex)
# Convert fields to down/up going components in current layer
Pjinv = 1./2*np.array([[1,zp],[1,-zp]])
Pjinv = 1./2*np.array([[1,zp],[1,-zp]],dtype=complex)
# Propagate down and up components through the current layer
elamh = np.array([[np.exp(-1j*k[lnr+1]*h),0],[0,np.exp(1j*k[lnr+1]*h)]])
@@ -48,7 +48,14 @@ 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]
# Scale the values such that 1 at the top
scaleVal = UDp[:,lnr+1::-1]/UDp[1,lnr+1]
if np.any(np.isnan(scaleVal)):
# If there is a nan (thickness very great), rebuild the move up cell
scaleVal = np.zeros_like(UDp[:,lnr+1::-1],dtype=complex)
scaleVal[1,0] = scaleValue
UDp[:,lnr+1::-1] = scaleVal
# Calculate the fields
Ed = np.empty((zd.size,),dtype=complex)
+5
View File
@@ -0,0 +1,5 @@
from MT1Dsolutions import get1DEfields # Add the names of the functions
from MT1Danalytic import getEHfields, getImpedance
from dataUtils import *
from ediFilesUtils import *
from testUtils import *
@@ -5,25 +5,25 @@ import numpy.lib.recfunctions as recFunc
from scipy.constants import mu_0
from scipy import interpolate as sciint
def getAppRes(MTdata):
def getAppRes(NSEMdata):
# Make impedance
zList = []
for src in MTdata.survey.srcList:
for src in NSEMdata.survey.srcList:
zc = [src.freq]
for rx in src.rxList:
if 'i' in rx.rxType:
m=1j
else:
m = 1
zc.append(m*MTdata[src,rx])
zc.append(m*NSEMdata[src,rx])
zList.append(zc)
return [appResPhs(zList[i][0],np.sum(zList[i][1:3])) for i in np.arange(len(zList))]
def rotateData(MTdata, rotAngle):
def rotateData(NSEMdata,rotAngle):
'''
Function that rotates clockwist by rotAngle (- negative for a counter-clockwise rotation)
'''
recData = MTdata.toRecArray('Complex')
recData = NSEMdata.toRecArray('Complex')
impData = rec2ndarr(recData[['zxx','zxy','zyx','zyy']],complex)
# Make the rotation matrix
# c,s,zxx,zxy,zyx,zyy = sympy.symbols('c,s,zxx,zxy,zyx,zyy')
@@ -40,27 +40,27 @@ def rotateData(MTdata, rotAngle):
for nr,comp in enumerate(['zxx','zxy','zyx','zyy']):
outRec[comp] = rotData[:,nr]
from SimPEG import MT
return MT.Data.fromRecArray(outRec)
from SimPEG import NSEM
return NSEM.Data.fromRecArray(outRec)
def appResPhs(freq, z):
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)
return app_res, app_phs
def skindepth(rho, freq):
def skindepth(rho,freq):
''' Function to calculate the skindepth of EM waves'''
return np.sqrt( (rho*((1/(freq * mu_0 * np.pi )))))
def rec2ndarr(x, dt=float):
def rec2ndarr(x,dt=float):
return x.view((dt, len(x.dtype.names)))
def makeAnalyticSolution(mesh, model, elev, freqs):
from SimPEG import MT
def makeAnalyticSolution(mesh,model,elev,freqs):
from SimPEG import NSEM
data1D = []
for freq in freqs:
anaEd, anaEu, anaHd, anaHu = MT.Utils.MT1Danalytic.getEHfields(mesh,model,freq,elev)
anaEd, anaEu, anaHd, anaHu = NSEM.Utils.MT1Danalytic.getEHfields(mesh,model,freq,elev)
anaE = anaEd+anaEu
anaH = anaHd+anaHu
@@ -70,8 +70,8 @@ def makeAnalyticSolution(mesh, model, elev, freqs):
dataRec = np.array(data1D,dtype=[('freq',float),('x',float),('y',float),('z',float),('zyx',complex)])
return dataRec
def plotMT1DModelData(problem, models, symList=None):
from SimPEG import MT
def plotMT1DModelData(problem,models,symList=None):
from SimPEG import NSEM
# Setup the figure
fontSize = 15
@@ -79,7 +79,7 @@ def plotMT1DModelData(problem, models, symList=None):
axM = fig.add_axes([0.075,.1,.25,.875])
axM.set_xlabel('Resistivity [Ohm*m]',fontsize=fontSize)
axM.set_xlim(1e-1,1e5)
axM.set_ylim(-10000,5000)
# axM.set_ylim(-10000,5000)
axM.set_ylabel('Depth [km]',fontsize=fontSize)
axR = fig.add_axes([0.42,.575,.5,.4])
axR.set_xscale('log')
@@ -132,38 +132,94 @@ def plotMT1DModelData(problem, models, symList=None):
freq = simpeg.mkvc(data1D['freq'],2)
res, phs = appResPhs(freq,allData)
stdCol = 'gray'
axRtw = axR.twinx()
axRtw.set_ylabel('Std of log10',color=stdCol)
[(t.set_color(stdCol), t.set_rotation(-45)) for t in axRtw.get_yticklabels()]
axPtw = axP.twinx()
axPtw.set_ylabel('Std ',color=stdCol)
[t.set_color(stdCol) for t in axPtw.get_yticklabels()]
axRtw.plot(freq, np.std(np.log10(res),1),'--',color=stdCol)
axPtw.plot(freq, np.std(phs,1),'--',color=stdCol)
if False:
stdCol = 'gray'
axRtw = axR.twinx()
axRtw.set_ylabel('Std of log10',color=stdCol)
[(t.set_color(stdCol), t.set_rotation(-45)) for t in axRtw.get_yticklabels()]
axPtw = axP.twinx()
axPtw.set_ylabel('Std ',color=stdCol)
[t.set_color(stdCol) for t in axPtw.get_yticklabels()]
axRtw.plot(freq, np.std(np.log10(res),1),'--',color=stdCol)
axPtw.plot(freq, np.std(phs,1),'--',color=stdCol)
# Fix labels and ticks
yMtick = [l/1000 for l in axM.get_yticks().tolist()]
axM.set_yticklabels(yMtick)
# yMtick = [l/1000 for l in axM.get_yticks().tolist()]
# axM.set_yticklabels(yMtick)
[ l.set_rotation(90) for l in axM.get_yticklabels()]
[ l.set_rotation(90) for l in axR.get_yticklabels()]
[(t.set_color(stdCol), t.set_rotation(-45)) for t in axRtw.get_yticklabels()]
[t.set_color(stdCol) for t in axPtw.get_yticklabels()]
# [(t.set_color(stdCol), t.set_rotation(-45)) for t in axRtw.get_yticklabels()]
# [t.set_color(stdCol) for t in axPtw.get_yticklabels()]
for ax in [axM,axR,axP]:
ax.xaxis.set_tick_params(labelsize=fontSize)
ax.yaxis.set_tick_params(labelsize=fontSize)
return fig
def plotImpAppRes(dataArrays,plotLoc,textStr=[]):
''' Plots amplitude impedance and phase'''
# fig = plt.figure(1,(7, 7))
import plotDataTypes as pDt
# axes = ImageGrid(fig, (0.05,0.05,0.875,0.875),nrows_ncols = (2, 2),axes_pad = 0.25,add_all=True,share_all=True,label_mode = "L")
# Make the figure and axes
fig,axT=plt.subplots(2,2,sharex=True)
axes = axT.ravel()
fig.set_size_inches((13.5,7.0))
fig.suptitle('{:s}\nStation at: {:.1f}x ; {:.1f}y'.format(textStr,plotLoc[0],plotLoc[1]))
# Have to deal with axes
# Set log
for ax in axes.ravel():
ax.set_xscale('log')
axes[0].invert_xaxis()
axes[0].set_yscale('log')
axes[2].set_yscale('log')
# Set labels
axes[2].set_xlabel('Frequency [Hz]')
axes[3].set_xlabel('Frequency [Hz]')
axes[0].set_ylabel('Apperent resistivity [Ohm m]')
axes[1].set_ylabel('Apperent phase [degrees]')
axes[1].set_ylim(-180,180)
axes[2].set_ylabel('Impedance amplitude [V/A]')
axes[3].set_ylim(-180,180)
axes[3].set_ylabel('Impedance angle [degrees]')
# Plot the data
for nr,dataArray in enumerate(dataArrays):
if nr==1:
parSym = '*'
else:
parSym = 's'
# app res
pDt.plotIsoStaImpedance(axes[0],plotLoc,dataArray,'zxy',par='res',pSym=parSym)
pDt.plotIsoStaImpedance(axes[0],plotLoc,dataArray,'zyx',par='res',pSym=parSym)
# app phs
pDt.plotIsoStaImpedance(axes[1],plotLoc,dataArray,'zxy',par='phs',pSym=parSym)
pDt.plotIsoStaImpedance(axes[1],plotLoc,dataArray,'zyx',par='phs',pSym=parSym)
# imp abs
pDt.plotIsoStaImpedance(axes[2],plotLoc,dataArray,'zxx',par='abs',pSym=parSym)
pDt.plotIsoStaImpedance(axes[2],plotLoc,dataArray,'zxy',par='abs',pSym=parSym)
pDt.plotIsoStaImpedance(axes[2],plotLoc,dataArray,'zyx',par='abs',pSym=parSym)
pDt.plotIsoStaImpedance(axes[2],plotLoc,dataArray,'zyy',par='abs',pSym=parSym)
# imp abs
pDt.plotIsoStaImpedance(axes[3],plotLoc,dataArray,'zxx',par='phs',pSym=parSym)
pDt.plotIsoStaImpedance(axes[3],plotLoc,dataArray,'zxy',par='phs',pSym=parSym)
pDt.plotIsoStaImpedance(axes[3],plotLoc,dataArray,'zyx',par='phs',pSym=parSym)
pDt.plotIsoStaImpedance(axes[3],plotLoc,dataArray,'zyy',par='phs',pSym=parSym)
return fig,axes
def printTime():
import time
print time.strftime("%a, %d %b %Y %H:%M:%S +0000", time.localtime())
def convert3Dto1Dobject(MTdata,rxType3D='zyx'):
from SimPEG import MT
def convert3Dto1Dobject(NSEMdata,rxType3D='zyx'):
from SimPEG import NSEM
# Find the unique locations
# Need to find the locations
recDataTemp = MTdata.toRecArray()
recDataTemp = NSEMdata.toRecArray()
# Check if survey.std has been assigned.
## NEED TO: write this...
# Calculte and add the DET of the tensor to the recArray
@@ -185,24 +241,24 @@ def convert3Dto1Dobject(MTdata,rxType3D='zyx'):
# Make the receiver list
rx1DList = []
for rxType in ['z1dr','z1di']:
rx1DList.append(MT.Rx(simpeg.mkvc(loc,2).T,rxType))
rx1DList.append(NSEM.Rx(simpeg.mkvc(loc,2).T,rxType))
# Source list
locrecData = recData[np.sqrt(np.sum( (rec2ndarr(recData[['x','y','z']]).data - loc )**2,axis=1)) < 1e-5]
dat1DList = []
src1DList = []
for freq in locrecData['freq']:
src1DList.append(MT.SrcMT.src_polxy_1Dprimary(rx1DList,freq))
src1DList.append(NSEM.SrcNSEM.src_polxy_1Dprimary(rx1DList,freq))
for comp in ['r','i']:
dat1DList.append( corr * locrecData[rxType3D+comp][locrecData['freq']== freq].data )
# Make the survey
sur1D = MT.Survey(src1DList)
sur1D = NSEM.Survey(src1DList)
# Make the data
dataVec = np.hstack(dat1DList)
dat1D = MT.Data(sur1D,dataVec)
dat1D = NSEM.Data(sur1D,dataVec)
sur1D.dobs = dataVec
# Need to take MTdata.survey.std and split it as well.
# Need to take NSEMdata.survey.std and split it as well.
std=0.05
sur1D.std = np.abs(sur1D.dobs*std) #+ 0.01*np.linalg.norm(sur1D.dobs)
mtData1DList.append(dat1D)
@@ -210,29 +266,29 @@ def convert3Dto1Dobject(MTdata,rxType3D='zyx'):
# Return the the list of data.
return mtData1DList
def resampleMTdataAtFreq(MTdata,freqs):
def resampleNSEMdataAtFreq(NSEMdata,freqs):
"""
Function to resample MTdata at set of frequencies
Function to resample NSEMdata at set of frequencies
"""
from SimPEG import MT
from SimPEG import NSEM
# Make a rec array
MTrec = MTdata.toRecArray().data
NSEMrec = NSEMdata.toRecArray().data
# Find unique locations
uniLoc = np.unique(MTrec[['x','y','z']])
uniFreq = MTdata.survey.freqs
uniLoc = np.unique(NSEMrec[['x','y','z']])
uniFreq = NSEMdata.survey.freqs
# Get the comps
dNames = MTrec.dtype
dNames = NSEMrec.dtype
# Loop over all the locations and interpolate
for loc in uniLoc:
# Find the index of the station
ind = np.sqrt(np.sum((rec2ndarr(MTrec[['x','y','z']]) - rec2ndarr(loc))**2,axis=1)) < 1. # Find dist of 1 m accuracy
ind = np.sqrt(np.sum((rec2ndarr(NSEMrec[['x','y','z']]) - rec2ndarr(loc))**2,axis=1)) < 1. # Find dist of 1 m accuracy
# Make a temporary recArray and interpolate all the components
tArrRec = np.concatenate((simpeg.mkvc(freqs,2),np.ones((len(freqs),1))*rec2ndarr(loc),np.nan*np.ones((len(freqs),12))),axis=1).view(dNames)
for comp in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi','tzxr','tzxi','tzyr','tzyi']:
int1d = sciint.interp1d(MTrec[ind]['freq'],MTrec[ind][comp],bounds_error=False)
int1d = sciint.interp1d(NSEMrec[ind]['freq'],NSEMrec[ind][comp],bounds_error=False)
tArrRec[comp] = simpeg.mkvc(int1d(freqs),2)
# Join together
@@ -241,5 +297,5 @@ def resampleMTdataAtFreq(MTdata,freqs):
except NameError as e:
outRecArr = tArrRec
# Make the MTdata and return
return MT.Data.fromRecArray(outRecArr)
# Make the NSEMdata and return
return NSEM.Data.fromRecArray(outRecArr)
@@ -2,7 +2,7 @@
from SimPEG import mkvc
from scipy.constants import mu_0
from numpy.lib import recfunctions as recFunc
from SimPEG.MT.Utils.dataUtils import rec2ndarr
from SimPEG.NSEM.Utils.dataUtils import rec2ndarr
# Import modules
import numpy as np
@@ -12,7 +12,7 @@ def homo1DModelSource(mesh,freq,sigma_1d):
'''
# import
from SimPEG.MT.Utils import get1DEfields
from SimPEG.NSEM.Utils import get1DEfields
# Get a 1d solution for a halfspace background
if mesh.dim == 1:
mesh1d = mesh
@@ -77,7 +77,7 @@ def analytic1DModelSource(mesh,freq,sigma_1d):
'''
# import
from SimPEG.MT.Utils import getEHfields
from SimPEG.NSEM.Utils import getEHfields
# Get a 1d solution for a halfspace background
if mesh.dim == 1:
mesh1d = mesh
+198
View File
@@ -0,0 +1,198 @@
import unittest
import sys
from scipy.constants import mu_0
import SimPEG as simpeg
from SimPEG.Utils import meshTensor
import numpy as np
np.random.seed(1100)
# Define the tolerances
TOLr = 5e-2
TOLp = 5e-2
def getAppResPhs(NSEMdata):
# Make impedance
from SimPEG.NSEM.Utils import appResPhs
zList = []
for src in NSEMdata.survey.srcList:
zc = [src.freq]
for rx in src.rxList:
if 'i' in rx.rxType:
m=1j
else:
m = 1
zc.append(m*NSEMdata[src,rx])
zList.append(zc)
return [appResPhs(zList[i][0],np.sum(zList[i][1:3])) for i in np.arange(len(zList))]
def setup1DSurvey(sigmaHalf,tD=True,structure=False):
from SimPEG import NSEM
# Frequency
nFreq = 33
freqs = np.logspace(3,-3,nFreq)
# Make the mesh
ct = 5
air = meshTensor([(ct,25,1.3)])
# coreT0 = meshTensor([(ct,15,1.2)])
# coreT1 = np.kron(meshTensor([(coreT0[-1],15,1.3)]),np.ones((7,)))
core = np.concatenate( ( np.kron(meshTensor([(ct,15,-1.2)]),np.ones((10,))) , meshTensor([(ct,20)]) ) )
bot = meshTensor([(core[0],20,-1.3)])
x0 = -np.array([np.sum(np.concatenate((core,bot)))])
m1d = simpeg.Mesh.TensorMesh([np.concatenate((bot,core,air))], x0=x0)
# Make the model
sigma = np.zeros(m1d.nC) + sigmaHalf
sigma[m1d.gridCC > 0 ] = 1e-8
sigmaBack = sigma.copy()
# Add structure
if structure:
shallow = (m1d.gridCC < -200) * (m1d.gridCC > -600)
deep = (m1d.gridCC < -3000) * (m1d.gridCC > -5000)
sigma[shallow] = 1
sigma[deep] = 0.1
rxList = []
for rxType in ['z1dr','z1di']:
rxList.append(NSEM.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType))
# Source list
srcList =[]
if tD:
for freq in freqs:
srcList.append(NSEM.SrcNSEM.polxy_1DhomotD(rxList,freq))
else:
for freq in freqs:
srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq))
survey = NSEM.Survey(srcList)
return survey, sigma, m1d
def setupSimpegNSEM_ePrimSec(inputSetup,comp='Imp',singleFreq=False,expMap=True):
from SimPEG import NSEM
M,freqs,sig,sigBG,rx_loc = inputSetup
# Make a receiver list
rxList = []
if comp == 'All':
for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi','tzxr','tzxi','tzyr','tzyi']:
rxList.append(NSEM.Rx(rx_loc,rxType))
elif comp == 'Imp':
for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi']:
rxList.append(NSEM.Rx(rx_loc,rxType))
elif comp == 'Tip':
for rxType in ['tzxr','tzxi','tzyr','tzyi']:
rxList.append(NSEM.Rx(rx_loc,rxType))
else:
rxList.append(NSEM.Rx(rx_loc,comp))
# Source list
srcList =[]
if singleFreq:
srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,singleFreq))
else:
for freq in freqs:
srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq))
# Survey NSEM
survey = NSEM.Survey(srcList)
## Setup the problem object
sigma1d = M.r(sigBG,'CC','CC','M')[0,0,:]
if expMap:
problem = NSEM.Problem3D_ePrimSec(M,sigmaPrimary= np.log(sigma1d) )
problem.mapping = simpeg.Maps.ExpMap(problem.mesh)
problem.curModel = np.log(sig)
else:
problem = NSEM.Problem3D_ePrimSec(M,sigmaPrimary= sigma1d)
problem.curModel = sig
problem.pair(survey)
problem.verbose = False
try:
from pymatsolver import MumpsSolver
problem.Solver = MumpsSolver
except:
pass
return (survey, problem)
def getInputs():
"""
Function that returns Mesh, freqs, rx_loc, elev.
"""
# Make a mesh
# M = simpeg.Mesh.TensorMesh([[(100,5,-1.5),(100.,10),(100,5,1.5)],[(100,5,-1.5),(100.,10),(100,5,1.5)],[(100,5,1.6),(100.,10),(100,3,2)]], x0=['C','C',-3529.5360])
# M = simpeg.Mesh.TensorMesh([[(1000,6,-1.5),(1000.,6),(1000,6,1.5)],[(1000,6,-1.5),(1000.,2),(1000,6,1.5)],[(1000,6,-1.3),(1000.,6),(1000,6,1.3)]], x0=['C','C','C'])# Setup the model
M = simpeg.Mesh.TensorMesh([[(200,6,-1.5),(200.,4),(200,6,1.5)],[(200,6,-1.5),(200.,4),(200,6,1.5)],[(200,8,-1.5),(200.,8),(200,8,1.5)]], x0=['C','C','C'])# Setup the model
# Set the frequencies
freqs = np.logspace(1,-3,5)
elev = 0
## Setup the the survey object
# Receiver locations
rx_x, rx_y = np.meshgrid(np.arange(-350,350,200),np.arange(-350,350,200))
rx_loc = np.hstack((simpeg.Utils.mkvc(rx_x,2),simpeg.Utils.mkvc(rx_y,2),elev+np.zeros((np.prod(rx_x.shape),1))))
return M, freqs, rx_loc, elev
def random(conds):
''' Returns a halfspace model based on the inputs'''
M, freqs, rx_loc, elev = getInputs()
# Backround
sigBG = np.ones(M.nC)*conds
# Add randomness to the model (10% of the value).
sig = np.exp( np.log(sigBG) + np.random.randn(M.nC)*(conds)*1e-1 )
return (M, freqs, sig, sigBG, rx_loc)
def halfSpace(conds):
''' Returns a halfspace model based on the inputs'''
M, freqs, rx_loc, elev = getInputs()
# Model
ccM = M.gridCC
# conds = [1e-2]
groundInd = ccM[:,2] < elev
sig = np.zeros(M.nC) + 1e-8
sig[groundInd] = conds
# Set the background, not the same as the model
sigBG = np.zeros(M.nC) + 1e-8
sigBG[groundInd] = conds
return (M, freqs, sig, sigBG, rx_loc)
def blockInhalfSpace(conds):
''' Returns a halfspace model based on the inputs'''
M, freqs, rx_loc, elev = getInputs()
# Model
ccM = M.gridCC
# conds = [1e-2]
groundInd = ccM[:,2] < elev
sig = simpeg.Utils.ModelBuilder.defineBlock(M.gridCC,np.array([-1000,-1000,-1500]),np.array([1000,1000,-1000]),conds)
sig[~groundInd] = 1e-8
# Set the background, not the same as the model
sigBG = np.zeros(M.nC) + 1e-8
sigBG[groundInd] = conds[1]
return (M, freqs, sig, sigBG, rx_loc)
def twoLayer(conds):
''' Returns a 2 layer model based on the conductivity values given'''
M, freqs, rx_loc, elev = getInputs()
# Model
ccM = M.gridCC
groundInd = ccM[:,2] < elev
botInd = ccM[:,2] < -3000
sig = np.zeros(M.nC) + 1e-8
sig[groundInd] = conds[1]
sig[botInd] = conds[0]
# Set the background, not the same as the model
sigBG = np.zeros(M.nC) + 1e-8
sigBG[groundInd] = conds[1]
return (M, freqs, sig, sigBG, rx_loc)
+5
View File
@@ -0,0 +1,5 @@
import Utils
from SurveyNSEM import Rx, Survey, Data
from FieldsNSEM import Fields1D_ePrimSec, Fields3D_ePrimSec
from ProblemNSEM import Problem1D_ePrimSec, Problem3D_ePrimSec
import SrcNSEM
+5 -5
View File
@@ -131,7 +131,7 @@ class Minimize(object):
Minimizes the function (evalFunction) starting at the location x0.
:param callable evalFunction: function handle that evaluates: f, g, H = F(x)
:param def evalFunction: function handle that evaluates: f, g, H = F(x)
:param numpy.ndarray x0: starting location
:rtype: numpy.ndarray
:return: x, the last iterate of the optimization algorithm
@@ -372,8 +372,8 @@ class Minimize(object):
Else, a modifySearchDirectionBreak call is preformed.
:param numpy.ndarray p: searchDirection
:rtype: tuple
:return: (xt, passLS) numpy.ndarray, bool
:rtype: numpy.ndarray,bool
:return: (xt, passLS)
"""
# Projected Armijo linesearch
self._LS_t = 1
@@ -408,8 +408,8 @@ class Minimize(object):
evalFunction returns a False indicating the break was not caught.
:param numpy.ndarray p: searchDirection
:rtype: tuple
:return: (xt, breakCaught) numpy.ndarray, bool
:rtype: numpy.ndarray,bool
:return: (xt, breakCaught)
"""
self.printDone(inLS=True)
print 'The linesearch got broken. Boo.'
+1 -1
View File
@@ -187,7 +187,7 @@ class _PropMapMetaClass(type):
attrs[attr + 'Model'] = prop._getModelProperty()
attrs[attr + 'Deriv'] = prop._getModelDerivProperty()
return type('PropModel', (PropModel, ), attrs)
return type(name.replace('PropMap', 'PropModel'), (PropModel, ), attrs)
class PropMap(object):
+6 -6
View File
@@ -10,7 +10,7 @@ class RegularizationMesh(object):
are not necessarily true differential operators, but are constructed from
a SimPEG Mesh.
:param BaseMesh mesh: problem mesh
:param Mesh mesh: problem mesh
:param numpy.array indActive: bool array, size nC, that is True where we have active cells. Used to reduce the operators so we regularize only on active cells
"""
@@ -383,8 +383,8 @@ class BaseRegularization(object):
:param numpy.array m: geophysical model
:param numpy.array v: vector to multiply
:rtype: scipy.sparse.csr_matrix
:return: WtW, or if v is supplied WtW*v (numpy.ndarray)
:rtype: scipy.sparse.csr_matrix or numpy.ndarray
:return: WtW or WtW*v
The regularization is:
@@ -650,8 +650,8 @@ class Tikhonov(Simple):
Note if the key word argument `mrefInSmooth` is False, then mref is not
included in the smoothness contribution.
:param BaseMesh mesh: SimPEG mesh
:param IdentityMap mapping: regularization mapping, takes the model from model space to the thing you want to regularize
:param Mesh mesh: SimPEG mesh
:param Maps mapping: regularization mapping, takes the model from model space to the thing you want to regularize
:param numpy.ndarray indActive: active cell indices for reducing the size of differential operators in the definition of a regularization mesh
:param bool mrefInSmooth: (default = False) put mref in the smoothness component?
:param float alpha_s: (default 1e-6) smallness weight
@@ -671,7 +671,7 @@ class Tikhonov(Simple):
alpha_yy = Utils.dependentProperty('_alpha_yy', 0.0, ['_W', '_Wyy'], "Weight for the second derivative in the y direction")
alpha_zz = Utils.dependentProperty('_alpha_zz', 0.0, ['_W', '_Wzz'], "Weight for the second derivative in the z direction")
def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
def __init__(self, mesh, mapping=None, indActive = None, **kwargs):
BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
@property
+3 -2
View File
@@ -311,6 +311,7 @@ class BaseSurvey(object):
if f is None: f = self.prob.fields(m)
return Utils.mkvc(self.eval(f))
@Utils.count
def eval(self, f):
"""eval(f)
@@ -321,7 +322,7 @@ class BaseSurvey(object):
d_\\text{pred} = \mathbf{P} f(m)
"""
raise NotImplementedError('eval is not yet implemented.')
raise NotImplemented('eval is not yet implemented.')
@Utils.count
def evalDeriv(self, f):
@@ -333,7 +334,7 @@ class BaseSurvey(object):
\\frac{\partial d_\\text{pred}}{\partial u} = \mathbf{P}
"""
raise NotImplementedError('eval is not yet implemented.')
raise NotImplemented('eval is not yet implemented.')
@Utils.count
def residual(self, m, f=None):
+1 -1
View File
@@ -237,7 +237,7 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tole
Compares error decay of 0th and 1st order Taylor approximation at point
x0 for a randomized search direction.
:param callable fctn: function handle
:param lambda fctn: function handle
:param numpy.array x0: point at which to check derivative
:param int num: number of times to reduce step length, h
:param bool plotIt: if you would like to plot
+13 -13
View File
@@ -7,11 +7,11 @@ def addBlock(gridCC, modelCC, p0, p1, blockProp):
"""
Add a block to an exsisting cell centered model, modelCC
:param numpy.array gridCC: mesh.gridCC is the cell centered grid
:param numpy.array modelCC: cell centered model
:param numpy.array p0: bottom, southwest corner of block
:param numpy.array p1: top, northeast corner of block
:blockProp float blockProp: property to assign to the model
:param numpy.array, gridCC: mesh.gridCC is the cell centered grid
:param numpy.array, modelCC: cell centered model
:param numpy.array, p0: bottom, southwest corner of block
:param numpy.array, p1: top, northeast corner of block
:blockProp float, blockProp: property to assign to the model
:return numpy.array, modelBlock: model with block
"""
@@ -147,7 +147,7 @@ def getIndicesSphere(center,radius,ccMesh):
if dimMesh == 1:
# Define the reference points
ind = np.abs(center[0] - ccMesh[:,0]) < radius
elif dimMesh == 2:
@@ -222,14 +222,14 @@ def layeredModel(ccMesh, layerTops, layerValues):
:param numpy.array ccMesh: cell-centered mesh
:param numpy.array layerTops: z-locations of the tops of each layer
:param numpy.array layerValue: values of the property to assign for each layer (starting at the top)
:param numpy.array layerValue: values of the property to assign for each layer (starting at the top)
:rtype: numpy.array
:return: M, layered model on the mesh
:return: M, layered model on the mesh
"""
descending = np.linalg.norm(sorted(layerTops, reverse=True) - layerTops) < 1e-20
# TODO: put an error check to make sure that there is an ordering... needs to work with inf elts
# TODO: put an error check to make sure that there is an ordering... needs to work with inf elts
# assert ascending or descending, "Layers must be listed in either ascending or descending order"
# start from bottom up
@@ -253,10 +253,10 @@ def layeredModel(ccMesh, layerTops, layerValues):
model = np.zeros(ccMesh.shape[0])
for i, top in enumerate(layerTops):
zind = z <= top
zind = z <= top
model[zind] = layerValues[i]
return model
return model
@@ -265,9 +265,9 @@ def randomModel(shape, seed=None, anisotropy=None, its=100, bounds=None):
Create a random model by convolving a kernel with a
uniformly distributed model.
:param tuple shape: shape of the model.
:param int,tuple shape: shape of the model.
:param int seed: pick which model to produce, prints the seed if you don't choose.
:param numpy.ndarray anisotropy: this is the (3 x n) blurring kernel that is used.
:param numpy.ndarray,list anisotropy: this is the (3 x n) blurring kernel that is used.
:param int its: number of smoothing iterations
:param list bounds: bounds on the model, len(list) == 2
:rtype: numpy.ndarray
+7 -7
View File
@@ -13,7 +13,7 @@ def _checkAccuracy(A, b, X, accuracyTol):
warnings.warn(msg, RuntimeWarning)
def SolverWrapD(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-6, name=None):
def SolverWrapD(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-6):
"""
Wraps a direct Solver.
@@ -72,11 +72,11 @@ def SolverWrapD(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-6, name=
if factorize and hasattr(self.solver, 'clean'):
return self.solver.clean()
return type(name if name is not None else fun.__name__, (object,), {"__init__": __init__, "clean": clean, "__mul__": __mul__})
return type(fun.__name__+'_Wrapped', (object,), {"__init__": __init__, "clean": clean, "__mul__": __mul__})
def SolverWrapI(fun, checkAccuracy=True, accuracyTol=1e-5, name=None):
def SolverWrapI(fun, checkAccuracy=True, accuracyTol=1e-5):
"""
Wraps an iterative Solver.
@@ -128,13 +128,13 @@ def SolverWrapI(fun, checkAccuracy=True, accuracyTol=1e-5, name=None):
def clean(self):
pass
return type(name if name is not None else fun.__name__, (object,), {"__init__": __init__, "clean": clean, "__mul__": __mul__})
return type(fun.__name__+'_Wrapped', (object,), {"__init__": __init__, "clean": clean, "__mul__": __mul__})
from scipy.sparse import linalg
Solver = SolverWrapD(linalg.spsolve, factorize=False, name="Solver")
SolverLU = SolverWrapD(linalg.splu, factorize=True, name="SolverLU")
SolverCG = SolverWrapI(linalg.cg, name="SolverCG")
Solver = SolverWrapD(linalg.spsolve, factorize=False)
SolverLU = SolverWrapD(linalg.splu, factorize=True)
SolverCG = SolverWrapI(linalg.cg)
class SolverDiag(object):
+1 -1
View File
@@ -25,7 +25,7 @@ def interpmat(locs, x, y=None, z=None):
:param numpy.ndarray x: Tensor vector of 1st dimension of grid.
:param numpy.ndarray y: Tensor vector of 2nd dimension of grid. None by default.
:param numpy.ndarray z: Tensor vector of 3rd dimension of grid. None by default.
:rtype: scipy.sparse.csr_matrix
:rtype: scipy.sparse.csr.csr_matrix
:return: Interpolation matrix
.. plot::
+6 -6
View File
@@ -27,7 +27,7 @@ def mkvc(x, numDims=1):
if isinstance(x, Zero):
return x
assert isinstance(x, np.ndarray), "Vector must be a numpy array"
if numDims == 1:
@@ -355,9 +355,9 @@ def diagEst(matFun, n, k=None, approach='Probing'):
2. Ones : random +/- 1 entries
3. Random : random vectors
:param callable matFun: takes a (numpy.array) and multiplies it by a matrix to estimate the diagonal
:param int n: size of the vector that should be used to compute matFun(v)
:param int k: number of vectors to be used to estimate the diagonal
:param lambda (numpy.array) matFun: matrix to estimate the diagonal of
:param int64 n: size of the vector that should be used to compute matFun(v)
:param int64 k: number of vectors to be used to estimate the diagonal
:param str approach: approach to be used for getting vectors
:rtype: numpy.array
:return: est_diag(A)
@@ -422,9 +422,9 @@ class Zero(object):
def __ge__(self, v):return 0 >= v
def __gt__(self, v):return 0 > v
@property
@property
def transpose(self): return Zero()
@property
def T(self): return Zero()
+14 -18
View File
@@ -83,7 +83,7 @@ def closestPoints(mesh, pts, gridLoc='CC'):
"""
Move a list of points to the closest points on a grid.
:param BaseMesh mesh: The mesh
:param simpeg.Mesh.BaseMesh mesh: The mesh
:param numpy.ndarray pts: Points to move
:param string gridLoc: ['CC', 'N', 'Fx', 'Fy', 'Fz', 'Ex', 'Ex', 'Ey', 'Ez']
:rtype: numpy.ndarray
@@ -104,20 +104,16 @@ def closestPoints(mesh, pts, gridLoc='CC'):
def ExtractCoreMesh(xyzlim, mesh, meshType='tensor'):
"""
Extracts Core Mesh from Global mesh
:param numpy.ndarray xyzlim: 2D array [ndim x 2]
:param BaseMesh mesh: The mesh
This function ouputs::
- actind: corresponding boolean index from global to core
- meshcore: core SimPEG mesh
Warning: 1D and 2D has not been tested
Extracts Core Mesh from Global mesh
xyzlim: 2D array [ndim x 2]
mesh: SimPEG mesh
This function ouputs:
- actind: corresponding boolean index from global to core
- meshcore: core SimPEG mesh
Warning: 1D and 2D has not been tested
"""
from SimPEG import Mesh
if mesh.dim == 1:
if mesh.dim ==1:
xyzlim = xyzlim.flatten()
xmin, xmax = xyzlim[0], xyzlim[1]
@@ -129,11 +125,11 @@ def ExtractCoreMesh(xyzlim, mesh, meshType='tensor'):
x0 = [xc[0]-hx[0]*0.5, yc[0]-hy[0]*0.5]
meshCore = Mesh.TensorMesh([hx, hy], x0=x0)
meshCore = Mesh.TensorMesh([hx, hy] ,x0=x0)
actind = (mesh.gridCC[:,0]>xmin) & (mesh.gridCC[:,0]<xmax)
elif mesh.dim == 2:
elif mesh.dim ==2:
xmin, xmax = xyzlim[0,0], xyzlim[0,1]
ymin, ymax = xyzlim[1,0], xyzlim[1,1]
@@ -148,12 +144,12 @@ def ExtractCoreMesh(xyzlim, mesh, meshType='tensor'):
x0 = [xc[0]-hx[0]*0.5, yc[0]-hy[0]*0.5]
meshCore = Mesh.TensorMesh([hx, hy], x0=x0)
meshCore = Mesh.TensorMesh([hx, hy] ,x0=x0)
actind = (mesh.gridCC[:,0]>xmin) & (mesh.gridCC[:,0]<xmax) \
& (mesh.gridCC[:,1]>ymin) & (mesh.gridCC[:,1]<ymax) \
elif mesh.dim == 3:
elif mesh.dim==3:
xmin, xmax = xyzlim[0,0], xyzlim[0,1]
ymin, ymax = xyzlim[1,0], xyzlim[1,1]
zmin, zmax = xyzlim[2,0], xyzlim[2,1]
@@ -172,7 +168,7 @@ def ExtractCoreMesh(xyzlim, mesh, meshType='tensor'):
x0 = [xc[0]-hx[0]*0.5, yc[0]-hy[0]*0.5, zc[0]-hz[0]*0.5]
meshCore = Mesh.TensorMesh([hx, hy, hz], x0=x0)
meshCore = Mesh.TensorMesh([hx, hy, hz] ,x0=x0)
actind = (mesh.gridCC[:,0]>xmin) & (mesh.gridCC[:,0]<xmax) \
& (mesh.gridCC[:,1]>ymin) & (mesh.gridCC[:,1]<ymax) \
+1 -1
View File
@@ -15,7 +15,7 @@ import Directives
import Inversion
import Tests
__version__ = '0.1.12'
__version__ = '0.1.10'
__author__ = 'Rowan Cockett'
__license__ = 'MIT'
__copyright__ = 'Copyright 2014 Rowan Cockett'

Before

Width:  |  Height:  |  Size: 49 KiB

After

Width:  |  Height:  |  Size: 49 KiB

Before

Width:  |  Height:  |  Size: 58 KiB

After

Width:  |  Height:  |  Size: 58 KiB

+1 -1
View File
@@ -2,7 +2,7 @@
#
# You can set these variables from the command line.
SPHINXOPTS = -n -w warnings.txt
SPHINXOPTS =
SPHINXBUILD = sphinx-build
PAPER =
BUILDDIR = _build

Before

Width:  |  Height:  |  Size: 30 KiB

After

Width:  |  Height:  |  Size: 30 KiB

-22
View File
@@ -1,22 +0,0 @@
{# Import the theme's layout. #}
{% extends "!layout.html" %}
{% block extrahead %}
{{ super() }}
<meta name="description" content="Simulation and Parameter Estimation in Geophysics">
<meta name="author" content="SimPEG Developers">
<meta name="keywords" content="python, geophysics, inversion, electromagnetics, magnetotellurics, magnetics, gravity, DC, flow inverse problems, open source, finite volume">
<script>
(function(i,s,o,g,r,a,m){i['GoogleAnalyticsObject']=r;i[r]=i[r]||function(){
(i[r].q=i[r].q||[]).push(arguments)},i[r].l=1*new Date();a=s.createElement(o),
m=s.getElementsByTagName(o)[0];a.async=1;a.src=g;m.parentNode.insertBefore(a,m)
})(window,document,'script','https://www.google-analytics.com/analytics.js','ga');
ga('create', 'UA-45185336-1', 'auto');
ga('send', 'pageview');
</script>
{% endblock %}
+9 -16
View File
@@ -1,3 +1,5 @@
.. _api_DC:
.. math::
\renewcommand{\div}{\nabla\cdot\,}
@@ -36,16 +38,8 @@
\renewcommand {\u} { {\vec u} }
\newcommand{\I}{\vec{I}}
Direct Current Resistivity
**************************
`SimPEG.DCIP` uses SimPEG as the framework for the forward and inverse
direct current (DC) resistivity and induced polarization (IP) geophysical problems.
DC resistivity survey
=====================
*********************
Electrical resistivity of subsurface materials is measured by causing an electrical current to flow in the earth between one pair of electrodes while the voltage across a second pair of electrodes is measured. The result is an "apparent" resistivity which is a value representing the weighted average resistivity over a volume of the earth. Variations in this measurement are caused by variations in the soil, rock, and pore fluid electrical resistivity. Surveys require contact with the ground, so they can be labour intensive. Results are sometimes interpreted directly, but more commonly, 1D, 2D or 3D models are estimated using inversion procedures (`GPG <http://www.eos.ubc.ca/courses/eosc350/content/>`_).
@@ -61,7 +55,7 @@ As direct current (DC) implies, in DC resistivity survey, we assume steady-state
\curl \e = 0
Then by taking \\(\\div\\) of the first equation, we have
Then by taking \\(\\curl\\) for the first equation, we have
.. math::
@@ -143,14 +137,13 @@ Comparing to the analytic function:
.. plot::
from SimPEG import Examples
Examples.DC_Analytic_Dipole.run(plotIt=True)
import simpegDC as DC
DC.Examples.Verification.run(plotIt=True)
API
===
API for DC codes
================
.. automodule:: SimPEG.DCIP.BaseDC
.. automodule:: simpegDC.BaseDC
:show-inheritance:
:members:
:undoc-members:
@@ -7,7 +7,7 @@ Examples
:maxdepth: 1
:glob:
../examples/*
examples/*
External Notebooks
+19
View File
@@ -0,0 +1,19 @@
.. _api_FiniteVolume:
Finite Volume
*************
Any numerical implementation requires the discretization of continuous functions into discrete approximations. These approximations are typically organized in a mesh, which defines boundaries, locations, and connectivity. Of specific interest to geophysical simulations, we require that averaging, interpolation and differential operators be defined for any mesh. In SimPEG, we have implemented a staggered mimetic finite volume approach (`Hyman and Shashkov, 1999 <http://math.lanl.gov/~mac/papers/numerics/HS99B.pdf>`_). This approach requires the definitions of variables at either cell-centers, nodes, faces, or edges as seen in the figure below.
.. image:: images/finitevolrealestate.png
:width: 400 px
:alt: FiniteVolume
:align: center
.. toctree::
:maxdepth: 2
api_Mesh
api_DiffOps
api_InnerProducts
@@ -52,15 +52,13 @@ We can take the derivative of the PDE:
\nabla_m c(m, u) \partial m + \nabla_u c(m, u) \partial u = 0
If the forward problem is invertible, then we can rearrange for
\\(\\frac{\\partial u}{\\partial m}\\):
If the forward problem is invertible, then we can rearrange for \\(\\frac{\\partial u}{\\partial m}\\):
.. math::
J = - P \left( \nabla_u c(m, u) \right)^{-1} \nabla_m c(m, u)
This can often be computed given a vector (i.e. \\(J(v)\\)) rather than
stored, as \\(J\\) is a large dense matrix.
This can often be computed given a vector (i.e. \\(J(v)\\)) rather than stored, as \\(J\\) is a large dense matrix.
@@ -69,45 +67,13 @@ The API
Problem
-------
.. autoclass:: SimPEG.Problem.BaseProblem
:members:
:undoc-members:
.. autoclass:: SimPEG.Problem.BaseTimeProblem
:members:
:undoc-members:
Fields
------
.. autoclass:: SimPEG.Fields.Fields
:members:
:undoc-members:
.. autoclass:: SimPEG.Fields.TimeFields
.. automodule:: SimPEG.Problem
:members:
:undoc-members:
Survey
------
.. autoclass:: SimPEG.Survey.BaseSurvey
.. automodule:: SimPEG.Survey
:members:
:undoc-members:
.. autoclass:: SimPEG.Survey.BaseSrc
:members:
:undoc-members:
.. autoclass:: SimPEG.Survey.BaseRx
:members:
:undoc-members:
.. autoclass:: SimPEG.Survey.BaseTimeRx
:members:
:undoc-members:
.. autoclass:: SimPEG.Survey.Data
:members:
:undoc-members:
@@ -4,10 +4,7 @@
Inner Products
**************
By using the weak formulation of many of the PDEs in geophysical applications,
we can rapidly develop discretizations. Much of this work, however, needs a
good understanding of how to approximate inner products on our discretized
meshes. We will define the inner product as:
By using the weak formulation of many of the PDEs in geophysical applications, we can rapidly develop discretizations. Much of this work, however, needs a good understanding of how to approximate inner products on our discretized meshes. We will define the inner product as:
.. math::
@@ -17,15 +14,12 @@ where a and b are either scalars or vectors.
.. note::
The InnerProducts class is a base class providing inner product matrices
for meshes and cannot run on its own.
The InnerProducts class is a base class providing inner product matrices for meshes and cannot run on its own.
Example problem for DC resistivity
----------------------------------
We will start with the formulation of the Direct Current (DC) resistivity
problem in geophysics.
We will start with the formulation of the Direct Current (DC) resistivity problem in geophysics.
.. math::
@@ -34,13 +28,12 @@ problem in geophysics.
\nabla\cdot \vec{j} = q
In the following discretization, :math:`\sigma` and :math:`\phi`
will be discretized on the cell-centers and the flux, :math:`\vec{j}`,
In the following discretization, \\\( \\sigma \\\) and \\\( \\phi \\\)
will be discretized on the cell-centers and the flux, \\\(\\vec{j}\\\),
will be on the faces. We will use the weak formulation to discretize
the DC resistivity equation.
We can define in weak form by integrating with a general face function
:math:`\vec{f}`:
We can define in weak form by integrating with a general face function \\\(\\vec{f}\\\):
.. math::
@@ -68,16 +61,9 @@ We can then discretize for every cell:
.. note::
We have discretized the dot product above, but remember that we do not
really have a single vector :math:`\mathbf{J}`, but approximations of
:math:`\vec{j}` on each face of our cell. In 2D that means 2
approximations of :math:`\mathbf{J}_x` and 2 approximations of
:math:`\mathbf{J}_y`. In 3D we also have 2 approximations of
:math:`\mathbf{J}_z`.
We have discretized the dot product above, but remember that we do not really have a single vector \\\(\\mathbf{J}\\\), but approximations of \\\(\\vec{j}\\\) on each face of our cell. In 2D that means 2 approximations of \\\(\\mathbf{J}_x\\\) and 2 approximations of \\\(\\mathbf{J}_y\\\). In 3D we also have 2 approximations of \\\(\\mathbf{J}_z\\\).
Regardless of how we choose to approximate this dot product, we can represent
this in vector form (again this is for every cell), and will generalize for
the case of anisotropic (tensor) sigma.
Regardless of how we choose to approximate this dot product, we can represent this in vector form (again this is for every cell), and will generalize for the case of anisotropic (tensor) sigma.
.. math::
@@ -85,17 +71,14 @@ the case of anisotropic (tensor) sigma.
-\phi^{\top} v_{\text{cell}} \mathbf{D}_{\text{cell}} \mathbf{F})
+ \text{BC}
We multiply by square-root of volume on each side of the tensor conductivity
to keep symmetry in the system. Here :math:`\mathbf{J}_c` is the Cartesian
:math:`\mathbf{J}` (on the faces that we choose to use in our approximation)
and must be calculated differently depending on the mesh:
We multiply by square-root of volume on each side of the tensor conductivity to keep symmetry in the system. Here \\\(\\mathbf{J}_c\\\) is the Cartesian \\\(\\mathbf{J}\\\) (on the faces that we choose to use in our approximation) and must be calculated differently depending on the mesh:
.. math::
\mathbf{J}_c = \mathbf{Q}_{(i)}\mathbf{J}_\text{TENSOR} \\
\mathbf{J}_c = \mathbf{N}_{(i)}^{-1}\mathbf{Q}_{(i)}\mathbf{J}_\text{Curv}
Here the :math:`i` index refers to where we choose to approximate this integral, as discussed in the note above.
We will approximate this integral by taking the fluxes clustered around every node of the cell, there are 8 combinations in 3D, and 4 in 2D. We will use a projection matrix :math:`\mathbf{Q}_{(i)}` to pick the appropriate fluxes. So, now that we have 8 approximations of this integral, we will just take the average. For the TensorMesh, this looks like:
Here the \\\(i\\\) index refers to where we choose to approximate this integral, as discussed in the note above.
We will approximate this integral by taking the fluxes clustered around every node of the cell, there are 8 combinations in 3D, and 4 in 2D. We will use a projection matrix \\\( \\mathbf{Q}_{(i)} \\\) to pick the appropriate fluxes. So, now that we have 8 approximations of this integral, we will just take the average. For the TensorMesh, this looks like:
.. math::
@@ -124,12 +107,10 @@ By defining the faceInnerProduct (8 combinations of fluxes in 3D, 4 in 2D, 2 in
\sum_{i=1}^{2^d}
\mathbf{P}_{(i)}^{\top} \Sigma^{-1} \mathbf{P}_{(i)}
Where :math:`d` is the dimension of the mesh.
The :math:`\mathbf{M}^f` is returned when given the input of :math:`\Sigma^{-1}`.
Where \\\(d\\\) is the dimension of the mesh.
The \\\( \\mathbf{M}^f \\\) is returned when given the input of \\\( \\Sigma^{-1} \\\).
Here each :math:`\mathbf{P} ~ \in ~ \mathbb{R}^{(d*nC, nF)}` is a combination
of the projection, volume, and any normalization to Cartesian coordinates
(where the dot product is well defined):
Here each \\( \\mathbf{P} \\in \\mathbb{R}^{(d*nC, nF)} \\\) is a combination of the projection, volume, and any normalization to Cartesian coordinates (where the dot product is well defined):
.. math::
@@ -148,10 +129,7 @@ If ``returnP=True`` is requested in any of these methods the projection matrices
# In 1D
P = [P0, P1]
The derivation for ``edgeInnerProducts`` is exactly the same, however, when we
approximate the integral using the fields around each node, the projection
matrices look a bit different because we have 12 edges in 3D instead of just 6
faces. The interface to the code is exactly the same.
The derivation for ``edgeInnerProducts`` is exactly the same, however, when we approximate the integral using the fields around each node, the projection matrices look a bit different because we have 12 edges in 3D instead of just 6 faces. The interface to the code is exactly the same.
Defining Tensor Properties
@@ -159,8 +137,7 @@ Defining Tensor Properties
**For 3D:**
Depending on the number of columns (either 1, 3, or 6) of mu, the material
property is interpreted as follows:
Depending on the number of columns (either 1, 3, or 6) of mu, the material property is interpreted as follows:
.. math::
@@ -211,16 +188,13 @@ Which is nice and easy to invert if necessary, however, in the fully anisotropic
Taking Derivatives
------------------
We will take the derivative of the fully anisotropic tensor for a 3D mesh, the
other cases are easier and will not be discussed here. Let us start with one
part of the sum which makes up :math:`\mathbf{M}^f_\Sigma` and take the
derivative when this is multiplied by some vector :math:`\mathbf{v}`:
We will take the derivative of the fully anisotropic tensor for a 3D mesh, the other cases are easier and will not be discussed here. Let us start with one part of the sum which makes up \\\(\\mathbf{M}^f_\\Sigma\\\) and take the derivative when this is multiplied by some vector \\\(\\mathbf{v}\\\):
.. math::
\mathbf{P}^\top \boldsymbol{\Sigma} \mathbf{Pv}
Here we will let :math:`\mathbf{Pv} = \mathbf{y}` and :math:`\mathbf{y}` will have the form:
Here we will let \\\( \\mathbf{Pv} = \\mathbf{y} \\\) and \\\(\\mathbf{y}\\\) will have the form:
.. math::
@@ -259,9 +233,7 @@ Here we will let :math:`\mathbf{Pv} = \mathbf{y}` and :math:`\mathbf{y}` will ha
\end{matrix}
\right]
Now it is easy to take the derivative with respect to any one of the
parameters, for example,
:math:`\frac{\partial}{\partial\boldsymbol{\sigma}_1}`
Now it is easy to take the derivative with respect to any one of the parameters, for example, \\\(\\frac{\\partial}{\\partial\\boldsymbol{\\sigma}_1}\\\)
.. math::
\frac{\partial}{\partial \boldsymbol{\sigma}_1}\left(\mathbf{P}^\top\Sigma\mathbf{y}\right)
@@ -275,8 +247,7 @@ parameters, for example,
\end{matrix}
\right]
Whereas :math:`\frac{\partial}{\partial\boldsymbol{\sigma}_4}`, for
example, is:
Whereas \\\(\\frac{\\partial}{\\partial\\boldsymbol{\\sigma}_4}\\\), for example, is:
.. math::
\frac{\partial}{\partial \boldsymbol{\sigma}_4}\left(\mathbf{P}^\top\Sigma\mathbf{y}\right)
@@ -290,12 +261,11 @@ example, is:
\end{matrix}
\right]
These are computed for each of the 8 projections, horizontally concatenated,
and returned.
These are computed for each of the 8 projections, horizontally concatenated, and returned.
The API
-------
.. autoclass:: SimPEG.Mesh.InnerProducts.InnerProducts
.. automodule:: SimPEG.Mesh.InnerProducts
:members:
:undoc-members:
@@ -3,7 +3,7 @@
InvProblem
**********
.. autoclass:: SimPEG.InvProblem.BaseInvProblem
.. automodule:: SimPEG.InvProblem
:show-inheritance:
:members:
:undoc-members:
@@ -12,7 +12,7 @@ InvProblem
Inversion
*********
.. autoclass:: SimPEG.Inversion.BaseInversion
.. automodule:: SimPEG.Inversion
:show-inheritance:
:members:
:undoc-members:
@@ -27,8 +27,7 @@ back to conductivity. This is a relatively trivial example (we are just taking
the exponential!) but by defining maps we can start to combine and manipulate
exactly what we think about as our model, \\\(m\\\). In code, this looks like
.. code-block:: python
:linenos:
::
M = Mesh.TensorMesh([100]) # Create a mesh
expMap = Maps.ExpMap(M) # Create a mapping
@@ -47,15 +46,14 @@ We will use an example where we want a 1D layered earth as
our model, but we want to map this to a 2D discretization to do our forward
modeling. We will also assume that we are working in log conductivity still,
so after the transformation we want to map to conductivity space.
To do this we will introduce the vertical 1D map (:class:`SimPEG.Maps.SurjectVertical1D`),
To do this we will introduce the vertical 1D map (:class:`SimPEG.Maps.Vertical1DMap`),
which does the first part of what we just described. The second part will be
done by the :class:`SimPEG.Maps.ExpMap` described above.
.. code-block:: python
:linenos:
::
M = Mesh.TensorMesh([7,5])
v1dMap = Maps.SurjectVertical1D(M)
v1dMap = Maps.Vertical1DMap(M)
expMap = Maps.ExpMap(M)
myMap = expMap * v1dMap
m = np.r_[0.2,1,0.1,2,2.9] # only 5 model parameters!
@@ -63,8 +61,26 @@ done by the :class:`SimPEG.Maps.ExpMap` described above.
.. plot::
from SimPEG import Examples
Examples.Maps_ComboMaps.run()
from SimPEG import *
import matplotlib.pyplot as plt
M = Mesh.TensorMesh([7,5])
v1dMap = Maps.Vertical1DMap(M)
expMap = Maps.ExpMap(M)
myMap = expMap * v1dMap
m = np.r_[0.2,1,0.1,2,2.9] # only 5 model parameters!
sig = myMap * m
figs, axs = plt.subplots(1,2)
axs[0].plot(m, M.vectorCCy, 'b-o')
axs[0].set_title('Model')
axs[0].set_ylabel('Depth, y')
axs[0].set_xlabel('Value, $m_i$')
axs[0].set_xlim(0,3)
axs[0].set_ylim(0,1)
clbar = plt.colorbar(M.plotImage(sig,ax=axs[1],grid=True,gridOpts=dict(color='grey'))[0])
axs[1].set_title('Physical Property')
axs[1].set_ylabel('Depth, y')
clbar.set_label('$\sigma = \exp(\mathbf{P}m)$')
plt.tight_layout()
If you noticed, it was pretty easy to combine maps. What is even cooler is
that the derivatives also are made for you (if everything goes right).
@@ -106,8 +122,6 @@ When these are used in the inverse problem, this is extremely important!!
The API
=======
The :code:`IdentityMap` is the base class for all mappings, and it does absolutely nothing.
.. autoclass:: SimPEG.Maps.IdentityMap
:members:
:undoc-members:
@@ -116,6 +130,7 @@ The :code:`IdentityMap` is the base class for all mappings, and it does absolute
Common Maps
===========
Exponential Map
---------------
@@ -133,7 +148,7 @@ lives (i.e. it varies logarithmically).
Vertical 1D Map
---------------
.. autoclass:: SimPEG.Maps.SurjectVertical1D
.. autoclass:: SimPEG.Maps.Vertical1DMap
:members:
:undoc-members:
@@ -149,10 +164,31 @@ Map 2D Cross-Section to 3D Model
Mesh to Mesh Map
----------------
.. plot::
from SimPEG import Examples
Examples.Maps_Mesh2Mesh.run()
from SimPEG import *
import matplotlib.pyplot as plt
M = Mesh.TensorMesh([100,100])
h1 = Utils.meshTensor([(6,7,-1.5),(6,10),(6,7,1.5)])
h1 = h1/h1.sum()
M2 = Mesh.TensorMesh([h1,h1])
V = Utils.ModelBuilder.randomModel(M.vnC, seed=79, its=50)
v = Utils.mkvc(V)
modh = Maps.Mesh2Mesh([M,M2])
modH = Maps.Mesh2Mesh([M2,M])
H = modH * v
h = modh * H
ax = plt.subplot(131)
M.plotImage(v, ax=ax)
ax.set_title('Fine Mesh (Original)')
ax = plt.subplot(132)
M2.plotImage(H,clim=[0,1],ax=ax)
ax.set_title('Course Mesh')
ax = plt.subplot(133)
M.plotImage(h,clim=[0,1],ax=ax)
ax.set_title('Fine Mesh (Interpolated)')
plt.show()
.. autoclass:: SimPEG.Maps.Mesh2Mesh
@@ -160,8 +196,8 @@ Mesh to Mesh Map
:undoc-members:
Under the Hood
==============
Some Extras
===========
Combo Map
---------
@@ -188,6 +188,6 @@ other types of meshes in this SimPEG framework.
The API
=======
.. autoclass:: SimPEG.Mesh.BaseMesh.BaseMesh
.. automodule:: SimPEG.Mesh.BaseMesh
:members:
:undoc-members:
+36
View File
@@ -0,0 +1,36 @@
.. _api_MeshCode:
Tensor Mesh
===========
.. automodule:: SimPEG.Mesh.TensorMesh
:show-inheritance:
:members:
:undoc-members:
Cylindrical Mesh
================
.. automodule:: SimPEG.Mesh.CylMesh
:show-inheritance:
:members:
:undoc-members:
Tree Mesh
=========
.. autoclass:: SimPEG.Mesh.TreeMesh.TreeMesh
:show-inheritance:
:members:
:undoc-members:
Curvilinear Mesh
================
.. automodule:: SimPEG.Mesh.CurvilinearMesh
:show-inheritance:
:members:
:undoc-members:
@@ -91,21 +91,10 @@ The API
:members:
:undoc-members:
.. autoclass:: SimPEG.Regularization.Simple
:show-inheritance:
:members:
.. autoclass:: SimPEG.Regularization.Tikhonov
:show-inheritance:
:members:
.. autoclass:: SimPEG.Regularization.Sparse
:show-inheritance:
:members:
.. autoclass:: SimPEG.Regularization.RegularizationMesh
:show-inheritance:
:members:
@@ -46,8 +46,6 @@ The API
=======
.. autofunction:: SimPEG.Utils.SolverUtils.SolverWrapD
:noindex:
.. autofunction:: SimPEG.Utils.SolverUtils.SolverWrapI
:noindex:
@@ -6,6 +6,5 @@ Utilities
api_Solver
api_Maps
api_PropMaps
api_Utils
api_Tests
@@ -21,7 +21,7 @@ Solver Utilities
:undoc-members:
Curv Utilities
==============
=============
.. automodule:: SimPEG.Utils.curvutils
:members:
@@ -51,9 +51,7 @@ Interpolation Utilities
Counter Utilities
=================
.. code-block:: python
:linenos:
::
class MyClass(object):
def __init__(self, url):
self.counter = Counter()
@@ -71,9 +69,7 @@ Counter Utilities
for i in range(300): c.MySecondMethod()
c.counter.summary()
.. code-block:: text
:linenos:
::
Counters:
MyClass.MyMethod : 100
@@ -81,8 +77,6 @@ Counter Utilities
Times: mean sum
MyClass.MySecondMethod : 1.70e-06, 5.10e-04, 300x
The API
-------
@@ -35,7 +35,7 @@ The Big Picture
Defining a well-posed inverse problem and solving it is a complex task that requires many components that must interact. It is helpful
to view this task as a workflow in which various elements are explicitly identified and integrated. The figure below outlines the inversion components that consists of inputs, implementation, and evaluation. The inputs are composed of the geophysical data, the equations which are a mathematical description of the governing physics, and prior knowledge or assumptions about the setting. The implementation consists of two broad categories: the forward simulation and the inversion. The **forward simulation** is the means by which we solve the governing equations given a model and the **inversion components** evaluate and update this model. We are considering a gradient based approach, which updates the model through an optimization routine. The output of this implementation is a model, which, prior to interpretation, must be evaluated. This requires considering, and often re-assessing, the choices and assumptions made in both the input and implementation stages.
.. image:: ../../images/InversionWorkflow-PreSimPEG.png
.. image:: InversionWorkflow-PreSimPEG.png
:width: 400 px
:alt: Components
:align: center
@@ -46,24 +46,24 @@ A Comprehensive Framework
There are an overwhelming amount of choices to be made as one works through the forward modeling and inversion process (see figure above). As a result, software implementations of this workflow often become complex and highly interdependent, making it difficult to interact with and to ask other scientists to pick up and change. Our approach to handling this complexity is to propose a framework, (see below), that compartmentalizes the implementation of inversions into various units. We present it in this specific modular style, as each unit contains a targeted subset of choices crucial to the inversion process.
.. image:: ../../images/InversionWorkflow.png
.. image:: InversionWorkflow.png
:width: 400 px
:alt: Framework
:align: center
The process of obtaining an acceptable model from an inversion generally requires the geophysicist to perform several iterations of the inversion workflow, rethinking and redesigning each piece of the framework to ensure it is appropriate in the current context. Inversions are experimental and empirical by nature and our software package is designed to facilitate this iterative process. To accomplish this, we have divided the inversion methodology into eight major components (See figure above). The :class:`SimPEG.Mesh.BaseMesh.BaseMesh` class handles the discretization of the earth and also provides numerical operators. The forward simulation is split into two classes, the :class:`SimPEG.Survey.BaseSurvey` and the :class:`SimPEG.Problem.BaseProblem`. The :class:`SimPEG.Survey.BaseSurvey` class handles the geometry of a geophysical problem as well as sources. The :class:`SimPEG.Problem.BaseProblem` class handles the simulation of the physics for the geophysical problem of interest. Although created independently, these two classes must be paired to form all of the components necessary for a geophysical forward simulation and calculation of the sensitivity. The :class:`SimPEG.Problem.BaseProblem` creates geophysical fields given a source from the :class:`SimPEG.Survey.BaseSurvey`. The :class:`SimPEG.Survey.BaseSurvey` interpolates these fields to the receiver locations and converts them to the appropriate data type, for example, by selecting only the measured components of the field. Each of these operations may have associated derivatives with respect to the model and the computed field; these are included in the calculation of the sensitivity. For the inversion, a :class:`SimPEG.DataMisfit.BaseDataMisfit` is chosen to capture the goodness of fit of the predicted data and a :class:`SimPEG.Regularization.BaseRegularization` is chosen to handle the non-uniqueness. These inversion elements and an Optimization routine are combined into an inverse problem class :class:`SimPEG.InvProblem.BaseInvProblem`. :class:`SimPEG.InvProblem.BaseInvProblem` is the mathematical statement that will be numerically solved by running an Inversion. The :class:`SimPEG.Inversion.BaseInversion` class handles organization and dispatch of directives between all of the various pieces of the framework.
The process of obtaining an acceptable model from an inversion generally requires the geophysicist to perform several iterations of the inversion workflow, rethinking and redesigning each piece of the framework to ensure it is appropriate in the current context. Inversions are experimental and empirical by nature and our software package is designed to facilitate this iterative process. To accomplish this, we have divided the inversion methodology into eight major components (See figure above). The (:class:`SimPEG.Mesh.BaseMesh`) class handles the discretization of the earth and also provides numerical operators. The forward simulation is split into two classes, the (:class:`SimPEG.Survey.BaseSurvey`) and the (:class:`SimPEG.Problem.BaseProblem`). The (:class:`SimPEG.Survey.BaseSurvey`) class handles the geometry of a geophysical problem as well as sources. The (:class:`SimPEG.Problem.BaseProblem`) class handles the simulation of the physics for the geophysical problem of interest. Although created independently, these two classes must be paired to form all of the components necessary for a geophysical forward simulation and calculation of the sensitivity. The (:class:`SimPEG.Problem.BaseProblem`) creates geophysical fields given a source from the (:class:`SimPEG.Survey.BaseSurvey`). The (:class:`SimPEG.Survey.BaseSurvey`) interpolates these fields to the receiver locations and converts them to the appropriate data type, for example, by selecting only the measured components of the field. Each of these operations may have associated derivatives with respect to the model and the computed field; these are included in the calculation of the sensitivity. For the inversion, a (:class:`SimPEG.DataMisfit.BaseDataMisfit`) is chosen to capture the goodness of fit of the predicted data and a (:class:`SimPEG.Regularization.BaseRegularization`) is chosen to handle the non-uniqueness. These inversion elements and an Optimization routine are combined into an inverse problem class (:class:`SimPEG.InvProblem.BaseInvProblem`). (:class:`SimPEG.InvProblem.BaseInvProblem`) is the mathematical statement that will be numerically solved by running an Inversion. The (:class:`SimPEG.Inversion.BaseInversion`) class handles organization and dispatch of directives between all of the various pieces of the framework.
The arrows in the figure above indicate what each class takes as a primary argument. For example, both the :class:`SimPEG.Problem.BaseProblem` and :class:`SimPEG.Regularization.BaseRegularization` classes take a :class:`SimPEG.Mesh.BaseMesh.BaseMesh` class as an argument. The diagram does not show class inheritance, as each of the base classes outlined have many subtypes that can be interchanged. The :class:`SimPEG.Mesh.BaseMesh.BaseMesh` class, for example, could be a regular Cartesian mesh :class:`SimPEG.Mesh.TensorMesh` or a cylindrical coordinate mesh :class:`SimPEG.Mesh.CylMesh`, which have many properties in common. These common features, such as both meshes being created from tensor products, can be exploited through inheritance of base classes, and differences can be expressed through subtype polymorphism. Please look at the documentation here for more in-depth information.
The arrows in the figure above indicate what each class takes as a primary argument. For example, both the (:class:`SimPEG.Problem.BaseProblem`) and (:class:`SimPEG.Regularization.BaseRegularization`) classes take a (:class:`SimPEG.Mesh.BaseMesh`) class as an argument. The diagram does not show class inheritance, as each of the base classes outlined have many subtypes that can be interchanged. The (:class:`SimPEG.Mesh.BaseMesh`) class, for example, could be a regular Cartesian mesh (:class:`SimPEG.Mesh.TensorMesh`) or a cylindrical coordinate mesh (:class:`SimPEG.Mesh.CylMesh`), which have many properties in common. These common features, such as both meshes being created from tensor products, can be exploited through inheritance of base classes, and differences can be expressed through subtype polymorphism. Please look at the documentation here for more in-depth information.
.. include:: ../../../CITATION.rst
.. include:: ../CITATION.rst
Authors
-------
.. include:: ../../../AUTHORS.rst
.. include:: ../AUTHORS.rst
License
-------
.. include:: ../../../LICENSE
.. include:: ../LICENSE
-95
View File
@@ -1,95 +0,0 @@
# application: simpegdocs
# version: 1
runtime: python27
api_version: 1
threadsafe: yes
handlers:
# favicon
- url: /images/logo-block\.ico
static_files: /images/logo-block.ico
upload: /images/logo-block\.ico
# all css
- url: /(.*\.css)
mime_type: text/css
static_files: _build/html/\1
upload: _build/html/(.*\.css)
# webfonts
- url: /(.*\.(eot|svg|ttf|woff|woff2|otf))
static_files: _build/html/\1
upload: _build/html/(.*\.(eot|svg|ttf|woff|woff2|otf))
# javascript
- url: /(.*\.js)
mime_type: text/javascript
static_files: _build/html/\1
upload: _build/html/(.*\.js)
# plain text source
- url: /(.*\.txt)
mime_type: text/plain
static_files: _build/html/\1
upload: _build/html/(.*\.txt)
# images
- url: /_images/(.*\.(gif|png|jpg|ico))
static_files: _build/html/_images/\1
upload: _build/html/_images/(.*\.(gif|png|jpg|ico))
# redirect en/latest traffic
- url: /en/latest/(.*\.html)
script: simpegdocs.app
# raw html
- url: /(.*\.html)
mime_type: text/html
static_files: _build/html/\1
upload: _build/html/(.*\.html)
# serve index files
- url: /(.+)/
static_files: _build/html/\1/index.html
upload: _build/html/(.+)/index.html
- url: /(.+)
static_files: _build/html/\1/index.html
upload: _build/html/(.+)/index.html
- url: /
static_files: _build/html/index.html
upload: _build/html/index.html
- url: .*
script: simpegdocs.app
# Recommended file skipping declaration from the GAE tutorials
skip_files:
- ^(.*/)?app\.yaml
- ^(.*/)?app\.yml
- ^(.*/)?#.*#
- ^(.*/)?.*~
- ^(.*/)?.*\.py[co]
- ^(.*/)?.*/RCS/.*
- ^(.*/)?\..*
- ^(.*/)?tests$
- ^(.*/)?test$
- ^test/(.*/)?
- ^COPYING.LESSER
- ^README\..*
- \.gitignore
- ^\.git/.*
- \.*\.lint$
- ^(.*/)?.*\.doctree$
libraries:
- name: webapp2
version: "2.5.2"
- name: PIL
version: "1.1.7"
- name: numpy
version: "latest"
- name: jinja2
version: "latest"
+6 -45
View File
@@ -28,7 +28,7 @@ sys.path.append('../')
# Add any Sphinx extension module names here, as strings. They can be extensions
# coming with Sphinx (named 'sphinx.ext.*') or your custom ones.
extensions = ['sphinx.ext.todo', 'sphinx.ext.mathjax', 'sphinx.ext.viewcode', 'sphinx.ext.autodoc', 'sphinx.ext.intersphinx', 'matplotlib.sphinxext.plot_directive']
extensions = ['sphinx.ext.todo', 'sphinx.ext.mathjax', 'sphinx.ext.viewcode', 'sphinx.ext.autodoc', 'matplotlib.sphinxext.plot_directive']
# Add any paths that contain templates here, relative to this directory.
templates_path = ['_templates']
@@ -44,16 +44,16 @@ master_doc = 'index'
# General information about the project.
project = u'SimPEG'
copyright = u'2013 - 2016, SimPEG Developers'
copyright = u'2013, SimPEG Developers'
# The version info for the project you're documenting, acts as replacement for
# |version| and |release|, also used in various other places throughout the
# built documents.
#
# The short X.Y version.
version = '0.1.12'
version = '0.1.10'
# The full version, including alpha/beta/rc tags.
release = '0.1.12'
release = '0.1.10'
# The language for content autogenerated by Sphinx. Refer to documentation
# for a list of supported languages.
@@ -124,12 +124,12 @@ except Exception, e:
# The name of an image file (within the static path) to use as favicon of the
# docs. This file should be a Windows icon file (.ico) being 16x16 or 32x32
# pixels large.
html_favicon = './images/logo-block.ico'
#html_favicon = None
# Add any paths that contain custom static files (such as style sheets) here,
# relative to this directory. They are copied after the builtin static files,
# so a file named "default.css" will overwrite the builtin "default.css".
html_static_path = []
html_static_path = ['_static']
# If not '', a 'Last updated on:' timestamp is inserted at every page bottom,
# using the given strftime format.
@@ -229,12 +229,6 @@ man_pages = [
# If true, show URL addresses after external links.
#man_show_urls = False
# Intersphinx
intersphinx_mapping = {'python': ('http://docs.python.org/2', None),
'numpy': ('http://docs.scipy.org/doc/numpy/', None),
'scipy': ('http://docs.scipy.org/doc/scipy/reference/', None),
'matplotlib': ('http://matplotlib.sourceforge.net/', None)}
# -- Options for Texinfo output ------------------------------------------------
@@ -257,36 +251,3 @@ texinfo_documents = [
#texinfo_show_urls = 'footnote'
autodoc_member_order = 'bysource'
def supress_nonlocal_image_warn():
import sphinx.environment
sphinx.environment.BuildEnvironment.warn_node = _supress_nonlocal_image_warn
def _supress_nonlocal_image_warn(self, msg, node):
from docutils.utils import get_source_line
if not msg.startswith('nonlocal image URI found:'):
self._warnfunc(msg, '%s:%s' % get_source_line(node))
supress_nonlocal_image_warn()
nitpick_ignore = [
('py:class', 'IdentityMap'),
('py:class', 'BaseSurvey'),
('py:class', 'BaseSrc'),
('py:class', 'BaseRx'),
('py:class', 'Survey'),
('py:class', 'FieldsFDEM'),
('py:class', 'Fields3D_e'),
('py:class', 'Fields3D_b'),
('py:class', 'Fields3D_j'),
('py:class', 'Fields3D_h'),
('py:class', 'SurveyTDEM'),
('py:class', 'SrcTDEM'),
('py:class', 'EMPropMap'),
('py:class', 'Data'),
('py:class', 'SurveyDC'),
('py:class', 'BaseMTFields'),
('py:class', 'SolverLU'),
]
@@ -1,27 +0,0 @@
.. _api_FiniteVolume:
Finite Volume
*************
Any numerical implementation requires the discretization of continuous
functions into discrete approximations. These approximations are typically
organized in a mesh, which defines boundaries, locations, and connectivity. Of
specific interest to geophysical simulations, we require that averaging,
interpolation and differential operators be defined for any mesh. In SimPEG,
we have implemented a staggered mimetic finite volume approach (`Hyman and
Shashkov, 1999 <http://math.lanl.gov/~mac/papers/numerics/HS99B.pdf>`_). This
approach requires the definitions of variables at either cell-centers, nodes,
faces, or edges as seen in the figure below.
.. image:: ../../images/finitevolrealestate.png
:width: 400 px
:alt: FiniteVolume
:align: center
.. toctree::
:maxdepth: 2
api_Mesh
api_DiffOps
api_InnerProducts
-68
View File
@@ -1,68 +0,0 @@
.. _api_MeshCode:
Tensor Mesh
===========
.. autoclass:: SimPEG.Mesh.TensorMesh
:members:
:undoc-members:
:show-inheritance:
Cylindrical Mesh
================
.. autoclass:: SimPEG.Mesh.CylMesh
:members:
:undoc-members:
:show-inheritance:
Tree Mesh
=========
.. autoclass:: SimPEG.Mesh.TreeMesh
:members:
:undoc-members:
:show-inheritance:
Curvilinear Mesh
================
.. autoclass:: SimPEG.Mesh.CurvilinearMesh
:members:
:undoc-members:
:show-inheritance:
Base Rectangular Mesh
=====================
.. autoclass:: SimPEG.Mesh.BaseMesh.BaseRectangularMesh
:members:
:undoc-members:
:show-inheritance:
Base Tensor Mesh
================
.. autoclass:: SimPEG.Mesh.TensorMesh.BaseTensorMesh
:members:
:undoc-members:
:show-inheritance:
Mesh IO
=======
.. automodule:: SimPEG.Mesh.MeshIO
:members:
:undoc-members:
:show-inheritance:
Mesh Viewing
============
.. automodule:: SimPEG.Mesh.View
:members:
:undoc-members:
:show-inheritance:
-29
View File
@@ -1,29 +0,0 @@
SimPEG PropMaps
***************
The API
=======
Property
--------
.. autoclass:: SimPEG.PropMaps.Property
:members:
:undoc-members:
PropMap
-------
.. autoclass:: SimPEG.PropMaps.PropMap
:members:
:undoc-members:
PropModel
---------
.. autoclass:: SimPEG.PropMaps.PropModel
:members:
:undoc-members:
-33
View File
@@ -1,33 +0,0 @@
Overview of Electromagnetics in SimPEG
**************************************
The API
=======
Physical Properties
-------------------
.. autoclass:: SimPEG.EM.Base.EMPropMap
:show-inheritance:
:members:
:undoc-members:
Problem
-------
.. autoclass:: SimPEG.EM.Base.BaseEMProblem
:show-inheritance:
:members:
:undoc-members:
Survey
------
.. autoclass:: SimPEG.EM.Base.BaseEMSurvey
:show-inheritance:
:members:
:undoc-members:
-48
View File
@@ -1,48 +0,0 @@
.. _examples_Maps_ComboMaps:
.. --------------------------------- ..
.. ..
.. THIS FILE IS AUTO GENEREATED ..
.. ..
.. SimPEG/Examples/__init__.py ..
.. ..
.. --------------------------------- ..
Maps: ComboMaps
===============
We will use an example where we want a 1D layered earth as
our model, but we want to map this to a 2D discretization to do our forward
modeling. We will also assume that we are working in log conductivity still,
so after the transformation we want to map to conductivity space.
To do this we will introduce the vertical 1D map (:class:`SimPEG.Maps.SurjectVertical1D`),
which does the first part of what we just described. The second part will be
done by the :class:`SimPEG.Maps.ExpMap` described above.
.. code-block:: python
:linenos:
M = Mesh.TensorMesh([7,5])
v1dMap = Maps.SurjectVertical1D(M)
expMap = Maps.ExpMap(M)
myMap = expMap * v1dMap
m = np.r_[0.2,1,0.1,2,2.9] # only 5 model parameters!
sig = myMap * m
If you noticed, it was pretty easy to combine maps. What is even cooler is
that the derivatives also are made for you (if everything goes right).
Just to be sure that the derivative is correct, you should always run the test
on the mapping that you create.
.. plot::
from SimPEG import Examples
Examples.Maps_ComboMaps.run()
.. literalinclude:: ../../../SimPEG/Examples/Maps_ComboMaps.py
:language: python
:linenos:

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