mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-08 14:24:21 +08:00
Jvec adjoint test not working, but all other derivatives adjoint tests are working.
This commit is contained in:
GudniRos
@@ -237,13 +237,13 @@
|
||||
"==================== checkDerivative ====================\n",
|
||||
"iter h |ft-f0| |ft-f0-h*J0*dx| Order\n",
|
||||
"---------------------------------------------------------\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",
|
||||
" 0 1.00e-01 1.076e-05 3.966e-07 nan\n",
|
||||
" 1 1.00e-02 1.040e-06 3.875e-09 2.010\n",
|
||||
" 2 1.00e-03 1.037e-07 3.866e-11 2.001\n",
|
||||
" 3 1.00e-04 1.037e-08 3.865e-13 2.000\n",
|
||||
" 4 1.00e-05 1.037e-09 3.865e-15 2.000\n",
|
||||
"========================= PASS! =========================\n",
|
||||
"Yay passed!\n",
|
||||
"The test be workin!\n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
@@ -277,7 +277,7 @@
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"array([ 0.00017028])"
|
||||
"array([ -2.07544500e-05])"
|
||||
]
|
||||
},
|
||||
"execution_count": 8,
|
||||
@@ -299,7 +299,7 @@
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"array([[ 0.0014539]])"
|
||||
"array([[-0.00014517]])"
|
||||
]
|
||||
},
|
||||
"execution_count": 9,
|
||||
@@ -338,13 +338,15 @@
|
||||
"---------------------------------------------------------\n",
|
||||
"Project at freq: 1.000e+02\n",
|
||||
"Project at freq: 1.000e+02\n",
|
||||
" 0 1.00e-01 8.466e-06 5.557e-08 nan\n",
|
||||
" 0 1.00e-01 1.142e-06 2.072e-08 nan\n",
|
||||
"Project at freq: 1.000e+02\n",
|
||||
" 1 1.00e-02 8.416e-07 5.519e-10 2.003\n",
|
||||
" 1 1.00e-02 1.161e-07 2.039e-10 2.007\n",
|
||||
"Project at freq: 1.000e+02\n",
|
||||
" 2 1.00e-03 8.411e-08 5.515e-12 2.000\n",
|
||||
" 2 1.00e-03 1.162e-08 2.036e-12 2.001\n",
|
||||
"Project at freq: 1.000e+02\n",
|
||||
" 3 1.00e-04 1.163e-09 2.036e-14 2.000\n",
|
||||
"========================= PASS! =========================\n",
|
||||
"Not just a pretty face Gudni\n",
|
||||
"Testing is important.\n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
@@ -374,7 +376,7 @@
|
||||
"def fun(x):\n",
|
||||
" \n",
|
||||
" return survey.dpred(x), lambda x: problem.Jvec(m0, x)\n",
|
||||
"simpeg.Tests.checkDerivative(fun, m0, num=3, plotIt=False)"
|
||||
"simpeg.Tests.checkDerivative(fun, m0, num=4, plotIt=False)"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -403,12 +405,195 @@
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Adjoint e formulation - projectFieldsDeriv\n",
|
||||
"-9.0191451349e-06 -9.0191451349e-06 3.38813178902e-21 1e-09 True\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"True"
|
||||
]
|
||||
},
|
||||
"execution_count": 13,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# Run a test\n",
|
||||
"TOL = 1e-4\n",
|
||||
"FLR = 1e-20\n",
|
||||
"\n",
|
||||
"def projectFieldsAdjointTest(fdemType, comp):\n",
|
||||
" print 'Adjoint %s formulation - %s' % (fdemType, comp)\n",
|
||||
"\n",
|
||||
" m = np.log(np.ones(problem.mesh.nC)*0.01)\n",
|
||||
" if True:\n",
|
||||
" m = m + np.random.randn(problem.mesh.nC)*0.01*1e-1 \n",
|
||||
"\n",
|
||||
" u = problem.fields(m)\n",
|
||||
" v = np.random.randn(1)#+np.random.randn(1)*1j\n",
|
||||
" # print prb.PropMap.PropModel.nP\n",
|
||||
" w = np.random.randn(m1d.nN)#+np.random.randn(m1d.nN)*1j\n",
|
||||
"\n",
|
||||
" vJw = v.dot(rx.projectFieldsDeriv(src,m1d,f0,w))\n",
|
||||
" wJtv = w.dot(rx.projectFieldsDeriv(src,m1d,f0,v,adjoint=True))\n",
|
||||
" tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR]) \n",
|
||||
" print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol\n",
|
||||
" return np.abs(vJw - wJtv) < tol\n",
|
||||
"projectFieldsAdjointTest('e','projectFieldsDeriv')"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 14,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Adjoint test e formulation - getADeriv_m\n",
|
||||
"(508851.474801-63591.6147417j) (508851.474801-63591.6147417j) (-1.7462298274e-10-5.82076609135e-11j) 10.0 True\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"True"
|
||||
]
|
||||
},
|
||||
"execution_count": 14,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# Run a test\n",
|
||||
"TOL = 1e-4\n",
|
||||
"FLR = 1e-20\n",
|
||||
"\n",
|
||||
"def getADeriv_mAdjointTest():\n",
|
||||
" print 'Adjoint test e formulation - getADeriv_m' \n",
|
||||
"\n",
|
||||
" m = np.log(np.ones(problem.mesh.nC)*0.01)\n",
|
||||
" if True:\n",
|
||||
" m = m + np.random.randn(problem.mesh.nC)*0.01*1e-1 \n",
|
||||
"\n",
|
||||
" u = problem.fields(m)\n",
|
||||
" v = np.random.randn(m1d.nN)#+np.random.randn(1)*1j\n",
|
||||
" # print prb.PropMap.PropModel.nP\n",
|
||||
" w = np.random.randn(m1d.nC)#+np.random.randn(m1d.nN)*1j\n",
|
||||
"\n",
|
||||
" vJw = v.dot(problem.getADeriv_m(freq,u0,w))\n",
|
||||
" wJtv = w.dot(problem.getADeriv_m(freq,u0,v,adjoint=True))\n",
|
||||
" tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR]) \n",
|
||||
" print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol\n",
|
||||
" return np.abs(vJw - wJtv) < tol\n",
|
||||
"getADeriv_mAdjointTest()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 15,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Adjoint test e formulation - getRHSDeriv_m\n",
|
||||
"(-29427.4754714-17701.4445852j) (-29427.4754714-17701.4445852j) (2.91038304567e-11+0j) 1.0 True\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"True"
|
||||
]
|
||||
},
|
||||
"execution_count": 15,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# Run a test\n",
|
||||
"TOL = 1e-4\n",
|
||||
"FLR = 1e-20\n",
|
||||
"\n",
|
||||
"def getRHSDeriv_mAdjointTest():\n",
|
||||
" print 'Adjoint test e formulation - getRHSDeriv_m'\n",
|
||||
"\n",
|
||||
" m = np.log(np.ones(problem.mesh.nC)*0.01)\n",
|
||||
" if True:\n",
|
||||
" m = m + np.random.randn(problem.mesh.nC)*0.01*1e-1 \n",
|
||||
"\n",
|
||||
" u = problem.fields(m)\n",
|
||||
" v = np.random.randn(m1d.nN)#+np.random.randn(1)*1j\n",
|
||||
" # print prb.PropMap.PropModel.nP\n",
|
||||
" w = np.random.randn(m1d.nC)#+np.random.randn(m1d.nN)*1j\n",
|
||||
"\n",
|
||||
" vJw = v.dot(problem.getRHSDeriv_m(freq,w))\n",
|
||||
" wJtv = w.dot(problem.getRHSDeriv_m(freq,v,adjoint=True))\n",
|
||||
" tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR]) \n",
|
||||
" print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol\n",
|
||||
" return np.abs(vJw - wJtv) < tol\n",
|
||||
"getRHSDeriv_mAdjointTest( )"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 16,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"2\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"simpeg.mkvc(np.random.randn(survey.nD)+np.random.randn(survey.nD)*1j,2)\n",
|
||||
"\n",
|
||||
"print survey.nD"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 17,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"TOL = 1e-4\n",
|
||||
"FLR = 1e-20\n",
|
||||
"\n",
|
||||
"def adjointTest(fdemType, comp):\n",
|
||||
"def JvecAdjointTest(fdemType, comp):\n",
|
||||
" print 'Adjoint %s formulation - %s' % (fdemType, comp)\n",
|
||||
"\n",
|
||||
" m = np.log(np.ones(problem.mesh.nC)*0.01)\n",
|
||||
@@ -430,7 +615,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 14,
|
||||
"execution_count": 18,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -440,7 +625,7 @@
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Adjoint E formulation - e\n",
|
||||
"-4.65700174631e-05 3.27055690471e-05 -7.92755865102e-05 1e-08 False\n"
|
||||
"-1.7867035849e-05 1.42750293876e-05 -3.21420652366e-05 1e-08 False\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -449,25 +634,63 @@
|
||||
"False"
|
||||
]
|
||||
},
|
||||
"execution_count": 14,
|
||||
"execution_count": 18,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"adjointTest('E','e')"
|
||||
"JvecAdjointTest('E','e')"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 17,
|
||||
"execution_count": 19,
|
||||
"metadata": {
|
||||
"collapsed": false,
|
||||
"scrolled": true
|
||||
},
|
||||
"outputs": [],
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stderr",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"ERROR: No traceback has been produced, nothing to debug.\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# %debug"
|
||||
"%debug"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 20,
|
||||
"metadata": {
|
||||
"collapsed": false,
|
||||
"scrolled": true
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Project at freq: 1.000e+02\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"array([ 9.80523303e-06, -1.98372645e-03])"
|
||||
]
|
||||
},
|
||||
"execution_count": 20,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"survey.dpred(sigma)"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -481,20 +704,23 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"execution_count": 21,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"<180x180 sparse matrix of type '<type 'numpy.complex128'>'\n",
|
||||
"\twith 180 stored elements in Compressed Sparse Row format>"
|
||||
]
|
||||
},
|
||||
"execution_count": 21,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"problem.mesh.getEdgeInnerProductDeriv(problem.curModel.sigma)(u0[1::])"
|
||||
]
|
||||
@@ -514,6 +740,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.9"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
|
||||
+3
-3
@@ -106,15 +106,15 @@ class BaseMTProblem(BaseFDEMProblem):
|
||||
dRHS_dmT = self.getRHSDeriv_m(freq, dA_duIT, adjoint=True)
|
||||
# Make du_dmT
|
||||
if dRHS_dmT is None:
|
||||
du_dmT = - dA_dmT
|
||||
du_dmT = -dA_dmT
|
||||
else:
|
||||
du_dmT = -dA_dmT + dRHS_dmT
|
||||
# Select the correct component
|
||||
real_or_imag = rx.projComp
|
||||
if real_or_imag == 'real':
|
||||
Jtv += du_dmT.real
|
||||
Jtv += du_dmT.real
|
||||
elif real_or_imag == 'imag':
|
||||
Jtv += - du_dmT.real
|
||||
Jtv += -du_dmT.real
|
||||
else:
|
||||
raise Exception('Must be real or imag')
|
||||
|
||||
|
||||
@@ -161,8 +161,8 @@ 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*v
|
||||
dP_dbTv = -mkvc(f._bDeriv_u(src,db_duv,adjoint=True),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*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':
|
||||
|
||||
Reference in New Issue
Block a user