mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-18 12:40:30 +08:00
Inversion problem working.
Fixed 1D problem to correct the phase quadrants.
This commit is contained in:
GudniRos
@@ -51,15 +51,7 @@
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Project at freq: 1.000e+02\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Setup the problem\n",
|
||||
"sigmaHalf = 1e-2\n",
|
||||
@@ -237,13 +229,13 @@
|
||||
"==================== checkDerivative ====================\n",
|
||||
"iter h |ft-f0| |ft-f0-h*J0*dx| Order\n",
|
||||
"---------------------------------------------------------\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",
|
||||
" 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",
|
||||
"========================= PASS! =========================\n",
|
||||
"You are awesome.\n",
|
||||
"That was easy!\n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
@@ -277,7 +269,7 @@
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"array([ 0.00127341])"
|
||||
"array([ 0.00052762])"
|
||||
]
|
||||
},
|
||||
"execution_count": 8,
|
||||
@@ -299,7 +291,7 @@
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"array([[ 0.00173178]])"
|
||||
"array([[ 0.00124017]])"
|
||||
]
|
||||
},
|
||||
"execution_count": 9,
|
||||
@@ -336,17 +328,12 @@
|
||||
"==================== checkDerivative ====================\n",
|
||||
"iter h |ft-f0| |ft-f0-h*J0*dx| Order\n",
|
||||
"---------------------------------------------------------\n",
|
||||
"Project at freq: 1.000e+02\n",
|
||||
"Project at freq: 1.000e+02\n",
|
||||
" 0 1.00e-01 2.646e-06 2.277e-08 nan\n",
|
||||
"Project at freq: 1.000e+02\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 2.627e-08 2.243e-12 2.001\n",
|
||||
"Project at freq: 1.000e+02\n",
|
||||
" 3 1.00e-04 2.627e-09 2.242e-14 2.000\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",
|
||||
"========================= PASS! =========================\n",
|
||||
"Yay passed!\n",
|
||||
"You get a gold star!\n",
|
||||
"\n"
|
||||
]
|
||||
},
|
||||
@@ -411,7 +398,7 @@
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Adjoint e formulation - projectFieldsDeriv\n",
|
||||
"-0.000991536643367 -0.000991536643367 2.16840434497e-19 1e-07 True\n"
|
||||
"-2.26989762698e-05 -2.26989762698e-05 0.0 1e-08 True\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -456,32 +443,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": 15,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Adjoint test e formulation - getADeriv_m\n",
|
||||
"(340193.379835-398622.996348j) (340193.379835-398622.996348j) (-2.32830643654e-10-3.49245965481e-10j) 10.0 True\n"
|
||||
"(-1977540.36505+2093781.70221j) (-1977540.36505+2093781.70221j) (-1.86264514923e-09+2.79396772385e-09j) 100.0 True\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -490,7 +458,7 @@
|
||||
"True"
|
||||
]
|
||||
},
|
||||
"execution_count": 15,
|
||||
"execution_count": 14,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
@@ -531,7 +499,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 16,
|
||||
"execution_count": 15,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -541,7 +509,7 @@
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Adjoint test e formulation - getRHSDeriv_m\n",
|
||||
"(-12351.1349263+433.4655562j) (-12351.1349263+433.4655562j) (-9.09494701773e-12-1.22781784739e-11j) 1.0 True\n"
|
||||
"(13201.2196403+13827.5790776j) (13201.2196403+13827.5790776j) (-5.45696821064e-12+3.63797880709e-12j) 1.0 True\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -550,7 +518,7 @@
|
||||
"True"
|
||||
]
|
||||
},
|
||||
"execution_count": 16,
|
||||
"execution_count": 15,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
@@ -582,7 +550,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 17,
|
||||
"execution_count": 16,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -603,7 +571,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 18,
|
||||
"execution_count": 17,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -634,7 +602,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 19,
|
||||
"execution_count": 18,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -644,7 +612,7 @@
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Adjoint e formulation - Jvec\n",
|
||||
"-3.61480355369e-05 -3.61480355369e-05 -2.71050543121e-20 1e-08 True\n"
|
||||
"1.96695386678e-05 1.96695386678e-05 3.38813178902e-21 1e-08 True\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -653,7 +621,7 @@
|
||||
"True"
|
||||
]
|
||||
},
|
||||
"execution_count": 19,
|
||||
"execution_count": 18,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
@@ -661,99 +629,6 @@
|
||||
"source": [
|
||||
"JvecAdjointTest()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 20,
|
||||
"metadata": {
|
||||
"collapsed": false,
|
||||
"scrolled": true
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stderr",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"ERROR: No traceback has been produced, nothing to debug.\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"%debug"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 21,
|
||||
"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": 21,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"survey.dpred(sigma)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 22,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"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::])"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
@@ -761,6 +636,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,
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
@@ -25,12 +25,12 @@
|
||||
"\n",
|
||||
"# Frequency\n",
|
||||
"nFreq = 33\n",
|
||||
"# freqs = np.logspace(3,-3,nFreq)\n",
|
||||
"freqs = np.array([100,10,1,0.1,0.01])\n",
|
||||
"freqs = np.logspace(3,-3,nFreq)\n",
|
||||
"# freqs = np.array([100,10,1,0.1,0.01])\n",
|
||||
"# Make the mesh\n",
|
||||
"ct = 5\n",
|
||||
"air = simpeg.Utils.meshTensor([(ct,25,1.3)])\n",
|
||||
"core = np.concatenate( ( np.kron(simpeg.Utils.meshTensor([(ct,15,-1.2)]),np.ones((10,))) , simpeg.Utils.meshTensor([(ct,20)]) ) )\n",
|
||||
"ct = 10\n",
|
||||
"air = simpeg.Utils.meshTensor([(ct,15,1.3)])\n",
|
||||
"core = np.concatenate( ( np.kron(simpeg.Utils.meshTensor([(ct,15,-1.2)]),np.ones((5,))) , simpeg.Utils.meshTensor([(ct,5)]) ) )\n",
|
||||
"bot = simpeg.Utils.meshTensor([(core[0],10,-1.3)])\n",
|
||||
"x0 = -np.array([np.sum(np.concatenate((core,bot)))])\n",
|
||||
"# Change to use no air\n",
|
||||
@@ -38,20 +38,24 @@
|
||||
"\n",
|
||||
"## Setup model varibles\n",
|
||||
"active = m1d.vectorCCx<0.\n",
|
||||
"layer = (m1d.vectorCCx<-200.) & (m1d.vectorCCx>=-600.)\n",
|
||||
"layer1 = (m1d.vectorCCx<-200.) & (m1d.vectorCCx>=-600.)\n",
|
||||
"layer2 = (m1d.vectorCCx<-2000.) & (m1d.vectorCCx>=-4000.)\n",
|
||||
"actMap = simpeg.Maps.ActiveCells(m1d, active, np.log(1e-8), nC=m1d.nCx)\n",
|
||||
"mappingExpAct = simpeg.Maps.ExpMap(m1d) * actMap\n",
|
||||
"sig_half = 2e-3\n",
|
||||
"sig_air = 1e-8\n",
|
||||
"sig_layer = 1e-3\n",
|
||||
"sig_layer1 = 1\n",
|
||||
"sig_layer2 = .1\n",
|
||||
"# Make the true model\n",
|
||||
"sigma_true = np.ones(m1d.nCx)*sig_air\n",
|
||||
"sigma_true[active] = sig_half\n",
|
||||
"sigma_true[layer] = sig_layer\n",
|
||||
"sigma_true[layer1] = sig_layer1\n",
|
||||
"sigma_true[layer2] = sig_layer2\n",
|
||||
"m_true = np.log(sigma_true[active])\n",
|
||||
"# Make the background model\n",
|
||||
"sigma_0 = np.ones(m_true.size)*sig_half\n",
|
||||
"m_0 = np.log10(np.ones(m_true.size)*sig_half)\n",
|
||||
"sigma_0 = np.ones(m1d.nCx)*sig_air\n",
|
||||
"sigma_0[active] = sig_half\n",
|
||||
"m_0 = np.log(sigma_0[active])\n",
|
||||
"\n",
|
||||
"# Receivers \n",
|
||||
"# 1D impedance at the surface (elevation 0)\n",
|
||||
@@ -69,9 +73,11 @@
|
||||
" srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freq,sigma_0))\n",
|
||||
"# Make the survey\n",
|
||||
"survey = simpegmt.SurveyMT.SurveyMT(srcList)\n",
|
||||
"survey.mtrue = sigma_true\n",
|
||||
"survey.mtrue = m_true\n",
|
||||
"# Set the problem\n",
|
||||
"problem = simpegmt.ProblemMT1D.eForm_psField(m1d,mapping=mappingExpAct)\n",
|
||||
"from pymatsolver import MumpsSolver\n",
|
||||
"problem.solver = MumpsSolver\n",
|
||||
"problem.pair(survey)"
|
||||
]
|
||||
},
|
||||
@@ -81,20 +87,9 @@
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"180"
|
||||
]
|
||||
},
|
||||
"execution_count": 3,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"np.sum(active)"
|
||||
"# problem.mapping.sigmaMap._transform(m_0)"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -105,29 +100,27 @@
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"ename": "Exception",
|
||||
"evalue": "Unexpected shape of tensor",
|
||||
"output_type": "error",
|
||||
"traceback": [
|
||||
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
|
||||
"\u001b[1;31mException\u001b[0m Traceback (most recent call last)",
|
||||
"\u001b[1;32m<ipython-input-4-785d1f579efc>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m 1\u001b[0m \u001b[1;31m## Make the observed data\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 2\u001b[0m \u001b[1;31m# Project the data\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 3\u001b[1;33m \u001b[0md_true\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0msurvey\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdpred\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mm_true\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 4\u001b[0m \u001b[0msurvey\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdtrue\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0md_true\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 5\u001b[0m \u001b[1;31m# Add noise\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/CounterUtils.pyc\u001b[0m in \u001b[0;36mwrapper\u001b[1;34m(self, *args, **kwargs)\u001b[0m\n\u001b[0;32m 81\u001b[0m \u001b[0mcounter\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mgetattr\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'counter'\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mNone\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 82\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0mtype\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mcounter\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mis\u001b[0m \u001b[0mCounter\u001b[0m\u001b[1;33m:\u001b[0m \u001b[0mcounter\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcount\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m__class__\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m__name__\u001b[0m\u001b[1;33m+\u001b[0m\u001b[1;34m'.'\u001b[0m\u001b[1;33m+\u001b[0m\u001b[0mf\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m__name__\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 83\u001b[1;33m \u001b[0mout\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mf\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 84\u001b[0m \u001b[1;32mreturn\u001b[0m \u001b[0mout\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 85\u001b[0m \u001b[1;32mreturn\u001b[0m \u001b[0mwrapper\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpeg/SimPEG/Utils/codeutils.pyc\u001b[0m in \u001b[0;36mrequiresVarWrapper\u001b[1;34m(self, *args, **kwargs)\u001b[0m\n\u001b[0;32m 224\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0mgetattr\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mvar\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[0;32m 225\u001b[0m \u001b[1;32mraise\u001b[0m \u001b[0mException\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mextra\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 226\u001b[1;33m \u001b[1;32mreturn\u001b[0m \u001b[0mf\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 227\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 228\u001b[0m \u001b[0mdoc\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mrequiresVarWrapper\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m__doc__\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;36mdpred\u001b[1;34m(self, m, u)\u001b[0m\n\u001b[0;32m 306\u001b[0m \u001b[0mWhere\u001b[0m \u001b[0mP\u001b[0m \u001b[1;32mis\u001b[0m \u001b[0ma\u001b[0m \u001b[0mprojection\u001b[0m \u001b[0mof\u001b[0m \u001b[0mthe\u001b[0m \u001b[0mfields\u001b[0m \u001b[0monto\u001b[0m \u001b[0mthe\u001b[0m \u001b[0mdata\u001b[0m \u001b[0mspace\u001b[0m\u001b[1;33m.\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 307\u001b[0m \"\"\"\n\u001b[1;32m--> 308\u001b[1;33m \u001b[1;32mif\u001b[0m \u001b[0mu\u001b[0m \u001b[1;32mis\u001b[0m \u001b[0mNone\u001b[0m\u001b[1;33m:\u001b[0m \u001b[0mu\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mprob\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfields\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mm\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 309\u001b[0m \u001b[1;32mreturn\u001b[0m \u001b[0mUtils\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmkvc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mprojectFields\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[0;32m 310\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpegmt/simpegMT/ProblemMT1D/Problems.pyc\u001b[0m in \u001b[0;36mfields\u001b[1;34m(self, m)\u001b[0m\n\u001b[0;32m 102\u001b[0m \u001b[0msys\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mstdout\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mflush\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 103\u001b[0m \u001b[0mA\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetA\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfreq\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 104\u001b[1;33m \u001b[0mrhs\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetRHS\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfreq\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 105\u001b[0m \u001b[0mAinv\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mSolver\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mA\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msolverOpts\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 106\u001b[0m \u001b[0me_s\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mAinv\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mrhs\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpegmt/simpegMT/ProblemMT1D/Problems.pyc\u001b[0m in \u001b[0;36mgetRHS\u001b[1;34m(self, freq)\u001b[0m\n\u001b[0;32m 74\u001b[0m \u001b[1;31m# Get sources for the frequncy(polarizations)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 75\u001b[0m \u001b[0mSrc\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[0mgetSrcByFreq\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfreq\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 76\u001b[1;33m \u001b[0mS_e\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mSrc\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mS_e\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 77\u001b[0m \u001b[1;32mreturn\u001b[0m \u001b[1;33m-\u001b[0m\u001b[1;36m1j\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0momega\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfreq\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mS_e\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 78\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpegmt/simpegMT/SurveyMT.pyc\u001b[0m in \u001b[0;36mS_e\u001b[1;34m(self, problem)\u001b[0m\n\u001b[0;32m 244\u001b[0m \u001b[0mGet\u001b[0m \u001b[0mthe\u001b[0m \u001b[0melectrical\u001b[0m \u001b[0mfield\u001b[0m \u001b[0msource\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 245\u001b[0m \"\"\"\n\u001b[1;32m--> 246\u001b[1;33m \u001b[0me_p\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mePrimary\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mproblem\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 247\u001b[0m \u001b[0mMap_sigma_p\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mMaps\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mVertical1DMap\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mproblem\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmesh\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 248\u001b[0m \u001b[0msigma_p\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mMap_sigma_p\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_transform\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msigma1d\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/SurveyMT.pyc\u001b[0m in \u001b[0;36mePrimary\u001b[1;34m(self, problem)\u001b[0m\n\u001b[0;32m 227\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0mePrimary\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mproblem\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 228\u001b[0m \u001b[1;31m# Get primary fields for both polarizations\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 229\u001b[1;33m \u001b[0meBG_bp\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mhomo1DModelSource\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mproblem\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmesh\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfreq\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msigma1d\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 230\u001b[0m \u001b[1;32mreturn\u001b[0m \u001b[0meBG_bp\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 231\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpegmt/simpegMT/Sources/backgroundModelSources.pyc\u001b[0m in \u001b[0;36mhomo1DModelSource\u001b[1;34m(mesh, freq, sigma_1d)\u001b[0m\n\u001b[0;32m 22\u001b[0m \u001b[0mmesh1d\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0msimpeg\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mMesh\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mTensorMesh\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mmesh\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mhz\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0marray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mmesh\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mx0\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m2\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[0m\n\u001b[0;32m 23\u001b[0m \u001b[1;31m# # Note: Everything is using e^iwt\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 24\u001b[1;33m \u001b[0me0_1d\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mget1DEfields\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmesh1d\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0msigma_1d\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mfreq\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 25\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0mmesh\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdim\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[0;32m 26\u001b[0m \u001b[0meBG_px\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m-\u001b[0m\u001b[0msimpeg\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmkvc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0me0_1d\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m2\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/Utils/MT1Dsolutions.pyc\u001b[0m in \u001b[0;36mget1DEfields\u001b[1;34m(m1d, sigma, freq, sourceAmp)\u001b[0m\n\u001b[0;32m 12\u001b[0m \u001b[0mMmu\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0msimpeg\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mUtils\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msdiag\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mm1d\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mvol\u001b[0m\u001b[1;33m*\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[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 13\u001b[0m \u001b[1;31m# Conductivity\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 14\u001b[1;33m \u001b[0mMsig\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mm1d\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetFaceInnerProduct\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msigma\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 15\u001b[0m \u001b[1;31m# Set up the solution matrix\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 16\u001b[0m \u001b[0mA\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mG\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mT\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mMmu\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mG\u001b[0m \u001b[1;33m+\u001b[0m \u001b[1;36m1j\u001b[0m\u001b[1;33m*\u001b[0m\u001b[1;36m2.\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mpi\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mfreq\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mMsig\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpeg/SimPEG/Mesh/InnerProducts.pyc\u001b[0m in \u001b[0;36mgetFaceInnerProduct\u001b[1;34m(self, prop, invProp, invMat, doFast)\u001b[0m\n\u001b[0;32m 20\u001b[0m \u001b[1;33m:\u001b[0m\u001b[1;32mreturn\u001b[0m\u001b[1;33m:\u001b[0m \u001b[0mM\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mthe\u001b[0m \u001b[0minner\u001b[0m \u001b[0mproduct\u001b[0m \u001b[0mmatrix\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mnF\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnF\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 21\u001b[0m \"\"\"\n\u001b[1;32m---> 22\u001b[1;33m \u001b[1;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_getInnerProduct\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'F'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mprop\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mprop\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0minvProp\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0minvProp\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0minvMat\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0minvMat\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdoFast\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mdoFast\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 23\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 24\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0mgetEdgeInnerProduct\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mprop\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0minvProp\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mFalse\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0minvMat\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mFalse\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdoFast\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[1;32m/media/gudni/ExtraDrive1/Codes/python/simpeg/SimPEG/Mesh/InnerProducts.pyc\u001b[0m in \u001b[0;36m_getInnerProduct\u001b[1;34m(self, projType, prop, invProp, invMat, doFast)\u001b[0m\n\u001b[0;32m 54\u001b[0m \u001b[0mprop\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0minvPropertyTensor\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mprop\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 55\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 56\u001b[1;33m \u001b[0mtensorType\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mTensorType\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mprop\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 57\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 58\u001b[0m \u001b[0mMu\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mmakePropertyTensor\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mprop\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/matutils.pyc\u001b[0m in \u001b[0;36m__init__\u001b[1;34m(self, M, tensor)\u001b[0m\n\u001b[0;32m 272\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_tts\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;34m'tensor'\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 273\u001b[0m \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 274\u001b[1;33m \u001b[1;32mraise\u001b[0m \u001b[0mException\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'Unexpected shape of tensor'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 275\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m__str__\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 276\u001b[0m \u001b[1;32mreturn\u001b[0m \u001b[1;34m'TensorType[%i]: %s'\u001b[0m \u001b[1;33m%\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_tt\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_tts\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
|
||||
"\u001b[1;31mException\u001b[0m: Unexpected shape of tensor"
|
||||
]
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"0.002"
|
||||
]
|
||||
},
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"sig_half"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 5,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"## Make the observed data \n",
|
||||
"# Project the data\n",
|
||||
@@ -150,19 +143,30 @@
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 6,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"## Setup the inversion proceedure\n",
|
||||
"C = simpeg.Utils.Counter()\n",
|
||||
"\n",
|
||||
"# Set the optimization\n",
|
||||
"opt = simpeg.Optimization.InexactGaussNewton(maxIter = 10)\n",
|
||||
"opt = simpeg.Optimization.InexactGaussNewton(maxIter = 30)\n",
|
||||
"opt.counter = C\n",
|
||||
"opt.LSshorten = 0.5\n",
|
||||
"opt.remember('xc')\n",
|
||||
"# Data misfit\n",
|
||||
"dmis = simpeg.DataMisfit.l2_DataMisfit(survey)\n",
|
||||
"# Regularization\n",
|
||||
"reg = simpeg.Regularization.Tikhonov(m1d)\n",
|
||||
"# regMesh = simpeg.Mesh.TensorMesh([m1d.hx[problem.mapping.sigmaMap.maps[-1].indActive]])\n",
|
||||
"# reg = simpeg.Regularization.Tikhonov(regMesh)\n",
|
||||
"reg = simpeg.Regularization.Tikhonov(m1d,mapping=mappingExpAct)\n",
|
||||
"reg.alpha_s = 1e-5\n",
|
||||
"reg.alpha_x = 1.\n",
|
||||
"\n",
|
||||
@@ -174,7 +178,355 @@
|
||||
"betaest = simpeg.Directives.BetaEstimate_ByEig(beta0_ratio=0.75)\n",
|
||||
"saveModel = simpeg.Directives.SaveModelEveryIteration()\n",
|
||||
"# Create an inversion object\n",
|
||||
"inv = simpeg.Inversion.BaseInversion(invProb, directiveList=[beta,betaest,saveModel]) \n"
|
||||
"inv = simpeg.Inversion.BaseInversion(invProb, directiveList=[beta,betaest])#,saveModel]) \n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 7,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"<SimPEG.Maps.ActiveCells at 0x7fcfc8df5650>"
|
||||
]
|
||||
},
|
||||
"execution_count": 7,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"problem.mapping.sigmaMap.maps[-1]"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 8,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"SimPEG.InvProblem will set Regularization.mref to m0.\n",
|
||||
"SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.\n",
|
||||
" ***Done using same solver as the problem***\n",
|
||||
"SimPEG.l2_DataMisfit is creating default weightings for Wd.\n",
|
||||
"============================ Inexact Gauss Newton ============================\n",
|
||||
" # beta phi_d phi_m f |proj(x-g)-x| LS Comment \n",
|
||||
"-----------------------------------------------------------------------------\n",
|
||||
" 0 4.83e+05 1.46e+06 6.90e-07 1.46e+06 2.09e+05 0 \n",
|
||||
" 1 4.83e+05 1.80e+05 5.76e-06 1.80e+05 2.64e+04 0 \n",
|
||||
" 2 4.83e+05 1.29e+05 7.58e-06 1.29e+05 2.05e+04 0 Skip BFGS \n",
|
||||
" 3 6.04e+04 1.10e+05 8.89e-06 1.10e+05 1.82e+04 0 Skip BFGS \n",
|
||||
" 4 6.04e+04 5.48e+04 1.97e-05 5.48e+04 1.05e+04 0 Skip BFGS \n",
|
||||
" 5 6.04e+04 4.57e+04 2.46e-05 4.57e+04 9.11e+03 0 Skip BFGS \n",
|
||||
" 6 7.55e+03 4.04e+04 2.88e-05 4.04e+04 8.25e+03 0 Skip BFGS \n",
|
||||
" 7 7.55e+03 2.22e+04 6.11e-05 2.22e+04 5.14e+03 0 Skip BFGS \n",
|
||||
" 8 7.55e+03 1.81e+04 7.86e-05 1.81e+04 4.38e+03 0 Skip BFGS \n",
|
||||
" 9 9.44e+02 1.58e+04 9.29e-05 1.58e+04 3.94e+03 0 Skip BFGS \n",
|
||||
" 10 9.44e+02 8.37e+03 1.97e-04 8.37e+03 2.43e+03 0 Skip BFGS \n",
|
||||
" 11 9.44e+02 6.71e+03 2.54e-04 6.71e+03 2.05e+03 0 Skip BFGS \n",
|
||||
" 12 1.18e+02 5.78e+03 3.01e-04 5.78e+03 1.83e+03 0 Skip BFGS \n",
|
||||
" 13 1.18e+02 3.04e+03 6.06e-04 3.04e+03 1.11e+03 0 Skip BFGS \n",
|
||||
" 14 1.18e+02 2.42e+03 7.75e-04 2.42e+03 9.25e+02 0 Skip BFGS \n",
|
||||
" 15 1.47e+01 2.09e+03 9.10e-04 2.09e+03 8.15e+02 0 Skip BFGS \n",
|
||||
" 16 1.47e+01 1.21e+03 1.66e-03 1.21e+03 4.81e+02 0 Skip BFGS \n",
|
||||
" 17 1.47e+01 9.97e+02 2.07e-03 9.97e+02 3.82e+02 0 Skip BFGS \n",
|
||||
" 18 1.84e+00 8.79e+02 2.36e-03 8.79e+02 3.30e+02 0 Skip BFGS \n",
|
||||
" 19 1.84e+00 6.16e+02 3.62e-03 6.16e+02 1.99e+02 0 Skip BFGS \n",
|
||||
" 20 1.84e+00 5.31e+02 4.24e-03 5.31e+02 1.64e+02 0 Skip BFGS \n",
|
||||
" 21 2.30e-01 4.79e+02 4.68e-03 4.79e+02 1.45e+02 0 Skip BFGS \n",
|
||||
" 22 2.30e-01 3.38e+02 6.36e-03 3.38e+02 1.01e+02 0 Skip BFGS \n",
|
||||
" 23 2.30e-01 2.80e+02 7.02e-03 2.80e+02 8.93e+01 0 Skip BFGS \n",
|
||||
" 24 2.88e-02 2.32e+02 7.55e-03 2.32e+02 8.14e+01 0 Skip BFGS \n",
|
||||
" 25 2.88e-02 1.53e+02 1.08e-02 1.53e+02 6.94e+01 0 Skip BFGS \n",
|
||||
" 26 2.88e-02 8.67e+01 1.27e-02 8.67e+01 3.84e+01 0 \n",
|
||||
" 27 3.60e-03 6.54e+01 1.40e-02 6.54e+01 3.13e+01 0 Skip BFGS \n",
|
||||
" 28 3.60e-03 3.48e+01 1.69e-02 3.48e+01 2.06e+01 0 Skip BFGS \n",
|
||||
" 29 3.60e-03 2.63e+01 1.86e-02 2.63e+01 1.07e+01 0 Skip BFGS \n",
|
||||
" 30 4.50e-04 2.19e+01 2.17e-02 2.19e+01 7.17e+00 0 Skip BFGS \n",
|
||||
"------------------------- STOP! -------------------------\n",
|
||||
"1 : |fc-fOld| = 4.3709e+00 <= tolF*(1+|f0|) = 1.4560e+05\n",
|
||||
"1 : |xc-x_last| = 2.1874e+00 <= tolX*(1+|x0|) = 5.9957e+00\n",
|
||||
"0 : |proj(x-g)-x| = 7.1673e+00 <= tolG = 1.0000e-01\n",
|
||||
"0 : |proj(x-g)-x| = 7.1673e+00 <= 1e3*eps = 1.0000e-02\n",
|
||||
"1 : maxIter = 30 <= iter = 30\n",
|
||||
"------------------------- DONE! -------------------------\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# Runn the inversion\n",
|
||||
"mopt = inv.run(m_0)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 18,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"[{'left': <function SimPEG.Optimization.<lambda>>,\n",
|
||||
" 'right': <function SimPEG.Optimization.<lambda>>,\n",
|
||||
" 'stopType': 'critical',\n",
|
||||
" 'str': '%d : maxIter = %3d <= iter = %3d'}]"
|
||||
]
|
||||
},
|
||||
"execution_count": 18,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": []
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 9,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"## Setup the inversion proceedure\n",
|
||||
"C = simpeg.Utils.Counter()\n",
|
||||
"\n",
|
||||
"# Set the optimization\n",
|
||||
"optc = simpeg.Optimization.InexactGaussNewton(maxIter = 20)\n",
|
||||
"optc.counter = C\n",
|
||||
"optc.LSshorten = 0.5\n",
|
||||
"optc.remember('xc')\n",
|
||||
"# Data misfit\n",
|
||||
"dmisc = simpeg.DataMisfit.l2_DataMisfit(survey)\n",
|
||||
"# Regularization\n",
|
||||
"# regMesh = simpeg.Mesh.TensorMesh([m1d.hx[problem.mapping.sigmaMap.maps[-1].indActive]])\n",
|
||||
"# reg = simpeg.Regularization.Tikhonov(regMesh)\n",
|
||||
"regc = simpeg.Regularization.Tikhonov(m1d,mapping=mappingExpAct)\n",
|
||||
"regc.alpha_s = 1e-5\n",
|
||||
"regc.alpha_x = 1.\n",
|
||||
"# Inversion problem\n",
|
||||
"invProbc = simpeg.InvProblem.BaseInvProblem(dmisc, regc, optc)\n",
|
||||
"invProbc.counter = C\n",
|
||||
"# Beta cooling\n",
|
||||
"betac = simpeg.Directives.BetaSchedule()\n",
|
||||
"betaestc = simpeg.Directives.BetaEstimate_ByEig(beta0_ratio=0.75)\n",
|
||||
"saveModel = simpeg.Directives.SaveModelEveryIteration()\n",
|
||||
"# Create an inversion object\n",
|
||||
"invc = simpeg.Inversion.BaseInversion(invProbc, directiveList=[betac,betaestc])"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 10,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"SimPEG.InvProblem will set Regularization.mref to m0.\n",
|
||||
"SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.\n",
|
||||
" ***Done using same solver as the problem***\n",
|
||||
"SimPEG.l2_DataMisfit is creating default weightings for Wd.\n",
|
||||
"============================ Inexact Gauss Newton ============================\n",
|
||||
" # beta phi_d phi_m f |proj(x-g)-x| LS Comment \n",
|
||||
"-----------------------------------------------------------------------------\n",
|
||||
" 0 9.87e+02 2.19e+01 2.35e-02 4.51e+01 3.23e+01 0 \n",
|
||||
" 1 9.87e+02 2.38e+01 1.74e-02 4.09e+01 7.61e+00 0 \n",
|
||||
" 2 9.87e+02 2.30e+01 1.79e-02 4.06e+01 2.68e+00 0 \n",
|
||||
" 3 1.23e+02 2.28e+01 1.79e-02 2.50e+01 1.54e+01 0 \n",
|
||||
" 4 1.23e+02 2.07e+01 2.24e-02 2.35e+01 5.36e+00 0 \n",
|
||||
" 5 1.23e+02 2.03e+01 2.22e-02 2.30e+01 6.78e+00 0 Skip BFGS \n",
|
||||
" 6 1.54e+01 2.00e+01 2.33e-02 2.03e+01 7.34e+00 0 Skip BFGS \n",
|
||||
" 7 1.54e+01 1.95e+01 3.01e-02 2.00e+01 7.09e+00 0 \n",
|
||||
" 8 1.54e+01 1.92e+01 3.30e-02 1.97e+01 5.00e+00 0 \n",
|
||||
" 9 1.93e+00 1.91e+01 3.72e-02 1.91e+01 5.68e+00 0 \n",
|
||||
" 10 1.93e+00 1.88e+01 6.62e-02 1.89e+01 8.43e+00 1 \n",
|
||||
" 11 1.93e+00 1.82e+01 1.12e-01 1.84e+01 5.37e+00 0 \n",
|
||||
" 12 2.41e-01 1.81e+01 1.21e-01 1.81e+01 5.14e+00 0 \n",
|
||||
" 13 2.41e-01 1.80e+01 1.18e-01 1.80e+01 3.93e+00 0 \n",
|
||||
" 14 2.41e-01 1.80e+01 1.43e-01 1.80e+01 6.08e+00 0 \n",
|
||||
" 15 3.01e-02 1.78e+01 1.19e-01 1.78e+01 4.97e+00 1 \n",
|
||||
" 16 3.01e-02 1.78e+01 1.13e-01 1.78e+01 2.41e+00 0 Skip BFGS \n",
|
||||
" 17 3.01e-02 1.78e+01 1.37e-01 1.78e+01 4.47e+00 0 \n",
|
||||
" 18 3.77e-03 1.76e+01 1.28e-01 1.76e+01 6.24e+00 0 Skip BFGS \n",
|
||||
" 19 3.77e-03 1.76e+01 1.07e-01 1.76e+01 7.34e+00 0 \n",
|
||||
" 20 3.77e-03 1.75e+01 1.25e-01 1.75e+01 6.14e+00 1 \n",
|
||||
"------------------------- STOP! -------------------------\n",
|
||||
"1 : |fc-fOld| = 4.7935e-02 <= tolF*(1+|f0|) = 4.6080e+00\n",
|
||||
"1 : |xc-x_last| = 3.4384e+00 <= tolX*(1+|x0|) = 4.0560e+00\n",
|
||||
"0 : |proj(x-g)-x| = 6.1413e+00 <= tolG = 1.0000e-01\n",
|
||||
"0 : |proj(x-g)-x| = 6.1413e+00 <= 1e3*eps = 1.0000e-02\n",
|
||||
"1 : maxIter = 20 <= iter = 20\n",
|
||||
"------------------------- DONE! -------------------------\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"mopt2 = invc.run(mopt)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 24,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"moptc=mopt2"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 26,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"## Setup the inversion proceedure\n",
|
||||
"C = simpeg.Utils.Counter()\n",
|
||||
"\n",
|
||||
"# Set the optimization\n",
|
||||
"optc1 = simpeg.Optimization.InexactGaussNewton(maxIter = 20)\n",
|
||||
"optc1.counter = C\n",
|
||||
"optc1.LSshorten = 0.1\n",
|
||||
"optc1.remember('xc')\n",
|
||||
"# Data misfit\n",
|
||||
"dmisc1 = simpeg.DataMisfit.l2_DataMisfit(survey)\n",
|
||||
"# Regularization\n",
|
||||
"# regMesh = simpeg.Mesh.TensorMesh([m1d.hx[problem.mapping.sigmaMap.maps[-1].indActive]])\n",
|
||||
"# reg = simpeg.Regularization.Tikhonov(regMesh)\n",
|
||||
"regc1 = simpeg.Regularization.Tikhonov(m1d,mapping=mappingExpAct)\n",
|
||||
"regc1.alpha_s = 1e-5\n",
|
||||
"regc1.alpha_x = 1.\n",
|
||||
"regc1.mref = reg.mref\n",
|
||||
"# Inversion problem\n",
|
||||
"invProbc1 = simpeg.InvProblem.BaseInvProblem(dmisc1, regc1, optc1)\n",
|
||||
"invProbc1.counter = C\n",
|
||||
"# Beta cooling\n",
|
||||
"betac1 = simpeg.Directives.BetaSchedule()\n",
|
||||
"betaestc1 = simpeg.Directives.BetaEstimate_ByEig(beta0_ratio=0.75)\n",
|
||||
"betaestc1.beta0 = 3.60e-03\n",
|
||||
"saveModel = simpeg.Directives.SaveModelEveryIteration()\n",
|
||||
"# Create an inversion object\n",
|
||||
"invc1 = simpeg.Inversion.BaseInversion(invProbc1, directiveList=[betac1,betaestc1])"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 27,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.\n",
|
||||
" ***Done using same solver as the problem***\n",
|
||||
"SimPEG.l2_DataMisfit is creating default weightings for Wd.\n",
|
||||
"============================ Inexact Gauss Newton ============================\n",
|
||||
" # beta phi_d phi_m f |proj(x-g)-x| LS Comment \n",
|
||||
"-----------------------------------------------------------------------------\n",
|
||||
" 0 1.79e-02 1.75e+01 1.23e-01 1.75e+01 6.14e+00 0 \n",
|
||||
" 1 1.79e-02 1.75e+01 1.23e-01 1.75e+01 6.14e+00 3 \n",
|
||||
" 2 1.79e-02 1.75e+01 1.23e-01 1.75e+01 6.15e+00 3 Skip BFGS \n",
|
||||
" 3 2.23e-03 1.75e+01 1.23e-01 1.75e+01 6.16e+00 3 Skip BFGS \n",
|
||||
" 4 2.23e-03 1.75e+01 1.23e-01 1.75e+01 6.17e+00 3 Skip BFGS \n",
|
||||
" 5 2.23e-03 1.75e+01 1.23e-01 1.75e+01 6.17e+00 3 Skip BFGS \n",
|
||||
" 6 2.79e-04 1.75e+01 1.23e-01 1.75e+01 6.17e+00 3 Skip BFGS \n",
|
||||
" 7 2.79e-04 1.75e+01 1.23e-01 1.75e+01 6.17e+00 3 Skip BFGS \n",
|
||||
" 8 2.79e-04 1.75e+01 1.23e-01 1.75e+01 6.17e+00 3 Skip BFGS \n",
|
||||
" 9 3.49e-05 1.75e+01 1.23e-01 1.75e+01 6.17e+00 3 Skip BFGS \n",
|
||||
" 10 3.49e-05 1.75e+01 1.23e-01 1.75e+01 6.17e+00 3 Skip BFGS \n",
|
||||
" 11 3.49e-05 1.75e+01 1.23e-01 1.75e+01 6.17e+00 3 Skip BFGS \n",
|
||||
" 12 4.36e-06 1.75e+01 1.23e-01 1.75e+01 6.16e+00 3 Skip BFGS \n",
|
||||
" 13 4.36e-06 1.75e+01 1.23e-01 1.75e+01 6.16e+00 3 Skip BFGS \n",
|
||||
" 14 4.36e-06 1.75e+01 1.23e-01 1.75e+01 6.29e+00 2 Skip BFGS \n",
|
||||
" 15 5.46e-07 1.75e+01 1.23e-01 1.75e+01 6.32e+00 2 Skip BFGS \n",
|
||||
" 16 5.46e-07 1.75e+01 1.23e-01 1.75e+01 6.31e+00 2 Skip BFGS \n",
|
||||
" 17 5.46e-07 1.75e+01 1.23e-01 1.75e+01 6.29e+00 2 Skip BFGS \n",
|
||||
" 18 6.82e-08 1.75e+01 1.23e-01 1.75e+01 6.26e+00 2 Skip BFGS \n",
|
||||
" 19 6.82e-08 1.75e+01 1.23e-01 1.75e+01 6.22e+00 2 Skip BFGS \n",
|
||||
" 20 6.82e-08 1.75e+01 1.23e-01 1.75e+01 6.18e+00 2 Skip BFGS \n",
|
||||
"------------------------- STOP! -------------------------\n",
|
||||
"1 : |fc-fOld| = 3.9817e-03 <= tolF*(1+|f0|) = 1.8547e+00\n",
|
||||
"0 : |xc-x_last| = 4.1858e+01 <= tolX*(1+|x0|) = 5.7798e+00\n",
|
||||
"0 : |proj(x-g)-x| = 6.1825e+00 <= tolG = 1.0000e-01\n",
|
||||
"0 : |proj(x-g)-x| = 6.1825e+00 <= 1e3*eps = 1.0000e-02\n",
|
||||
"1 : maxIter = 20 <= iter = 20\n",
|
||||
"------------------------- DONE! -------------------------\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"moptc1 = invc1.run(mopt2)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 11,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Counters:\n",
|
||||
" InexactGaussNewton.doEndIteration : 30\n",
|
||||
" InexactGaussNewton.doStartIteration : 31\n",
|
||||
" InexactGaussNewton.scaleSearchDirection : 30\n",
|
||||
"\n",
|
||||
"Times: mean sum\n",
|
||||
" BaseInvProblem.evalFunction : 3.35e+00, 2.04e+02, 61x\n",
|
||||
" InexactGaussNewton.findSearchDirection : 2.11e+01, 6.34e+02, 30x\n",
|
||||
" InexactGaussNewton.minimize : 8.39e+02, 8.39e+02, 1x\n",
|
||||
" InexactGaussNewton.modifySearchDirection: 1.90e+00, 5.71e+01, 30x\n",
|
||||
" InexactGaussNewton.projection : 4.65e-05, 5.86e-03, 126x\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"opt.counter.summary()\n",
|
||||
"xc = opt.recall('xc')"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 12,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# import matplotlib.pyplot as plt\n",
|
||||
"# # plt.figure(1)\n",
|
||||
"# # for i in range(problem.G.shape[0]):\n",
|
||||
"# # plt.plot(problem.G[i,:])\n",
|
||||
"# meshPts = np.concatenate((mesh.gridN[0:1],np.kron(mesh.gridN[1::],np.ones(2))[:-1]))\n",
|
||||
"# modelPts = np.kron(1./model,np.ones(2,))\n",
|
||||
"# axM.semilogx(modelPts,meshPts,color=col)\n",
|
||||
"# plt.figure(2)\n",
|
||||
"# plt.plot(m1d.vectorCCx[active], np.log10(mappingExpAct*survey.mtrue)[active], 'b-')\n",
|
||||
"# plt.plot(m1d.vectorCCx[active], np.log10(mappingExpAct*mopt)[active], 'r-')\n",
|
||||
"# plt.show()\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -188,56 +540,209 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"execution_count": 13,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Runn the inversion\n",
|
||||
"mopt = inv.run(m_0)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"opt.counter.summary()\n",
|
||||
"xc = opt.recall('xc')"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"# plt.figure(1)\n",
|
||||
"# for i in range(problem.G.shape[0]):\n",
|
||||
"# plt.plot(problem.G[i,:])\n",
|
||||
"def plotMT1DModelData(problem,models,symList=None):\n",
|
||||
" # Make the analytic solution\n",
|
||||
" # \tdef makeAnalyticSolution(mesh,model,elev,freqs):\n",
|
||||
" # \t\tdata1D = []\n",
|
||||
" # \t\tfor freq in freqs:\n",
|
||||
" # \t\t\tanaEd, anaEu, anaHd, anaHu = simpegmt.Utils.MT1Danalytic.getEHfields(mesh,model,freq,elev)\n",
|
||||
" # \t\t\tanaE = anaEd+anaEu\n",
|
||||
" # \t\t\tanaH = anaHd+anaHu\n",
|
||||
" # \t\t\t# Scale the solution\n",
|
||||
" # \t\t\t# anaE = (anaEtemp/anaEtemp[-1])#.conj()\n",
|
||||
" # \t\t\t# anaH = (anaHtemp/anaEtemp[-1])#.conj()\n",
|
||||
" # \t\t\tanaZ = anaE/anaH\n",
|
||||
" # \t\t\t# Add to the list\n",
|
||||
" # \t\t\tdata1D.append((freq,0,0,elev,anaZ[0]))\n",
|
||||
" # \t\tdataRec = np.array(data1D,dtype=[('freq',float),('x',float),('y',float),('z',float),('zxx',complex)])\n",
|
||||
" # \t\treturn dataRec\n",
|
||||
" def appResPhs(freq,z):\n",
|
||||
" fr = simpeg.mkvc(freq,2)*np.ones(z.shape)\n",
|
||||
" app_res = ((1./(8e-7*np.pi**2))/fr)*np.abs(z)**2\n",
|
||||
" app_phs = np.arctan2(z.imag,z.real)*(180/np.pi)\n",
|
||||
" return app_res, app_phs\n",
|
||||
" \n",
|
||||
" # Setup the figure\n",
|
||||
" fontSize = 15\n",
|
||||
"\n",
|
||||
"plt.figure(2)\n",
|
||||
"plt.plot(m1d.vectorCCx, np.log10(survey.mtrue), 'b-')\n",
|
||||
"plt.plot(m1d.vectorCCx, np.log10(mopt), 'r-')\n",
|
||||
" fig = plt.figure(figsize=[9,7])\n",
|
||||
" axM = fig.add_axes([0.075,.1,.25,.875])\n",
|
||||
" axM.set_xlabel('Resistivity [Ohm*m]',fontsize=fontSize)\n",
|
||||
" axM.set_xlim(1e-1,1e5)\n",
|
||||
" axM.set_ylim(-10000,5000)\n",
|
||||
" axM.set_ylabel('Depth [km]',fontsize=fontSize)\n",
|
||||
" axR = fig.add_axes([0.42,.575,.5,.4])\n",
|
||||
" axR.set_xscale('log')\n",
|
||||
" axR.set_yscale('log')\n",
|
||||
" axR.invert_xaxis()\n",
|
||||
" # axR.set_xlabel('Frequency [Hz]')\n",
|
||||
" axR.set_ylabel('Apparent resistivity [Ohm m]',fontsize=fontSize)\n",
|
||||
"\n",
|
||||
" axP = fig.add_axes([0.42,.1,.5,.4])\n",
|
||||
" axP.set_xscale('log')\n",
|
||||
" axP.invert_xaxis()\n",
|
||||
" axP.set_ylim(0,90)\n",
|
||||
" axP.set_xlabel('Frequency [Hz]',fontsize=fontSize)\n",
|
||||
" axP.set_ylabel('Apparent phase [deg]',fontsize=fontSize)\n",
|
||||
"\n",
|
||||
" # if not symList:\n",
|
||||
" # \tsymList = ['x']*len(models)\n",
|
||||
" sys.path.append('/home/gudni/Dropbox/code/python/MTview')\n",
|
||||
" import plotDataTypes as pDt\n",
|
||||
" # Loop through the models.\n",
|
||||
" modelList = [problem.survey.mtrue]\n",
|
||||
" modelList.extend(models)\n",
|
||||
" if False:\n",
|
||||
" modelList = [problem.mapping.sigmaMap*mod for mod in modelList]\n",
|
||||
" for nr, model in enumerate(modelList):\n",
|
||||
" # Calculate the data\n",
|
||||
" if nr==0:\n",
|
||||
" data1D = problem.dataPair(problem.survey,problem.survey.dobs).toRecArray('Complex')\n",
|
||||
" else:\n",
|
||||
" data1D = problem.dataPair(problem.survey,problem.survey.dpred(model)).toRecArray('Complex')\n",
|
||||
" # Plot the data and the model \n",
|
||||
" colRat = nr/((len(modelList)-2)*1.)\n",
|
||||
" if colRat > 1.:\n",
|
||||
" col = 'k'\n",
|
||||
" else:\n",
|
||||
" col = plt.cm.seismic(1-colRat)\n",
|
||||
" # The model - make the pts to plot\n",
|
||||
" meshPts = np.concatenate((problem.mesh.gridN[0:1],np.kron(problem.mesh.gridN[1::],np.ones(2))[:-1]))\n",
|
||||
" modelPts = np.kron(1./(problem.mapping.sigmaMap*model),np.ones(2,))\n",
|
||||
" axM.semilogx(modelPts,meshPts,color=col)\n",
|
||||
"\n",
|
||||
" ## Data\n",
|
||||
" # Appres\n",
|
||||
" pDt.plotIsoStaImpedance(axR,np.array([0,0]),data1D,'zyx','res',pColor=col)\n",
|
||||
" # Appphs\n",
|
||||
" pDt.plotIsoStaImpedance(axP,np.array([0,0]),data1D,'zyx','phs',pColor=col)\n",
|
||||
" try:\n",
|
||||
" allData = np.concatenate((allData,mkvc(data1D['zyx'],2)),1)\n",
|
||||
" except:\n",
|
||||
" allData = simpeg.mkvc(data1D['zyx'],2)\n",
|
||||
" freq = data1D['freq']\n",
|
||||
" res, phs = appResPhs(freq,allData)\n",
|
||||
"\n",
|
||||
" stdCol = 'gray'\n",
|
||||
" axRtw = axR.twinx()\n",
|
||||
" axRtw.set_ylabel('Std of log10',color=stdCol)\n",
|
||||
" [(t.set_color(stdCol), t.set_rotation(-45)) for t in axRtw.get_yticklabels()]\n",
|
||||
" axPtw = axP.twinx()\n",
|
||||
" axPtw.set_ylabel('Std ',color=stdCol)\n",
|
||||
" [t.set_color(stdCol) for t in axPtw.get_yticklabels()]\n",
|
||||
" axRtw.plot(freq, np.std(np.log10(res),1),'--',color=stdCol)\n",
|
||||
" axPtw.plot(freq, np.std(phs,1),'--',color=stdCol)\n",
|
||||
"\n",
|
||||
" # Fix labels and ticks\n",
|
||||
"\n",
|
||||
" yMtick = [l/1000 for l in axM.get_yticks().tolist()]\n",
|
||||
" axM.set_yticklabels(yMtick)\n",
|
||||
" [ l.set_rotation(90) for l in axM.get_yticklabels()]\n",
|
||||
" [ l.set_rotation(90) for l in axR.get_yticklabels()]\n",
|
||||
" [(t.set_color(stdCol), t.set_rotation(-45)) for t in axRtw.get_yticklabels()]\n",
|
||||
" [t.set_color(stdCol) for t in axPtw.get_yticklabels()]\n",
|
||||
" for ax in [axM,axR,axP]:\n",
|
||||
" ax.xaxis.set_tick_params(labelsize=fontSize)\n",
|
||||
" ax.yaxis.set_tick_params(labelsize=fontSize)\n",
|
||||
" return fig\n",
|
||||
"# plotMT1DModelData(problem,[mopt])"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 14,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stderr",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"/home/gudni/anaconda/lib/python2.7/site-packages/numpy/ma/core.py:2834: FutureWarning: Numpy has detected that you (may be) writing to an array returned\n",
|
||||
"by numpy.diagonal or by selecting multiple fields in a record\n",
|
||||
"array. This code will likely break in a future numpy release --\n",
|
||||
"see numpy.diagonal or arrays.indexing reference docs for details.\n",
|
||||
"The quick fix is to make an explicit copy (e.g., do\n",
|
||||
"arr.diagonal().copy() or arr[['f0','f1']].copy()).\n",
|
||||
" if (obj.__array_interface__[\"data\"][0]\n",
|
||||
"/home/gudni/anaconda/lib/python2.7/site-packages/numpy/ma/core.py:2835: FutureWarning: Numpy has detected that you (may be) writing to an array returned\n",
|
||||
"by numpy.diagonal or by selecting multiple fields in a record\n",
|
||||
"array. This code will likely break in a future numpy release --\n",
|
||||
"see numpy.diagonal or arrays.indexing reference docs for details.\n",
|
||||
"The quick fix is to make an explicit copy (e.g., do\n",
|
||||
"arr.diagonal().copy() or arr[['f0','f1']].copy()).\n",
|
||||
" != self.__array_interface__[\"data\"][0]):\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"plotMT1DModelData(problem,[mopt,mopt2])\n",
|
||||
"plt.show()\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"execution_count": 15,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"mopt"
|
||||
"modelList = [problem.survey.mtrue]\n",
|
||||
"modelList.extend([mopt])\n",
|
||||
"# problem.mapping.sigmaMap*mopt"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 16,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"[['zyxr'], ['zyxi']]"
|
||||
]
|
||||
},
|
||||
"execution_count": 16,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"src = survey.srcList[0]\n",
|
||||
"[[rx.rxType.replace('z1d','zyx')] for rx in src.rxList ]"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 17,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"'zyxr'"
|
||||
]
|
||||
},
|
||||
"execution_count": 17,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"'z1dr'.replace('z1d','zyx')"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -255,6 +760,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,
|
||||
|
||||
+10
-10
@@ -24,7 +24,7 @@ class BaseMTProblem(BaseFDEMProblem):
|
||||
# Use the forward and devs from BaseFDEMProblem
|
||||
# Might need to add more stuff here.
|
||||
|
||||
def Jvec(self, m, v, f=None):
|
||||
def Jvec(self, m, v, u=None):
|
||||
"""
|
||||
Function to calculate the data sensitivities dD/dm times a vector.
|
||||
|
||||
@@ -36,8 +36,8 @@ class BaseMTProblem(BaseFDEMProblem):
|
||||
"""
|
||||
|
||||
# Calculate the fields
|
||||
if f is None:
|
||||
f = self.fields(m)
|
||||
if u is None:
|
||||
u = self.fields(m)
|
||||
# Set current model
|
||||
self.curModel = m
|
||||
# Initiate the Jv object
|
||||
@@ -52,7 +52,7 @@ class BaseMTProblem(BaseFDEMProblem):
|
||||
# We need fDeriv_m = df/du*du/dm + df/dm
|
||||
# Construct du/dm, it requires a solve
|
||||
ftype = self._fieldType + 'Solution'
|
||||
u_src = f[src, ftype]
|
||||
u_src = u[src, ftype]
|
||||
dA_dm = self.getADeriv_m(freq, u_src, v)
|
||||
dRHS_dm = self.getRHSDeriv_m(freq, v)
|
||||
if dRHS_dm is None:
|
||||
@@ -62,14 +62,14 @@ class BaseMTProblem(BaseFDEMProblem):
|
||||
# Calculate the projection derivatives
|
||||
for rx in src.rxList:
|
||||
# Get the projection derivative
|
||||
PDeriv = lambda v: rx.projectFieldsDeriv(src, self.mesh, f, v) # wrt u, also have wrt m
|
||||
PDeriv = lambda v: rx.projectFieldsDeriv(src, self.mesh, u, v) # wrt u, also have wrt m
|
||||
Jv[src, rx] = PDeriv(du_dm)
|
||||
# Return the vectorized sensitivities
|
||||
return mkvc(Jv)
|
||||
|
||||
def Jtvec(self, m, v, f=None):
|
||||
if f is None:
|
||||
f = self.fields(m)
|
||||
def Jtvec(self, m, v, u=None):
|
||||
if u is None:
|
||||
u = self.fields(m)
|
||||
|
||||
self.curModel = m
|
||||
|
||||
@@ -85,11 +85,11 @@ class BaseMTProblem(BaseFDEMProblem):
|
||||
|
||||
for src in self.survey.getSrcByFreq(freq):
|
||||
ftype = self._fieldType + 'Solution'
|
||||
u_src = f[src, ftype]
|
||||
u_src = u[src, ftype]
|
||||
|
||||
for rx in src.rxList:
|
||||
# Get the adjoint projectFieldsDeriv
|
||||
PTv = rx.projectFieldsDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt u, need possibility wrt m
|
||||
PTv = rx.projectFieldsDeriv(src, self.mesh, u, v[src, rx], adjoint=True) # wrt u, need possibility wrt m
|
||||
# Get the
|
||||
dA_duIT = ATinv * PTv
|
||||
dA_dmT = self.getADeriv_m(freq, u_src, dA_duIT, adjoint=True)
|
||||
|
||||
+9
-5
@@ -8,9 +8,9 @@ from numpy.lib import recfunctions as recFunc
|
||||
############
|
||||
class DataMT(Survey.Data):
|
||||
'''
|
||||
Data class for MTdata
|
||||
Data class for MTdata
|
||||
|
||||
:param SimPEG survey object survey:
|
||||
:param SimPEG survey object survey:
|
||||
:param v vector with data
|
||||
|
||||
'''
|
||||
@@ -37,14 +37,18 @@ class DataMT(Survey.Data):
|
||||
# Note: needs to be written more generally, using diffterent rxTypes and not all the data at the locaitons
|
||||
# Assume the same locs for all RX
|
||||
locs = src.rxList[0].locs
|
||||
if locs.shape[1] == 1:
|
||||
locs = np.hstack((np.array([[0.0,0.0]]),locs))
|
||||
elif locs.shape[1] == 2:
|
||||
locs = np.hstack((np.array([[0.0]]),locs))
|
||||
tArrRec = np.concatenate((src.freq*np.ones((locs.shape[0],1)),locs,np.nan*np.ones((locs.shape[0],8))),axis=1).view(dtRI)
|
||||
# np.array([(src.freq,rx.locs[0,0],rx.locs[0,1],rx.locs[0,2],np.nan ,np.nan ,np.nan ,np.nan ,np.nan ,np.nan ,np.nan ,np.nan ) for rx in src.rxList],dtype=dtRI)
|
||||
# Get the type and the value for the DataMT object as a list
|
||||
typeList = [[rx.rxType,self[src,rx]] for rx in src.rxList]
|
||||
typeList = [[rx.rxType.replace('z1d','zyx'),self[src,rx]] for rx in src.rxList]
|
||||
# Insert the values to the temp array
|
||||
for nr,(key,val) in enumerate(typeList):
|
||||
tArrRec[key] = mkvc(val,2)
|
||||
# Masked array
|
||||
# Masked array
|
||||
mArrRec = np.ma.MaskedArray(rec2ndarr(tArrRec),mask=np.isnan(rec2ndarr(tArrRec))).view(dtype=tArrRec.dtype)
|
||||
# Unique freq and loc of the masked array
|
||||
uniFLmarr = np.unique(mArrRec[['freq','x','y','z']])
|
||||
@@ -66,5 +70,5 @@ class DataMT(Survey.Data):
|
||||
for comp in ['zxx','zxy','zyx','zyy']:
|
||||
outArr[comp] = outTemp[comp+'r'].copy() + 1j*outTemp[comp+'i'].copy()
|
||||
|
||||
# Return
|
||||
# Return
|
||||
return outArr
|
||||
|
||||
+15
-1
@@ -118,4 +118,18 @@ class FieldsMT_1D(FieldsMT):
|
||||
|
||||
This function stacks the fields derivatives appropriately
|
||||
"""
|
||||
return None
|
||||
return None
|
||||
|
||||
class FieldsMT_3D(FieldsMT):
|
||||
"""
|
||||
Fields storage for the 3D MT solution.
|
||||
"""
|
||||
knownFields = {'e_px':'E','e_py':'E','b_px':'F','b_py':'F'}
|
||||
aliasFields = { }
|
||||
# 'e_1d' : ['e_1dSolution','F','_e'],
|
||||
# 'e_1dPrimary' : ['e_1dSolution','F','_ePrimary'],
|
||||
# 'e_1dSecondary' : ['e_1dSolution','F','_eSecondary'],
|
||||
# 'b_1d' : ['e_1dSolution','E','_b'],
|
||||
# 'b_1dPrimary' : ['e_1dSolution','E','_bPrimary'],
|
||||
# 'b_1dSecondary' : ['e_1dSolution','E','_bSecondary']
|
||||
# }
|
||||
@@ -108,7 +108,7 @@ class eForm_psField(BaseMTProblem):
|
||||
# Store the fields
|
||||
Src = self.survey.getSrcByFreq(freq)[0]
|
||||
# NOTE: only store the e_solution(secondary), all other components calculated in the fields object
|
||||
F[Src, 'e_1dSolution'] = e_s[:,1] # Only storing the yx polarization as 1d
|
||||
F[Src, 'e_1dSolution'] = e_s[:,-1] # Only storing the yx polarization as 1d
|
||||
|
||||
# Note curl e = -iwb so b = -curl e /iw
|
||||
# b = -( self.mesh.nodalGrad * e )/( 1j*omega(freq) )
|
||||
@@ -195,7 +195,7 @@ class eForm_TotalField(BaseMTProblem):
|
||||
eBC = np.r_[Etot[0],Etot[-1]]
|
||||
# The right hand side
|
||||
|
||||
return Aio*eBC, eBC
|
||||
return -Aio*eBC, eBC
|
||||
|
||||
def getRHSderiv(self, freq, backSigma, u, v, adjoint=False):
|
||||
raise NotImplementedError('getRHSDeriv not implemented yet')
|
||||
@@ -211,7 +211,7 @@ class eForm_TotalField(BaseMTProblem):
|
||||
self.curModel = m
|
||||
# RHS, CalcFields = self.getRHS(freq,m_back), self.calcFields
|
||||
|
||||
F = FieldsMT(self.mesh, self.survey)
|
||||
F = FieldsMT_1D(self.mesh, self.survey)
|
||||
for freq in self.survey.freqs:
|
||||
if self.verbose:
|
||||
startTime = time.time()
|
||||
@@ -224,14 +224,8 @@ class eForm_TotalField(BaseMTProblem):
|
||||
e = mkvc(np.r_[e_o[0], e_i, e_o[1]],2)
|
||||
# Store the fields
|
||||
Src = self.survey.getSrcByFreq(freq)
|
||||
# Store the fields
|
||||
# NOTE: only store
|
||||
F[Src, 'e_1d'] = e
|
||||
# F[Src, 'e_py'] = 0*e[:,0]
|
||||
# Note curl e = -iwb so b = -curl e /iw
|
||||
b = -( self.mesh.nodalGrad * e )/( 1j*omega(freq) )
|
||||
# F[Src, 'b_px'] = 0*b[:,0]
|
||||
F[Src, 'b_1d'] = b
|
||||
# NOTE: only store e fields
|
||||
F[Src, 'e_1dSolution'] = e[:,0]
|
||||
if self.verbose:
|
||||
print 'Ran for {:f} seconds'.format(time.time()-startTime)
|
||||
sys.stdout.flush()
|
||||
|
||||
@@ -3,7 +3,7 @@ from simpegEM.Utils.EMUtils import omega
|
||||
from scipy.constants import mu_0
|
||||
from simpegMT.BaseMT import BaseMTProblem
|
||||
from simpegMT.SurveyMT import SurveyMT
|
||||
from simpegMT.FieldsMT import FieldsMT
|
||||
from simpegMT.FieldsMT import FieldsMT_3D
|
||||
from simpegMT.DataMT import DataMT
|
||||
import multiprocessing, sys, time
|
||||
|
||||
@@ -22,7 +22,7 @@ class eForm_ps(BaseMTProblem):
|
||||
# From FDEMproblem: Used to project the fields. Currently not used for MTproblem.
|
||||
_fieldType = 'e'
|
||||
_eqLocs = 'FE'
|
||||
|
||||
fieldsPair = FieldsMT_3D
|
||||
# Need to add the src ....
|
||||
|
||||
|
||||
@@ -112,7 +112,7 @@ class eForm_ps(BaseMTProblem):
|
||||
# Set the current model
|
||||
self.curModel = m
|
||||
|
||||
F = FieldsMT(self.mesh, self.survey)
|
||||
F = FieldsMT_3D(self.mesh, self.survey)
|
||||
for freq in self.survey.freqs:
|
||||
if self.verbose:
|
||||
startTime = time.time()
|
||||
@@ -149,7 +149,7 @@ class eForm_Tp(BaseMTProblem):
|
||||
|
||||
_fieldType = 'e'
|
||||
_eqLocs = 'FE'
|
||||
fieldsPair = FieldsMT
|
||||
fieldsPair = FieldsMT_3D
|
||||
|
||||
# Set new properties
|
||||
# Background model
|
||||
@@ -242,7 +242,7 @@ class eForm_Tp(BaseMTProblem):
|
||||
self.backModel = m_back
|
||||
# RHS, CalcFields = self.getRHS(freq,m_back), self.calcFields
|
||||
|
||||
F = FieldsMT(self.mesh, self.survey)
|
||||
F = FieldsMT_3D(self.mesh, self.survey)
|
||||
for freq in self.survey.freqs:
|
||||
if self.verbose:
|
||||
startTime = time.time()
|
||||
|
||||
@@ -23,8 +23,8 @@ def homo1DModelSource(mesh,freq,sigma_1d):
|
||||
# # Note: Everything is using e^iwt
|
||||
e0_1d = get1DEfields(mesh1d,sigma_1d,freq)
|
||||
if mesh.dim == 1:
|
||||
eBG_px = -simpeg.mkvc(e0_1d,2)
|
||||
eBG_py = simpeg.mkvc(e0_1d,2)
|
||||
eBG_px = simpeg.mkvc(e0_1d,2)
|
||||
eBG_py = -simpeg.mkvc(e0_1d,2) # added a minus to make the results in the correct quadrents.
|
||||
elif mesh.dim == 2:
|
||||
ex_px = np.zeros(mesh.vnEx,dtype=complex)
|
||||
ey_px = np.zeros((mesh.nEy,1),dtype=complex)
|
||||
|
||||
@@ -93,7 +93,9 @@ 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
|
||||
f_part_complex = ex/bx
|
||||
# Note: Has a minus sign in front, to comply with quadrant calculations.
|
||||
# Can be derived from zyx case for the 3D case.
|
||||
f_part_complex = -ex/bx
|
||||
# elif self.projType is 'Z2D':
|
||||
elif self.projType is 'Z3D':
|
||||
# Get the projection
|
||||
@@ -146,8 +148,8 @@ 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 = 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)
|
||||
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)
|
||||
elif self.projType is 'Z2D':
|
||||
raise NotImplementedError('Has not be implement for 2D impedance tensor')
|
||||
@@ -162,8 +164,8 @@ 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_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_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_real = np.sum(np.hstack((dP_deTv,dP_dbTv)),1)
|
||||
elif self.projType is 'Z2D':
|
||||
@@ -323,7 +325,6 @@ class SurveyMT(Survey.BaseSurvey):
|
||||
def projectFields(self, u):
|
||||
data = DataMT(self)
|
||||
for src in self.srcList:
|
||||
print 'Project at freq: {:.3e}'.format(src.freq)
|
||||
sys.stdout.flush()
|
||||
for rx in src.rxList:
|
||||
data[src, rx] = rx.projectFields(src, self.mesh, u)
|
||||
|
||||
@@ -45,7 +45,7 @@ def getAppResPhs(MTdata):
|
||||
# Make impedance
|
||||
def appResPhs(freq,z):
|
||||
app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
|
||||
app_phs = np.arctan2(-z.imag,z.real)*(180/np.pi)
|
||||
app_phs = np.arctan2(z.imag,z.real)*(180/np.pi)
|
||||
return app_res, app_phs
|
||||
zList = []
|
||||
for src in MTdata.survey.srcList:
|
||||
@@ -93,7 +93,7 @@ def appPhs_TotalFieldNorm(sigmaHalf):
|
||||
# Calculate the app phs
|
||||
app_p = np.array(getAppResPhs(data))[:,1]
|
||||
|
||||
return np.linalg.norm(np.abs(app_p - np.ones(survey.nFreq)*135)/ 135)
|
||||
return np.linalg.norm(np.abs(app_p - np.ones(survey.nFreq)*45)/ 45)
|
||||
|
||||
def appRes_psFieldNorm(sigmaHalf):
|
||||
|
||||
@@ -129,7 +129,7 @@ def appPhs_psFieldNorm(sigmaHalf):
|
||||
# Calculate the app phs
|
||||
app_p = np.array(getAppResPhs(data))[:,1]
|
||||
|
||||
return np.linalg.norm(np.abs(app_p - np.ones(survey.nFreq)*135)/ 135)
|
||||
return np.linalg.norm(np.abs(app_p - np.ones(survey.nFreq)*45)/ 45)
|
||||
|
||||
class TestAnalytics(unittest.TestCase):
|
||||
|
||||
|
||||
@@ -0,0 +1,155 @@
|
||||
import unittest
|
||||
import SimPEG as simpeg
|
||||
import simpegMT as simpegmt
|
||||
from SimPEG.Utils import meshTensor
|
||||
import numpy as np
|
||||
# Define the tolerances
|
||||
TOLr = 5e-2
|
||||
TOLp = 5e-2
|
||||
|
||||
|
||||
def setupSurvey(sigmaHalf,tD=True):
|
||||
|
||||
# Frequency
|
||||
nFreq = 33
|
||||
freqs = np.logspace(3,-3,nFreq)
|
||||
# Make the mesh
|
||||
ct = 5
|
||||
air = meshTensor([(ct,25,1.3)])
|
||||
# coreT0 = meshTensor([(ct,15,1.2)])
|
||||
# coreT1 = np.kron(meshTensor([(coreT0[-1],15,1.3)]),np.ones((7,)))
|
||||
core = np.concatenate( ( np.kron(meshTensor([(ct,15,-1.2)]),np.ones((10,))) , meshTensor([(ct,20)]) ) )
|
||||
bot = meshTensor([(core[0],15,-1.3)])
|
||||
x0 = -np.array([np.sum(np.concatenate((core,bot)))])
|
||||
m1d = simpeg.Mesh.TensorMesh([np.concatenate((bot,core,air))], x0=x0)
|
||||
# Make the model
|
||||
sigma = np.zeros(m1d.nC) + sigmaHalf
|
||||
sigma[m1d.gridCC > 0 ] = 1e-8
|
||||
sigmaBack = sigma.copy()
|
||||
# Add structure
|
||||
shallow = (m1d.gridCC < -200) * (m1d.gridCC > -600)
|
||||
deep = (m1d.gridCC < -3000) * (m1d.gridCC > -5000)
|
||||
sigma[shallow] = 1
|
||||
sigma[deep] = 0.1
|
||||
|
||||
rxList = []
|
||||
for rxType in ['z1dr','z1di']:
|
||||
rxList.append(simpegmt.SurveyMT.RxMT(simpeg.mkvc(np.array([0.0]),2).T,rxType))
|
||||
# Source list
|
||||
srcList =[]
|
||||
if tD:
|
||||
for freq in freqs:
|
||||
srcList.append(simpegmt.SurveyMT.srcMT_polxy_1DhomotD(rxList,freq))
|
||||
else:
|
||||
for freq in freqs:
|
||||
srcList.append(simpegmt.SurveyMT.srcMT_polxy_1Dprimary(rxList,freq,sigmaBack))
|
||||
|
||||
survey = simpegmt.SurveyMT.SurveyMT(srcList)
|
||||
return survey, sigma, m1d
|
||||
|
||||
def getAppResPhs(MTdata):
|
||||
# Make impedance
|
||||
def appResPhs(freq,z):
|
||||
app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
|
||||
app_phs = np.arctan2(z.imag,z.real)*(180/np.pi)
|
||||
return app_res, app_phs
|
||||
zList = []
|
||||
for src in MTdata.survey.srcList:
|
||||
zc = [src.freq]
|
||||
for rx in src.rxList:
|
||||
if 'i' in rx.rxType:
|
||||
m=1j
|
||||
else:
|
||||
m = 1
|
||||
zc.append(m*MTdata[src,rx])
|
||||
zList.append(zc)
|
||||
return [appResPhs(zList[i][0],np.sum(zList[i][1:3])) for i in np.arange(len(zList))]
|
||||
|
||||
def calculateAnalyticSolution(srcList,mesh,model):
|
||||
surveyAna = simpegmt.SurveyMT.SurveyMT(srcList)
|
||||
data1D = simpegmt.DataMT.DataMT(surveyAna)
|
||||
for src in surveyAna.srcList:
|
||||
elev = src.rxList[0].locs[0]
|
||||
anaEd, anaEu, anaHd, anaHu = simpegmt.Utils.MT1Danalytic.getEHfields(mesh,model,src.freq,elev)
|
||||
anaE = anaEd+anaEu
|
||||
anaH = anaHd+anaHu
|
||||
# Scale the solution
|
||||
# anaE = (anaEtemp/anaEtemp[-1])#.conj()
|
||||
# anaH = (anaHtemp/anaEtemp[-1])#.conj()
|
||||
anaZ = anaE/anaH
|
||||
for rx in src.rxList:
|
||||
data1D[src,rx] = getattr(anaZ, rx.projComp)
|
||||
return data1D
|
||||
|
||||
def dataMis_AnalyticTotalDomain(sigmaHalf):
|
||||
|
||||
# Make the survey
|
||||
|
||||
# Total domain solution
|
||||
surveyTD, sigma, mesh = setupSurvey(sigmaHalf)
|
||||
problemTD = simpegmt.ProblemMT1D.eForm_TotalField(mesh)
|
||||
problemTD.pair(surveyTD)
|
||||
# Analytic data
|
||||
dataAnaObj = calculateAnalyticSolution(surveyTD.srcList,mesh,sigma)
|
||||
# dataTDObj = simpegmt.DataMT.DataMT(surveyTD, surveyTD.dpred(sigma))
|
||||
dataTD = surveyTD.dpred(sigma)
|
||||
dataAna = simpeg.mkvc(dataAnaObj)
|
||||
return np.all((dataTD - dataAna)/dataAna < 2.)
|
||||
# surveyTD.dtrue = -simpeg.mkvc(dataAna,2)
|
||||
# surveyTD.dobs = -simpeg.mkvc(dataAna,2)
|
||||
# surveyTD.Wd = np.ones(surveyTD.dtrue.shape) #/(np.abs(surveyTD.dtrue)*0.01)
|
||||
# # Setup the data misfit
|
||||
# dmis = simpeg.DataMisfit.l2_DataMisfit(surveyTD)
|
||||
# dmis.Wd = surveyTD.Wd
|
||||
# return dmis.eval(sigma)
|
||||
|
||||
|
||||
def dataMis_AnalyticPrimarySecondary(sigmaHalf):
|
||||
|
||||
# Make the survey
|
||||
# Primary secondary
|
||||
surveyPS, sigmaPS, mesh = setupSurvey(sigmaHalf,False)
|
||||
problemPS = simpegmt.ProblemMT1D.eForm_psField(mesh)
|
||||
problemPS.pair(surveyPS)
|
||||
# Analytic data
|
||||
dataAna = calculateAnalyticSolution(surveyPS.srcList,mesh,sigma)
|
||||
|
||||
surveyPS.dtrue = dataAna
|
||||
# Project the data
|
||||
data = surveyPS.dpred(sigmaPS)
|
||||
|
||||
# Setup the data misfit
|
||||
dmis = simpeg.DataMisfit.l2_DataMisfit(survey)
|
||||
return dmis.eval(sigma)
|
||||
|
||||
|
||||
|
||||
class TestNumericVsAnalytics(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
pass
|
||||
# Total Fields
|
||||
# def test_appRes2en1(self):self.assertLess(appRes_TotalFieldNorm(2e-1), TOLr)
|
||||
# def test_appPhs2en1(self):self.assertLess(appPhs_TotalFieldNorm(2e-1), TOLp)
|
||||
|
||||
def test_appRes2en2(self):self.assertTrue(dataMis_AnalyticTotalDomain(2e-2))
|
||||
# def test_appPhs2en2(self):self.assert(appPhs_TotalFieldNorm(2e-2), TOLp)
|
||||
|
||||
# def test_appRes2en3(self):self.assertLess(appRes_TotalFieldNorm(2e-3), TOLr)
|
||||
# def test_appPhs2en3(self):self.assertLess(appPhs_TotalFieldNorm(2e-3), TOLp)
|
||||
|
||||
# def test_appRes2en4(self):self.assertLess(appRes_TotalFieldNorm(2e-4), TOLr)
|
||||
# def test_appPhs2en4(self):self.assertLess(appPhs_TotalFieldNorm(2e-4), TOLp)
|
||||
|
||||
# def test_appRes2en5(self):self.assertLess(appRes_TotalFieldNorm(2e-5), TOLr)
|
||||
# def test_appPhs2en5(self):self.assertLess(appPhs_TotalFieldNorm(2e-5), TOLp)
|
||||
|
||||
# def test_appRes2en6(self):self.assertLess(appRes_TotalFieldNorm(2e-6), TOLr)
|
||||
# def test_appPhs2en6(self):self.assertLess(appPhs_TotalFieldNorm(2e-6), TOLp)
|
||||
|
||||
# Primary/secondary
|
||||
# def test_appRes2en2_ps(self):self.assertLess(appRes_psFieldNorm(2e-2), TOLr)
|
||||
# def test_appPhs2en2_ps(self):self.assertLess(appPhs_psFieldNorm(2e-2), TOLp)
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
@@ -53,7 +53,7 @@ def getEHfields(m1d,sigma,freq,zd):
|
||||
|
||||
# Loop over the layers and calculate the fields
|
||||
# In the halfspace below the mesh
|
||||
dup = m1d.vectorNx[0]
|
||||
dup = m1d.vectorNx[0]
|
||||
dind = dup >= zd
|
||||
Ed[dind] = UDp[1,0]*np.exp(-1j*k[0]*(dup-zd[dind]))
|
||||
Eu[dind] = UDp[0,0]*np.exp(1j*k[0]*(dup-zd[dind]))
|
||||
|
||||
Reference in New Issue
Block a user