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

projection with sdiag for all the elements. Not sure if the is correct, but archiving the work

Browse Source
This commit is contained in:
GudniRos
2015-11-01 13:05:24 -08:00
parent 9d6a1dcac6
commit 42bc4404ba
6 changed files with 360 additions and 118 deletions
Download Patch File Download Diff File Expand all files Collapse all files
notebooks/Derivative test MT1D.ipynb
+314 -79
View File
@@ -84,12 +84,13 @@
" srcList.append(simpegmt.SurveyMT.srcMT_polxy_1DhomotD(rxList,freq))\n",
"else:\n",
" for freq in freqs:\n",
" srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freq,sigma))\n",
" srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freq))\n",
"# Make the survey\n",
"survey = simpegmt.SurveyMT.SurveyMT(srcList)\n",
"\n",
"# Set the problem\n",
"problem = simpegmt.ProblemMT1D.eForm_psField(m1d)\n",
"problem.sigmaPrimary = sigma\n",
"problem.pair(survey)\n",
"\n",
"# Get the fields\n",
@@ -101,12 +102,22 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": []
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"ERROR: No traceback has been produced, nothing to debug.\n"
]
}
],
"source": [
"%debug\n"
]
},
{
"cell_type": "markdown",
@@ -136,7 +147,7 @@
},
{
"cell_type": "code",
"execution_count": 4,
"execution_count": 5,
"metadata": {
"collapsed": true
},
@@ -168,7 +179,7 @@
},
{
"cell_type": "code",
"execution_count": 5,
"execution_count": 6,
"metadata": {
"collapsed": true
},
@@ -192,7 +203,7 @@
},
{
"cell_type": "code",
"execution_count": 6,
"execution_count": 7,
"metadata": {
"collapsed": false
},
@@ -216,7 +227,7 @@
},
{
"cell_type": "code",
"execution_count": 7,
"execution_count": 8,
"metadata": {
"collapsed": false,
"scrolled": true
@@ -229,13 +240,13 @@
"==================== checkDerivative ====================\n",
"iter h |ft-f0| |ft-f0-h*J0*dx| Order\n",
"---------------------------------------------------------\n",
" 0 1.00e-01 1.884e-05 1.227e-07 nan\n",
" 1 1.00e-02 1.873e-06 1.265e-09 1.987\n",
" 2 1.00e-03 1.872e-07 1.269e-11 1.999\n",
" 3 1.00e-04 1.872e-08 1.269e-13 2.000\n",
" 4 1.00e-05 1.872e-09 1.269e-15 2.000\n",
" 0 1.00e-01 3.089e-04 2.730e-05 nan\n",
" 1 1.00e-02 2.841e-05 2.582e-07 2.024\n",
" 2 1.00e-03 2.818e-06 2.568e-09 2.002\n",
" 3 1.00e-04 2.816e-07 2.567e-11 2.000\n",
" 4 1.00e-05 2.816e-08 2.567e-13 2.000\n",
"========================= PASS! =========================\n",
"That was easy!\n",
"The test be workin!\n",
"\n"
]
},
@@ -245,7 +256,7 @@
"True"
]
},
"execution_count": 7,
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
@@ -261,7 +272,7 @@
},
{
"cell_type": "code",
"execution_count": 8,
"execution_count": 9,
"metadata": {
"collapsed": false
},
@@ -269,10 +280,10 @@
{
"data": {
"text/plain": [
"array([ 0.00052762])"
"array([-0.00285083])"
]
},
"execution_count": 8,
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
@@ -283,7 +294,7 @@
},
{
"cell_type": "code",
"execution_count": 9,
"execution_count": 10,
"metadata": {
"collapsed": false
},
@@ -291,10 +302,10 @@
{
"data": {
"text/plain": [
"array([[ 0.00124017]])"
"array([[-0.00130618]])"
]
},
"execution_count": 9,
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
@@ -305,7 +316,7 @@
},
{
"cell_type": "code",
"execution_count": 10,
"execution_count": 11,
"metadata": {
"collapsed": false
},
@@ -316,7 +327,7 @@
},
{
"cell_type": "code",
"execution_count": 11,
"execution_count": 12,
"metadata": {
"collapsed": false
},
@@ -328,12 +339,12 @@
"==================== checkDerivative ====================\n",
"iter h |ft-f0| |ft-f0-h*J0*dx| Order\n",
"---------------------------------------------------------\n",
" 0 1.00e-01 4.417e-08 4.873e-09 nan\n",
" 1 1.00e-02 4.132e-09 4.832e-11 2.004\n",
" 2 1.00e-03 4.105e-10 4.828e-13 2.000\n",
" 3 1.00e-04 4.103e-11 4.827e-15 2.000\n",
" 0 1.00e-01 1.537e-05 1.005e-07 nan\n",
" 1 1.00e-02 1.541e-06 1.008e-09 1.999\n",
" 2 1.00e-03 1.541e-07 1.008e-11 2.000\n",
" 3 1.00e-04 1.541e-08 1.008e-13 2.000\n",
"========================= PASS! =========================\n",
"You get a gold star!\n",
"The test be workin!\n",
"\n"
]
},
@@ -343,7 +354,7 @@
"True"
]
},
"execution_count": 11,
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
@@ -368,7 +379,26 @@
},
{
"cell_type": "code",
"execution_count": 12,
"execution_count": 13,
"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": 14,
"metadata": {
"collapsed": true
},
@@ -388,7 +418,7 @@
},
{
"cell_type": "code",
"execution_count": 13,
"execution_count": 15,
"metadata": {
"collapsed": false
},
@@ -398,7 +428,7 @@
"output_type": "stream",
"text": [
"Adjoint e formulation - projectFieldsDeriv\n",
"-2.26989762698e-05 -2.26989762698e-05 0.0 1e-08 True\n"
"0.000177264079929 0.000177264079929 3.52365706058e-19 1e-07 True\n"
]
},
{
@@ -407,7 +437,7 @@
"True"
]
},
"execution_count": 13,
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
@@ -439,7 +469,7 @@
},
{
"cell_type": "code",
"execution_count": 14,
"execution_count": 21,
"metadata": {
"collapsed": false
},
@@ -448,19 +478,25 @@
"name": "stdout",
"output_type": "stream",
"text": [
"Adjoint test e formulation - getADeriv_m\n",
"(-1977540.36505+2093781.70221j) (-1977540.36505+2093781.70221j) (-1.86264514923e-09+2.79396772385e-09j) 100.0 True\n"
"Adjoint test e formulation - getADeriv_m\n"
]
},
{
"data": {
"text/plain": [
"True"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
"ename": "KeyError",
"evalue": "'Source does not have a uid: e_1dSolution'",
"output_type": "error",
"traceback": [
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[1;31mKeyError\u001b[0m Traceback (most recent call last)",
"\u001b[1;32m<ipython-input-21-4c7e7058b8e9>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\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[0;32m 21\u001b[0m \u001b[1;32mreturn\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---> 22\u001b[1;33m \u001b[0mgetADeriv_mAdjointTest\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[1;32m<ipython-input-21-4c7e7058b8e9>\u001b[0m in \u001b[0;36mgetADeriv_mAdjointTest\u001b[1;34m()\u001b[0m\n\u001b[0;32m 15\u001b[0m \u001b[0mw\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mrandom\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mrandn\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mm1d\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mnC\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;31m#+np.random.randn(m1d.nN)*1j\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 16\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 17\u001b[1;33m \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[0mgetADeriv_m\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfreq\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mf0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mw\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 18\u001b[0m \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[0mgetADeriv_m\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfreq\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mf0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mv\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[1;33m\u001b[0m\u001b[0m\n\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[1;32m/media/gudni/ExtraDrive1/Codes/python/simpegmt/simpegMT/ProblemMT1D/Problems.py\u001b[0m in \u001b[0;36mgetADeriv_m\u001b[1;34m(self, freq, u, v, adjoint)\u001b[0m\n\u001b[0;32m 67\u001b[0m \u001b[0mMeMui\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[0mgetEdgeInnerProduct\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m1.0\u001b[0m\u001b[1;33m/\u001b[0m\u001b[0mmu_0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 68\u001b[0m \u001b[1;31m#\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 69\u001b[1;33m \u001b[0mu_src\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mu\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'e_1dSolution'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 70\u001b[0m \u001b[0mdMf_dsig\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[0mgetFaceInnerProductDeriv\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcurModel\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msigma\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mu_src\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcurModel\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msigmaDeriv\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 71\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0madjoint\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/Fields.pyc\u001b[0m in \u001b[0;36m__getitem__\u001b[1;34m(self, key)\u001b[0m\n\u001b[0;32m 150\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 151\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m__getitem__\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mkey\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 152\u001b[1;33m \u001b[0mind\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mname\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_indexAndNameFromKey\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m'get'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 153\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0mname\u001b[0m \u001b[1;32mis\u001b[0m \u001b[0mNone\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 154\u001b[0m \u001b[0mout\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m{\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/Fields.pyc\u001b[0m in \u001b[0;36m_indexAndNameFromKey\u001b[1;34m(self, key, accessType)\u001b[0m\n\u001b[0;32m 131\u001b[0m \u001b[0msrcTestList\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mname\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mkey\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 132\u001b[0m \u001b[0mname\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_nameIndex\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mname\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maccessType\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 133\u001b[1;33m \u001b[0mind\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_srcIndex\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msrcTestList\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 134\u001b[0m \u001b[1;32mreturn\u001b[0m \u001b[0mind\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mname\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 135\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpeg/SimPEG/Fields.pyc\u001b[0m in \u001b[0;36m_srcIndex\u001b[1;34m(self, srcTestList)\u001b[0m\n\u001b[0;32m 100\u001b[0m \u001b[0mind\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0msrcTestList\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 101\u001b[0m \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 102\u001b[1;33m \u001b[0mind\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msurvey\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetSourceIndex\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msrcTestList\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 103\u001b[0m \u001b[1;32mreturn\u001b[0m \u001b[0mind\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 104\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpeg/SimPEG/Survey.pyc\u001b[0m in \u001b[0;36mgetSourceIndex\u001b[1;34m(self, sources)\u001b[0m\n\u001b[0;32m 339\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0msrc\u001b[0m \u001b[1;32min\u001b[0m \u001b[0msources\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 340\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0mgetattr\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'uid'\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mNone\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mis\u001b[0m \u001b[0mNone\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 341\u001b[1;33m \u001b[1;32mraise\u001b[0m \u001b[0mKeyError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'Source does not have a uid: %s'\u001b[0m\u001b[1;33m%\u001b[0m\u001b[0mstr\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msrc\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 342\u001b[0m \u001b[0minds\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mmap\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;32mlambda\u001b[0m \u001b[0msrc\u001b[0m\u001b[1;33m:\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_sourceOrder\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mget\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0muid\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mNone\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0msources\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 343\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0mNone\u001b[0m \u001b[1;32min\u001b[0m \u001b[0minds\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;31mKeyError\u001b[0m: 'Source does not have a uid: e_1dSolution'"
]
}
],
"source": [
@@ -480,8 +516,8 @@
" # 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",
" vJw = v.dot(problem.getADeriv_m(freq,f0,w))\n",
" wJtv = w.dot(problem.getADeriv_m(freq,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",
@@ -490,16 +526,226 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 17,
"metadata": {
"collapsed": true
"collapsed": false
},
"outputs": [],
"source": []
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"> \u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpegmt/simpegMT/ProblemMT1D/Problems.py\u001b[0m(69)\u001b[0;36mgetADeriv_m\u001b[1;34m()\u001b[0m\n",
"\u001b[1;32m 68 \u001b[1;33m \u001b[1;31m#\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[0m\u001b[1;32m---> 69 \u001b[1;33m \u001b[0mu_src\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mu\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'e_1dSolution'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[0m\u001b[1;32m 70 \u001b[1;33m \u001b[0mdMf_dsig\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[0mgetFaceInnerProductDeriv\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcurModel\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msigma\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mu_src\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcurModel\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msigmaDeriv\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[0m\n",
"ipdb> u\n",
"> \u001b[1;32m<ipython-input-16-0d25a7e16094>\u001b[0m(17)\u001b[0;36mgetADeriv_mAdjointTest\u001b[1;34m()\u001b[0m\n",
"\u001b[1;32m 16 \u001b[1;33m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[0m\u001b[1;32m---> 17 \u001b[1;33m \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[0mgetADeriv_m\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfreq\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mu0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mw\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[0m\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[0mgetADeriv_m\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfreq\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mu0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mv\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[1;33m\u001b[0m\u001b[0m\n",
"\u001b[0m\n",
"ipdb> d\n",
"> \u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpegmt/simpegMT/ProblemMT1D/Problems.py\u001b[0m(69)\u001b[0;36mgetADeriv_m\u001b[1;34m()\u001b[0m\n",
"\u001b[1;32m 68 \u001b[1;33m \u001b[1;31m#\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[0m\u001b[1;32m---> 69 \u001b[1;33m \u001b[0mu_src\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mu\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'e_1dSolution'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[0m\u001b[1;32m 70 \u001b[1;33m \u001b[0mdMf_dsig\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[0mgetFaceInnerProductDeriv\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcurModel\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msigma\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mu_src\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcurModel\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msigmaDeriv\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[0m\n",
"ipdb> p u\n",
"array([ 7.09826705e-01 -3.42993400e-01j,\n",
" 2.06495724e-02 -1.13236591e+00j,\n",
" 7.19940065e-01 +6.18299116e-01j,\n",
" -3.76573540e-01 +9.09687373e-01j,\n",
" -5.96166731e-01 +2.23309643e-01j,\n",
" 3.30533731e-01 +7.98068064e-01j,\n",
" -1.41356795e+00 -1.74137927e-01j,\n",
" 1.51554175e+00 -2.16363375e-01j,\n",
" 1.59469627e+00 -1.87971566e+00j,\n",
" -5.54637698e-01 -1.01936458e+00j,\n",
" 1.01869865e+00 +8.33443055e-01j,\n",
" -4.92616671e-01 -1.97036224e+00j,\n",
" 1.00389102e+00 -1.52567261e-01j,\n",
" 7.90254779e-01 -5.68065046e-01j,\n",
" 9.17329964e-01 -1.18879097e+00j,\n",
" -1.38317362e+00 -2.17789273e+00j,\n",
" 2.58324368e-01 -1.51003598e+00j,\n",
" -2.04803634e+00 -1.44299943e-01j,\n",
" 3.33011285e+00 +1.52905384e+00j,\n",
" 2.91429972e-01 -1.91849967e+00j,\n",
" -1.10273784e+00 -8.65794290e-01j,\n",
" -2.80103156e-01 -1.42653192e+00j,\n",
" -9.81404885e-01 +9.06845144e-02j,\n",
" 6.37579045e-01 +3.84148006e-01j,\n",
" -6.83597988e-01 +6.91187965e-01j,\n",
" -8.98575981e-01 +3.80347476e-02j,\n",
" 1.02647319e+00 -1.13876879e-01j,\n",
" 8.28006255e-01 -9.95903406e-01j,\n",
" 3.53931083e-01 -2.80188697e-01j,\n",
" -1.09339690e+00 +1.43646692e+00j,\n",
" 2.59276626e+00 +4.09662561e-01j,\n",
" -5.97873868e-01 -5.23852065e-01j,\n",
" 8.87054287e-01 -1.69872744e+00j,\n",
" 2.29767481e-01 -5.88143928e-01j,\n",
" -6.97610315e-01 +5.68810360e-01j,\n",
" -4.47049016e-01 +3.29354463e-01j,\n",
" -3.31498311e-01 -8.27120276e-01j,\n",
" 1.57689998e-01 +6.68378872e-01j,\n",
" 9.47018631e-02 +2.60180446e-01j,\n",
" -4.59177257e-01 -7.66896049e-01j,\n",
" -1.16791204e+00 +4.20585804e-01j,\n",
" 2.72405147e-01 +3.48079787e-01j,\n",
" -2.45880514e-02 +1.90226778e-01j,\n",
" -7.03428363e-01 +1.20531568e+00j,\n",
" 1.60362217e+00 +6.88472259e-01j,\n",
" -5.49144335e-01 +6.96208292e-01j,\n",
" -9.83200346e-01 -1.15726575e+00j,\n",
" 1.22826524e+00 -8.58070855e-01j,\n",
" -5.99020439e-01 -2.19182124e+00j,\n",
" -5.69992871e-01 +1.68586900e+00j,\n",
" 9.26752644e-01 +1.05767829e+00j,\n",
" 8.64247444e-01 +1.51282356e-01j,\n",
" -6.24547862e-01 -1.54517937e+00j,\n",
" -1.72043526e+00 +1.91659834e+00j,\n",
" -1.84501951e+00 -2.40165623e-01j,\n",
" 1.15210121e+00 -1.19627222e+00j,\n",
" -3.93036395e-01 -1.94918396e+00j,\n",
" -1.18266273e+00 -5.24652722e-01j,\n",
" 2.67461628e-01 -8.07279138e-01j,\n",
" -1.56363641e+00 +1.59257615e-03j,\n",
" -3.17211883e-01 -3.23379690e+00j,\n",
" -1.08658939e-01 -8.66079039e-01j,\n",
" -9.30350013e-01 -1.22668219e+00j,\n",
" -1.07263137e+00 -5.56847000e-02j,\n",
" 1.27878013e+00 -1.34173537e+00j,\n",
" 7.91441492e-01 +6.42111414e-01j,\n",
" -2.60121575e+00 -1.08140886e+00j,\n",
" -1.15027808e+00 -4.47110142e-01j,\n",
" -5.69747800e-01 +2.34065916e-01j,\n",
" 1.52155089e+00 -5.80340775e-01j,\n",
" 2.50052750e-01 -6.58793179e-01j,\n",
" 2.96944135e-01 +1.27455206e+00j,\n",
" 9.18939766e-01 -4.58122761e-01j,\n",
" 4.90015531e-01 -8.86288179e-01j,\n",
" -1.65155083e+00 -6.57733967e-01j,\n",
" 9.72148169e-01 -1.19446820e+00j,\n",
" 2.08365350e+00 -1.53585595e-01j,\n",
" -4.11101072e-01 +1.01418536e+00j,\n",
" 9.72580562e-01 -2.20470840e-01j,\n",
" 1.72922147e-01 -6.86360934e-01j,\n",
" 3.80255034e-01 +1.96862619e-01j,\n",
" -1.53628572e+00 +5.31946011e-01j,\n",
" -4.65939167e-01 +9.11725611e-01j,\n",
" -1.06390365e+00 -3.07658816e-02j,\n",
" -2.40734296e+00 -7.60530882e-01j,\n",
" 1.28242272e+00 +3.32085973e+00j,\n",
" 3.17004751e+00 +1.62654811e+00j,\n",
" -8.49758476e-01 +9.86658359e-01j,\n",
" 4.99578557e-01 +4.01392488e-01j,\n",
" -9.48695998e-01 +4.04322381e-01j,\n",
" 1.75371116e-01 +7.85458813e-01j,\n",
" 4.90232886e-01 -1.39280042e-01j,\n",
" -2.31430252e-01 -8.44065839e-01j,\n",
" -1.02319133e+00 +1.90560784e+00j,\n",
" 1.08940667e+00 +1.63046186e+00j,\n",
" -4.34233324e-01 +2.90439312e-01j,\n",
" 4.47622371e-01 +1.26507976e+00j,\n",
" -2.07317470e+00 -2.21558532e+00j,\n",
" -3.19553487e-01 -9.19210427e-01j,\n",
" 6.09870390e-01 +3.05692898e-01j,\n",
" 4.92081860e-01 -2.18326059e-01j,\n",
" 1.48544627e+00 -2.02054490e-01j,\n",
" -1.59612367e-01 +7.90853777e-02j,\n",
" 2.85854411e-01 +2.16567814e-01j,\n",
" 5.28063071e-01 +2.08838949e+00j,\n",
" 6.17580558e-01 +9.90748539e-01j,\n",
" -2.64539211e-01 +1.16598665e+00j,\n",
" -1.07357542e+00 -7.58002064e-01j,\n",
" 6.13399713e-02 +2.52338250e+00j,\n",
" -4.84180135e-01 -1.45046020e+00j,\n",
" -6.57519374e-01 -1.52719059e+00j,\n",
" -1.54949084e-01 -7.64796563e-01j,\n",
" -2.01172915e-02 +1.07377565e+00j,\n",
" 2.19840636e-01 +6.89622678e-01j,\n",
" -2.64578500e+00 +1.46933966e+00j,\n",
" 1.22168696e+00 +2.33458454e-01j,\n",
" 2.71303177e+00 +2.77858685e-01j,\n",
" 4.39426665e-01 +1.79924695e-01j,\n",
" -1.73858667e+00 -3.03473636e-01j,\n",
" 2.68736300e-01 +8.32359634e-01j,\n",
" 2.57517784e+00 +1.88097576e-01j,\n",
" 2.76775469e-01 -6.49891901e-01j,\n",
" -9.87790301e-01 +4.50328269e-01j,\n",
" 9.19033669e-01 +1.01453748e+00j,\n",
" -1.35706747e+00 +5.90240590e-01j,\n",
" 6.60118738e-01 +5.10361208e-01j,\n",
" -2.17652642e-01 -4.26632365e-01j,\n",
" -1.39163540e+00 +9.12739752e-01j,\n",
" -1.76241829e+00 -1.04835908e+00j,\n",
" -7.28279529e-01 +2.75961750e-01j,\n",
" -1.32043513e+00 -3.28790183e-01j,\n",
" 5.88227093e-01 +9.41865970e-01j,\n",
" 1.97695607e+00 +1.04535450e+00j,\n",
" -3.02599128e+00 +1.02031744e+00j,\n",
" -1.34308763e+00 -4.48222977e-01j,\n",
" -1.86687861e+00 +2.26767568e+00j,\n",
" 5.61472023e-01 -5.28652462e-01j,\n",
" -5.26747946e-01 +2.75545054e-01j,\n",
" 6.16116514e-01 +3.41765793e-01j,\n",
" 1.38709750e+00 +5.87537003e-01j,\n",
" -2.27923203e-01 +7.19317677e-01j,\n",
" -8.66404713e-01 -1.80932081e-01j,\n",
" -4.58146668e-01 +1.75955171e+00j,\n",
" -5.12526647e-01 -1.91939928e-01j,\n",
" -6.26532699e-01 -5.71301355e-01j,\n",
" -3.35566040e-01 -5.06173544e-01j,\n",
" -3.86124497e-01 -8.00560567e-01j,\n",
" 2.47252265e+00 +3.83928843e-01j,\n",
" 3.35561169e-01 +1.19309445e+00j,\n",
" -6.19075444e-01 -1.19021628e+00j,\n",
" -1.17335297e+00 +5.80500783e-01j,\n",
" -2.35111248e-01 -1.59785495e+00j,\n",
" -1.24819583e-01 -8.05281256e-01j,\n",
" -9.18606131e-01 -1.25312332e-01j,\n",
" 3.14400972e-01 +1.18650079e+00j,\n",
" 6.50630049e-01 -6.99853497e-01j,\n",
" 1.28269722e+00 -1.29832850e+00j,\n",
" 2.07880202e-01 +2.56385692e-02j,\n",
" -1.45089461e+00 -9.35093103e-01j,\n",
" 6.84146951e-01 +7.74422758e-01j,\n",
" 9.45439196e-01 -3.50695577e-01j,\n",
" 1.00057580e+00 -1.27412936e+00j,\n",
" -2.01665427e-01 +1.73591867e+00j,\n",
" 1.76623768e+00 -1.42964638e+00j,\n",
" 1.12282181e+00 +9.02767375e-01j,\n",
" -2.22663798e+00 +7.34672692e-02j,\n",
" 8.26230180e-02 -7.46601441e-01j,\n",
" -4.92810641e-01 -5.49544995e-01j,\n",
" 1.06703148e-01 +1.70461156e+00j,\n",
" -6.84949746e-01 -1.46950185e-01j,\n",
" 1.47509281e+00 +8.77421590e-01j,\n",
" -9.83722735e-03 -6.30597166e-01j,\n",
" 2.36006569e-03 +2.12883413e-01j,\n",
" -3.67899785e-01 -9.95508241e-02j,\n",
" 1.26752410e+00 +2.49370098e-01j,\n",
" -1.19439605e+00 +2.94935402e-01j,\n",
" -2.28617888e+00 +5.70564449e-02j,\n",
" 3.95617530e-02 -2.80995512e-01j,\n",
" -7.51027805e-01 -8.14567233e-01j,\n",
" -1.24160227e+00 -4.13553574e-01j, -4.76626111e-01 +4.16480574e-01j])\n",
"ipdb> u_src\n",
"*** NameError: name 'u_src' is not defined\n",
"ipdb> c\n"
]
}
],
"source": [
"%debug"
]
},
{
"cell_type": "code",
"execution_count": 15,
"execution_count": 18,
"metadata": {
"collapsed": false
},
@@ -509,7 +755,7 @@
"output_type": "stream",
"text": [
"Adjoint test e formulation - getRHSDeriv_m\n",
"(13201.2196403+13827.5790776j) (13201.2196403+13827.5790776j) (-5.45696821064e-12+3.63797880709e-12j) 1.0 True\n"
"(-10077.7224119-28916.4723594j) (-10077.7224119-28916.4723594j) (1.81898940355e-12-3.63797880709e-12j) 1.0 True\n"
]
},
{
@@ -518,7 +764,7 @@
"True"
]
},
"execution_count": 15,
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
@@ -550,19 +796,11 @@
},
{
"cell_type": "code",
"execution_count": 16,
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"2\n"
]
}
],
"outputs": [],
"source": [
"simpeg.mkvc(np.random.randn(survey.nD)+np.random.randn(survey.nD)*1j,2)\n",
"\n",
@@ -571,7 +809,7 @@
},
{
"cell_type": "code",
"execution_count": 17,
"execution_count": 19,
"metadata": {
"collapsed": false
},
@@ -602,7 +840,7 @@
},
{
"cell_type": "code",
"execution_count": 18,
"execution_count": 20,
"metadata": {
"collapsed": false
},
@@ -612,7 +850,7 @@
"output_type": "stream",
"text": [
"Adjoint e formulation - Jvec\n",
"1.96695386678e-05 1.96695386678e-05 3.38813178902e-21 1e-08 True\n"
"1.09508355274e-05 1.09508355274e-05 -1.01643953671e-20 1e-08 True\n"
]
},
{
@@ -621,7 +859,7 @@
"True"
]
},
"execution_count": 18,
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
@@ -629,6 +867,15 @@
"source": [
"JvecAdjointTest()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
@@ -636,18 +883,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
+2 -2
View File
@@ -92,8 +92,8 @@ class BaseMTProblem(BaseFDEMProblem):
# Calculate the projection derivatives
for rx in src.rxList:
# Get the projection derivative
# v should be of size nE,2 (each column for 2 polarizations)
PDeriv_u = lambda v: rx.projectFieldsDeriv(src, self.mesh, u, v) # wrt u, we don't have have PDeriv wrt m
# v should be of size 2*nE (for 2 polarizations)
PDeriv_u = lambda t: rx.projectFieldsDeriv(src, self.mesh, u, t) # wrt u, we don't have have PDeriv wrt m
Jv[src, rx] = PDeriv_u(mkvc(du_dm))
# Return the vectorized sensitivities
return mkvc(Jv)
simpegMT/FieldsMT.py
+2 -2
View File
@@ -233,13 +233,13 @@ class FieldsMT_3D(FieldsMT):
# NOTE: v needs to be length 2*nE to account for both polarizations
def _b_pxSecondaryDeriv_u(self, src, v, adjoint = False):
C = sp.kron(self.mesh.edgeCurl,[[1,0],[0,1]])
C = sp.kron(self.mesh.edgeCurl,[[1,0],[0,0]])
if adjoint:
return - 1./(1j*omega(src.freq)) * (C.T * v)
return - 1./(1j*omega(src.freq)) * (C * v)
def _b_pySecondaryDeriv_u(self, src, v, adjoint = False):
C = sp.kron(self.mesh.edgeCurl,[[1,0],[0,1]])
C = sp.kron(self.mesh.edgeCurl,[[0,0],[0,1]])
if adjoint:
return - 1./(1j*omega(src.freq)) * (C.T * v)
return - 1./(1j*omega(src.freq)) * (C * v)
simpegMT/ProblemMT1D/Problems.py
+3 -2
View File
@@ -65,8 +65,9 @@ class eForm_psField(BaseMTProblem):
dsig_dm = self.curModel.sigmaDeriv
MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
# Need to make the dMf_dsig symmetirc (nN,nN), don't know how to do this
dMf_dsig = self.mesh.getFaceInnerProductDeriv(self.curModel.sigma)(u) * self.curModel.sigmaDeriv
#
u_src = u['e_1dSolution']
dMf_dsig = self.mesh.getFaceInnerProductDeriv(self.curModel.sigma)(u_src) * self.curModel.sigmaDeriv
if adjoint:
return 1j * omega(freq) * ( dMf_dsig.T * v )
# Note: output has to be nN/nF, not nC/nE.
simpegMT/SurveyMT.py
+27 -25
View File
@@ -188,46 +188,48 @@ class RxMT(Survey.BaseRx):
Pby = mesh.getInterpolationMat(bFLocs,'Fy')
# Get the fields at location
# px: x-polaration and py: y-polaration.
ex_px = Pex*mkvc(f[src,'e_px'],2)
ey_px = Pey*mkvc(f[src,'e_px'],2)
ex_py = Pex*mkvc(f[src,'e_py'],2)
ey_py = Pey*mkvc(f[src,'e_py'],2)
hx_px = Pbx*mkvc(f[src,'b_px']/mu_0,2)
hy_px = Pby*mkvc(f[src,'b_px']/mu_0,2)
hx_py = Pbx*mkvc(f[src,'b_py']/mu_0,2)
hy_py = Pby*mkvc(f[src,'b_py']/mu_0,2)
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))
# Derivatives as lambda functions
ex_px_u = lambda vec: sp.hstack((Pex,Pex))*f._e_pxDeriv_u(src,vec)
ey_px_u = lambda vec: sp.hstack((Pey,Pey))*f._e_pxDeriv_u(src,vec)
ex_py_u = lambda vec: sp.hstack((Pex,Pex))*f._e_pyDeriv_u(src,vec)
ey_py_u = lambda vec: sp.hstack((Pey,Pey))*f._e_pyDeriv_u(src,vec)
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))
# NOTE: Think b_p?Deriv_u should return a 2*nF size matrix
hx_px_u = lambda vec: sp.hstack((Pbx,Pbx))*f._b_pxDeriv_u(src,vec)/mu_0
hy_px_u = lambda vec: sp.hstack((Pby,Pby))*f._b_pxDeriv_u(src,vec)/mu_0
hx_py_u = lambda vec: sp.hstack((Pbx,Pbx))*f._b_pyDeriv_u(src,vec)/mu_0
hy_py_u = lambda vec: sp.hstack((Pby,Pby))*f._b_pyDeriv_u(src,vec)/mu_0
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))
# Update the input vector
v = mkvc(v,2) # Make v into a column vector
# v = mkvc(v,2) # Make v into a column vector
# Define the components of the derivative
Hd = Utils.sdiag(1/(hx_px*hy_py - hx_py*hy_px))
Hd_uV = hx_px_u(hy_py*v) - hx_py*hy_px_u(v) + hx_px*hy_py_u(v) - hx_py_u(hy_px*v)
Hd = Utils.sdiag(mkvc(1./(hx_px*hy_py - hx_py*hy_px).data,2))
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:
Zij = ( ex_px*hy_py - ex_py*hy_px)*Hd
ZijN_uV = ex_px_u(hy_py*v) - ex_py*hy_px_u(v) + ex_px*hy_py_u(v) - ex_py_u(hy_px*v)
ZijN_uV = ex_px_u(v)*hy_py + ex_px*hy_py_u(v) - ex_py*hy_px_u(v) - ex_py_u(v)*hy_px
elif 'zxy' in self.rxType:
Zij = (-ex_px*hx_py + ex_py*hx_px)*Hd
ZijN_uV = -ex_px_u(hx_py*v) + ex_py*hx_px_u(v) - ex_px*hx_py_u(v) + ex_py_u(hx_px*v)
ZijN_uV = -ex_px_u(v)*hx_py - ex_px*hx_py_u(v) + ex_py*hx_px_u(v) + ex_py_u(v)*hx_px
elif 'zyx' in self.rxType:
Zij = ( ey_px*hy_py - ey_py*hy_px)*Hd
ZijN_uV = ey_px_u(hy_py*v) - ey_py*hy_px_u(v) + ey_px*hy_py_u(v) - ey_py_u(hy_px*v)
ZijN_uV = ey_px_u(v)*hy_py + ey_px*hy_py_u(v) - ey_py*hy_px_u(v) - ey_py_u(v)*hy_px
elif 'zyy' in self.rxType:
Zij = (-ey_px*hx_py + ey_py*hx_px)*Hd
ZijN_uV = -ey_px_u(hx_py*v) + ey_py*hx_px_u(v) - ey_px*hx_py_u(v) +ey_py_u(hx_px*v)
ZijN_uV = -ey_px_u(v)*hx_py - ey_px*hx_py_u(v) + ey_py*hx_px_u(v) + ey_py_u(v)*hx_px
# Calculate the complex derivative
PDeriv_complex = ZijN_uV*Hd.toarray() - Zij * (Hd_uV*Hd.toarray())
PDeriv_complex = Hd * (ZijN_uV - Zij * Hd_uV )
# Extract the real number for the real/imag components.
Pv = np.array(getattr(PDeriv_complex, real_or_imag))
simpegMT/Tests/test_Problem3D_againstAnalytic.py
+12 -8
View File
@@ -173,18 +173,22 @@ def test_DerivProjfields(inputSetup,comp='All',freq=False):
survey, problem = setupSimpegMTfwd_eForm_ps(inputSetup,comp,freq)
print 'Derivative test of data projection for eFormulation primary/secondary\n\n'
# problem.mapping = simpeg.Maps.ExpMap(problem.mesh)
# Initate things for the derivs Test
src = survey.srcList[0]
rx = src.rxList[0]
# Define a src and rx
src = survey.srcList[-1]
rx = src.rxList[1]
u0 = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j
u0x = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j
u0y = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j
u0 = np.vstack((simpeg.mkvc(u0x,2),simpeg.mkvc(u0y,2)))
f0 = problem.fieldsPair(survey.mesh,survey)
f0[src,'e_pxSolution'] = u0
f0[src,'e_pySolution'] = u0
# u0 = np.hstack((simpeg.mkvc(u0_px,2),simpeg.mkvc(u0_py,2)))
f0[src,'e_pxSolution'] = u0[:len(u0)/2]#u0x
f0[src,'e_pySolution'] = u0[len(u0)/2::]#u0y
def fun(u):
f = problem.fieldsPair(survey.mesh,survey)
f[src,'e_pxSolution'] = u.ravel()
f[src,'e_pySolution'] = u.ravel()
f[src,'e_pxSolution'] = u[:len(u)/2]
f[src,'e_pySolution'] = u[len(u)/2::]
return rx.projectFields(src,survey.mesh,f), lambda t: rx.projectFieldsDeriv(src,survey.mesh,f0,simpeg.mkvc(t,2))
return simpeg.Tests.checkDerivative(fun, u0, num=3, plotIt=False, eps=FLR)