mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-07 12:34:04 +08:00
Adding a MT1D inversion test notebook
This commit is contained in:
GudniRos
@@ -0,0 +1,262 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import SimPEG as simpeg\n",
|
||||
"import simpegMT as simpegmt\n",
|
||||
"import numpy as np"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"## Setup the problem\n",
|
||||
"\n",
|
||||
"# Frequency\n",
|
||||
"nFreq = 33\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",
|
||||
"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",
|
||||
"m1d = simpeg.Mesh.TensorMesh([np.concatenate((bot,core,air))], x0=x0)\n",
|
||||
"\n",
|
||||
"## Setup model varibles\n",
|
||||
"active = m1d.vectorCCx<0.\n",
|
||||
"layer = (m1d.vectorCCx<-200.) & (m1d.vectorCCx>=-600.)\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",
|
||||
"# 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",
|
||||
"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",
|
||||
"\n",
|
||||
"# Receivers \n",
|
||||
"# 1D impedance at the surface (elevation 0)\n",
|
||||
"rxList = []\n",
|
||||
"for rxType in ['z1dr','z1di']:\n",
|
||||
" rxList.append(simpegmt.SurveyMT.RxMT(simpeg.mkvc(np.array([0.0]),2).T,rxType))\n",
|
||||
"# Source list\n",
|
||||
"srcList =[]\n",
|
||||
"tD = False\n",
|
||||
"if tD:\n",
|
||||
" for freq in freqs:\n",
|
||||
" 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_0))\n",
|
||||
"# Make the survey\n",
|
||||
"survey = simpegmt.SurveyMT.SurveyMT(srcList)\n",
|
||||
"survey.mtrue = sigma_true\n",
|
||||
"# Set the problem\n",
|
||||
"problem = simpegmt.ProblemMT1D.eForm_psField(m1d,mapping=mappingExpAct)\n",
|
||||
"problem.pair(survey)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"180"
|
||||
]
|
||||
},
|
||||
"execution_count": 3,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"np.sum(active)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"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"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"## Make the observed data \n",
|
||||
"# Project the data\n",
|
||||
"d_true = survey.dpred(m_true)\n",
|
||||
"survey.dtrue = d_true\n",
|
||||
"# Add noise\n",
|
||||
"std = 0.05 # 5% std\n",
|
||||
"noise = std*abs(survey.dtrue)*np.random.randn(*survey.dtrue.shape)\n",
|
||||
"# Assign the dobs\n",
|
||||
"survey.dobs = survey.dtrue + noise\n",
|
||||
"survey.std = survey.dobs*0 + std\n",
|
||||
"# Assign the data weight\n",
|
||||
"survey.Wd = 1/(abs(survey.dobs)*std)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"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.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",
|
||||
"reg.alpha_s = 1e-5\n",
|
||||
"reg.alpha_x = 1.\n",
|
||||
"\n",
|
||||
"# Inversion problem\n",
|
||||
"invProb = simpeg.InvProblem.BaseInvProblem(dmis, reg, opt)\n",
|
||||
"invProb.counter = C\n",
|
||||
"# Beta cooling\n",
|
||||
"beta = simpeg.Directives.BetaSchedule()\n",
|
||||
"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"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"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",
|
||||
"\n",
|
||||
"plt.figure(2)\n",
|
||||
"plt.plot(m1d.vectorCCx, np.log10(survey.mtrue), 'b-')\n",
|
||||
"plt.plot(m1d.vectorCCx, np.log10(mopt), 'r-')\n",
|
||||
"plt.show()\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"mopt"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Python 2",
|
||||
"language": "python",
|
||||
"name": "python2"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 0
|
||||
}
|
||||
Reference in New Issue
Block a user