removedFWI example to merge into master.

This commit is contained in:
Rowan Cockett committed 2013-11-12 13:55:03 -08:00
1 parent 5248ac81ce
commit 7d3e41d591
4 files changed
+43 -318

No files matched your search

-180
View File
@@ -1,180 +0,0 @@
{
"metadata": {
"name": "FWI"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "code",
"collapsed": false,
"input": [
"import sys\n",
"sys.path.append('/Users/haber/dropbox/simpegMaster/')\n",
"import SimPEG\n"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"Warning: you may not have the required pubsub installed, use pypubsub. You will not be able to listen to events.\n"
]
}
],
"prompt_number": 6
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import scipy.sparse as sp"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 9
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"n1 = 128; n2 = 128;\n",
"h1 = 1.0/n1; h2 = 1.0/n2;\n",
"w = 100; m = np.ones(n1*n2)\n",
"mbc = 1;\n",
"k = w*np.sqrt(w**2 * mbc)\n",
"import numpy as np\n",
"import SimPEG.utils as utils\n",
"\n",
"mesh = SimPEG.mesh.TensorMesh([h1*np.ones(n1),h2*np.ones(n2)])"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 86
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"L1 = (1.0/h1**2) * sp.spdiags((np.ones((n1+1,1))*[1, -2, 1]).T,[-1,0,1],n1+1,n1+1);\n",
"L2 = (1.0/h2**2) * sp.spdiags((np.ones((n2+1,1))*[1, -2, 1]).T,[-1,0,1],n2+1,n2+1);\n",
"\n",
"av1 = sp.spdiags((np.ones((n1+1,1))*[0.5, 0.5]).T,[0,1],n1,n1+1);\n",
"av2 = sp.spdiags((np.ones((n2+1,1))*[0.5, 0.5]).T,[0,1],n2,n2+1);\n"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 87
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"B1 = utils.spzeros(n1+1,n1+1);\n",
"B2 = utils.spzeros(n2+1,n2+1);\n",
"B1.dtype = complex\n",
"B2.dtype = complex\n"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 93
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"B1[0,0] = 1+2*1j*k*h1; B1[-1,-1] = 1+2*1j*k*h1\n",
"B2[0,0] = 1+2*1j*k*h2; B2[-1,-1] = 1+2*1j*k*h2"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 94
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# generate the 2D Laplacian\n",
"L = sp.kron(sp.identity(n2+1),L1+B1) + sp.kron(L2+B2,sp.identity(n1+1)) \n",
"\n",
"# Generate the Mass matrix\n",
"Av = sp.kron(av2,av1)\n",
"M = utils.sdiag(Av.T*utils.mkvc(m))\n",
"A = L + w**2 * M"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 95
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"q = np.zeros((n1+1,n2+1)); q[64,64] = 1; q = utils.mkvc(q);\n",
"u = sp.linalg.spsolve(A,q);\n",
"u = u.reshape((n1+1,n2+1),order='F')\n",
"\n"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 96
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"imshow(u.real)\n"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "pyout",
"prompt_number": 97,
"text": [
"<matplotlib.image.AxesImage at 0x113625ed0>"
]
},
{
"output_type": "display_data",
"png": "iVBORw0KGgoAAAANSUhEUgAAAQIAAAD8CAYAAACcoKqNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvNmPZWl67vX71jzseYq9Y55yrKzMyqqu6nb70K42lm1h\nhEEgYRA6XPrO4saXB/e549LG/wAgBEI60oED4iD5iDZN2z3UPOcQmTFH7Ig9T2te6+Ni7QrryO4j\nkMoqS+QrhSIjcsWavvd9vvd7nufbQkopeRWv4lX8/zqUb/sGXsWreBXffrwCglfxKl7FKyB4Fa/i\nVbwCglfxKl4Fr4DgVbyKV8ErIHgVr+JV8PcABD/5yU+4d+8et27d4s///M+/6dO/ilfxKv4eQnzT\nPoLHjx/zZ3/2Z2xtbfE7v/M7/PSnP6XRaHyTl3gVr+JVfMPxjXYEk8kEgB/84AdsbW3x27/92/zi\nF7/4Ji/xKl7Fq/h7CO2bPNl7773H3bt3b36+f/8+P//5z/m93/u9m98JIb7JS76KV/Eq/j/G37UI\n+EaB4P9t/Jd/IPiniWT2JSwuobwB1luC7D8S/OLh2/yz1d/nJ5/9kI/eewf5lwL5REAE7ErE27D2\nbx9y53tf8Pv8C359+DP2nxzjHnmIEwklyNoK/lsaBxvb/Fj5IX8ZvctfBu8SflAk/sKCoQAB1CXi\nvX+C/Wf/BW9b7/OW8QFv8DF3sqdsJccUggDdk3glnZnj0GWFLm0u6TCgwYgqARYpKiopGgkqKQKJ\nROBjs8BlTIUJZSaUCTFIUSkwp86QDU7Z4ZA9XvAXP/olP/pRjJmEpFKlr9V5Lm7xCY/4mfw+7y/e\nZvGTCuFHDvIC0EGsQOHdEZ1fO+Xf55/zW8Mf8+bnn1H95RjlryQsABXkY8HonTLPf7DF/1r7d/lv\n+c8Z/TctFv9TCZ5E0AMw4JEA9Ufs/Vd/wBvff5//hP+B7x2+R+PHE4xnCZxJ2BBEtzV6Pyzz8+13\n+B/5T/n4r77Dy//9NvxYwKcZEENTwD0d9w+m1P7xFf+Y/45/b/C/sf+TI6q/nCI+kpABBch+XTB8\np8KHr73Ov6r9Jv8z/wHdv15n9k/+a3j9TyAGsQbmGx7ub4x52/klvyb+mkd8wr48oJEMUEVKqJmc\nsc4xWxywzyE7nLLBkCpziqgkmETL0ZhQYYzLAhv/ZtxSVBK0m3G1CPjgR/+S//hH+3S4pLPMhJLn\nYU8TYkcws2xOtE2eKHf5mDf4IHqL94O3iD50Sb60oL+cAGsS/UGA+eaMH1p/ybvGj/lh9mP2To+x\n309QrjKYgdwQLHYcDu5u8dPq9/lf+H2e/fw+5/9qG/nf/1MwfgQGiHsS8a7kzbd/wQ9e+z/5Dy//\nBd/95H2UfyYJPpBMTsFdg+J9SP8zBe33s7+zJr9RIHj77bf54z/+45ufv/jiC373d3/3bx0np5LZ\nERydw+UMdoB2XeIMJJmv4EmbaGCQHQh4GcJxAilg6siOSTZTkSg4eBTDGUY3Rj3K4AiogZAp1kxi\nRimZoRJFFsGkSHplIE8UGC2fXAoIFaLYAAMMIhwWFJIFzjzC9BPUQGLYoMkEKRQiTBa4TCkyoUyM\nDoBJiEqKToxKgkKOujE6ChkpKgEWXuoQZCaJqqEpKVVGROj532YJ5cUCbZaQZQqiPGBo1ihpM7Qs\nIUk0sqGCPBPwgrzAhxDsOMxuVYmLJjIWKJMMcSnJDiCYQ6yCU5Zomwl27KMkGYvYJe7qcBTD8BwW\nC8CBYQVMkHOBkJJisqA4nqMdpsivJOkhqDOJpqWUHs8pxguElh/PGTAcwGIMeKC6cLhO3NXxfAdV\nT7EjD62bkLyUeAegp2AVQexIlImERBCHFrNZleDYgT7wRIAEmUG2pZAkGlqWUFRn1JMhzXBAYeKj\nqBlWIWFgBihaSozOApcJZcZZlXnqYikhqAtSVBQyTEJsfBwWZKg3IBCjE2AhkISYN+MeYZAh0GSC\nHiXoiwQlFWRJiFvwcI0FBhFkEMUm6UxD9hW4BBJgIUiaBkwKRIpJpqmYUYo1jVEvM8QpeX5mEt2O\nKW7PsfHIEGQzBXkOTNN8YFWQpoZcN4l2DXwcpK/AIGN+CN1TOJxCR8K2AOv87wYB+IaBoFwuA7ly\nsLm5yV/8xV/wJ3/yJ3/ruMzPO4HLGbzIoDCD0gjMAOJEY4FLPDHgVEI3gn6Yz+ATCRMDQlDI8gGM\nfJRRlidMb3mBCig+iEQh0XXiyCKZ2cixyF/yBDCBAMhAZHmhf50UVhqi+xlKAFksSFOFRGr42MxE\nkRFVRtQYUkMibroAkxCdGJMQjRiBJEZHX/47QcNLHRaxS2qoaCJlJooE2KQoKFmGHUQoI4mMUzQt\npqaOKalTLAKULAOf/Bm6+b3jQ3xh4fVLBLpDnOnISCDnkA0hmIKvgt4HZZphJBEizYgCg3SiwjCB\nqAcMgAJE+cwlIokqU5w4wJkHKFeQnkNyAsIAtZ3hLALsOEBVU0QoYQpEE+AcmENUh2GHdKISBQZC\nkehpjDKVxAOY98CWYMSgLEBGgjjTCQIbr1civrDysUoBlRzkfYmSZlgElOWUWjKispii9UHooGsZ\nmpqRaQo+NlNKTGSZSVrGixykPsdQIwQSnRgbH5cFLoub4pcIIgw8HFJUBBIfmxFVqozxpUMiddJM\nIYsFCqCTYdkBFj4WAZpM8vGJl3k2I+9qgWyskcwU4oJFInVEoqB8fcxomccFUFcy7NjDJsgnlhAY\nA2EKYx8kUDLh1CCaGCykS5xoJEGOxdfzvL7EDJoXoA1/de1+40uDP/3TP+UP//APieOYP/qjP/o7\nFYPffAPKA9glB4GNNSjeAfbAbxmMRYUAi7z6bVAM0IGyCmsCqxhQZIZGkv++BLTIE6YN2ZbKuObS\ntypMlQKhYoAuoSygvTzeATbBWvsu9eI1VX2EywKVlEyDuChICyoxGn2rwYVo80Ls8ZIdjthmQJ0p\nZXRiLAJUEgrMMYhuWk1jCQABFnMKGEQoWUYSa/jSYZomjIwqA6XOkDr3f7OO53hYeoIWpWg+uEZA\nw+jTUq5pGVd06ybRqg3X5AmmAwlkC4VxWmFo14g3dcRtUO6COwYrA3MT5k2FQLfINBXb9UlWbNJ1\nEw53IVzLK6nswt676KUYU4RIXZK4Cko9Q22DCEBpg6xD4iqgS0wRoJdj2BJw1ATFBRmDY8K6jrYS\nYrkBmaoQ6hZZU8HagOoAVAWUKojbkGzoDKw646RC5in5DNp4F4w8DeiAVQ9YMa5oKdc05AAnCNF8\nEBISXSVwNCZ6kSF1htQYZVUmUQk/sUkSDaHmgF1gToUxFcaUmeDg4WPnMy8CH4sJZQIsYnTUd7/P\nIW0kClIIFDI67hXNTh+NhFQIMk2gkeKwoGqMWCleMe40mW8bOZB5y7wr5+MWKgZTpUDfqtCsFahs\nzVFllj9ra5mnGujElJhiFQNYE9D4TZgX8/eDCggCLEZKBb9lIvagdAc2RyAuoFXMl9/Gv0G8+8aB\n4Dd+4zf46quv/o3HfP+3TMxJxEpDUhxD8TaI7xr0tkpclZqMlSqBZUFVQEcHRQcX1Psp6gOfUnNC\nlREqKZGlEa+poGXQEkQNHb9tc1lpcaZ1mFAh0VVM1ydd08gyBSXJUJ0UvRXT2LxN2zqhzgAbH4nA\nU21GdolUUfFUhzPWOWEzX3NmO5ymG0xlCV86ONqCkjrFZYFE3ABBkRkWIRoJEQZzCjl4iBIIiGKD\neVpgqNXoKU26tHnzN3fop3NqzpSC9FGkxEpDaozoiEs29DO8jRLePZfE15ETJdd9LEhTlWFW48Lp\ncL3axHoQ4Cw89FGGEUvSOyrBtsncLCJVQVmdEG07hN9xodjIGwJNoDxK0d79Hk7zY2zhEaoGi4qD\nettHU1LUJsg1SPdUFmWHUDWwhYfTXGA8DkjGDhlFSCTUgTsCc2tC2RiToTCzCvjbJs7cx3ZSpCFI\nawreA4er1SaXVpvhvEaaqHlB3H0379oqGdq9mPL6iHXjjI64pCqH2GmIkJA5Cp5rMTSLXC/f5zUt\nhlmNRVggyvKOxBDRDQjUGVBnQIkpJiEKGREGIAixmFJimpbxEgf7+/8OB/GEUDWIFINEaPimTWhq\nOKmPmqX4qoUEHHwaap8N9RStDSISxBWd1FfJNAVlNUN1E2JNYywqnGmrlCpTsv0eTsFH78RQkUQr\nKrGVLxurjCg1JxgPAtLrH5Jaas7/tIEqhJbFWKlwVWrS36pT/d6UuhJjPAWnCtaOINg0yNuKvx3f\nClnYe63KqnKNM5CYAbAH/c0iH7cf8YX2GkOqhHUT7gFF8ntvgPGaT/GdESvlS9p0AZiVCtgPfPTb\nWl4MWpWe0eTYygu3T4PMEpT1PqFhEe9oGFqEqy2oGmPa5iVrnLJCF4uAEIu+0mQkasxxGVPllA2O\n2eIlu5wl63QXbYLYJks1KuURwpHUGQA5z1BgTpURCdpyfZe3qWMqzLQiY6tCMHMIYpeBXeeCVSqM\nqTKiqMyQjWPUUowVJRh6REWMWeeMff0Zk7slFmWXaaVKcm3ks8wqpIZKX23wwtjlo8YDwncU9u69\nwJlE6IFk0TAZVYv03RoZCh0u8B8WmRZqcCZgDlhg7Pq4r02oFAcUmOfJtdLCfPcc9fUUMQJZg7Cm\ncdVqMlYqFJlT2R5QrA1YrFcI/pELgYACsCEp7M7pcEmGYODWab7Zx9rxKX7PJ7YVvLLBi+IWXziv\ncaDv0lfqZKYCq8AjwAFtJaH4cEy7c84t7RnrnFERE3QjIi0JAkNlYJQ5FpscssMxW1zSYZjVCQIX\nVInpeBS1GXUGrHBFmy4rXOGyQCNBIljgAhBiMqPIIGwwnlRR1ATL8Jk6ZaZGOf8/6vRoUFLine truncated
"text": [
"<matplotlib.figure.Figure at 0x10109a5d0>"
]
}
],
"prompt_number": 97
},
{
"cell_type": "code",
"collapsed": false,
"input": [],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 97
},
{
"cell_type": "code",
"collapsed": false,
"input": [],
"language": "python",
"metadata": {},
"outputs": []
}
],
"metadata": {}
}
]
}
-101
View File
@@ -1,101 +0,0 @@
try:
import sys
sys.path.append('/Users/haber/dropbox/simpegMaster/')
except Exception, e:
pass
import SimPEG
import scipy.sparse as sp
import numpy as np
import SimPEG.utils as utils
import matplotlib.pylab as plt
from SimPEG.mesh import TensorMesh
def HelmholtzSol(model,mesh,w,mbc,q,P):
n1 = mesh.nCx; n2 = mesh.nCy;
h1 = mesh.hx; h2 = mesh.hy;
k = np.sqrt(w**2 * mbc)
# D1 = utils.sdiag(1./h1) * utils.ddx(mesh.nCx)
# D2 = utils.sdiag(1./h2) * utils.ddx(mesh.nCy)
# L1 = - D1.T*D1
# L2 = - D2.T*D2
Av = mesh.aveN2CC
B1 = utils.spzeros(n1+1,n1+1);
B2 = utils.spzeros(n2+1,n2+1);
B1.dtype = complex
B2.dtype = complex
B1[0,0] = 2*1j*k*h1[0]; B1[-1,-1] = 2*1j*k*h1[-1]
B2[0,0] = 2*1j*k*h2[0]; B2[-1,-1] = 2*1j*k*h2[-1]
# generate the 2D Laplacian
# L = sp.kron(sp.identity(n2+1),L1) + sp.kron(L2,sp.identity(n1+1))
L = mesh.nodalLaplacian
B = sp.kron(sp.identity(n2+1),B1) + sp.kron(B2,sp.identity(n1+1))
L = L+B
#plt.spy(L)
#plt.show()
# Generate the Mass matrix
M = utils.sdiag(Av.T*utils.mkvc(model))
A = - L - w**2 * M
mesh.ForModMat = A
u = sp.linalg.spsolve(A,q);
d = P*u
return u, d
def HelmholtzJmatVec(v,model,u,mesh,w,mbc,P):
Cm = -w**2 * utils.sdiag(u)*mesh.aveN2CC.T
Cu = mesh.ForModMat
Cmv = Cm*v;
lam = sp.linalg.spsolve(Cu,Cmv);
return -P*lam
def HelmholtzJTmatVec(v,model,u,mesh,w,mbc,P):
Cm = -w**2 * utils.sdiag(u)*mesh.aveN2CC.T
Cu = mesh.ForModMat
Pv = -P.T*v
z = sp.linalg.spsolve(Cu.T,Pv);
return Cm.T*z
if __name__ == '__main__':
# odel,mesh,w,mbc,q
n1 = 128; n2 = 128
h1 = np.ones(n1); h2 = np.ones(n2);
P = sp.identity((n1+1)*(n2+1))
mesh = TensorMesh([h1,h2])
model = np.ones(mesh.nC)
w = 1
mbc = 1
q = np.zeros((mesh.nNx,mesh.nNy))
q[n1/2,n2/2] = 1.0
q = q.reshape(mesh.nN,order = 'F')
u, d = HelmholtzSol(model,mesh,w,mbc,q,P)
u = u.reshape((mesh.nCx+1,mesh.nCy+1),order = 'F')
plt.imshow(u.real)
plt.show()
dm = np.random.rand(mesh.nC)*1e-1+2
u1, d1 = HelmholtzSol(model+dm,mesh,w,mbc,q,P)
dd = HelmholtzJmatVec(dm,model,u,mesh,w,mbc,P)
print np.linalg.norm(d1-d)
print np.linalg.norm(d1-d-dd)
-1
View File
@@ -5,7 +5,6 @@ from SimPEG.utils import ModelBuilder, sdiag, mkvc
from SimPEG import Solver
import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as linalg
class DCProblem(ModelTransforms.LogModel, Problem):
+43 -36
View File
@@ -4,47 +4,54 @@ import unittest
import HTMLTestRunner
# This code will run all tests in directory named test_*.py
def main(html=False):
TITLE = 'Test Results'
test_file_strings = glob.glob('test_*.py')
module_strings = [str[0:len(str)-3] for str in test_file_strings]
suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str
in module_strings]
testSuite = unittest.TestSuite(suites)
TITLE = 'Test Results'
test_file_strings = glob.glob('test_*.py')
module_strings = [str[0:len(str)-3] for str in test_file_strings]
suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str
in module_strings]
testSuite = unittest.TestSuite(suites)
unittest.TextTestRunner(verbosity=2).run(testSuite)
if not html:
unittest.TextTestRunner(verbosity=2).run(testSuite)
return
outfile = open("report.html", "w")
runner = HTMLTestRunner.HTMLTestRunner(
stream=outfile,
title=TITLE,
description='SimPEG Test Report was automatically generated.'
)
outfile = open("report.html", "w")
runner = HTMLTestRunner.HTMLTestRunner(
stream=outfile,
title=TITLE,
description='SimPEG Test Report was automatically generated.',
verbosity=2
)
runner.run(testSuite)
outfile.close()
runner.run(testSuite)
outfile.close()
reader = open("report.html", "r")
writer = open("../../docs/api_TestResults.rst", "w")
reader = open("report.html", "r")
writer = open("../../docs/api_TestResults.rst", "w")
writer.write('.. _api_TestResults:\n\nTest Results\n============\n\n.. raw:: html\n\n')
writer.write('.. _api_TestResults:\n\nTest Results\n============\n\n.. raw:: html\n\n')
go = False
for line in reader:
skip = False
if line == '<style type="text/css" media="screen">\n':
go = True
elif line == "<div id='ending'>&nbsp;</div>\n":
go = False
elif line == '</head>\n':
skip = True
elif line == '<h1>'+TITLE+'</h1>\n':
skip = True
elif line == '<body>\n':
skip = True
if go and not skip:
writer.write(' '+line)
go = False
for line in reader:
skip = False
if line == '<style type="text/css" media="screen">\n':
go = True
elif line == "<div id='ending'>&nbsp;</div>\n":
go = False
elif line == '</head>\n':
skip = True
elif line == '<h1>'+TITLE+'</h1>\n':
skip = True
elif line == '<body>\n':
skip = True
if go and not skip:
writer.write(' '+line)
writer.close()
reader.close()
os.remove("report.html")
writer.close()
reader.close()
os.remove("report.html")
if __name__ == '__main__':
main(True)