wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-19 11:28:00 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #4 from simpeg/Tx2Src

Browse Source 
Tx2 src
This commit is contained in:
Rowan Cockett
2015-05-15 13:07:25 -07:00
parent 59c2216ec0 bcc3394079
commit 1455901474
14 changed files with 657 additions and 213 deletions
Download Patch File Download Diff File Expand all files Collapse all files
LICENSE
+1 -1
View File
@@ -1,6 +1,6 @@
The MIT License (MIT)
Copyright (c) 2013-2014 SimPEG Developers
Copyright (c) 2013-2015 SimPEG Developers
Permission is hereby granted, free of charge, to any person obtaining a copy of
this software and associated documentation files (the "Software"), to deal in
notebooks/DC_schumberger_FWD.ipynb
+45 -35
View File
File diff suppressed because one or more lines are too long
notebooks/DC_schumberger_Inv.ipynb
+121 -138
View File
File diff suppressed because one or more lines are too long
notebooks/DCexample.ipynb
+194 -30
View File
File diff suppressed because one or more lines are too long
notebooks/DCfwd.png
BIN
View File
Binary file not shown.
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 105 KiB

After

Width:  |  Height:  |  Size: 115 KiB

notebooks/DCsurvey.png
BIN
View File
Binary file not shown.
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 57 KiB

After

Width:  |  Height:  |  Size: 55 KiB

notebooks/ThreeMesh.png
BIN
View File
Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 206 KiB

notebooks/ThreeMeshDC.ipynb
+278
View File
@@ -0,0 +1,278 @@
{
"metadata": {
"name": "",
"signature": "sha256:4ee97c21ba8374164213c87af95095c4694b605af46490529f6e0f5ff4051dea"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "code",
"collapsed": false,
"input": [
"from SimPEG import Mesh, Utils, np, SolverLU\n",
"import matplotlib.pyplot as plt\n",
"# %pylab inline"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 5
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import matplotlib\n",
"from matplotlib.mlab import griddata"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 17
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"matplotlib.rcParams.update({'font.size': 16, 'text.usetex': True})"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 18
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"sz = [15,15]\n",
"tM = Mesh.TensorMesh(sz)\n",
"qM = Mesh.TreeMesh(sz)\n",
"qM.refine(lambda X: 1 if np.sqrt(((X-0.5)**2).sum()) < 0.45 else 0)\n",
"rM = Mesh.LogicallyRectMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate'))"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 57
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"fun = lambda X: 1 if np.sqrt(((X-0.5)**2).sum()) < 0.45 else 0"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 58
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"fun(np.ones(30))"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 59,
"text": [
"0"
]
}
],
"prompt_number": 59
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def DCfun(mesh, pts):\n",
" D = mesh.faceDiv\n",
" G = D.T\n",
" sigma = 1e-2*np.ones(mesh.nC)\n",
" Msigi = mesh.getFaceInnerProduct(1./sigma)\n",
" MsigI = Utils.sdInv(Msigi)\n",
" A = D*MsigI*G\n",
" A[-1,-1] /= mesh.vol[-1] # Remove null space\n",
" rhs = np.zeros(mesh.nC)\n",
" txind = Utils.meshutils.closestPoints(mesh, pts)\n",
" rhs[txind] = np.r_[1,-1]\n",
" return A, rhs"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 60
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"tM.vectorCCy"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 61,
"text": [
"array([ 0.03333333, 0.1 , 0.16666667, 0.23333333, 0.3 ,\n",
" 0.36666667, 0.43333333, 0.5 , 0.56666667, 0.63333333,\n",
" 0.7 , 0.76666667, 0.83333333, 0.9 , 0.96666667])"
]
}
],
"prompt_number": 61
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"pts = np.vstack((np.r_[0.25, 0.5], np.r_[0.75, 0.5]))"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 62
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"AtM, rhstM = DCfun(tM, pts)\n",
"AinvtM = SolverLU(AtM)\n",
"phitM = AinvtM*rhstM"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 63
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"AqM, rhsqM = DCfun(qM, pts)\n",
"AinvqM = SolverLU(AqM)\n",
"phiqM = AinvqM*rhsqM"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 64
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"ArM, rhsrM = DCfun(rM, pts)\n",
"AinvrM = SolverLU(ArM)\n",
"phirM = AinvrM*rhsrM"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 65
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"coreind = (qM.gridCC[:,0]-0.5)**2+(qM.gridCC[:,1]-0.5)**2 >0.43**2\n",
"phiqM[coreind] = 0."
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 66
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"Xi = tM.gridCC[:,0].reshape(sz[0], sz[1], order='F')\n",
"Yi = tM.gridCC[:,1].reshape(sz[0], sz[1], order='F')\n",
"PHItM = griddata(tM.gridCC[:,0], tM.gridCC[:,1], phitM, Xi, Yi, interp='linear')\n",
"PHIqM = griddata(qM.gridCC[:,0], qM.gridCC[:,1], phiqM, Xi, Yi, interp='linear')\n",
"PHIrM = griddata(rM.gridCC[:,0], rM.gridCC[:,1], phirM, Xi, Yi, interp='linear')"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 67
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"fig, axes = plt.subplots(1,3,figsize=(14*1.2,4*1.2))\n",
"label = [\"(a)\", \"(b)\", \"(c)\"]\n",
"opts = {}\n",
"vmin, vmax = PHItM.min(), PHItM.max()\n",
"dat = axes[0].contourf(Xi, Yi, PHItM, 100)\n",
"tM.plotGrid(ax=axes[0], **opts)\n",
"axes[0].set_title('TensorMesh')\n",
"axes[1].contourf(Xi, Yi, PHIqM, 100)\n",
"qM.plotGrid(ax=axes[1], **opts)\n",
"axes[1].set_title('TreeMesh')\n",
"axes[2].contourf(Xi, Yi, PHIrM, 100)\n",
"rM.plotGrid(ax=axes[2], **opts)\n",
"axes[2].set_title('CurvilinearMesh')\n",
"for i in range(3):\n",
" axes[i].set_xlim(0.025, 0.975)\n",
" axes[i].set_ylim(0.025, 0.975)\n",
" axes[i].text(0., 1.0, label[i], fontsize=24)\n",
" if i==0: \n",
" axes[i].set_ylabel(\"y\")\n",
" else:\n",
" axes[i].set_ylabel(\" \")\n",
" axes[i].set_xlabel(\"x\")\n",
"plt.show()\n",
"# fig.savefig(\"./ThreeMesh.png\", dpi=100)"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 68
},
{
"cell_type": "code",
"collapsed": false,
"input": [],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 68
},
{
"cell_type": "code",
"collapsed": false,
"input": [],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 68
},
{
"cell_type": "code",
"collapsed": false,
"input": [],
"language": "python",
"metadata": {},
"outputs": []
}
],
"metadata": {}
}
]
}
notebooks/comp_dc.png
BIN
View File
Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 11 KiB

notebooks/mappingDC.png
BIN
View File
Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 103 KiB

notebooks/obspredDC.png
BIN
View File
Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 146 KiB

simpegDC/BaseDC.py
+8 -5
View File
@@ -1,17 +1,20 @@
from SimPEG import *
class DipoleSrc(Survey.BaseSrc):
class SrcDipole(Survey.BaseSrc):
"""A dipole source, locA and locB are moved to the closest cell-centers"""
current = 1
loc = None
_rhsDict = None
def __init__(self, rxList, locA, locB, **kwargs):
super(DipoleSrc, self).__init__(rxList, **kwargs)
self.loc = (locA, locB)
self._rhsDict = {}
super(SrcDipole, self).__init__(rxList, **kwargs)
def getRhs(self, mesh):
if getattr(self, '_rhsDict', None) is None:
self._rhsDict = {}
if mesh not in self._rhsDict:
pts = [self.loc[0], self.loc[1]]
inds = Utils.closestPoints(mesh, pts)
@@ -21,12 +24,12 @@ class DipoleSrc(Survey.BaseSrc):
return self._rhsDict[mesh]
class DipoleRx(Survey.BaseRx):
class RxDipole(Survey.BaseRx):
"""A dipole source, locA and locB are moved to the closest cell-centers"""
def __init__(self, locsM, locsN, **kwargs):
locs = (locsM, locsN)
assert locsM.shape == locsN.shape, 'locs must be the same shape.'
super(DipoleRx, self).__init__(locs, 'dipole', storeProjections=False, **kwargs)
super(RxDipole, self).__init__(locs, 'dipole', storeProjections=False, **kwargs)
@property
def nD(self):
simpegDC/Examples/Verification.py
+8 -2
View File
@@ -2,6 +2,7 @@ from SimPEG import *
import simpegDC as DC
import matplotlib.pyplot as plt
def run(plotIt=False):
cs = 25.
hx = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
@@ -22,11 +23,16 @@ def run(plotIt=False):
# ax.plot(xyz_rxP[:,0],xyz_rxP[:,1], 'w.')
# ax.plot(xyz_rxN[:,0],xyz_rxN[:,1], 'r.', ms = 3)
rx = DC.DipoleRx(xyz_rxP, xyz_rxN)
src = DC.DipoleSrc([rx], [-200, 0, -12.5],[+200, 0, -12.5])
rx = DC.RxDipole(xyz_rxP, xyz_rxN)
src = DC.SrcDipole([rx], [-200, 0, -12.5], [+200, 0, -12.5])
survey = DC.SurveyDC([src])
problem = DC.ProblemDC(mesh)
problem.pair(survey)
try:
from pymatsolver import MumpsSolver
problem.Solver = MumpsSolver
except Exception, e:
pass
data = survey.dpred(sigma)
def DChalf(srclocP, srclocN, rxloc, sigma, I=1.):
simpegDC/Examples/WennerArray.py
+2 -2
View File
@@ -32,8 +32,8 @@ def getSrcList(nElecs, aSpacing, in2D=False, plotIt=False):
getLoc = lambda ii, abmn: np.r_[elocs[WENNER[ii,abmn]],0, 0]
srcList = []
for i in range(WENNER.shape[0]):
rx = DC.DipoleRx(getLoc(i,1),getLoc(i,2))
src = DC.DipoleSrc([rx],getLoc(i,0),getLoc(i,3))
rx = DC.RxDipole(getLoc(i,1),getLoc(i,2))
src = DC.SrcDipole([rx], getLoc(i,0),getLoc(i,3))
srcList += [src]
return srcList