wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-10 03:30:57 +08:00
Code Issues Packages Projects Releases Wiki Activity

Jvec adjoint test is working for MT1D primary/secondary formulation.

Browse Source
This commit is contained in:
GudniRos
2015-06-25 10:36:10 -07:00
parent d233e40a95
commit 4a39602ca4
3 changed files with 89 additions and 82 deletions
Download Patch File Download Diff File Expand all files Collapse all files
notebooks/Derivative test MT1D.ipynb
+74 -65
View File
@@ -237,13 +237,13 @@
"==================== checkDerivative ====================\n",
"iter h |ft-f0| |ft-f0-h*J0*dx| Order\n",
"---------------------------------------------------------\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",
" 0 1.00e-01 1.094e-05 3.087e-08 nan\n",
" 1 1.00e-02 1.097e-06 3.093e-10 1.999\n",
" 2 1.00e-03 1.097e-07 3.093e-12 2.000\n",
" 3 1.00e-04 1.097e-08 3.093e-14 2.000\n",
" 4 1.00e-05 1.097e-09 3.095e-16 2.000\n",
"========================= PASS! =========================\n",
"The test be workin!\n",
"You are awesome.\n",
"\n"
]
},
@@ -277,7 +277,7 @@
{
"data": {
"text/plain": [
"array([ -2.07544500e-05])"
"array([ 0.00127341])"
]
},
"execution_count": 8,
@@ -299,7 +299,7 @@
{
"data": {
"text/plain": [
"array([[-0.00014517]])"
"array([[ 0.00173178]])"
]
},
"execution_count": 9,
@@ -338,15 +338,15 @@
"---------------------------------------------------------\n",
"Project at freq: 1.000e+02\n",
"Project at freq: 1.000e+02\n",
" 0 1.00e-01 1.142e-06 2.072e-08 nan\n",
" 0 1.00e-01 2.646e-06 2.277e-08 nan\n",
"Project at freq: 1.000e+02\n",
" 1 1.00e-02 1.161e-07 2.039e-10 2.007\n",
" 1 1.00e-02 2.629e-07 2.246e-10 2.006\n",
"Project at freq: 1.000e+02\n",
" 2 1.00e-03 1.162e-08 2.036e-12 2.001\n",
" 2 1.00e-03 2.627e-08 2.243e-12 2.001\n",
"Project at freq: 1.000e+02\n",
" 3 1.00e-04 1.163e-09 2.036e-14 2.000\n",
" 3 1.00e-04 2.627e-09 2.242e-14 2.000\n",
"========================= PASS! =========================\n",
"Testing is important.\n",
"Yay passed!\n",
"\n"
]
},
@@ -411,7 +411,7 @@
"output_type": "stream",
"text": [
"Adjoint e formulation - projectFieldsDeriv\n",
"-9.0191451349e-06 -9.0191451349e-06 3.38813178902e-21 1e-09 True\n"
"-0.000991536643367 -0.000991536643367 2.16840434497e-19 1e-07 True\n"
]
},
{
@@ -440,10 +440,10 @@
" 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",
" 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",
" wJtv = w.dot(rx.projectFieldsDeriv(src,m1d,f0,v,adjoint=True)).real\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",
@@ -456,13 +456,32 @@
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"ERROR: No traceback has been produced, nothing to debug.\n"
]
}
],
"source": [
"%debug"
]
},
{
"cell_type": "code",
"execution_count": 15,
"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"
"(340193.379835-398622.996348j) (340193.379835-398622.996348j) (-2.32830643654e-10-3.49245965481e-10j) 10.0 True\n"
]
},
{
@@ -471,7 +490,7 @@
"True"
]
},
"execution_count": 14,
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
@@ -512,7 +531,7 @@
},
{
"cell_type": "code",
"execution_count": 15,
"execution_count": 16,
"metadata": {
"collapsed": false
},
@@ -522,7 +541,7 @@
"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"
"(-12351.1349263+433.4655562j) (-12351.1349263+433.4655562j) (-9.09494701773e-12-1.22781784739e-11j) 1.0 True\n"
]
},
{
@@ -531,7 +550,7 @@
"True"
]
},
"execution_count": 15,
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
@@ -563,7 +582,7 @@
},
{
"cell_type": "code",
"execution_count": 16,
"execution_count": 17,
"metadata": {
"collapsed": false
},
@@ -584,7 +603,7 @@
},
{
"cell_type": "code",
"execution_count": 17,
"execution_count": 18,
"metadata": {
"collapsed": false
},
@@ -593,8 +612,8 @@
"TOL = 1e-4\n",
"FLR = 1e-20\n",
"\n",
"def JvecAdjointTest(fdemType, comp):\n",
" print 'Adjoint %s formulation - %s' % (fdemType, comp)\n",
"def JvecAdjointTest():\n",
" print 'Adjoint e formulation - Jvec' \n",
"\n",
" m = np.log(np.ones(problem.mesh.nC)*0.01)\n",
" if True:\n",
@@ -615,7 +634,7 @@
},
{
"cell_type": "code",
"execution_count": 18,
"execution_count": 19,
"metadata": {
"collapsed": false
},
@@ -624,28 +643,28 @@
"name": "stdout",
"output_type": "stream",
"text": [
"Adjoint E formulation - e\n",
"-1.7867035849e-05 1.42750293876e-05 -3.21420652366e-05 1e-08 False\n"
"Adjoint e formulation - Jvec\n",
"-3.61480355369e-05 -3.61480355369e-05 -2.71050543121e-20 1e-08 True\n"
]
},
{
"data": {
"text/plain": [
"False"
"True"
]
},
"execution_count": 18,
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"JvecAdjointTest('E','e')"
"JvecAdjointTest()"
]
},
{
"cell_type": "code",
"execution_count": 19,
"execution_count": 20,
"metadata": {
"collapsed": false,
"scrolled": true
@@ -665,7 +684,7 @@
},
{
"cell_type": "code",
"execution_count": 20,
"execution_count": 21,
"metadata": {
"collapsed": false,
"scrolled": true
@@ -684,7 +703,7 @@
"array([ 9.80523303e-06, -1.98372645e-03])"
]
},
"execution_count": 20,
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
@@ -695,32 +714,34 @@
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 21,
"execution_count": 22,
"metadata": {
"collapsed": false
},
"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"
"ename": "NameError",
"evalue": "name 'r' is not defined",
"output_type": "error",
"traceback": [
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[1;31mNameError\u001b[0m Traceback (most recent call last)",
"\u001b[1;32m<ipython-input-22-72cfd272ace1>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mr\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[1;31mNameError\u001b[0m: name 'r' is not defined"
]
}
],
"source": [
"r"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"problem.mesh.getEdgeInnerProductDeriv(problem.curModel.sigma)(u0[1::])"
]
@@ -740,18 +761,6 @@
"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,
simpegMT/BaseMT.py
-10
View File
@@ -61,16 +61,6 @@ class BaseMTProblem(BaseFDEMProblem):
du_dm = dA_duI * ( - dA_dm + dRHS_dm )
# Calculate the projection derivatives
for rx in src.rxList:
# Get the stacked derivative
# df_duFun = getattr(f, '_fDeriv_u', None)
# df_dmFun = getattr(f, '_fDeriv_m', None)
# df_dm = df_dmFun(src,v,adjoint=False)
# if df_dm is None:
# fDeriv_m = df_duFun(src, du_dm, adjoint=False)
# else:
# fDeriv_m = df_duFun(src, du_dm, adjoint=False) + df_dm
# Not needed for now. Since PDeriv does this currently.
# Get the projection derivative
PDeriv = lambda v: rx.projectFieldsDeriv(src, self.mesh, f, v) # wrt u, also have wrt m
Jv[src, rx] = PDeriv(du_dm)
simpegMT/SurveyMT.py
+15 -7
View File
@@ -149,12 +149,14 @@ class RxMT(Survey.BaseRx):
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 'Z2D':
raise NotImplementedError('Has not be implement for 2D impedance tensor')
elif self.projType is 'Z3D':
raise NotImplementedError('Has not be implement for full impedance tensor')
raise NotImplementedError('Has not be implement for full 3D impedance tensor')
# Extract the real number for the real/imag components.
Pv = np.array(getattr(PDeriv_complex, real_or_imag))
elif adjoint:
# Note: The v vector is real and the return should be complex
if self.projType is 'Z1D':
Pex = mesh.getInterpolationMat(self.locs,'Fx')
Pbx = mesh.getInterpolationMat(self.locs,'Ex')
@@ -163,11 +165,17 @@ class RxMT(Survey.BaseRx):
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)
PDeriv_complex = np.sum(np.hstack((dP_deTv,dP_dbTv)),1)
# raise Exception('Debug error')
PDeriv_real = np.sum(np.hstack((dP_deTv,dP_dbTv)),1)
elif self.projType is 'Z2D':
raise NotImplementedError('Has not be implement for 2D impedance tensor')
elif self.projType is 'Z3D':
raise NotImplementedError('must be real or imag')
Pv = np.array(getattr(PDeriv_complex, real_or_imag))
raise NotImplementedError('Has not be implement for full 3D impedance tensor')
# Extract the data
if real_or_imag == 'imag':
Pv = 1j*PDeriv_real
elif real_or_imag == 'real':
Pv = PDeriv_real.astype(complex)
return Pv