mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-29 05:35:13 +08:00
JTvec is working but not converging.
This commit is contained in:
GudniRos
@@ -237,11 +237,11 @@
|
||||
"==================== checkDerivative ====================\n",
|
||||
"iter h |ft-f0| |ft-f0-h*J0*dx| Order\n",
|
||||
"---------------------------------------------------------\n",
|
||||
" 0 1.00e-01 3.893e-05 5.627e-08 nan\n",
|
||||
" 1 1.00e-02 3.898e-06 5.632e-10 2.000\n",
|
||||
" 2 1.00e-03 3.898e-07 5.633e-12 2.000\n",
|
||||
" 3 1.00e-04 3.898e-08 5.633e-14 2.000\n",
|
||||
" 4 1.00e-05 3.898e-09 5.637e-16 2.000\n",
|
||||
" 0 1.00e-01 9.779e-05 3.219e-06 nan\n",
|
||||
" 1 1.00e-02 1.007e-05 3.145e-08 2.010\n",
|
||||
" 2 1.00e-03 1.010e-06 3.137e-10 2.001\n",
|
||||
" 3 1.00e-04 1.010e-07 3.136e-12 2.000\n",
|
||||
" 4 1.00e-05 1.010e-08 3.136e-14 2.000\n",
|
||||
"========================= PASS! =========================\n",
|
||||
"Yay passed!\n",
|
||||
"\n"
|
||||
@@ -277,7 +277,7 @@
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"array(-0.005683351098616839)"
|
||||
"array([ 0.00017028])"
|
||||
]
|
||||
},
|
||||
"execution_count": 8,
|
||||
@@ -299,7 +299,7 @@
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"array([[ 0.00374478]])"
|
||||
"array([[ 0.0014539]])"
|
||||
]
|
||||
},
|
||||
"execution_count": 9,
|
||||
@@ -338,13 +338,13 @@
|
||||
"---------------------------------------------------------\n",
|
||||
"Project at freq: 1.000e+02\n",
|
||||
"Project at freq: 1.000e+02\n",
|
||||
" 0 1.00e-01 7.283e-06 9.326e-09 nan\n",
|
||||
" 0 1.00e-01 8.466e-06 5.557e-08 nan\n",
|
||||
"Project at freq: 1.000e+02\n",
|
||||
" 1 1.00e-02 7.290e-07 9.338e-11 1.999\n",
|
||||
" 1 1.00e-02 8.416e-07 5.519e-10 2.003\n",
|
||||
"Project at freq: 1.000e+02\n",
|
||||
" 2 1.00e-03 7.290e-08 9.339e-13 2.000\n",
|
||||
" 2 1.00e-03 8.411e-08 5.515e-12 2.000\n",
|
||||
"========================= PASS! =========================\n",
|
||||
"Testing is important.\n",
|
||||
"Not just a pretty face Gudni\n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
@@ -439,22 +439,19 @@
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Adjoint E formulation - e\n"
|
||||
"Adjoint E formulation - e\n",
|
||||
"-4.65700174631e-05 3.27055690471e-05 -7.92755865102e-05 1e-08 False\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"ename": "IndexError",
|
||||
"evalue": "tuple index out of range",
|
||||
"output_type": "error",
|
||||
"traceback": [
|
||||
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
|
||||
"\u001b[1;31mIndexError\u001b[0m Traceback (most recent call last)",
|
||||
"\u001b[1;32m<ipython-input-14-250b6ff497c4>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0madjointTest\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'E'\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'e'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
|
||||
"\u001b[1;32m<ipython-input-13-3fac9635b176>\u001b[0m in \u001b[0;36madjointTest\u001b[1;34m(fdemType, comp)\u001b[0m\n\u001b[0;32m 16\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 17\u001b[0m \u001b[0mvJw\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mv\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mproblem\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mJvec\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mm\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mw\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mu\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 18\u001b[1;33m \u001b[0mwJtv\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mw\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mproblem\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mJtvec\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mm\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mv\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mu\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 19\u001b[0m \u001b[0mtol\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmax\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mTOL\u001b[0m\u001b[1;33m*\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m10\u001b[0m\u001b[1;33m**\u001b[0m\u001b[0mint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mlog10\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mabs\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mvJw\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mFLR\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 20\u001b[0m \u001b[1;32mprint\u001b[0m \u001b[0mvJw\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mwJtv\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mvJw\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mwJtv\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtol\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mabs\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mvJw\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mwJtv\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m<\u001b[0m \u001b[0mtol\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpegmt/simpegMT/BaseMT.pyc\u001b[0m in \u001b[0;36mJtvec\u001b[1;34m(self, m, v, f)\u001b[0m\n\u001b[0;32m 102\u001b[0m \u001b[0mPTv\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mrx\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mprojectFieldsDeriv\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmesh\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mf\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mv\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mrx\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0madjoint\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mTrue\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;31m# wrt u, need possibility wrt m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 103\u001b[0m \u001b[1;31m# Get the\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 104\u001b[1;33m \u001b[0mdA_duIT\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mATinv\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mPTv\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 105\u001b[0m \u001b[0mdA_dmT\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetADeriv_m\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfreq\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mu_src\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdA_duIT\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0madjoint\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mTrue\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 106\u001b[0m \u001b[0mdRHS_dmT\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetRHSDeriv_m\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdA_duIT\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0madjoint\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mTrue\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpeg/SimPEG/Utils/SolverUtils.pyc\u001b[0m in \u001b[0;36m__mul__\u001b[1;34m(self, b)\u001b[0m\n\u001b[0;32m 35\u001b[0m \u001b[1;32mraise\u001b[0m \u001b[0mTypeError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'Can only multiply by a numpy array.'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 36\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 37\u001b[1;33m \u001b[1;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mb\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;36m1\u001b[0m \u001b[1;32mor\u001b[0m \u001b[0mb\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 38\u001b[0m \u001b[0mb\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mb\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mflatten\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 39\u001b[0m \u001b[1;31m# Just one RHS\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[1;31mIndexError\u001b[0m: tuple index out of range"
|
||||
]
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"False"
|
||||
]
|
||||
},
|
||||
"execution_count": 14,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
@@ -463,34 +460,14 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"execution_count": 17,
|
||||
"metadata": {
|
||||
"collapsed": false,
|
||||
"scrolled": true
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"> \u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpeg/SimPEG/Utils/SolverUtils.py\u001b[0m(37)\u001b[0;36m__mul__\u001b[1;34m()\u001b[0m\n",
|
||||
"\u001b[1;32m 36 \u001b[1;33m\u001b[1;33m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0m\u001b[1;32m---> 37 \u001b[1;33m \u001b[1;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mb\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;36m1\u001b[0m \u001b[1;32mor\u001b[0m \u001b[0mb\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0m\u001b[1;32m 38 \u001b[1;33m \u001b[0mb\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mb\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mflatten\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0m\n",
|
||||
"ipdb> u\n",
|
||||
"> \u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpegmt/simpegMT/BaseMT.py\u001b[0m(104)\u001b[0;36mJtvec\u001b[1;34m()\u001b[0m\n",
|
||||
"\u001b[1;32m 103 \u001b[1;33m \u001b[1;31m# Get the\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0m\u001b[1;32m--> 104 \u001b[1;33m \u001b[0mdA_duIT\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mATinv\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mPTv\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0m\u001b[1;32m 105 \u001b[1;33m \u001b[0mdA_dmT\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetADeriv_m\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfreq\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mu_src\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdA_duIT\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0madjoint\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mTrue\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[0m\n",
|
||||
"ipdb> PTv\n",
|
||||
"array(-6.9926641874700175e-06)\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"%debug"
|
||||
"# %debug"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -500,10 +477,7 @@
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"v = np.random.randn(1,1) + 1j* np.random.randn(1,1)\n",
|
||||
"rx.projectFieldsDeriv(src, problem.mesh, fields, v, adjoint=True)"
|
||||
]
|
||||
"source": []
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
@@ -512,21 +486,7 @@
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"dMf_dsig = problem.mesh.getFaceInnerProductDeriv(problem.curModel.sigma)(u0) * problem.curModel.sigmaDeriv\n",
|
||||
"dsig_dm = self.curModel.sigmaDeriv\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"problem.mesh.getFaceInnerProductDeriv(problem.curModel.sigma)(u0)"
|
||||
]
|
||||
"source": []
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
|
||||
+1
-1
@@ -103,7 +103,7 @@ class BaseMTProblem(BaseFDEMProblem):
|
||||
# Get the
|
||||
dA_duIT = ATinv * PTv
|
||||
dA_dmT = self.getADeriv_m(freq, u_src, dA_duIT, adjoint=True)
|
||||
dRHS_dmT = self.getRHSDeriv_m(src, dA_duIT, adjoint=True)
|
||||
dRHS_dmT = self.getRHSDeriv_m(freq, dA_duIT, adjoint=True)
|
||||
# Make du_dmT
|
||||
if dRHS_dmT is None:
|
||||
du_dmT = - dA_dmT
|
||||
|
||||
@@ -146,9 +146,10 @@ class RxMT(Survey.BaseRx):
|
||||
Pbx = mesh.getInterpolationMat(self.locs,'Ex')
|
||||
# ex = Pex*mkvc(f[src,'e_1d'],2)
|
||||
# bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0
|
||||
dP_de = Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))*(Pex*v)
|
||||
dP_db = - Utils.sdiag(Pex*mkvc(f[src,'e_1d'],2))*(Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)).T*Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)))*(Pbx*f._bDeriv_u(src,v)/mu_0)
|
||||
PDeriv_complex = np.sum((dP_de,dP_db))
|
||||
dP_de = mkvc(Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))*(Pex*v),2)
|
||||
dP_db = mkvc(- Utils.sdiag(Pex*mkvc(f[src,'e_1d'],2))*(Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)).T*Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)))*(Pbx*f._bDeriv_u(src,v)/mu_0),2)
|
||||
PDeriv_complex = np.sum(np.hstack((dP_de,dP_db)),1)
|
||||
# raise Exception('Debug error')
|
||||
# elif self.projType is 'Z2D
|
||||
elif self.projType is 'Z3D':
|
||||
raise NotImplementedError('Has not be implement for full impedance tensor')
|
||||
@@ -160,9 +161,10 @@ class RxMT(Survey.BaseRx):
|
||||
# ex = Pex*mkvc(f[src,'e_1d'],2)
|
||||
# bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0
|
||||
dP_deTv = mkvc(Pex.T*Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)).T*v,2)
|
||||
db_duv = Pbx.T/mu_0*Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))*(Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))).T*Utils.sdiag(Pex*mkvc(f[src,'e_1d'],2)).T
|
||||
dP_dbTv = -mkvc(f._bDeriv_u(src,db_duv,adjoint=True)*v,2)
|
||||
PDeriv_complex = np.sum((dP_deTv,dP_dbTv))
|
||||
db_duv = Pbx.T/mu_0*Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))*(Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))).T*Utils.sdiag(Pex*mkvc(f[src,'e_1d'],2)).T*v
|
||||
dP_dbTv = -mkvc(f._bDeriv_u(src,db_duv,adjoint=True),2)
|
||||
PDeriv_complex = np.sum(np.hstack((dP_deTv,dP_dbTv)),1)
|
||||
# raise Exception('Debug error')
|
||||
elif self.projType is 'Z3D':
|
||||
raise NotImplementedError('must be real or imag')
|
||||
Pv = np.array(getattr(PDeriv_complex, real_or_imag))
|
||||
|
||||
Reference in New Issue
Block a user