Compare commits

...
Author SHA1 Message Date
gitea 714afa860f misc fixes
- 2to3 on setup.py to allow python 3 import
- removed uneeded circular import in Maps.py
- added missing old_div o mathutils
2016-07-21 14:20:37 +08:00
Brendan Smithyman f7a70aa6a7 Possibly deal with a good chunk of the old_divs 2016-07-17 16:02:43 -05:00
Brendan Smithyman 7189ec5b2f Test Python3. 2016-07-16 14:37:45 -05:00
Brendan Smithyman 4c2a61bf54 Fix deps. 2016-07-16 14:35:56 -05:00
Brendan Smithyman ca8d8f8c2d Futurize 1, futurize 2, pasteurize. 2016-07-16 14:17:02 -05:00
Lindsey Heagy 362975d2bd Merge pull request #359 from simpeg/dev
Dev
2016-07-15 15:26:39 -05:00
Lindsey Heagy cf34415ae8 Merge pull request #358 from simpeg/feat/codecov
Feat/codecov
2016-07-14 08:51:57 -06:00
Lindsey Heagy 09ec5621ae remove coveralls 2016-07-13 19:46:27 -07:00
Lindsey Heagy 4e583fc566 fix indentation level 2016-07-13 13:40:38 -07:00
Lindsey Heagy ea62998250 use codecov.io 2016-07-13 12:15:15 -07:00
SEOGI KANG 4796b0f91f Merge pull request #357 from simpeg/analytics
Analytics
2016-06-30 00:25:28 -07:00
seogi_macbook 52b25e2dc5 Merge branch 'dev' of https://github.com/simpeg/simpeg into analytics 2016-06-30 00:23:12 -07:00
seogi_macbook a289b656cd Fixes for kwargs variables in FDEMDipolarfields.py 2016-06-29 13:09:11 -07:00
Lindsey Heagy 334cd8e454 Bump version: 0.1.11 → 0.1.12 2016-06-29 09:49:29 -07:00
Lindsey Heagy ecbdd90f63 Merge pull request #354 from simpeg/dev
Two new examples.
2016-06-29 09:46:33 -07:00
dfournier 394dc9106a Merge pull request #332 from simpeg/ref/regularization
Automate the epsilon picking based on percentile of model values for …
2016-06-29 08:39:29 -07:00
seogi_macbook eda2394411 fix bug for omega. 2016-06-27 13:04:30 -07:00
Rowan Cockett 3deca9ed77 Merge pull request #351 from simpeg/example/mesh2mesh
Mesh2Mesh and Combo Map examples.
2016-06-26 21:22:28 -06:00
Rowan Cockett ba173674ec Mesh2Mesh and Combo Map examples.
Also fixed plotting codes to show the plots by default.
2016-06-26 17:07:07 -06:00
Rowan Cockett 303da372aa Merged branch master into dev 2016-06-26 16:31:24 -06:00
Rowan Cockett 6d6e7fc8bd Merge pull request #350 from simpeg/fix/docs-images
Update index.rst
2016-06-26 16:30:29 -06:00
Rowan Cockett 8ed3ec18fa Update README.rst 2016-06-26 16:29:20 -06:00
Rowan Cockett 3960cfc313 Update index.rst 2016-06-26 16:17:30 -06:00
Rowan Cockett 2eba0b841f Merge pull request #338 from simpeg/dev
Dev
2016-06-26 14:01:03 -06:00
sgkang c79bb998cb Merge pull request #348 from simpeg/analytics
fix minor bugs in analytics (just for binder deploy)
2016-06-23 14:14:29 -07:00
seogi_macbook 0763925743 fix minor bugs in analytics 2016-06-23 14:12:49 -07:00
sgkang 1a0b81a206 Merge pull request #341 from simpeg/analytics
Analytics
2016-06-23 11:33:33 -07:00
seogi_macbook 8b44f8d96b change FDEM_fields.py to FDEMDipolarfield.py 2016-06-23 10:14:13 -07:00
seogi_macbook 2bfd01ed7c Merge branch 'dev' of https://github.com/simpeg/simpeg into analytics 2016-06-23 10:09:56 -07:00
seogi_macbook e1ba80883d Incorporate Lindsey's suggestoins 2016-06-23 09:10:50 -07:00
Rowan Cockett a54713f546 Change to NotImplementedError. 2016-06-22 11:32:54 -06:00
Lindsey Heagy ef382aed85 Merge branch 'master' into dev 2016-06-21 18:43:28 -06:00
seogi_macbook 1b33804e5a Merge branch 'dev' of https://github.com/simpeg/simpeg into analytics 2016-06-21 11:15:23 -07:00
seogi_macbook c161c5eab2 Merge branch 'master' of https://github.com/simpeg/simpeg into analytics 2016-06-21 11:14:56 -07:00
micmitch f0944362c8 Silly mistake... needed zero arrays instead of scalars. 2016-06-14 15:22:36 -07:00
micmitch e815ddaec7 Removed d typos from the end of function names. 2016-06-14 15:15:53 -07:00
Lindsey 64b0b4561f Merge pull request #337 from simpeg/feat/docs-gae-travis-deploy
Feat/docs gae travis deploy
2016-06-13 19:16:22 -06:00
micmitch b9d30af4a8 Changed \sigma to \hat{\sigma} = \sigma + i \omega \epsilon in the E field calculations. 2016-06-13 17:47:46 -07:00
micmitch 093f441331 Added functions to split electric field into "galvanic" and "inductive" portions. 2016-06-13 17:35:52 -07:00
seogi_macbook 5e3c1da8e2 add place holder for galvanic and inductive electric fields... 2016-06-13 16:41:59 -07:00
seogi_macbook a3a5c86008 Fix couple bugs in FDEM analytics 2016-06-11 21:58:03 +02:00
Lindsey Heagy 7a82f57367 - doc the simple regularization, sparse regularization and regularization mesh
- typo fix in MT_3D_Forward example - Solver
- use EM.Static.DC for DC example
2016-06-11 08:16:22 -07:00
Lindsey Heagy 845c3c10ab version on travis commit 2016-06-11 07:48:24 -07:00
Lindsey Heagy c4d86b4a29 Merge branch 'dev' into feat/docs-gae-travis-deploy 2016-06-11 07:14:00 -07:00
Lindsey Heagy 3fc855f3c9 only push on master 2016-06-10 20:50:41 -07:00
Lindsey Heagy fea508a507 fix path to app.yaml 2016-06-10 19:39:02 -07:00
Lindsey Heagy 5ca52cf49f linting travis deploy of docs 2016-06-10 19:14:24 -07:00
Lindsey Heagy 7a06453c42 editing travis deploy 2016-06-10 18:54:31 -07:00
Lindsey Heagy 291da78b97 try deploying only from docs branch 2016-06-10 15:16:56 -07:00
Lindsey Heagy 504592c8de travis typo fixes 2016-06-01 08:12:21 -07:00
Lindsey Heagy 75647f8fc3 first pass at gae deploy 2016-06-01 00:00:27 -07:00
Lindsey Heagy e3462666bd working on travis decrypt 2016-05-31 23:27:32 -07:00
Lindsey Heagy fb37bf0fe2 move unpacking of credentials to after success, point credentials to docs folder 2016-05-31 23:15:50 -07:00
Lindsey Heagy b023adbb33 add encrypted credentials for gae site and decrypt on travis 2016-05-31 23:03:55 -07:00
Lindsey Heagy 279bd49b4c Bump version: 0.1.10 → 0.1.11 2016-05-31 15:38:05 -07:00
Lindsey Heagy 3398d5ab4d use templates in conf.py, add google analytics to docs 2016-05-31 15:33:19 -07:00
Lindsey 3e26bb7de9 Merge pull request #277 from simpeg/dev
Dev
2016-05-31 15:11:29 -07:00
Lindsey Heagy 01e19e4227 start of gae site 2016-05-30 22:01:15 -07:00
Lindsey Heagy 9fbdaaf0a5 fix typo in examples path 2016-05-30 21:12:31 -07:00
Lindsey Heagy 2a159e20b9 update docs for examples to point to the correct path 2016-05-30 20:34:10 -07:00
Lindsey Heagy 5e8d3fbc78 organizing the docs - put the content in a content folder. put the SimPEG core api docs in core_api 2016-05-30 17:06:29 -07:00
Lindsey Heagy 414418a996 Merge branch 'dev' into feat/docs-deploy
# Conflicts:
#	SimPEG/Mesh/View.py
2016-05-30 15:56:24 -07:00
Rowan Cockett 40f0874dfb Doc testing, I think that is most of them! 2016-05-29 22:18:22 -07:00
Rowan Cockett 231e6dbc93 Suppress image warning. 2016-05-29 19:22:15 -07:00
Rowan Cockett 18f98b2ecd Merge branch 'docs' of https://github.com/simpeg/simpeg into feat/docs-deploy
# Conflicts:
#	SimPEG/Utils/meshutils.py
#	docs/api_Utils.rst
#	docs/conf.py
#	docs/flow/index.rst
2016-05-29 19:18:36 -07:00
Rowan Cockett bc073e49b5 Updates to docs errors. 2016-05-29 18:57:38 -07:00
Rowan Cockett a131383dae Correct solver location. 2016-05-29 18:35:41 -07:00
Rowan Cockett ad4a0240d1 Remove Vertical1DMap from tests. 2016-05-29 18:34:17 -07:00
Rowan Cockett 74f5395573 Surject1D updates. 2016-05-29 18:25:42 -07:00
Rowan Cockett 6f7a0b1279 Rename Vertical1DMap to SurjectVertical1D due to depreciation. 2016-05-29 18:14:42 -07:00
Rowan Cockett e2bb9c8d8e rename flow example. 2016-05-29 18:12:52 -07:00
Rowan Cockett de693adaa7 Minor updates to get it 'working'
There still seems to be a problem with this example:

	- The line search breaks.
	- The plots are not informative.
	- There are a lot of errors in the structured array codes.
2016-05-29 18:04:09 -07:00
Rowan Cockett fc07993006 Spacings in functions. 2016-05-29 18:03:03 -07:00
Rowan Cockett 12a12c7b5a updates to FDEM docs. 2016-05-29 17:51:41 -07:00
Rowan Cockett 0b4215f33e Add DC and IP docs. 2016-05-29 17:45:49 -07:00
Rowan Cockett 5e2a8232a3 Minor updates to titles in examples. 2016-05-29 17:21:37 -07:00
Rowan Cockett feba384911 Add solver parameter to the Casing example. 2016-05-29 17:08:11 -07:00
Rowan Cockett 4844b7230a TOC updates for docs index. 2016-05-29 17:05:42 -07:00
micmitch 1960b52dfd First stab at analytic functions for the fields from a harmonic electric dipole source. Not sure about the exception that I try to throw if multiple frequencies and multiple evaluation locations are both specified. 2016-05-27 14:56:48 -07:00
seogi_macbook 28d67e3112 Start of ED !! 2016-05-27 11:26:47 -07:00
Lindsey Heagy 756b738ef2 Merge branch 'dev' into docs 2016-04-05 17:47:17 -07:00
Lindsey Heagy 021e7c794c Merge branch 'master' into docs
# Conflicts:
#	SimPEG/Mesh/TensorMesh.py
2016-03-29 22:53:16 -07:00
Lindsey Heagy b5b70390cb Merge branch 'dev' into docs 2016-03-06 21:50:27 -08:00
Lindsey Heagy 6acaa81faf fixed merge conflicts in FDEM docs that I missed 2016-02-09 09:03:20 -08:00
Lindsey Heagy 312b5d79c5 resolved merge conflicts in TensorMesh 2016-02-09 08:53:01 -08:00
Lindsey Heagy 999a37547e Merge branch 'dev' into docs
# Conflicts:
#	.travis.yml
#	SimPEG/EM/FDEM/FDEM.py
#	SimPEG/Mesh/TensorMesh.py
2016-02-09 08:32:41 -08:00
Lindsey Heagy cbe8758465 corrected scipy.sparse.csr_matrix, move size descriptions to :return: instead of :type: 2016-02-01 08:22:00 -08:00
Lindsey Heagy 6b359f49b5 docs clean-up (using autoclass is more stable than automodule) 2016-01-31 15:21:46 -08:00
Lindsey Heagy 2254eedbac indentations clean up in FDEM.py 2016-01-31 13:54:39 -08:00
Lindsey Heagy 012d2cadf1 use intersphinx mapping to get numpy, scipy, matplotlib, python inventories 2016-01-31 13:54:24 -08:00
Lindsey Heagy 841ba61006 clean up the html build 2016-01-31 13:22:36 -08:00
Lindsey Heagy 2874e204ee exclude _static from warnings 2016-01-31 12:52:42 -08:00
Lindsey Heagy adca273565 ignore nonlocal images in sphinx build 2016-01-31 12:46:24 -08:00
Lindsey Heagy d9d6f70958 better description of paths in test_docs 2016-01-31 12:45:51 -08:00
Lindsey Heagy f4ef767764 seperate out docs test so it runs independently (not on every test) 2016-01-31 12:08:29 -08:00
Lindsey Heagy e314bdb740 add sphinx to travis conda install 2016-01-31 11:06:47 -08:00
Lindsey Heagy e005ed8f5f use cd to get into docs directories for testing 2016-01-31 10:51:55 -08:00
Lindsey Heagy ac2e38e89d test docs first 2016-01-31 10:47:00 -08:00
Lindsey Heagy ade37fb493 add travis to docs. nit-picky testing on html, latex, link check 2016-01-31 09:51:57 -08:00
270 changed files with 3847 additions and 1550 deletions
+1 -1
View File
@@ -1,4 +1,4 @@
[bumpversion]
current_version = 0.1.10
current_version = 0.1.12
files = setup.py SimPEG/__init__.py docs/conf.py
+2
View File
@@ -39,3 +39,5 @@ nosetests.xml
*.sublime-workspace
docs/_build/
Makefile
docs/warnings.txt
.DS_Store
+28 -5
View File
@@ -1,6 +1,7 @@
language: python
python:
- 2.7
- 3.4
sudo: false
@@ -24,18 +25,25 @@ env:
- TEST_DIR=tests/examples
- TEST_DIR=tests/em/fdem/inverse/adjoint
- TEST_DIR=tests/em/fdem/forward
- TEST_DIR=tests/docs;
GAE_PYTHONPATH=${HOME}/.cache/google_appengine;
PATH=$PATH:${HOME}/google-cloud-sdk/bin;
PYTHONPATH=${PYTHONPATH}:${GAE_PYTHONPATH};
CLOUDSDK_CORE_DISABLE_PROMPTS=1
# Setup anaconda
before_install:
- 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
# 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
- chmod +x miniconda.sh
- ./miniconda.sh -b
- export PATH=/home/travis/anaconda/bin:/home/travis/miniconda/bin:$PATH
- export PATH=/home/travis/anaconda/bin:/home/travis/anaconda3/bin:/home/travis/miniconda/bin:/home/travis/miniconda3/bin:$PATH
- conda update --yes conda
# Install packages
install:
- conda install --yes pip python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib cython ipython nose vtk
- conda install --yes pip python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib cython ipython nose vtk sphinx
- pip install nose-cov python-coveralls
- git clone https://github.com/rowanc1/pymatsolver.git
@@ -46,11 +54,26 @@ 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
- bash <(curl -s https://codecov.io/bash)
- 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:
+6 -2
View File
@@ -1,4 +1,4 @@
.. image:: https://raw.github.com/simpeg/simpeg/master/docs/simpeg-logo.png
.. image:: https://raw.github.com/simpeg/simpeg/master/docs/images/simpeg-logo.png
:alt: SimPEG Logo
======
@@ -28,7 +28,11 @@ SimPEG
.. image:: http://img.shields.io/badge/GITTER-JOIN_CHAT-brightgreen.svg?style=flat-square
:alt: gitter chat room at https://gitter.im/simpeg/simpeg
:target: https://gitter.im/simpeg/simpeg
.. image:: https://codecov.io/gh/simpeg/simpeg/branch/master/graph/badge.svg
  :target: https://codecov.io/gh/simpeg/simpeg
Simulation and Parameter Estimation in Geophysics - A python package for simulation and gradient based parameter estimation in the context of geophysical applications.
The vision is to create a package for finite volume simulation with applications to geophysical imaging and subsurface flow. To enable the understanding of the many different components, this package has the following features:
+11 -3
View File
@@ -1,3 +1,11 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from builtins import super
from future import standard_library
standard_library.install_aliases()
from builtins import range
from SimPEG import *
class FieldsDC_CC(Problem.Fields):
@@ -61,7 +69,7 @@ class SrcDipole(Survey.BaseSrc):
pts = [self.loc[0], self.loc[1]]
inds = Utils.closestPoints(prob.mesh, pts)
q = np.zeros(prob.mesh.nC)
q[inds] = - self.current * ( np.r_[1., -1.] / prob.mesh.vol[inds] )
q[inds] = - self.current * (np.r_[1., -1.] / prob.mesh.vol[inds])
# self._rhsDict[mesh] = q
# return self._rhsDict[mesh]
return q
@@ -162,8 +170,8 @@ class ProblemDC_CC(Problem.BaseProblem):
"""
Makes the matrix A(m) for the DC resistivity problem.
:param numpy.array m: model
:rtype: scipy.csc_matrix
:param numpy.ndarray m: model
:rtype: scipy.sparse.csc_matrix
:return: A(m)
.. math::
+10 -3
View File
@@ -1,5 +1,12 @@
from __future__ import absolute_import
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from future import standard_library
standard_library.install_aliases()
from builtins import range
from SimPEG import *
from BaseDC import SurveyDC, FieldsDC_CC
from .BaseDC import SurveyDC, FieldsDC_CC
class SurveyIP(SurveyDC):
"""
@@ -52,7 +59,7 @@ class ProblemIP(Problem.BaseProblem):
# sigma = self.curModel.transform
sigma = self.sigma
Av = self.mesh.aveF2CC
self._Msig = Utils.sdiag(1/(self.mesh.dim * Av.T * (1/sigma)))
self._Msig = Utils.sdiag(1//(self.mesh.dim * Av.T * (1/sigma)))
return self._Msig
@property
@@ -71,7 +78,7 @@ class ProblemIP(Problem.BaseProblem):
Makes the matrix A(m) for the DC resistivity problem.
:param numpy.array m: model
:rtype: scipy.csc_matrix
:rtype: scipy.sparse.csc_matrix
:return: A(m)
.. math::
+31 -21
View File
@@ -1,6 +1,16 @@
from __future__ import print_function
from __future__ import absolute_import
from __future__ import division
from __future__ import unicode_literals
from builtins import open
from builtins import int
from future import standard_library
standard_library.install_aliases()
from builtins import map
from builtins import range
from SimPEG import np, Utils
import BaseDC as DC
import BaseDC as IP
from . import BaseDC as DC
from . import BaseDC as IP
import warnings
def getActiveindfromTopo(mesh, topo):
@@ -67,7 +77,7 @@ def readUBC_DC3Dobstopo(filename,mesh,topo,probType="CC"):
if "!" in line.split(): continue
elif line == '\n': continue
elif line == ' \n': continue
temp = map(float, line.split())
temp = list(map(float, line.split()))
# Read a line for the current electrode
if len(temp) == 5: # SRC: Only X and Y are provided (assume no topography)
#TODO consider topography and assign the closest cell center in the earth
@@ -228,10 +238,10 @@ def plot_pseudoSection(DCsurvey, axs, surveyType='dipole-dipole', unitType='volt
elif surveyType == 'dipole-dipole':
leg = data * 2*np.pi / ( 1/MA - 1/MB - 1/NB + 1/NA )
leg = data * 2*np.pi / (1/MA - 1/MB - 1/NB + 1/NA)
else:
print """unitType must be 'pole-dipole' | 'dipole-dipole' """
print("""unitType must be 'pole-dipole' | 'dipole-dipole' """)
break
@@ -246,11 +256,11 @@ def plot_pseudoSection(DCsurvey, axs, surveyType='dipole-dipole', unitType='volt
rho = np.hstack([rho,leg])
else:
print """unitType must be 'appResistivity' | 'appConductivity' | 'volt' """
print("""unitType must be 'appResistivity' | 'appConductivity' | 'volt' """)
break
midx = np.hstack([midx, ( Cmid + Pmid )/2 ])
midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + (Tx[0][2] + Tx[1][2])/2 ])
midx = np.hstack([midx, (Cmid + Pmid)/2])
midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + (Tx[0][2] + Tx[1][2])/2])
# Grid points
grid_x, grid_z = np.mgrid[np.min(midx):np.max(midx), np.min(midz):np.max(midz)]
@@ -340,11 +350,11 @@ def gen_DCIPsurvey(endl, mesh, surveyType, AM_sep, MN_sep, nrx):
dl_x = ( endl[1,0] - endl[0,0] ) / dl_len
dl_y = ( endl[1,1] - endl[0,1] ) / dl_len
nstn = np.floor( dl_len / AM_sep )
nstn = np.floor(dl_len / AM_sep)
# Compute discrete pole location along line
stn_x = endl[0,0] + np.array(range(int(nstn)))*dl_x*AM_sep
stn_y = endl[0,1] + np.array(range(int(nstn)))*dl_y*AM_sep
stn_x = endl[0,0] + np.array(list(range(int(nstn))))*dl_x*AM_sep
stn_y = endl[0,1] + np.array(list(range(int(nstn))))*dl_y*AM_sep
# Create line of P1 locations
M = np.c_[stn_x, stn_y, np.ones(nstn).T*mesh.vectorNz[-1]]
@@ -376,15 +386,15 @@ def gen_DCIPsurvey(endl, mesh, surveyType, AM_sep, MN_sep, nrx):
AB = xy_2_r(tx[0,1],endl[1,0],tx[1,1],endl[1,1])
# Number of receivers to fit
nstn = np.min([np.floor( (AB - MN_sep) / AM_sep ) , nrx])
nstn = np.min([(AB - MN_sep) // AM_sep, nrx])
# Check if there is enough space, else break the loop
if nstn <= 0:
continue
# Compute discrete pole location along line
stn_x = N[ii,0] + dl_x*MN_sep + np.array(range(int(nstn)))*dl_x*AM_sep
stn_y = N[ii,1] + dl_y*MN_sep + np.array(range(int(nstn)))*dl_y*AM_sep
stn_x = N[ii,0] + dl_x*MN_sep + np.array(list(range(int(nstn))))*dl_x*AM_sep
stn_y = N[ii,1] + dl_y*MN_sep + np.array(list(range(int(nstn))))*dl_y*AM_sep
# Create receiver poles
# Create line of P1 locations
@@ -419,15 +429,15 @@ def gen_DCIPsurvey(endl, mesh, surveyType, AM_sep, MN_sep, nrx):
box_l = np.sqrt( (min_x - max_x)**2 + (min_y - max_y)**2 )
box_w = box_l/2.
nstn = np.floor( box_l / AM_sep )
nstn = box_l // AM_sep
# Compute discrete pole location along line
stn_x = min_x + np.array(range(int(nstn)))*dl_x*AM_sep
stn_y = min_y + np.array(range(int(nstn)))*dl_y*AM_sep
stn_x = min_x + np.array(list(range(int(nstn))))*dl_x*AM_sep
stn_y = min_y + np.array(list(range(int(nstn))))*dl_y*AM_sep
# Define number of cross lines
nlin = int(np.floor( box_w / AM_sep ))
lind = range(-nlin,nlin+1)
nlin = int(box_w // AM_sep)
lind = list(range(-nlin,nlin+1))
ngrad = nstn * len(lind)
@@ -449,7 +459,7 @@ def gen_DCIPsurvey(endl, mesh, surveyType, AM_sep, MN_sep, nrx):
srcClass = DC.SrcDipole([rxClass], M[0,:], N[-1,:])
SrcList.append(srcClass)
else:
print """surveyType must be either 'pole-dipole', 'dipole-dipole' or 'gradient'. """
print("""surveyType must be either 'pole-dipole', 'dipole-dipole' or 'gradient'. """)
survey = DC.SurveyDC(SrcList)
return survey, Tx, Rx
@@ -668,7 +678,7 @@ def readUBC_DC3Dobs(fileName, rtype = 'DC'):
obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!')
else:
print "rtype must be 'DC'(default) | 'IP'"
print("rtype must be 'DC'(default) | 'IP'")
# Pre-allocate
srcLists = []
+8 -1
View File
@@ -1,3 +1,10 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
import numpy as np
def WennerSrcList(nElecs, aSpacing, in2D=False, plotIt=False):
@@ -5,7 +12,7 @@ def WennerSrcList(nElecs, aSpacing, in2D=False, plotIt=False):
import SimPEG.DCIP as DC
elocs = np.arange(0,aSpacing*nElecs,aSpacing)
elocs -= (nElecs*aSpacing - aSpacing)/2
elocs -= (nElecs*aSpacing - aSpacing) / 2
space = 1
WENNER = np.zeros((0,),dtype=int)
for ii in range(nElecs):
+10 -4
View File
@@ -1,4 +1,10 @@
from BaseDC import *
from BaseIP import *
from DCIPUtils import *
import Utils
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from .BaseDC import *
from .BaseIP import *
from .DCIPUtils import *
from . import Utils
+13 -6
View File
@@ -1,7 +1,16 @@
import Utils, Survey, Problem, numpy as np, scipy.sparse as sp, gc
from __future__ import print_function
from __future__ import division
from __future__ import unicode_literals
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import object
from . import Utils, Survey, Problem
import numpy as np, scipy.sparse as sp, gc
from future.utils import with_metaclass
class BaseDataMisfit(object):
class BaseDataMisfit(with_metaclass(Utils.SimPEGMetaClass, object)):
"""BaseDataMisfit
.. note::
@@ -9,8 +18,6 @@ class BaseDataMisfit(object):
You should inherit from this class to create your own data misfit term.
"""
__metaclass__ = Utils.SimPEGMetaClass
debug = False #: Print debugging information
counter = None #: Set this to a SimPEG.Utils.Counter() if you want to count things
@@ -93,11 +100,11 @@ class l2_DataMisfit(BaseDataMisfit):
survey = self.survey
if getattr(survey,'std', None) is None:
print 'SimPEG.DataMisfit.l2_DataMisfit assigning default std of 5%'
print('SimPEG.DataMisfit.l2_DataMisfit assigning default std of 5%')
survey.std = 0.05
if getattr(survey, 'eps', None) is None:
print 'SimPEG.DataMisfit.l2_DataMisfit assigning default eps of 1e-5 * ||dobs||'
print('SimPEG.DataMisfit.l2_DataMisfit assigning default eps of 1e-5 * ||dobs||')
survey.eps = np.linalg.norm(Utils.mkvc(survey.dobs),2)*1e-5
self._Wd = Utils.sdiag(1/(abs(survey.dobs)*survey.std+survey.eps))
+35 -24
View File
@@ -1,4 +1,15 @@
import Utils, numpy as np
from __future__ import print_function
from __future__ import division
from __future__ import unicode_literals
from __future__ import absolute_import
from builtins import open
from builtins import int
from future import standard_library
standard_library.install_aliases()
from builtins import str
from builtins import object
from . import Utils
import numpy as np
class InversionDirective(object):
"""InversionDirective"""
@@ -15,7 +26,7 @@ class InversionDirective(object):
@inversion.setter
def inversion(self, i):
if getattr(self,'_inversion',None) is not None:
print 'Warning: InversionDirective %s has switched to a new inversion.' % self.__name__
print('Warning: InversionDirective %s has switched to a new inversion.' % self.__name__)
self._inversion = i
@property
@@ -68,7 +79,7 @@ class DirectiveList(object):
def inversion(self, i):
if self.inversion is i: return
if getattr(self,'_inversion',None) is not None:
print 'Warning: %s has switched to a new inversion.' % self.__name__
print('Warning: %s has switched to a new inversion.' % self.__name__)
for d in self.dList:
d.inversion = i
self._inversion = i
@@ -120,7 +131,7 @@ class BetaEstimate_ByEig(InversionDirective):
:return: beta0
"""
if self.debug: print 'Calculating the beta0 parameter.'
if self.debug: print('Calculating the beta0 parameter.')
m = self.invProb.curModel
f = self.invProb.getFields(m, store=True, deleteWarmstart=False)
@@ -141,7 +152,7 @@ class BetaSchedule(InversionDirective):
def endIter(self):
if self.opt.iter > 0 and self.opt.iter % self.coolingRate == 0:
if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter
if self.debug: print('BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter)
self.invProb.beta /= self.coolingFactor
@@ -167,7 +178,7 @@ class TargetMisfit(InversionDirective):
class _SaveEveryIteration(InversionDirective):
class SaveEveryIteration(InversionDirective):
@property
def name(self):
if getattr(self, '_name', None) is None:
@@ -188,21 +199,21 @@ class _SaveEveryIteration(InversionDirective):
self._fileName = value
class SaveModelEveryIteration(_SaveEveryIteration):
class SaveModelEveryIteration(SaveEveryIteration):
"""SaveModelEveryIteration"""
def initialize(self):
print "SimPEG.SaveModelEveryIteration will save your models as: '###-%s.npy'"%self.fileName
print("SimPEG.SaveModelEveryIteration will save your models as: '###-%s.npy'"%self.fileName)
def endIter(self):
np.save('%03d-%s' % (self.opt.iter, self.fileName), self.opt.xc)
class SaveOutputEveryIteration(_SaveEveryIteration):
class SaveOutputEveryIteration(SaveEveryIteration):
"""SaveModelEveryIteration"""
def initialize(self):
print "SimPEG.SaveOutputEveryIteration will save your inversion progress as: '###-%s.txt'"%self.fileName
print("SimPEG.SaveOutputEveryIteration will save your inversion progress as: '###-%s.txt'"%self.fileName)
f = open(self.fileName+'.txt', 'w')
f.write(" # beta phi_d phi_m f\n")
f.close()
@@ -212,11 +223,11 @@ 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):
class SaveOutputDictEveryIteration(SaveEveryIteration):
"""SaveOutputDictEveryIteration"""
def initialize(self):
print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '###-%s.npz'"%self.fileName
print("SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '###-%s.npz'"%self.fileName)
def endIter(self):
# Save the data.
@@ -294,7 +305,7 @@ class Update_IRLS(InversionDirective):
# After reaching target misfit with l2-norm, switch to IRLS (mode:2)
if self.invProb.phi_d < self.target and self.mode == 1:
print "Convergence with smooth l2-norm regularization: Start IRLS steps..."
print("Convergence with smooth l2-norm regularization: Start IRLS steps...")
self.mode = 2
@@ -302,17 +313,17 @@ class Update_IRLS(InversionDirective):
# model values
if getattr(self, 'reg.eps', None) is None:
self.reg.eps_p = np.percentile(np.abs(self.invProb.curModel),self.prctile)
else:
else:
self.reg.eps_p = self.eps[0]
if getattr(self, 'reg.eps', None) is None:
self.reg.eps_q = np.percentile(np.abs(self.reg.regmesh.cellDiffxStencil*(self.reg.mapping * self.invProb.curModel)),self.prctile)
else:
else:
self.reg.eps_q = self.eps[1]
print "L[p qx qy qz]-norm : " + str(self.reg.norms)
print "eps_p: " + str(self.reg.eps_p) + " eps_q: " + str(self.reg.eps_q)
print("L[p qx qy qz]-norm : " + str(self.reg.norms))
print("eps_p: " + str(self.reg.eps_p) + " eps_q: " + str(self.reg.eps_q))
self.reg.norms = self.norms
self.coolingFactor = 1.
self.coolingRate = 1
@@ -328,7 +339,7 @@ class Update_IRLS(InversionDirective):
# Beta Schedule
if self.opt.iter > 0 and self.opt.iter % self.coolingRate == 0:
if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter
if self.debug: print('BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter)
self.invProb.beta /= self.coolingFactor
@@ -340,17 +351,17 @@ class Update_IRLS(InversionDirective):
phim_new = self.reg.eval(self.invProb.curModel)
self.f_change = np.abs(self.f_old - phim_new) / self.f_old
print "Regularization decrease: %6.3e" % (self.f_change)
print("Regularization decrease: %6.3e" % (self.f_change))
# Check for maximum number of IRLS cycles
if self.IRLSiter == self.maxIRLSiter:
print "Reach maximum number of IRLS cycles: %i" % self.maxIRLSiter
print("Reach maximum number of IRLS cycles: %i" % self.maxIRLSiter)
self.opt.stopNextIteration = True
return
# Check if the function has changed enough
if self.f_change < self.f_min_change and self.IRLSiter > 1:
print "Minimum decrease in regularization. End of IRLS"
print("Minimum decrease in regularization. End of IRLS")
self.opt.stopNextIteration = True
return
else:
+7
View File
@@ -1,3 +1,10 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
import numpy as np
from scipy.constants import mu_0, pi
from scipy import special
+5
View File
@@ -1,4 +1,9 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
import numpy as np
from scipy.constants import mu_0, pi
from scipy.special import erf
+307
View File
@@ -0,0 +1,307 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
import numpy as np
from scipy.constants import mu_0, pi, epsilon_0
from scipy.special import erf
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
+7 -1
View File
@@ -1,3 +1,9 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import Utils, np
from scipy.constants import mu_0, epsilon_0
from SimPEG.EM.Utils.EMUtils import k
@@ -34,7 +40,7 @@ def _getCasingHertzMagDipoleDeriv_r(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones
sqrtr2z2 = np.sqrt(r2 + dxyz[:,2]**2)
k2 = k(freq,sigma[2],mu[2],eps)
return -HertzZ * np.sqrt(r2) / sqrtr2z2 * (1j*k2 + 1./ sqrtr2z2)
return -HertzZ * np.sqrt(r2) / sqrtr2z2 * (1j*k2 + 1. / sqrtr2z2)
def _getCasingHertzMagDipoleDeriv_z(srcloc,obsloc,freq,sigma,a,b,mu=mu_0*np.ones(3),eps=epsilon_0,moment=1.):
+6
View File
@@ -1,3 +1,9 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
import numpy as np
from scipy.constants import mu_0, pi
from scipy.special import erf
+11 -4
View File
@@ -1,4 +1,11 @@
from TDEM import hzAnalyticDipoleT
from FDEM import hzAnalyticDipoleF
from FDEMcasing import *
from DC import DCAnalyticHalf, DCAnalyticSphere
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from .TDEM import hzAnalyticDipoleT
from .FDEM import hzAnalyticDipoleF
from .FDEMcasing import *
from .DC import DCAnalyticHalf, DCAnalyticSphere
from .FDEMDipolarfields import *
+10 -3
View File
@@ -1,3 +1,9 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import Survey, Problem, Utils, Models, Maps, PropMaps, np, sp, Solver as SimpegSolver
from scipy.constants import mu_0
@@ -20,10 +26,10 @@ class BaseEMProblem(Problem.BaseProblem):
Problem.BaseProblem.__init__(self, mesh, **kwargs)
surveyPair = Survey.BaseSurvey
dataPair = Survey.Data
surveyPair = Survey.BaseSurvey #: The survey to pair with.
dataPair = Survey.Data #: The data to pair with.
PropMap = EMPropMap
PropMap = EMPropMap #: The property mapping
Solver = SimpegSolver
solverOpts = {}
@@ -217,6 +223,7 @@ class BaseEMSurvey(Survey.BaseSurvey):
def eval(self, f):
"""
Project fields to receiver locations
:param Fields u: fields object
:rtype: numpy.ndarray
:return: data
+54 -59
View File
@@ -1,3 +1,10 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from builtins import int
from future import standard_library
standard_library.install_aliases()
import numpy as np
import scipy.sparse as sp
import SimPEG
@@ -6,11 +13,11 @@ from SimPEG.EM.Utils import omega
from SimPEG.Utils import Zero, Identity, sdiag
class Fields(SimPEG.Problem.Fields):
class FieldsFDEM(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 obejct acts like an array and is indexed by
each problem, the rest are computed. The fields object acts like an array and is indexed by
.. code-block:: python
@@ -42,7 +49,7 @@ class Fields(SimPEG.Problem.Fields):
:return: total electric field
"""
if getattr(self, '_ePrimary', None) is None or getattr(self, '_eSecondary', None) is None:
raise NotImplementedError ('Getting e from %s is not implemented' %self.knownFields.keys()[0])
raise NotImplementedError ('Getting e from %s is not implemented' %list(self.knownFields.keys())[0])
return self._ePrimary(solution,srcList) + self._eSecondary(solution,srcList)
@@ -56,7 +63,7 @@ class Fields(SimPEG.Problem.Fields):
:return: total magnetic flux density
"""
if getattr(self, '_bPrimary', None) is None or getattr(self, '_bSecondary', None) is None:
raise NotImplementedError ('Getting b from %s is not implemented' %self.knownFields.keys()[0])
raise NotImplementedError ('Getting b from %s is not implemented' %list(self.knownFields.keys())[0])
return self._bPrimary(solution, srcList) + self._bSecondary(solution, srcList)
@@ -70,7 +77,7 @@ class Fields(SimPEG.Problem.Fields):
:return: total magnetic field
"""
if getattr(self, '_hPrimary', None) is None or getattr(self, '_hSecondary', None) is None:
raise NotImplementedError ('Getting h from %s is not implemented' %self.knownFields.keys()[0])
raise NotImplementedError ('Getting h from %s is not implemented' %list(self.knownFields.keys())[0])
return self._hPrimary(solution, srcList) + self._hSecondary(solution, srcList)
@@ -84,7 +91,7 @@ class Fields(SimPEG.Problem.Fields):
:return: total current density
"""
if getattr(self, '_jPrimary', None) is None or getattr(self, '_jSecondary', None) is None:
raise NotImplementedError ('Getting j from %s is not implemented' %self.knownFields.keys()[0])
raise NotImplementedError ('Getting j from %s is not implemented' %list(self.knownFields.keys())[0])
return self._jPrimary(solution, srcList) + self._jSecondary(solution, srcList)
@@ -92,7 +99,7 @@ class Fields(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 Src src: sorce
:param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: source
: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?
@@ -100,7 +107,7 @@ class Fields(SimPEG.Problem.Fields):
:return: derivative times a vector (or tuple for adjoint)
"""
if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None:
raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0])
raise NotImplementedError ('Getting eDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
if adjoint:
return self._eDeriv_u(src, v, adjoint), self._eDeriv_m(src, v, adjoint)
@@ -110,7 +117,7 @@ class Fields(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 Src src: sorce
:param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: source
: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?
@@ -118,7 +125,7 @@ class Fields(SimPEG.Problem.Fields):
:return: derivative times a vector (or tuple for adjoint)
"""
if getattr(self, '_bDeriv_u', None) is None or getattr(self, '_bDeriv_m', None) is None:
raise NotImplementedError ('Getting bDerivs from %s is not implemented' %self.knownFields.keys()[0])
raise NotImplementedError ('Getting bDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
if adjoint:
return self._bDeriv_u(src, v, adjoint), self._bDeriv_m(src, v, adjoint)
@@ -128,7 +135,7 @@ class Fields(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 Src src: sorce
:param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: source
: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?
@@ -136,7 +143,7 @@ class Fields(SimPEG.Problem.Fields):
:return: derivative times a vector (or tuple for adjoint)
"""
if getattr(self, '_hDeriv_u', None) is None or getattr(self, '_hDeriv_m', None) is None:
raise NotImplementedError ('Getting hDerivs from %s is not implemented' %self.knownFields.keys()[0])
raise NotImplementedError ('Getting hDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
if adjoint:
return self._hDeriv_u(src, v, adjoint), self._hDeriv_m(src, v, adjoint)
@@ -146,7 +153,7 @@ class Fields(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 Src src: sorce
:param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: source
: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?
@@ -154,18 +161,18 @@ class Fields(SimPEG.Problem.Fields):
:return: derivative times a vector (or tuple for adjoint)
"""
if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None:
raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0])
raise NotImplementedError ('Getting jDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
if adjoint:
return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint)
return np.array(self._jDeriv_u(src, du_dm_v, adjoint) + self._jDeriv_m(src, v, adjoint), dtype = complex)
class Fields3D_e(Fields):
class Fields3D_e(FieldsFDEM):
"""
Fields object for Problem3D_e.
:param Mesh mesh: mesh
:param Survey survey: survey
:param BaseMesh mesh: mesh
:param SimPEG.EM.FDEM.SurveyFDEM.Survey survey: survey
"""
knownFields = {'eSolution':'E'}
@@ -180,9 +187,6 @@ class Fields3D_e(Fields):
'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
@@ -288,7 +292,7 @@ class Fields3D_e(Fields):
C = self._edgeCurl
b = (C * eSolution)
for i, src in enumerate(srcList):
b[:,i] *= - 1./(1j*omega(src.freq))
b[:,i] *= -1./(1j*omega(src.freq))
s_m, _ = src.eval(self.prob)
b[:,i] = b[:,i]+ 1./(1j*omega(src.freq)) * s_m
return b
@@ -348,7 +352,7 @@ class Fields3D_e(Fields):
:rtype: numpy.ndarray
:return: product of the derivative of the current density with respect to the field we solved for with a vector
"""
n = int(self._aveE2CCV.shape[0] / self._nC) # number of components (instead of checking if cyl or not)
n = int(self._aveE2CCV.shape[0] // self._nC) # number of components (instead of checking if cyl or not)
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
@@ -385,8 +389,8 @@ class Fields3D_e(Fields):
:rtype: numpy.ndarray
:return: magnetic field
"""
n = int(self._aveF2CCV.shape[0] / self._nC) # Number of Components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
n = int(self._aveF2CCV.shape[0] // self._nC) # Number of Components
VI = sdiag(np.kron(np.ones(n), 1. // self.prob.mesh.vol))
return VI * (self._aveF2CCV * (self._MfMui * self._b(eSolution, srcList)))
@@ -400,7 +404,7 @@ class Fields3D_e(Fields):
:rtype: numpy.ndarray
:return: product of the derivative of the magnetic field with respect to the field we solved for with a vector
"""
n = int(self._aveF2CCV.shape[0] / self._nC) # Number of Components
n = int(self._aveF2CCV.shape[0] // self._nC) # Number of Components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
v = self._MfMui.T * (self._aveF2CCV.T * (VI.T * du_dm_v))
@@ -417,7 +421,7 @@ class Fields3D_e(Fields):
:rtype: numpy.ndarray
:return: product of the magnetic field derivative with respect to the inversion model with a vector
"""
n = int(self._aveF2CCV.shape[0] / self._nC) # Number of Components
n = int(self._aveF2CCV.shape[0] // self._nC) # Number of Components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
v = self._MfMui.T * (self._aveF2CCV.T * (VI.T * v))
@@ -426,12 +430,12 @@ class Fields3D_e(Fields):
class Fields3D_b(Fields):
class Fields3D_b(FieldsFDEM):
"""
Fields object for Problem3D_b.
:param Mesh mesh: mesh
:param Survey survey: survey
:param BaseMesh mesh: mesh
:param SimPEG.EM.FDEM.SurveyFDEM.Survey survey: survey
"""
knownFields = {'bSolution':'F'}
@@ -446,9 +450,6 @@ class Fields3D_b(Fields):
'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
@@ -613,7 +614,7 @@ class Fields3D_b(Fields):
:return: primary current density
"""
n = int(self._aveE2CCV.shape[0] / self._nC) # number of components
n = int(self._aveE2CCV.shape[0] // self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (self._aveE2CCV * ( self._MeSigma * self._e(bSolution,srcList ) ) )
@@ -630,7 +631,7 @@ class Fields3D_b(Fields):
:rtype: numpy.ndarray
:return: product of the derivative of the current density with respect to the field we solved for with a vector
"""
n = int(self._aveE2CCV.shape[0] / self._nC) # number of components
n = int(self._aveE2CCV.shape[0] // self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
return self._MfMui.T * ( self._edgeCurl * ( self._aveE2CCV.T * (VI.T * du_dm_v) ) )
@@ -658,7 +659,7 @@ class Fields3D_b(Fields):
:rtype: numpy.ndarray
:return: magnetic field
"""
n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
n = int(self._aveF2CCV.shape[0] // self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (self._aveF2CCV * (self._MfMui * self._b(bSolution, srcList)))
@@ -673,7 +674,7 @@ class Fields3D_b(Fields):
:rtype: numpy.ndarray
:return: product of the derivative of the magnetic field with respect to the field we solved for with a vector
"""
n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
n = int(self._aveF2CCV.shape[0] // self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
@@ -693,12 +694,12 @@ class Fields3D_b(Fields):
return Zero()
class Fields3D_j(Fields):
class Fields3D_j(FieldsFDEM):
"""
Fields object for Problem3D_j.
:param Mesh mesh: mesh
:param Survey survey: survey
:param BaseMesh mesh: mesh
:param SimPEG.EM.FDEM.SurveyFDEM.Survey survey: survey
"""
knownFields = {'jSolution':'F'}
@@ -713,9 +714,6 @@ class Fields3D_j(Fields):
'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
@@ -899,7 +897,7 @@ class Fields3D_j(Fields):
:rtype: numpy.ndarray
:return: electric field
"""
n = int(self._aveF2CCV.shape[0] / self._nC) # number of components
n = int(self._aveF2CCV.shape[0] // self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (self._aveF2CCV * (self._MfRho * self._j(jSolution, srcList)))
@@ -913,7 +911,7 @@ class Fields3D_j(Fields):
:rtype: numpy.ndarray
:return: product of the derivative of the electric field with respect to the field we solved for with a vector
"""
n = int(self._aveF2CCV.shape[0] / self._nC) # number of components
n = int(self._aveF2CCV.shape[0] // self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
return self._MfRho.T * ( self._aveF2CCV.T * ( VI.T * du_dm_v ) )
@@ -930,7 +928,7 @@ class Fields3D_j(Fields):
:return: product of the derivative of the electric field with respect to the model with a vector
"""
jSolution = Utils.mkvc(self[src,'jSolution'])
n = int(self._aveF2CCV.shape[0] / self._nC) # number of components
n = int(self._aveF2CCV.shape[0] // self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
return self._MfRhoDeriv(jSolution).T * ( self._aveF2CCV.T * ( VI.T * v ) )
@@ -945,7 +943,7 @@ class Fields3D_j(Fields):
:rtype: numpy.ndarray
:return: secondary magnetic flux density
"""
n = int(self._aveE2CCV.shape[0] / self._nC) # number of components
n = int(self._aveE2CCV.shape[0] // self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (self._aveE2CCV * ( self._MeMu * self._h(jSolution,srcList)) )
@@ -960,7 +958,7 @@ class Fields3D_j(Fields):
:rtype: numpy.ndarray
:return: product of the derivative of the magnetic flux density with respect to the field we solved for with a vector
"""
n = int(self._aveF2CCV.shape[0] / self._nC) # number of components
n = int(self._aveF2CCV.shape[0] // self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
@@ -978,7 +976,7 @@ class Fields3D_j(Fields):
:return: product of the derivative of the magnetic flux density with respect to the model with a vector
"""
jSolution = self[src,'jSolution']
n = int(self._aveE2CCV.shape[0] / self._nC) # number of components
n = int(self._aveE2CCV.shape[0] // self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
s_mDeriv,_ = src.evalDeriv(self.prob, adjoint = adjoint)
@@ -988,12 +986,12 @@ class Fields3D_j(Fields):
return 1./(1j * omega(src.freq)) * VI * (self._aveE2CCV * ( s_mDeriv(v) - self._edgeCurl.T * ( self._MfRhoDeriv(jSolution) * v ) ) )
class Fields3D_h(Fields):
class Fields3D_h(FieldsFDEM):
"""
Fields object for Problem3D_h.
:param Mesh mesh: mesh
:param Survey survey: survey
:param BaseMesh mesh: mesh
:param SimPEG.EM.FDEM.SurveyFDEM.Survey survey: survey
"""
knownFields = {'hSolution':'E'}
@@ -1008,9 +1006,6 @@ class Fields3D_h(Fields):
'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
@@ -1163,7 +1158,7 @@ class Fields3D_h(Fields):
:rtype: numpy.ndarray
:return: electric field
"""
n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
n = int(self._aveF2CCV.shape[0] // self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (self._aveF2CCV * (self._MfRho * self._j(hSolution, srcList)))
@@ -1177,7 +1172,7 @@ class Fields3D_h(Fields):
:rtype: numpy.ndarray
:return: product of the derivative of the electric field with respect to the field we solved for with a vector
"""
n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
n = int(self._aveF2CCV.shape[0] // self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
return self._edgeCurl.T * ( self._MfRho.T * ( self._aveF2CCV.T * ( VI.T * du_dm_v ) ) )
@@ -1194,7 +1189,7 @@ class Fields3D_h(Fields):
:return: product of the electric field derivative with respect to the inversion model with a vector
"""
hSolution = Utils.mkvc(self[src,'hSolution'])
n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
n = int(self._aveF2CCV.shape[0] // self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
return ( self._MfRhoDeriv(self._edgeCurl * hSolution).T * ( self._aveF2CCV.T * (VI.T * v) ) )
@@ -1210,7 +1205,7 @@ class Fields3D_h(Fields):
:return: magnetic flux density
"""
h = self._h(hSolution, srcList)
n = int(self._aveE2CCV.shape[0] / self._nC) #number of components
n = int(self._aveE2CCV.shape[0] // self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (self._aveE2CCV * (self._MeMu * h))
@@ -1225,7 +1220,7 @@ class Fields3D_h(Fields):
:rtype: numpy.ndarray
:return: product of the derivative of the magnetic flux density with respect to the field we solved for with a vector
"""
n = int(self._aveE2CCV.shape[0] / self._nC) #number of components
n = int(self._aveE2CCV.shape[0] // self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
return self._MeMu.T * (self._aveE2CCV.T * ( VI.T * du_dm_v ))
+28 -17
View File
@@ -1,7 +1,13 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from SimPEG import Problem, Utils, np, sp, Solver as SimpegSolver
from scipy.constants import mu_0
from SurveyFDEM import Survey as SurveyFDEM
from FieldsFDEM import Fields, Fields3D_e, Fields3D_b, Fields3D_h, Fields3D_j
from .SurveyFDEM import Survey as SurveyFDEM
from .FieldsFDEM import FieldsFDEM, Fields3D_e, Fields3D_b, Fields3D_h, Fields3D_j
from SimPEG.EM.Base import BaseEMProblem
from SimPEG.EM.Utils import omega
@@ -31,10 +37,11 @@ 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 = Fields
fieldsPair = FieldsFDEM
def fields(self, m):
"""
@@ -64,7 +71,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.Fields u: fields object
:param SimPEG.EM.FDEM.FieldsFDEM.FieldsFDEM u: fields object
:rtype numpy.array:
:return: Jv (ndata,)
"""
@@ -99,7 +106,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.Fields u: fields object
:param SimPEG.EM.FDEM.FieldsFDEM.FieldsFDEM u: fields object
:rtype numpy.array:
:return: Jv (ndata,)
"""
@@ -153,8 +160,8 @@ class BaseFDEMProblem(BaseEMProblem):
Evaluates the sources for a given frequency and puts them in matrix form
:param float freq: Frequency
:rtype: (numpy.ndarray, numpy.ndarray)
:return: s_m, s_e (nE or nF, nSrc)
:rtype: tuple
:return: (s_m, s_e) (nE or nF, nSrc)
"""
Srcs = self.survey.getSrcByFreq(freq)
if self._formulation is 'EB':
@@ -194,7 +201,7 @@ class Problem3D_e(BaseFDEMProblem):
which we solve for :math:`\mathbf{e}`.
:param SimPEG.Mesh mesh: mesh
:param SimPEG.Mesh.BaseMesh.BaseMesh mesh: mesh
"""
_solutionType = 'eSolution'
@@ -269,7 +276,7 @@ class Problem3D_e(BaseFDEMProblem):
Derivative of the right hand side with respect to the model
:param float freq: frequency
:param SimPEG.EM.FDEM.Src src: FDEM source
:param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
@@ -305,7 +312,7 @@ class Problem3D_b(BaseFDEMProblem):
.. note ::
The inverse problem will not work with full anisotropy
:param SimPEG.Mesh mesh: mesh
:param SimPEG.Mesh.BaseMesh.BaseMesh mesh: mesh
"""
_solutionType = 'bSolution'
@@ -400,7 +407,7 @@ class Problem3D_b(BaseFDEMProblem):
Derivative of the right hand side with respect to the model
:param float freq: frequency
:param SimPEG.EM.FDEM.Src src: FDEM source
:param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
@@ -444,6 +451,7 @@ 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 ::
@@ -453,7 +461,7 @@ class Problem3D_j(BaseFDEMProblem):
.. note::
This implementation does not yet work with full anisotropy!!
:param SimPEG.Mesh mesh: mesh
:param SimPEG.Mesh.BaseMesh.BaseMesh mesh: mesh
"""
_solutionType = 'jSolution'
@@ -529,8 +537,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 (nE, nSrc)
:return: RHS
:rtype: numpy.ndarray
:return: RHS (nE, nSrc)
"""
s_m, s_e = self.getSourceTerm(freq)
@@ -549,7 +557,7 @@ class Problem3D_j(BaseFDEMProblem):
Derivative of the right hand side with respect to the model
:param float freq: frequency
:param SimPEG.EM.FDEM.Src src: FDEM source
:param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
@@ -591,7 +599,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 mesh: mesh
:param SimPEG.Mesh.BaseMesh.BaseMesh mesh: mesh
"""
_solutionType = 'hSolution'
@@ -608,9 +616,11 @@ 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
@@ -653,6 +663,7 @@ class Problem3D_h(BaseFDEMProblem):
:param float freq: Frequency
:rtype: numpy.ndarray
:return: RHS (nE, nSrc)
"""
s_m, s_e = self.getSourceTerm(freq)
@@ -666,7 +677,7 @@ class Problem3D_h(BaseFDEMProblem):
Derivative of the right hand side with respect to the model
:param float freq: frequency
:param SimPEG.EM.FDEM.Src src: FDEM source
:param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
+12 -5
View File
@@ -1,3 +1,10 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from builtins import super
from future import standard_library
standard_library.install_aliases()
import SimPEG
from SimPEG import sp
@@ -25,10 +32,10 @@ class BaseRx(SimPEG.Survey.BaseRx):
def eval(self, src, mesh, f):
"""
Project fields to recievers to get data.
Project fields to receivers to get data.
:param Source src: FDEM source
:param Mesh mesh: mesh used
:param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
:param BaseMesh mesh: mesh used
:param Fields f: fields object
:rtype: numpy.ndarray
:return: fields projected to recievers
@@ -44,8 +51,8 @@ class BaseRx(SimPEG.Survey.BaseRx):
"""
Derivative of projected fields with respect to the inversion model times a vector.
:param Source src: FDEM source
:param Mesh mesh: mesh used
:param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
:param BaseMesh mesh: mesh used
:param Fields f: fields object
:param numpy.ndarray v: vector to multiply
:rtype: numpy.ndarray
+34 -28
View File
@@ -1,3 +1,9 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import Survey, Problem, Utils, np, sp
from scipy.constants import mu_0
from SimPEG.EM.Utils import *
@@ -23,8 +29,8 @@ class BaseSrc(Survey.BaseSrc):
- :math:`s_m` : magnetic source term
- :math:`s_e` : electric source term
:param Problem prob: FDEM Problem
:rtype: (numpy.ndarray, numpy.ndarray)
:param BaseFDEMProblem prob: FDEM Problem
:rtype: tuple
:return: tuple with magnetic source term and electric source term
"""
s_m = self.s_m(prob)
@@ -37,10 +43,10 @@ class BaseSrc(Survey.BaseSrc):
- :code:`s_mDeriv` : derivative of the magnetic source term
- :code:`s_eDeriv` : derivative of the electric source term
:param Problem prob: FDEM Problem
:param BaseFDEMProblem prob: FDEM Problem
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: (numpy.ndarray, numpy.ndarray)
:rtype: tuple
:return: tuple with magnetic source term and electric source term derivatives times a vector
"""
if v is not None:
@@ -52,7 +58,7 @@ class BaseSrc(Survey.BaseSrc):
"""
Primary magnetic flux density
:param Problem prob: FDEM Problem
:param BaseFDEMProblem prob: FDEM Problem
:rtype: numpy.ndarray
:return: primary magnetic flux density
"""
@@ -64,7 +70,7 @@ class BaseSrc(Survey.BaseSrc):
"""
Primary magnetic field
:param Problem prob: FDEM Problem
:param BaseFDEMProblem prob: FDEM Problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
@@ -76,7 +82,7 @@ class BaseSrc(Survey.BaseSrc):
"""
Primary electric field
:param Problem prob: FDEM Problem
:param BaseFDEMProblem prob: FDEM Problem
:rtype: numpy.ndarray
:return: primary electric field
"""
@@ -88,7 +94,7 @@ class BaseSrc(Survey.BaseSrc):
"""
Primary current density
:param Problem prob: FDEM Problem
:param BaseFDEMProblem prob: FDEM Problem
:rtype: numpy.ndarray
:return: primary current density
"""
@@ -100,7 +106,7 @@ class BaseSrc(Survey.BaseSrc):
"""
Magnetic source term
:param Problem prob: FDEM Problem
:param BaseFDEMProblem prob: FDEM Problem
:rtype: numpy.ndarray
:return: magnetic source term on mesh
"""
@@ -110,7 +116,7 @@ class BaseSrc(Survey.BaseSrc):
"""
Electric source term
:param Problem prob: FDEM Problem
:param BaseFDEMProblem prob: FDEM Problem
:rtype: numpy.ndarray
:return: electric source term on mesh
"""
@@ -120,7 +126,7 @@ class BaseSrc(Survey.BaseSrc):
"""
Derivative of magnetic source term with respect to the inversion model
:param Problem prob: FDEM Problem
:param BaseFDEMProblem prob: FDEM Problem
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
@@ -133,7 +139,7 @@ class BaseSrc(Survey.BaseSrc):
"""
Derivative of electric source term with respect to the inversion model
:param Problem prob: FDEM Problem
:param BaseFDEMProblem prob: FDEM Problem
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: numpy.ndarray
@@ -162,7 +168,7 @@ class RawVec_e(BaseSrc):
"""
Electric source term
:param Problem prob: FDEM Problem
:param BaseFDEMProblem prob: FDEM Problem
:rtype: numpy.ndarray
:return: electric source term on mesh
"""
@@ -191,7 +197,7 @@ class RawVec_m(BaseSrc):
"""
Magnetic source term
:param Problem prob: FDEM Problem
:param BaseFDEMProblem prob: FDEM Problem
:rtype: numpy.ndarray
:return: magnetic source term on mesh
"""
@@ -220,7 +226,7 @@ class RawVec(BaseSrc):
"""
Magnetic source term
:param Problem prob: FDEM Problem
:param BaseFDEMProblem prob: FDEM Problem
:rtype: numpy.ndarray
:return: magnetic source term on mesh
"""
@@ -232,7 +238,7 @@ class RawVec(BaseSrc):
"""
Electric source term
:param Problem prob: FDEM Problem
:param BaseFDEMProblem prob: FDEM Problem
:rtype: numpy.ndarray
:return: electric source term on mesh
"""
@@ -301,7 +307,7 @@ class MagDipole(BaseSrc):
"""
The primary magnetic flux density from a magnetic vector potential
:param Problem prob: FDEM problem
:param BaseFDEMProblem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
@@ -339,7 +345,7 @@ class MagDipole(BaseSrc):
"""
The primary magnetic field from a magnetic vector potential
:param Problem prob: FDEM problem
:param BaseFDEMProblem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
@@ -350,7 +356,7 @@ class MagDipole(BaseSrc):
"""
The magnetic source term
:param Problem prob: FDEM problem
:param BaseFDEMProblem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
@@ -364,7 +370,7 @@ class MagDipole(BaseSrc):
"""
The electric source term
:param Problem prob: FDEM problem
:param BaseFDEMProblem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
@@ -416,7 +422,7 @@ class MagDipole_Bfield(BaseSrc):
"""
The primary magnetic flux density from the analytic solution for magnetic fields from a dipole
:param Problem prob: FDEM problem
:param BaseFDEMProblem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
@@ -455,7 +461,7 @@ class MagDipole_Bfield(BaseSrc):
"""
The primary magnetic field from a magnetic vector potential
:param Problem prob: FDEM problem
:param BaseFDEMProblem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
@@ -466,7 +472,7 @@ class MagDipole_Bfield(BaseSrc):
"""
The magnetic source term
:param Problem prob: FDEM problem
:param BaseFDEMProblem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
@@ -479,7 +485,7 @@ class MagDipole_Bfield(BaseSrc):
"""
The electric source term
:param Problem prob: FDEM problem
:param BaseFDEMProblem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
@@ -530,7 +536,7 @@ class CircularLoop(BaseSrc):
"""
The primary magnetic flux density from a magnetic vector potential
:param Problem prob: FDEM problem
:param BaseFDEMProblem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
@@ -567,7 +573,7 @@ class CircularLoop(BaseSrc):
"""
The primary magnetic field from a magnetic vector potential
:param Problem prob: FDEM problem
:param BaseFDEMProblem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
@@ -578,7 +584,7 @@ class CircularLoop(BaseSrc):
"""
The magnetic source term
:param Problem prob: FDEM problem
:param BaseFDEMProblem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
@@ -591,7 +597,7 @@ class CircularLoop(BaseSrc):
"""
The electric source term
:param Problem prob: FDEM problem
:param BaseFDEMProblem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
"""
+8 -2
View File
@@ -1,10 +1,16 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
import SimPEG
from SimPEG.EM.Utils import *
from SimPEG.EM.Base import BaseEMSurvey
from scipy.constants import mu_0
from SimPEG.Utils import Zero, Identity
import SrcFDEM as Src
import RxFDEM as Rx
from . import SrcFDEM as Src
from . import RxFDEM as Rx
from SimPEG import sp
class Survey(BaseEMSurvey):
+11 -5
View File
@@ -1,5 +1,11 @@
from SurveyFDEM import Survey
import SrcFDEM as Src
import RxFDEM as Rx
from ProblemFDEM import Problem3D_e, Problem3D_b, Problem3D_j, Problem3D_h
from FieldsFDEM import Fields3D_e, Fields3D_b, Fields3D_j, Fields3D_h
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from .SurveyFDEM import Survey
from . import SrcFDEM as Src
from . import RxFDEM as Rx
from .ProblemFDEM import Problem3D_e, Problem3D_b, Problem3D_j, Problem3D_h
from .FieldsFDEM import Fields3D_e, Fields3D_b, Fields3D_j, Fields3D_h
+6
View File
@@ -1,3 +1,9 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
import numpy as np
def getxBCyBC_CC(mesh, alpha, beta, gamma):
+9 -3
View File
@@ -1,3 +1,9 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
import SimPEG
from SimPEG.Utils import Identity, Zero
import numpy as np
@@ -9,7 +15,7 @@ class Fields(SimPEG.Problem.Fields):
def _phiDeriv(self, src, du_dm_v, v, adjoint=False):
if getattr(self, '_phiDeriv_u', None) is None or getattr(self, '_phiDeriv_m', None) is None:
raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %self.knownFields.keys()[0])
raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
if adjoint:
return self._phiDeriv_u(src, v, adjoint=adjoint), self._phiDeriv_m(src, v, adjoint=adjoint)
@@ -18,7 +24,7 @@ class Fields(SimPEG.Problem.Fields):
def _eDeriv(self, src, du_dm_v, v, adjoint=False):
if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None:
raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0])
raise NotImplementedError ('Getting eDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
if adjoint:
return self._eDeriv_u(src, v, adjoint), self._eDeriv_m(src, v, adjoint)
@@ -26,7 +32,7 @@ class Fields(SimPEG.Problem.Fields):
def _jDeriv(self, src, du_dm_v, v, adjoint=False):
if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None:
raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0])
raise NotImplementedError ('Getting jDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
if adjoint:
return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint)
+9 -3
View File
@@ -1,3 +1,9 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
import SimPEG
from SimPEG.Utils import Identity, Zero
import numpy as np
@@ -32,7 +38,7 @@ class Fields_ky(SimPEG.Problem.TimeFields):
def _phiDeriv(self,kyInd, src, du_dm_v, v, adjoint=False):
if getattr(self, '_phiDeriv_u', None) is None or getattr(self, '_phiDeriv_m', None) is None:
raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %self.knownFields.keys()[0])
raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
if adjoint:
return self._phiDeriv_u(kyInd, src, v, adjoint=adjoint), self._phiDeriv_m(kyInd, src, v, adjoint=adjoint)
@@ -41,7 +47,7 @@ class Fields_ky(SimPEG.Problem.TimeFields):
def _eDeriv(self,kyInd, src, du_dm_v, v, adjoint=False):
if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None:
raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0])
raise NotImplementedError ('Getting eDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
if adjoint:
return self._eDeriv_u(kyInd, src, v, adjoint), self._eDeriv_m(kyInd, src, v, adjoint)
@@ -49,7 +55,7 @@ class Fields_ky(SimPEG.Problem.TimeFields):
def _jDeriv(self,kyInd, src, du_dm_v, v, adjoint=False):
if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None:
raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0])
raise NotImplementedError ('Getting jDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
if adjoint:
return self._jDeriv_u(kyInd, src, v, adjoint), self._jDeriv_m(kyInd, src, v, adjoint)
+9 -3
View File
@@ -1,11 +1,17 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from SimPEG import Problem, Utils
from SimPEG.EM.Base import BaseEMProblem
from SurveyDC import Survey
from FieldsDC import Fields, Fields_CC, Fields_N
from .SurveyDC import Survey
from .FieldsDC import Fields, Fields_CC, Fields_N
from SimPEG.Utils import sdiag
import numpy as np
from SimPEG.Utils import Zero
from BoundaryUtils import getxBCyBC_CC
from .BoundaryUtils import getxBCyBC_CC
class BaseDCProblem(BaseEMProblem):
+10 -3
View File
@@ -1,11 +1,18 @@
from __future__ import absolute_import
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from future import standard_library
standard_library.install_aliases()
from builtins import range
from SimPEG import Problem, Utils
from SimPEG.EM.Base import BaseEMProblem
from SurveyDC import Survey, Survey_ky
from FieldsDC_2D import Fields_ky, Fields_ky_CC, Fields_ky_N
from .SurveyDC import Survey, Survey_ky
from .FieldsDC_2D import Fields_ky, Fields_ky_CC, Fields_ky_N
from SimPEG.Utils import sdiag
import numpy as np
from SimPEG.Utils import Zero
from BoundaryUtils import getxBCyBC_CC
from .BoundaryUtils import getxBCyBC_CC
class BaseDCProblem_2D(BaseEMProblem):
+7
View File
@@ -1,3 +1,10 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
import SimPEG
import numpy as np
from SimPEG.Utils import Zero, closestPoints
+6
View File
@@ -1,3 +1,9 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
import SimPEG
# from SimPEG.EM.Base import BaseEMSurvey
from SimPEG.Utils import Zero, closestPoints, mkvc
+8 -2
View File
@@ -1,9 +1,15 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
import SimPEG
from SimPEG.EM.Base import BaseEMSurvey
from SimPEG import sp, Survey
from SimPEG.Utils import Zero, Identity
from RxDC import BaseRx
from SrcDC import BaseSrc
from .RxDC import BaseRx
from .SrcDC import BaseSrc
class Survey(BaseEMSurvey):
rxPair = BaseRx
+7
View File
@@ -1,3 +1,10 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
import numpy as np
def WennerSrcList(nElecs, aSpacing, in2D=False, plotIt=False):
+14 -8
View File
@@ -1,8 +1,14 @@
from ProblemDC import Problem3D_CC, Problem3D_N
from ProblemDC_2D import Problem2D_CC, Problem2D_N
from SurveyDC import Survey, Survey_ky
import SrcDC as Src #Pole
import RxDC as Rx
from FieldsDC import Fields_CC
from BoundaryUtils import getxBCyBC_CC
import Utils
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from .ProblemDC import Problem3D_CC, Problem3D_N
from .ProblemDC_2D import Problem2D_CC, Problem2D_N
from .SurveyDC import Survey, Survey_ky
from . import SrcDC as Src #Pole
from . import RxDC as Rx
from .FieldsDC import Fields_CC
from .BoundaryUtils import getxBCyBC_CC
from . import Utils
+7 -1
View File
@@ -1,3 +1,9 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from SimPEG import Problem, Utils, Maps, Mesh
from SimPEG.EM.Base import BaseEMProblem
from SimPEG.EM.Static.DC.FieldsDC import Fields, Fields_CC, Fields_N
@@ -5,7 +11,7 @@ from SimPEG.Utils import sdiag
import numpy as np
from SimPEG.Utils import Zero
from SimPEG.EM.Static.DC import getxBCyBC_CC
from SurveyIP import Survey
from .SurveyIP import Survey
class IPPropMap(Maps.PropMap):
"""
+6
View File
@@ -1,3 +1,9 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
import SimPEG
from SimPEG.EM.Base import BaseEMSurvey
from SimPEG import sp, Survey
+8 -2
View File
@@ -1,2 +1,8 @@
from ProblemIP import Problem3D_CC, Problem3D_N
from SurveyIP import Survey
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from .ProblemIP import Problem3D_CC, Problem3D_N
from .SurveyIP import Survey
+10 -2
View File
@@ -1,3 +1,11 @@
from __future__ import absolute_import
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from builtins import int
from future import standard_library
standard_library.install_aliases()
from builtins import range
from SimPEG import Problem, Utils, Maps, Mesh
from SimPEG.EM.Base import BaseEMProblem
from SimPEG.EM.Static.DC.FieldsDC import Fields, Fields_CC, Fields_N
@@ -5,7 +13,7 @@ from SimPEG.Utils import sdiag
import numpy as np
from SimPEG.Utils import Zero
from SimPEG.EM.Static.DC import getxBCyBC_CC
from SurveySIP import Survey, Data
from .SurveySIP import Survey, Data
class ColeColePropMap(Maps.PropMap):
"""
@@ -105,7 +113,7 @@ class BaseSIPProblem(BaseEMProblem):
JvAll = []
#Assume only eta and tau (eta first then tau)
# v = [2*Mx1]
v = v.reshape((int(v.size/2), 2), order='F')
v = v.reshape((v.size//2), 2), order='F')
for tind in range(len(self.survey.times)):
t = self.survey.times[tind]
+7
View File
@@ -1,3 +1,10 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
from SimPEG import Utils, Maps, Mesh, sp, np
from SimPEG.Regularization import BaseRegularization, Simple
+6
View File
@@ -1,3 +1,9 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
import SimPEG
import numpy as np
from SimPEG.Utils import Zero, closestPoints
+6
View File
@@ -1,3 +1,9 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
import SimPEG
# from SimPEG.EM.Base import BaseEMSurvey
from SimPEG.Utils import Zero, closestPoints, mkvc
+7
View File
@@ -1,3 +1,10 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import str
import SimPEG
from SimPEG.EM.Base import BaseEMSurvey
from SimPEG import np, sp, Survey, Utils
+11 -5
View File
@@ -1,5 +1,11 @@
from ProblemSIP import Problem3D_CC, Problem3D_N
from SurveySIP import Survey, Data
import SrcSIP as Src #Pole
import RxSIP as Rx
from Regularization import MultiRegularization
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from .ProblemSIP import Problem3D_CC, Problem3D_N
from .SurveySIP import Survey, Data
from . import SrcSIP as Src #Pole
from . import RxSIP as Rx
from .Regularization import MultiRegularization
+27 -19
View File
@@ -1,3 +1,11 @@
from __future__ import print_function
from __future__ import division
from __future__ import unicode_literals
from __future__ import absolute_import
from builtins import int
from future import standard_library
standard_library.install_aliases()
from builtins import range
from SimPEG import np
from SimPEG.EM.Static import DC, IP
@@ -85,10 +93,10 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None):
elif stype == 'dpdp':
leg = data * 2*np.pi / ( 1/MA - 1/MB + 1/NB - 1/NA )
LEG.append(1./(2*np.pi) *( 1/MA - 1/MB + 1/NB - 1/NA ))
leg = data * 2*np.pi / (1/MA - 1/MB + 1/NB - 1/NA)
LEG.append(1./(2*np.pi) * (1/MA - 1/MB + 1/NB - 1/NA))
else:
print """dtype must be 'pdp'(pole-dipole) | 'dpdp' (dipole-dipole) """
print("""dtype must be 'pdp'(pole-dipole) | 'dpdp' (dipole-dipole) """)
break
@@ -103,7 +111,7 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None):
rho = np.hstack([rho,leg])
else:
print """dtype must be 'appr' | 'appc' | 'volt' """
print("""dtype must be 'appr' | 'appc' | 'volt' """)
break
@@ -184,14 +192,14 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
# Mesure survey length and direction
dl_len = xy_2_r(endl[0,0],endl[1,0],endl[0,1],endl[1,1])
dl_x = ( endl[1,0] - endl[0,0] ) / dl_len
dl_y = ( endl[1,1] - endl[0,1] ) / dl_len
dl_x = (endl[1,0] - endl[0,0]) / dl_len
dl_y = (endl[1,1] - endl[0,1]) / dl_len
nstn = np.floor( dl_len / a )
nstn = np.floor(dl_len / a)
# Compute discrete pole location along line
stn_x = endl[0,0] + np.array(range(int(nstn)))*dl_x*a
stn_y = endl[0,1] + np.array(range(int(nstn)))*dl_y*a
stn_x = endl[0,0] + np.array(list(range(int(nstn))))*dl_x*a
stn_y = endl[0,1] + np.array(list(range(int(nstn))))*dl_y*a
if mesh.dim==2:
ztop = mesh.vectorNy[-1]
@@ -230,15 +238,15 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
AB = xy_2_r(tx[0,1],endl[1,0],tx[1,1],endl[1,1])
# Number of receivers to fit
nstn = np.min([np.floor( (AB - b) / a ) , n])
nstn = np.min([(AB - b) // a, n])
# Check if there is enough space, else break the loop
if nstn <= 0:
continue
# Compute discrete pole location along line
stn_x = N[ii,0] + dl_x*b + np.array(range(int(nstn)))*dl_x*a
stn_y = N[ii,1] + dl_y*b + np.array(range(int(nstn)))*dl_y*a
stn_x = N[ii,0] + dl_x*b + np.array(list(range(int(nstn))))*dl_x*a
stn_y = N[ii,1] + dl_y*b + np.array(list(range(int(nstn))))*dl_y*a
# Create receiver poles
@@ -275,17 +283,17 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
max_y = endl[1,1] - dl_y * b
box_l = np.sqrt( (min_x - max_x)**2 + (min_y - max_y)**2 )
box_w = box_l/2.
box_w = box_l / 2.
nstn = np.floor( box_l / a )
nstn = np.floor(box_l / a)
# Compute discrete pole location along line
stn_x = min_x + np.array(range(int(nstn)))*dl_x*a
stn_y = min_y + np.array(range(int(nstn)))*dl_y*a
stn_x = min_x + np.array(list(range(int(nstn))))*dl_x*a
stn_y = min_y + np.array(list(range(int(nstn))))*dl_y*a
# Define number of cross lines
nlin = int(np.floor( box_w / a ))
lind = range(-nlin,nlin+1)
nlin = int(box_w // a)
lind = list(range(-nlin,nlin+1))
ngrad = nstn * len(lind)
@@ -310,7 +318,7 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
srcClass = DC.Src.Dipole([rxClass], M[0,:], N[-1,:])
SrcList.append(srcClass)
else:
print """stype must be either 'pdp', 'dpdp' or 'gradient'. """
print("""stype must be either 'pdp', 'dpdp' or 'gradient'. """)
return SrcList
+7 -1
View File
@@ -1 +1,7 @@
from StaticUtils import *
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from .StaticUtils import *
+9 -3
View File
@@ -1,3 +1,9 @@
import DC
import IP
import SIP
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from . import DC
from . import IP
from . import SIP
+24 -17
View File
@@ -1,3 +1,10 @@
from __future__ import print_function
from __future__ import unicode_literals
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
from SimPEG import Solver, Problem
from SimPEG.Problem import BaseTimeProblem
from SimPEG.EM.Utils import *
@@ -47,7 +54,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
self.waveformType = "GENERAL"
def fields(self, m):
if self.verbose: print '%s\nCalculating fields(m)\n%s'%('*'*50,'*'*50)
if self.verbose: print('%s\nCalculating fields(m)\n%s'%('*'*50,'*'*50))
self.curModel = m
# Create a fields storage object
F = self._FieldsForward_pair(self.mesh, self.survey)
@@ -55,7 +62,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
# Set the initial conditions
F[src,:,0] = src.getInitialFields(self.mesh)
F = self.forward(m, self.getRHS, F=F)
if self.verbose: print '%s\nDone calculating fields(m)\n%s'%('*'*50,'*'*50)
if self.verbose: print('%s\nDone calculating fields(m)\n%s'%('*'*50,'*'*50))
return F
def forward(self, m, RHS, F=None):
@@ -70,13 +77,13 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
if Ainv is not None:
Ainv.clean()
A = self.getA(tInd)
if self.verbose: print 'Factoring... (dt = %e)'%dt
if self.verbose: print('Factoring... (dt = %e)'%dt)
Ainv = self.Solver(A, **self.solverOpts)
if self.verbose: print 'Done'
if self.verbose: print('Done')
rhs = RHS(tInd, F)
if self.verbose: print ' Solving... (tInd = %d)'%tInd
if self.verbose: print(' Solving... (tInd = %d)'%tInd)
sol = Ainv * rhs
if self.verbose: print ' Done...'
if self.verbose: print(' Done...')
if sol.ndim == 1:
sol.shape = (sol.size,1)
F[:,self.solType,tInd+1] = sol
@@ -95,13 +102,13 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
if Ainv is not None:
Ainv.clean()
A = self.getA(tInd)
if self.verbose: print 'Factoring (Adjoint)... (dt = %e)'%dt
if self.verbose: print('Factoring (Adjoint)... (dt = %e)'%dt)
Ainv = self.Solver(A, **self.solverOpts)
if self.verbose: print 'Done'
if self.verbose: print('Done')
rhs = RHS(tInd, F)
if self.verbose: print ' Solving (Adjoint)... (tInd = %d)'%tInd
if self.verbose: print(' Solving (Adjoint)... (tInd = %d)'%tInd)
sol = Ainv * rhs
if self.verbose: print ' Done...'
if self.verbose: print(' Done...')
if sol.ndim == 1:
sol.shape = (sol.size,1)
F[:,self.solType,tInd+1] = sol
@@ -112,7 +119,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
"""
:param numpy.array m: Conductivity model
:param numpy.ndarray v: vector (model object)
:param simpegEM.TDEM.FieldsTDEM f: Fields resulting from m
:param FieldsTDEM f: Fields resulting from m
:rtype: numpy.ndarray
:return: w (data object)
@@ -123,21 +130,21 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
* Compute \\\(\\\\vec{w} = -\\\mathbf{Q} \\\\vec{y}\\\)
"""
if self.verbose: print '%s\nCalculating J(v)\n%s'%('*'*50,'*'*50)
if self.verbose: print('%s\nCalculating J(v)\n%s'%('*'*50,'*'*50))
self.curModel = m
if f is None:
f = self.fields(m)
p = self.Gvec(m, v, f)
y = self.solveAh(m, p)
Jv = self.survey.evalDeriv(f, v=y)
if self.verbose: print '%s\nDone calculating J(v)\n%s'%('*'*50,'*'*50)
if self.verbose: print('%s\nDone calculating J(v)\n%s'%('*'*50,'*'*50))
return - mkvc(Jv)
def Jtvec(self, m, v, f=None):
"""
:param numpy.array m: Conductivity model
:param numpy.ndarray,SimPEG.Survey.Data v: vector (data object)
:param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m
:param numpy.ndarray v: vector (or a :class:`SimPEG.Survey.Data` object)
:param FieldsTDEM u: Fields resulting from m
:rtype: numpy.ndarray
:return: w (model object)
@@ -148,7 +155,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
* Compute \\\(\\\\vec{w} = -\\\mathbf{G}^\\\\top y\\\)
"""
if self.verbose: print '%s\nCalculating J^T(v)\n%s'%('*'*50,'*'*50)
if self.verbose: print('%s\nCalculating J^T(v)\n%s'%('*'*50,'*'*50))
self.curModel = m
if f is None:
f = self.fields(m)
@@ -159,6 +166,6 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
p = self.survey.evalDeriv(f, v=v, adjoint=True)
y = self.solveAht(m, p)
w = self.Gtvec(m, y, f)
if self.verbose: print '%s\nDone calculating J^T(v)\n%s'%('*'*50,'*'*50)
if self.verbose: print('%s\nDone calculating J^T(v)\n%s'%('*'*50,'*'*50))
return - mkvc(w)
+11 -5
View File
@@ -1,7 +1,13 @@
from __future__ import print_function
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from SimPEG import Utils, Survey, np
from SimPEG.Survey import BaseSurvey
from SimPEG.EM.Utils import *
from BaseTDEM import FieldsTDEM
from .BaseTDEM import FieldsTDEM
class RxTDEM(Survey.BaseTimeRx):
@@ -87,7 +93,7 @@ class SrcTDEM_VMD_MVP(SrcTDEM):
def getInitialFields(self, mesh):
"""Vertical magnetic dipole, magnetic vector potential"""
if self.waveformType == "STEPOFF":
print ">> Step waveform: Non-zero initial condition"
print(">> Step waveform: Non-zero initial condition")
if mesh._meshType is 'CYL':
if mesh.isSymmetric:
MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
@@ -99,7 +105,7 @@ class SrcTDEM_VMD_MVP(SrcTDEM):
raise Exception('Unknown mesh for VMD')
return {"b": mesh.edgeCurl*MVP}
elif self.waveformType == "GENERAL":
print ">> General waveform: Zero initial condition"
print(">> General waveform: Zero initial condition")
return {"b": np.zeros(mesh.nF)}
else:
raise NotImplementedError("Only use STEPOFF or GENERAL")
@@ -127,7 +133,7 @@ class SrcTDEM_CircularLoop_MVP(SrcTDEM):
def getInitialFields(self, mesh):
"""Circular Loop, magnetic vector potential"""
if self.waveformType == "STEPOFF":
print ">> Step waveform: Non-zero initial condition"
print(">> Step waveform: Non-zero initial condition")
if mesh._meshType is 'CYL':
if mesh.isSymmetric:
MVP = MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
@@ -139,7 +145,7 @@ class SrcTDEM_CircularLoop_MVP(SrcTDEM):
raise Exception('Unknown mesh for CircularLoop')
return {"b": mesh.edgeCurl*MVP}
elif self.waveformType == "GENERAL":
print ">> General waveform: Zero initial condition"
print(">> General waveform: Zero initial condition")
return {"b": np.zeros(mesh.nF)}
else:
raise NotImplementedError("Only use STEPOFF or GENERAL")
+22 -15
View File
@@ -1,7 +1,14 @@
from BaseTDEM import BaseTDEMProblem, FieldsTDEM
from __future__ import absolute_import
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from future import standard_library
standard_library.install_aliases()
from builtins import range
from .BaseTDEM import BaseTDEMProblem, FieldsTDEM
from SimPEG.Utils import mkvc, sdiag
import numpy as np
from SurveyTDEM import SurveyTDEM
from .SurveyTDEM import SurveyTDEM
class FieldsTDEM_e_from_b(FieldsTDEM):
@@ -87,8 +94,8 @@ class ProblemTDEM_b(BaseTDEMProblem):
"""
:param numpy.array m: Conductivity model
:param numpy.array vec: vector (like a model)
:param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m
:rtype: simpegEM.TDEM.FieldsTDEM
:param FieldsTDEM u: Fields resulting from m
:rtype: FieldsTDEM
:return: f
Multiply G by a vector
@@ -125,9 +132,9 @@ class ProblemTDEM_b(BaseTDEMProblem):
"""
:param numpy.array m: Conductivity model
:param numpy.array vec: vector (like a fields)
:param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m
:rtype: np.ndarray (like a model)
:return: p
:param FieldsTDEM u: Fields resulting from m
:rtype: numpy.ndarray
:return: p (like a model)
Multiply G.T by a vector
"""
@@ -153,8 +160,8 @@ class ProblemTDEM_b(BaseTDEMProblem):
def solveAh(self, m, p):
"""
:param numpy.array m: Conductivity model
:param simpegEM.TDEM.FieldsTDEM p: Fields object
:rtype: simpegEM.TDEM.FieldsTDEM
:param FieldsTDEM p: Fields object
:rtype: FieldsTDEM
:return: y
Solve the block-matrix system \\\(\\\hat{A} \\\hat{y} = \\\hat{p}\\\):
@@ -200,8 +207,8 @@ class ProblemTDEM_b(BaseTDEMProblem):
def solveAht(self, m, p):
"""
:param numpy.array m: Conductivity model
:param simpegEM.TDEM.FieldsTDEM p: Fields object
:rtype: simpegEM.TDEM.FieldsTDEM
:param FieldsTDEM p: Fields object
:rtype: FieldsTDEM
:return: y
Solve the block-matrix system \\\(\\\hat{A}^\\\\top \\\hat{y} = \\\hat{p}\\\):
@@ -270,8 +277,8 @@ class ProblemTDEM_b(BaseTDEMProblem):
def _AhVec(self, m, vec):
"""
:param numpy.array m: Conductivity model
:param simpegEM.TDEM.FieldsTDEM vec: Fields object
:rtype: simpegEM.TDEM.FieldsTDEM
:param FieldsTDEM vec: Fields object
:rtype: FieldsTDEM
:return: f
Multiply the matrix \\\(\\\hat{A}\\\) by a fields vector where
@@ -315,8 +322,8 @@ class ProblemTDEM_b(BaseTDEMProblem):
def _AhtVec(self, m, vec):
"""
:param numpy.array m: Conductivity model
:param simpegEM.TDEM.FieldsTDEM vec: Fields object
:rtype: simpegEM.TDEM.FieldsTDEM
:param FieldsTDEM vec: Fields object
:rtype: FieldsTDEM
:return: f
Multiply the matrix \\\(\\\hat{A}\\\) by a fields vector where
+9 -3
View File
@@ -1,3 +1,9 @@
from SurveyTDEM import * #SurveyTDEM, RxTDEM, SrcTDEM
from BaseTDEM import BaseTDEMProblem, FieldsTDEM
from TDEM_b import ProblemTDEM_b
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from .SurveyTDEM import * #SurveyTDEM, RxTDEM, SrcTDEM
from .BaseTDEM import BaseTDEMProblem, FieldsTDEM
from .TDEM_b import ProblemTDEM_b
+12 -5
View File
@@ -1,3 +1,10 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
from SimPEG import *
from scipy.special import ellipk, ellipe
from scipy.constants import mu_0, pi
@@ -17,7 +24,7 @@ def MagneticDipoleVectorPotential(srcLoc, obsLoc, component, moment=1., dipoleMo
#TODO: break this out!
if type(component) in [list, tuple]:
out = range(len(component))
out = list(range(len(component)))
for i, comp in enumerate(component):
out[i] = MagneticDipoleVectorPotential(srcLoc, obsLoc, comp, dipoleMoment=dipoleMoment)
return np.concatenate(out)
@@ -118,7 +125,7 @@ def MagneticLoopVectorPotential(srcLoc, obsLoc, component, radius, mu=mu_0):
"""
if type(component) in [list, tuple]:
out = range(len(component))
out = list(range(len(component)))
for i, comp in enumerate(component):
out[i] = MagneticLoopVectorPotential(srcLoc, obsLoc, comp, radius, mu)
return np.concatenate(out)
@@ -158,11 +165,11 @@ def MagneticLoopVectorPotential(srcLoc, obsLoc, component, radius, mu=mu_0):
# % 1/r singular at r = 0 and K(m) singular at m = 1
Aphi = np.zeros(n)
# % Common factor is (mu * I) / pi with I = 1 and mu = 4e-7 * pi.
Aphi[ind] = 4e-7 / np.sqrt(m[ind]) * np.sqrt(radius / r[ind]) *((1. - m[ind] / 2.) * K[ind] - E[ind])
Aphi[ind] = 4e-7 / np.sqrt(m[ind]) * np.sqrt(radius/ r[ind]) *((1. - m[ind] / 2.) * K[ind] - E[ind])
if component == 'x':
A[ind, i] = Aphi[ind] * (-y[ind] / r[ind] )
A[ind, i] = Aphi[ind] * (-y[ind] / r[ind])
elif component == 'y':
A[ind, i] = Aphi[ind] * ( x[ind] / r[ind] )
A[ind, i] = Aphi[ind] * (x[ind] / r[ind])
else:
raise ValueError('Invalid component')
+8 -2
View File
@@ -1,3 +1,9 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
import numpy as np
from scipy.constants import mu_0, epsilon_0
@@ -9,8 +15,8 @@ def omega(freq):
def k(freq, sigma, mu=mu_0, eps=epsilon_0):
""" Eq 1.47 - 1.49 in Ward and Hohmann """
w = omega(freq)
alp = w * np.sqrt( mu*eps/2 * ( np.sqrt(1. + (sigma / (eps*w))**2 ) + 1) )
beta = w * np.sqrt( mu*eps/2 * ( np.sqrt(1. + (sigma / (eps*w))**2 ) - 1) )
alp = w * np.sqrt( mu*eps/2 * ( np.sqrt(1. + (sigma / (eps*w)))**2 ) + 1)
beta = w * np.sqrt( mu*eps/2 * ( np.sqrt(1. + (sigma / (eps*w)))**2 ) - 1)
return alp - 1j*beta
+8 -2
View File
@@ -1,2 +1,8 @@
from EMUtils import omega, k
from AnalyticUtils import MagneticDipoleFields, MagneticDipoleVectorPotential, MagneticLoopVectorPotential
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from .EMUtils import omega, k
from .AnalyticUtils import MagneticDipoleFields, MagneticDipoleVectorPotential, MagneticLoopVectorPotential
+13 -6
View File
@@ -1,3 +1,10 @@
from __future__ import print_function
from __future__ import division
from __future__ import unicode_literals
from __future__ import absolute_import
from builtins import int
from future import standard_library
standard_library.install_aliases()
import unittest
from SimPEG import *
from SimPEG import EM
@@ -58,7 +65,7 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False):
Src.append(EM.FDEM.Src.RawVec([rx0], freq, mesh.getEdgeInnerProduct()*S_m, S_e))
if verbose:
print ' Fetching %s problem' % (fdemType)
print(' Fetching %s problem' % (fdemType))
if fdemType == 'e':
survey = EM.FDEM.Survey(Src)
@@ -83,7 +90,7 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False):
try:
from pymatsolver import MumpsSolver
prb.Solver = MumpsSolver
except ImportError, e:
except ImportError as e:
prb.Solver = SolverLU
return prb
@@ -94,7 +101,7 @@ def crossCheckTest(SrcList, fdemType1, fdemType2, comp, addrandoms = False, useM
prb1 = getFDEMProblem(fdemType1, comp, SrcList, freq, useMu, verbose)
mesh = prb1.mesh
print 'Cross Checking Forward: %s, %s formulations - %s' % (fdemType1, fdemType2, comp)
print('Cross Checking Forward: %s, %s formulations - %s' % (fdemType1, fdemType2, comp))
logsig = np.log(np.ones(mesh.nC)*CONDUCTIVITY)
mu = np.ones(mesh.nC)*MU
@@ -112,7 +119,7 @@ def crossCheckTest(SrcList, fdemType1, fdemType2, comp, addrandoms = False, useM
d1 = survey1.dpred(m)
if verbose:
print ' Problem 1 solved'
print(' Problem 1 solved')
prb2 = getFDEMProblem(fdemType2, comp, SrcList, freq, useMu, verbose)
@@ -121,11 +128,11 @@ def crossCheckTest(SrcList, fdemType1, fdemType2, comp, addrandoms = False, useM
d2 = survey2.dpred(m)
if verbose:
print ' Problem 2 solved'
print(' Problem 2 solved')
r = d2-d1
l2r = l2norm(r)
tol = np.max([TOL*(10**int(np.log10(0.5* (l2norm(d1) + l2norm(d2)) ))),FLR])
print l2norm(d1), l2norm(d2), l2r , tol, l2r < tol
print(l2norm(d1), l2norm(d2), l2r , tol, l2r < tol)
return l2r < tol
+12 -6
View File
@@ -1,7 +1,13 @@
import TDEM
import FDEM
import Static
import Base
import Analytics
import Utils
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from . import TDEM
from . import FDEM
from . import Static
from . import Base
from . import Analytics
from . import Utils
from scipy.constants import mu_0, epsilon_0
+15 -9
View File
@@ -1,7 +1,13 @@
from __future__ import print_function
from __future__ import division
from __future__ import unicode_literals
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import *
import SimPEG.DCIP as DC
import SimPEG.EM.Static.DC as DC
def run(plotIt=False):
def run(plotIt=True):
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,15 +27,15 @@ def run(plotIt=False):
# ax.plot(xyz_rxP[:,0],xyz_rxP[:,1], 'w.')
# ax.plot(xyz_rxN[:,0],xyz_rxN[:,1], 'r.', ms = 3)
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)
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)
problem.pair(survey)
try:
from pymatsolver import MumpsSolver
problem.Solver = MumpsSolver
except Exception, e:
except Exception as e:
pass
data = survey.dpred(sigma)
@@ -61,8 +67,8 @@ def run(plotIt=False):
ax[0].set_title('Computed')
plt.show()
return np.linalg.norm(data-data_ana)/np.linalg.norm(data_ana)
return np.linalg.norm(data-data_ana) / np.linalg.norm(data_ana)
if __name__ == '__main__':
print run(plotIt=True)
print(run())
+14 -6
View File
@@ -1,3 +1,11 @@
from __future__ import print_function
from __future__ import division
from __future__ import unicode_literals
from __future__ import absolute_import
from builtins import int
from future import standard_library
standard_library.install_aliases()
from builtins import range
from SimPEG import Mesh, Utils, np, sp
import SimPEG.DCIP as DC
import time
@@ -57,7 +65,7 @@ def run(loc=None, sig=None, radi=None, param=None, surveyType='dipole-dipole', u
model[ind] = sig[2]
# Get index of the center
indy = int(mesh.nCy/2)
indy = int(mesh.nCy // 2)
# Plot the model for reference
# Define core mesh extent
@@ -124,10 +132,10 @@ def run(loc=None, sig=None, radi=None, param=None, surveyType='dipole-dipole', u
tx = np.squeeze(Tx[ii][:,0:1])
tinf = tx + np.array([dl_x,dl_y,0])*dl_len*2
inds = Utils.closestPoints(mesh, np.c_[tx,tinf].T)
RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*( [-1] / mesh.vol[inds] )
RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*([-1] / mesh.vol[inds])
else:
inds = Utils.closestPoints(mesh, np.asarray(Tx[ii]).T )
RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*( [-1,1] / mesh.vol[inds] )
RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*([-1,1] / mesh.vol[inds])
# Iterative Solve
Ainvb = sp.linalg.bicgstab(P*A,P*RHS, tol=1e-5)
@@ -143,10 +151,10 @@ def run(loc=None, sig=None, radi=None, param=None, surveyType='dipole-dipole', u
dtemp = (P1*phi - P2*phi)*np.pi
data.append( dtemp )
print '\rTransmitter {0} of {1} -> Time:{2} sec'.format(ii,len(Tx),time.time()- start_time),
print('\rTransmitter {0} of {1} -> Time:{2} sec'.format(ii,len(Tx),time.time()- start_time), end=' ')
print 'Transmitter {0} of {1}'.format(ii,len(Tx))
print 'Forward completed'
print('Transmitter {0} of {1}'.format(ii,len(Tx)))
print('Forward completed')
# Let's just convert the 3D format into 2D (distance along line) and plot
survey2D = DC.convertObs_DC3D_to_2D(survey, np.ones(survey.nSrc) , 'Xloc')
+7 -1
View File
@@ -1,3 +1,9 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import *
import SimPEG.EM as EM
from SimPEG.EM import mu_0
@@ -56,7 +62,7 @@ def run(plotIt=True):
try:
from pymatsolver import MumpsSolver
prb.Solver = MumpsSolver
except ImportError, e:
except ImportError as e:
prb.Solver = SolverLU
prb.pair(survey)
@@ -1,3 +1,9 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import *
import SimPEG.EM as EM
+26 -16
View File
@@ -1,3 +1,10 @@
from __future__ import print_function
from __future__ import division
from __future__ import unicode_literals
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
from SimPEG import *
from SimPEG.EM import FDEM, Analytics, mu_0
import time
@@ -19,10 +26,13 @@ def run(plotIt=True):
Morrison Casing Model, and the results are used in a 2016 SEG abstract by
Yang et al.
- 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.
.. 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.
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
@@ -64,8 +74,8 @@ def run(plotIt=True):
casing_l = 300 # length of the casing
casing_r = 0.1
casing_a = casing_r - casing_t/2. # inner radius
casing_b = casing_r + casing_t/2. # outer radius
casing_a = casing_r - casing_t / 2. # inner radius
casing_b = casing_r + casing_t / 2. # outer radius
casing_z = np.r_[-casing_l,0.]
@@ -75,25 +85,25 @@ def run(plotIt=True):
src_loc = np.r_[0.,0.,dsz]
inf_loc = np.r_[0.,0.,1e4]
print 'Skin Depth: ', [(500./np.sqrt(sigmaback*_)) for _ in freqs]
print('Skin Depth: ', [(500. / np.sqrt(sigmaback*_)) for _ in freqs])
# ------------------ MESH ------------------
# fine cells near well bore
csx1, csx2 = 2e-3, 60.
pfx1, pfx2 = 1.3, 1.3
ncx1 = np.ceil(casing_b/csx1+2)
ncx1 = np.ceil(casing_b/csx1)+2
# pad nicely to second cell size
npadx1 = np.floor(np.log(csx2/csx1) / np.log(pfx1))
npadx1 = np.log(csx2/csx1) // np.log(pfx1)
hx1a,hx1b = Utils.meshTensor([(csx1,ncx1)]),Utils.meshTensor([(csx1,npadx1,pfx1)])
dx1 = sum(hx1a)+sum(hx1b)
dx1 = np.floor(dx1/csx2)
hx1b *= (dx1*csx2 - sum(hx1a))/sum(hx1b)
dx1 = dx1 // csx2
hx1b *= (dx1*csx2 - sum(hx1a)) / sum(hx1b)
# second chunk of mesh
dx2 = 300. # uniform mesh out to here
ncx2 = np.ceil((dx2 - dx1)/csx2)
ncx2 = np.ceil((dx2 - dx1) / csx2)
npadx2 = 45
hx2a, hx2b = Utils.meshTensor([(csx2,ncx2)]), Utils.meshTensor([(csx2,npadx2,pfx2)])
hx = np.hstack([hx1a,hx1b,hx2a,hx2b])
@@ -107,8 +117,8 @@ def run(plotIt=True):
# Mesh
mesh = Mesh.CylMesh([hx,1.,hz], [0.,0.,-np.sum(hz[:npadzu+ncz-nza])])
print 'Mesh Extent xmax: %f,: zmin: %f, zmax: %f'%(mesh.vectorCCx.max(), mesh.vectorCCz.min(), mesh.vectorCCz.max())
print 'Number of cells', mesh.nC
print('Mesh Extent xmax: %f,: zmin: %f, zmax: %f'%(mesh.vectorCCx.max(), mesh.vectorCCz.min(), mesh.vectorCCz.max()))
print('Number of cells', mesh.nC)
if plotIt is True:
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
@@ -215,13 +225,13 @@ 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)
problem = FDEM.Problem3D_h(mesh, mapping=mapping, Solver=solver)
problem.pair(survey)
# ------------- Solve ---------------------------
t0 = time.time()
fieldsCasing = problem.fields(sigCasing)
print 'Time to solve 2 sources', time.time() - t0
print('Time to solve 2 sources', time.time() - t0)
# Plot current
@@ -248,9 +258,9 @@ def run(plotIt=True):
in1_in = in1[np.r_[inds]]
z_in = mesh.gridFz[inds_fz,2]
in0_in = in0_in.reshape([in0_in.shape[0]/3,3])
in1_in = in1_in.reshape([in1_in.shape[0]/3,3])
z_in = z_in.reshape([z_in.shape[0]/3,3])
in0_in = in0_in.reshape([in0_in.shape[0]//3,3])
in1_in = in1_in.reshape([in1_in.shape[0]//3,3])
z_in = z_in.reshape([z_in.shape[0]//3,3])
I0 = in0_in.sum(1).real
I1 = in1_in.sum(1).real
+6
View File
@@ -1,3 +1,9 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import *
import SimPEG.EM as EM
from SimPEG.EM import mu_0
@@ -1,3 +1,9 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import *
from SimPEG.FLOW import Richards
+13 -5
View File
@@ -1,3 +1,11 @@
from __future__ import print_function
from __future__ import division
from __future__ import unicode_literals
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import str
from builtins import range
from SimPEG import *
@@ -47,7 +55,7 @@ def run(N=100, plotIt=True):
# Distance weighting
wr = np.sum(prob.G**2.,axis=0)**0.5
wr = ( wr/np.max(wr) )
wr = ( wr/np.max(wr))
dmis = DataMisfit.l2_DataMisfit(survey)
dmis.Wd = 1./wd
@@ -59,15 +67,15 @@ def run(N=100, plotIt=True):
reg.cell_weights = wr
reg.mref = np.zeros(mesh.nC)
opt = Optimization.ProjectedGNCG(maxIter=100 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 10, tolCG = 1e-3)
invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
update_Jacobi = Directives.Update_lin_PreCond()
# Set the IRLS directive, penalize the lowest 25 percentile of model values
# Start with an l2-l2, then switch to lp-norms
norms = [0., 0., 2., 2.]
norms = [0., 0., 2., 2.]
IRLS = Directives.Update_IRLS( norms=norms, prctile = 25, maxIRLSiter = 15, minGNiter=3)
inv = Inversion.BaseInversion(invProb, directiveList=[IRLS,betaest,update_Jacobi])
@@ -75,7 +83,7 @@ def run(N=100, plotIt=True):
# Run inversion
mrec = inv.run(m0)
print "Final misfit:" + str(invProb.dmisfit.eval(mrec))
print("Final misfit:" + str(invProb.dmisfit.eval(mrec)))
if plotIt:
+7
View File
@@ -1,3 +1,10 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
from SimPEG import *
+9 -3
View File
@@ -1,3 +1,9 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
import SimPEG as simpeg
import numpy as np
import SimPEG.MT as MT
@@ -7,7 +13,7 @@ import matplotlib.pyplot as plt
def run(plotIt=True):
"""
MT: 1D: Inversion
=======================
=================
Forward model 1D MT data.
Setup and run a MT 1D inversion.
@@ -50,7 +56,7 @@ def run(plotIt=True):
m_0 = np.log(sigma_0[active])
# Set the mapping
actMap = simpeg.Maps.ActiveCells(m1d, active, np.log(1e-8), nC=m1d.nCx)
actMap = simpeg.Maps.InjectActiveCells(m1d, active, np.log(1e-8), nC=m1d.nCx)
mappingExpAct = simpeg.Maps.ExpMap(m1d) * actMap
## Setup the layout of the survey, set the sources and the connected receivers
@@ -76,7 +82,7 @@ def run(plotIt=True):
survey.dobs = survey.dtrue + 0.025*abs(survey.dtrue)*np.random.randn(*survey.dtrue.shape)
if plotIt:
fig = MT.Utils.dataUtils.plotMT1DModelData(problem)
fig = MT.Utils.dataUtils.plotMT1DModelData(problem, [m_0])
fig.suptitle('Target - smooth true')
+9 -4
View File
@@ -1,3 +1,9 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
# Test script to use SimPEG.MT platform to forward model synthetic data.
# Import
@@ -12,7 +18,7 @@ except:
def run(plotIt=True, nFreq=1):
"""
MT: 3D: Forward
=======================
===============
Forward model 3D MT data.
@@ -46,16 +52,15 @@ def run(plotIt=True, nFreq=1):
survey = MT.Survey(srcList)
## Setup the problem object
problem = MT.Problem3D.eForm_ps(M, sigmaPrimary=sigBG)
problem = MT.Problem3D.eForm_ps(M, sigmaPrimary=sigBG, Solver=Solver)
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 = MT.Data(survey, dataVec)
# Add plots
if plotIt:
pass
+69
View File
@@ -0,0 +1,69 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from builtins import dict
from future import standard_library
standard_library.install_aliases()
from SimPEG import Mesh, Maps, np
def run(plotIt=True):
"""
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()
+47
View File
@@ -0,0 +1,47 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import Mesh, Maps, Utils
def run(plotIt=True):
"""
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()
+6
View File
@@ -1,3 +1,9 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import Mesh, Utils, np, SolverLU
def run(plotIt=True):
+6
View File
@@ -1,3 +1,9 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import *
def run(plotIt=True):
+6
View File
@@ -1,3 +1,9 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import *
def run(plotIt=True):
@@ -1,3 +1,11 @@
from __future__ import print_function
from __future__ import unicode_literals
from __future__ import division
from __future__ import absolute_import
from builtins import int
from future import standard_library
standard_library.install_aliases()
from builtins import zip
from SimPEG import *
def run(plotIt=True, n=60):
@@ -87,7 +95,7 @@ def run(plotIt=True, n=60):
if elapsed > capture[jj]:
PHIS += [(elapsed, phi.copy())]
jj += 1
if ii % 10 == 0: print ii, elapsed
if ii % 10 == 0: print(ii, elapsed)
ii += 1
if plotIt:
@@ -1,3 +1,10 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
from SimPEG import *
def run(plotIt=True):
+11 -3
View File
@@ -1,3 +1,11 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import zip
from builtins import range
from SimPEG import *
def run(plotIt=True, n=60):
@@ -28,15 +36,15 @@ def run(plotIt=True, n=60):
axes[0].set_xlim([-1,17])
axes[0].set_ylim([-1,17])
for ii, loc in zip(range(M.nC),M.gridCC):
for ii, loc in zip(list(range(M.nC)),M.gridCC):
axes[0].text(loc[0]+0.2,loc[1],'%d'%ii, color='r')
axes[0].plot(M.gridFx[:,0],M.gridFx[:,1], 'g>')
for ii, loc in zip(range(M.nFx),M.gridFx):
for ii, loc in zip(list(range(M.nFx)),M.gridFx):
axes[0].text(loc[0]+0.2,loc[1],'%d'%ii, color='g')
axes[0].plot(M.gridFy[:,0],M.gridFy[:,1], 'm^')
for ii, loc in zip(range(M.nFy),M.gridFy):
for ii, loc in zip(list(range(M.nFy)),M.gridFy):
axes[0].text(loc[0]+0.2,loc[1]+0.2,'%d'%(ii+M.nFx), color='m')
axes[1].spy(M.faceDiv)
@@ -1,3 +1,9 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import *
def run(plotIt=True):
+6
View File
@@ -1,3 +1,9 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import *
def run(plotIt=True):
+15 -7
View File
@@ -1,9 +1,19 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import *
from SimPEG.Utils import surface2ind_topo
def run(plotIt=False, nx = 5, ny = 5):
def run(plotIt=True, nx=5, ny=5):
"""
Utils: surface2ind_topo
=======================
Here we show how to use :code:`Utils.surface2ind_topo` to identify cells below
a topographic surface.
@@ -13,27 +23,25 @@ def run(plotIt=False, 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()
+34 -25
View File
@@ -1,28 +1,37 @@
from __future__ import print_function
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import division
from builtins import open
from future import standard_library
standard_library.install_aliases()
# Run this file to add imports.
##### AUTOIMPORTS #####
import DC_Analytic_Dipole
import DC_Forward_PseudoSection
import EM_FDEM_1D_Inversion
import EM_FDEM_Analytic_MagDipoleWholespace
import EM_Schenkel_Morrison_Casing
import EM_TDEM_1D_Inversion
import FLOW_Richards_1D_Celia1990
import Inversion_IRLS
import Inversion_Linear
import Mesh_Basic_ForwardDC
import Mesh_Basic_PlotImage
import Mesh_Basic_Types
import Mesh_Operators_CahnHilliard
import Mesh_QuadTree_Creation
import Mesh_QuadTree_FaceDiv
import Mesh_QuadTree_HangingNodes
import Mesh_Tensor_Creation
import MT_1D_ForwardAndInversion
import MT_3D_Foward
import Utils_surface2ind_topo
from . import DC_Analytic_Dipole
from . import DC_Forward_PseudoSection
from . import EM_FDEM_1D_Inversion
from . import EM_FDEM_Analytic_MagDipoleWholespace
from . import EM_Schenkel_Morrison_Casing
from . import EM_TDEM_1D_Inversion
from . import FLOW_Richards_1D_Celia1990
from . import Inversion_IRLS
from . import Inversion_Linear
from . import Maps_ComboMaps
from . import Maps_Mesh2Mesh
from . import Mesh_Basic_ForwardDC
from . import Mesh_Basic_PlotImage
from . import Mesh_Basic_Types
from . import Mesh_Operators_CahnHilliard
from . import Mesh_QuadTree_Creation
from . import Mesh_QuadTree_FaceDiv
from . import Mesh_QuadTree_HangingNodes
from . import Mesh_Tensor_Creation
from . import MT_1D_ForwardAndInversion
from . import MT_3D_Foward
from . import Utils_surface2ind_topo
__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Inversion_IRLS", "Inversion_Linear", "Mesh_Basic_ForwardDC", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward", "Utils_surface2ind_topo"]
__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Inversion_IRLS", "Inversion_Linear", "Maps_ComboMaps", "Maps_Mesh2Mesh", "Mesh_Basic_ForwardDC", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward", "Utils_surface2ind_topo"]
##### AUTOIMPORTS #####
@@ -38,7 +47,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', 'examples'])
docExamplesDir = os.path.sep.join(fName.split(os.path.sep)[:-3] + ['docs', 'content', 'examples'])
shutil.rmtree(docExamplesDir)
os.makedirs(docExamplesDir)
@@ -95,14 +104,14 @@ 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', 'examples', name + '.rst']))
rst = os.path.sep.join((filePath.split(os.path.sep)[:-3] + ['docs', 'content', 'examples', name + '.rst']))
print 'Creating: %s.rst'%name
print('Creating: %s.rst'%name)
f = open(rst, 'w')
f.write(out)
f.close()
+14 -8
View File
@@ -1,14 +1,20 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import object
from SimPEG import Mesh, Maps, Utils, np
from future.utils import with_metaclass
class NonLinearMap(object):
class NonLinearMap(with_metaclass(Utils.SimPEGMetaClass, object)):
"""
SimPEG NonLinearMap
"""
__metaclass__ = Utils.SimPEGMetaClass
counter = None #: A SimPEG.Utils.Counter object
mesh = None #: A SimPEG Mesh
@@ -31,7 +37,7 @@ class NonLinearMap(object):
"""
:param numpy.array u: fields
:param numpy.array m: model
:rtype: scipy.csr_matrix
:rtype: scipy.sparse.csr_matrix
:return: derivative of transformed model
The *transform* changes the model into the physical property.
@@ -44,7 +50,7 @@ class NonLinearMap(object):
"""
:param numpy.array u: fields
:param numpy.array m: model
:rtype: scipy.csr_matrix
:rtype: scipy.sparse.csr_matrix
:return: derivative of transformed model
The *transform* changes the model into the physical property.
@@ -186,7 +192,7 @@ class _haverkamp_theta(NonLinearMap):
def transformDerivU(self, u, m):
self.setModel(m)
g = (self.alpha*((self.theta_s - self.theta_r)/
g = (self.alpha*((self.theta_s - self.theta_r) /
(self.alpha + abs(u)**self.beta)**2)
*(-self.beta*abs(u)**(self.beta-1)*np.sign(u)))
g[u >= 0] = 0
@@ -273,7 +279,7 @@ class _vangenuchten_theta(NonLinearMap):
def transform(self, u, m):
self.setModel(m)
m = 1 - 1.0/self.n
f = (( self.theta_s - self.theta_r )/
f = (( self.theta_s - self.theta_r ) /
((1+abs(self.alpha*u)**self.n)**m) + self.theta_r)
if Utils.isScalar(self.theta_s):
f[u >= 0] = self.theta_s
@@ -343,7 +349,7 @@ class _vangenuchten_k(NonLinearMap):
Ks = self.Ks
m = 1.0 - 1.0/n
g = I*alpha*n*np.exp(Ks)*abs(alpha*u)**(n - 1.0)*np.sign(alpha*u)*(1.0/n - 1.0)*((abs(alpha*u)**n + 1)**(1.0/n - 1))**(I - 1)*((1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)**2*(abs(alpha*u)**n + 1)**(1.0/n - 2) - (2*alpha*n*np.exp(Ks)*abs(alpha*u)**(n - 1)*np.sign(alpha*u)*(1.0/n - 1)*((abs(alpha*u)**n + 1)**(1.0/n - 1))**I*((1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)*(abs(alpha*u)**n + 1)**(1.0/n - 2))/(((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1) + 1)*(1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1.0/n))
g = I*alpha*n*np.exp(Ks)*abs(alpha*u)**(n - 1.0)*np.sign(alpha*u)*(1.0/n - 1.0)*((abs(alpha*u)**n + 1)**(1.0/n - 1))**(I - 1)*((1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)**2*(abs(alpha*u)**n + 1)**(1.0/n - 2) - (2*alpha*n*np.exp(Ks)*abs(alpha*u)**(n - 1)*np.sign(alpha*u)*(1.0/n - 1)*((abs(alpha*u)**n + 1)**(1.0/n - 1))**I*((1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1 - 1.0/n) - 1)*(abs(alpha*u)**n + 1)**(1.0/n - 2))/(((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1))+ 1)*(1 - 1.0/((abs(alpha*u)**n + 1)**(1.0/n - 1))**(1.0/(1.0/n - 1)))**(1.0/n)
g[u >= 0] = 0
g = Utils.sdiag(g)
return g
+15 -8
View File
@@ -1,5 +1,12 @@
from __future__ import print_function
from __future__ import absolute_import
from __future__ import division
from __future__ import unicode_literals
from future import standard_library
standard_library.install_aliases()
from builtins import range
from SimPEG import *
from Empirical import RichardsMap
from .Empirical import RichardsMap
import time
@@ -61,7 +68,7 @@ class RichardsSurvey(Survey.BaseSurvey):
@Utils.requires('prob')
def eval(self, U, m):
Ds = range(len(self.rxList))
Ds = list(range(len(self.rxList)))
for ii, rx in enumerate(self.rxList):
Ds[ii] = rx.eval(U, m,
self.prob.mapping,
@@ -73,7 +80,7 @@ class RichardsSurvey(Survey.BaseSurvey):
@Utils.requires('prob')
def evalDeriv(self, U, m):
"""The Derivative with respect to the fields."""
Ds = range(len(self.rxList))
Ds = list(range(len(self.rxList)))
for ii, rx in enumerate(self.rxList):
Ds[ii] = rx.evalDeriv(U, m,
self.prob.mapping,
@@ -135,12 +142,12 @@ class RichardsProblem(Problem.BaseTimeProblem):
@Utils.timeIt
def fields(self, m):
tic = time.time()
u = range(self.nT+1)
u = list(range(self.nT+1))
u[0] = self.initialConditions
for ii, dt in enumerate(self.timeSteps):
bc = self.getBoundaryConditions(ii, u[ii])
u[ii+1] = self.rootFinder.root(lambda hn1m, return_g=True: self.getResidual(m, u[ii], hn1m, dt, bc, return_g=return_g), u[ii])
if self.debug: print "Solving Fields (%4d/%d - %3.1f%% Done) %d Iterations, %4.2f seconds"%(ii+1, self.nT, 100.0*(ii+1)/self.nT, self.rootFinder.iter, time.time() - tic)
if self.debug: print("Solving Fields (%4d/%d - %3.1f%% Done) %d Iterations, %4.2f seconds"%(ii+1, self.nT, 100.0*(ii+1)/self.nT, self.rootFinder.iter, time.time() - tic))
return u
@Utils.timeIt
@@ -238,7 +245,7 @@ class RichardsProblem(Problem.BaseTimeProblem):
f = self.fields(m)
nn = len(f)-1
Asubs, Adiags, Bs = range(nn), range(nn), range(nn)
Asubs, Adiags, Bs = list(range(nn)), list(range(nn)), list(range(nn))
for ii in range(nn):
dt = self.timeSteps[ii]
bc = self.getBoundaryConditions(ii, f[ii])
@@ -263,7 +270,7 @@ class RichardsProblem(Problem.BaseTimeProblem):
if f is None:
f = self.fields(m)
JvC = range(len(f)-1) # Cell to hold each row of the long vector.
JvC = list(range(len(f)-1)) # Cell to hold each row of the long vector.
# This is done via forward substitution.
bc = self.getBoundaryConditions(0, f[0])
@@ -295,7 +302,7 @@ class RichardsProblem(Problem.BaseTimeProblem):
bc = self.getBoundaryConditions(ii-1, f[ii-1])
Asub, Adiag, B = self.diagsJacobian(m, f[ii-1], f[ii], self.timeSteps[ii-1], bc)
#select the correct part of v
vpart = range((ii)*Adiag.shape[0], (ii+1)*Adiag.shape[0])
vpart = list(range((ii)*Adiag.shape[0], (ii+1)*Adiag.shape[0]))
AdiaginvT = self.Solver(Adiag.T, **self.solverOpts)
JTvC = AdiaginvT * (PTv[vpart] - minus)
minus = Asub.T*JTvC # this is now the super diagonal.
+8 -2
View File
@@ -1,2 +1,8 @@
import Empirical
from RichardsProblem import *
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from . import Empirical
from .RichardsProblem import *
+7 -1
View File
@@ -1 +1,7 @@
import Richards
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from . import Richards
+11 -2
View File
@@ -1,4 +1,13 @@
import Utils, numpy as np, scipy.sparse as sp
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
from builtins import object
from . import Utils
import numpy as np, scipy.sparse as sp
class Fields(object):
"""Fancy Field Storage
@@ -244,7 +253,7 @@ class TimeFields(Fields):
out = func(pointerFields, srcII, timeII)
else: #loop over the time steps
nT = pointerShape[2]
out = range(nT)
out = list(range(nT))
for i, TIND_i in enumerate(timeII):
fieldI = pointerFields[:,:,i]
if fieldI.shape[0] == fieldI.size:
+19 -12
View File
@@ -1,14 +1,21 @@
import Utils, Survey, Problem, numpy as np, scipy.sparse as sp, gc
from Utils.SolverUtils import *
import DataMisfit
import Regularization
from __future__ import print_function
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from builtins import object
from . import Utils, Survey, Problem
import numpy as np, scipy.sparse as sp, gc
from .Utils.SolverUtils import *
from . import DataMisfit
from . import Regularization
from future.utils import with_metaclass
class BaseInvProblem(object):
class BaseInvProblem(with_metaclass(Utils.SimPEGMetaClass, object)):
"""BaseInvProblem(dmisfit, reg, opt)"""
__metaclass__ = Utils.SimPEGMetaClass
beta = 1.0 #: Trade-off parameter
debug = False #: Print debugging information
@@ -54,10 +61,10 @@ class BaseInvProblem(object):
Called when inversion is first starting.
"""
if self.debug: print 'Calling InvProblem.startup'
if self.debug: print('Calling InvProblem.startup')
if self.reg.mref is None:
print 'SimPEG.InvProblem will set Regularization.mref to m0.'
print('SimPEG.InvProblem will set Regularization.mref to m0.')
self.reg.mref = m0
self.phi_d = np.nan
@@ -65,8 +72,8 @@ class BaseInvProblem(object):
self.curModel = m0
print """SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
***Done using same Solver and solverOpts as the problem***"""
print("""SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
***Done using same Solver and solverOpts as the problem***""")
self.opt.bfgsH0 = self.prob.Solver(self.reg.eval2Deriv(self.curModel), **self.prob.solverOpts)
@property
@@ -87,7 +94,7 @@ class BaseInvProblem(object):
for mtest, u_ofmtest in self.warmstart:
if m is mtest:
f = u_ofmtest
if self.debug: print 'InvProb is Warm Starting!'
if self.debug: print('InvProb is Warm Starting!')
break
if f is None:
+11 -5
View File
@@ -1,18 +1,24 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from builtins import object
import SimPEG
from SimPEG import Utils, sp, np
from Optimization import Remember, IterationPrinters, StoppingCriteria
import Directives
from .Optimization import Remember, IterationPrinters, StoppingCriteria
from . import Directives
from future.utils import with_metaclass
class BaseInversion(object):
class BaseInversion(with_metaclass(Utils.SimPEGMetaClass, object)):
"""
Inversion Class.
"""
__metaclass__ = Utils.SimPEGMetaClass
name = 'BaseInversion'
debug = False #: Print debugging information
+8 -2
View File
@@ -1,7 +1,13 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from SimPEG import SolverLU as SimpegSolver, PropMaps, Utils, mkvc, sp, np
from SimPEG.EM.FDEM.ProblemFDEM import BaseFDEMProblem
from SurveyMT import Survey, Data
from FieldsMT import BaseMTFields
from .SurveyMT import Survey, Data
from .FieldsMT import BaseMTFields
class BaseMTProblem(BaseFDEMProblem):
+12 -6
View File
@@ -1,3 +1,9 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import Survey, Utils, Problem, np, sp, mkvc
from scipy.constants import mu_0
import sys
@@ -63,7 +69,7 @@ class Fields1D_e(BaseMTFields):
C = self.mesh.nodalGrad
b = (C * eSolution)
for i, src in enumerate(srcList):
b[:,i] *= - 1./(1j*omega(src.freq))
b[:,i] *= -1./(1j*omega(src.freq))
# There is no magnetic source in the MT problem
# S_m, _ = src.eval(self.survey.prob)
# if S_m is not None:
@@ -188,7 +194,7 @@ class Fields3D_e(BaseMTFields):
# adjoint: returns a 2*nE long vector with zero's for py
return np.vstack((v,np.zeros_like(v)))
# Not adjoint: return only the px part of the vector
return v[:len(v)/2]
return v[:len(v)//2]
def _e_pyDeriv_u(self, src, v, adjoint = False):
'''
@@ -198,7 +204,7 @@ class Fields3D_e(BaseMTFields):
# adjoint: returns a 2*nE long vector with zero's for px
return np.vstack((np.zeros_like(v),v))
# Not adjoint: return only the px part of the vector
return v[len(v)/2::]
return v[len(v)//2::]
def _e_pxDeriv_m(self, src, v, adjoint = False):
# assuming primary does not depend on the model
@@ -227,7 +233,7 @@ class Fields3D_e(BaseMTFields):
C = self.mesh.edgeCurl
b = (C * e_pxSolution)
for i, src in enumerate(srcList):
b[:,i] *= - 1./(1j*omega(src.freq))
b[:,i] *= -1./(1j*omega(src.freq))
# There is no magnetic source in the MT problem
# S_m, _ = src.eval(self.survey.prob)
# if S_m is not None:
@@ -238,7 +244,7 @@ class Fields3D_e(BaseMTFields):
C = self.mesh.edgeCurl
b = (C * e_pySolution)
for i, src in enumerate(srcList):
b[:,i] *= - 1./(1j*omega(src.freq))
b[:,i] *= -1./(1j*omega(src.freq))
# There is no magnetic source in the MT problem
# S_m, _ = src.eval(self.survey.prob)
# if S_m is not None:
@@ -348,4 +354,4 @@ class Fields3D_e(BaseMTFields):
This function stacks the fields derivatives appropriately
"""
# The fields have no dependance to the model.
return None
return None
+13 -7
View File
@@ -1,3 +1,9 @@
from __future__ import print_function
from __future__ import division
from __future__ import unicode_literals
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG.EM.Utils import omega
from SimPEG import mkvc
from scipy.constants import mu_0
@@ -46,7 +52,7 @@ class eForm_psField(BaseMTProblem):
Edge inner product matrix
"""
if getattr(self, '_MeMui', None) is None:
self._MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
self._MeMui = self.mesh.getEdgeInnerProduct(old_div(1.0,mu_0))
return self._MeMui
@property
@@ -142,7 +148,7 @@ class eForm_psField(BaseMTProblem):
for freq in self.survey.freqs:
if self.verbose:
startTime = time.time()
print 'Starting work for {:.3e}'.format(freq)
print('Starting work for {:.3e}'.format(freq))
sys.stdout.flush()
A = self.getA(freq)
rhs = self.getRHS(freq)
@@ -158,7 +164,7 @@ class eForm_psField(BaseMTProblem):
# b = -( self.mesh.nodalGrad * e )/( 1j*omega(freq) )
# F[Src, 'b_1d'] = b[:,1]
if self.verbose:
print 'Ran for {:f} seconds'.format(time.time()-startTime)
print('Ran for {:f} seconds'.format(time.time()-startTime))
sys.stdout.flush()
return F
@@ -191,7 +197,7 @@ class eForm_TotalField(BaseMTProblem):
Edge inner product matrix
"""
if getattr(self, '_MeMui', None) is None:
self._MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
self._MeMui = self.mesh.getEdgeInnerProduct(old_div(1.0,mu_0))
return self._MeMui
@property
@@ -249,7 +255,7 @@ class eForm_TotalField(BaseMTProblem):
Ed, Eu, Hd, Hu = getEHfields(self.mesh,self.curModel.sigma,freq,self.mesh.vectorNx)
Etot = (Ed + Eu)
sourceAmp = 1.0
Etot = ((Etot/Etot[-1])*sourceAmp) # Scale the fields to be equal to sourceAmp at the top
Etot = ((old_div(Etot,Etot[-1]))*sourceAmp) # Scale the fields to be equal to sourceAmp at the top
## Note: The analytic solution is derived with e^iwt
eBC = np.r_[Etot[0],Etot[-1]]
# The right hand side
@@ -274,7 +280,7 @@ class eForm_TotalField(BaseMTProblem):
for freq in self.survey.freqs:
if self.verbose:
startTime = time.time()
print 'Starting work for {:.3e}'.format(freq)
print('Starting work for {:.3e}'.format(freq))
sys.stdout.flush()
A = self.getA(freq)
rhs, e_o = self.getRHS(freq)
@@ -286,6 +292,6 @@ class eForm_TotalField(BaseMTProblem):
# NOTE: only store e fields
F[Src, 'e_1dSolution'] = e[:,0]
if self.verbose:
print 'Ran for {:f} seconds'.format(time.time()-startTime)
print('Ran for {:f} seconds'.format(time.time()-startTime))
sys.stdout.flush()
return F
+7 -1
View File
@@ -1 +1,7 @@
from Probs import eForm_TotalField, eForm_psField
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from .Probs import eForm_TotalField, eForm_psField
+6
View File
@@ -1 +1,7 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
pass
+8 -2
View File
@@ -1,3 +1,9 @@
from __future__ import print_function
from __future__ import unicode_literals
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import Survey, Problem, Utils, Models, np, sp, mkvc, SolverLU as SimpegSolver
from SimPEG.EM.Utils import omega
from scipy.constants import mu_0
@@ -115,7 +121,7 @@ class eForm_ps(BaseMTProblem):
for freq in self.survey.freqs:
if self.verbose:
startTime = time.time()
print 'Starting work for {:.3e}'.format(freq)
print('Starting work for {:.3e}'.format(freq))
sys.stdout.flush()
A = self.getA(freq)
rhs = self.getRHS(freq)
@@ -131,7 +137,7 @@ class eForm_ps(BaseMTProblem):
# Note curl e = -iwb so b = -curl/iw
if self.verbose:
print 'Ran for {:f} seconds'.format(time.time()-startTime)
print('Ran for {:f} seconds'.format(time.time()-startTime))
sys.stdout.flush()
Ainv.clean()
return F
+7 -1
View File
@@ -1 +1,7 @@
from Probs import eForm_ps
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from .Probs import eForm_ps
+12 -6
View File
@@ -1,10 +1,16 @@
from __future__ import absolute_import
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from future import standard_library
standard_library.install_aliases()
from SimPEG import Utils, Problem, Maps, np, sp, mkvc
from SimPEG.EM.FDEM.SrcFDEM import BaseSrc as FDEMBaseSrc
from SimPEG.EM.Utils import omega
from scipy.constants import mu_0
from numpy.lib import recfunctions as recFunc
from Utils.sourceUtils import homo1DModelSource
from Utils import rec2ndarr
from .Utils.sourceUtils import homo1DModelSource
from .Utils import rec2ndarr
import sys
#################
@@ -78,7 +84,7 @@ class polxy_1Dprimary(BaseMTSrc):
C = problem.mesh.nodalGrad
elif problem.mesh.dim == 3:
C = problem.mesh.edgeCurl
bBG_bp = (- C * self.ePrimary(problem) )*(1/( 1j*omega(self.freq) ))
bBG_bp = (- C * self.ePrimary(problem) )*(1/(1j*omega(self.freq)))
return bBG_bp
def S_e(self,problem):
@@ -86,7 +92,7 @@ class polxy_1Dprimary(BaseMTSrc):
Get the electrical field source
"""
e_p = self.ePrimary(problem)
Map_sigma_p = Maps.Vertical1DMap(problem.mesh)
Map_sigma_p = Maps.SurjectVertical1D(problem.mesh)
sigma_p = Map_sigma_p._transform(self.sigma1d)
# Make mass matrix
# Note: M(sig) - M(sig_p) = M(sig - sig_p)
@@ -155,7 +161,7 @@ class polxy_3Dprimary(BaseMTSrc):
C = problem.mesh.nodalGrad
elif problem.mesh.dim == 3:
C = problem.mesh.edgeCurl
bBG_bp = (- C * self.ePrimary(problem) )*(1/( 1j*omega(self.freq) ))
bBG_bp = (- C * self.ePrimary(problem) )*(1/(1j*omega(self.freq)))
return bBG_bp
def S_e(self,problem):
@@ -163,7 +169,7 @@ class polxy_3Dprimary(BaseMTSrc):
Get the electrical field source
"""
e_p = self.ePrimary(problem)
Map_sigma_p = Maps.Vertical1DMap(problem.mesh)
Map_sigma_p = Maps.SurjectVertical1D(problem.mesh)
sigma_p = Map_sigma_p._transform(self.sigma1d)
# Make mass matrix
# Note: M(sig) - M(sig_p) = M(sig - sig_p)
+27 -21
View File
@@ -1,10 +1,16 @@
from __future__ import absolute_import
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from future import standard_library
standard_library.install_aliases()
from SimPEG import Survey as SimPEGsurvey, Utils, Problem, Maps, np, sp, mkvc
from SimPEG.EM.FDEM.SrcFDEM import BaseSrc as FDEMBaseSrc
from SimPEG.EM.Utils import omega
from scipy.constants import mu_0
from numpy.lib import recfunctions as recFunc
from Utils import rec2ndarr
import SrcMT
from .Utils import rec2ndarr
from . import SrcMT
import sys
#################
@@ -76,7 +82,7 @@ class Rx(SimPEGsurvey.BaseRx):
bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0
# Note: Has a minus sign in front, to comply with quadrant calculations.
# Can be derived from zyx case for the 3D case.
f_part_complex = -ex/bx
f_part_complex = old_div(-ex,bx)
# elif self.projType is 'Z2D':
elif self.projType is 'Z3D':
## NOTE: Assumes that e is on edges and b on the faces. Need to generalize that or use a prop of fields to determine that.
@@ -103,13 +109,13 @@ class Rx(SimPEGsurvey.BaseRx):
hy_py = Pby*f[src,'b_py']/mu_0
# Make the complex data
if 'zxx' in self.rxType:
f_part_complex = ( ex_px*hy_py - ex_py*hy_px)/(hx_px*hy_py - hx_py*hy_px)
f_part_complex = old_div(( ex_px*hy_py - ex_py*hy_px),(hx_px*hy_py - hx_py*hy_px))
elif 'zxy' in self.rxType:
f_part_complex = (-ex_px*hx_py + ex_py*hx_px)/(hx_px*hy_py - hx_py*hy_px)
f_part_complex = old_div((-ex_px*hx_py + ex_py*hx_px),(hx_px*hy_py - hx_py*hy_px))
elif 'zyx' in self.rxType:
f_part_complex = ( ey_px*hy_py - ey_py*hy_px)/(hx_px*hy_py - hx_py*hy_px)
f_part_complex = old_div(( ey_px*hy_py - ey_py*hy_px),(hx_px*hy_py - hx_py*hy_px))
elif 'zyy' in self.rxType:
f_part_complex = (-ey_px*hx_py + ey_py*hx_px)/(hx_px*hy_py - hx_py*hy_px)
f_part_complex = old_div((-ey_px*hx_py + ey_py*hx_px),(hx_px*hy_py - hx_py*hy_px))
elif self.projType is 'T3D':
if self.locs.ndim == 3:
horLoc = self.locs[:,:,0]
@@ -127,9 +133,9 @@ class Rx(SimPEGsurvey.BaseRx):
by_py = Pby*f[src,'b_py']
bz_py = Pbz*f[src,'b_py']
if 'tzx' in self.rxType:
f_part_complex = (- by_px*bz_py + by_py*bz_px)/(bx_px*by_py - bx_py*by_px)
f_part_complex = old_div((- by_px*bz_py + by_py*bz_px),(bx_px*by_py - bx_py*by_px))
if 'tzy' in self.rxType:
f_part_complex = ( bx_px*bz_py - bx_py*bz_px)/(bx_px*by_py - bx_py*by_px)
f_part_complex = old_div(( bx_px*bz_py - bx_py*bz_px),(bx_px*by_py - bx_py*by_px))
else:
NotImplementedError('Projection of {:s} receiver type is not implemented.'.format(self.rxType))
@@ -157,8 +163,8 @@ class Rx(SimPEGsurvey.BaseRx):
Pbx = mesh.getInterpolationMat(self.locs[:,-1],'Ex')
# ex = Pex*mkvc(f[src,'e_1d'],2)
# bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0
dP_de = -mkvc(Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))*(Pex*v),2)
dP_db = mkvc( Utils.sdiag(Pex*mkvc(f[src,'e_1d'],2))*(Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)).T*Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)))*(Pbx*f._bDeriv_u(src,v)/mu_0),2)
dP_de = -mkvc(Utils.sdiag(old_div(1.,(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)))*(Pex*v),2)
dP_db = mkvc( Utils.sdiag(Pex*mkvc(f[src,'e_1d'],2))*(Utils.sdiag(old_div(1.,(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))).T*Utils.sdiag(old_div(1.,(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))))*(Pbx*f._bDeriv_u(src,v)/mu_0),2)
PDeriv_complex = np.sum(np.hstack((dP_de,dP_db)),1)
elif self.projType is 'Z2D':
raise NotImplementedError('Has not been implement for 2D impedance tensor')
@@ -198,7 +204,7 @@ class Rx(SimPEGsurvey.BaseRx):
# Update the input vector
sDiag = lambda t: Utils.sdiag(mkvc(t,2))
# Define the components of the derivative
Hd = sDiag(1./(sDiag(hx_px)*hy_py - sDiag(hx_py)*hy_px))
Hd = sDiag(old_div(1.,(sDiag(hx_px)*hy_py - sDiag(hx_py)*hy_px)))
Hd_uV = sDiag(hy_py)*hx_px_u(v) + sDiag(hx_px)*hy_py_u(v) - sDiag(hx_py)*hy_px_u(v) - sDiag(hy_px)*hx_py_u(v)
# Calculate components
if 'zxx' in self.rxType:
@@ -247,7 +253,7 @@ class Rx(SimPEGsurvey.BaseRx):
# Update the input vector
sDiag = lambda t: Utils.sdiag(mkvc(t,2))
# Define the components of the derivative
Hd = sDiag(1./(sDiag(bx_px)*by_py - sDiag(bx_py)*by_px))
Hd = sDiag(old_div(1.,(sDiag(bx_px)*by_py - sDiag(bx_py)*by_px)))
Hd_uV = sDiag(by_py)*bx_px_u(v) + sDiag(bx_px)*by_py_u(v) - sDiag(bx_py)*by_px_u(v) - sDiag(by_px)*bx_py_u(v)
if 'tzx' in self.rxType:
Tij = sDiag(Hd*( - sDiag(by_px)*bz_py + sDiag(by_py)*bz_px ))
@@ -267,8 +273,8 @@ class Rx(SimPEGsurvey.BaseRx):
Pbx = mesh.getInterpolationMat(self.locs[:,-1],'Ex')
# ex = Pex*mkvc(f[src,'e_1d'],2)
# bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0
dP_deTv = -mkvc(Pex.T*Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)).T*v,2)
db_duv = Pbx.T/mu_0*Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))*(Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))).T*Utils.sdiag(Pex*mkvc(f[src,'e_1d'],2)).T*v
dP_deTv = -mkvc(Pex.T*Utils.sdiag(old_div(1.,(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))).T*v,2)
db_duv = Pbx.T/mu_0*Utils.sdiag(old_div(1.,(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)))*(Utils.sdiag(old_div(1.,(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)))).T*Utils.sdiag(Pex*mkvc(f[src,'e_1d'],2)).T*v
dP_dbTv = mkvc(f._bDeriv_u(src,db_duv,adjoint=True),2)
PDeriv_real = np.sum(np.hstack((dP_deTv,dP_dbTv)),1)
elif self.projType is 'Z2D':
@@ -300,17 +306,17 @@ class Rx(SimPEGsurvey.BaseRx):
aey_px_u = lambda vec: f._e_pxDeriv_u(src,Pey.T*vec,adjoint=True)
aex_py_u = lambda vec: f._e_pyDeriv_u(src,Pex.T*vec,adjoint=True)
aey_py_u = lambda vec: f._e_pyDeriv_u(src,Pey.T*vec,adjoint=True)
ahx_px_u = lambda vec: f._b_pxDeriv_u(src,Pbx.T*vec,adjoint=True)/mu_0
ahy_px_u = lambda vec: f._b_pxDeriv_u(src,Pby.T*vec,adjoint=True)/mu_0
ahx_py_u = lambda vec: f._b_pyDeriv_u(src,Pbx.T*vec,adjoint=True)/mu_0
ahy_py_u = lambda vec: f._b_pyDeriv_u(src,Pby.T*vec,adjoint=True)/mu_0
ahx_px_u = lambda vec: old_div(f._b_pxDeriv_u(src,Pbx.T*vec,adjoint=True),mu_0)
ahy_px_u = lambda vec: old_div(f._b_pxDeriv_u(src,Pby.T*vec,adjoint=True),mu_0)
ahx_py_u = lambda vec: old_div(f._b_pyDeriv_u(src,Pbx.T*vec,adjoint=True),mu_0)
ahy_py_u = lambda vec: old_div(f._b_pyDeriv_u(src,Pby.T*vec,adjoint=True),mu_0)
# Update the input vector
# Define shortcuts
sDiag = lambda t: Utils.sdiag(mkvc(t,2))
sVec = lambda t: Utils.sp.csr_matrix(mkvc(t,2))
# Define the components of the derivative
aHd = sDiag(1./(sDiag(ahx_px)*ahy_py - sDiag(ahx_py)*ahy_px))
aHd = sDiag(old_div(1.,(sDiag(ahx_px)*ahy_py - sDiag(ahx_py)*ahy_px)))
aHd_uV = lambda x: ahx_px_u(sDiag(ahy_py)*x) + ahx_px_u(sDiag(ahy_py)*x) - ahy_px_u(sDiag(ahx_py)*x) - ahx_py_u(sDiag(ahy_px)*x)
# Need to fix this to reflect the adjoint
if 'zxx' in self.rxType:
@@ -362,7 +368,7 @@ class Rx(SimPEGsurvey.BaseRx):
sDiag = lambda t: Utils.sdiag(mkvc(t,2))
sVec = lambda t: Utils.sp.csr_matrix(mkvc(t,2))
# Define the components of the derivative
aHd = sDiag(1./(sDiag(abx_px)*aby_py - sDiag(abx_py)*aby_px))
aHd = sDiag(old_div(1.,(sDiag(abx_px)*aby_py - sDiag(abx_py)*aby_px)))
aHd_uV = lambda x: abx_px_u(sDiag(aby_py)*x) + abx_px_u(sDiag(aby_py)*x) - aby_px_u(sDiag(abx_py)*x) - abx_py_u(sDiag(aby_px)*x)
# Need to fix this to reflect the adjoint
if 'tzx' in self.rxType:
+17 -10
View File
@@ -1,3 +1,10 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import zip
# Analytic solution of EM fields due to a plane wave
import numpy as np, SimPEG as simpeg
@@ -33,8 +40,8 @@ def getEHfields(m1d,sigma,freq,zd,scaleUD=True):
# Loop over all the layers, starting at the bottom layer
for lnr, h in enumerate(m1d.hx): # lnr-number of layer, h-thickness of the layer
# Calculate
yp1 = k[lnr]/(w*mu[lnr]) # Admittance of the layer below the current layer
zp = (w*mu[lnr+1])/k[lnr+1] # Impedance in the current layer
yp1 = old_div(k[lnr],(w*mu[lnr])) # Admittance of the layer below the current layer
zp = old_div((w*mu[lnr+1]),k[lnr+1]) # Impedance in the current layer
# Build the propagation matrix
# Convert fields to down/up going components in layer below current layer
@@ -48,7 +55,7 @@ def getEHfields(m1d,sigma,freq,zd,scaleUD=True):
UDp[:,lnr+1] = elamh.dot(Pjinv.dot(Pj1)).dot(UDp[:,lnr])
if scaleUD:
UDp[:,lnr+1::-1] = UDp[:,lnr+1::-1]/UDp[1,lnr+1]
UDp[:,lnr+1::-1] = old_div(UDp[:,lnr+1::-1],UDp[1,lnr+1])
# Calculate the fields
Ed = np.empty((zd.size,),dtype=complex)
@@ -62,14 +69,14 @@ def getEHfields(m1d,sigma,freq,zd,scaleUD=True):
dind = dup >= zd
Ed[dind] = UDp[1,0]*np.exp(-1j*k[0]*(dup-zd[dind]))
Eu[dind] = UDp[0,0]*np.exp(1j*k[0]*(dup-zd[dind]))
Hd[dind] = (k[0]/(w*mu[0]))*UDp[1,0]*np.exp(-1j*k[0]*(dup-zd[dind]))
Hu[dind] = -(k[0]/(w*mu[0]))*UDp[0,0]*np.exp(1j*k[0]*(dup-zd[dind]))
Hd[dind] = (old_div(k[0],(w*mu[0])))*UDp[1,0]*np.exp(-1j*k[0]*(dup-zd[dind]))
Hu[dind] = -(old_div(k[0],(w*mu[0])))*UDp[0,0]*np.exp(1j*k[0]*(dup-zd[dind]))
for ki,mui,epsi,dlow,dup,Up,Dp in zip(k[1::],mu[1::],eps[1::],m1d.vectorNx[:-1],m1d.vectorNx[1::],UDp[0,1::],UDp[1,1::]):
dind = np.logical_and(dup >= zd, zd > dlow)
Ed[dind] = Dp*np.exp(-1j*ki*(dup-zd[dind]))
Eu[dind] = Up*np.exp(1j*ki*(dup-zd[dind]))
Hd[dind] = (ki/(w*mui))*Dp*np.exp(-1j*ki*(dup-zd[dind]))
Hu[dind] = -(ki/(w*mui))*Up*np.exp(1j*ki*(dup-zd[dind]))
Hd[dind] = (old_div(ki,(w*mui)))*Dp*np.exp(-1j*ki*(dup-zd[dind]))
Hu[dind] = -(old_div(ki,(w*mui)))*Up*np.exp(1j*ki*(dup-zd[dind]))
# Return return the fields
return Ed, Eu, Hd, Hu
@@ -92,15 +99,15 @@ def getImpedance(m1d,sigma,freq):
om = 2*np.pi*fr
Zall = np.empty(len(h)+1,dtype='complex')
# Calculate the impedance for the bottom layer
Zall[0] = (mu_0*om)/np.sqrt(mu_0*eps_0*(om)**2 - 1j*mu_0*sigma[0]*om)
Zall[0] = old_div((mu_0*om),np.sqrt(mu_0*eps_0*(om)**2 - 1j*mu_0*sigma[0]*om))
for nr,hi in enumerate(h):
# Calculate the wave number
# print nr,sigma[nr]
k = np.sqrt(mu_0*eps_0*om**2 - 1j*mu_0*sigma[nr]*om)
Z = (mu_0*om)/k
Z = old_div((mu_0*om),k)
Zall[nr+1] = Z *((Zall[nr] + Z*np.tanh(1j*k*hi))/(Z + Zall[nr]*np.tanh(1j*k*hi)))
Zall[nr+1] = Z *(old_div((Zall[nr] + Z*np.tanh(1j*k*hi)),(Z + Zall[nr]*np.tanh(1j*k*hi))))
#pdb.set_trace()
Z1d[nrFr] = Zall[-1]
+9 -3
View File
@@ -1,5 +1,11 @@
from __future__ import absolute_import
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from future import standard_library
standard_library.install_aliases()
import numpy as np, SimPEG as simpeg
from MT1Danalytic import getEHfields
from .MT1Danalytic import getEHfields
from scipy.constants import mu_0
def get1DEfields(m1d,sigma,freq,sourceAmp=1.0):
@@ -9,7 +15,7 @@ def get1DEfields(m1d,sigma,freq,sourceAmp=1.0):
G = m1d.nodalGrad
# Mass matrices
# Magnetic permeability
Mmu = simpeg.Utils.sdiag(m1d.vol*(1.0/mu_0))
Mmu = simpeg.Utils.sdiag(m1d.vol*(old_div(1.0,mu_0)))
# Conductivity
Msig = m1d.getFaceInnerProduct(sigma)
# Set up the solution matrix
@@ -23,7 +29,7 @@ def get1DEfields(m1d,sigma,freq,sourceAmp=1.0):
Ed, Eu, Hd, Hu = getEHfields(m1d,sigma,freq,m1d.vectorNx)
Etot = (Ed + Eu)
if sourceAmp is not None:
Etot = ((Etot/Etot[-1])*sourceAmp) # Scale the fields to be equal to sourceAmp at the top
Etot = ((old_div(Etot,Etot[-1]))*sourceAmp) # Scale the fields to be equal to sourceAmp at the top
## Note: The analytic solution is derived with e^iwt
bc = np.r_[Etot[0],Etot[-1]]
# The right hand side
+10 -4
View File
@@ -1,4 +1,10 @@
from MT1Dsolutions import * # Add the names of the functions
from MT1Danalytic import *
from dataUtils import *
from ediFilesUtils import *
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from .MT1Dsolutions import * # Add the names of the functions
from .MT1Danalytic import *
from .dataUtils import *
from .ediFilesUtils import *
+21 -15
View File
@@ -1,3 +1,9 @@
from __future__ import print_function
from __future__ import absolute_import
from __future__ import division
from __future__ import unicode_literals
from future import standard_library
standard_library.install_aliases()
# Utils used for the data,
import numpy as np, matplotlib.pyplot as plt, sys
import SimPEG as simpeg
@@ -19,7 +25,7 @@ def getAppRes(MTdata):
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(MTdata, rotAngle):
'''
Function that rotates clockwist by rotAngle (- negative for a counter-clockwise rotation)
'''
@@ -44,19 +50,19 @@ def rotateData(MTdata,rotAngle):
return MT.Data.fromRecArray(outRec)
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)
def appResPhs(freq, z):
app_res = (old_div((old_div(1.,(8e-7*np.pi**2))),freq))*np.abs(z)**2
app_phs = np.arctan2(z.imag,z.real)*(old_div(180,np.pi))
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 )))))
return np.sqrt( (rho*((old_div(1,(freq * mu_0 * np.pi ))))))
def rec2ndarr(x,dt=float):
def rec2ndarr(x, dt=float):
return x.view((dt, len(x.dtype.names)))
def makeAnalyticSolution(mesh,model,elev,freqs):
def makeAnalyticSolution(mesh, model, elev, freqs):
from SimPEG import MT
data1D = []
for freq in freqs:
@@ -64,13 +70,13 @@ def makeAnalyticSolution(mesh,model,elev,freqs):
anaE = anaEd+anaEu
anaH = anaHd+anaHu
anaZ = anaE/anaH
anaZ = old_div(anaE,anaH)
# Add to the list
data1D.append((freq,0,0,elev,anaZ[0]))
dataRec = np.array(data1D,dtype=[('freq',float),('x',float),('y',float),('z',float),('zyx',complex)])
return dataRec
def plotMT1DModelData(problem,models,symList=None):
def plotMT1DModelData(problem, models, symList=None):
from SimPEG import MT
# Setup the figure
fontSize = 15
@@ -97,7 +103,7 @@ def plotMT1DModelData(problem,models,symList=None):
# if not symList:
# symList = ['x']*len(models)
import plotDataTypes as pDt
from . import plotDataTypes as pDt
# Loop through the models.
modelList = [problem.survey.mtrue]
modelList.extend(models)
@@ -110,14 +116,14 @@ def plotMT1DModelData(problem,models,symList=None):
else:
data1D = problem.dataPair(problem.survey,problem.survey.dpred(model)).toRecArray('Complex')
# Plot the data and the model
colRat = nr/((len(modelList)-1.999)*1.)
colRat = old_div(nr,((len(modelList)-1.999)*1.))
if colRat > 1.:
col = 'k'
else:
col = plt.cm.seismic(1-colRat)
# The model - make the pts to plot
meshPts = np.concatenate((problem.mesh.gridN[0:1],np.kron(problem.mesh.gridN[1::],np.ones(2))[:-1]))
modelPts = np.kron(1./(problem.mapping.sigmaMap*model),np.ones(2,))
modelPts = np.kron(old_div(1.,(problem.mapping.sigmaMap*model)),np.ones(2,))
axM.semilogx(modelPts,meshPts,color=col)
## Data
@@ -144,7 +150,7 @@ def plotMT1DModelData(problem,models,symList=None):
# Fix labels and ticks
yMtick = [l/1000 for l in axM.get_yticks().tolist()]
yMtick = [old_div(l,1000) for l in axM.get_yticks().tolist()]
axM.set_yticklabels(yMtick)
[ l.set_rotation(90) for l in axM.get_yticklabels()]
[ l.set_rotation(90) for l in axR.get_yticklabels()]
@@ -157,7 +163,7 @@ def plotMT1DModelData(problem,models,symList=None):
def printTime():
import time
print time.strftime("%a, %d %b %Y %H:%M:%S +0000", time.localtime())
print(time.strftime("%a, %d %b %Y %H:%M:%S +0000", time.localtime()))
def convert3Dto1Dobject(MTdata,rxType3D='zyx'):
from SimPEG import MT
+13 -4
View File
@@ -1,3 +1,12 @@
from __future__ import print_function
from __future__ import division
from __future__ import unicode_literals
from __future__ import absolute_import
from builtins import open
from builtins import int
from future import standard_library
standard_library.install_aliases()
from builtins import object
# Functions to import and export MT EDI files.
from SimPEG import mkvc
from scipy.constants import mu_0
@@ -9,7 +18,7 @@ import numpy as np
import os, sys, re
class EDIimporter:
class EDIimporter(object):
"""
A class to import EDIfiles.
@@ -18,7 +27,7 @@ class EDIimporter:
# Define data converters
_impUnitEDI2SI = 4*np.pi*1e-4 # Convert Z[mV/km/nT] (as in EDI)to Z[V/A] SI unit
_impUnitSI2EDI = 1./_impUnitEDI2SI # ConvertZ[V/A] SI unit to Z[mV/km/nT] (as in EDI)
_impUnitSI2EDI = old_div(1.,_impUnitEDI2SI) # ConvertZ[V/A] SI unit to Z[mV/km/nT] (as in EDI)
# Properties
filesList = None
@@ -116,7 +125,7 @@ class EDIimporter:
try:
import osr
except ImportError as e:
print 'Could not import osr, missing the gdal package\nCan not project coordinates'
print('Could not import osr, missing the gdal package\nCan not project coordinates')
raise e
# Coordinates convertor
if self._2out is None:
@@ -126,7 +135,7 @@ class EDIimporter:
if self._outEPSG is None:
# Find the UTM EPSG number
Nnr = 700 if latD < 0.0 else 600
utmZ = int(1+(longD+180.0)/6.0)
utmZ = int(1+old_div((longD+180.0),6.0))
self._outEPSG = 32000 + Nnr + utmZ
out.ImportFromEPSG(self._outEPSG)
self._2out = osr.CoordinateTransformation(src,out)
+55 -49
View File
@@ -1,3 +1,9 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from matplotlib import pyplot as plt, colors, numpy as np
@@ -35,7 +41,7 @@ def plotIsoFreqNSimpedance(ax,freq,array,flag,par='abs',colorbar=True,colorNorm=
cmap = plt.get_cmap('RdYlBu')
level = np.concatenate((-np.logspace(0,-10,31),np.logspace(-10,0,31)))
clevel = np.concatenate((-np.logspace(0,-8,5),np.logspace(-8,0,5)))
plotNorm = colors.SymLogNorm(1e-10,linscale=2)
plotNorm = colors.SymLogNorm(1e-10,linscale=2)
if cLevel:
level = np.concatenate((-np.logspace(0,-10,31),np.logspace(-10,0,31)))
clevel = np.concatenate((-np.logspace(0,-8,5),np.logspace(-8,0,5)))
@@ -65,45 +71,45 @@ def plotIsoFreqNSDiff(ax,freq,arrayList,flag,par='abs',colorbar=True,cLevel=True
x, y = arrayList[0]['x'][indUniFreq0],arrayList[0]['y'][indUniFreq0]
if par == 'abs':
if useLog:
zPlot = (np.log10(np.abs(arrayList[0][flag][indUniFreq0])) - np.log10(np.abs(arrayList[1][flag][indUniFreq1])))/np.log10(np.abs(arrayList[1][flag][indUniFreq1]))
zPlot = old_div((np.log10(np.abs(arrayList[0][flag][indUniFreq0])) - np.log10(np.abs(arrayList[1][flag][indUniFreq1]))),np.log10(np.abs(arrayList[1][flag][indUniFreq1])))
else:
zPlot = (np.abs(arrayList[0][flag][indUniFreq0]) - np.abs(arrayList[1][flag][indUniFreq1]))/np.abs(arrayList[1][flag][indUniFreq1])
zPlot = old_div((np.abs(arrayList[0][flag][indUniFreq0]) - np.abs(arrayList[1][flag][indUniFreq1])),np.abs(arrayList[1][flag][indUniFreq1]))
if mask:
maskInd = np.logical_or(np.abs(arrayList[0][flag][indUniFreq0])< 1e-3,np.abs(arrayList[1][flag][indUniFreq1]) < 1e-3)
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
if cLevel:
level = np.arange(-200,201,10)
level = np.arange(-200,201,10)
clevel = np.arange(-200,201,25)
else:
level = np.linspace(zPlot.min(),zPlot.max(),100)
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
elif par == 'real':
if useLog:
zPlot = (np.log10(np.real(arrayList[0][flag][indUniFreq0])) -np.log10(np.real(arrayList[1][flag][indUniFreq1])))/np.log10(np.abs((np.real(arrayList[1][flag][indUniFreq1]))))
zPlot = old_div((np.log10(np.real(arrayList[0][flag][indUniFreq0])) -np.log10(np.real(arrayList[1][flag][indUniFreq1]))),np.log10(np.abs((np.real(arrayList[1][flag][indUniFreq1])))))
else:
zPlot = (np.real(arrayList[0][flag][indUniFreq0]) -np.real(arrayList[1][flag][indUniFreq1]))/np.abs((np.real(arrayList[1][flag][indUniFreq1])))
zPlot = old_div((np.real(arrayList[0][flag][indUniFreq0]) -np.real(arrayList[1][flag][indUniFreq1])),np.abs((np.real(arrayList[1][flag][indUniFreq1]))))
if mask:
maskInd = np.logical_or(np.abs(np.real(arrayList[0][flag][indUniFreq0])) < 1e-3,np.abs(np.real(arrayList[1][flag][indUniFreq1])) < 1e-3)
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
if cLevel:
level = np.arange(-200,201,10)
level = np.arange(-200,201,10)
clevel = np.arange(-200,201,25)
else:
level = np.linspace(zPlot.min(),zPlot.max(),100)
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
elif par == 'imag':
if useLog:
zPlot = (np.log10(np.imag(arrayList[0][flag][indUniFreq0])) -np.log10(np.imag(arrayList[1][flag][indUniFreq1])))/np.log10(np.abs((np.imag(arrayList[1][flag][indUniFreq1]))))
zPlot = old_div((np.log10(np.imag(arrayList[0][flag][indUniFreq0])) -np.log10(np.imag(arrayList[1][flag][indUniFreq1]))),np.log10(np.abs((np.imag(arrayList[1][flag][indUniFreq1])))))
else:
zPlot = (np.imag(arrayList[0][flag][indUniFreq0]) -np.imag(arrayList[1][flag][indUniFreq1]))/np.abs((np.imag(arrayList[1][flag][indUniFreq1])))
zPlot = old_div((np.imag(arrayList[0][flag][indUniFreq0]) -np.imag(arrayList[1][flag][indUniFreq1])),np.abs((np.imag(arrayList[1][flag][indUniFreq1]))))
if mask:
maskInd = np.logical_or(np.abs(np.imag(arrayList[0][flag][indUniFreq0])) < 1e-3,np.abs(np.imag(arrayList[1][flag][indUniFreq1])) < 1e-3)
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
if cLevel:
level = np.arange(-200,201,10)
level = np.arange(-200,201,10)
clevel = np.arange(-200,201,25)
else:
level = np.linspace(zPlot.min(),zPlot.max(),100)
@@ -174,13 +180,13 @@ def plotIsoFreqNStipper(ax,freq,array,flag,par='abs',colorbar=True,colorNorm='Sy
plt.colorbar(cs,cax=ax.cax,ticks=clevel,format='%1.2e')
ax.set_title(flag+' '+par,fontsize=8)
def plotIsoStaImpedance(ax,loc,array,flag,par='abs',pSym='s',pColor=None):
appResFact = 1/(8*np.pi**2*10**(-7))
def plotIsoStaImpedance(ax,loc,array,flag,par='abs',pSym='s',pColor=None):
appResFact = old_div(1,(8*np.pi**2*10**(-7)))
treshold = 1.0 # 1 meter
indUniSta = np.sqrt(np.sum((rec2nd(array[['x','y']])-loc)**2,axis=1)) < treshold
freq = array['freq'][indUniSta]
if par == 'abs':
zPlot = np.abs(array[flag][indUniSta])
elif par == 'real':
@@ -188,9 +194,9 @@ def plotIsoStaImpedance(ax,loc,array,flag,par='abs',pSym='s',pColor=None):
elif par == 'imag':
zPlot = np.imag(array[flag][indUniSta])
elif par == 'res':
zPlot = (appResFact/freq)*np.abs(array[flag][indUniSta])**2
zPlot = (old_div(appResFact,freq))*np.abs(array[flag][indUniSta])**2
elif par == 'phs':
zPlot = np.arctan2(array[flag][indUniSta].imag,array[flag][indUniSta].real)*(180/np.pi)
zPlot = np.arctan2(array[flag][indUniSta].imag,array[flag][indUniSta].real)*(old_div(180,np.pi))
if not pColor:
if 'xx' in flag:
@@ -198,7 +204,7 @@ def plotIsoStaImpedance(ax,loc,array,flag,par='abs',pSym='s',pColor=None):
pColor = 'g'
elif 'xy' in flag:
lab = 'XY'
pColor = 'r'
pColor = 'r'
elif 'yx' in flag:
lab = 'YX'
pColor = 'b'
@@ -207,14 +213,14 @@ def plotIsoStaImpedance(ax,loc,array,flag,par='abs',pSym='s',pColor=None):
pColor = 'y'
ax.plot(freq,zPlot,color=pColor,marker=pSym,label=flag)
def plotPsudoSectNSimpedance(ax,sectDict,array,flag,par='abs',colorbar=True,colorNorm='None',cLevel=None,contour=True):
indSect = np.where(sectDict.values()[0]==array[sectDict.keys()[0]])
indSect = np.where(list(sectDict.values())[0]==array[list(sectDict.keys())[0]])
# Define the plot axes
if 'x' in sectDict.keys()[0]:
if 'x' in list(sectDict.keys())[0]:
x = array['y'][indSect]
else:
x = array['x'][indSect]
@@ -229,9 +235,9 @@ def plotPsudoSectNSimpedance(ax,sectDict,array,flag,par='abs',colorbar=True,colo
else:
level = np.linspace(zPlot.min(),zPlot.max(),100,endpoint=True)
clevel = np.linspace(zPlot.min(),zPlot.max(),10,endpoint=True)
elif par == 'ares':
zPlot = np.abs(array[flag][indSect])**2/(8*np.pi**2*10**(-7)*array['freq'][indSect])
elif par == 'ares':
zPlot = old_div(np.abs(array[flag][indSect])**2,(8*np.pi**2*10**(-7)*array['freq'][indSect]))
cmap = plt.get_cmap('RdYlBu')#seismic)
if cLevel:
zMax = np.log10(cLevel[1])
@@ -244,7 +250,7 @@ def plotPsudoSectNSimpedance(ax,sectDict,array,flag,par='abs',colorbar=True,colo
plotNorm = colors.LogNorm()
elif par == 'aphs':
zPlot = np.arctan2(array[flag][indSect].imag,array[flag][indSect].real)*(180/np.pi)
zPlot = np.arctan2(array[flag][indSect].imag,array[flag][indSect].real)*(old_div(180,np.pi))
cmap = plt.get_cmap('RdYlBu')#seismic)
if cLevel:
zMax = cLevel[1]
@@ -271,7 +277,7 @@ def plotPsudoSectNSimpedance(ax,sectDict,array,flag,par='abs',colorbar=True,colo
elif par == 'imag':
zPlot = np.imag(array[flag][indSect])
cmap = plt.get_cmap('Spectral') #('RdYlBu')
if cLevel:
zMax = np.log10(cLevel[1])
zMin = np.log10(cLevel[0])
@@ -281,7 +287,7 @@ def plotPsudoSectNSimpedance(ax,sectDict,array,flag,par='abs',colorbar=True,colo
level = np.concatenate((-np.logspace(zMax,zMin-.125,(zMax-zMin)*8+1,endpoint=True),np.logspace(zMin-.125,zMax,(zMax-zMin)*8+1,endpoint=True)))
clevel = np.concatenate((-np.logspace(zMax,zMin,(zMax-zMin)*1+1,endpoint=True),np.logspace(zMin,zMax,(zMax-zMin)*1+1,endpoint=True)))
plotNorm = colors.SymLogNorm(np.abs(level).min(),linscale=0.1)
if colorNorm=='SymLog':
plotNorm = colors.SymLogNorm(np.abs(level).min(),linscale=0.1)
elif colorNorm=='Lin':
@@ -307,14 +313,14 @@ def plotPsudoSectNSDiff(ax,sectDict,arrayList,flag,par='abs',colorbar=True,color
def sortInArr(arr):
return np.sort(arr,order=['freq','x','y','z'])
# Find the index for the slice
indSect0 = np.where(sectDict.values()[0]==arrayList[0][sectDict.keys()[0]])
indSect1 = np.where(sectDict.values()[0]==arrayList[1][sectDict.keys()[0]])
indSect0 = np.where(list(sectDict.values())[0]==arrayList[0][list(sectDict.keys())[0]])
indSect1 = np.where(list(sectDict.values())[0]==arrayList[1][list(sectDict.keys())[0]])
# Extract and sort the mats
arr0 = sortInArr(arrayList[0][indSect0])
arr1 = sortInArr(arrayList[1][indSect1])
# Define the plot axes
if 'x' in sectDict.keys()[0]:
if 'x' in list(sectDict.keys())[0]:
x0 = arr0['y']
x1 = arr1['y']
else:
@@ -322,35 +328,35 @@ def plotPsudoSectNSDiff(ax,sectDict,arrayList,flag,par='abs',colorbar=True,color
x1 = arr1['x']
y0 = arr0['freq']
y1 = arr1['freq']
if par == 'abs':
if useLog:
zPlot = (np.log10(np.abs(arr0[flag])) - np.log10(np.abs(arr1[flag])))/np.log10(np.abs(arr1[flag]))
zPlot = old_div((np.log10(np.abs(arr0[flag])) - np.log10(np.abs(arr1[flag]))),np.log10(np.abs(arr1[flag])))
else:
zPlot = (np.abs(arr0[flag]) - np.abs(arr1[flag]))/np.abs(arr1[flag])
zPlot = old_div((np.abs(arr0[flag]) - np.abs(arr1[flag])),np.abs(arr1[flag]))
if mask:
maskInd = np.logical_or(np.abs(arr0[flag])< 1e-3,np.abs(arr1[flag]) < 1e-3)
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
cmap = plt.get_cmap('RdYlBu')#seismic)
elif par == 'ares':
arF = 1/(8*np.pi**2*10**(-7))
arF = old_div(1,(8*np.pi**2*10**(-7)))
if useLog:
zPlot = (np.log10((arF/arr0['freq'])*np.abs(arr0[flag])**2) - np.log10((arF/arr1['freq'])*np.abs(arr1[flag])**2))/np.log10((arF/arr1['freq'])*np.abs(arr1[flag])**2)
zPlot = old_div((np.log10((old_div(arF,arr0['freq']))*np.abs(arr0[flag])**2) - np.log10((old_div(arF,arr1['freq']))*np.abs(arr1[flag])**2)),np.log10((old_div(arF,arr1['freq']))*np.abs(arr1[flag])**2))
else:
zPlot = ((arF/arr0['freq'])*np.abs(arr0[flag])**2 - (arF/arr1['freq'])*np.abs(arr1[flag])**2)/((arF/arr1['freq'])*np.abs(arr1[flag])**2)
zPlot = old_div(((old_div(arF,arr0['freq']))*np.abs(arr0[flag])**2 - (old_div(arF,arr1['freq']))*np.abs(arr1[flag])**2),((old_div(arF,arr1['freq']))*np.abs(arr1[flag])**2))
if mask:
maskInd = np.logical_or(np.abs(arr0[flag])< 1e-3,np.abs(arr1[flag]) < 1e-3)
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
cmap = plt.get_cmap('Spectral')#seismic)
elif par == 'aphs':
if useLog:
zPlot = (np.log10(np.arctan2(arr0[flag].imag,arr0[flag].real)*(180/np.pi)) - np.log10(np.arctan2(arr1[flag].imag,arr1[flag].real)*(180/np.pi)) )/np.log10(np.arctan2(arr1[flag].imag,arr1[flag].real)*(180/np.pi))
zPlot = old_div((np.log10(np.arctan2(arr0[flag].imag,arr0[flag].real)*(old_div(180,np.pi))) - np.log10(np.arctan2(arr1[flag].imag,arr1[flag].real)*(old_div(180,np.pi))) ),np.log10(np.arctan2(arr1[flag].imag,arr1[flag].real)*(old_div(180,np.pi))))
else:
zPlot = ( np.arctan2(arr0[flag].imag,arr0[flag].real)*(180/np.pi) - np.arctan2(arr1[flag].imag,arr1[flag].real)*(180/np.pi) )/(np.arctan2(arr1[flag].imag,arr1[flag].real)*(180/np.pi))
zPlot = old_div(( np.arctan2(arr0[flag].imag,arr0[flag].real)*(old_div(180,np.pi)) - np.arctan2(arr1[flag].imag,arr1[flag].real)*(old_div(180,np.pi)) ),(np.arctan2(arr1[flag].imag,arr1[flag].real)*(old_div(180,np.pi))))
if mask:
maskInd = np.logical_or(np.abs(arr0[flag])< 1e-3,np.abs(arr1[flag]) < 1e-3)
zPlot = np.ma.array(zPlot)
@@ -358,26 +364,26 @@ def plotPsudoSectNSDiff(ax,sectDict,arrayList,flag,par='abs',colorbar=True,color
cmap = plt.get_cmap('Spectral')#seismic)
elif par == 'real':
if useLog:
zPlot = (np.log10(arr0[flag].real) - np.log10(arr1[flag].real))/np.log10(arr1[flag].real)
zPlot = old_div((np.log10(arr0[flag].real) - np.log10(arr1[flag].real)),np.log10(arr1[flag].real))
else:
zPlot = (arr0[flag].real - arr1[flag].real)/arr1[flag].real
zPlot = old_div((arr0[flag].real - arr1[flag].real),arr1[flag].real)
if mask:
maskInd = np.logical_or(arr0[flag].real< 1e-3,arr1[flag].real < 1e-3)
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
cmap = plt.get_cmap('Spectral') #('Spectral')
elif par == 'imag':
if useLog:
zPlot = (np.log10(arr0[flag].imag) - np.log10(arr1[flag].imag))/np.log10(arr1[flag].imag)
zPlot = old_div((np.log10(arr0[flag].imag) - np.log10(arr1[flag].imag)),np.log10(arr1[flag].imag))
else:
zPlot = (arr0[flag].imag - arr1[flag].imag)/arr1[flag].imag
zPlot = old_div((arr0[flag].imag - arr1[flag].imag),arr1[flag].imag)
if mask:
maskInd = np.logical_or(arr0[flag].imag< 1e-3,arr1[flag].imag < 1e-3)
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
cmap = plt.get_cmap('Spectral') #('RdYlBu')
if cLevel:
zMax = np.log10(cLevel[1])
zMin = np.log10(cLevel[0])
@@ -385,18 +391,18 @@ def plotPsudoSectNSDiff(ax,sectDict,arrayList,flag,par='abs',colorbar=True,color
zMax = (np.ceil(np.log10(np.abs(zPlot).max())))
zMin = (np.floor(np.log10(np.abs(zPlot).min())))
if colorNorm=='SymLog':
level = np.concatenate((-np.logspace(zMax,zMin-.125,(zMax-zMin)*8+1,endpoint=True),np.logspace(zMin-.125,zMax,(zMax-zMin)*8+1,endpoint=True)))
clevel = np.concatenate((-np.logspace(zMax,zMin,(zMax-zMin)*1+1,endpoint=True),np.logspace(zMin,zMax,(zMax-zMin)*1+1,endpoint=True)))
plotNorm = colors.SymLogNorm(np.abs(level).min(),linscale=0.1)
elif colorNorm=='Lin':
if cLevel:
level = np.arange(cLevel[0],cLevel[1]+.1,(cLevel[1] - cLevel[0])/50.)
clevel = np.arange(cLevel[0],cLevel[1]+.1,(cLevel[1] - cLevel[0])/10.)
level = np.arange(cLevel[0],cLevel[1]+.1,old_div((cLevel[1] - cLevel[0]),50.))
clevel = np.arange(cLevel[0],cLevel[1]+.1,old_div((cLevel[1] - cLevel[0]),10.))
else:
level = np.arange(zPlot.min(),zPlot.max(),(zPlot.max() - zPlot.min())/50.)
clevel = np.arange(zPlot.min(),zPlot.max(),(zPlot.max() - zPlot.min())/10.)
level = np.arange(zPlot.min(),zPlot.max(),old_div((zPlot.max() - zPlot.min()),50.))
clevel = np.arange(zPlot.min(),zPlot.max(),old_div((zPlot.max() - zPlot.min()),10.))
plotNorm = colors.Normalize()
elif colorNorm=='Log':
level = np.logspace(zMin-.125,zMax,(zMax-zMin)*8+1,endpoint=True)

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