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

Added a Derivative test notebook

Browse Source
This commit is contained in:
GudniRos
2015-11-02 11:23:06 -08:00
parent 42bc4404ba
commit 960cb0a3e5
2 changed files with 453 additions and 17 deletions
Download Patch File Download Diff File Expand all files Collapse all files
notebooks/MT3D derivatives test.ipynb
+436
View File
@@ -0,0 +1,436 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import SimPEG as simpeg\n",
"import simpegEM as simpegem, simpegMT as simpegmt\n",
"from SimPEG.Utils import meshTensor\n",
"import numpy as np\n",
"import simpegMT.Tests.test_Problem3D_againstAnalytic as t3Dmt"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"sigmaHalf = 0.01\n",
"survey, problem = t3Dmt.setupSimpegMTfwd_eForm_ps(t3Dmt.halfSpace(sigmaHalf),comp='zxyi',singleFreq=1)\n",
"if False:\n",
" fields = problem.fields(problem.curModel.sigma)\n",
" data = survey.dpred(problem.curModel.sigma)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"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.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",
"========================= PASS! =========================\n",
"Happy little convergence test!\n",
"\n"
]
},
{
"data": {
"text/plain": [
"True"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"def testDeriv_dA_dm(prb,cond):\n",
" TOL = 1e-4\n",
" FLR = 1e-20\n",
" x0 = np.log(np.ones(prb.mesh.nC)*cond)\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",
" freq = src.freq\n",
" v1 = np.random.randn(prb.mesh.nE,1)\n",
" v2 = np.random.randn(prb.mesh.nE,1)\n",
" v = np.hstack(( simpeg.mkvc(v1,2), simpeg.mkvc(v2,2)))\n",
" u_0 = prb.fieldsPair(prb.mesh,prb.survey)\n",
" u_0[src,'e_pxSolution'] = v1\n",
" u_0[src,'e_pySolution'] = v2\n",
" def fun(x):\n",
" prb.curModel = x\n",
" A = prb.getA(freq) #\n",
"# return simpeg.mkvc(A*v1)+simpeg.mkvc(A*v2), lambda t: simpeg.mkvc(prb.getADeriv_m(freq, u_0[src], t))\n",
" return A*v, lambda t: (prb.getADeriv_m(freq, u_0[src], t))\n",
" return simpeg.Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=FLR)\n",
"testDeriv_dA_dm(problem,0.1)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"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.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",
"========================= PASS! =========================\n",
"Yay passed!\n",
"\n"
]
},
{
"data": {
"text/plain": [
"True"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"def testDeriv_dRHS_dm(prb,cond):\n",
" TOL = 1e-4\n",
" FLR = 1e-20\n",
" x0 = np.log(np.ones(prb.mesh.nC)*cond)\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",
" u0 = prb.fields(x0)\n",
" freq = src.freq\n",
" A = prb.getA(freq) #\n",
" A_I = prb.Solver(A, **prb.solverOpts)\n",
"\n",
" ftype = prb._fieldType + 'Solution'\n",
" u0_src = u0[src, ftype]\n",
" v = np.random.randn(prb.mesh.nE,1)\n",
" \n",
" def fun(x):\n",
" prb.curModel = x\n",
" return simpeg.mkvc(np.sum(prb.getRHS(freq))), lambda x: simpeg.mkvc(prb.getRHSDeriv_m(freq, x))\n",
"# return simpeg.mkvc(prb.fields(x)[src,prb._fieldType + 'Solution']), lambda x: simpeg.mkvc(prb.getADeriv_m(freq, u0_src, x))\n",
" return simpeg.Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=FLR)\n",
"testDeriv_dA_dm(problem,0.1)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"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 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",
"========================= PASS! =========================\n",
"That was easy!\n",
"\n"
]
}
],
"source": [
"def testDeriv_S_e(prb,cond):\n",
" # Initate things for the derivs Test\n",
" TOL = 1e-4\n",
" FLR = 1e-20\n",
" \n",
" src = prb.survey.srcList[0]\n",
" rx = src.rxList[0]\n",
"\n",
" x0 = np.log(np.ones(prb.mesh.nC)*cond)\n",
"# prb.mapping = simpeg.Maps.ExpMap(prb.mesh)\n",
" if True:\n",
" x0 = x0 + np.random.randn(prb.mesh.nC)*cond*1e-1 \n",
" def fun(x):\n",
" prb.curModel = x\n",
" return src.S_e(prb), lambda t: src.S_eDeriv_m(prb,t)\n",
" simpeg.Tests.checkDerivative(fun,x0,num=3,plotIt=False)\n",
"survey, problem = t3Dmt.setupSimpegMTfwd_eForm_ps(t3Dmt.random(1),comp='All',singleFreq=1)\n",
"testDeriv_S_e(problem,0.1)"
]
},
{
"cell_type": "code",
"execution_count": 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",
"output_type": "stream",
"text": [
"==================== 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",
"\n"
]
}
],
"source": [
"def testDeriv_ProjFields(prb):\n",
" # Initate things for the derivs Test\n",
" src = 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 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.random(.1),comp='zyxr',singleFreq=.01)\n",
"testDeriv_ProjFields(problem)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"src = survey.srcList[0]\n",
"rx = src.rxList[0]"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"rx.locs = np.array([[0, 0, 0],[1,1,1]])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"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",
"src = survey.srcList[0]\n",
"rx = src.rxList[0]\n",
"f0 = problem.fieldsPair(survey.mesh,survey)\n",
"u0 = np.vstack((simpeg.mkvc(u0x,2),simpeg.mkvc(u0y,2)))\n",
"f0[src,'e_pxSolution'] = u0x\n",
"f0[src,'e_pySolution'] = u0y"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"(2, 1)"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"rx.projectFieldsDeriv(src,survey.mesh,f0,u0).shape"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"simpeg.Utils.spzeros?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"rx.projectFields(src,survey.mesh,f0)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%debug"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
simpegMT/SurveyMT.py
+17 -17
View File
@@ -188,31 +188,31 @@ class RxMT(Survey.BaseRx):
Pby = mesh.getInterpolationMat(bFLocs,'Fy')
# Get the fields at location
# px: x-polaration and py: y-polaration.
ex_px = Utils.sdiag(mkvc(Pex*f[src,'e_px'],2))
ey_px = Utils.sdiag(mkvc(Pey*f[src,'e_px'],2))
ex_py = Utils.sdiag(mkvc(Pex*f[src,'e_py'],2))
ey_py = Utils.sdiag(mkvc(Pey*f[src,'e_py'],2))
hx_px = Utils.sdiag(mkvc(Pbx*f[src,'b_px']/mu_0,2))
hy_px = Utils.sdiag(mkvc(Pby*f[src,'b_px']/mu_0,2))
hx_py = Utils.sdiag(mkvc(Pbx*f[src,'b_py']/mu_0,2))
hy_py = Utils.sdiag(mkvc(Pby*f[src,'b_py']/mu_0,2))
ex_px = Pex*f[src,'e_px']
ey_px = Pey*f[src,'e_px']
ex_py = Pex*f[src,'e_py']
ey_py = Pey*f[src,'e_py']
hx_px = Pbx*f[src,'b_px']/mu_0
hy_px = Pby*f[src,'b_px']/mu_0
hx_py = Pbx*f[src,'b_py']/mu_0
hy_py = Pby*f[src,'b_py']/mu_0
# Derivatives as lambda functions
spPe = Utils.spzeros(self.nD,mesh.nE)
spPb = Utils.spzeros(self.nD,mesh.nF)
# The size of the diratives should be nD,nU
ex_px_u = lambda vec: Utils.sdiag(mkvc(sp.hstack((Pex,spPe))*f._e_pxDeriv_u(src,vec),2))
ey_px_u = lambda vec: Utils.sdiag(mkvc(sp.hstack((Pey,spPe))*f._e_pxDeriv_u(src,vec),2))
ex_py_u = lambda vec: Utils.sdiag(mkvc(sp.hstack((spPe,Pex))*f._e_pyDeriv_u(src,vec),2))
ey_py_u = lambda vec: Utils.sdiag(mkvc(sp.hstack((spPe,Pey))*f._e_pyDeriv_u(src,vec),2))
ex_px_u = lambda vec: sp.hstack((Pex,spPe))*f._e_pxDeriv_u(src,vec)
ey_px_u = lambda vec: sp.hstack((Pey,spPe))*f._e_pxDeriv_u(src,vec)
ex_py_u = lambda vec: sp.hstack((spPe,Pex))*f._e_pyDeriv_u(src,vec)
ey_py_u = lambda vec: sp.hstack((spPe,Pey))*f._e_pyDeriv_u(src,vec)
# NOTE: Think b_p?Deriv_u should return a 2*nF size matrix
hx_px_u = lambda vec: Utils.sdiag(mkvc(sp.hstack((Pbx,spPb))*f._b_pxDeriv_u(src,vec)/mu_0,2))
hy_px_u = lambda vec: Utils.sdiag(mkvc(sp.hstack((Pby,spPb))*f._b_pxDeriv_u(src,vec)/mu_0,2))
hx_py_u = lambda vec: Utils.sdiag(mkvc(sp.hstack((spPb,Pbx))*f._b_pyDeriv_u(src,vec)/mu_0,2))
hy_py_u = lambda vec: Utils.sdiag(mkvc(sp.hstack((spPb,Pby))*f._b_pyDeriv_u(src,vec)/mu_0,2))
hx_px_u = lambda vec: sp.hstack((Pbx,spPb))*f._b_pxDeriv_u(src,vec)/mu_0
hy_px_u = lambda vec: sp.hstack((Pby,spPb))*f._b_pxDeriv_u(src,vec)/mu_0
hx_py_u = lambda vec: sp.hstack((spPb,Pbx))*f._b_pyDeriv_u(src,vec)/mu_0
hy_py_u = lambda vec: sp.hstack((spPb,Pby))*f._b_pyDeriv_u(src,vec)/mu_0
# Update the input vector
# v = mkvc(v,2) # Make v into a column vector
# Define the components of the derivative
Hd = Utils.sdiag(mkvc(1./(hx_px*hy_py - hx_py*hy_px).data,2))
Hd = 1./(hx_px*hy_py - hx_py*hy_px)
Hd_uV = hx_px_u(v)*hy_py + hx_px*hy_py_u(v) - hx_py*hy_px_u(v) - hx_py_u(v)*hy_px
# Calculate components
if 'zxx' in self.rxType: