Compare commits

..
139 Commits
Author SHA1 Message Date
D Fournier 71098a01b5 Allow to forward G data without storing the kernel. Add figure-axes options to plot data function 2016-07-12 06:36:25 -07:00
D Fournier b13b36dcb8 Small updates to the function call 2016-07-01 13:51:47 -07:00
D Fournier 66ee0f0422 Merge branch 'dev' into pf/dev
Conflicts:
	.gitignore
	SimPEG/Examples/__init__.py
2016-06-30 13:55:57 -07:00
Rowan Cockett 8e944078c8 Merged branch dev into pf/dev 2016-06-26 16:25:13 -06:00
Rowan Cockett 64108eb36d Updates to the asset folder location. 2016-06-22 12:38:19 -06:00
Rowan Cockett a79fc7456d Change headers in docs. 2016-06-22 12:27:05 -06:00
Rowan Cockett 4074c67bcf Initial merge and minor refactor of simpegPF. 2016-06-22 12:25:08 -06:00
Rowan Cockett 0d8fa20257 Minor path cleanup. 2016-06-22 11:30:46 -06:00
Rowan Cockett 8d3a710d92 Updates to testing the boundary conditions. 2016-06-22 11:29:04 -06:00
Rowan Cockett bf2035425c Add ploting code. 2016-06-22 11:26:17 -06:00
Rowan Cockett 31e7bcb3f0 Analytic test 2016-06-22 11:24:46 -06:00
Rowan Cockett f216b6a748 Name change for problem. 2016-06-22 11:19:11 -06:00
Rowan Cockett 9336b821da Delete MagSurveyIE.py 2016-06-22 11:16:51 -06:00
Rowan Cockett d7dad2fe60 Delete Ex_sir_Wathen.py 2016-06-22 11:15:54 -06:00
D Fournier 7dd0b1b2df Latest updates in preparation for merge with SimPEG 2016-06-22 09:59:50 -07:00
D Fournier 375a18c53a Add Active cell input file. Allow for Inactive, Active-Static and Active-Dynamic cells [0, 1, -1]
Fix MAG and GRAV test scripts.
2016-06-06 12:29:23 -07:00
D Fournier 7f7ce5e25b Add driver for gravity. Adapt example for new Regularization formulation 2016-05-30 17:35:47 -07:00
D Fournier f1eaac7cd0 Small changes to Example 2016-05-29 14:14:00 -07:00
D Fournier 90fb125521 Adapt MAG code for refactor IRLS iterations 2016-05-28 11:25:27 -07:00
D Fournier e4998d0bb9 Update example for constrained and compact-constrained MAG inversion 2016-05-05 12:37:03 -07:00
D Fournier fd5e56202a New Example script to have inactive cells in the inversion 2016-05-04 16:02:28 -07:00
Rowan Cockett c32abd215c Merge pull request #7 from simpeg/ref/merge
Driver File
2016-04-29 17:16:53 -07:00
Rowan Cockett 64610d3a0a Merge branch 'Dom' of https://github.com/simpeg/simpegPF into ref/merge
# Conflicts:
#	simpegPF/Dev/MAG/Intgrl_MAG_Inv_Driver.py
2016-04-29 17:12:57 -07:00
Rowan Cockett b01b14efb9 Rename to mad driver 2016-04-29 16:58:09 -07:00
Rowan Cockett 8bbc70228d Rename to Problem3D_Integral 2016-04-29 16:06:08 -07:00
Rowan Cockett 47cf1151e3 updates to headers. 2016-04-29 16:03:14 -07:00
Rowan Cockett 7ad1b6f5aa put main script at end. 2016-04-29 15:53:08 -07:00
Rowan Cockett 73b6316962 Separate out the input files. 2016-04-29 15:49:49 -07:00
D Fournier 6a3d0e4ffe Updated drivers. Clean version 2016-04-29 15:08:12 -07:00
Rowan Cockett a45e1fb0fd Remove gocad and surfaces, use mesh_io utils in SimPEG 2016-04-29 14:22:00 -07:00
D Fournier 363236394b Update Directives in Tutotial MAG 2016-04-29 13:04:57 -07:00
D Fournier c40b9a9add Merge branch 'Dom'
Conflicts:
	simpegPF/Dev/model.png
	simpegPF/Dev/obspred.png
2016-04-29 11:05:45 -07:00
D Fournier 4c1a893109 Add Update_Beta directive. Minor changes to the example 2016-04-29 11:03:05 -07:00
D Fournier 5bd2805e45 Add gravity example. Add Mag_parametric script 2016-04-22 15:24:22 -07:00
D Fournier 3bb71f0fd4 Update grav code and example 2016-04-18 12:53:11 -07:00
D Fournier c61da61061 Change u for f (field) in Jvec and Jtvec.
Update the Notebook in preperation for example.
2016-04-07 10:23:40 -07:00
D Fournier 0578e1ff6c Change to source parameters
Add color axis handles for data plot
2016-04-03 17:23:24 -07:00
D Fournier 03dd9915b7 First commit for GRAV3D
Add example for mag
Add example for grav
2016-03-16 17:10:30 -07:00
D Fournier c70256d028 Minor change to lp norms 2016-03-16 11:54:56 -07:00
D Fournier b55897453e Create example for Craig and add Driver script 2016-03-16 11:53:27 -07:00
seogi_macbook 241dc83354 Merge branch 'master' of https://github.com/simpeg/simpegpf
Conflicts:
	simpegPF/notebooks/MagInversion.ipynb
	simpegPF/notebooks/tutorials/Tutorial_1_Mag forward modeling.ipynb
2016-03-15 00:13:52 -07:00
seogi_macbook 0cbda152ec temporary commit 2016-03-15 00:13:07 -07:00
D Fournier ceffc0b53f Major overhaul of Survey and Problem class.
Create SrcField and RxObs
retrotfit I/O to use survey
Forward matrix as property
Inversion running again but needs cleanup
2016-03-14 22:06:24 -07:00
D Fournier 000e227ea4 Minor change 2016-03-06 16:56:39 -08:00
D Fournier d2db56337a Updated notebook example 2016-03-06 16:53:46 -08:00
D Fournier c977a16b9e Adapt Mag Problem for topography 2016-02-25 08:57:55 -08:00
D Fournier 6ba2cb549c Clean up examples 2016-02-09 08:13:37 -08:00
D Fournier 516522deda Save example for next time 2016-01-31 15:51:08 -08:00
D Fournier bb9624f4ff Implement sparse norm on mag problem.
Add Nutcracker example to folder.
2016-01-31 15:33:57 -08:00
D Fournier 90838aa74b Implement sparse norm and test on nutcracker 2016-01-29 00:56:39 -08:00
D Fournier d718ffd3a3 Progress on 3D inversion using simple regularization. 2016-01-28 18:39:13 -08:00
D Fournier 45a5ab556e no message 2016-01-19 15:26:14 -08:00
D Fournier 6f508b0d27 Investigating problem with GNCG 2016-01-19 15:19:06 -08:00
D Fournier b2fc35ee75 Test script for multiple *.ts surfaces to block model (Lalor model) 2016-01-18 17:10:31 -08:00
D Fournier 6a7d72bf85 Work on Example and mesh refinement for MinSIM 2016-01-17 13:39:10 -08:00
D Fournier c3c96a5bd8 VTK example for surface to mesh model - Thanks Gudni! 2016-01-15 16:49:27 -08:00
D Fournier 44b7dd184f Fix merge conflicts and add read GOCAD triangulated surface file 2016-01-15 12:38:08 -08:00
D Fournier dc7b8607eb Fix merge conflicts 2016-01-15 10:56:30 -08:00
D Fournier e4e19ceac1 Merge branch 'Dom' of https://github.com/simpeg/simpegpf into Dom
Conflicts:
	simpegPF/Magnetics.py
2016-01-15 00:03:37 -08:00
D Fournier c0445e5db5 Modify Inv Integral function
Start example in a notebook
2016-01-15 00:01:24 -08:00
Rowan Cockett dbcd57bc4e Minor clean ups. 2016-01-11 18:35:19 -08:00
D Fournier 27900498c6 Add surface to mesh test
Fix forward function call
2016-01-11 16:40:28 -08:00
D Fournier 958c4d151a Test interpolation methods.
Attempt at using FFT to resample grid.
Need to can an FFT method for 2 grid ( pad, taper, resample, extract)
2016-01-04 12:29:28 -08:00
D Fournier 201fe6bf46 Implement active cell from topo and topocheck function 2015-12-29 16:11:16 -08:00
D Fournier 1a7656c34d Implement distance weighting for potential fields in BaseMag 2015-12-27 11:29:01 -08:00
D Fournier cbe3b3de42 Move sub functions to PF.Magnetics
Add Mag Integral forward operator
Begin Inversion script for magnetic integral
2015-12-21 13:11:46 -08:00
D Fournier d6a2490a7b Create function to generate linear forward operator (F)
Clean project
2015-12-21 08:42:50 -08:00
D Fournier e9f02e1d77 Finish test with analytical answer. O(h^2) error pass
Add function MagSphereFreeSpace in MagAnalytic. Same results as MagSphereAnaFunA but with z-observation flipped.
2015-11-27 16:08:21 -08:00
D Fournier 291cce9d52 Finish forward modeling in integral
Start test function.
2015-11-13 12:49:19 -08:00
D Fournier 6c71237328 Scripts added for the magnetic integral formulation 2015-11-10 14:23:29 -08:00
D Fournier 6174c6fece Start working on mag integral 2015-10-23 12:08:19 -07:00
D Fournier 7d516d2879 Change functions in the header and start a development branch. 2015-09-27 09:25:53 -07:00
Lindsey 80c839390c updated travis 2015-08-07 09:51:22 -07:00
Lindsey Heagy f416dfacf8 added myself to travis 2015-08-06 20:45:18 -07:00
Lindsey 3ae311c7b1 merge conflicts resolved 2015-08-06 16:42:38 -07:00
Lindsey 203640c803 for the sphere problem, the magnetic field values need to be calculated wrt the center of the sphere. analytics abbreviate to Ana 2015-08-06 16:40:19 -07:00
seogi_macbook 97d00effb6 Merge branch 'master' of https://github.com/simpeg/simpegpf 2015-05-15 16:04:20 -07:00
seogi_macbook e05221c6ea Modification for meshes 2015-05-15 16:03:41 -07:00
Lindsey Heagy 6afff26e04 abbreviate Analytic Ana 2015-05-10 17:04:34 -07:00
seogi 668e435cce Tutorial for Mag forward problem:
TOODs:

- Change code that we can implement SimPEG's solver class (Now use jacobi-bicgstab)
- Put linear forward problem
2015-01-28 10:58:07 -08:00
SEOGI KANG c769ac6d8d modify docs 2014-07-08 14:17:38 -07:00
SEOGI KANG ab1e639fbd Documentations 2014-07-04 22:29:28 -07:00
seogi d85ff18bd8 ss 2014-07-04 16:54:03 -07:00
seogi 0812f14acd modification for no congrous bc 2014-07-04 16:53:22 -07:00
seogi d1c34da5f0 archive 2014-07-04 16:50:00 -07:00
seogi d43caf092b Modifications for inversion 2014-07-04 14:34:39 -07:00
rowanc1 83a7a7bc12 name it correctly! 2014-07-03 13:45:19 -07:00
rowanc1 124108f84d add coverage in readme 2014-07-03 13:44:50 -07:00
rowanc1 64c9e8fce0 Coverage updates 2014-05-19 11:40:02 -07:00
Rowan Cockett 640457cf0b Merge pull request #5 from simpeg/MapChanges
Updates
2014-05-19 11:35:20 -07:00
rowanc1 e334d3efa9 Updates 2014-05-19 11:34:14 -07:00
Rowan Cockett 55f9f43c0e Merge pull request #4 from simpeg/MapChanges
Map changes
2014-04-15 16:03:18 -07:00
rowanc1 95f0cb9972 docs updates 2014-04-15 15:41:44 -07:00
rowanc1 1d5c3ced39 Fix derivatives. 2014-04-15 15:07:00 -07:00
rowanc1 cbefdd2e84 derivatives not yet working. 2014-04-15 14:59:25 -07:00
Lindsey Heagy e6e9616197 a couple typo fixes in the documentation 2014-04-05 11:20:14 -07:00
Lindsey Heagy a5f57d3a56 fixes a couple of minor typos 2014-04-05 10:32:50 -07:00
seogi 4b9304e476 Mag inversion 2014-03-13 12:33:18 -07:00
Rowan Cockett 4b7d5e1d32 Update api_PF.rst
Get rid of conflict
2014-03-12 23:08:58 -07:00
SEOGI KANG 0f63b0d6b6 Mag example 2014-03-12 22:34:20 -07:00
rowanc1 68021e0a42 Changed Data --> Survey 2014-03-07 14:07:47 -08:00
seogi 3f7a521833 Working codes 2014-03-03 10:50:26 -08:00
seogi 9c22a43759 commit test 2014-03-03 09:51:07 -08:00
seogi fbc69d4709 Commit test 2014-03-03 09:50:55 -08:00
seogi d3fc3d2058 test commit 2014-03-03 09:47:14 -08:00
SEOGI KANG 699c91f6d7 Initial work Mag inversion 2014-03-03 08:53:19 -08:00
seogi 18546d038b resolved conflicts 2014-02-28 16:44:12 -08:00
rowanc1 5816f75728 Hopefully fix problems with read the docs. 2014-02-27 15:26:10 -08:00
rowanc1 f83c84da81 Try fixing RTD 2014-02-27 09:33:48 -08:00
seogi ea0c4719fe added play notebook! 2014-02-27 09:15:41 -08:00
seogi 4c8741cd2e added git ignore 2014-02-27 09:15:41 -08:00
seogi 751d7098b6 Fix examples folder 2014-02-27 09:15:41 -08:00
SEOGI KANG 3cba7845ca Little work for jacobian 2014-02-27 08:10:37 -08:00
seogi bcfefde2da modified: simpegPF/notebooks/Jacobian.ipynb
modified:   simpegPF/Tests/test_forward_PFproblem.py
2014-02-26 17:28:48 -08:00
seogi 8232f94fcf Working for jacobian 2014-02-26 17:26:21 -08:00
seogi 004e225844 Documentation Mag 2014-02-26 16:56:11 -08:00
rowanc1 a253adcf4e slightly more complicated... lots of ToDos 2014-02-25 23:22:37 -08:00
rowanc1 c84a112ebc initial writing of the sensitivity (J) 2014-02-25 22:14:35 -08:00
rowanc1 fcb63d5078 make a positive definite 2014-02-25 20:57:48 -08:00
rowanc1 1fc9c0407c fixed plotting in documentation 2014-02-24 22:07:58 -08:00
rowanc1 a576c566d0 fixed imports 2014-02-24 22:00:56 -08:00
seogi e385b55674 First documentation 2014-02-24 16:59:00 -08:00
rowanc1 e662276cc1 mag work, forward problem running 2014-02-24 10:40:17 -08:00
rowanc1 5c8b6f4a21 updates to initial structure 2014-02-23 13:23:16 -08:00
rowanc1 3b97872601 Initial work on bringing mag into simpeg framework 2014-02-23 13:18:02 -08:00
SEOGI KANG 53146016ff Merge branch 'master' of https://github.com/simpeg/simpegpf 2014-02-23 12:03:54 -08:00
SEOGI KANG 0a9ed614a8 Some analytics and set up for forward problem 2014-02-23 12:00:14 -08:00
SEOGI KANG 2e6164c62e MagForwrad secondary and total 2014-02-23 11:57:51 -08:00
rowanc1 c9fa9aa23f update to develop branch 2014-02-21 21:55:27 -08:00
rowanc1 454a5f9f8c fixed travis. 2014-02-21 10:49:56 -08:00
rowanc1 6aebe56596 Merge branch 'master' of https://github.com/simpeg/simpegpf 2014-02-21 10:44:55 -08:00
rowanc1 544e69b928 update Travis 2014-02-21 10:44:49 -08:00
SEOGI KANG 35d132f52b Secondary field approach 2014-02-21 08:26:02 -08:00
seogi cc33d148b6 initial set up 2014-02-20 17:30:31 -08:00
seogi 3bc66134aa Total field approach 2014-02-20 17:29:39 -08:00
SEOGI KANG 37bca3b3e9 Some functions for neumann BC weak form 2014-02-17 11:19:44 -08:00
rowanc1 d896e99bb2 added email notifications to the build testing 2014-02-13 16:09:55 -08:00
seogi bf4e2a774f Testing inhomogenous neumann boundary conditions, but quite general
since we can put different boundary condtions on each side.

Testing has been done for homogenous case with some fictitious source test.
See: Boundary_conditions.ipynb

two functions are written.

- ddxFaceDivBC
- faceDivBC

Need to be more tested and after that they can be merged in to main
2014-02-13 13:48:57 -08:00
rowanc1 575c340130 (boiler plate) 2014-02-12 23:13:10 -08:00
225 changed files with 115551 additions and 2376 deletions
+4 -3
View File
@@ -1,7 +1,6 @@
language: python
python:
- 2.7
- 3.4
sudo: false
@@ -19,6 +18,7 @@ env:
- TEST_DIR="tests/mesh tests/base tests/utils"
- TEST_DIR=tests/em/fdem/inverse/derivs
- TEST_DIR=tests/em/tdem
- TEST_DIR=tests/pf
- TEST_DIR=tests/dcip
- TEST_DIR=tests/flow
- TEST_DIR=tests/mt
@@ -39,7 +39,7 @@ before_install:
-O miniconda.sh; fi
- chmod +x miniconda.sh
- ./miniconda.sh -b
- export PATH=/home/travis/anaconda/bin:/home/travis/anaconda3/bin:/home/travis/miniconda/bin:/home/travis/miniconda3/bin:$PATH
- export PATH=/home/travis/anaconda/bin:/home/travis/miniconda/bin:$PATH
- conda update --yes conda
install:
@@ -59,7 +59,8 @@ script:
# Calculate coverage
after_success:
- bash <(curl -s https://codecov.io/bash)
- coveralls --config_file .coveragerc
- if [ "$TRAVIS_BRANCH" = "master" -a "$TRAVIS_PULL_REQUEST" = "false" ]; then
if [ ${TEST_DIR} == "tests/docs" ]; then
python scripts/fetch_gae_sdk.py $(dirname "${GAE_PYTHONPATH}");
+1 -5
View File
@@ -28,11 +28,7 @@ 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:
+1 -9
View File
@@ -1,11 +1,3 @@
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):
@@ -69,7 +61,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
+2 -9
View File
@@ -1,12 +1,5 @@
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):
"""
@@ -59,7 +52,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
+21 -31
View File
@@ -1,16 +1,6 @@
from __future__ import print_function
from __future__ import absolute_import
from __future__ import division
from __future__ import unicode_literals
from builtins import open
from builtins import int
from future import standard_library
standard_library.install_aliases()
from builtins import map
from builtins import range
from SimPEG import np, Utils
from . import BaseDC as DC
from . import BaseDC as IP
import BaseDC as DC
import BaseDC as IP
import warnings
def getActiveindfromTopo(mesh, topo):
@@ -77,7 +67,7 @@ def readUBC_DC3Dobstopo(filename,mesh,topo,probType="CC"):
if "!" in line.split(): continue
elif line == '\n': continue
elif line == ' \n': continue
temp = list(map(float, line.split()))
temp = 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
@@ -238,10 +228,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
@@ -256,11 +246,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)]
@@ -350,11 +340,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(list(range(int(nstn))))*dl_x*AM_sep
stn_y = endl[0,1] + np.array(list(range(int(nstn))))*dl_y*AM_sep
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
# Create line of P1 locations
M = np.c_[stn_x, stn_y, np.ones(nstn).T*mesh.vectorNz[-1]]
@@ -386,15 +376,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([(AB - MN_sep) // AM_sep, nrx])
nstn = np.min([np.floor( (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(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
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
# Create receiver poles
# Create line of P1 locations
@@ -429,15 +419,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 = box_l // AM_sep
nstn = np.floor( box_l / AM_sep )
# Compute discrete pole location along line
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
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
# Define number of cross lines
nlin = int(box_w // AM_sep)
lind = list(range(-nlin,nlin+1))
nlin = int(np.floor( box_w / AM_sep ))
lind = range(-nlin,nlin+1)
ngrad = nstn * len(lind)
@@ -459,7 +449,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
@@ -678,7 +668,7 @@ def readUBC_DC3Dobs(fileName, rtype = 'DC'):
obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!')
else:
print("rtype must be 'DC'(default) | 'IP'")
print "rtype must be 'DC'(default) | 'IP'"
# Pre-allocate
srcLists = []
+1 -8
View File
@@ -1,10 +1,3 @@
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):
@@ -12,7 +5,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):
+4 -10
View File
@@ -1,10 +1,4 @@
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
from BaseDC import *
from BaseIP import *
from DCIPUtils import *
import Utils
+6 -13
View File
@@ -1,16 +1,7 @@
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
import Utils, Survey, Problem, numpy as np, scipy.sparse as sp, gc
class BaseDataMisfit(with_metaclass(Utils.SimPEGMetaClass, object)):
class BaseDataMisfit(object):
"""BaseDataMisfit
.. note::
@@ -18,6 +9,8 @@ class BaseDataMisfit(with_metaclass(Utils.SimPEGMetaClass, 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
@@ -100,11 +93,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))
+18 -29
View File
@@ -1,15 +1,4 @@
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
import Utils, numpy as np
class InversionDirective(object):
"""InversionDirective"""
@@ -26,7 +15,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
@@ -79,7 +68,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
@@ -131,7 +120,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)
@@ -152,7 +141,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
@@ -203,7 +192,7 @@ class SaveModelEveryIteration(SaveEveryIteration):
"""SaveModelEveryIteration"""
def initialize(self):
print("SimPEG.SaveModelEveryIteration will save your models as: '###-%s.npy'"%self.fileName)
print "SimPEG.SaveModelEveryIteration will save your models as: '###-%s.npy'"%self.fileName
def endIter(self):
np.save('%03d-%s' % (self.opt.iter, self.fileName), self.opt.xc)
@@ -213,7 +202,7 @@ class SaveOutputEveryIteration(SaveEveryIteration):
"""SaveModelEveryIteration"""
def initialize(self):
print("SimPEG.SaveOutputEveryIteration will save your inversion progress as: '###-%s.txt'"%self.fileName)
print "SimPEG.SaveOutputEveryIteration will save your inversion progress as: '###-%s.txt'"%self.fileName
f = open(self.fileName+'.txt', 'w')
f.write(" # beta phi_d phi_m f\n")
f.close()
@@ -227,7 +216,7 @@ class SaveOutputDictEveryIteration(SaveEveryIteration):
"""SaveOutputDictEveryIteration"""
def initialize(self):
print("SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '###-%s.npz'"%self.fileName)
print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '###-%s.npz'"%self.fileName
def endIter(self):
# Save the data.
@@ -305,25 +294,22 @@ 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
# Either use the supplied epsilon, or fix base on distribution of
# model values
if getattr(self, 'reg.eps', None) is None:
if getattr(self, 'eps', None) is None:
self.reg.eps_p = np.percentile(np.abs(self.invProb.curModel),self.prctile)
else:
self.reg.eps_p = self.eps[0]
if getattr(self, 'reg.eps', None) is None:
if getattr(self, 'eps', None) is None:
self.reg.eps_q = np.percentile(np.abs(self.reg.regmesh.cellDiffxStencil*(self.reg.mapping * self.invProb.curModel)),self.prctile)
else:
self.reg.eps_q = self.eps[1]
print("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
@@ -334,12 +320,15 @@ class Update_IRLS(InversionDirective):
self.reg.l2model = self.invProb.curModel
self.reg.curModel = self.invProb.curModel
print "L[p qx qy qz]-norm : " + str(self.reg.norms)
print "eps_p: " + str(self.reg.eps_p) + " eps_q: " + str(self.reg.eps_q)
if getattr(self, 'f_old', None) is None:
self.f_old = self.reg.eval(self.invProb.curModel)#self.invProb.evalFunction(self.invProb.curModel, return_g=False, return_H=False)
# 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
@@ -351,17 +340,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,10 +1,3 @@
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,9 +1,4 @@
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
-5
View File
@@ -1,9 +1,4 @@
from __future__ import division
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
import numpy as np
from scipy.constants import mu_0, pi, epsilon_0
from scipy.special import erf
+1 -7
View File
@@ -1,9 +1,3 @@
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
@@ -40,7 +34,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,9 +1,3 @@
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
+5 -11
View File
@@ -1,11 +1,5 @@
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 *
from TDEM import hzAnalyticDipoleT
from FDEM import hzAnalyticDipoleF
from FDEMcasing import *
from DC import DCAnalyticHalf, DCAnalyticSphere
from FDEMDipolarfields import *
-6
View File
@@ -1,9 +1,3 @@
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
+29 -36
View File
@@ -1,10 +1,3 @@
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
@@ -49,7 +42,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
:return: total electric field
"""
if getattr(self, '_ePrimary', None) is None or getattr(self, '_eSecondary', None) is None:
raise NotImplementedError ('Getting e from %s is not implemented' %list(self.knownFields.keys())[0])
raise NotImplementedError ('Getting e from %s is not implemented' %self.knownFields.keys()[0])
return self._ePrimary(solution,srcList) + self._eSecondary(solution,srcList)
@@ -63,7 +56,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
:return: total magnetic flux density
"""
if getattr(self, '_bPrimary', None) is None or getattr(self, '_bSecondary', None) is None:
raise NotImplementedError ('Getting b from %s is not implemented' %list(self.knownFields.keys())[0])
raise NotImplementedError ('Getting b from %s is not implemented' %self.knownFields.keys()[0])
return self._bPrimary(solution, srcList) + self._bSecondary(solution, srcList)
@@ -77,7 +70,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
:return: total magnetic field
"""
if getattr(self, '_hPrimary', None) is None or getattr(self, '_hSecondary', None) is None:
raise NotImplementedError ('Getting h from %s is not implemented' %list(self.knownFields.keys())[0])
raise NotImplementedError ('Getting h from %s is not implemented' %self.knownFields.keys()[0])
return self._hPrimary(solution, srcList) + self._hSecondary(solution, srcList)
@@ -91,7 +84,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
:return: total current density
"""
if getattr(self, '_jPrimary', None) is None or getattr(self, '_jSecondary', None) is None:
raise NotImplementedError ('Getting j from %s is not implemented' %list(self.knownFields.keys())[0])
raise NotImplementedError ('Getting j from %s is not implemented' %self.knownFields.keys()[0])
return self._jPrimary(solution, srcList) + self._jSecondary(solution, srcList)
@@ -107,7 +100,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
:return: derivative times a vector (or tuple for adjoint)
"""
if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None:
raise NotImplementedError ('Getting eDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0])
if adjoint:
return self._eDeriv_u(src, v, adjoint), self._eDeriv_m(src, v, adjoint)
@@ -125,7 +118,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
:return: derivative times a vector (or tuple for adjoint)
"""
if getattr(self, '_bDeriv_u', None) is None or getattr(self, '_bDeriv_m', None) is None:
raise NotImplementedError ('Getting bDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
raise NotImplementedError ('Getting bDerivs from %s is not implemented' %self.knownFields.keys()[0])
if adjoint:
return self._bDeriv_u(src, v, adjoint), self._bDeriv_m(src, v, adjoint)
@@ -143,7 +136,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
:return: derivative times a vector (or tuple for adjoint)
"""
if getattr(self, '_hDeriv_u', None) is None or getattr(self, '_hDeriv_m', None) is None:
raise NotImplementedError ('Getting hDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
raise NotImplementedError ('Getting hDerivs from %s is not implemented' %self.knownFields.keys()[0])
if adjoint:
return self._hDeriv_u(src, v, adjoint), self._hDeriv_m(src, v, adjoint)
@@ -161,7 +154,7 @@ class FieldsFDEM(SimPEG.Problem.Fields):
:return: derivative times a vector (or tuple for adjoint)
"""
if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None:
raise NotImplementedError ('Getting jDerivs from %s is not implemented' %list(self.knownFields.keys())[0])
raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0])
if adjoint:
return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint)
@@ -292,7 +285,7 @@ class Fields3D_e(FieldsFDEM):
C = self._edgeCurl
b = (C * eSolution)
for i, src in enumerate(srcList):
b[:,i] *= -1./(1j*omega(src.freq))
b[:,i] *= - 1./(1j*omega(src.freq))
s_m, _ = src.eval(self.prob)
b[:,i] = b[:,i]+ 1./(1j*omega(src.freq)) * s_m
return b
@@ -352,7 +345,7 @@ class Fields3D_e(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the derivative of the current density with respect to the field we solved for with a vector
"""
n = int(self._aveE2CCV.shape[0] // self._nC) # number of components (instead of checking if cyl or not)
n = int(self._aveE2CCV.shape[0] / self._nC) # number of components (instead of checking if cyl or not)
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
@@ -389,8 +382,8 @@ class Fields3D_e(FieldsFDEM):
:rtype: numpy.ndarray
:return: magnetic field
"""
n = int(self._aveF2CCV.shape[0] // self._nC) # Number of Components
VI = sdiag(np.kron(np.ones(n), 1. // self.prob.mesh.vol))
n = int(self._aveF2CCV.shape[0] / self._nC) # Number of Components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (self._aveF2CCV * (self._MfMui * self._b(eSolution, srcList)))
@@ -404,7 +397,7 @@ class Fields3D_e(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the derivative of the magnetic field with respect to the field we solved for with a vector
"""
n = int(self._aveF2CCV.shape[0] // self._nC) # Number of Components
n = int(self._aveF2CCV.shape[0] / self._nC) # Number of Components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
v = self._MfMui.T * (self._aveF2CCV.T * (VI.T * du_dm_v))
@@ -421,7 +414,7 @@ class Fields3D_e(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the magnetic field derivative with respect to the inversion model with a vector
"""
n = int(self._aveF2CCV.shape[0] // self._nC) # Number of Components
n = int(self._aveF2CCV.shape[0] / self._nC) # Number of Components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
v = self._MfMui.T * (self._aveF2CCV.T * (VI.T * v))
@@ -614,7 +607,7 @@ class Fields3D_b(FieldsFDEM):
:return: primary current density
"""
n = int(self._aveE2CCV.shape[0] // self._nC) # number of components
n = int(self._aveE2CCV.shape[0] / self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (self._aveE2CCV * ( self._MeSigma * self._e(bSolution,srcList ) ) )
@@ -631,7 +624,7 @@ class Fields3D_b(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the derivative of the current density with respect to the field we solved for with a vector
"""
n = int(self._aveE2CCV.shape[0] // self._nC) # number of components
n = int(self._aveE2CCV.shape[0] / self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
return self._MfMui.T * ( self._edgeCurl * ( self._aveE2CCV.T * (VI.T * du_dm_v) ) )
@@ -659,7 +652,7 @@ class Fields3D_b(FieldsFDEM):
:rtype: numpy.ndarray
:return: magnetic field
"""
n = int(self._aveF2CCV.shape[0] // self._nC) #number of components
n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (self._aveF2CCV * (self._MfMui * self._b(bSolution, srcList)))
@@ -674,7 +667,7 @@ class Fields3D_b(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the derivative of the magnetic field with respect to the field we solved for with a vector
"""
n = int(self._aveF2CCV.shape[0] // self._nC) #number of components
n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
@@ -897,7 +890,7 @@ class Fields3D_j(FieldsFDEM):
:rtype: numpy.ndarray
:return: electric field
"""
n = int(self._aveF2CCV.shape[0] // self._nC) # number of components
n = int(self._aveF2CCV.shape[0] / self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (self._aveF2CCV * (self._MfRho * self._j(jSolution, srcList)))
@@ -911,7 +904,7 @@ class Fields3D_j(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the derivative of the electric field with respect to the field we solved for with a vector
"""
n = int(self._aveF2CCV.shape[0] // self._nC) # number of components
n = int(self._aveF2CCV.shape[0] / self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
return self._MfRho.T * ( self._aveF2CCV.T * ( VI.T * du_dm_v ) )
@@ -928,7 +921,7 @@ class Fields3D_j(FieldsFDEM):
:return: product of the derivative of the electric field with respect to the model with a vector
"""
jSolution = Utils.mkvc(self[src,'jSolution'])
n = int(self._aveF2CCV.shape[0] // self._nC) # number of components
n = int(self._aveF2CCV.shape[0] / self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
return self._MfRhoDeriv(jSolution).T * ( self._aveF2CCV.T * ( VI.T * v ) )
@@ -943,7 +936,7 @@ class Fields3D_j(FieldsFDEM):
:rtype: numpy.ndarray
:return: secondary magnetic flux density
"""
n = int(self._aveE2CCV.shape[0] // self._nC) # number of components
n = int(self._aveE2CCV.shape[0] / self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (self._aveE2CCV * ( self._MeMu * self._h(jSolution,srcList)) )
@@ -958,7 +951,7 @@ class Fields3D_j(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the derivative of the magnetic flux density with respect to the field we solved for with a vector
"""
n = int(self._aveF2CCV.shape[0] // self._nC) # number of components
n = int(self._aveF2CCV.shape[0] / self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
@@ -976,7 +969,7 @@ class Fields3D_j(FieldsFDEM):
:return: product of the derivative of the magnetic flux density with respect to the model with a vector
"""
jSolution = self[src,'jSolution']
n = int(self._aveE2CCV.shape[0] // self._nC) # number of components
n = int(self._aveE2CCV.shape[0] / self._nC) # number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
s_mDeriv,_ = src.evalDeriv(self.prob, adjoint = adjoint)
@@ -1158,7 +1151,7 @@ class Fields3D_h(FieldsFDEM):
:rtype: numpy.ndarray
:return: electric field
"""
n = int(self._aveF2CCV.shape[0] // self._nC) #number of components
n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (self._aveF2CCV * (self._MfRho * self._j(hSolution, srcList)))
@@ -1172,7 +1165,7 @@ class Fields3D_h(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the derivative of the electric field with respect to the field we solved for with a vector
"""
n = int(self._aveF2CCV.shape[0] // self._nC) #number of components
n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
return self._edgeCurl.T * ( self._MfRho.T * ( self._aveF2CCV.T * ( VI.T * du_dm_v ) ) )
@@ -1189,7 +1182,7 @@ class Fields3D_h(FieldsFDEM):
:return: product of the electric field derivative with respect to the inversion model with a vector
"""
hSolution = Utils.mkvc(self[src,'hSolution'])
n = int(self._aveF2CCV.shape[0] // self._nC) #number of components
n = int(self._aveF2CCV.shape[0] / self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
return ( self._MfRhoDeriv(self._edgeCurl * hSolution).T * ( self._aveF2CCV.T * (VI.T * v) ) )
@@ -1205,7 +1198,7 @@ class Fields3D_h(FieldsFDEM):
:return: magnetic flux density
"""
h = self._h(hSolution, srcList)
n = int(self._aveE2CCV.shape[0] // self._nC) #number of components
n = int(self._aveE2CCV.shape[0] / self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
return VI * (self._aveE2CCV * (self._MeMu * h))
@@ -1220,7 +1213,7 @@ class Fields3D_h(FieldsFDEM):
:rtype: numpy.ndarray
:return: product of the derivative of the magnetic flux density with respect to the field we solved for with a vector
"""
n = int(self._aveE2CCV.shape[0] // self._nC) #number of components
n = int(self._aveE2CCV.shape[0] / self._nC) #number of components
VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol))
if adjoint:
return self._MeMu.T * (self._aveE2CCV.T * ( VI.T * du_dm_v ))
+2 -8
View File
@@ -1,13 +1,7 @@
from __future__ import absolute_import
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from future import standard_library
standard_library.install_aliases()
from SimPEG import Problem, Utils, np, sp, Solver as SimpegSolver
from scipy.constants import mu_0
from .SurveyFDEM import Survey as SurveyFDEM
from .FieldsFDEM import FieldsFDEM, Fields3D_e, Fields3D_b, Fields3D_h, Fields3D_j
from SurveyFDEM import Survey as SurveyFDEM
from FieldsFDEM import FieldsFDEM, Fields3D_e, Fields3D_b, Fields3D_h, Fields3D_j
from SimPEG.EM.Base import BaseEMProblem
from SimPEG.EM.Utils import omega
-7
View File
@@ -1,10 +1,3 @@
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
-6
View File
@@ -1,9 +1,3 @@
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 *
+2 -8
View File
@@ -1,16 +1,10 @@
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
from . import SrcFDEM as Src
from . import RxFDEM as Rx
import SrcFDEM as Src
import RxFDEM as Rx
from SimPEG import sp
class Survey(BaseEMSurvey):
+5 -11
View File
@@ -1,11 +1,5 @@
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
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
-6
View File
@@ -1,9 +1,3 @@
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):
+3 -9
View File
@@ -1,9 +1,3 @@
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
@@ -15,7 +9,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' %list(self.knownFields.keys())[0])
raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %self.knownFields.keys()[0])
if adjoint:
return self._phiDeriv_u(src, v, adjoint=adjoint), self._phiDeriv_m(src, v, adjoint=adjoint)
@@ -24,7 +18,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' %list(self.knownFields.keys())[0])
raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0])
if adjoint:
return self._eDeriv_u(src, v, adjoint), self._eDeriv_m(src, v, adjoint)
@@ -32,7 +26,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' %list(self.knownFields.keys())[0])
raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0])
if adjoint:
return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint)
+3 -9
View File
@@ -1,9 +1,3 @@
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
@@ -38,7 +32,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' %list(self.knownFields.keys())[0])
raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %self.knownFields.keys()[0])
if adjoint:
return self._phiDeriv_u(kyInd, src, v, adjoint=adjoint), self._phiDeriv_m(kyInd, src, v, adjoint=adjoint)
@@ -47,7 +41,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' %list(self.knownFields.keys())[0])
raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0])
if adjoint:
return self._eDeriv_u(kyInd, src, v, adjoint), self._eDeriv_m(kyInd, src, v, adjoint)
@@ -55,7 +49,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' %list(self.knownFields.keys())[0])
raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0])
if adjoint:
return self._jDeriv_u(kyInd, src, v, adjoint), self._jDeriv_m(kyInd, src, v, adjoint)
+3 -9
View File
@@ -1,17 +1,11 @@
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):
+3 -10
View File
@@ -1,18 +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()
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,10 +1,3 @@
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,9 +1,3 @@
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
+2 -8
View File
@@ -1,15 +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()
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,10 +1,3 @@
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):
+8 -14
View File
@@ -1,14 +1,8 @@
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
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
+1 -7
View File
@@ -1,9 +1,3 @@
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
@@ -11,7 +5,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,9 +1,3 @@
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
+2 -8
View File
@@ -1,8 +1,2 @@
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
from ProblemIP import Problem3D_CC, Problem3D_N
from SurveyIP import Survey
+2 -10
View File
@@ -1,11 +1,3 @@
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
@@ -13,7 +5,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):
"""
@@ -113,7 +105,7 @@ class BaseSIPProblem(BaseEMProblem):
JvAll = []
#Assume only eta and tau (eta first then tau)
# v = [2*Mx1]
v = v.reshape((v.size//2), 2), order='F')
v = v.reshape((int(v.size/2), 2), order='F')
for tind in range(len(self.survey.times)):
t = self.survey.times[tind]
-7
View File
@@ -1,10 +1,3 @@
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,9 +1,3 @@
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,9 +1,3 @@
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,10 +1,3 @@
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
+5 -11
View File
@@ -1,11 +1,5 @@
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
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
+19 -27
View File
@@ -1,11 +1,3 @@
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
@@ -93,10 +85,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
@@ -111,7 +103,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
@@ -192,14 +184,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(list(range(int(nstn))))*dl_x*a
stn_y = endl[0,1] + np.array(list(range(int(nstn))))*dl_y*a
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
if mesh.dim==2:
ztop = mesh.vectorNy[-1]
@@ -238,15 +230,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([(AB - b) // a, n])
nstn = np.min([np.floor( (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(list(range(int(nstn))))*dl_x*a
stn_y = N[ii,1] + dl_y*b + np.array(list(range(int(nstn))))*dl_y*a
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
# Create receiver poles
@@ -283,17 +275,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(list(range(int(nstn))))*dl_x*a
stn_y = min_y + np.array(list(range(int(nstn))))*dl_y*a
stn_x = min_x + np.array(range(int(nstn)))*dl_x*a
stn_y = min_y + np.array(range(int(nstn)))*dl_y*a
# Define number of cross lines
nlin = int(box_w // a)
lind = list(range(-nlin,nlin+1))
nlin = int(np.floor( box_w / a ))
lind = range(-nlin,nlin+1)
ngrad = nstn * len(lind)
@@ -318,7 +310,7 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
srcClass = DC.Src.Dipole([rxClass], M[0,:], N[-1,:])
SrcList.append(srcClass)
else:
print("""stype must be either 'pdp', 'dpdp' or 'gradient'. """)
print """stype must be either 'pdp', 'dpdp' or 'gradient'. """
return SrcList
+1 -7
View File
@@ -1,7 +1 @@
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 *
from StaticUtils import *
+3 -9
View File
@@ -1,9 +1,3 @@
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
import DC
import IP
import SIP
+14 -21
View File
@@ -1,10 +1,3 @@
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 *
@@ -54,7 +47,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)
@@ -62,7 +55,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):
@@ -77,13 +70,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
@@ -102,13 +95,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
@@ -130,14 +123,14 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
* Compute \\\(\\\\vec{w} = -\\\mathbf{Q} \\\\vec{y}\\\)
"""
if self.verbose: print('%s\nCalculating J(v)\n%s'%('*'*50,'*'*50))
if self.verbose: print '%s\nCalculating J(v)\n%s'%('*'*50,'*'*50)
self.curModel = m
if f is None:
f = self.fields(m)
p = self.Gvec(m, v, f)
y = self.solveAh(m, p)
Jv = self.survey.evalDeriv(f, v=y)
if self.verbose: print('%s\nDone calculating J(v)\n%s'%('*'*50,'*'*50))
if self.verbose: print '%s\nDone calculating J(v)\n%s'%('*'*50,'*'*50)
return - mkvc(Jv)
def Jtvec(self, m, v, f=None):
@@ -155,7 +148,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)
@@ -166,6 +159,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)
+5 -11
View File
@@ -1,13 +1,7 @@
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):
@@ -93,7 +87,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')
@@ -105,7 +99,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")
@@ -133,7 +127,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)
@@ -145,7 +139,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")
+2 -9
View File
@@ -1,14 +1,7 @@
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 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):
+3 -9
View File
@@ -1,9 +1,3 @@
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
from SurveyTDEM import * #SurveyTDEM, RxTDEM, SrcTDEM
from BaseTDEM import BaseTDEMProblem, FieldsTDEM
from TDEM_b import ProblemTDEM_b
+5 -12
View File
@@ -1,10 +1,3 @@
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
@@ -24,7 +17,7 @@ def MagneticDipoleVectorPotential(srcLoc, obsLoc, component, moment=1., dipoleMo
#TODO: break this out!
if type(component) in [list, tuple]:
out = list(range(len(component)))
out = range(len(component))
for i, comp in enumerate(component):
out[i] = MagneticDipoleVectorPotential(srcLoc, obsLoc, comp, dipoleMoment=dipoleMoment)
return np.concatenate(out)
@@ -125,7 +118,7 @@ def MagneticLoopVectorPotential(srcLoc, obsLoc, component, radius, mu=mu_0):
"""
if type(component) in [list, tuple]:
out = list(range(len(component)))
out = range(len(component))
for i, comp in enumerate(component):
out[i] = MagneticLoopVectorPotential(srcLoc, obsLoc, comp, radius, mu)
return np.concatenate(out)
@@ -165,11 +158,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')
+2 -8
View File
@@ -1,9 +1,3 @@
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
@@ -15,8 +9,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
+2 -8
View File
@@ -1,8 +1,2 @@
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
from EMUtils import omega, k
from AnalyticUtils import MagneticDipoleFields, MagneticDipoleVectorPotential, MagneticLoopVectorPotential
+6 -13
View File
@@ -1,10 +1,3 @@
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
@@ -65,7 +58,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)
@@ -90,7 +83,7 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False):
try:
from pymatsolver import MumpsSolver
prb.Solver = MumpsSolver
except ImportError as e:
except ImportError, e:
prb.Solver = SolverLU
return prb
@@ -101,7 +94,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
@@ -119,7 +112,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)
@@ -128,11 +121,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
+6 -12
View File
@@ -1,13 +1,7 @@
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
import TDEM
import FDEM
import Static
import Base
import Analytics
import Utils
from scipy.constants import mu_0, epsilon_0
+3 -9
View File
@@ -1,9 +1,3 @@
from __future__ import print_function
from __future__ import division
from __future__ import unicode_literals
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import *
import SimPEG.EM.Static.DC as DC
@@ -35,7 +29,7 @@ def run(plotIt=True):
try:
from pymatsolver import MumpsSolver
problem.Solver = MumpsSolver
except Exception as e:
except Exception, e:
pass
data = survey.dpred(sigma)
@@ -67,8 +61,8 @@ def run(plotIt=True):
ax[0].set_title('Computed')
plt.show()
return np.linalg.norm(data-data_ana) / np.linalg.norm(data_ana)
return np.linalg.norm(data-data_ana)/np.linalg.norm(data_ana)
if __name__ == '__main__':
print(run())
print run()
+6 -14
View File
@@ -1,11 +1,3 @@
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
@@ -65,7 +57,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
@@ -132,10 +124,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)
@@ -151,10 +143,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), end=' ')
print '\rTransmitter {0} of {1} -> Time:{2} sec'.format(ii,len(Tx),time.time()- start_time),
print('Transmitter {0} of {1}'.format(ii,len(Tx)))
print('Forward completed')
print 'Transmitter {0} of {1}'.format(ii,len(Tx))
print 'Forward completed'
# Let's just convert the 3D format into 2D (distance along line) and plot
survey2D = DC.convertObs_DC3D_to_2D(survey, np.ones(survey.nSrc) , 'Xloc')
+1 -7
View File
@@ -1,9 +1,3 @@
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
@@ -62,7 +56,7 @@ def run(plotIt=True):
try:
from pymatsolver import MumpsSolver
prb.Solver = MumpsSolver
except ImportError as e:
except ImportError, e:
prb.Solver = SolverLU
prb.pair(survey)
@@ -1,9 +1,3 @@
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
+14 -21
View File
@@ -1,10 +1,3 @@
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
@@ -74,8 +67,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.]
@@ -85,25 +78,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.log(csx2/csx1) // np.log(pfx1)
npadx1 = np.floor(np.log(csx2/csx1) / np.log(pfx1))
hx1a,hx1b = Utils.meshTensor([(csx1,ncx1)]),Utils.meshTensor([(csx1,npadx1,pfx1)])
dx1 = sum(hx1a)+sum(hx1b)
dx1 = dx1 // csx2
hx1b *= (dx1*csx2 - sum(hx1a)) / sum(hx1b)
dx1 = np.floor(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])
@@ -117,8 +110,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))
@@ -231,7 +224,7 @@ def run(plotIt=True):
# ------------- Solve ---------------------------
t0 = time.time()
fieldsCasing = problem.fields(sigCasing)
print('Time to solve 2 sources', time.time() - t0)
print 'Time to solve 2 sources', time.time() - t0
# Plot current
@@ -258,9 +251,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,9 +1,3 @@
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,9 +1,3 @@
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
+5 -13
View File
@@ -1,11 +1,3 @@
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 *
@@ -55,7 +47,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
@@ -67,15 +59,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])
@@ -83,7 +75,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,10 +1,3 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from builtins import range
from SimPEG import *
@@ -1,9 +1,3 @@
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
-6
View File
@@ -1,9 +1,3 @@
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
-7
View File
@@ -1,10 +1,3 @@
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):
-6
View File
@@ -1,9 +1,3 @@
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):
-6
View File
@@ -1,9 +1,3 @@
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,9 +1,3 @@
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):
+5 -10
View File
@@ -1,10 +1,5 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from SimPEG import *
from SimPEG import Mesh, Utils, np
def run(plotIt=True):
"""
@@ -14,15 +9,15 @@ def run(plotIt=True):
Here we show SimPEG used to create three different types of meshes.
"""
sz = [16,16]
sz = [16, 16]
tM = Mesh.TensorMesh(sz)
qM = Mesh.TreeMesh(sz)
qM.refine(lambda cell: 4 if np.sqrt(((np.r_[cell.center]-0.5)**2).sum()) < 0.4 else 3)
rM = Mesh.CurvilinearMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate'))
rM = Mesh.CurvilinearMesh(Utils.meshutils.exampleLrmGrid(sz, 'rotate'))
if plotIt:
import matplotlib.pyplot as plt
fig, axes = plt.subplots(1,3,figsize=(14,5))
fig, axes = plt.subplots(1, 3, figsize=(14, 5))
opts = {}
tM.plotGrid(ax=axes[0], **opts)
axes[0].set_title('TensorMesh')
@@ -1,11 +1,3 @@
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):
@@ -95,7 +87,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,10 +1,3 @@
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):
+3 -11
View File
@@ -1,11 +1,3 @@
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):
@@ -36,15 +28,15 @@ def run(plotIt=True, n=60):
axes[0].set_xlim([-1,17])
axes[0].set_ylim([-1,17])
for ii, loc in zip(list(range(M.nC)),M.gridCC):
for ii, loc in zip(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(list(range(M.nFx)),M.gridFx):
for ii, loc in zip(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(list(range(M.nFy)),M.gridFy):
for ii, loc in zip(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,9 +1,3 @@
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,9 +1,3 @@
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):
+65
View File
@@ -0,0 +1,65 @@
from SimPEG import Mesh, np, PF
def run(plotIt=True):
"""
PF: Magnetics: Analytics
========================
Comparing the magnetics field in Vancouver to Seoul
"""
xr = np.linspace(-300, 300, 41)
yr = np.linspace(-300, 300, 41)
X, Y = np.meshgrid(xr, yr)
Z = np.ones((np.size(xr), np.size(yr)))*150
# Bz component in Korea
inckr = -8. + 3./60
deckr = 54. + 9./60
btotkr = 50898.6
Bokr = PF.MagAnalytics.IDTtoxyz(inckr, deckr, btotkr)
bx, by, bz = PF.MagAnalytics.MagSphereAnaFunA(
X, Y, Z, 100., 0., 0., 0., 0.01, Bokr, 'secondary'
)
Bzkr = np.reshape(bz, (np.size(xr), np.size(yr)), order='F')
# Bz component in Canada
incca = 16. + 49./60
decca = 70. + 19./60
btotca = 54692.1
Boca = PF.MagAnalytics.IDTtoxyz(incca, decca, btotca)
bx, by, bz = PF.MagAnalytics.MagSphereAnaFunA(
X, Y, Z, 100., 0., 0., 0., 0.01, Boca, 'secondary'
)
Bzca = np.reshape(bz, (np.size(xr), np.size(yr)), order='F')
if plotIt:
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
fig = plt.figure(figsize=(14, 5))
ax1 = plt.subplot(121)
dat1 = plt.imshow(Bzkr, extent=[min(xr), max(xr), min(yr), max(yr)])
divider = make_axes_locatable(ax1)
cax1 = divider.append_axes("right", size="5%", pad=0.05)
ax1.set_xlabel('East-West (m)')
ax1.set_ylabel('South-North (m)')
plt.colorbar(dat1, cax=cax1)
ax1.set_title('$B_z$ field at Seoul, South Korea')
ax2 = plt.subplot(122)
dat2 = plt.imshow(Bzca, extent=[min(xr), max(xr), min(yr), max(yr)])
divider = make_axes_locatable(ax2)
cax2 = divider.append_axes("right", size="5%", pad=0.05)
ax2.set_xlabel('East-West (m)')
ax2.set_ylabel('South-North (m)')
plt.colorbar(dat2, cax=cax2)
ax2.set_title('$B_z$ field at Vancouver, Canada')
plt.show()
if __name__ == '__main__':
run()
@@ -1,9 +1,3 @@
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
+23 -30
View File
@@ -1,35 +1,28 @@
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 #####
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
import DC_Analytic_Dipole
import DC_Forward_PseudoSection
import EM_FDEM_1D_Inversion
import EM_FDEM_Analytic_MagDipoleWholespace
import EM_Schenkel_Morrison_Casing
import EM_TDEM_1D_Inversion
import FLOW_Richards_1D_Celia1990
import Inversion_IRLS
import Inversion_Linear
import Maps_ComboMaps
import Maps_Mesh2Mesh
import Mesh_Basic_ForwardDC
import Mesh_Basic_PlotImage
import Mesh_Basic_Types
import Mesh_Operators_CahnHilliard
import Mesh_QuadTree_Creation
import Mesh_QuadTree_FaceDiv
import Mesh_QuadTree_HangingNodes
import Mesh_Tensor_Creation
import MT_1D_ForwardAndInversion
import MT_3D_Foward
import Utils_surface2ind_topo
__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"]
@@ -111,7 +104,7 @@ if __name__ == '__main__':
rst = os.path.sep.join((filePath.split(os.path.sep)[:-3] + ['docs', 'content', 'examples', name + '.rst']))
print('Creating: %s.rst'%name)
print 'Creating: %s.rst'%name
f = open(rst, 'w')
f.write(out)
f.close()
+6 -12
View File
@@ -1,20 +1,14 @@
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(with_metaclass(Utils.SimPEGMetaClass, object)):
class NonLinearMap(object):
"""
SimPEG NonLinearMap
"""
__metaclass__ = Utils.SimPEGMetaClass
counter = None #: A SimPEG.Utils.Counter object
mesh = None #: A SimPEG Mesh
@@ -192,7 +186,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
@@ -279,7 +273,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
@@ -349,7 +343,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
+8 -15
View File
@@ -1,12 +1,5 @@
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
@@ -68,7 +61,7 @@ class RichardsSurvey(Survey.BaseSurvey):
@Utils.requires('prob')
def eval(self, U, m):
Ds = list(range(len(self.rxList)))
Ds = range(len(self.rxList))
for ii, rx in enumerate(self.rxList):
Ds[ii] = rx.eval(U, m,
self.prob.mapping,
@@ -80,7 +73,7 @@ class RichardsSurvey(Survey.BaseSurvey):
@Utils.requires('prob')
def evalDeriv(self, U, m):
"""The Derivative with respect to the fields."""
Ds = list(range(len(self.rxList)))
Ds = range(len(self.rxList))
for ii, rx in enumerate(self.rxList):
Ds[ii] = rx.evalDeriv(U, m,
self.prob.mapping,
@@ -142,12 +135,12 @@ class RichardsProblem(Problem.BaseTimeProblem):
@Utils.timeIt
def fields(self, m):
tic = time.time()
u = list(range(self.nT+1))
u = 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
@@ -245,7 +238,7 @@ class RichardsProblem(Problem.BaseTimeProblem):
f = self.fields(m)
nn = len(f)-1
Asubs, Adiags, Bs = list(range(nn)), list(range(nn)), list(range(nn))
Asubs, Adiags, Bs = range(nn), range(nn), range(nn)
for ii in range(nn):
dt = self.timeSteps[ii]
bc = self.getBoundaryConditions(ii, f[ii])
@@ -270,7 +263,7 @@ class RichardsProblem(Problem.BaseTimeProblem):
if f is None:
f = self.fields(m)
JvC = list(range(len(f)-1)) # Cell to hold each row of the long vector.
JvC = 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])
@@ -302,7 +295,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 = list(range((ii)*Adiag.shape[0], (ii+1)*Adiag.shape[0]))
vpart = 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.
+2 -8
View File
@@ -1,8 +1,2 @@
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 *
import Empirical
from RichardsProblem import *
+1 -7
View File
@@ -1,7 +1 @@
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
import Richards
+2 -11
View File
@@ -1,13 +1,4 @@
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
import Utils, numpy as np, scipy.sparse as sp
class Fields(object):
"""Fancy Field Storage
@@ -253,7 +244,7 @@ class TimeFields(Fields):
out = func(pointerFields, srcII, timeII)
else: #loop over the time steps
nT = pointerShape[2]
out = list(range(nT))
out = range(nT)
for i, TIND_i in enumerate(timeII):
fieldI = pointerFields[:,:,i]
if fieldI.shape[0] == fieldI.size:
+12 -19
View File
@@ -1,21 +1,14 @@
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
import Utils, Survey, Problem, numpy as np, scipy.sparse as sp, gc
from Utils.SolverUtils import *
import DataMisfit
import Regularization
class BaseInvProblem(with_metaclass(Utils.SimPEGMetaClass, object)):
class BaseInvProblem(object):
"""BaseInvProblem(dmisfit, reg, opt)"""
__metaclass__ = Utils.SimPEGMetaClass
beta = 1.0 #: Trade-off parameter
debug = False #: Print debugging information
@@ -61,10 +54,10 @@ class BaseInvProblem(with_metaclass(Utils.SimPEGMetaClass, 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
@@ -72,8 +65,8 @@ class BaseInvProblem(with_metaclass(Utils.SimPEGMetaClass, 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
@@ -94,7 +87,7 @@ class BaseInvProblem(with_metaclass(Utils.SimPEGMetaClass, 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:
+5 -11
View File
@@ -1,24 +1,18 @@
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
from . import Directives
from future.utils import with_metaclass
from Optimization import Remember, IterationPrinters, StoppingCriteria
import Directives
class BaseInversion(with_metaclass(Utils.SimPEGMetaClass, object)):
class BaseInversion(object):
"""
Inversion Class.
"""
__metaclass__ = Utils.SimPEGMetaClass
name = 'BaseInversion'
debug = False #: Print debugging information
+2 -8
View File
@@ -1,13 +1,7 @@
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):
+6 -12
View File
@@ -1,9 +1,3 @@
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
@@ -69,7 +63,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:
@@ -194,7 +188,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):
'''
@@ -204,7 +198,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
@@ -233,7 +227,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:
@@ -244,7 +238,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:
@@ -354,4 +348,4 @@ class Fields3D_e(BaseMTFields):
This function stacks the fields derivatives appropriately
"""
# The fields have no dependance to the model.
return None
return None
+7 -13
View File
@@ -1,9 +1,3 @@
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
@@ -52,7 +46,7 @@ class eForm_psField(BaseMTProblem):
Edge inner product matrix
"""
if getattr(self, '_MeMui', None) is None:
self._MeMui = self.mesh.getEdgeInnerProduct(old_div(1.0,mu_0))
self._MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
return self._MeMui
@property
@@ -148,7 +142,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)
@@ -164,7 +158,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
@@ -197,7 +191,7 @@ class eForm_TotalField(BaseMTProblem):
Edge inner product matrix
"""
if getattr(self, '_MeMui', None) is None:
self._MeMui = self.mesh.getEdgeInnerProduct(old_div(1.0,mu_0))
self._MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
return self._MeMui
@property
@@ -255,7 +249,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 = ((old_div(Etot,Etot[-1]))*sourceAmp) # Scale the fields to be equal to sourceAmp at the top
Etot = ((Etot/Etot[-1])*sourceAmp) # Scale the fields to be equal to sourceAmp at the top
## Note: The analytic solution is derived with e^iwt
eBC = np.r_[Etot[0],Etot[-1]]
# The right hand side
@@ -280,7 +274,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)
@@ -292,6 +286,6 @@ class eForm_TotalField(BaseMTProblem):
# NOTE: only store e fields
F[Src, 'e_1dSolution'] = e[:,0]
if self.verbose:
print('Ran for {:f} seconds'.format(time.time()-startTime))
print 'Ran for {:f} seconds'.format(time.time()-startTime)
sys.stdout.flush()
return F
+1 -7
View File
@@ -1,7 +1 @@
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
from Probs import eForm_TotalField, eForm_psField
-6
View File
@@ -1,7 +1 @@
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
+2 -8
View File
@@ -1,9 +1,3 @@
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
@@ -121,7 +115,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)
@@ -137,7 +131,7 @@ class eForm_ps(BaseMTProblem):
# Note curl e = -iwb so b = -curl/iw
if self.verbose:
print('Ran for {:f} seconds'.format(time.time()-startTime))
print 'Ran for {:f} seconds'.format(time.time()-startTime)
sys.stdout.flush()
Ainv.clean()
return F
+1 -7
View File
@@ -1,7 +1 @@
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
from Probs import eForm_ps
+4 -10
View File
@@ -1,16 +1,10 @@
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
#################
@@ -84,7 +78,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):
@@ -161,7 +155,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):
+21 -27
View File
@@ -1,16 +1,10 @@
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
from . import SrcMT
from Utils import rec2ndarr
import SrcMT
import sys
#################
@@ -82,7 +76,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 = old_div(-ex,bx)
f_part_complex = -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.
@@ -109,13 +103,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 = old_div(( ex_px*hy_py - ex_py*hy_px),(hx_px*hy_py - hx_py*hy_px))
f_part_complex = ( ex_px*hy_py - ex_py*hy_px)/(hx_px*hy_py - hx_py*hy_px)
elif 'zxy' in self.rxType:
f_part_complex = old_div((-ex_px*hx_py + ex_py*hx_px),(hx_px*hy_py - hx_py*hy_px))
f_part_complex = (-ex_px*hx_py + ex_py*hx_px)/(hx_px*hy_py - hx_py*hy_px)
elif 'zyx' in self.rxType:
f_part_complex = old_div(( ey_px*hy_py - ey_py*hy_px),(hx_px*hy_py - hx_py*hy_px))
f_part_complex = ( ey_px*hy_py - ey_py*hy_px)/(hx_px*hy_py - hx_py*hy_px)
elif 'zyy' in self.rxType:
f_part_complex = old_div((-ey_px*hx_py + ey_py*hx_px),(hx_px*hy_py - hx_py*hy_px))
f_part_complex = (-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]
@@ -133,9 +127,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 = old_div((- by_px*bz_py + by_py*bz_px),(bx_px*by_py - bx_py*by_px))
f_part_complex = (- by_px*bz_py + by_py*bz_px)/(bx_px*by_py - bx_py*by_px)
if 'tzy' in self.rxType:
f_part_complex = old_div(( bx_px*bz_py - bx_py*bz_px),(bx_px*by_py - bx_py*by_px))
f_part_complex = ( 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))
@@ -163,8 +157,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(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)
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)
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')
@@ -204,7 +198,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(old_div(1.,(sDiag(hx_px)*hy_py - sDiag(hx_py)*hy_px)))
Hd = sDiag(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:
@@ -253,7 +247,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(old_div(1.,(sDiag(bx_px)*by_py - sDiag(bx_py)*by_px)))
Hd = sDiag(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 ))
@@ -273,8 +267,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(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_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_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':
@@ -306,17 +300,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: 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)
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
# 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(old_div(1.,(sDiag(ahx_px)*ahy_py - sDiag(ahx_py)*ahy_px)))
aHd = sDiag(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:
@@ -368,7 +362,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(old_div(1.,(sDiag(abx_px)*aby_py - sDiag(abx_py)*aby_px)))
aHd = sDiag(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:
+10 -17
View File
@@ -1,10 +1,3 @@
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
@@ -40,8 +33,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 = 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
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
# Build the propagation matrix
# Convert fields to down/up going components in layer below current layer
@@ -55,7 +48,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] = old_div(UDp[:,lnr+1::-1],UDp[1,lnr+1])
UDp[:,lnr+1::-1] = UDp[:,lnr+1::-1]/UDp[1,lnr+1]
# Calculate the fields
Ed = np.empty((zd.size,),dtype=complex)
@@ -69,14 +62,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] = (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]))
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]))
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] = (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]))
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]))
# Return return the fields
return Ed, Eu, Hd, Hu
@@ -99,15 +92,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] = old_div((mu_0*om),np.sqrt(mu_0*eps_0*(om)**2 - 1j*mu_0*sigma[0]*om))
Zall[0] = (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 = old_div((mu_0*om),k)
Z = (mu_0*om)/k
Zall[nr+1] = Z *(old_div((Zall[nr] + Z*np.tanh(1j*k*hi)),(Z + Zall[nr]*np.tanh(1j*k*hi))))
Zall[nr+1] = Z *((Zall[nr] + Z*np.tanh(1j*k*hi))/(Z + Zall[nr]*np.tanh(1j*k*hi)))
#pdb.set_trace()
Z1d[nrFr] = Zall[-1]
+3 -9
View File
@@ -1,11 +1,5 @@
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):
@@ -15,7 +9,7 @@ def get1DEfields(m1d,sigma,freq,sourceAmp=1.0):
G = m1d.nodalGrad
# Mass matrices
# Magnetic permeability
Mmu = simpeg.Utils.sdiag(m1d.vol*(old_div(1.0,mu_0)))
Mmu = simpeg.Utils.sdiag(m1d.vol*(1.0/mu_0))
# Conductivity
Msig = m1d.getFaceInnerProduct(sigma)
# Set up the solution matrix
@@ -29,7 +23,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 = ((old_div(Etot,Etot[-1]))*sourceAmp) # Scale the fields to be equal to sourceAmp at the top
Etot = ((Etot/Etot[-1])*sourceAmp) # Scale the fields to be equal to sourceAmp at the top
## Note: The analytic solution is derived with e^iwt
bc = np.r_[Etot[0],Etot[-1]]
# The right hand side
+4 -10
View File
@@ -1,10 +1,4 @@
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 *
from MT1Dsolutions import * # Add the names of the functions
from MT1Danalytic import *
from dataUtils import *
from ediFilesUtils import *
+9 -15
View File
@@ -1,9 +1,3 @@
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
@@ -51,13 +45,13 @@ def rotateData(MTdata, rotAngle):
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))
app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
app_phs = np.arctan2(z.imag,z.real)*(180/np.pi)
return app_res, app_phs
def skindepth(rho, freq):
''' Function to calculate the skindepth of EM waves'''
return np.sqrt( (rho*((old_div(1,(freq * mu_0 * np.pi ))))))
return np.sqrt( (rho*((1/(freq * mu_0 * np.pi )))))
def rec2ndarr(x, dt=float):
return x.view((dt, len(x.dtype.names)))
@@ -70,7 +64,7 @@ def makeAnalyticSolution(mesh, model, elev, freqs):
anaE = anaEd+anaEu
anaH = anaHd+anaHu
anaZ = old_div(anaE,anaH)
anaZ = 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)])
@@ -103,7 +97,7 @@ def plotMT1DModelData(problem, models, symList=None):
# if not symList:
# symList = ['x']*len(models)
from . import plotDataTypes as pDt
import plotDataTypes as pDt
# Loop through the models.
modelList = [problem.survey.mtrue]
modelList.extend(models)
@@ -116,14 +110,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 = old_div(nr,((len(modelList)-1.999)*1.))
colRat = 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(old_div(1.,(problem.mapping.sigmaMap*model)),np.ones(2,))
modelPts = np.kron(1./(problem.mapping.sigmaMap*model),np.ones(2,))
axM.semilogx(modelPts,meshPts,color=col)
## Data
@@ -150,7 +144,7 @@ def plotMT1DModelData(problem, models, symList=None):
# Fix labels and ticks
yMtick = [old_div(l,1000) for l in axM.get_yticks().tolist()]
yMtick = [l/1000 for l in axM.get_yticks().tolist()]
axM.set_yticklabels(yMtick)
[ l.set_rotation(90) for l in axM.get_yticklabels()]
[ l.set_rotation(90) for l in axR.get_yticklabels()]
@@ -163,7 +157,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
+4 -13
View File
@@ -1,12 +1,3 @@
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
@@ -18,7 +9,7 @@ import numpy as np
import os, sys, re
class EDIimporter(object):
class EDIimporter:
"""
A class to import EDIfiles.
@@ -27,7 +18,7 @@ class EDIimporter(object):
# Define data converters
_impUnitEDI2SI = 4*np.pi*1e-4 # Convert Z[mV/km/nT] (as in EDI)to Z[V/A] SI unit
_impUnitSI2EDI = old_div(1.,_impUnitEDI2SI) # ConvertZ[V/A] SI unit to Z[mV/km/nT] (as in EDI)
_impUnitSI2EDI = 1./_impUnitEDI2SI # ConvertZ[V/A] SI unit to Z[mV/km/nT] (as in EDI)
# Properties
filesList = None
@@ -125,7 +116,7 @@ class EDIimporter(object):
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:
@@ -135,7 +126,7 @@ class EDIimporter(object):
if self._outEPSG is None:
# Find the UTM EPSG number
Nnr = 700 if latD < 0.0 else 600
utmZ = int(1+old_div((longD+180.0),6.0))
utmZ = int(1+(longD+180.0)/6.0)
self._outEPSG = 32000 + Nnr + utmZ
out.ImportFromEPSG(self._outEPSG)
self._2out = osr.CoordinateTransformation(src,out)
+49 -55
View File
@@ -1,9 +1,3 @@
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
@@ -41,7 +35,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)))
@@ -71,45 +65,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 = 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])))
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]))
else:
zPlot = old_div((np.abs(arrayList[0][flag][indUniFreq0]) - np.abs(arrayList[1][flag][indUniFreq1])),np.abs(arrayList[1][flag][indUniFreq1]))
zPlot = (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 = 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])))))
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]))))
else:
zPlot = old_div((np.real(arrayList[0][flag][indUniFreq0]) -np.real(arrayList[1][flag][indUniFreq1])),np.abs((np.real(arrayList[1][flag][indUniFreq1]))))
zPlot = (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 = 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])))))
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]))))
else:
zPlot = old_div((np.imag(arrayList[0][flag][indUniFreq0]) -np.imag(arrayList[1][flag][indUniFreq1])),np.abs((np.imag(arrayList[1][flag][indUniFreq1]))))
zPlot = (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)
@@ -180,13 +174,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 = old_div(1,(8*np.pi**2*10**(-7)))
def plotIsoStaImpedance(ax,loc,array,flag,par='abs',pSym='s',pColor=None):
appResFact = 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':
@@ -194,9 +188,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 = (old_div(appResFact,freq))*np.abs(array[flag][indUniSta])**2
zPlot = (appResFact/freq)*np.abs(array[flag][indUniSta])**2
elif par == 'phs':
zPlot = np.arctan2(array[flag][indUniSta].imag,array[flag][indUniSta].real)*(old_div(180,np.pi))
zPlot = np.arctan2(array[flag][indUniSta].imag,array[flag][indUniSta].real)*(180/np.pi)
if not pColor:
if 'xx' in flag:
@@ -204,7 +198,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'
@@ -213,14 +207,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(list(sectDict.values())[0]==array[list(sectDict.keys())[0]])
indSect = np.where(sectDict.values()[0]==array[sectDict.keys()[0]])
# Define the plot axes
if 'x' in list(sectDict.keys())[0]:
if 'x' in sectDict.keys()[0]:
x = array['y'][indSect]
else:
x = array['x'][indSect]
@@ -235,9 +229,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 = old_div(np.abs(array[flag][indSect])**2,(8*np.pi**2*10**(-7)*array['freq'][indSect]))
elif par == 'ares':
zPlot = 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])
@@ -250,7 +244,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)*(old_div(180,np.pi))
zPlot = np.arctan2(array[flag][indSect].imag,array[flag][indSect].real)*(180/np.pi)
cmap = plt.get_cmap('RdYlBu')#seismic)
if cLevel:
zMax = cLevel[1]
@@ -277,7 +271,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])
@@ -287,7 +281,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':
@@ -313,14 +307,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(list(sectDict.values())[0]==arrayList[0][list(sectDict.keys())[0]])
indSect1 = np.where(list(sectDict.values())[0]==arrayList[1][list(sectDict.keys())[0]])
indSect0 = np.where(sectDict.values()[0]==arrayList[0][sectDict.keys()[0]])
indSect1 = np.where(sectDict.values()[0]==arrayList[1][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 list(sectDict.keys())[0]:
if 'x' in sectDict.keys()[0]:
x0 = arr0['y']
x1 = arr1['y']
else:
@@ -328,35 +322,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 = old_div((np.log10(np.abs(arr0[flag])) - np.log10(np.abs(arr1[flag]))),np.log10(np.abs(arr1[flag])))
zPlot = (np.log10(np.abs(arr0[flag])) - np.log10(np.abs(arr1[flag])))/np.log10(np.abs(arr1[flag]))
else:
zPlot = old_div((np.abs(arr0[flag]) - np.abs(arr1[flag])),np.abs(arr1[flag]))
zPlot = (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 = old_div(1,(8*np.pi**2*10**(-7)))
arF = 1/(8*np.pi**2*10**(-7))
if useLog:
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))
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)
else:
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))
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)
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 = 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))))
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))
else:
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))))
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))
if mask:
maskInd = np.logical_or(np.abs(arr0[flag])< 1e-3,np.abs(arr1[flag]) < 1e-3)
zPlot = np.ma.array(zPlot)
@@ -364,26 +358,26 @@ def plotPsudoSectNSDiff(ax,sectDict,arrayList,flag,par='abs',colorbar=True,color
cmap = plt.get_cmap('Spectral')#seismic)
elif par == 'real':
if useLog:
zPlot = old_div((np.log10(arr0[flag].real) - np.log10(arr1[flag].real)),np.log10(arr1[flag].real))
zPlot = (np.log10(arr0[flag].real) - np.log10(arr1[flag].real))/np.log10(arr1[flag].real)
else:
zPlot = old_div((arr0[flag].real - arr1[flag].real),arr1[flag].real)
zPlot = (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 = old_div((np.log10(arr0[flag].imag) - np.log10(arr1[flag].imag)),np.log10(arr1[flag].imag))
zPlot = (np.log10(arr0[flag].imag) - np.log10(arr1[flag].imag))/np.log10(arr1[flag].imag)
else:
zPlot = old_div((arr0[flag].imag - arr1[flag].imag),arr1[flag].imag)
zPlot = (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])
@@ -391,18 +385,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,old_div((cLevel[1] - cLevel[0]),50.))
clevel = np.arange(cLevel[0],cLevel[1]+.1,old_div((cLevel[1] - cLevel[0]),10.))
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.)
else:
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.))
level = np.arange(zPlot.min(),zPlot.max(),(zPlot.max() - zPlot.min())/50.)
clevel = np.arange(zPlot.min(),zPlot.max(),(zPlot.max() - zPlot.min())/10.)
plotNorm = colors.Normalize()
elif colorNorm=='Log':
level = np.logspace(zMin-.125,zMax,(zMax-zMin)*8+1,endpoint=True)
+1 -9
View File
@@ -1,11 +1,3 @@
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from builtins import dict
from future import standard_library
standard_library.install_aliases()
from builtins import zip
import SimPEG as simpeg, numpy as np
def homo1DModelSource(mesh,freq,sigma_1d):
@@ -98,7 +90,7 @@ def analytic1DModelSource(mesh,freq,sigma_1d):
Eu, Ed, _, _ = getEHfields(mesh1d,sigma_1d,freq,mesh.vectorNz)
# Make the fields into a dictionary of location and the fields
e0_1d = Eu+Ed
E1dFieldDict = dict(list(zip(mesh.vectorNz,e0_1d)))
E1dFieldDict = dict(zip(mesh.vectorNz,e0_1d))
if mesh.dim == 1:
eBG_px = simpeg.mkvc(e0_1d,2)
eBG_py = -simpeg.mkvc(e0_1d,2) # added a minus to make the results in the correct quadrents.

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