Compare commits

..
29 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
173 changed files with 2465 additions and 117267 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 -27
View File
@@ -18,32 +18,24 @@ env:
- TEST_DIR="tests/mesh tests/base tests/utils"
- TEST_DIR=tests/em/fdem/inverse/derivs
- TEST_DIR=tests/em/tdem
- TEST_DIR=tests/pf
- TEST_DIR=tests/dcip
- TEST_DIR=tests/flow
- TEST_DIR=tests/mt
- 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
@@ -54,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,19 +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, '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, '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 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
@@ -320,9 +301,6 @@ class Update_IRLS(InversionDirective):
self.reg.l2model = self.invProb.curModel
self.reg.curModel = self.invProb.curModel
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)
if getattr(self, 'f_old', None) is None:
self.f_old = self.reg.eval(self.invProb.curModel)#self.invProb.evalFunction(self.invProb.curModel, return_g=False, return_H=False)
@@ -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()
+4 -5
View File
@@ -1,5 +1,4 @@
from SimPEG import Mesh, Utils, np
from SimPEG import *
def run(plotIt=True):
"""
@@ -9,15 +8,15 @@ def run(plotIt=True):
Here we show SimPEG used to create three different types of meshes.
"""
sz = [16, 16]
sz = [16,16]
tM = Mesh.TensorMesh(sz)
qM = Mesh.TreeMesh(sz)
qM.refine(lambda cell: 4 if np.sqrt(((np.r_[cell.center]-0.5)**2).sum()) < 0.4 else 3)
rM = Mesh.CurvilinearMesh(Utils.meshutils.exampleLrmGrid(sz, 'rotate'))
rM = Mesh.CurvilinearMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate'))
if plotIt:
import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 3, figsize=(14, 5))
fig, axes = plt.subplots(1,3,figsize=(14,5))
opts = {}
tM.plotGrid(ax=axes[0], **opts)
axes[0].set_title('TensorMesh')
-65
View File
@@ -1,65 +0,0 @@
from SimPEG import Mesh, np, PF
def run(plotIt=True):
"""
PF: Magnetics: Analytics
========================
Comparing the magnetics field in Vancouver to Seoul
"""
xr = np.linspace(-300, 300, 41)
yr = np.linspace(-300, 300, 41)
X, Y = np.meshgrid(xr, yr)
Z = np.ones((np.size(xr), np.size(yr)))*150
# Bz component in Korea
inckr = -8. + 3./60
deckr = 54. + 9./60
btotkr = 50898.6
Bokr = PF.MagAnalytics.IDTtoxyz(inckr, deckr, btotkr)
bx, by, bz = PF.MagAnalytics.MagSphereAnaFunA(
X, Y, Z, 100., 0., 0., 0., 0.01, Bokr, 'secondary'
)
Bzkr = np.reshape(bz, (np.size(xr), np.size(yr)), order='F')
# Bz component in Canada
incca = 16. + 49./60
decca = 70. + 19./60
btotca = 54692.1
Boca = PF.MagAnalytics.IDTtoxyz(incca, decca, btotca)
bx, by, bz = PF.MagAnalytics.MagSphereAnaFunA(
X, Y, Z, 100., 0., 0., 0., 0.01, Boca, 'secondary'
)
Bzca = np.reshape(bz, (np.size(xr), np.size(yr)), order='F')
if plotIt:
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
fig = plt.figure(figsize=(14, 5))
ax1 = plt.subplot(121)
dat1 = plt.imshow(Bzkr, extent=[min(xr), max(xr), min(yr), max(yr)])
divider = make_axes_locatable(ax1)
cax1 = divider.append_axes("right", size="5%", pad=0.05)
ax1.set_xlabel('East-West (m)')
ax1.set_ylabel('South-North (m)')
plt.colorbar(dat1, cax=cax1)
ax1.set_title('$B_z$ field at Seoul, South Korea')
ax2 = plt.subplot(122)
dat2 = plt.imshow(Bzca, extent=[min(xr), max(xr), min(yr), max(yr)])
divider = make_axes_locatable(ax2)
cax2 = divider.append_axes("right", size="5%", pad=0.05)
ax2.set_xlabel('East-West (m)')
ax2.set_ylabel('South-North (m)')
plt.colorbar(dat2, cax=cax2)
ax2.set_title('$B_z$ field at Vancouver, Canada')
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
+7 -9
View File
@@ -41,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)
@@ -86,7 +86,7 @@ class IdentityMap(object):
The derivative of the transformation.
:param numpy.array m: model
:rtype: scipy.sparse.csr_matrix
:rtype: scipy.csr_matrix
:return: derivative of transformed model
"""
@@ -216,7 +216,7 @@ class ExpMap(IdentityMap):
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.
@@ -366,7 +366,7 @@ class SurjectVertical1D(IdentityMap):
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()
@@ -427,7 +427,7 @@ class Surject2Dto3D(IdentityMap):
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,9 +502,7 @@ 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]
+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.'
-53
View File
@@ -1,53 +0,0 @@
from SimPEG import Maps, Survey, Utils, np, sp
from scipy.constants import mu_0
import re
class LinearSurvey(Survey.BaseSurvey):
"""Base Magnetics Survey"""
rxLoc = None #: receiver locations
rxType = None #: receiver type
def __init__(self, srcField, **kwargs):
self.srcField = srcField
Survey.BaseSurvey.__init__(self, **kwargs)
def eval(self, u):
return u
@property
def nD(self):
return self.prob.G.shape[0]
@property
def nRx(self):
return self.srcField.rxList[0].locs.shape[0]
# def setBackgroundField(self, SrcField):
# if getattr(self, 'B0', None) is None:
# self._B0 = SrcField.param[0] * dipazm_2_xyz(SrcField.param[1],SrcField.param[2])
# return self._B0
class SrcField(Survey.BaseSrc):
""" Define the inducing field """
param = None #: Inducing field param (Amp, Incl, Decl)
def __init__(self, rxList, **kwargs):
super(SrcField, self).__init__(rxList, **kwargs)
class RxObs(Survey.BaseRx):
"""A station location must have be located in 3-D"""
def __init__(self, locsXYZ, **kwargs):
locs = locsXYZ
assert locsXYZ.shape[1] == 3, 'locs must in 3-D (x,y,z).'
super(RxObs, self).__init__(locs, 'tmi', storeProjections=False, **kwargs)
@property
def nD(self):
"""Number of data in the receiver."""
return self.locs[0].shape[0]
-194
View File
@@ -1,194 +0,0 @@
from SimPEG import Maps, Survey, Utils, np, sp
from scipy.constants import mu_0
import re
class BaseMagSurvey(Survey.BaseSurvey):
"""Base Magnetics Survey"""
rxLoc = None #: receiver locations
rxType = None #: receiver type
def __init__(self, **kwargs):
Survey.BaseSurvey.__init__(self, **kwargs)
def setBackgroundField(self, Inc, Dec, Btot):
Bx = Btot*np.cos(Inc/180.*np.pi)*np.sin(Dec/180.*np.pi)
By = Btot*np.cos(Inc/180.*np.pi)*np.cos(Dec/180.*np.pi)
Bz = -Btot*np.sin(Inc/180.*np.pi)
self.B0 = np.r_[Bx, By, Bz]
@property
def Qfx(self):
if getattr(self, '_Qfx', None) is None:
self._Qfx = self.prob.mesh.getInterpolationMat(self.rxLoc, 'Fx')
return self._Qfx
@property
def Qfy(self):
if getattr(self, '_Qfy', None) is None:
self._Qfy = self.prob.mesh.getInterpolationMat(self.rxLoc, 'Fy')
return self._Qfy
@property
def Qfz(self):
if getattr(self, '_Qfz', None) is None:
self._Qfz = self.prob.mesh.getInterpolationMat(self.rxLoc, 'Fz')
return self._Qfz
def projectFields(self, u):
"""
This function projects the fields onto the data space.
Especially, here for we use total magnetic intensity (TMI) data,
which is common in practice.
First we project our B on to data location
.. math::
\mathbf{B}_{rec} = \mathbf{P} \mathbf{B}
then we take the dot product between B and b_0
.. math ::
\\text{TMI} = \\vec{B}_s \cdot \hat{B}_0
"""
# TODO: There can be some different tyes of data like |B| or B
bfx = self.Qfx*u['B']
bfy = self.Qfy*u['B']
bfz = self.Qfz*u['B']
# Generate unit vector
B0 = self.prob.survey.B0
Bot = np.sqrt(B0[0]**2+B0[1]**2+B0[2]**2)
box = B0[0]/Bot
boy = B0[1]/Bot
boz = B0[2]/Bot
# return bfx*box + bfx*boy + bfx*boz
return bfx*box + bfy*boy + bfz*boz
@Utils.count
def projectFieldsDeriv(self, B):
"""
This function projects the fields onto the data space.
.. math::
\\frac{\partial d_\\text{pred}}{\partial \mathbf{B}} = \mathbf{P}
Especially, this function is for TMI data type
"""
# Generate unit vector
B0 = self.prob.survey.B0
Bot = np.sqrt(B0[0]**2+B0[1]**2+B0[2]**2)
box = B0[0]/Bot
boy = B0[1]/Bot
boz = B0[2]/Bot
return self.Qfx*box+self.Qfy*boy+self.Qfz*boz
def projectFieldsAsVector(self, B):
bfx = self.Qfx*B
bfy = self.Qfy*B
bfz = self.Qfz*B
return np.r_[bfx, bfy, bfz]
class LinearSurvey(Survey.BaseSurvey):
"""Base Magnetics Survey"""
rxLoc = None #: receiver locations
rxType = None #: receiver type
def __init__(self, srcField, **kwargs):
self.srcField = srcField
Survey.BaseSurvey.__init__(self, **kwargs)
def eval(self, u):
return u
@property
def nD(self):
return self.prob.G.shape[0]
@property
def nRx(self):
return self.srcField.rxList[0].locs.shape[0]
# def setBackgroundField(self, SrcField):
# if getattr(self, 'B0', None) is None:
# self._B0 = SrcField.param[0] * dipazm_2_xyz(SrcField.param[1],SrcField.param[2])
# return self._B0
class SrcField(Survey.BaseSrc):
""" Define the inducing field """
param = None #: Inducing field param (Amp, Incl, Decl)
def __init__(self, rxList, **kwargs):
super(SrcField, self).__init__(rxList, **kwargs)
class RxObs(Survey.BaseRx):
"""A station location must have be located in 3-D"""
def __init__(self, locsXYZ, **kwargs):
locs = locsXYZ
assert locsXYZ.shape[1] == 3, 'locs must in 3-D (x,y,z).'
super(RxObs, self).__init__(locs, 'tmi', storeProjections=False, **kwargs)
@property
def nD(self):
"""Number of data in the receiver."""
return self.locs[0].shape[0]
class MagSurveyBx(object):
"""docstring for MagSurveyBx"""
def __init__(self, **kwargs):
Survey.BaseData.__init__(self, **kwargs)
def projectFields(self, B):
bfx = self.Qfx*B
return bfx
class BaseMagMap(Maps.IdentityMap):
"""BaseMagMap"""
def __init__(self, mesh, **kwargs):
Maps.IdentityMap.__init__(self, mesh)
def _transform(self, m):
return mu_0*(1 + m)
def deriv(self, m):
return mu_0*sp.identity(self.nP)
class WeightMap(Maps.IdentityMap):
"""Weighted Map for distributed parameters"""
def __init__(self, nP, weight, **kwargs):
Maps.IdentityMap.__init__(self, nP)
self.mesh = None
self.weight = weight
def _transform(self, m):
return m*self.weight
def deriv(self, m):
return Utils.sdiag(self.weight)
-576
View File
@@ -1,576 +0,0 @@
from SimPEG import *
import BaseGrav as GRAV
import re
class GravityIntegral(Problem.BaseProblem):
# surveyPair = Survey.LinearSurvey
forwardOnly = False #: Determine if the forward matrix is stored (defaut:yes)
actInd = None #: Active cell indices provided
rtype = 'z'
def __init__(self, mesh, mapping=None, **kwargs):
Problem.BaseProblem.__init__(self, mesh, mapping=mapping, **kwargs)
def fwr_op(self):
# Add forward function
# kappa = self.curModel.kappa TODO
rho = self.mapping*self.curModel
if self.forwardOnly:
if getattr(self, 'actInd', None) is not None:
if self.actInd.dtype=='bool':
inds = np.asarray([inds for inds, elem in enumerate(self.actInd, 1) if elem], dtype = int) - 1
else:
inds = self.actInd
else:
inds = np.asarray(range(self.mesh.nC))
nC = len(inds)
# Create active cell projector
P = sp.csr_matrix(
(np.ones(nC), (inds, range(nC))),
shape=(self.mesh.nC, nC)
)
# Create vectors of nodal location (lower and upper corners for each cell)
xn = self.mesh.vectorNx
yn = self.mesh.vectorNy
zn = self.mesh.vectorNz
yn2, xn2, zn2 = np.meshgrid(yn[1:], xn[1:], zn[1:])
yn1, xn1, zn1 = np.meshgrid(yn[0:-1], xn[0:-1], zn[0:-1])
Yn = P.T*np.c_[mkvc(yn1), mkvc(yn2)]
Xn = P.T*np.c_[mkvc(xn1), mkvc(xn2)]
Zn = P.T*np.c_[mkvc(zn1), mkvc(zn2)]
rxLoc = self.survey.srcField.rxList[0].locs
ndata = rxLoc.shape[0]
# Pre-allocate space and create magnetization matrix if required
# Pre-allocate space
if self.rtype == 'z':
fwr_d = np.zeros(self.survey.nRx)
elif self.rtype == 'xyz':
fwr_d = np.zeros(3*self.survey.nRx)
else:
print """Flag must be either 'z' | 'xyz', please revised"""
return
# Add counter to dsiplay progress. Good for large problems
count = -1;
for ii in range(ndata):
tx, ty, tz = get_T_mat(Xn, Yn, Zn, rxLoc[ii, :])
if self.rtype =='z':
fwr_d[ii] =tz.dot(rho)
elif self.rtype =='xyz':
fwr_d[ii] = tx.dot(rho)
fwr_d[ii+ndata] = ty.dot(rho)
fwr_d[ii+2*ndata] = tz.dot(rho)
# Display progress
count = progress(ii,count,ndata)
print "Done 100% ...forward operator completed!!\n"
return fwr_d
else:
return self.G.dot(rho)
def fields(self, m):
self.curModel = m
fields = self.fwr_op()
return fields
# return self.G.dot(self.mapping*(m))
def Jvec(self, m, v, f=None):
dmudm = self.mapping.deriv(m)
return self.G.dot(dmudm*v)
def Jtvec(self, m, v, f=None):
dmudm = self.mapping.deriv(m)
return dmudm.T * (self.G.T.dot(v))
@property
def G(self):
if not self.ispaired:
raise Exception('Need to pair!')
if getattr(self, '_G', None) is None:
self._G = self.Intrgl_Fwr_Op( 'z' )
return self._G
def Intrgl_Fwr_Op(self, flag):
"""
Gravity forward operator in integral form
flag = 'z' | 'xyz'
Return
_G = Linear forward modeling operation
Created on March, 15th 2016
@author: dominiquef
"""
# Find non-zero cells
# inds = np.nonzero(actv)[0]
if getattr(self, 'actInd', None) is not None:
if self.actInd.dtype=='bool':
inds = np.asarray([inds for inds, elem in enumerate(self.actInd, 1) if elem], dtype = int) - 1
else:
inds = self.actInd
else:
inds = np.asarray(range(self.mesh.nC))
nC = len(inds)
# Create active cell projector
P = sp.csr_matrix(
(np.ones(nC), (inds, range(nC))),
shape=(self.mesh.nC, nC)
)
# Create vectors of nodal location (lower and upper corners for each cell)
xn = self.mesh.vectorNx
yn = self.mesh.vectorNy
zn = self.mesh.vectorNz
yn2, xn2, zn2 = np.meshgrid(yn[1:], xn[1:], zn[1:])
yn1, xn1, zn1 = np.meshgrid(yn[0:-1], xn[0:-1], zn[0:-1])
Yn = P.T*np.c_[mkvc(yn1), mkvc(yn2)]
Xn = P.T*np.c_[mkvc(xn1), mkvc(xn2)]
Zn = P.T*np.c_[mkvc(zn1), mkvc(zn2)]
rxLoc = self.survey.srcField.rxList[0].locs
ndata = rxLoc.shape[0]
# Pre-allocate space and create magnetization matrix if required
# Pre-allocate space
if flag == 'z':
G = np.zeros((ndata, nC))
elif flag == 'xyz':
G = np.zeros((int(3*ndata), nC))
else:
print """Flag must be either 'z' | 'xyz', please revised"""
return
# Loop through all observations and create forward operator (ndata-by-nC)
print "Begin calculation of forward operator: " + flag
# Add counter to dsiplay progress. Good for large problems
count = -1;
for ii in range(ndata):
if flag=='z':
tt = get_T_mat(Xn, Yn, Zn, rxLoc[ii, :])
G[ii, :] = tt
elif flag == 'xyz':
print "Sorry 3-component not implemented yet"
# Display progress
count = progress(ii, count, ndata)
print "Done 100% ...forward operator completed!!\n"
return G
def get_T_mat(Xn, Yn, Zn, rxLoc):
"""
Load in the active nodes of a tensor mesh and computes the gravity tensor
for a given observation location rxLoc[obsx, obsy, obsz]
INPUT:
Xn, Yn, Zn: Node location matrix for the lower and upper most corners of
all cells in the mesh shape[nC,2]
M
OUTPUT:
Tx = [Txx Txy Txz]
Ty = [Tyx Tyy Tyz]
Tz = [Tzx Tzy Tzz]
where each elements have dimension 1-by-nC.
Only the upper half 5 elements have to be computed since symetric.
Currently done as for-loops but will eventually be changed to vector
indexing, once the topography has been figured out.
"""
NewtG=6.6738e-3
eps = 1e-10 # add a small value to the locations to avoid /0
nC = Xn.shape[0]
# Pre-allocate space for 1D array
tx = np.zeros((1,nC))
ty = np.zeros((1,nC))
tz = np.zeros((1,nC))
dz = rxLoc[2] - Zn + eps
dy = Yn - rxLoc[1] + eps
dx = Xn - rxLoc[0] + eps
# Compute contribution from each corners
for aa in range(2):
for bb in range(2):
for cc in range(2):
r = (
dx[:, aa] ** 2 +
dy[:, bb] ** 2 +
dz[:, cc] ** 2
) ** (0.50)
tx = tx - NewtG * (-1) ** aa * (-1) ** bb * (-1) ** cc * (
dy[:, bb] * np.log(dz[:, cc] + r) +
dz[:, cc] * np.log(dy[:, bb] + r) -
dx[:, aa] * np.arctan(dy[:, bb] * dz[:, cc] / (dx[:, aa] * r)))
ty = ty - NewtG * (-1) ** aa * (-1) ** bb * (-1) ** cc * (
dx[:, aa] * np.log(dz[:, cc] + r) +
dz[:, cc] * np.log(dx[:, aa] + r) -
dy[:, bb] * np.arctan(dx[:, aa] * dz[:, cc] / (dy[:, bb] * r)))
tz = tz - NewtG * (-1) ** aa * (-1) ** bb * (-1) ** cc * (
dx[:, aa] * np.log(dy[:, bb] + r) +
dy[:, bb] * np.log(dx[:, aa] + r) -
dz[:, cc] * np.arctan(dx[:, aa] * dy[:, bb] / (dz[:, cc] * r)))
return tx,ty,tz
def progress(iter, prog, final):
"""
progress(iter,prog,final)
Function measuring the progress of a process and print to screen the %.
Useful to estimate the remaining runtime of a large problem.
Created on Dec, 20th 2015
@author: dominiquef
"""
arg = np.floor(float(iter)/float(final)*10.)
if arg > prog:
strg = "Done " + str(arg*10) + " %"
print strg
prog = arg
return prog
def writeUBCobs(filename, survey, d):
"""
writeUBCobs(filename,survey,d)
Function writing an observation file in UBC-GRAV3D format.
INPUT
filename : Name of out file including directory
survey
flag : dobs | dpred
OUTPUT
Obsfile
"""
rxLoc = survey.srcField.rxList[0].locs
wd = survey.std
data = np.c_[rxLoc, d, wd]
with file(filename, 'w') as fid:
fid.write('%i\n' % len(d))
np.savetxt(fid, data, fmt='%e', delimiter=' ', newline='\n')
print "Observation file saved to: " + filename
def getActiveTopo(mesh, topo, flag):
"""
getActiveTopo(mesh,topo)
Function creates an active cell model from topography
INPUT
mesh : Mesh in SimPEG format
topo : Scatter points defining topography [x,y,z]
OUTPUT
actv : Active cell model
"""
import scipy.interpolate as interpolation
if flag == 'N':
Zn = np.zeros((mesh.nNx, mesh.nNy))
# wght = np.zeros((mesh.nNx,mesh.nNy))
cx = mesh.vectorNx
cy = mesh.vectorNy
F = interpolation.NearestNDInterpolator(topo[:, 0:2], topo[:, 2])
[Y, X] = np.meshgrid(cy, cx)
Zn = F(X, Y)
actv = np.zeros((mesh.nCx, mesh.nCy, mesh.nCz))
if flag == 'N':
Nz = mesh.vectorNz[1:]
for jj in range(mesh.nCy):
for ii in range(mesh.nCx):
temp = [kk for kk in range(len(Nz)) if np.all(Zn[ii:(ii+2), jj:(jj+2)] > Nz[kk]) ]
actv[ii, jj, temp] = 1
actv = mkvc(actv == 1)
inds = np.asarray([inds for inds, elem in enumerate(actv, 1) if elem], dtype = int) - 1
return inds
def plot_obs_2D(survey,varstr, fig = None):
""" Function plot_obs(rxLoc,d,wd)
Generate a 2d interpolated plot from scatter points of data
INPUT
rxLoc : Observation locations [x,y,z]
d : Data vector
wd : Uncertainty vector
OUTPUT
figure()
Created on Dec, 27th 2015
@author: dominiquef
"""
from scipy.interpolate import griddata
import pylab as plt
rxLoc = survey.srcField.rxList[0].locs
d = survey.dobs
wd = survey.std
# Create grid of points
x = np.linspace(rxLoc[:,0].min(), rxLoc[:,0].max(), 100)
y = np.linspace(rxLoc[:,1].min(), rxLoc[:,1].max(), 100)
X, Y = np.meshgrid(x,y)
# Interpolate
d_grid = griddata(rxLoc[:,0:2],d,(X,Y), method ='linear')
# Plot result
if fig is None:
fig = plt.figure()
ax = plt.subplot()
plt.imshow(d_grid, extent=[x.min(), x.max(), y.min(), y.max()],origin = 'lower', cmap='plasma')
plt.colorbar(fraction=0.02)
plt.contour(X,Y, d_grid,10)
plt.scatter(rxLoc[:,0],rxLoc[:,1], c=d, s=20)
plt.title(varstr)
plt.gca().set_aspect('equal', adjustable='box')
def readUBCgravObs(obs_file):
"""
Read UBC grav file format
INPUT:
:param fileName, path to the UBC obs grav file
OUTPUT:
:param survey
"""
fid = open(obs_file,'r')
# First line has the number of rows
line = fid.readline()
ndat = np.array(line.split(),dtype=int)
# Pre-allocate space for obsx, obsy, obsz, data, uncert
line = fid.readline()
temp = np.array(line.split(),dtype=float)
d = np.zeros(ndat, dtype=float)
wd = np.zeros(ndat, dtype=float)
locXYZ = np.zeros( (ndat,3), dtype=float)
for ii in range(ndat):
temp = np.array(line.split(),dtype=float)
locXYZ[ii,:] = temp[:3]
d[ii] = temp[3]
wd[ii] = temp[4]
line = fid.readline()
rxLoc = GRAV.RxObs(locXYZ)
srcField = GRAV.SrcField([rxLoc])
survey = GRAV.LinearSurvey(srcField)
survey.dobs = d
survey.std = wd
return survey
def read_GRAVinv_inp(input_file):
"""Read input files for forward modeling MAG data with integral form
INPUT:
input_file: File name containing the forward parameter
OUTPUT:
mshfile
obsfile
topofile
start model
ref model
weightfile
chi_target
as, ax ,ay, az
upper, lower bounds
lp, lqx, lqy, lqz
# All files should be in the working directory, otherwise the path must
# be specified.
Created on Dec 21th, 2015
@author: dominiquef
"""
fid = open(input_file,'r')
# Line 1
line = fid.readline()
l_input = line.split('!')
mshfile = l_input[0].rstrip()
# Line 2
line = fid.readline()
l_input = line.split('!')
obsfile = l_input[0].rstrip()
# Line 3
line = fid.readline()
l_input = re.split('[!\s]',line)
if l_input=='null':
topofile = []
else:
topofile = l_input[0].rstrip()
# Line 4
line = fid.readline()
l_input = re.split('[!\s]',line)
if l_input[0]=='VALUE':
mstart = float(l_input[1])
else:
mstart = l_input[0].rstrip()
# Line 5
line = fid.readline()
l_input = re.split('[!\s]',line)
if l_input[0]=='VALUE':
mref = float(l_input[1])
else:
mref = l_input[0].rstrip()
# Line 7
line = fid.readline()
l_input = re.split('[!\s]',line)
if l_input[0]=='DEFAULT':
wgtfile = None
else:
wgtfile = l_input[0].rstrip()
# Line 8
line = fid.readline()
l_input = re.split('[!\s]',line)
chi = float(l_input[0])
# Line 9
line = fid.readline()
l_input = re.split('[!\s]',line)
val = np.array(l_input[0:4])
alphas = val.astype(np.float)
# Line 10
line = fid.readline()
l_input = re.split('[!\s]',line)
if l_input[0]=='VALUE':
val = np.array(l_input[1:3])
bounds = val.astype(np.float)
else:
bounds = l_input[0].rstrip()
# Line 11
line = fid.readline()
l_input = re.split('[!\s]',line)
if l_input[0]=='VALUE':
val = np.array(l_input[1:6])
lpnorms = val.astype(np.float)
else:
lpnorms = l_input[0].rstrip()
return mshfile, obsfile, topofile, mstart, mref, wgtfile, chi, alphas, bounds, lpnorms
-295
View File
@@ -1,295 +0,0 @@
import re, os
from SimPEG import Mesh, np, Utils
import BaseGrav, Gravity
class GravityDriver_Inv(object):
"""docstring for GravityDriver_Inv"""
def __init__(self, input_file=None):
if input_file is not None:
self.basePath = os.path.sep.join(input_file.split(os.path.sep)[:-1])
if len(self.basePath) > 0:
self.basePath += os.path.sep
self.readDriverFile(input_file.split(os.path.sep)[-1])
def readDriverFile(self, input_file):
"""
Read input files for forward modeling GRAV data with integral form
INPUT:
input_file: File name containing the forward parameter
OUTPUT:
mshfile
obsfile
topofile
start model
ref model
active cells model
weightfile
chi_target
as, ax ,ay, az
upper, lower bounds
lp, lqx, lqy, lqz
eps_p, eps_q
# All files should be in the working directory, otherwise the path must
# be specified.
"""
fid = open(self.basePath + input_file, 'r')
# Line 1
line = fid.readline()
l_input = line.split('!')
mshfile = l_input[0].rstrip()
# Line 2
line = fid.readline()
l_input = line.split('!')
obsfile = l_input[0].rstrip()
# Line 3
line = fid.readline()
l_input = re.split('[!\s]', line)
if l_input=='null':
topofile = []
else:
topofile = l_input[0].rstrip()
# Line 4
line = fid.readline()
l_input = re.split('[!\s]', line)
if l_input[0]=='VALUE':
mstart = float(l_input[1])
else:
mstart = l_input[0].rstrip()
# Line 5
line = fid.readline()
l_input = re.split('[!\s]',line)
if l_input[0]=='VALUE':
mref = float(l_input[1])
else:
mref = l_input[0].rstrip()
# Line 6
line = fid.readline()
l_input = re.split('[!\s]', line)
if l_input[0]=='VALUE':
staticInput = float(l_input[1])
elif l_input[0]=='DEFAULT':
staticInput = None
else:
staticInput = l_input[0].rstrip()
# Line 7
line = fid.readline()
l_input = re.split('[!\s]', line)
if l_input=='DEFAULT':
wgtfile = []
else:
wgtfile = l_input[0].rstrip()
# Line 8
line = fid.readline()
l_input = re.split('[!\s]', line)
chi = float(l_input[0])
# Line 9
line = fid.readline()
l_input = re.split('[!\s]', line)
val = np.array(l_input[0:4])
alphas = val.astype(np.float)
# Line 10
line = fid.readline()
l_input = re.split('[!\s]', line)
if l_input[0]=='VALUE':
val = np.array(l_input[1:3])
bounds = val.astype(np.float)
else:
bounds = l_input[0].rstrip()
# Line 11
line = fid.readline()
l_input = re.split('[!\s]', line)
if l_input[0]=='VALUE':
val = np.array(l_input[1:6])
lpnorms = val.astype(np.float)
else:
lpnorms = l_input[0].rstrip()
# Line 12
line = fid.readline()
l_input = re.split('[!\s]', line)
if l_input[0]=='VALUE':
val = np.array(l_input[1:3])
eps = val.astype(np.float)
else:
eps = [None, None]
self.mshfile = mshfile
self.obsfile = obsfile
self.topofile = topofile
self.mstart = mstart
self._mrefInput = mref
self._staticInput = staticInput
self.wgtfile = wgtfile
self.chi = chi
self.alphas = alphas
self.bounds = bounds
self.lpnorms = lpnorms
self.eps = eps
@property
def mesh(self):
if getattr(self, '_mesh', None) is None:
self._mesh = Mesh.TensorMesh.readUBC(self.basePath + self.mshfile)
return self._mesh
@property
def survey(self):
if getattr(self, '_survey', None) is None:
self._survey = self.readGravityObservations(self.basePath + self.obsfile)
return self._survey
@property
def activeCells(self):
if getattr(self, '_activeCells', None) is None:
if self.topofile == 'null':
self._activeCells = np.arange(mesh.nC)
else:
topo = np.genfromtxt(self.basePath + self.topofile, skip_header=1)
# Find the active cells
active = Utils.surface2ind_topo(self.mesh,topo,'N')
inds = np.asarray([inds for inds, elem in enumerate(active, 1) if elem], dtype = int) - 1
self._activeCells = inds
return self._activeCells
@property
def staticCells(self):
if getattr(self, '_staticCells', None) is None:
if getattr(self, '_staticInput', None) is None:
# All cells are dynamic: 1's
self._dynamicCells = np.arange(len(self.m0))
self._staticCells = []
# Cells with specific value are static: 0's
else:
if isinstance(self._staticInput, float):
staticCells = self.m0 == self._staticInput
else:
# Read from file active cells with 0:air, 1:dynamic, -1 static
staticCells = Mesh.TensorMesh.readModelUBC(self.mesh, self.basePath + self._staticInput)
staticCells = staticCells[self.activeCells] == -1
inds = np.asarray([inds for inds, elem in enumerate(staticCells, 1) if elem], dtype = int) - 1
self._staticCells = inds
return self._staticCells
@property
def dynamicCells(self):
if getattr(self, '_dynamicCells', None) is None:
if getattr(self, '_staticInput', None) is None:
# All cells are dynamic: 1's
self._dynamicCells = np.arange(len(self.m0))
# Cells with specific value are static: 0's
else:
if isinstance(self._staticInput, float):
dynamicCells = self.m0 != self._staticInput
else:
# Read from file active cells with 0:air, 1:dynamic, -1 static
dynamicCells = Mesh.TensorMesh.readModelUBC(self.mesh, self.basePath + self._staticInput)
dynamicCells = dynamicCells[self.activeCells] == 1
inds = np.asarray([inds for inds, elem in enumerate(dynamicCells, 1) if elem], dtype = int) - 1
self._dynamicCells = inds
return self._dynamicCells
@property
def nC(self):
if getattr(self, '_nC', None) is None:
self._nC = len(self.activeCells)
return self._nC
@property
def m0(self):
if getattr(self, '_m0', None) is None:
if isinstance(self.mstart, float):
self._m0 = np.ones(self.nC) * self.mstart
else:
self._m0 = Mesh.TensorMesh.readModelUBC(self.mesh, self.basePath + self.mstart)
self._m0 = self._m0[self.activeCells]
return self._m0
@property
def mref(self):
if getattr(self, '_mref', None) is None:
if isinstance(self._mrefInput, float):
self._mref = np.ones(self.nC) * self._mrefInput
else:
self._mref = Mesh.TensorMesh.readModelUBC(self.mesh, self.basePath + self._mrefInput)
self._mref = self._mref[self.activeCells]
return self._mref
def readGravityObservations(self, obs_file):
"""
Read UBC grav file format
INPUT:
:param fileName, path to the UBC obs grav file
OUTPUT:
:param survey
"""
fid = open(obs_file,'r')
# First line has the number of rows
line = fid.readline()
ndat = np.array(line.split(),dtype=int)
# Pre-allocate space for obsx, obsy, obsz, data, uncert
line = fid.readline()
temp = np.array(line.split(),dtype=float)
d = np.zeros(ndat, dtype=float)
wd = np.zeros(ndat, dtype=float)
locXYZ = np.zeros( (ndat,3), dtype=float)
for ii in range(ndat):
temp = np.array(line.split(),dtype=float)
locXYZ[ii,:] = temp[:3]
d[ii] = temp[3]
wd[ii] = temp[4]
line = fid.readline()
rxLoc = BaseGrav.RxObs(locXYZ)
srcField = BaseGrav.SrcField([rxLoc])
survey = BaseGrav.LinearSurvey(srcField)
survey.dobs = d
survey.std = wd
return survey
-278
View File
@@ -1,278 +0,0 @@
from scipy.constants import mu_0
from SimPEG import *
from SimPEG.Utils import kron3, speye, sdiag
import matplotlib.pyplot as plt
def spheremodel(mesh, x0, y0, z0, r):
"""
Generate model indicies for sphere
- (x0, y0, z0 ): is the center location of sphere
- r: is the radius of the sphere
- it returns logical indicies of cell-center model
"""
ind = np.sqrt( (mesh.gridCC[:,0]-x0)**2+(mesh.gridCC[:,1]-y0)**2+(mesh.gridCC[:,2]-z0)**2 ) < r
return ind
def MagSphereAnaFun(x, y, z, R, x0, y0, z0, mu1, mu2, H0, flag='total'):
"""
test
Analytic function for Magnetics problem. The set up here is
magnetic sphere in whole-space assuming that the inducing field is oriented in the x-direction.
* (x0,y0,z0)
* (x0, y0, z0 ): is the center location of sphere
* r: is the radius of the sphere
.. math::
\mathbf{H}_0 = H_0\hat{x}
"""
if (~np.size(x)==np.size(y)==np.size(z)):
print "Specify same size of x, y, z"
return
dim = x.shape
x = Utils.mkvc(x)
y = Utils.mkvc(y)
z = Utils.mkvc(z)
ind = np.sqrt((x-x0)**2+(y-y0)**2+(z-z0)**2 ) < R
r = Utils.mkvc(np.sqrt((x-x0)**2+(y-y0)**2+(z-z0)**2 ))
Bx = np.zeros(x.size)
By = np.zeros(x.size)
Bz = np.zeros(x.size)
# Inside of the sphere
rf2 = 3*mu1/(mu2+2*mu1)
if flag is 'total' and any(ind):
Bx[ind] = mu2*H0*(rf2)
elif (flag == 'secondary'):
Bx[ind] = mu2*H0*(rf2)-mu1*H0
By[ind] = 0.
Bz[ind] = 0.
# Outside of the sphere
rf1 = (mu2-mu1)/(mu2+2*mu1)
if (flag == 'total'):
Bx[~ind] = mu1*(H0+H0/r[~ind]**5*(R**3)*rf1*(2*(x[~ind]-x0)**2-(y[~ind]-y0)**2-(z[~ind]-z0)**2))
elif (flag == 'secondary'):
Bx[~ind] = mu1*(H0/r[~ind]**5*(R**3)*rf1*(2*(x[~ind]-x0)**2-(y[~ind]-y0)**2-(z[~ind]-z0)**2))
By[~ind] = mu1*(H0/r[~ind]**5*(R**3)*rf1*(3*(x[~ind]-x0)*(y[~ind]-y0)))
Bz[~ind] = mu1*(H0/r[~ind]**5*(R**3)*rf1*(3*(x[~ind]-x0)*(z[~ind]-z0)))
return np.reshape(Bx, x.shape, order='F'), np.reshape(By, x.shape, order='F'), np.reshape(Bz, x.shape, order='F')
def CongruousMagBC(mesh, Bo, chi):
"""
Computing boundary condition using Congrous sphere method.
This is designed for secondary field formulation.
>> Input
* mesh: Mesh class
* Bo: np.array([Box, Boy, Boz]): Primary magnetic flux
* chi: susceptibility at cell volume
.. math::
\\vec{B}(r) = \\frac{\mu_0}{4\pi} \\frac{m}{ \| \\vec{r} - \\vec{r}_0\|^3}[3\hat{m}\cdot\hat{r}-\hat{m}]
"""
ind = chi > 0.
V = mesh.vol[ind].sum()
gamma = 1/V*(chi*mesh.vol).sum() # like a mass!
Bot = np.sqrt(sum(Bo**2))
mx = Bo[0]/Bot
my = Bo[1]/Bot
mz = Bo[2]/Bot
mom = 1/mu_0*Bot*gamma*V/(1+gamma/3)
xc = sum(chi[ind]*mesh.gridCC[:,0][ind])/sum(chi[ind])
yc = sum(chi[ind]*mesh.gridCC[:,1][ind])/sum(chi[ind])
zc = sum(chi[ind]*mesh.gridCC[:,2][ind])/sum(chi[ind])
indxd, indxu, indyd, indyu, indzd, indzu = mesh.faceBoundaryInd
const = mu_0/(4*np.pi)*mom
rfun = lambda x: np.sqrt((x[:,0]-xc)**2 + (x[:,1]-yc)**2 + (x[:,2]-zc)**2)
mdotrx = (mx*(mesh.gridFx[(indxd|indxu),0]-xc)/rfun(mesh.gridFx[(indxd|indxu),:]) +
my*(mesh.gridFx[(indxd|indxu),1]-yc)/rfun(mesh.gridFx[(indxd|indxu),:]) +
mz*(mesh.gridFx[(indxd|indxu),2]-zc)/rfun(mesh.gridFx[(indxd|indxu),:]))
Bbcx = const/(rfun(mesh.gridFx[(indxd|indxu),:])**3)*(3*mdotrx*(mesh.gridFx[(indxd|indxu),0]-xc)/rfun(mesh.gridFx[(indxd|indxu),:])-mx)
mdotry = (mx*(mesh.gridFy[(indyd|indyu),0]-xc)/rfun(mesh.gridFy[(indyd|indyu),:]) +
my*(mesh.gridFy[(indyd|indyu),1]-yc)/rfun(mesh.gridFy[(indyd|indyu),:]) +
mz*(mesh.gridFy[(indyd|indyu),2]-zc)/rfun(mesh.gridFy[(indyd|indyu),:]))
Bbcy = const/(rfun(mesh.gridFy[(indyd|indyu),:])**3)*(3*mdotry*(mesh.gridFy[(indyd|indyu),1]-yc)/rfun(mesh.gridFy[(indyd|indyu),:])-my)
mdotrz = (mx*(mesh.gridFz[(indzd|indzu),0]-xc)/rfun(mesh.gridFz[(indzd|indzu),:]) +
my*(mesh.gridFz[(indzd|indzu),1]-yc)/rfun(mesh.gridFz[(indzd|indzu),:]) +
mz*(mesh.gridFz[(indzd|indzu),2]-zc)/rfun(mesh.gridFz[(indzd|indzu),:]))
Bbcz = const/(rfun(mesh.gridFz[(indzd|indzu),:])**3)*(3*mdotrz*(mesh.gridFz[(indzd|indzu),2]-zc)/rfun(mesh.gridFz[(indzd|indzu),:])-mz)
return np.r_[Bbcx, Bbcy, Bbcz], (1/gamma-1/(3+gamma))*1/V
def MagSphereAnaFunA(x, y, z, R, xc, yc, zc, chi, Bo, flag):
"""
Computing boundary condition using Congrous sphere method.
This is designed for secondary field formulation.
>> Input
mesh: Mesh class
Bo: np.array([Box, Boy, Boz]): Primary magnetic flux
Chi: susceptibility at cell volume
.. math::
\\vec{B}(r) = \\frac{\mu_0}{4\pi}\\frac{m}{\| \\vec{r}-\\vec{r}_0\|^3}[3\hat{m}\cdot\hat{r}-\hat{m}]
"""
if (~np.size(x)==np.size(y)==np.size(z)):
print "Specify same size of x, y, z"
return
dim = x.shape
x = Utils.mkvc(x)
y = Utils.mkvc(y)
z = Utils.mkvc(z)
Bot = np.sqrt(sum(Bo**2))
mx = Bo[0]/Bot
my = Bo[1]/Bot
mz = Bo[2]/Bot
ind = np.sqrt((x-xc)**2+(y-yc)**2+(z-zc)**2 ) < R
Bx = np.zeros(x.size)
By = np.zeros(x.size)
Bz = np.zeros(x.size)
# Inside of the sphere
rf2 = 3/(chi+3)*(1+chi)
if (flag == 'total'):
Bx[ind] = Bo[0]*(rf2)
By[ind] = Bo[1]*(rf2)
Bz[ind] = Bo[2]*(rf2)
elif (flag == 'secondary'):
Bx[ind] = Bo[0]*(rf2)-Bo[0]
By[ind] = Bo[1]*(rf2)-Bo[1]
Bz[ind] = Bo[2]*(rf2)-Bo[2]
r = Utils.mkvc(np.sqrt((x-xc)**2+(y-yc)**2+(z-zc)**2 ))
V = 4*np.pi*R**3/3
mom = Bot/mu_0*chi/(1+chi/3)*V
const = mu_0/(4*np.pi)*mom
mdotr = (mx*(x[~ind]-xc)/r[~ind] + my*(y[~ind]-yc)/r[~ind] + mz*(z[~ind]-zc)/r[~ind])
Bx[~ind] = const/(r[~ind]**3)*(3*mdotr*(x[~ind]-xc)/r[~ind]-mx)
By[~ind] = const/(r[~ind]**3)*(3*mdotr*(y[~ind]-yc)/r[~ind]-my)
Bz[~ind] = const/(r[~ind]**3)*(3*mdotr*(z[~ind]-zc)/r[~ind]-mz)
return Bx, By, Bz
def IDTtoxyz(Inc, Dec, Btot):
"""
Convert from Inclination, Declination, Total intensity of earth field to x, y, z
"""
Bx = Btot*np.cos(Inc/180.*np.pi)*np.sin(Dec/180.*np.pi)
By = Btot*np.cos(Inc/180.*np.pi)*np.cos(Dec/180.*np.pi)
Bz = -Btot*np.sin(Inc/180.*np.pi)
return np.r_[Bx, By, Bz]
def MagSphereFreeSpace(x, y, z, R, xc, yc, zc, chi, Bo):
"""
Computing boundary condition using Congrous sphere method.
This is designed for secondary field formulation.
>> Input
mesh: Mesh class
Bo: np.array([Box, Boy, Boz]): Primary magnetic flux
Chi: susceptibility at cell volume
.. math::
\\vec{B}(r) = \\frac{\mu_0}{4\pi}\\frac{m}{\| \\vec{r}-\\vec{r}_0\|^3}[3\hat{m}\cdot\hat{r}-\hat{m}]
"""
if (~np.size(x)==np.size(y)==np.size(z)):
print "Specify same size of x, y, z"
return
x = Utils.mkvc(x)
y = Utils.mkvc(y)
z = Utils.mkvc(z)
nobs = len(x)
Bot = np.sqrt(sum(Bo**2))
mx = np.ones([nobs]) * Bo[0,0] * R**3 / 3. * chi
my = np.ones([nobs]) * Bo[0,1] * R**3 / 3. * chi
mz = np.ones([nobs]) * Bo[0,2] * R**3 / 3. * chi
M = np.c_[mx, my, mz]
rx = (x - xc)
ry = (y - yc)
rz = (zc - z)
rvec = np.c_[rx, ry, rz]
r = np.sqrt((rx)**2+(ry)**2+(rz)**2 )
B = -Utils.sdiag(1./r**3)*M + Utils.sdiag((3 * np.sum(M*rvec,axis=1))/r**5)*rvec
Bx = B[:,0]
By = B[:,1]
Bz = B[:,2]
return Bx, By, Bz
if __name__ == '__main__':
hxind = [(0,25,1.3),(21, 12.5),(0,25,1.3)]
hyind = [(0,25,1.3),(21, 12.5),(0,25,1.3)]
hzind = [(0,25,1.3),(20, 12.5),(0,25,1.3)]
# hx, hy, hz = Utils.meshTensors(hxind, hyind, hzind)
M3 = Mesh.TensorMesh([hxind, hyind, hzind], "CCC")
indxd, indxu, indyd, indyu, indzd, indzu = M3.faceBoundaryInd
mu0 = 4*np.pi*1e-7
chibkg = 0.
chiblk = 0.01
chi = np.ones(M3.nC)*chibkg
sph_ind = spheremodel(M3, 0, 0, 0, 100)
chi[sph_ind] = chiblk
mu = (1.+chi)*mu0
Bbc, const = CongruousMagBC(M3, np.array([1., 0., 0.]), chi)
flag = 'secondary'
Box = 1.
H0 = Box/mu_0
Bbcxx, Bbcxy, Bbcxz = MagSphereAnaFun(M3.gridFx[(indxd|indxu),0], M3.gridFx[(indxd|indxu),1], M3.gridFx[(indxd|indxu),2], 100, 0., 0., 0., mu_0, mu_0*(1+chiblk), H0, flag)
Bbcyx, Bbcyy, Bbcyz = MagSphereAnaFun(M3.gridFy[(indyd|indyu),0], M3.gridFy[(indyd|indyu),1], M3.gridFy[(indyd|indyu),2], 100, 0., 0., 0., mu_0, mu_0*(1+chiblk), H0, flag)
Bbczx, Bbczy, Bbczz = MagSphereAnaFun(M3.gridFz[(indzd|indzu),0], M3.gridFz[(indzd|indzu),1], M3.gridFz[(indzd|indzu),2], 100, 0., 0., 0., mu_0, mu_0*(1+chiblk), H0, flag)
Bbc_ana = np.r_[Bbcxx, Bbcyy, Bbczz]
# fig, ax = plt.subplots(1,1, figsize = (10, 10))
# ax.plot(Bbc_ana)
# ax.plot(Bbc)
# plt.show()
err = np.linalg.norm(Bbc-Bbc_ana)/np.linalg.norm(Bbc_ana)
if err < 0.1:
print 'Mag Boundary computation is valid, err = ', err
else:
print 'Mag Boundary computation is wrong!!, err = ', err
pass
File diff suppressed because it is too large Load Diff
-334
View File
@@ -1,334 +0,0 @@
import re, os
from SimPEG import Mesh, np, Utils
import BaseMag, Magnetics
class MagneticsDriver_Inv(object):
"""docstring for MagneticsDriver_Inv"""
def __init__(self, input_file=None):
if input_file is not None:
self.basePath = os.path.sep.join(input_file.split(os.path.sep)[:-1])
if len(self.basePath) > 0:
self.basePath += os.path.sep
self.readDriverFile(input_file.split(os.path.sep)[-1])
def readDriverFile(self, input_file):
"""
Read input files for forward modeling MAG data with integral form
INPUT:
input_file: File name containing the forward parameter
OUTPUT:
mshfile
obsfile
topofile
start model
ref model
mag model
weightfile
chi_target
as, ax ,ay, az
upper, lower bounds
lp, lqx, lqy, lqz
# All files should be in the working directory, otherwise the path must
# be specified.
"""
fid = open(self.basePath + input_file,'r')
# Line 1
line = fid.readline()
l_input = line.split('!')
mshfile = l_input[0].rstrip()
# Line 2
line = fid.readline()
l_input = line.split('!')
obsfile = l_input[0].rstrip()
# Line 3
line = fid.readline()
l_input = re.split('[!\s]',line)
if l_input=='null':
topofile = []
else:
topofile = l_input[0].rstrip()
# Line 4
line = fid.readline()
l_input = re.split('[!\s]',line)
if l_input[0]=='VALUE':
mstart = float(l_input[1])
else:
mstart = l_input[0].rstrip()
# Line 5
line = fid.readline()
l_input = re.split('[!\s]',line)
if l_input[0]=='VALUE':
mref = float(l_input[1])
else:
mref = l_input[0].rstrip()
# Line 6
line = fid.readline()
l_input = re.split('[!\s]',line)
if l_input[0]=='VALUE':
staticInput = float(l_input[1])
elif l_input[0]=='DEFAULT':
staticInput = None
else:
staticInput = l_input[0].rstrip()
# Line 7
line = fid.readline()
l_input = re.split('[!\s]',line)
if l_input=='DEFAULT':
magfile = []
else:
magfile = l_input[0].rstrip()
# Line 8
line = fid.readline()
l_input = re.split('[!\s]',line)
if l_input=='DEFAULT':
wgtfile = []
else:
wgtfile = l_input[0].rstrip()
# Line 9
line = fid.readline()
l_input = re.split('[!\s]',line)
chi = float(l_input[0])
# Line 10
line = fid.readline()
l_input = re.split('[!\s]',line)
val = np.array(l_input[0:4])
alphas = val.astype(np.float)
# Line 11
line = fid.readline()
l_input = re.split('[!\s]',line)
if l_input[0]=='VALUE':
val = np.array(l_input[1:3])
bounds = val.astype(np.float)
else:
bounds = l_input[0].rstrip()
# Line 12
line = fid.readline()
l_input = re.split('[!\s]',line)
if l_input[0]=='VALUE':
val = np.array(l_input[1:6])
lpnorms = val.astype(np.float)
else:
lpnorms = l_input[0].rstrip()
# Line 13
line = fid.readline()
l_input = re.split('[!\s]',line)
if l_input[0]=='VALUE':
val = np.array(l_input[1:3])
eps = val.astype(np.float)
else:
eps = [None,None]
self.mshfile = mshfile
self.obsfile = obsfile
self.topofile = topofile
self.mstart = mstart
self._mrefInput = mref
self._staticInput = staticInput
self.magfile = magfile
self.wgtfile = wgtfile
self.chi = chi
self.alphas = alphas
self.bounds = bounds
self.lpnorms = lpnorms
self.eps = eps
@property
def mesh(self):
if getattr(self, '_mesh', None) is None:
self._mesh = Mesh.TensorMesh.readUBC(self.basePath + self.mshfile)
return self._mesh
@property
def survey(self):
if getattr(self, '_survey', None) is None:
self._survey = self.readMagneticsObservations(self.obsfile)
return self._survey
@property
def activeCells(self):
if getattr(self, '_activeCells', None) is None:
if self.topofile == 'null':
self._activeCells = np.arange(self.mesh.nC)
else:
topo = np.genfromtxt(self.basePath + self.topofile, skip_header=1)
# Find the active cells
active = Utils.surface2ind_topo(self.mesh,topo,'N')
inds = np.asarray([inds for inds, elem in enumerate(active, 1) if elem], dtype = int) - 1
self._activeCells = inds
return self._activeCells
@property
def staticCells(self):
if getattr(self, '_staticCells', None) is None:
if getattr(self, '_staticInput', None) is None:
# All cells are dynamic: 1's
self._dynamicCells = np.arange(len(self.m0))
self._staticCells = []
# Cells with specific value are static: 0's
else:
if isinstance(self._staticInput, float):
staticCells = self.m0 == self._staticInput
else:
# Read from file active cells with 0:air, 1:dynamic, -1 static
staticCells = Mesh.TensorMesh.readModelUBC(self.mesh, self.basePath + self._staticInput)
staticCells = staticCells[self.activeCells] == -1
inds = np.asarray([inds for inds, elem in enumerate(staticCells, 1) if elem], dtype = int) - 1
self._staticCells = inds
return self._staticCells
@property
def dynamicCells(self):
if getattr(self, '_dynamicCells', None) is None:
if getattr(self, '_staticInput', None) is None:
# All cells are dynamic: 1's
self._dynamicCells = np.arange(len(self.m0))
# Cells with specific value are static: 0's
else:
if isinstance(self._staticInput, float):
dynamicCells = self.m0 != self._staticInput
else:
# Read from file active cells with 0:air, 1:dynamic, -1 static
dynamicCells = Mesh.TensorMesh.readModelUBC(self.mesh, self.basePath + self._staticInput)
dynamicCells = dynamicCells[self.activeCells] == 1
inds = np.asarray([inds for inds, elem in enumerate(dynamicCells, 1) if elem], dtype = int) - 1
self._dynamicCells = inds
return self._dynamicCells
@property
def nC(self):
if getattr(self, '_nC', None) is None:
self._nC = len(self.activeCells)
return self._nC
@property
def m0(self):
if getattr(self, '_m0', None) is None:
if isinstance(self.mstart, float):
self._m0 = np.ones(self.nC) * self.mstart
else:
self._m0 = Mesh.TensorMesh.readModelUBC(self.mesh,self.basePath + self.mstart)
self._m0 = self._m0[self.activeCells]
return self._m0
@property
def mref(self):
if getattr(self, '_mref', None) is None:
if isinstance(self._mrefInput, float):
self._mref = np.ones(self.nC) * self._mrefInput
else:
self._mref = Mesh.TensorMesh.readModelUBC(self.mesh,self.basePath + self._mrefInput)
self._mref = self._mref[self.activeCells]
return self._mref
@property
def magnetizationModel(self):
"""
magnetization vector
"""
if self.magfile == 'DEFAULT':
return Magnetics.dipazm_2_xyz(np.ones(self.nC) * self.survey.srcField.param[1], np.ones(self.nC) * self.survey.srcField.param[2])
else:
raise NotImplementedError("this will require you to read in a three column vector model")
self._mref = Utils.meshutils.readUBCTensorModel(self.basePath + self._mrefInput, self.mesh)
return np.genfromtxt(self.magfile,delimiter=' \n',dtype=np.str,comments='!')
def readMagneticsObservations(self, obs_file):
"""
Read and write UBC mag file format
INPUT:
:param fileName, path to the UBC obs mag file
OUTPUT:
:param survey
:param M, magnetization orentiaton (MI, MD)
"""
fid = open(self.basePath + obs_file,'r')
# First line has the inclination,declination and amplitude of B0
line = fid.readline()
B = np.array(line.split(),dtype=float)
# Second line has the magnetization orientation and a flag
line = fid.readline()
M = np.array(line.split(),dtype=float)
# Third line has the number of rows
line = fid.readline()
ndat = np.array(line.split(),dtype=int)
# Pre-allocate space for obsx, obsy, obsz, data, uncert
line = fid.readline()
temp = np.array(line.split(),dtype=float)
d = np.zeros(ndat, dtype=float)
wd = np.zeros(ndat, dtype=float)
locXYZ = np.zeros( (ndat,3), dtype=float)
for ii in range(ndat):
temp = np.array(line.split(),dtype=float)
locXYZ[ii,:] = temp[:3]
if len(temp) > 3:
d[ii] = temp[3]
if len(temp)==5:
wd[ii] = temp[4]
line = fid.readline()
rxLoc = BaseMag.RxObs(locXYZ)
srcField = BaseMag.SrcField([rxLoc],param=(B[2],B[0],B[1]))
survey = BaseMag.LinearSurvey(srcField)
survey.dobs = d
survey.std = wd
return survey
-7
View File
@@ -1,7 +0,0 @@
import MagAnalytics
import BaseMag
import Magnetics
import BaseGrav
import Gravity
import MagneticsDriver
import GravityDriver
+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

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