Skip to content

Instantly share code, notes, and snippets.

@sgkang

sgkang/DCsensTest.ipynb

Created March 3, 2015 20:02
Show Gist options
  • Save sgkang/96f9140785ebcc8e165c to your computer and use it in GitHub Desktop.
Save sgkang/96f9140785ebcc8e165c to your computer and use it in GitHub Desktop.
Download ZIP
DCSenseTest
Display the source blob
Display the rendered blob
Raw
DCsensTest.ipynb
{
"metadata": {
"language": "Julia",
"name": "",
"signature": "sha256:aefd00f6dc6e17fe33ad06592c57164c36e31178df25943802eb2c115cf1f1ca"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "code",
"collapsed": false,
"input": [
"using DCIP\n",
"using Mesh\n",
"using Base.Test\n",
"using Utils\n",
"using LinearSolvers\n",
"using MUMPS"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 1
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# using PyPlot"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 2
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"cs, npad, rate = 200., 2, 1.3\n",
"ncx, ncy, ncz = 15, 15, 10\n",
"hpad = cs*(1:npad).^rate;"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 3
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"hx = [hpad, cs*ones(ncx), hpad[end:-1:1]]; \n",
"hy = [hpad, cs*ones(ncy), hpad[end:-1:1]]; \n",
"hz = [hpad, cs*ones(ncz)]; \n",
"x0 = [-sum(hx)/2., -sum(hy)/2., -sum(hz)];\n",
"M = getTensorMesh3D(hx, hy, hz, x0);"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 4
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"xn, yn, zn = getNodalAxes(M);"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 5
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"XYZN = getNodalGrid(M)\n",
"xc1_p, yc1_p, zc1_p = -1500., 0., 0.; \n",
"xc1_n, yc1_n, zc1_n = 1500., 0., 0.; \n",
"xc2_p, yc2_p, zc2_p = 0., -1500., 0.; \n",
"xc2_n, yc2_n, zc2_n = 0., 1500., 0.; \n",
"(dum, indp1) = findmin(abs(XYZN[:,1]-xc1_p)+abs(XYZN[:,2]-yc1_p)+abs(XYZN[:,3]-zc1_p));\n",
"(dum, indn1) = findmin(abs(XYZN[:,1]-xc1_n)+abs(XYZN[:,2]-yc1_n)+abs(XYZN[:,3]-zc1_n));\n",
"(dum, indp2) = findmin(abs(XYZN[:,1]-xc2_p)+abs(XYZN[:,2]-yc2_p)+abs(XYZN[:,3]-zc2_p));\n",
"(dum, indn2) = findmin(abs(XYZN[:,1]-xc2_n)+abs(XYZN[:,2]-yc2_n)+abs(XYZN[:,3]-zc2_n));"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 6
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"Sources = zeros(Float64, M.nn, 2);\n",
"Sources[indp1,1] = 1./3;\n",
"Sources[indn1,1] = -1./3;\n",
"Sources[indp2,2] = 1./3;\n",
"Sources[indn2,2] = -1./3;"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 7
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"xr = linspace(-1000., 1000., 11);\n",
"(XrP, YrP, ZrP) = ndgrid(xr, xr, -1e-4*ones(1));\n",
"(XrN, YrN, ZrN) = ndgrid(xr+100., xr, -1e0*ones(1));\n",
"xyzrP = [XrP[:] YrP[:] ZrP[:]];\n",
"xyzrN = [XrN[:] YrN[:] ZrN[:]];\n",
"QP = getNodalInterpolationMatrix(M, xyzrP);\n",
"QN = getNodalInterpolationMatrix(M, xyzrN);"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 8
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"Q = blkdiag(QP-QN, QP-QN);"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 9
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"mumps = getMUMPSsolver([],1);\n",
"param = getDCIPParam(M, Sources, Q', [], mumps);"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 10
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"sigma = ones(M.nc);\n",
"dm = sigma*10;\n",
"sigma0 = sigma+dm\n",
"dobs, param = DCIPfwd(sigma, param);\n",
"dini, param = DCIPfwd(sigma0, param);"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"elapsed time: 0."
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"044402391 seconds\n",
"elapsed time: 0.064166271 seconds\n"
]
}
],
"prompt_number": 19
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"dr = dini-dobs\n",
"Jtdr = DCIPSensTMatVec(dr, sigma0, param);"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 31
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"Step1: Test Jtv"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"function testDCJtv(x, v, Jtdr, dobs, param) \n",
" pred, param = DCIPfwd(x, param)\n",
" misfit = norm(pred-dobs)^2*0.5\n",
" dmisfit = dot(Jtdr, v)\n",
" return misfit, dmisfit\n",
"end"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 32,
"text": [
"testDCJtv (generic function with 1 method)"
]
}
],
"prompt_number": 32
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"testDCJtvfun(x, v=zeros(size(x))) = testDCJtv(x, v, Jtdr, dobs, param)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 33,
"text": [
"testDCJtvfun (generic function with 2 methods)"
]
}
],
"prompt_number": 33
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"pass, err, order = checkDerivative(testDCJtvfun, sigma0)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
" h\t E0\t E1\t O1\t O2\t OK?\n"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"elapsed time: 0.058913509 seconds\n",
"elapsed time: "
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"0.060815606 seconds\n",
"1.000e-01\t2.863e-04\t2.526e-06\t0.000e+00\t0.000e+00\t 1\n",
"elapsed time: "
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"0.048162417 seconds\n",
"1.000e-02\t2.841e-05\t2.519e-08\t1.003e+00\t2.001e+00\t 1\n",
"elapsed time: 0.045833286 seconds\n",
"1.000e-03\t2.838e-06\t2.519e-10\t1.000e+00\t2.000e+00\t 1\n",
"elapsed time: 0.030228583 seconds\n",
"1.000e-04\t2.838e-07\t2.517e-12\t1.000e+00\t2.000e+00\t 1\n",
"elapsed time: "
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"0.031033603 seconds\n",
"1.000e-05\t2.838e-08\t2.451e-14\t1.000e+00\t2.012e+00\t 1\n",
"elapsed time: 0.044765998 seconds\n",
"1.000e-06\t2.838e-09\t3.860e-16\t1.000e+00\t1.803e+00\t 1\n",
"elapsed time: 0.030751767 seconds\n",
"1.000e-07\t2.838e-10\t5.789e-16\t1.000e+00\t-1.761e-01\t 1\n",
"elapsed time: "
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"0.032226345 seconds\n",
"1.000e-08\t2.838e-11\t9.649e-16\t1.000e+00\t-2.218e-01\t 1\n",
"elapsed time: 0.042841435 seconds\n",
"1.000e-09\t2.838e-12\t5.789e-16\t1.000e+00\t2.218e-01\t 1\n",
"elapsed time: "
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"0.031701324 seconds\n",
"1.000e-10\t2.831e-13\t7.719e-16\t1.001e+00\t-1.249e-01\t 1\n"
]
},
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 34,
"text": [
"(true,\n",
"10x2 Array{Float64,2}:\n",
" 0.000286347 2.52603e-6 \n",
" 2.84073e-5 2.51931e-8 \n",
" 2.83846e-6 2.51864e-10\n",
" 2.83824e-7 2.51742e-12\n",
" 2.83821e-8 2.45085e-14\n",
" 2.83821e-9 3.8596e-16 \n",
" 2.8382e-10 5.7894e-16 \n",
" 2.83812e-11 9.649e-16 \n",
" 2.83758e-12 5.7894e-16 \n",
" 2.83102e-13 7.7192e-16 ,\n",
"\n",
"10x2 Array{Float64,2}:\n",
" 0.0 0.0 \n",
" 1.00346 2.00116 \n",
" 1.00035 2.00012 \n",
" 1.00003 2.00021 \n",
" 1.0 2.01164 \n",
" 1.0 1.80277 \n",
" 1.0 -0.176091\n",
" 1.00001 -0.221849\n",
" 1.00008 0.221849\n",
" 1.00101 -0.124939)"
]
}
],
"prompt_number": 34
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"Step2: Test Jvec"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"sigma = ones(M.nc);\n",
"dm = sigma*10;\n",
"sigma0 = sigma+dm\n",
"dobs, param = DCIPfwd(sigma, param);\n",
"dini, param = DCIPfwd(sigma0, param);"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"elapsed time: 0."
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"074214266 seconds\n",
"elapsed time: 0.072871131 seconds\n"
]
}
],
"prompt_number": 35
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"function testDCJv(x, v, param) \n",
" pred, param = DCIPfwd(x, param)\n",
" Jdm = DCIPSensMatVec(v, x, param);\n",
" return pred, Jdm\n",
"end"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 46,
"text": [
"testDCJv (generic function with 2 methods)"
]
}
],
"prompt_number": 46
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"testDCJvfun(x, v=zeros(size(x))) = testDCJv(x, v, param)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 47,
"text": [
"testDCJvfun (generic function with 2 methods)"
]
}
],
"prompt_number": 47
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"pass, err, order = checkDerivative(testDCJvfun, sigma0)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
" h\t E0\t E1\t O1\t O2\t OK?\n"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"elapsed time: 0.037673534 seconds\n",
"elapsed time: "
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"0.036428875 seconds\n",
"1.000e-01\t2.902e-03\t1.790e-05\t0.000e+00\t0.000e+00\t 1"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"elapsed time: 0.043053445 seconds\n",
"1.000e-02\t2.900e-04\t1.788e-07\t1.000e+00\t2.001e+00\t 1\n",
"elapsed time: "
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"0.032310843 seconds\n",
"1.000e-03\t2.900e-05\t1.788e-09\t1.000e+00\t2.000e+00\t 1\n",
"elapsed time: "
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"0.031032602 seconds\n",
"1.000e-04\t2.900e-06\t1.788e-11\t1.000e+00\t2.000e+00\t 1\n",
"elapsed time: 0.031137021 seconds\n",
"1.000e-05\t2.900e-07\t1.759e-13\t1.000e+00\t2.007e+00\t 1"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"elapsed time: 0.042587977 seconds\n",
"1.000e-06\t2.900e-08\t6.131e-15\t1.000e+00\t1.458e+00\t 1\n",
"elapsed time: "
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"0.031519125 seconds\n",
"1.000e-07\t2.900e-09\t7.374e-15\t1.000e+00\t-8.014e-02\t 1\n",
"elapsed time: 0."
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"030675628 seconds\n",
"1.000e-08\t2.900e-10\t6.150e-15\t1.000e+00\t7.882e-02\t 1\n",
"elapsed time: 0.031992284 seconds\n",
"1.000e-09\t2.899e-11\t6.475e-15\t1.000e+00\t-2.238e-02\t 1"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"\n",
"elapsed time: 0.039846114 seconds\n",
"1.000e-10\t2.899e-12\t6.173e-15\t1.000e+00\t2.078e-02\t 1\n"
]
},
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 48,
"text": [
"(true,\n",
"10x2 Array{Float64,2}:\n",
" 0.0029019 1.79028e-5 \n",
" 0.000289973 1.78815e-7 \n",
" 2.89953e-5 1.78795e-9 \n",
" 2.89951e-6 1.78765e-11\n",
" 2.89951e-7 1.75874e-13\n",
" 2.89951e-8 6.13141e-15\n",
" 2.8995e-9 7.37395e-15\n",
" 2.8995e-10 6.15009e-15\n",
" 2.89946e-11 6.47525e-15\n",
" 2.8992e-12 6.17268e-15,\n",
"\n",
"10x2 Array{Float64,2}:\n",
" 0.0 0.0 \n",
" 1.00032 2.00052 \n",
" 1.00003 2.00005 \n",
" 1.0 2.00007 \n",
" 1.0 2.00708 \n",
" 1.0 1.45764 \n",
" 1.0 -0.0801397\n",
" 1.0 0.0788191\n",
" 1.00001 -0.0223753\n",
" 1.00004 0.0207826)"
]
}
],
"prompt_number": 48
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"Step3: Adjoint test"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"#vtJtJv\n",
"v = randn(size(sigma0))\n",
"Jv = DCIPSensMatVec(v, sigma0, param);\n",
"JtJv = DCIPSensTMatVec(Jv, sigma0, param);\n",
"vtJtJv = dot(v, JtJv)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 51,
"text": [
"9.439497936091597e-13"
]
}
],
"prompt_number": 51
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"#vJJtv\n",
"v = randn(size(dobs))\n",
"Jtv = DCIPSensTMatVec(v, sigma0, param);\n",
"JJtv = DCIPSensMatVec(Jtv, sigma0, param);\n",
"vtJJtv = dot(v, JJtv)"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "pyout",
"prompt_number": 58,
"text": [
"9.476487711215165e-13"
]
}
],
"prompt_number": 58
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"println(norm(vtJtJv-vtJJtv)/norm(vtJJtv))"
],
"language": "python",
"metadata": {},
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"0.0"
]
},
{
"output_type": "stream",
"stream": "stdout",
"text": [
"03903321172441546\n"
]
}
],
"prompt_number": 60
}
],
"metadata": {}
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment