mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-06 05:16:51 +08:00
Corrected the notebook for MT3D derivative tests
This commit is contained in:
GudniRos
@@ -44,11 +44,11 @@
|
||||
"==================== checkDerivative ====================\n",
|
||||
"iter h |ft-f0| |ft-f0-h*J0*dx| Order\n",
|
||||
"---------------------------------------------------------\n",
|
||||
" 0 1.00e-01 3.058e+11 3.238e+10 nan\n",
|
||||
" 1 1.00e-02 3.096e+10 3.246e+08 1.999\n",
|
||||
" 2 1.00e-03 3.102e+09 3.251e+06 1.999\n",
|
||||
" 0 1.00e-01 3.111e+11 3.481e+10 nan\n",
|
||||
" 1 1.00e-02 3.112e+10 3.421e+08 2.008\n",
|
||||
" 2 1.00e-03 3.116e+09 3.420e+06 2.000\n",
|
||||
"========================= PASS! =========================\n",
|
||||
"Happy little convergence test!\n",
|
||||
"Once upon a time, a happy little test passed.\n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
@@ -103,11 +103,11 @@
|
||||
"==================== checkDerivative ====================\n",
|
||||
"iter h |ft-f0| |ft-f0-h*J0*dx| Order\n",
|
||||
"---------------------------------------------------------\n",
|
||||
" 0 1.00e-01 3.151e+11 3.512e+10 nan\n",
|
||||
" 1 1.00e-02 3.103e+10 3.483e+08 2.004\n",
|
||||
" 2 1.00e-03 3.101e+09 3.486e+06 2.000\n",
|
||||
" 0 1.00e-01 3.187e+11 3.398e+10 nan\n",
|
||||
" 1 1.00e-02 3.176e+10 3.369e+08 2.004\n",
|
||||
" 2 1.00e-03 3.178e+09 3.372e+06 2.000\n",
|
||||
"========================= PASS! =========================\n",
|
||||
"Yay passed!\n",
|
||||
"That was easy!\n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
@@ -164,11 +164,11 @@
|
||||
"==================== checkDerivative ====================\n",
|
||||
"iter h |ft-f0| |ft-f0-h*J0*dx| Order\n",
|
||||
"---------------------------------------------------------\n",
|
||||
" 0 1.00e-01 2.848e+10 2.846e+09 nan\n",
|
||||
" 1 1.00e-02 2.837e+09 2.801e+07 2.007\n",
|
||||
" 2 1.00e-03 2.838e+08 2.800e+05 2.000\n",
|
||||
" 0 1.00e-01 2.873e+10 3.211e+09 nan\n",
|
||||
" 1 1.00e-02 2.879e+09 3.204e+07 2.001\n",
|
||||
" 2 1.00e-03 2.882e+08 3.207e+05 2.000\n",
|
||||
"========================= PASS! =========================\n",
|
||||
"That was easy!\n",
|
||||
"You are awesome.\n",
|
||||
"\n"
|
||||
]
|
||||
}
|
||||
@@ -200,64 +200,6 @@
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"src = survey.srcList[0]\n",
|
||||
"u = problem.fields(problem.curModel)\n",
|
||||
"u_src = u[src,:]\n",
|
||||
"w = np.random.rand(problem.mesh.nC)\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 7,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"dA_dm = problem.getADeriv_m(src.freq,u_src,w)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 8,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"(20536, 2)"
|
||||
]
|
||||
},
|
||||
"execution_count": 8,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"dA_dm.shape"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 9,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"a, b, sig, d, e = t3Dmt.halfSpace(.01)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 21,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
@@ -266,26 +208,28 @@
|
||||
"==================== checkDerivative ====================\n",
|
||||
"iter h |ft-f0| |ft-f0-h*J0*dx| Order\n",
|
||||
"---------------------------------------------------------\n",
|
||||
" 0 1.00e-01 3.699e-07 6.864e-07 nan\n",
|
||||
" 1 1.00e-02 3.693e-08 6.856e-08 1.001\n",
|
||||
" 2 1.00e-03 3.693e-09 6.855e-09 1.000\n",
|
||||
" 3 1.00e-04 3.693e-10 6.855e-10 1.000\n",
|
||||
" 4 1.00e-05 3.693e-11 6.855e-11 1.000\n",
|
||||
" 5 1.00e-06 3.693e-12 6.855e-12 1.000\n",
|
||||
" 6 1.00e-07 3.693e-13 6.855e-13 1.000\n",
|
||||
" 7 1.00e-08 3.693e-14 6.855e-14 1.000\n",
|
||||
"*********************************************************\n",
|
||||
"<<<<<<<<<<<<<<<<<<<<<<<<< FAIL! >>>>>>>>>>>>>>>>>>>>>>>>>\n",
|
||||
"*********************************************************\n",
|
||||
"You had so much promise Gudni, oh well...\n",
|
||||
" 0 1.00e-01 1.436e+01 1.792e-11 nan\n",
|
||||
" 1 1.00e-02 1.436e+00 1.792e-11 0.000\n",
|
||||
" 2 1.00e-03 1.436e-01 1.730e-11 0.015\n",
|
||||
"========================= PASS! =========================\n",
|
||||
"That was easy!\n",
|
||||
"\n",
|
||||
"==================== checkDerivative ====================\n",
|
||||
"iter h |ft-f0| |ft-f0-h*J0*dx| Order\n",
|
||||
"---------------------------------------------------------\n",
|
||||
" 0 1.00e-01 1.433e+01 1.792e-11 nan\n",
|
||||
" 1 1.00e-02 1.433e+00 1.776e-11 0.004\n",
|
||||
" 2 1.00e-03 1.433e-01 1.763e-11 0.003\n",
|
||||
"========================= PASS! =========================\n",
|
||||
"The test be workin!\n",
|
||||
"\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"def testDeriv_ProjFields(prb):\n",
|
||||
"def testDeriv_epx(prb):\n",
|
||||
" # Initate things for the derivs Test\n",
|
||||
" src = survey.srcList[0]\n",
|
||||
" src = prb.survey.srcList[0]\n",
|
||||
" rx = src.rxList[0]\n",
|
||||
" rx.locs = np.array([[0.,0.,0,],[1.,1.,1.]])\n",
|
||||
" u0x = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j\n",
|
||||
@@ -303,10 +247,379 @@
|
||||
" f[src,'e_pxSolution'] = u[:len(u)/2]\n",
|
||||
" f[src,'e_pySolution'] = u[len(u)/2::]\n",
|
||||
" # f[src,'b_1d'] = -(m1d.nodalGrad*u)/(1j*simpegem.Utils.EMUtils.omega(src.freq))\n",
|
||||
" return rx.projectFields(src,survey.mesh,f), lambda t: rx.projectFieldsDeriv(src,survey.mesh,f0,t)\n",
|
||||
" simpeg.Tests.checkDerivative(fun,u0,num=8,plotIt=False)\n",
|
||||
" return f._e_px(f[src,'e_pxSolution'],[src]).ravel(), lambda t: f._e_pxDeriv_u([src],t)[:len(u0)/2].ravel()\n",
|
||||
" simpeg.Tests.checkDerivative(fun,u0,num=3,plotIt=False)\n",
|
||||
" \n",
|
||||
"def testDeriv_epy(prb):\n",
|
||||
" # Initate things for the derivs Test\n",
|
||||
" src = prb.survey.srcList[0]\n",
|
||||
" rx = src.rxList[0]\n",
|
||||
" rx.locs = np.array([[0.,0.,0,],[1.,1.,1.]])\n",
|
||||
" u0x = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j\n",
|
||||
" u0y = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j\n",
|
||||
"# u0 = np.vstack((simpeg.mkvc(u0x,2),simpeg.mkvc(u0y,2)))\n",
|
||||
" u0 = np.r_[u0x, u0y]\n",
|
||||
" f0 = prb.fieldsPair(survey.mesh,survey)\n",
|
||||
" # u0 = np.hstack((simpeg.mkvc(u0_px,2),simpeg.mkvc(u0_py,2)))\n",
|
||||
" f0[src,'e_pxSolution'] = u0[:len(u0)/2]#u0x\n",
|
||||
" f0[src,'e_pySolution'] = u0[len(u0)/2::]#u0y\n",
|
||||
" # f0[src,'b_1d'] = -1/(1j*simpegem.Utils.EMUtils.omega(src.freq))*m1d.nodalGrad*u0\n",
|
||||
" # Run a test\n",
|
||||
" def fun(u):\n",
|
||||
" f = prb.fieldsPair(survey.mesh,survey)\n",
|
||||
" f[src,'e_pxSolution'] = u[:len(u)/2]\n",
|
||||
" f[src,'e_pySolution'] = u[len(u)/2::]\n",
|
||||
" return f._e_py(f[src,'e_pySolution'],[src]).ravel(), lambda t: f._e_pyDeriv_u([src],t).ravel()\n",
|
||||
" simpeg.Tests.checkDerivative(fun,u0,num=3,plotIt=False)\n",
|
||||
" \n",
|
||||
"survey, problem = t3Dmt.setupSimpegMTfwd_eForm_ps(t3Dmt.random(.1),comp='zxyr',singleFreq=.01)\n",
|
||||
"testDeriv_epx(problem)\n",
|
||||
"testDeriv_epy(problem)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 7,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"==================== checkDerivative ====================\n",
|
||||
"iter h |ft-f0| |ft-f0-h*J0*dx| Order\n",
|
||||
"---------------------------------------------------------\n",
|
||||
" 0 1.00e-01 3.995e-01 1.107e-13 nan\n",
|
||||
" 1 1.00e-02 3.995e-02 1.105e-13 0.001\n",
|
||||
" 2 1.00e-03 3.995e-03 1.110e-13 -0.002\n",
|
||||
"========================= PASS! =========================\n",
|
||||
"The test be workin!\n",
|
||||
"\n",
|
||||
"==================== checkDerivative ====================\n",
|
||||
"iter h |ft-f0| |ft-f0-h*J0*dx| Order\n",
|
||||
"---------------------------------------------------------\n",
|
||||
" 0 1.00e-01 4.075e-01 1.086e-13 nan\n",
|
||||
" 1 1.00e-02 4.075e-02 1.049e-13 0.015\n",
|
||||
" 2 1.00e-03 4.075e-03 1.025e-13 0.010\n",
|
||||
"========================= PASS! =========================\n",
|
||||
"Happy little convergence test!\n",
|
||||
"\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"def testDeriv_bpx(prb):\n",
|
||||
" # Initate things for the derivs Test\n",
|
||||
" src = prb.survey.srcList[0]\n",
|
||||
" rx = src.rxList[0]\n",
|
||||
" rx.locs = np.array([[0.,0.,0,],[1.,1.,1.]])\n",
|
||||
" u0x = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j\n",
|
||||
" u0y = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j\n",
|
||||
"# u0 = np.vstack((simpeg.mkvc(u0x,2),simpeg.mkvc(u0y,2)))\n",
|
||||
" u0 = np.r_[u0x, u0y]\n",
|
||||
" f0 = prb.fieldsPair(survey.mesh,survey)\n",
|
||||
" # u0 = np.hstack((simpeg.mkvc(u0_px,2),simpeg.mkvc(u0_py,2)))\n",
|
||||
" f0[src,'e_pxSolution'] = u0[:len(u0)/2]#u0x\n",
|
||||
" f0[src,'e_pySolution'] = u0[len(u0)/2::]#u0y\n",
|
||||
" # f0[src,'b_1d'] = -1/(1j*simpegem.Utils.EMUtils.omega(src.freq))*m1d.nodalGrad*u0\n",
|
||||
" # Run a test\n",
|
||||
" def fun(u):\n",
|
||||
" f = prb.fieldsPair(survey.mesh,survey)\n",
|
||||
" f[src,'e_pxSolution'] = u[:len(u)/2]\n",
|
||||
" f[src,'e_pySolution'] = u[len(u)/2::]\n",
|
||||
" # f[src,'b_1d'] = -(m1d.nodalGrad*u)/(1j*simpegem.Utils.EMUtils.omega(src.freq))\n",
|
||||
" return f._b_px(f[src,'e_pxSolution'],[src]).ravel(), lambda t: f._b_pxDeriv_u(src,t).ravel()\n",
|
||||
" simpeg.Tests.checkDerivative(fun,u0,num=3,plotIt=False)\n",
|
||||
" \n",
|
||||
"def testDeriv_bpy(prb):\n",
|
||||
" # Initate things for the derivs Test\n",
|
||||
" src = prb.survey.srcList[0]\n",
|
||||
" rx = src.rxList[0]\n",
|
||||
" rx.locs = np.array([[0.,0.,0,],[1.,1.,1.]])\n",
|
||||
" u0x = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j\n",
|
||||
" u0y = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j\n",
|
||||
"# u0 = np.vstack((simpeg.mkvc(u0x,2),simpeg.mkvc(u0y,2)))\n",
|
||||
" u0 = np.r_[u0x, u0y]\n",
|
||||
" f0 = prb.fieldsPair(survey.mesh,survey)\n",
|
||||
" # u0 = np.hstack((simpeg.mkvc(u0_px,2),simpeg.mkvc(u0_py,2)))\n",
|
||||
" f0[src,'e_pxSolution'] = u0[:len(u0)/2]#u0x\n",
|
||||
" f0[src,'e_pySolution'] = u0[len(u0)/2::]#u0y\n",
|
||||
" # f0[src,'b_1d'] = -1/(1j*simpegem.Utils.EMUtils.omega(src.freq))*m1d.nodalGrad*u0\n",
|
||||
" # Run a test\n",
|
||||
" def fun(u):\n",
|
||||
" f = prb.fieldsPair(survey.mesh,survey)\n",
|
||||
" f[src,'e_pxSolution'] = u[:len(u)/2]\n",
|
||||
" f[src,'e_pySolution'] = u[len(u)/2::]\n",
|
||||
" return f._b_py(f[src,'e_pySolution'],[src]).ravel(), lambda t: f._b_pyDeriv_u(src,t).ravel()\n",
|
||||
" simpeg.Tests.checkDerivative(fun,u0,num=3,plotIt=False)\n",
|
||||
" \n",
|
||||
"survey, problem = t3Dmt.setupSimpegMTfwd_eForm_ps(t3Dmt.random(.1),comp='zyxr',singleFreq=.01)\n",
|
||||
"testDeriv_ProjFields(problem)"
|
||||
"testDeriv_bpx(problem)\n",
|
||||
"testDeriv_bpy(problem)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 8,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# %debug"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 9,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"==================== checkDerivative ====================\n",
|
||||
"iter h |ft-f0| |ft-f0-h*J0*dx| Order\n",
|
||||
"---------------------------------------------------------\n",
|
||||
" 0 1.00e-01 5.819e+11 7.675e+10 nan\n",
|
||||
" 1 1.00e-02 5.440e+10 8.852e+08 1.938\n",
|
||||
" 2 1.00e-03 5.403e+09 4.494e+07 1.294\n",
|
||||
" 3 1.00e-04 5.399e+08 4.428e+06 1.006\n",
|
||||
" 4 1.00e-05 5.399e+07 4.428e+05 1.000\n",
|
||||
" 5 1.00e-06 5.399e+06 4.428e+04 1.000\n",
|
||||
"*********************************************************\n",
|
||||
"<<<<<<<<<<<<<<<<<<<<<<<<< FAIL! >>>>>>>>>>>>>>>>>>>>>>>>>\n",
|
||||
"*********************************************************\n",
|
||||
"You break it, you fix it.\n",
|
||||
"\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"from scipy.constants import mu_0\n",
|
||||
"def testDeriv_Hd(prb):\n",
|
||||
" src = prb.survey.srcList[0]\n",
|
||||
" rx = src.rxList[0]\n",
|
||||
" rx.locs = np.array([[0.,0.,0,],[1.,1.,1.]])\n",
|
||||
" u0x = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j\n",
|
||||
" u0y = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j\n",
|
||||
"# u0 = np.vstack((simpeg.mkvc(u0x,2),simpeg.mkvc(u0y,2)))\n",
|
||||
" u0 = np.r_[u0x, u0y]\n",
|
||||
" f0 = prb.fieldsPair(survey.mesh,survey)\n",
|
||||
" # u0 = np.hstack((simpeg.mkvc(u0_px,2),simpeg.mkvc(u0_py,2)))\n",
|
||||
" f0[src,'e_pxSolution'] = u0[:len(u0)/2]#u0x\n",
|
||||
" f0[src,'e_pySolution'] = u0[len(u0)/2::]#u0y\n",
|
||||
" # Run a testdef testDeriv_ZijN(rx):\n",
|
||||
" def getHdcomp(rx,f,v=None):\n",
|
||||
" if rx.locs.ndim == 3:\n",
|
||||
" eFLocs = rx.locs[:,:,0]\n",
|
||||
" bFLocs = rx.locs[:,:,1]\n",
|
||||
" else:\n",
|
||||
" eFLocs = rx.locs\n",
|
||||
" bFLocs = rx.locs\n",
|
||||
" # Get the projection\n",
|
||||
" Pbx = prb.mesh.getInterpolationMat(bFLocs,'Fx')\n",
|
||||
" Pby = prb.mesh.getInterpolationMat(bFLocs,'Fy')\n",
|
||||
" # Get the fields at location\n",
|
||||
" # px: x-polaration and py: y-polaration.\n",
|
||||
" hx_px = Pbx*f[src,'b_px']/mu_0\n",
|
||||
" hy_px = Pby*f[src,'b_px']/mu_0\n",
|
||||
" hx_py = Pbx*f[src,'b_py']/mu_0\n",
|
||||
" hy_py = Pby*f[src,'b_py']/mu_0\n",
|
||||
" # Derivatives as lambda functions\n",
|
||||
" spPe = simpeg.Utils.spzeros(rx.nD,prb.mesh.nE)\n",
|
||||
" spPb = simpeg.Utils.spzeros(rx.nD,prb.mesh.nF)\n",
|
||||
"\n",
|
||||
" # NOTE: Think b_p?Deriv_u should return a 2*nF size matrix\n",
|
||||
" hx_px_u = lambda vec: Pbx*f._b_pxDeriv_u(src,vec)/mu_0\n",
|
||||
" hy_px_u = lambda vec: Pby*f._b_pxDeriv_u(src,vec)/mu_0\n",
|
||||
" hx_py_u = lambda vec: Pbx*f._b_pyDeriv_u(src,vec)/mu_0\n",
|
||||
" hy_py_u = lambda vec: Pby*f._b_pyDeriv_u(src,vec)/mu_0\n",
|
||||
" # Update the input vector\n",
|
||||
" sdiag = lambda t: simpeg.Utils.sdiag(simpeg.mkvc(t,2))\n",
|
||||
" # Define the components of the derivative\n",
|
||||
" if v is not None:\n",
|
||||
" return (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))\n",
|
||||
" else:\n",
|
||||
" return (sdiag(hx_px)*hy_py) - (sdiag(hx_py)*hy_px)\n",
|
||||
" def fun(u):\n",
|
||||
" f = prb.fieldsPair(survey.mesh,survey)\n",
|
||||
" f[src,'e_pxSolution'] = u[:len(u)/2]\n",
|
||||
" f[src,'e_pySolution'] = u[len(u)/2::]\n",
|
||||
" return getHdcomp(rx,f), lambda t: getHdcomp(rx,f0,t)\n",
|
||||
" simpeg.Tests.checkDerivative(fun,u0,num=6,plotIt=False)\n",
|
||||
"survey, problem = t3Dmt.setupSimpegMTfwd_eForm_ps(t3Dmt.random(.01),comp='zxyr',singleFreq=.0001)\n",
|
||||
"testDeriv_Hd(problem)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 10,
|
||||
"metadata": {
|
||||
"collapsed": false,
|
||||
"scrolled": true
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"==================== checkDerivative ====================\n",
|
||||
"iter h |ft-f0| |ft-f0-h*J0*dx| Order\n",
|
||||
"---------------------------------------------------------\n",
|
||||
" 0 1.00e-01 1.733e+05 7.391e+02 nan\n",
|
||||
" 1 1.00e-02 1.734e+04 6.892e+01 1.030\n",
|
||||
" 2 1.00e-03 1.734e+03 6.896e+00 1.000\n",
|
||||
" 3 1.00e-04 1.734e+02 6.896e-01 1.000\n",
|
||||
" 4 1.00e-05 1.734e+01 6.897e-02 1.000\n",
|
||||
" 5 1.00e-06 1.734e+00 6.897e-03 1.000\n",
|
||||
"*********************************************************\n",
|
||||
"<<<<<<<<<<<<<<<<<<<<<<<<< FAIL! >>>>>>>>>>>>>>>>>>>>>>>>>\n",
|
||||
"*********************************************************\n",
|
||||
"You break it, you fix it.\n",
|
||||
"\n",
|
||||
"==================== checkDerivative ====================\n",
|
||||
"iter h |ft-f0| |ft-f0-h*J0*dx| Order\n",
|
||||
"---------------------------------------------------------\n",
|
||||
" 0 1.00e-01 3.075e+05 3.072e+03 nan\n",
|
||||
" 1 1.00e-02 3.057e+04 1.202e+02 1.407\n",
|
||||
" 2 1.00e-03 3.055e+03 1.169e+01 1.012\n",
|
||||
" 3 1.00e-04 3.055e+02 1.168e+00 1.000\n",
|
||||
" 4 1.00e-05 3.055e+01 1.168e-01 1.000\n",
|
||||
" 5 1.00e-06 3.055e+00 1.168e-02 1.000\n",
|
||||
"*********************************************************\n",
|
||||
"<<<<<<<<<<<<<<<<<<<<<<<<< FAIL! >>>>>>>>>>>>>>>>>>>>>>>>>\n",
|
||||
"*********************************************************\n",
|
||||
"Did you put your clever trousers on today?\n",
|
||||
"\n",
|
||||
"==================== checkDerivative ====================\n",
|
||||
"iter h |ft-f0| |ft-f0-h*J0*dx| Order\n",
|
||||
"---------------------------------------------------------\n",
|
||||
" 0 1.00e-01 4.726e+05 1.575e+04 nan\n",
|
||||
" 1 1.00e-02 4.791e+04 1.217e+03 1.112\n",
|
||||
" 2 1.00e-03 4.797e+03 1.216e+02 1.001\n",
|
||||
" 3 1.00e-04 4.798e+02 1.216e+01 1.000\n",
|
||||
" 4 1.00e-05 4.798e+01 1.216e+00 1.000\n",
|
||||
" 5 1.00e-06 4.798e+00 1.216e-01 1.000\n",
|
||||
"*********************************************************\n",
|
||||
"<<<<<<<<<<<<<<<<<<<<<<<<< FAIL! >>>>>>>>>>>>>>>>>>>>>>>>>\n",
|
||||
"*********************************************************\n",
|
||||
"Coffee break?\n",
|
||||
"\n",
|
||||
"==================== checkDerivative ====================\n",
|
||||
"iter h |ft-f0| |ft-f0-h*J0*dx| Order\n",
|
||||
"---------------------------------------------------------\n",
|
||||
" 0 1.00e-01 4.323e+05 2.689e+04 nan\n",
|
||||
" 1 1.00e-02 4.326e+04 4.536e+02 1.773\n",
|
||||
" 2 1.00e-03 4.326e+03 3.671e+01 1.092\n",
|
||||
" 3 1.00e-04 4.327e+02 3.660e+00 1.001\n",
|
||||
" 4 1.00e-05 4.327e+01 3.660e-01 1.000\n",
|
||||
" 5 1.00e-06 4.327e+00 3.660e-02 1.000\n",
|
||||
"*********************************************************\n",
|
||||
"<<<<<<<<<<<<<<<<<<<<<<<<< FAIL! >>>>>>>>>>>>>>>>>>>>>>>>>\n",
|
||||
"*********************************************************\n",
|
||||
"No gold star for you.\n",
|
||||
"\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"def testDeriv_ZijN(prb):\n",
|
||||
" src = prb.survey.srcList[0]\n",
|
||||
" rx = src.rxList[0]\n",
|
||||
" rx.locs = np.array([[0.,0.,0,],[1.,1.,1.]])\n",
|
||||
" u0x = np.random.randn(survey.mesh.nE)+np.random.randn(prb.mesh.nE)*1j\n",
|
||||
" u0y = np.random.randn(survey.mesh.nE)+np.random.randn(prb.mesh.nE)*1j\n",
|
||||
"# u0 = np.vstack((simpeg.mkvc(u0x,2),simpeg.mkvc(u0y,2)))\n",
|
||||
" u0 = np.r_[u0x, u0y]\n",
|
||||
" f0 = prb.fieldsPair(prb.mesh,prb.survey)\n",
|
||||
" # u0 = np.hstack((simpeg.mkvc(u0_px,2),simpeg.mkvc(u0_py,2)))\n",
|
||||
" f0[src,'e_pxSolution'] = u0[:len(u0)/2]#u0x\n",
|
||||
" f0[src,'e_pySolution'] = u0[len(u0)/2::]#u0y\n",
|
||||
" def getZijNcomp(rx,f,v=None):\n",
|
||||
" if rx.locs.ndim == 3:\n",
|
||||
" eFLocs = rx.locs[:,:,0]\n",
|
||||
" bFLocs = rx.locs[:,:,1]\n",
|
||||
" else:\n",
|
||||
" eFLocs = rx.locs\n",
|
||||
" bFLocs = rx.locs\n",
|
||||
" # Get the projection\n",
|
||||
" Pex = prb.mesh.getInterpolationMat(eFLocs,'Ex')\n",
|
||||
" Pey = prb.mesh.getInterpolationMat(eFLocs,'Ey')\n",
|
||||
" Pbx = prb.mesh.getInterpolationMat(bFLocs,'Fx')\n",
|
||||
" Pby = prb.mesh.getInterpolationMat(bFLocs,'Fy')\n",
|
||||
" # Get the fields at location\n",
|
||||
" # px: x-polaration and py: y-polaration.\n",
|
||||
" ex_px = Pex*f[src,'e_px']\n",
|
||||
" ey_px = Pey*f[src,'e_px']\n",
|
||||
" ex_py = Pex*f[src,'e_py']\n",
|
||||
" ey_py = Pey*f[src,'e_py']\n",
|
||||
" hx_px = Pbx*f[src,'b_px']/mu_0\n",
|
||||
" hy_px = Pby*f[src,'b_px']/mu_0\n",
|
||||
" hx_py = Pbx*f[src,'b_py']/mu_0\n",
|
||||
" hy_py = Pby*f[src,'b_py']/mu_0\n",
|
||||
" # Derivatives as lambda functions\n",
|
||||
" # The size of the deratives should be nD,nU\n",
|
||||
" ex_px_u = lambda vec: Pex*f._e_pxDeriv_u(src,vec)\n",
|
||||
" ey_px_u = lambda vec: Pey*f._e_pxDeriv_u(src,vec)\n",
|
||||
" ex_py_u = lambda vec: Pex*f._e_pyDeriv_u(src,vec)\n",
|
||||
" ey_py_u = lambda vec: Pey*f._e_pyDeriv_u(src,vec)\n",
|
||||
" # NOTE: Think b_p?Deriv_u should return a 2*nF size matrix\n",
|
||||
" hx_px_u = lambda vec: Pbx*f._b_pxDeriv_u(src,vec)/mu_0\n",
|
||||
" hy_px_u = lambda vec: Pby*f._b_pxDeriv_u(src,vec)/mu_0\n",
|
||||
" hx_py_u = lambda vec: Pbx*f._b_pyDeriv_u(src,vec)/mu_0\n",
|
||||
" hy_py_u = lambda vec: Pby*f._b_pyDeriv_u(src,vec)/mu_0\n",
|
||||
"\n",
|
||||
" # Update the input vector\n",
|
||||
" # Define the components of the derivative# Calculate components\n",
|
||||
" # Update the input vector\n",
|
||||
" sdiag = lambda t: simpeg.Utils.sdiag(simpeg.mkvc(t,2))\n",
|
||||
" # Define the components of the derivative\n",
|
||||
" if 'zxx' in rx.rxType:\n",
|
||||
" if v is not None:\n",
|
||||
" return sdiag(hy_py)*ex_px_u(v) + sdiag(ex_px)*hy_py_u(v) - sdiag(ex_py)*hy_px_u(v) - sdiag(hy_px)*ex_py_u(v)\n",
|
||||
" return ( sdiag(ex_px)*hy_py - sdiag(ex_py)*hy_px)\n",
|
||||
" elif 'zxy' in rx.rxType:\n",
|
||||
" if v is not None:\n",
|
||||
" return -sdiag(hx_py)*ex_px_u(v) - sdiag(ex_px)*hx_py_u(v) + sdiag(ex_py)*hx_px_u(v) + sdiag(hx_px)*ex_py_u(v)\n",
|
||||
" return (-sdiag(ex_px)*hx_py + sdiag(ex_py)*hx_px)\n",
|
||||
" elif 'zyx' in rx.rxType:\n",
|
||||
" if v is not None:\n",
|
||||
" return ey_px_u(v)*hy_py + ey_px*hy_py_u(v) - ey_py*hy_px_u(v) - ey_py_u(v)*hy_px\n",
|
||||
" return ( ey_px*hy_py - ey_py*hy_px)\n",
|
||||
" elif 'zyy' in rx.rxType:\n",
|
||||
" if v is not None:\n",
|
||||
" return -ey_px_u(v)*hx_py - ey_px*hx_py_u(v) + ey_py*hx_px_u(v) + ey_py_u(v)*hx_px\n",
|
||||
" return (-ey_px*hx_py + ey_py*hx_px)\n",
|
||||
" \n",
|
||||
" def fun(u):\n",
|
||||
" f = prb.fieldsPair(survey.mesh,prb.survey)\n",
|
||||
" f[src,'e_pxSolution'] = u[:len(u)/2]\n",
|
||||
" f[src,'e_pySolution'] = u[len(u)/2::]\n",
|
||||
" return getZijNcomp(rx,f), lambda t: getZijNcomp(rx,f0,t)\n",
|
||||
" simpeg.Tests.checkDerivative(fun,u0,num=6,plotIt=False)\n",
|
||||
" \n",
|
||||
"surveyxx, problemxx = t3Dmt.setupSimpegMTfwd_eForm_ps(t3Dmt.random(.01),comp='zxxr',singleFreq=.001)\n",
|
||||
"testDeriv_ZijN(problemxx)\n",
|
||||
"surveyxy, problemxy = t3Dmt.setupSimpegMTfwd_eForm_ps(t3Dmt.random(.01),comp='zxyr',singleFreq=.0001)\n",
|
||||
"testDeriv_ZijN(problemxy)\n",
|
||||
"surveyyx, problemyx = t3Dmt.setupSimpegMTfwd_eForm_ps(t3Dmt.random(.01),comp='zyxr',singleFreq=.0001)\n",
|
||||
"testDeriv_ZijN(problemyx)\n",
|
||||
"surveyyy, problemyy = t3Dmt.setupSimpegMTfwd_eForm_ps(t3Dmt.random(.01),comp='zyyr',singleFreq=.0001)\n",
|
||||
"testDeriv_ZijN(problemyy)"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -316,6 +629,127 @@
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# %debug"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 12,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# dA_dm = problem.getADeriv_m(src.freq,u_src,w)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 13,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# dA_dm.shape"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 14,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# a, b, sig, d, e = t3Dmt.halfSpace(.01)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 40,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"==================== checkDerivative ====================\n",
|
||||
"iter h |ft-f0| |ft-f0-h*J0*dx| Order\n",
|
||||
"---------------------------------------------------------\n",
|
||||
" 0 1.00e-01 5.050e-06 1.906e-09 nan\n",
|
||||
" 1 1.00e-02 5.049e-07 1.949e-11 1.990\n",
|
||||
" 2 1.00e-03 5.048e-08 1.969e-13 1.996\n",
|
||||
" 3 1.00e-04 5.048e-09 1.377e-15 2.155\n",
|
||||
" 4 1.00e-05 5.048e-10 1.073e-15 0.108\n",
|
||||
" 5 1.00e-06 5.049e-11 4.179e-15 -0.590\n",
|
||||
" 6 1.00e-07 5.059e-12 1.079e-14 -0.412\n",
|
||||
" 7 1.00e-08 5.092e-13 4.411e-15 0.388\n",
|
||||
"========================= PASS! =========================\n",
|
||||
"Yay passed!\n",
|
||||
"\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"def testDeriv_ProjFields(prb):\n",
|
||||
" # Initate things for the derivs Test\n",
|
||||
" src = prb.survey.srcList[0]\n",
|
||||
" rx = src.rxList[0]\n",
|
||||
"# rx.locs = np.array([[0.,0.,0,],[1.,1.,1.]])\n",
|
||||
"# u0x = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j\n",
|
||||
"# u0y = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j\n",
|
||||
"# u0 = np.vstack((simpeg.mkvc(u0x,2),simpeg.mkvc(u0y,2)))\n",
|
||||
"# u0 = np.r_[u0x, u0y]\n",
|
||||
"# f0 = prb.fieldsPair(survey.mesh,survey)\n",
|
||||
" # u0 = np.hstack((simpeg.mkvc(u0_px,2),simpeg.mkvc(u0_py,2)))\n",
|
||||
"# f0[src,'e_pxSolution'] = u0[:len(u0)/2]#u0x\n",
|
||||
"# f0[src,'e_pySolution'] = u0[len(u0)/2::]#u0y\n",
|
||||
"\n",
|
||||
"# prb.mapping = simpeg.Maps.ExpMap(prb.mesh)\n",
|
||||
"# if True:\n",
|
||||
"# x0 = x0 + np.random.randn(prb.mesh.nC)*cond*1e-1 \n",
|
||||
" survey = prb.survey\n",
|
||||
" src = survey.srcList[0]\n",
|
||||
" \n",
|
||||
" f0 = prb.fields(prb.curModel)\n",
|
||||
" u0 = np.r_[f0[src,'e_pxSolution'],f0[src,'e_pySolution']]\n",
|
||||
"# def fun(x):\n",
|
||||
"# prb.curModel = x\n",
|
||||
" # f0[src,'b_1d'] = -1/(1j*simpegem.Utils.EMUtils.omega(src.freq))*m1d.nodalGrad*u0\n",
|
||||
" # Run a test\n",
|
||||
" def fun(u):\n",
|
||||
" f = prb.fieldsPair(survey.mesh,survey)\n",
|
||||
" f[src,'e_pxSolution'] = u[:len(u)/2]\n",
|
||||
" f[src,'e_pySolution'] = u[len(u)/2::]\n",
|
||||
" # f[src,'b_1d'] = -(m1d.nodalGrad*u)/(1j*simpegem.Utils.EMUtils.omega(src.freq))\n",
|
||||
" return rx.projectFields(src,survey.mesh,f), lambda t: rx.projectFieldsDeriv(src,survey.mesh,f0,t)\n",
|
||||
" simpeg.Tests.checkDerivative(fun,u0,num=8,plotIt=False)\n",
|
||||
"survey, problem = t3Dmt.setupSimpegMTfwd_eForm_ps(t3Dmt.halfSpace(1),comp='zxyr',singleFreq=.1)\n",
|
||||
"testDeriv_ProjFields(problem)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 36,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# %debug"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 17,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"src = survey.srcList[0]\n",
|
||||
"rx = src.rxList[0]"
|
||||
@@ -323,7 +757,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 12,
|
||||
"execution_count": 18,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -334,16 +768,30 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"execution_count": 30,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"array([ 19.30432593, 8.59301637, 7.75725432, ..., 2.18228759,\n",
|
||||
" 16.72231401, 7.67029076])"
|
||||
]
|
||||
},
|
||||
"execution_count": 30,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"problem.curModel.sigma"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 14,
|
||||
"execution_count": 19,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -361,7 +809,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 17,
|
||||
"execution_count": 20,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -369,47 +817,89 @@
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"(2, 1)"
|
||||
"(20536, 1)"
|
||||
]
|
||||
},
|
||||
"execution_count": 17,
|
||||
"execution_count": 20,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"rx.projectFieldsDeriv(src,survey.mesh,f0,u0).shape"
|
||||
"f0._e_px(f0[src,'e_pxSolution'],[src]).shape"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 16,
|
||||
"execution_count": 21,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"(0, 1)"
|
||||
]
|
||||
},
|
||||
"execution_count": 21,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"f0._e_pxDeriv_u([src],u0)[len(u0)/2::].shape"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 22,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"simpeg.Utils.spzeros?"
|
||||
"simpeg.Utils.spzeros"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"execution_count": 23,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"array([[ 0.00031425],\n",
|
||||
" [ 0.00030875]])"
|
||||
]
|
||||
},
|
||||
"execution_count": 23,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"rx.projectFields(src,survey.mesh,f0)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"execution_count": 24,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stderr",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"ERROR: No traceback has been produced, nothing to debug.\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"%debug"
|
||||
]
|
||||
@@ -429,6 +919,18 @@
|
||||
"display_name": "Python 2",
|
||||
"language": "python",
|
||||
"name": "python2"
|
||||
},
|
||||
"language_info": {
|
||||
"codemirror_mode": {
|
||||
"name": "ipython",
|
||||
"version": 2
|
||||
},
|
||||
"file_extension": ".py",
|
||||
"mimetype": "text/x-python",
|
||||
"name": "python",
|
||||
"nbconvert_exporter": "python",
|
||||
"pygments_lexer": "ipython2",
|
||||
"version": "2.7.10"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
|
||||
Reference in New Issue
Block a user