mirror of
https://github.com/wassname/simpeg.git
synced 2026-10-28 13:30:34 +08:00
Wrote 1D codes, both analytic and FV solutions. \n Made test notebooks for 1D and worked on the 3D notebook
This commit is contained in:
GudniRos committed
1 parent
828960c29f
commit
47e7132f26
19 files changed
+586
-52
No files matched your search
@@ -0,0 +1,191 @@
|
||||
{
|
||||
"metadata": {
|
||||
"name": "",
|
||||
"signature": "sha256:f84d7f6fac432ae58377d9da7313463a2286c599c51cd1ac6b9a6ad5cb6ebd91"
|
||||
},
|
||||
"nbformat": 3,
|
||||
"nbformat_minor": 0,
|
||||
"worksheets": [
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"# Test 1D solution of MT problem and compare to a analytic solution\n"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": 1
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"# import the simpegMT module\n",
|
||||
"from simpegMT.Utils import MT1Danalytic, MT1Dsolutions\n",
|
||||
"\n",
|
||||
"def omega(freq):\n",
|
||||
" \"\"\"Change frequency to angular frequency, omega\"\"\"\n",
|
||||
" return 2.*np.pi*freq"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": 2
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"# Set up the mesh.\n",
|
||||
"freq = 10.0\n",
|
||||
"hz = [(100.,50)]\n",
|
||||
"M = simpeg.Mesh.TensorMesh([hz],'C')\n",
|
||||
"sig = np.zeros(M.nC) + 1e-8\n",
|
||||
"sig[M.vectorCCx<=0] = 1.0\n"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": 3
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"# Get the fields\n",
|
||||
"anaEd, anaEu, anaHd, anaHu = MT1Danalytic.getEHfields(M,sig,freq,M.vectorNx)\n",
|
||||
"anaEtemp = anaEd+anaEu\n",
|
||||
"# Scale the solution\n",
|
||||
"anaE = anaEtemp/anaEtemp[-1]\n",
|
||||
"\n",
|
||||
"solE = MT1Dsolutions.get1DEfields(M,sig,freq)\n"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": 4
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"plot(solE.real,M.vectorNx,'r*--',anaE.real,M.vectorNx,'b+:')"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"metadata": {},
|
||||
"output_type": "pyout",
|
||||
"prompt_number": 5,
|
||||
"text": [
|
||||
"[<matplotlib.lines.Line2D at 0x7f5b370f1a50>,\n",
|
||||
" <matplotlib.lines.Line2D at 0x7f5b370f1cd0>]"
|
||||
]
|
||||
},
|
||||
{
|
||||
"metadata": {},
|
||||
"output_type": "display_data",
|
||||
"png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEACAYAAACpoOGTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VOX1x/HPYYeyiwWJIFaIJahgAUGxmogKWBWqVVDA\nBSoWfi7V1gVtBWoVl4pof2LjUkRQEcWyuLBoCVp+Cq4oAjqo2BggIqCgICRwfn/MnTDAELLPTfi+\nX6+8uPPce2eeKzgnz3Oee665OyIiInurluwOiIhIOClAiIhIQgoQIiKSkAKEiIgkpAAhIiIJKUCI\niEhCpQoQZlbHzBab2QdmttzMxgbtTc1svpl9ambzzKxx3DkjzSxiZivN7My49s5m9lGw74HS9EtE\nREqvVAHC3X8EMty9E3AckGFmJwM3A/PdPRV4LXiNmaUB/YE0oDcwwcwseLuHgaHu3g5oZ2a9S9M3\nEREpnVJPMbn71mCzFlAd2AScC0wK2icB/YLtvsAz7p7n7quBVUA3MzsMaODuS4Ljnow7R0REkqDU\nAcLMqpnZB0AusMDdPwaau3tucEgu0DzYbgl8FXf6V0BKgvacoF1ERJKkRmnfwN13AZ3MrBEw18wy\n9trvZqZ6HiIilUypA0SMu39nZi8BnYFcM2vh7uuC6aOvg8NygFZxpx1OdOSQE2zHt+fs/RkKNCIi\nJePuduCj9lTaVUzNYiuUzKwucAbwPjALuDQ47FJgRrA9CxhgZrXM7EigHbDE3dcBm82sW5C0Hhx3\nzh7cvcr+jBo1Kul90PXp+g62azsYrq+kSjuCOAyYZGbViAabye7+mpm9D0wzs6HAauDC4Mt9uZlN\nA5YD+cAI3937EcATQF3gZXefU8q+iYhIKZQqQLj7R8AvErRvBE7fzzl3AncmaH8XOLY0/RERkbKj\nO6lDJD09PdldKFe6vsqrKl8bVP3rKykrzfxURTMzr0z9FREJAzPDKzpJLSIiVZcChIiIJKQAISIi\nCSlAiIhIQgoQIiKSkAKEiIgkpAAhIiIJKUCIiEhCChAiIpKQAoSISBVWmuoTChAiIlXY3OnTS3xu\nmT0wSEREwmNKZiZTH3yQn21pUOL30AhCRKSK2bYNUlKH0ShtGvPW3Fji91GAEBGpZNydKy9+ZM/8\nQn4+AG++Cf36QUaGMbj/CvrUu6zEn6MAISJSycydPp2l0zcwb+JEfhj/KBccmsWugYMBOOEEePHF\n6HHZkQi9J04s8efoeRAiIpXElMxMpo4fT401p/LzzUdQ3W5laYMGdPjl7dw57Uqq16ud8LySPg9C\nSWoRkZDbujU6gxTNK5zEfz7dwEzS6VG/Nm169ab38PZUr1fs7/8D0hSTiEgI7JNX+PHHaLYZuPXW\n6LRRRoYx+KJPOO8n59K92UN0YTSD+68gI6PsgwNoBCEiEgoFeYU2t5L/1iH85/+qMfb5dnD22Ywb\nBxbEgFheoeGH53HScS3IjkTKrU/KQYiIJNGUzEwm3fk41b8+mW4/NiS/zt283eBIzv7tVVxz54j9\nnpeVBenpRfsM5SBERELG3bl35Ei69hq7xzSQOyxbBsceG80rNGlzDG/n7GAuGfSo1Yg26b049oz2\nhb53UYNDaSgHISJSTuZOn87aCRN4/OEV4I5/8ikQDRBXXQWbNkXzCkOuXku/en2jeQUfVa55heLQ\nCEJEpIzFylx0zMtj3JYtnDJ7NinVl3JJo+mMXfcU1WrXZuHC3cdXZF6hOJSDEBEpY+7O5V2ms+S9\nfC5kJWMYzXm/XEKH9K6cdprtd3qoOHmF4ihpDkJTTCIixZSo1MXGjfD229FtM+P4o7PpVPt/2dx+\nGt1qjWXYtdn85S/7Dw5QMXmF4lCAEBEpptiS1Bf/MRE++giAzz+Hf/1r9zH1jv2RS5+6jvs+/pi2\nffuGZtqoOEo1xWRmrYAngZ8CDjzi7g+aWVPgWeAIYDVwobt/G5wzEhgC7ASucfd5QXtn4AmgDvCy\nu1+b4PM0xSQiSTMlM5Op48bR8dtvyf/6Oh7hQk46tAcX3T6aQVdeud/zymvqqKhKOsVU2gDRAmjh\n7h+YWX3gXaAfcDnwjbvfY2Y3AU3c/WYzSwOeBroCKcCrQDt3dzNbAlzl7kvM7GXgQXefs9fnKUCI\nSHKsX89ZaV9Qd9Micm0ri/Jv5cT6D/KzPqczdHj7UKw62p+k5CDcfZ27fxBsfw+sIPrFfy4wKThs\nEtGgAdAXeMbd89x9NbAK6GZmhwEN3H1JcNyTceeIiJS7RHmF11+PTh0BcMgh/GH4NgZOaEHXunfT\nvdlDnGB/Cs2S1PJQZjkIM2sDHA8sBpq7e26wKxdoHmy3BL6KO+0rogFl7/acoF1EpELMnT6dD57f\nyPRrx8C6dQCsXAlr1wYHVKtGz7+cyoYNq+k9cSJnjhhBn4kTK2VuoajK5D6IYHppOnCtu28xi79j\n0N3MNC8kIqE0ZcIEpt59Nx03beLQvKcY9/BKnnipBwNuvJFhCfIKV4wcCUDtQyA9/fyK7m6FKnWA\nMLOaRIPDZHefETTnmlkLd18XTB99HbTnAK3iTj+c6MghJ9iOb89J9HmjR48u2E5PTyc9bOvCRKTS\nWDl+DpNuSOOnDf/Iwh0/sIhz6FFvFY07zyYlNfmlLkoqKyuLrKysUr9PaZPURjTHsMHdr4trvydo\nu9vMbgYa75WkPoHdSeq2wShjMXANsAR4CSWpRaSU3J3fDXyUfzx1BWZGXh488AD84Q/R6qjbV2Wz\n8tNqrN36JnOHDOGt2mPptn0kfSZOpNf5VWd0kKxifT2AQcCHZvZ+0DYSuAuYZmZDCZa5Arj7cjOb\nBiwH8oERcd/4I4guc61LdJnrHsFBRKS45j7/PEuey2N29Ss5d/Ij1KgRfcTCtm1Qrx7UbtuKjm1h\nydhwlrpINpXaEJEqZ/Ko0Tw74SE6fvcdM/PepVvDa8hNyWXAtdeG+n6F8qJSGyJyUEi0HDXe3WlP\nsPhvrWjU5AEW1rqTjzmWTziXxsc8T0rqsELfuyoGh9JQgBCRSqXgyWsvvADABx/AY4/t3n/55J78\nbf3lDL6zFl2r/SV0JbQrEwUIEakUpmRmcnaHDiy44c+k7fgpr191FWd36MC/X3qa5s13H/fTzq2o\nU69aQQntg+F+hfKi50GISOjl5W4kZdXhNFr/Pyxcv4XFDKXH5o20+eWvOP6k9mRk7HvOwXS/QnnR\nCEJEks7duefmm1mwYN+8gn/+Bce1XM/Ry2cwePBnnNjgjui0UfXbizRtpLxCySlAiEjS7fFoTuDe\ne+Gtt6L77Mg2LMluScuXHiW7WTNNG1UgLXMVkaSZkpnJ1AceIOWbZlyx/geuazSYRimP8otf3cGV\n1/YjZT8V2arqctTykqwb5UREii8/n29nv8HAxYvJ+eIIJm2/lPXU5D/fnc/AM8/E6rYnEmG/AULB\noWJoiklEykVheYXFV09hwGV1sGOOoeP9felVfxht0m6jW62xDO6/gjFjCn80p1QMBQgRKRfxeYWt\nW6FXL8jPj+7r+vdLeHnTiXD99WRv2kTviRO5b9mySvtozqpKOQgRKVNTMjOZes89+Lpf8eTWKfRr\nOoZGLf5B57Pv4LY7+1G9+v7PVW6hfKjUhohUiP2Vuti0JMKGm+5h4P3388tvjmHF9nO4nT/zn41X\n07jDc1CnL2+8Ufh7KziEi5LUIlIsBaUuzn9hj5LY9170LselNmfAE0/QMTubvkPPx1q1otuqHxnc\nP5Ve56clsddSEppiEpEimZKZydQHH6RjXh41Ihez4rA1LPnxVP469vt9KqQ+OnYsrVNTOfO88xjc\nfwXpv5jFb2++OUk9Fy1zFZHys349dd9pTc7qR/k2fyGLGEm37x+ia0YXUlJT9zk8VuYC4Lcj0khP\n1+ihMlIOQkQS5hV88xYW3fISfmYvaNeOPt9O5aQTc+laeyzdmz1Ed0by20uWqdRFFaYAISIFeYW5\n019g166g8euvuWviT9nYfzjk5FDvuUl06rlSpS4OIspBiBzEpmRmMnX8eDru3EmNyMU8V/8EGjR9\ngatv6XpQPnmtqlIOQkSKbscOePVVPp9wCN9+egcL637OIv5IN5vAUSdcT0pq+0JPV3A4OGiKSaSK\nSpRX+Hru+yw852/QsiXccQf9+uxg4J92FTx5rTs3c8kAPXlNohQgRKqoWF5hxlOzCtpyl64ja/uJ\n8O67sGgRx911MTXq6MlrkphyECJVzJTMTKbedx8dgfzIEB6pMZCT2p3FRddepbzCQUo5CJGD3bJl\nMHUqT9zXmwZ1b2Hhzk0s4jpOqjueJsc8q7yCFJummEQqkX1KaLsz77KnWdHuXDjrLNi+nTv+cQiX\nPtiQrj6K7s0eoiu3FenRnCJ70whCpBKZO3062Q89xuLlQ8jISAUz1uc3oeEfb4crjoVq1egGPDp2\nBr0nTqThh+dx0nEtlFeQElEOQqQSmHLvvUx99FE6Am0jJzGqzgiO+9nlDLjmGuUV5IBU7lukktrv\nk9c2bIBHHmFl98uYesux/DJtBAtzh7OaNmT/eAKNOzxHZM0wsrL2/94KDlIammISSbKCJ699fgkZ\nGWnkfbSS23/9HqO/HkG1Pr1o+8eLGN8+g1Urvmfdv4ewJa0V3VbVVgltKXcKECJJEl8+e+iWVgx7\n9XXO7nABA4YMoelJGWy79yt+0rw+NYC2wIJZwf0KQQnt7MisA32ESKmUOgdhZv8EfgV87e7HBm1N\ngWeBI4DVwIXu/m2wbyQwBNgJXOPu84L2zsATQB3gZXe/NsFnKQchVcP27ex8eQ5/e78Bsx9YyrbN\nv+Q9ujDwguUc1b49GRlW6PSQcgtSHMnMQUwEeu/VdjMw391TgdeC15hZGtAfSAvOmWBmsULine truncated
|
||||
"text": [
|
||||
"<matplotlib.figure.Figure at 0x7f5b379ae790>"
|
||||
]
|
||||
}
|
||||
],
|
||||
"prompt_number": 5
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"plot(solE.imag,M.vectorNx,'r*--',-anaE.imag,M.vectorNx,'b+:')"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"metadata": {},
|
||||
"output_type": "pyout",
|
||||
"prompt_number": 6,
|
||||
"text": [
|
||||
"[<matplotlib.lines.Line2D at 0x7f5b36fded50>,\n",
|
||||
" <matplotlib.lines.Line2D at 0x7f5b36fdefd0>]"
|
||||
]
|
||||
},
|
||||
{
|
||||
"metadata": {},
|
||||
"output_type": "display_data",
|
||||
"png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEACAYAAACkvpHUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4FFXWx/HvYQlbIIIg+zYOKHEQFRE3JHkVAUVUXMAB\n3FAcUNFhFkHUJDMq4z6io4OiCEZFFAdEFFCHFnVEXBBBEBsQDYGEPSxhz33/qErShBBI0iGd5vd5\nnn5SXV3VfVJon9Q9t06Zcw4REZEjVam8AxARkYpFiUNERIpFiUNERIpFiUNERIpFiUNERIpFiUNE\nRIqlVInDzKqb2Zdm9p2ZLTGz0f76emb2oZn9ZGazzey4kH1GmlnQzH40s4tD1nc0s0X+a0+XJi4R\nESk7pUoczrldQKJz7jTgVCDRzM4HRgAfOufaAh/7zzGzeKAvEA/0AJ4zM/Pf7nlgkHOuDdDGzHqU\nJjYRESkbpR6qcs5l+4sxQGVgM9AbmOCvnwBc4S9fDrzhnNvrnFsFLAc6m1ljoLZzbr6/3cSQfURE\nJIKUOnGYWSUz+w7IBOY4534AGjrnMv1NMoGG/nITYHXI7quBpoWsT/fXi4hIhKlS2jdwzuUAp5lZ\nHDDLzBILvO7MTH1NRESiRKkTRy7nXJaZzQA6Aplm1sg5l+EPQ63zN0sHmofs1gzvTCPdXw5dn17w\nM5SARERKxjlnh9/qyJR2VlX93BlTZlYD6AYsAN4FbvA3uwGY6i+/C/Qzsxgzaw20AeY75zKArWbW\n2S+WDwzZ5wDOuYh/JCUllXsMilNxVtQYFWf4H+FW2jOOxsAEM6uEl4Redc59bGYLgMlmNghYBVwL\n4JxbYmaTgSXAPmCoy/+thgKvADWA951zM0sZm4iIlIFSJQ7n3CLgjELWbwIuOsQ+DwMPF7L+G6B9\naeIREZGypyvHy0BCQkJ5h3BEFGd4VYQ4K0KMoDgjnZXF+FdZMTNXkeIVEYkEZoaLlOK4iIgce5Q4\nRESkWJQ4RESkWJQ4RESkWJQ4RESkWJQ4RESkWJQ4RESkWJQ4RESkWJQ4RESkWJQ4wsw5x22/f6FM\nOlKKiEQCJY4wmzVlCgunbGT2O++UdygiImVCvarCJHXsWCaNGUOHvXupEvw9+9q8zsKqVek3bBgD\nbrutvMMTkWNYuHtVhe0OgMe6pm0HExd/Pp/Mns3n/JHzMuOofvqVxDVpWd6hiYiElYaqwiQx0RjY\ndymdXBJn1/8XZ7okatXcx6mn5if5nH05B+wzZ47j0REjVA8RkQpFiSOM0oJBeowfz8VDh9Jz/Hgu\nu+BtWvonHFu2wMmNNrPvzLPhiScgLY2Xnl/K2ueeUz1ERCoU1TjKQCAAhd3fZWPmPo5fFCA1KYlx\nX2wincf4yV3GfSeeyMJq1VQPEZEyEe4ahxLHURYIeENUC+b9xPTZJ5FEMh/ViiXhD+cy6u/nUKNG\n2P5tRUQA3cipwktIgJQUY+iti+gcM5qt8W9xZqW/sSajNuPHh/y7OgdZWYCXbJxTPUREIoMSRzlJ\nCwb57eWX88TixfQcP55zfzeDIUPyXx/UdzvLmv4f9O5N4JlFzEpNVT1ERCKChqrK0aFqIQDz5sHp\nJ27lrb/eRcprbWi/9zsmM4Wkpo1ZGBeneoiIHDENVUWRQyUNgF27YPS/6hBs8TLL997L/6rcx9+5\nn5+3dWTIAylcd8vgg/YJBLyfGtYSkbKkxBGhEhIgOdmrh/S/ZgnX1TifbfFvcYKbw7ffN2bgwAJ/\nPASDBD7eD3htTzSsJSJlRYmjAsjauJEe48fn1UMaxX7K88/nv/7OOxC49TW+f+QRetWty6d3382T\n27Yxd+RIep1yCqljx5Zf8CISdVTjqACKqoUEAvDKK1CtGrzwAvy+/TssW7yGR9w7fBT7FV1feonu\n11yDmab5ihyrVOM4BhVVC0lI8BLH2LGQlAQDH8hhU6UTmdg8hp379mGVKx+UNHJrIaB6iIgUnxJH\nlEkLBnn+zWxe+eUDeqam8tOiX+jWDfbvz98m8J/NsHAhOKd6iIgUW6mGqsysOTAROAFwwAvOuTFm\nVg94E2gJrAKudc5t8fcZCdwM7AeGOedm++s7Aq8A1YH3nXN3FfJ5x+RQ1ZEqbEhr3z4vR3Ts6D1P\nT4eRA9K4+PvTmLR9Ox1iY3lw0ybua9NGbeBFolREtRwxs0ZAI+fcd2YWC3wDXAHcBGxwzj1qZvcA\ndZ1zI8wsHngd6AQ0BT4C2jjnnJnNB+5wzs03s/eBMc65mQU+T4mjhAIB77FmDbz4IjzwgGP5p5+y\n74vHeXPXdEZWrUrXlBS6jxiheohIlImoGodzLsM5952/vB1YipcQegMT/M0m4CUTgMuBN5xze51z\nq4DlQGczawzUds7N97ebGLKPhEHu9N4XXvBqISkpxsCh61iRk8hFjZLYWbUq1rDhAUkjtBYCqoeI\niCdsNQ4zawWcDnwJNHTOZfovZQIN/eUmwOqQ3VbjJZqC69P99VKG0oJB/p7agndXJNFz4kTS1q3j\n4Yfhl1+81wMBvOLI2LGQkaF6iIgAYUoc/jDVFOAu59y20Nf8sSX9iRpBcusgt44cSc9rrqJmTaP7\nVVdxy4gRNG8O9erlb7to3g5efXk8vZo25dPrr/euD7nnHl0fInIMK/WtY82sKl7SeNU5N9VfnWlm\njZxzGf4w1Dp/fTrQPGT3ZnhnGun+cuj69MI+Lzk5OW85ISGBhKLmqkqhijpkzZt795kCSEmB8ePr\ncMMNX9Dloo/Y8sLvsZ07yVm5kjsuu4zugw9ueyIi5S8QCBAoONYcTs65Ej8Aw6tHPFVg/aPAPf7y\nCOAf/nI88B0QA7QGVpBfoP8S6Oy/5/tAj0I+z8nRk5SUv/zBW2+5PjUuc2fUe8XdVbu2m/nss3mv\nzZlz4H45OTnukXvucTk5OUclThEpmv/dWarv+9BHac84zgMGAN+b2QJ/3UjgH8BkMxuEPx3X/9Zf\nYmaTgSXAPmCo/0sBDMWbjlsDbzruATOqpHylBYNcP+5mTu54OasWx5IWDPL55xAXV2Aa8IsvMuvH\nH1n74ovM7tSJ7lddVY5Ri0hZUMsROaSiWp0ATJoE9evDZ595M7ZefmYc7zySTIfNm3kwO5v7atdm\n4XHH0e+eexhw++1HKWoRKSiipuNKdCsqaQQC8OOPXtJISfESxwOPDKLt5e+Qc/zxGJBTpQp3xMbS\nPyUFsrMPmt4LmuIrUhGVujgux6aEhAMTS3IyjBhhzJn+K7Nf3cKwk08jsOJ2uiTXxS5KhJo1Cz2D\nyZviq2EtkQpDZxwSNtWrw+rlQXqMH8/DX31Lrzvas2ZlEOrVY/Nm2LQpf9vUsWPp1aoVn95yi1rA\ni1QwqnFIqR2u7XsgAKtWwYQJ3lXrAF27Onb9+G/mjhzJ6KwsRlaqRNdLLqH7/fdjnTqB2p6IhE24\naxwaqpJSO1zb99zXW7XyhrQAkpONXVvOYHdODsPj48n55ResWjVswAAC1z5HwoMXHfRezjkeGzmS\nv4werX5aIuVIQ1VSLu68E06o9UXenQ1r3PoxP5x8ASxbRqDS/xW6j1qeiEQGJQ45akLPTI4/HoY/\ndDfdr7oKMyOmTmeuv3uYN0RVqRLbQhrXpI4dS6/4eD4dONCrh/z5z6qHiJQjDVXJUVPUkFZiIjz7\nrLeckgJjxnhnJYmJ0H/wYI6vW5e5d96J7dpFzq+/ckerVnTfuhV+/RVatDgq8YuIR4lDIkLB6b2j\nRkHVqt7yV18Z/5lxErV27vTqIWlp2IABfBKAhMX3e1X3AlQPESk7GqqSiJSbNAAaN4Y6MYvy6iHt\nH3qLH/YfT6DTXwpNGgCz3n5b9RCRMqLEIRGn4JBWs2bw+IsD8uoh23O607T9sEL3TR07ll6nnJJ/\nfcjtt9Pr5JNVDxEJIyUOiTiHa3WyeTMsXZrf6uTEE+HVV73X+w8ezO3JyeTUqeO1Pdm6lQt+bkn/\n996D116DXbsOek+1PREpHtU4pEIprNXJjTfm18f37zemvt+GGlu25tVD0ro/ivX5Ht56C3r39i5x\nD6G2JyLFozMOqfBatYJK/n/J27dDxuqt9HzlZZ5YvJgLx05g85YdMHAgTJ0KtWvn7Zc3rDVypNqe\niBSDWo5IhXUkrU6CQXj99fxWJ6FnLM45Zr79NnNvv53R69czsmZNut59N92TkrCYmLIOX+SoUVt1\nEd/hWp0kJ3tljfvv95aTk+Gnn2Cmf4swM8PM2LVrF8PbtCG49zxs8mSsaVO47Tb4/vuD3lf1EBEl\nDjkGVAr5r/yss6BNm/znn3y8g4R/T+SJZcuofsU/SRs0CL76yqu4b9hw0Hup7YmIhqrkGFDUkNb1\n18NDD0Hz5t4ZSVJS4Y15U8eOZdKYMXTYu5cHg0Hua9OGhVWr0m/YMAbcdlsZRi9SeuEeqlLikGNa\nbi0EvOm9DRt6o1SJiQcmm7x6yJ/+xOi0NEZWrkzX3r3p/sADWIcOagMvEU1t1UXCqOD03ltu8S44\nBFi2DBYuhGuvDamHbNnCNS1voOn6t7CcHOyKK7zpvYMGwV/+csB7q+2JRCvVOERC5CYNgL17Yf/+\n/OdLF66m27jxxN8wnp4TJ5J29tnw88/wyivQqNFB76V6iEQrDVWJ+IqqhYDXeLFlS1izJv+GVIU5\nqB7SujULa9RQPUTKjabjipSRw7U6qVoV0tPzW5106gSTJh28bV7bk127vLYnq1Zxx/799N+9GzIy\nyiZ4kaNIiUPkCOReF5KS4s28Sk72ZmP16eO97px3fUhOTn49JLjhdwyPj2dnbCzWpw/29dfQrh38\n3/9BWtoB76/rQ6QiUeIQKaGLL4bcC8w3bvQuNsytgf/6U5DqvR7nicWL6Tl+PGl16sDEibLine truncated
|
||||
"text": [
|
||||
"<matplotlib.figure.Figure at 0x7f5b379c6810>"
|
||||
]
|
||||
}
|
||||
],
|
||||
"prompt_number": 6
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"semilogx(abs(solE),M.vectorNx,'r*--',abs(anaE),M.vectorNx,'b+:')"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"metadata": {},
|
||||
"output_type": "pyout",
|
||||
"prompt_number": 7,
|
||||
"text": [
|
||||
"[<matplotlib.lines.Line2D at 0x7f5b36fcda50>,\n",
|
||||
" <matplotlib.lines.Line2D at 0x7f5b36f1fd50>]"
|
||||
]
|
||||
},
|
||||
{
|
||||
"metadata": {},
|
||||
"output_type": "display_data",
|
||||
"png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEFCAYAAAD9mKAdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VOX1x/HPYVOUTaDKqlSFSkBRwaJiNbFVQlFRUYMG\n1IoFwbp2M1YlWC1qKyr6E2LVgERBNnGpsrgEXCqodYGCErRg2HfBBQLk/P6YmzDEMITMJDNJvu/X\nKy/uPHPnzkkyzMl9znPPmLsjIiKyL7XiHYCIiCQ2JQoREYlIiUJERCJSohARkYiUKEREJCIlChER\niSiqRGFmB5vZPDP7xMwWmdmIYLypmc02syVmNsvMmoQ9JsPM8szsczM7N2y8q5ktCO57JJq4REQk\ndqJKFO6+HUhx9xOBE4AUMzsDuA2Y7e4dgDeC25hZEpAGJAGpwONmZsHhRgMD3b090N7MUqOJTURE\nYiPqqSd3/z7YrAfUBjYDFwDjgvFxwIXBdh9ggrvvdPdlwFKgu5m1BBq6+/xgv2fCHiMiInEUdaIw\ns1pm9gmwFnjL3f8LHOHua4Nd1gJHBNutgBVhD18BtC5lfGUwLiIicVYn2gO4eyFwopk1BmaaWUqJ\n+93M1CdERKSKijpRFHH3b8zsX0BXYK2ZtXD3NcG00rpgt5VA27CHtSF0JrEy2A4fX1nyOZRwRETK\nx91t/3uVLtpVT82LVjSZWX3gHOBj4CXgqmC3q4DpwfZLQD8zq2dmPwXaA/PdfQ2w1cy6B8XtAWGP\n2Yu7l+lr2LBhZd430v6ljYePxep5YhlTrJ7nQPfdXwyKqWy/v4qKKdL+ep2X//dXFWKKVrRnFC2B\ncWZWi1DSGe/ub5jZx8AkMxsILAMuA3D3RWY2CVgE7AKG+p7vYigwFqgPvOruM6IJLDk5OSb7lzZ+\noMcu72OrUkyl3aeYyn5feePS67zsatJrKpqYSnUgfyXE+ysUbmIZNmxYvEMoVSLGpZjKRjGVXSLG\nlYgxBe+d5X7v1ZXZUYp55o6RRIxLMZWNYiq7RIwrEWOKlnkM5q8qi5l5VYpXRCQRmBker2K2iIhU\nf0oUIiISkRKFiIhEpEQhIiIRKVGIiEhEShQiIhKREoWIiESkRCEiIhEpUYiISERKFCIi1Vgsulko\nUYiIVGMzp06N+hgx++AiERFJHDlZWUwcNYouO3dGfSydUYiIVEOtOwyicdIk5qwdEvWxlChERKqh\nlBSj/2WLaPbtK1EfS4lCRKSaenP6dD4svCHq4yhRiIhUMzlZWZzXqRN1582nNS2jPp6K2SIi1Uyo\nPnEGc197jQ/oHvXxlChERKqZlBRjx8bFzHjlTlbXOoKvCqM7nqaeRESqofy8PDa2vI5vCztGfSwl\nChGRaiYnK4sXc3LYtm033VkZ9fGUKEREqpmiayg2FPyUl+kT9fFUoxARqWZSUoztGxYxZvqrdG9e\nj3kbojuezihERKqh5576lI8Kh3BsSkrUx1KiEBGpRoquoTh48SecWfgDvPe3qI+pRCEiUo0U1Sfe\nXnU5E0hn2bauUR9TiUJEpBpJSTE6H/k0DXblcXKDB+nmw6I+ZlSJwszamtlbZvZfM1toZjcG403N\nbLaZLTGzWWbWJOwxGWaWZ2afm9m5YeNdzWxBcN8j0cQlIlITFU07vffUdL7jXA6pW5ePGzeO+rjR\nnlHsBG5x907AqcD1ZtYRuA2Y7e4dgDeC25hZEpAGJAGpwONmZsGxRgMD3b090N7MUqOMTUSkRima\ndlpfcBOfczpbdrSk7Wkzoz5uVInC3de4+yfB9rfAYqA1cAEwLthtHHBhsN0HmODuO919GbAU6G5m\nLYGG7j4/2O+ZsMeIiEgZpKQYA9IWc2rhnzi61njOrj2QAWmLoz5uzGoUZtYOOAmYBxzh7muDu9YC\nRwTbrYAVYQ9bQSixlBxfGYyLiMgByM/L44jWrdnhR9Fy8GDy8/KiPmZMLrgzswbAVOAmd9+2ZzYJ\n3N3NLPpP9xYRkYiKWncs+9/dnOn5bHvxRd6pWzfq40adKMysLqEkMd7dpwfDa82shbuvCaaV1gXj\nK4G2YQ9vQ+hMYmWwHT5eaoOSzMzM4u3k5GSSk5Oj/RZERKqFovbi25a+wISCXbRd/lOaHHsssCiq\n45p7+f/YDwrR44CN7n5L2PgDwdj9ZnYb0MTdbwuK2c8BPyc0tfQ6cGxw1jEPuBGYD/wLGOXuM0o8\nn0cTr4hIdffa5MmMueIJ1ja5kFN3ZNArO5vUSy7B3W3/jy5dtGcUPYD+wGdm9nEwlgHcB0wys4HA\nMuAyAHdfZGaTCKW3XcDQsHf+ocBYoD7waskkISIi+/fm9Oks9bM4KSWFXmnZMalRRHVGUdl0RiEi\nUrqcrCwmjhrFihWZnLp1A9+2fpctjT+m3403MuC666I6o9CV2SIi1UDRNRTffFefLIawbFtXmnSa\nTOsOg6I+ttqMi4hUAykpxvtTR/PF7uY0bbiEbp5Jr7QjSUlJivrYOqMQEaniilp3fJSdzXHAIXXq\n8HHjxryYkxOT4+uMQkSkiitaFvvm4u2soSudd0ymy69mcumQ6D8vG5QoRESqvKJPtNs8ZRzrm/Sm\nx+4/0ystOybTTqCpJxGRauHN8eNZTjd+9qsz6ZUdm2WxRZQoRESqsKL6RO15/+ESgH+P4NG77uLg\nww6L2XNo6klEpAorqk+Mf+F7VnEKPbY9RLvTb6d1h9jUJ0CJQkSkSiuqT2x74Z+0bnYx3QqGxWxZ\nbBFNPYmIVHFvTp/Ol5xJh7NjX58AJQoRkSqrqD6x6/2lXLTb4b2/xbw+AZp6EhGpsorqE9Ne2sLX\n9KiQ+gQoUYiIVFlF9Ylvpz9By+YX0m1H7OsToKknEZEq7c2cHJYWnsmxKSkVUp8AJQoRkSqpqD6x\n7e119Co8uMLqE6CpJxGRKqm4PrHoe9ZX0PUTRZQoRESqoJQU4/0XsviCprRtMJduPrxC6hOgqScR\nkSqnuK342LGhtuJ168a0rXhJOqMQEaliiqadPvliIYtIi3lb8ZJ0RiEiUsWkpBj9L1tErcL1HN/o\nWX5ZeyAD0haTklLuj8WOSIlCRKQKenP6dAprbeGEnidV2LLYIkoUIiJVSFF9ot4HH3DJ7l0Vuiy2\niGoUIiJVSHFb8Ze2ks9pFbostogShYhIFZKSYvywfhH/nZZPy+afVljbjnCaehIRqWLeemEaO+17\n2qdUTFvxkpQoRESqiOL6xAsvVFp9AjT1JCJSZbRq/1sWrf8Z7+1oxWY6VEp9ApQoRESqjIJN0zh4\n/ef0rDWVZU2vrZT6BChRiIgkvJysLCaOGkWXnTs5hiuo23QFG4BDLr+8wusTEIMahZk9bWZrzWxB\n2FhTM5ttZkvMbJaZNQm7L8PM8szsczM7N2y8q5ktCO57JNq4RESqi9CS2Em8tXIor5DJVwW3065T\nb87pN4Zrb7utwp8/FsXsbCC1xNhtwGx37wC8EdzGzJKANCApeMzjZlZ0zfloYKC7twfam1nJY4qI\n1EhFLTu2fH8Ibeq+TTcfxh9u+KjCWnaUFHWicPe3gc0lhi8AxgXb44ALg+0+wAR33+nuy4ClQHcz\nawk0dPf5wX7PhD1GRKTGe3P6dKz2Os44v1mlLIkNV1HLY49w97XB9lrgiGC7FbAibL8VQOtSxlcG\n4yIiNVrJlh2151XOkthwFV7Mdnc3M4/V8TIzM4u3k5OTSU5OjtWhRUQSTlHLjhdf/Zz/0rdMS2Jz\nc3PJzc2NWQzmHv17uJm1A1529+OD258Dye6+JphWesvdjzOz2wDc/b5gvxnAMGB5sE/HYPxy4Cx3\nv67E83gs4hURqSrcnfN6DGPJByfRtPFaTiv4E72ys+nZt2+Zj2FmuHu5CxoVdUbxEnAVcH/w7/Sw\n8efMbCShqaX2wPzgrGOrmXUH5gMDgFEVFJuISJUxc+pUls6ry8Bb3+OHQx7g9BN+Uqn1CYjBGYWZ\nTQDOApoTqkfcBbwITAKOBJYBl7n7lmD/24FrgF3ATe4+MxjvCowF6gOvuvuNpTyXzihEpEYovnZi\n+3bqfDWAXe2f482C07k+4zT6Dx58QMeK9owiJlNPlUWJQkRqirfecp58fBEfvPAlebsvoEejh2jX\nsycDh3Q84GWxiTr1JCIiUUhJMbbnf8SKKW9Tv9FWunnltOsojbrHiogkqDcffpgNtOH4nidX+rUT\n4ZQoREQSTE5WFr2PS4JPl3IJXmntxPdFU08iIgmmVfvf8tm6E/m3H1ap7cT3RWcUIiIJpmDTNBpt\nnU3aQd05tfn/0c2HMSBtcaX1dipJZxQiIgkivJ340buv4LufHMyGXZXXTnxflChERBJE6w6DOLTD\nmTz32jqWcRY9vm9Mu5N6c06/oXE7mwAlChGRhODufDAzgw7Hns86XqFF84WhT7C74UhSUtrFNTYl\nChGRBDBz6lRWP/YY89p1Y9izp/LeZxdz+gkt4jrlVESJQkQkjorqEifsKODO7+rSZ/lyHr3rabr0\nrHtAjf8qklp4iIjEkbvzwPA5TPzbl2zZeTbL+Cnply7imI6hVh2x+CQFtfAQEani3p0+nbN2PoIf\n25F5Xw9gQFoHevat/FYd+6LrKERE4mjmhAl88WkTWvTuzcNL/suxffokRF0inBKFiEgcFH3E6WPX\nzeMYTmHbkiWc37kzR7VZyrW33Rbv8PaiqScRkUrm7kx8viWNkyax9Ou5zKc3W9cuoV3Pnvzq/Pi0\n6YhEiUJEpJLNnDqVDe/m0eOGuRxuT1CreWHomok4tRHfHyUKEZFKUrQUttMO2FkwjvefvpQtjRvT\nokt3ev0mfm3E90eJQkSkkrTuMIjGSWfw7sxZ/Idu1N91I+1O60nX0zrSs2+3eIe3T0oUIiLine truncated
|
||||
"text": [
|
||||
"<matplotlib.figure.Figure at 0x7f5b370013d0>"
|
||||
]
|
||||
}
|
||||
],
|
||||
"prompt_number": 7
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"anaE[1].conj()"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"metadata": {},
|
||||
"output_type": "pyout",
|
||||
"prompt_number": 8,
|
||||
"text": [
|
||||
"(-1.2220068301411678e-08-1.5509245597740526e-09j)"
|
||||
]
|
||||
}
|
||||
],
|
||||
"prompt_number": 8
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": 8
|
||||
}
|
||||
],
|
||||
"metadata": {}
|
||||
}
|
||||
]
|
||||
}
|
||||
+105
-51
@@ -1,6 +1,7 @@
|
||||
{
|
||||
"metadata": {
|
||||
"name": ""
|
||||
"name": "",
|
||||
"signature": "sha256:5ff50687aa6b5eff8d9f077fee29dfb14e6eaa12563cbe8511b1dddc2ba76e23"
|
||||
},
|
||||
"nbformat": 3,
|
||||
"nbformat_minor": 0,
|
||||
@@ -22,7 +23,7 @@
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": 35
|
||||
"prompt_number": 1
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
@@ -33,7 +34,7 @@
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": 18
|
||||
"prompt_number": 2
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
@@ -61,12 +62,13 @@
|
||||
]
|
||||
}
|
||||
],
|
||||
"prompt_number": 25
|
||||
"prompt_number": 3
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"# Setup the model\n",
|
||||
"conds = [1e-2,1]\n",
|
||||
"sig = Utils.ModelBuilder.defineBlock(M.gridCC,[0,0,-300],[500,500,-100],conds)\n",
|
||||
"sig[M.gridCC[:,2]>0] = 1e-8\n",
|
||||
@@ -80,26 +82,28 @@
|
||||
{
|
||||
"metadata": {},
|
||||
"output_type": "pyout",
|
||||
"prompt_number": 34,
|
||||
"prompt_number": 4,
|
||||
"text": [
|
||||
"<matplotlib.colorbar.Colorbar instance at 0x10b57b7e8>"
|
||||
"<matplotlib.colorbar.Colorbar instance at 0x7f55faf27290>"
|
||||
]
|
||||
},
|
||||
{
|
||||
"metadata": {},
|
||||
"output_type": "display_data",
|
||||
"png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEKCAYAAAAPVd6lAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XtcVHX++PHXACNpiiYpCKNiiAmIiBe02lZbQw1X80Kk\ntuK90lyvpXYV259CtV00pXK/WGYmWrsJ3tg0b5WKF9RUTKcE5V6KIISCjOf3BzJJDDDAXNn38/E4\nj5hzPucz7zdDb898zjmfo1IURUEIIYTdcrB2AEIIIRpGCrkQQtg5KeRCCGHnpJALIYSdk0IuhBB2\nTgq5EELYOSnkQghh56SQC5v3+eef07t3b1q0aIGHhwehoaF8//33AJw/f54nnniCNm3a0KpVKwID\nA3n33Xe5deuWlaMWwnKkkAub9s477zB37lxeeeUVfvnlF9LT03nuuedISEjg559/pm/fvnTs2JHT\np0+Tn5/PF198wbFjxygsLLR26EJYjEru7BS2qqCgAI1GwyeffMLo0aOrbP/b3/5GQUEBW7ZssUJ0\nQtgOOSIXNuvgwYPcuHGDkSNHGtz+zTffEBYWZuGohLA9UsiFzbpy5Qr33nsvDg6G/0yvXLlCu3bt\nLByVELZHCrmwWa6urly+fLnaE5eurq5kZWVZOCohbI8UcmGzHnjgAZydnfnqq68Mbn/00Uf597//\nbeGohLA9UsiFzWrZsiWvv/46zz33HPHx8RQXF3Pz5k127NjBwoULWbJkCQcOHGDBggXk5uYC8NNP\nPzF+/HgKCgqsHL0QliOFXNi0efPm8c477/D//t//o23btnTo0IGYmBhGjhzJfffdx8GDB0lLS8Pf\n359WrVoRFhZGnz59aNGihbVDF8Ji5PJDIYQwscTERObMmYNOp2Pq1KksXLiwSptZs2axY8cOmjVr\nxieffEJQUFC930+OyIUQwoR0Oh0zZ84kMTGRlJQUNmzYwNmzZyu12b59Oz/99BNarZbVq1czffr0\nBr2nFHIhhDChw4cP07lzZ7y8vFCr1YwZM4b4+PhKbRISEpgwYQIAffv2JT8/X3+epz6kkAshhAll\nZmbSvn17/WuNRkNmZmatbTIyMur9nlLIhRDChFQqlVHt/nh60tj9DHGq954W1KNHD06ePGntMIQQ\ndiAwMJATJ040qI+mKhU3jGzbvHnzSpO0eXp6kp6ern+dnp6ORqOptM8f22RkZODp6VnveO2ikJcX\n8UhrhwHsAR4xSU+KshgAlWqJ1WJomKpx1D8n08VgDhV5VScyMpLIyEizx2HrMdhKHA05sq1wA/h/\nRrZ9paio0uvevXuj1WpJS0vDw8ODjRs3smHDhkpthg8fzsqVKxkzZgyHDh2iVatWuLm51Tteuyjk\nQghhaep67ufk5MTKlSsZPHgwOp2OKVOm4Ovry0cffQTAM888Q2hoKNu3b6dz587cfffdfPzxxw2K\nVQq5EEIY0JDi+Nhjj/HYY49VWvfMM89Uer1y5coGvENlUsjrxMvaAWAbMYBtxOFl7QAAGDBggLVD\nsIkYwHbiMIWm1g6gDqSQ10knaweAbcQAthGHLcRgG8XLFmIA24nDFOo7tGINUsiFEMIAeyqO9hSr\nEEJYjByRCyGEnbOn4mhPsQohhMXIEbkQQtg5KeRCCGHn5PJDIYSwc/ZUHO0pViGEsBgZWhFCCDtn\nT8XRnmIVQgiLkSNyIYSwc/ZUHO0pViGEsBh7OiKXR70JIYQBTY1c6iIvL4+QkBC6dOnCoEGDyM/P\nN9guKioKf39/AgICGDduHCUlJTX2K4VcCCEMUBu51EV0dDQhISGcP3+egQMHEh0dXaVNWloa//rX\nv0hOTubUqVPodDri4uJq7FcKuRBCGOBk5FIXCQkJTJgwAYAJEyawefPmKm1cXFxQq9UUFxdTVlZG\ncXFxrc/zrLGQT548GTc3NwICAvTrXnjhBXx9fQkMDGTUqFEUFBTot0VFReHj40PXrl35+uuv9euP\nHTtGQEAAPj4+zJ49W7++pKSEJ598Eh8fH/r168fFixdr+TUIIYRlqJ2MW+oiNzdX/2xONzc3cnNz\nq7Rp3bo18+fPp0OHDnh4eNCqVSseffTRGvutsZBPmjSJxMTESusGDRrEmTNnOHnyJF26dCEqKgqA\nlJQUNm7cSEpKComJicyYMQNFUQCYPn06sbGxaLVatFqtvs/Y2FhcXV3RarXMnTuXhQsXGvnrEEII\n83JyMrwcVMFbyu/LH4WEhBAQEFBlSUhIqNROpVIZfFD0zz//zHvvvUdaWhpZWVkUFRWxfv36mmOt\naePDDz9MWlpalSAr9O3bl3//+98AxMfHM3bsWNRqNV5eXnTu3JmkpCQ6duxIYWEhwcHBAERERLB5\n82aGDBlCQkICS5aUP3F99OjRzJw5s8ZghRDCUtSOhtf/xRH+csfr6N8qb9+5c2e1fbq5uZGTk4O7\nuzvZ2dm0bdu2SpujR4/y4IMP4urqCsCoUaM4cOAATz31VLX9NmiMfM2aNYSGhgKQlZWFRqPRb9No\nNGRmZlZZ7+npSWZmJgCZmZm0b98eKH/ydMuWLcnLy2tISEIIYRLVHZH/camL4cOHs3btWgDWrl3L\niBEjqrTp2rUrhw4d4vr16yiKwq5du/Dz86ux33oX8qVLl9KkSRPGjRtX3y6EEMJmqZ2NW+pi0aJF\n7Ny5ky5durB7924WLVoElB8IDx06FIDAwEAiIiLo3bs33bt3B+Dpp5+usd963RD0ySefsH37dr75\n5hv9Ok9PT9LT0/WvMzIy0Gg0eHp6kpGRUWV9xT6XLl3Cw8ODsrIyCgoKaN26dTXvuueOn72wlQfv\nCiGsa+/evezdu9f0HZvhdsnWrVuza9euKus9PDzYtm2b/vWCBQtYsGCB0f3WOdTExETeeust9u3b\nx1133aVfP3z4cMaNG8e8efPIzMxEq9USHByMSqXCxcWFpKQkgoODWbduHbNmzdLvs3btWvr168eX\nX37JwIEDa3jnR+oaqhDif8CAAQMYMGCA/nXFebcGs6P73msMdezYsezbt4/Lly/Tvn17lixZQlRU\nFKWlpfqTng888AAxMTH4+fkRHh6On58fTk5OxMTE6M/IxsTEMHHiRK5fv05oaChDhgwBYMqUKYwf\nPx4fHx9cXV1rvehdCCEsxo4KuUqpuEbQhpX/gxBp7TBMSlEWA6BSmejowQY0xpzg97yEfVCpVDS0\nrKlUKhQfI9tqafD7NZQd/ZsjhBAWZEfV0Y5CFUIIC6rjFSnWJIVcCCEMsaPqaEehCiGEBdlRdbSj\nUIUQwoKquUXfFkkhF0IIQ+yoOtpRqEIIYUF2VB3tKFQhhLAgO6qOdhSqEEJYkFx+KIQQds6OqqM8\ns1MIIQxxNHKpgy+++AJ/f38cHR1JTk6utl1+fj5hYWH4+vri5+fHoUOHauxXCrkBOTnz6dmzHQD7\n9k1kzJhu+m1+fm3YtCmMc+dmUlb2KqtXD6uyf//+HdHpXquyTJrUw2I5/FFDcwJwcFCxcOFD/Pjj\nc1y//jI5OfNZtSrUIvFXp6F5ffzx4wY/q7KyV3F1bWqxPIQNMsPTlwMCAvjqq6/485//XGO72bNn\nExoaytmzZ/nhhx/w9fWtNVRxB2/ve2jWTM3x49mo1Q707u3Bd99d0m9v2tSJtLQC4uPPMW/eAzVO\nlhMU9BHZ2YX619eulZg19uqYKqdPPnmcvn01LFiwkxMncmjRwhkvr1aWSqMKU+Q1a9YOFiz4/dFc\nKpWKzZufpKiolCtXrlskD2GjzFAdu3btWmubgoICvv32W/2ThCqenlYTKeR/8NBDHUhKykRRoE8f\nT65cKSYj45p++7Fj2Rw7lg3AlClBNfZ1+XIxv/5abNZ4jWGKnAYM8GLMmG507/4hP/54Wb/+9Olf\nzBt8DUyRV2FhKYWFpfrXPj6t6dtXwxNPfGHe4IXts9INQampqbRp04ZJkyZx8uRJevXqxfLly2nW\nrFm1+0ghv+3q1YUoioKzsxMODiry8hagVjvi7OxIXt4CFAVcXd+sU5/ffTeJZs3U/PRTHh99dIx1\n634wU/SGmTKn0aN9uXDhKoMGebNly1iaNHHk4MF0nn9+Z6XiaQnm+KwqPPNMb3Jyiti8+UcTRy3s\nTjXVcW8W7M2ufreQkBBycnKqrF+2bBnDhhketrxTWVkZycnJrFy5kj59+jBnzhyio6N5/fXX6xrq\n/57u3T9ApVJx6NAUnn12GydO5BAXN5rPPz9NfHzd/qfOyipkxoxtHD2axa1bCqGhPqxePYzOnVuz\nePFe8yRggClz8va+hw4dWvK3vwUwZUoCpaU6li79C7t3R9Ct2weUlurMlEVVpszrTk2aODJhQiAf\nfXSMW7dsfpp+YW53GV494L7ypcKS45W379y5k4bQaDRoNBr69OkDQFhYGNHR0TXuI4X8tvT0awQE\ntEWtdmTLlnM0b96EHj3cGT48jsuX6zY8otXmodXm6V8fP56Do6MDzz//AEuW7LNYkTBlTg4OKpyd\nnYiI2KwfWnnyyS/Jzp5PaKiPRY9gTZnXncLC/LjnnrtYvfqYCaMVdsvMQyvVnYtyd3enffv2nD9/\nni5durBr1y78/f1r7EsKOXD69HQ6dGiJk5MDarUjBQWL9IXrwoXy54v6+q4iM7Owlp6ql5SUwd13\nN6FNm2bk5v5mqtCrZeqcsrOLUBSl0vj45cvFXL5cTIcONZ+IMSVzflbPPtuL//73Zy5dKjB12MIe\nmaE6fvXVV8yaNYvLly8zdOhQgoKC2LFjB1lZWUybNk3/AOb333+fp556itLSUry9vfn4448tHar9\nGTJkPU2aOLJmzXB27PiJTZvOsHhxf0pKdERHfweUF7KG6NmzHcXFNxt0xFgXps5p//6LRELine truncated
|
||||
"png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEKCAYAAAAPVd6lAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8FeX1+PHPySrIJqhhF43Ysi8iqHxb0IKNWFG+btgW\nqeCKijuCqAS/v1rF1h0XqlixIKJVBBEVF6AuLCIg+1IIJCFBIYSlLIFwfn/MJNyEuTc3yV3jeb9e\n9+XcZ56Ze06CJ3OfeWZGVBVjjDHxKyHaARhjjKkeK+TGGBPnrJAbY0ycs0JujDFxzgq5McbEOSvk\nxhgT56yQG2NMnLNCbmKeiPxeRL4Tkb0isk1EPhKRnu66s0TkHRH5SUQKRWS5iNwtIvZv2/xs2D92\nE9NE5B7gaeD/AacCLYDxQH8RSQcWAluA9qraALgKOBuoG52IjYk8sSs7TawSkfpADvAnVf2Xx/p/\nAvVV9dKIB2dMDLEjchPLzgNOAN73s/43wLuRC8eY2GSF3MSyRsAOVT0aYH1eBOMxJiZZITexbCdw\ncoATlzuBphGMx5iYZIXcxLJvgUPAAD/rPwOuiFw4xsQmK+QmZqnqbuARYLyIXCYitUUkWUQuFpEn\ngDHA+SIyTkTSAETkTBF50z1RaszPghVyE9NU9SngHuAh4EdgKzAMeF9VN+GcEG0FrBKRQpyTn4uB\nvVEJ2JgosOmHxhgTYiKSATwDJAKvquoTHn2eAy4G9uNMsV1a1c+zI3JjjAkhEUkEXgAygLbAtSLS\nplyffsCZqtoauAl4qTqfaYXcGGNCqzuwUVWzVPUwMBW4rFyf/sAbAKq6EGhQcp6nKqyQG2NMaDUD\nsn3e57htFfVpXtUPtEJujDGhFeyJR6nidsdJquqGkSQidkbWGBM0VS1fJCulsjWn3Ofl4tzcrUQL\nnCNuAvRp7rZVSVwUckdmtAMAvgQuCMmeVMcAIDI2ajFUz/FxVD2n0MUQDiV5+ZOZmUlmZmbY44j1\nGGIlDpFq1fBS/y/Ifg8d3/Qd0FpEWgHbgGuAa8v1mQHcDkwVkXOBQlXdXsVQ46mQG2NM5CRXcTtV\nPSIitwOf4Ew/fE1V14jIze76V1T1IxHpJyIbgf8C11cnVivkxhjjoTrFUVVnA7PLtb1S7v3t1fiI\nMqyQV0qraAdAbMQAsRFHq2gHAEDv3r2jHUJMxACxE0co1Ip2AJVghbxSTo92AMRGDBAbccRCDLFR\nvGIhBoidOEKhqkMr0WCF3BhjPMRTcYynWI0xJmLsiNwYY+JcPBXHeIrVGGMixo7IjTEmzlkhN8aY\nOGfTD40xJs7FU3GMp1iNMSZibGjFGGPiXDwVx3iK1RhjIsaOyI0xJs7FU3GMp1iNMSZi4umI3B71\nZowxHmoF+aoMEWkoInNEZL2IfCoiDfz0GyUiq0RkhYhMEZHUQPu1Qm6MMR6Sg3xV0khgjqqeBXzu\nvi/DfbLQjUBXVe2A83CKgYF2aoXcGGM8JAX5qqT+wBvu8hvA5R599gCHgdoikgTUpoLneQYs5CIy\nUUS2i8gKn7YnRWSNiCwXkfdEpL7PulEiskFE1orIRT7tZ7tfETaIyLM+7aki8rbbvkBETgsUjzHG\nREpyUnCvSkrzeTbndiCtfAdVLQD+BmzFeeZnoap+FminFR2Rvw5klGv7FGinqp2A9cAoABFpi/OQ\n0bbuNi/KsaegvgQMVdXWOA8lLdnnUGCn2/408EQF8RhjTEQkJXm/vhV4Uo+9ynPHwFd4vPr79lNV\nBY7bg4ikA3fhPAKrKVBHRP4QMNZAK1X13+54jW/bHJ+3C4Er3OXLgLdU9TCQ5T5UtIeIbAHqquoi\nt98knK8TH+N8zSh5RPm/gBcCxWOMMZGSnOjdfmEiXOjz/vH/ll2vqn397dMd4Wisqvki0gT40aNb\nN+AbVd3pbvMecD4w2d9+qztGPgT4yF1uCuT4rMsBmnm057rtuP/NBufJ08BuEWlYzZiMMaba/B2R\nl39V0gxgsLs8GJju0WctcK6I1HJHNfoAqwPttMqFXERGA0WqOqWq+zDGmFiVnBrcq5IeB/qKyHqc\nA/vHAUSkqYjMAlDV5TgjF98BP7jbTQi00ypdECQifwL6Ab/xac4FWvi8b45zJJ7rLpdvL9mmJbDN\nPTtb3x3o9/Clz3IrYuXBu8aY6Jo7dy5z584N/Y7DcLmkW9/6eLRvAy7xeT8OGBfsfisdqnui8n6g\nl6oe9Fk1A5giIk/hDJm0BhapqorIHhHpASwCBgHP+WwzGFgAXIkzr9KPCyobqjHmZ6B379707t27\n9P3YsWNDs+M4uu49YKgi8hbQCzhZRLJxTkyOAlKAOe6klG9VdZiqrhaRaThjOUeAYe5ZWYBhwD9w\nLoT6SFU/dttfA94UkQ3ATiqY9G6MMRFTUwq5ql7r0TwxQP/HgMc82pcAHTzaDwFXVxymMcZEmJ9Z\nK7Eojv7mGGNMBMVRdYyjUI0xJoIqPyMlaqyQG2OMlziqjnEUqjHGRFAcVcc4CtUYYyLITnYaY0yc\ni6PqGEehGmNMBMVRdYyjUI0xJoLiqDrGUajGGBNBNv3QGGPiXBxVR3tmpzHGeEkM8lUJInKViKwS\nkWIR6RqgXwMRedd9rOZqETk30H6tkHvIz7+Xrl2bADBv3p8YOLB96bq2bU9h2rQrWbfudo4ceZgJ\nEy49bvtevU6juPiR417XX985YjmUV92cABIShAce6Mnatbdx4MBo8vPvZfz4fhGJ35/q5vX665d5\n/q6OHHmYRo1qRSwPE4PC8/TlFcAAYH4F/Z7FucFgG6AjsKaiUI2P9PSTqF07maVL80hOTqBbt6Z8\n9dXW0vW1aiWRlbWbDz5Yxz33nMexGzwer0uXV8jL21v6fs+eQ2GN3Z9Q5fSPf1xGjx7NGTFiDsuW\n5VO3biqtWjWIVBrHCUVew4fPZsSIY08vFBGmT7+GffuK2LnzQETyMDEqPPcjXwvOvzN/3Afa/0pV\nB7vbHAF2B9qvFfJyevZsycKFuajCOec0Y+fO/eTk7Cldv2RJHkuW5AEwdGiXgPvasWM/P/20P6zx\nBiMUOfXu3YqBA9vTsePLrF27o7R95UqvRw5GRijy2ru3iL17i0rft27dkB49mnPVVe+EN3gT+6J3\nQdDpwE8i8jrQCVgC3KmqfouJFXLXrl0PoKqkpiaRkCAUFIwgOTmR1NRECgpGoAqNGgX9wA4Avvrq\nemrXTmbjxgJeeWUJb775Q8UbhVAoc7riijZs2rSLiy5KZ+bMa0lJSeTbb7O57745ZYpnJITjd1Xi\n5pu7kZ+/j+nT14Y4ahN3/FTHudtgbp7/zURkDtDYY9WDqjozyE/uCtyuqotF5BlgJPBIJUP9+enY\n8SVEhAULhnLLLbNYtiyfqVOvYMqUlXzwQeX+p962bS/Dhs3iu++2cfSo0q9fayZMuJQzz2zImDFz\nw5OAh1DmlJ5+Ei1b1uePf+zA0KEzKCoq5s9/vpAvvriO9u1foqioOExZHC+UeflKSUlk8OBOvPLK\nEo4e9T9kZn4mTvBu7n2G8yoxdmnZ9arat5qfnAPkqOpi9/27OIXcLyvkruzsPXTocCrJyYnMnLmO\nOnVS6Ny5Mf37T2XHjsoNj2zYUMCGDccePbp0aT6JiQncd995jB07L2JFIpQ5JSQIqalJXHfd9NKh\nlWuueZe8vHvp1691RI9gQ5mXryuvbMtJJ53AhAlLQhitiVvhH1rxHChX1XwRyRaRs1R1Pc4zPlcF\n2pEVcmDlyltp2bI+SUkJJCcnsnv3yNLCtWnTcADatBlPbu7eCvbk38KFOZx4YgqnnFKb7dv/G6rQ\n/Qp1Tnl5+1DVMuPjO3bsZ8eO/bRsWT8sOXgJ5+/qllvO5pNP/sPWrQHPK5mfizBURxEZgPPM4pOB\nWSKyVFUvFpGmwN9VteQBzHcAk0UkBfgPcH2EQ40/GRmTSUlJZOLE/syevZFp01YxZkwvDh0q5vHH\nvwKcQlYdXbs2Yf/+w9U6YqyMUOc0f/4WrruuE2ed1Yj163cC0LBhLU4+uTZZWYVhycFLuH5Xbdqc\nTM+eLRkw4O1Qh2ziVXhmrbwPvO/Rvg24xOf9cuCcYPdrhRzIydlDQoLQsWMaN930IZs3F9KhQxqZ\nmXPZvLlskUpKSqBdu1MAqFs3lUaNatGpUxpFRcWsWeMcrd5117ls2VLI6tU/oQq//W06o0f/ihde\nWExxcWSGVUKd01tvrWT06F8xcWJ/hg//mMOHi3niiT5s2LCT2bM3RCSncORV4uabu7Ft215mzlwX\nsVxMjLPb2MafLl0ac+hQMevX76RevVTatTuF+fO3HNevWbO6fP/9zQCoKl27NmHAgDZkZRWSnv4c\nAImJwmOP/YYWLepx+PBRNmzYyfDhHzNx4tLj9hcvOR08eIQ+fd7kmWd+y9y5g9m//zBz52bRp8+b\nHD58NG7zAjjhhCQGDerI888vIsBlAebnJo6qowS6oCVWiIhCZrTDCCnVMQCIjI1yJKFTE3OCY3mZ\n+CAiqKr/K26C24fqo0H2fYRqf151xdHfHGOMiSC7+6ExxsS5OKqOAYdWRGQizpnUH1W1g9vWEHgb\nOA3IAq5W1UJ33ShgCFAMDFfVT932s4F/4Eyx/0hV73TbU4FJOFcx7QSuUdXjBjtFRONhCMgYE30h\nG1p5Nsi+d0Z/aKWiux++DmSUaxsJzFHVs4DP3feISFvgGqCtu82LcuzOMC8BQ1W1NdBaREr2ORTY\n6bY/DTxRzXyMMSY0wnAb23AJ+OVBVf8tIq3KNfcHernLbwBzcYr5ZcBbqnoYyBKRjUAPEdkC1FXV\nRe42k4DLgY/dfZWcSfoX8IK/WGrqCbSalFdNzAlqZl41MScI8YnpOBpaqUqoaaq63V3eDqLine truncated
|
||||
"text": [
|
||||
"<matplotlib.figure.Figure at 0x10ae16e90>"
|
||||
"<matplotlib.figure.Figure at 0x7f55de17de90>"
|
||||
]
|
||||
}
|
||||
],
|
||||
"prompt_number": 34
|
||||
"prompt_number": 4
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"# Get the mass matrix \n",
|
||||
"# The model\n",
|
||||
"Msig = M.getEdgeInnerProduct(sig)\n",
|
||||
"MsigBG = M.getEdgeInnerProduct(sigBG)\n",
|
||||
"\n",
|
||||
@@ -108,7 +112,7 @@
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": 38
|
||||
"prompt_number": 5
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
@@ -122,78 +126,128 @@
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": 39
|
||||
"prompt_number": 6
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"def get1Dfields(sigHalf, freq):\n",
|
||||
" hz = [(100.,18)]\n",
|
||||
" M = Mesh.TensorMesh([hz],'C')\n",
|
||||
" sig = np.zeros(M.nC) + 1e-8\n",
|
||||
" sig[M.vectorCCx<=0] = sigHalf\n",
|
||||
" G = M1D.nodalGrad\n",
|
||||
" Mmu = Utils.sdiag(M.vol*(1.0/mu_0))\n",
|
||||
" Msig = M.getFaceInnerProduct(sig)\n",
|
||||
" A = G.T*Mmu*G - 1j*omega(freq)*Msig\n",
|
||||
" Aii = A[1:-1,1:-1]\n",
|
||||
" Aio = A[1:-1,[0,-1]]\n",
|
||||
" bc = np.r_[0.0,1.0]\n",
|
||||
" rhs = -Aio*bc\n",
|
||||
" eii = Solver(Aii).solve(rhs)\n",
|
||||
" e = np.r_[bc[0],eii,bc[1]]\n",
|
||||
" return e"
|
||||
"# Need to solve x and y polarizations of the source.\n",
|
||||
"from simpegMT.Utils import get1DEfields\n",
|
||||
"# Get a 1d solution for a halfspace background\n",
|
||||
"mesh1d = simpeg.Mesh.TensorMesh([M.hz],np.array([M.x0[2]]))\n",
|
||||
"e0_1d = get1DEfields(mesh1d,M.r(sigBG,'CC','CC','M')[0,0,:],freq)\n",
|
||||
"# Setup x (east) polarization (_x)\n",
|
||||
"ex_x = np.zeros(M.vnEx,dtype=complex)\n",
|
||||
"ey_x = np.zeros((M.nEy,1),dtype=complex)\n",
|
||||
"ez_x = np.zeros((M.nEz,1),dtype=complex)\n",
|
||||
"# Assign the source to ex_x\n",
|
||||
"for i in arange(M.vnEx[0]):\n",
|
||||
" for j in arange(M.vnEx[2]):\n",
|
||||
" ex_x[i,j,:] = e0_1d\n",
|
||||
"eBG_x = np.vstack((simpeg.Utils.mkvc(M.r(ex_x,'Ex','Ex','V'),2),ey_x,ez_x))\n",
|
||||
"rhs_x = ABG.dot(eBG_x)"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": 89
|
||||
"prompt_number": 50
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"x0 = -h1.sum()/2"
|
||||
"# Setup y (north) polarization (_y)\n",
|
||||
"ex_y = np.zeros((M.vnEx,1))\n",
|
||||
"ey_y = np.zeros(M.vnEy)\n",
|
||||
"ez_y = np.zeros((M.vnEz,1))\n",
|
||||
"# Assign the source to ex_x\n",
|
||||
"for i in arange(M.vnEy[0]):\n",
|
||||
" for j in arange(M.vnEy[2]):\n",
|
||||
" ey_y[i,j,:] = e0_1d\n",
|
||||
"eBG_y = np.vstack((ex_y,simpeg.Utils.mkvc(M.r(ey_y,'Ey','Ey','V'),2),ez_y))\n",
|
||||
"rhs_y = ABG.dot(eBG_y)"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": 87
|
||||
"prompt_number": 12
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"sigs = np.logspace(-3,1,4)\n",
|
||||
"for s in sigs:\n",
|
||||
" e = get1Dfields(s,1e1)\n",
|
||||
" plot(e.imag,M.vectorNx)\n",
|
||||
"ylim([-900,0])\n",
|
||||
"ylabel('Depth, m')\n",
|
||||
"legend(['$\\sigma = 1e%d$'%log10(s) for s in sigs],'best')"
|
||||
"# Solve the systems for each polarization\n",
|
||||
"Ainv = simpeg.SolverCG(A)\n",
|
||||
"e_x = Ainv*rhs_x\n",
|
||||
"\n",
|
||||
"e_y = Ainv*rhs_x"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": []
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"prompt_number": 28
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"\n"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"metadata": {},
|
||||
"output_type": "pyout",
|
||||
"prompt_number": 101,
|
||||
"text": [
|
||||
"<matplotlib.legend.Legend at 0x10c40d750>"
|
||||
]
|
||||
},
|
||||
{
|
||||
"metadata": {},
|
||||
"output_type": "display_data",
|
||||
"png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEACAYAAACUMoD1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XlYVOX7x/H3KLiL+1IzKSYogqikorZShrmkqZWGFVJZ\n5l6m2e5Siqa2mJFWal/J3EuslMCFckUFd1xQsWAAV1REZD2/P47wcwEZYGbODNyv6+JKD2fO3Egz\nn3mW8zw6RVEUhBBCiFKqoHUBQgghygYJFCGEEGYhgSKEEMIsJFCEEEKYhQSKEEIIs5BAEUIIYRZ2\nFyihoaG4ubnh6urKjBkztC5HCCHEDTp7ug8lJyeHli1bsmHDBvR6PR07dmTp0qW0atVK69KEEKLc\ns6sWyq5du3BxccHZ2RlHR0deeOEFQkJCtC5LCCEEdhYoRqOR++67L//vBoMBo9GoYUVCCCHy2FWg\n6HQ6rUsQQghRCAetCygOvV5PfHx8/t/j4+MxGAy3nOPi4sLJkyetXZoQQti15s2bc+LEiVJdw65a\nKB06dCA2NpbTp0+TmZnJ8uXL6dOnzy3nnDx5EkVR7PZr4sSJxTrfyUkhJUX7uu9We25uLuEXLvD0\ngQPU37qV906e5L/0dM3rLe2/va19Sf12UH92NtmvvsLe+6vzw6ZZmtd885c5PojbVQvFwcGBuXPn\n8tRTT5GTk8Nrr71W7md4VasG165B7dpaV3Knazk5BJ85w5yEBCrodIzW61nu7k61ihW1Lk0I68vO\nhldeoWJCAjU2beX9Fb50aNWVdo3baV2Z2dhVoAD06NGDHj16aF2GzcgLFFvy7/XrBBmNLExO5kEn\nJ75xdeXx2rVlDEyUX5mZ8OKLkJoKf/6JS7VqfPXUV7yw6gX2vLGHGpVqaF2hWdhVl1d54OPjU6zz\nbSVQFEWhbocOPHfoEA/s2UOWohD5wAOEeHryRJ06dhEmxf23tzVSv7YKrf/6dXjuOTVUQkLUFy3w\nYpsX6XJfF0auG2m9Ii3Mrm5sNIVOp6OM/Uh35e0N33wDnTpp8/zXc3JYdvYsc4xG0nJyGKXXM7hx\nY2o62F3jVwjzu3YN+vWDWrVgyRJwdLzl22mZabT/vj0fPvIhL7d9WaMiVeZ475RXvZ3TqoWSmJHB\nd4mJfJ+YiFeNGkxt1oyn6talgh20RISwitRU6N0bmjSBhQuhgA9Z1StVZ8XzK+i6uCudDJ1oUa+F\nBoWaj3R52TlrB8quK1d4MSYGj927uZiVxd/t2hHati096tWTMBEiz+XL8NRT0KIF/PRTgWGSp02j\nNkzxmcLAVQPJyM6wXo0WIIFi56wRKFm5uSw9c4bOUVEMjImhfc2axHXqxLctWuBWvbpln1wIe3Ph\nAnTtCh06wPz5UKHot9k3O7xJ8zrNGR8+3goFWo50edk5SwbKucxM5icm8l1iIi2qVeO9Jk3oXb8+\nFaUlIkTBzpwBX1/o0QOmTwcTXys6nY4fev+A13wvujbryjNuz1i4UMuQQLFz1apBerp5r7kvNZU5\nRiO/nT/Ps/Xrs75NG9rUKBvTGoWwGKMRnnwSBg6EiRNNDpM8darWYemzS+m7vC9e93jRpFYTCxVq\nORIodq5qVfO0ULJzc1l74QJfJyRwMj2dEXo9sd7e1K9UqfQXF6Ks+/dftZvr9ddhwoQSX6bLfV0Y\n23ksg1YPIiIgAocK9vUWLWModq60XV4pWVnM/O8/XCIjmR0fz3C9nrjOnXm/aVMJEyFMceIEPPYY\njBpVqjDJM/6h8VRzrMakiEmlr83K7Cv+xB2qVYO0tOI/LiYtjW+MRpadPcvT9eqx0sODjk5O5i9Q\niLLs6FF1zOTjj+GNN8xyyQq6CgT3C8ZrvhePOz9O1/u7muW61iAtFDtXnBZKrqLw54ULdNu/nyf2\n7aORoyMxHTsS3KqVhIkQxXXgADzxBEydarYwydOoRiMW91uM/xp/zqadNeu1LUlaKHbOlEH5K9nZ\n/JSczDdGI04VKzLGYGCgpyeVTZjOKIQoQFQU9OoFX3+tDsJbwJP3P8ngtoPx/82fdS+uo4LO9l+v\ntl+huKu7tVBOXLvGmNhYnHfuZOvly/zk5sae9u3xb9xYwkSIktqxQ50WPG+excIkz2SfyaRmpjJ7\n+2yLPo+5SAvFzt0+y0tRFDakpDDHaGTnlSsMuece9nfowH1VqmhXpBBlRUQEDBgAixdD9+4WfzrH\nio780v8XvH/05pGmj9DZ0Nniz1kaEih2Lq+FkpaTw8+y94gQlhMWpi5Bv3y5OnZiJU1rN2Ver3n4\nrfZj79C91K5ig5sf3SCrDdu5ZRHXGb/VSLpPEg/VqsUYg0H2HhHC3H7/HV57DX77DR56SJMSRq4b\nyZm0M6x4boVFXt/meO/UpCN95cqVeHh4ULFiRaKjo2/5XmBgIK6urri5uREWFpZ/PCoqCk9PT1xd\nXRkzZoy1S7YpiqLwz6VLPHvoEEPZQ1auwq727e1q7xEh7MbKlTBkCPz5p2ZhAjCr2yxiL8TyfdT3\nmtVQJEUDR44cUY4dO6b4+PgoUVFR+ccPHz6stG3bVsnMzFTi4uKU5s2bK7m5uYqiKErHjh2VyMhI\nRVEUpUePHsr69esLvLZGP5JVpGdnK4sSE5V2u3crLXbuVL6Jj1d27s9SPDy0rkyIMio4WFEaN1aU\nffu0rkRRFEU5eu6oUv/z+sqB5ANmv7Y53js1aaG4ubnRosWd6/6HhITg5+eHo6Mjzs7OuLi4EBkZ\nSVJSEqmpqXh7ewPg7+/PmjVrrF22ZhIzMvg4Lo6mO3ey7OxZpjVrxhFvb0YaDDSo4WATOzYKUeb8\n+KN65/vGjdC2rdbVANCyfktmd5vNwFUDScsswR3NFmZTc0cTExMxGAz5fzcYDBiNxjuO6/V6jEaj\nFiVa1a4rVxhUxN4j5lrLSwhxk7lzYcoUdVaXu7vW1dzCv60/He7twOj1o7Uu5Q4Wm+Xl6+tLcnLy\nHcenTZtG7969LfW0AEyaNCn/zz4+Pna5V/XRtDS6HzjAh02bEuTqSu3btg7NU6WK+VcbFqJcungR\nli6FBQvUT2l//w3NmmldVYGCegXhMseFExdP4FLXpUTXiIiIICIiwqx1WSxQwsPDi/0YvV5PfHx8\n/t8TEhIwGAzo9XoSEhJuOa7X6wu9zs2BYq+irl7Ft04d3rnvvrueV6EClKNJbUKYV26u2qW1cCGs\nX6/esDhjhjot2Ian3NeoVINmdZpxNu1siQPl9g/bkydPLnVdmnd5KTe9G/bp04dly5aRmZlJXFwc\nsbGxeHt707hxY5ycnIiMjERRFIKDg+nbt6+GVVvegatX8TRhN0QJFCFKIC5O3bOkWTN47z14+GE4\ndUptofj62nSY5KlTpQ4p6Slal3ELTQLlt99+47777mPnzp306tWLHj16AODu7s6AAQNwd3enR48e\nBAUF5U+BDQoKYsiQIbi6uuLi4kJ3K9ylqqWDaWl4mrCplU6nfsgSQhQhPR2WLFH3LenYEVJSICRE\nXZdrxAioW1frCoulTtU6pFy3rUCRGxtt1H07dvB3u3bcX7XqXc9LS4MGDWRgXogCKYoaGAsXqne4\nd+wIr74KffqoA5B2bOS6kbSs15JRnUaZ5XrmeO+UpVdsUEpWFpeys3E24X94nU66vIS4w7lzamtk\n4UK4elUNkX37oIgxSXtSp4rttVAkUGzQobQ0PKpVy58afDcSKELckJ2trre1cCFs2KC2QubMgUcf\nVQcby5g6VesQfzm+6BOtSALFBpk6fgISKEJw4oQaIv/7HxgMamtkwQKoVUvryiyqTpU6HDhzQOsy\nbiGBYoMOpqWZNMMLZJaXKKfS0mDVKjVIjh6Fl16Cv/6C1q21rsxqbHFQXgLFBh28epUBDRqYdK7M\n8hLlhqLAzp1qiKxapU71festdefESpW0rs7qbHHasASKjVEUhUPFaKFIl5co85KTIThYDZKcHLVL\n6/BhuPderSvTlLRQRJHiMzKoWrEi9U38xCWBIsqkrCz1zvWFC9UlUPr1gx9+UJePl+0ZAGmhCBMU\nZ/wEJFBEGXPkCCxapG6x6+KitkaCg6FmTa0rsznSQhFFOmjikit55MOasHtXrsCKFWprJC4OBg9W\nWyUtW2pdmU2r7lidrJwsMnMyqVTRNsaQJFBszMG0NHzr1DH5/LxAURQJF2FHFAW2blWn965ZA48/\nDu+/ry7O6CBvS6bQ6XTUrlKblPQUGtVopHU5gA0sDiluVZx7UPLITC9hN4xGCAyEFi1g6FDw9IRj\nx9S92nv3ljApJlvr9pLfng3Jys0lNj0d92rVivU4GUcRNi0zE37/Xe3S2rEDnn8efv4ZvL2lWV1K\ndarU4WL6Ra3LyCeBYkOOXbtGk8qVqVrMpbMlUIRNOnRIDZGff1Z3PXztNXWspBhjhOLu6lS1rZle\nEig2pLgzvPJIoAibcekSLFumBkliIgQEwPbt6owtYXa2tkCkBIoNKcn4CUigCI3l5qp7ry9cCH/8\nAd26qfux28lGVfbM1u5FkUCxIQfT0nilceNiP07W8xKa+O8/dUHGRYvU+0Reew2++grq19e6snLD\n1gblNZnlNX78eFq1akXbtm3p378/ly9fzv9eYGAgrq6uuLm5ERYWln88KioKT09PXF1dGTNmjBZl\nW1xx70HJI7O8hNVcv65uVPXUU+DlpS6LsnKlutfI6NESJlZmay0UTQKlW7duHD58mP3799OiRQsC\nAwMBiImJYfny5cTExBAaGsrw4cPzdxAbNmwYCxYsIDY2ltjYWEJDQ7Uo3WKuZGdzLiuryB0aCyJd\nXsLi9u6FUaPU5eF//FEdG0lIgG+/hfbtZbaWRqSFAvj6+lLhxoY3nTp1IiEhAYCQkBD8/PLine truncated
|
||||
"text": [
|
||||
"<matplotlib.figure.Figure at 0x10c969890>"
|
||||
"ename": "NameError",
|
||||
"evalue": "name 'e0_1d' is not defined",
|
||||
"output_type": "pyerr",
|
||||
"traceback": [
|
||||
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m\n\u001b[1;31mNameError\u001b[0m Traceback (most recent call last)",
|
||||
"\u001b[1;32m<ipython-input-36-4c5017fec386>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m 5\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[1;32min\u001b[0m \u001b[0marange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mM\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mvnEx\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\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 6\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mj\u001b[0m \u001b[1;32min\u001b[0m \u001b[0marange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mM\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mvnEx\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[0m\n\u001b[1;32m----> 7\u001b[1;33m \u001b[0mex_x\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0me0_1d\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
|
||||
"\u001b[1;31mNameError\u001b[0m: name 'e0_1d' is not defined"
|
||||
]
|
||||
}
|
||||
],
|
||||
"prompt_number": 101
|
||||
"prompt_number": 36
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"collapsed": false,
|
||||
"input": [
|
||||
"# Get a 1d solution for a halfspace background\n",
|
||||
"mesh1d = simpeg.Mesh.TensorMesh(M.hz,M.x0[2])\n",
|
||||
"e0_1d = get1DEfields(mesh1d,M.r(sigBG,'CC','CC','M')[0,0,:],freq)\n",
|
||||
"# Setup x (east) polarization (_x)\n",
|
||||
"ex_x = np.zeros(M.vnEx)\n",
|
||||
"ey_x = np.zeros(M.vnEy)\n",
|
||||
"ez_x = np.zeros(M.vnEz)\n",
|
||||
"# Assign the source to ex_x\n",
|
||||
"for i in arange(M.vnEx[0]):\n",
|
||||
" for j in arange(M.vnEx[2]):\n",
|
||||
" ex_x[i,j,:] = e0_1d\n",
|
||||
"eBG_x = np.vstack(M.r(ex_x,'Ex','Ex','V'),M.r(ex_y)"
|
||||
],
|
||||
"language": "python",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"ename": "SyntaxError",
|
||||
"evalue": "invalid syntax (<ipython-input-34-f14ef29654fd>, line 12)",
|
||||
"output_type": "pyerr",
|
||||
"traceback": [
|
||||
"\u001b[1;36m File \u001b[1;32m\"<ipython-input-34-f14ef29654fd>\"\u001b[1;36m, line \u001b[1;32m12\u001b[0m\n\u001b[1;33m eBG_x = np.vstack(M.r(ex_x,'Ex','Ex','V'),M.r(ex_y)\u001b[0m\n\u001b[1;37m ^\u001b[0m\n\u001b[1;31mSyntaxError\u001b[0m\u001b[1;31m:\u001b[0m invalid syntax\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"prompt_number": 34
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
|
||||
@@ -0,0 +1,101 @@
|
||||
# Analytic solution of EM fields due to a plane wave
|
||||
|
||||
import numpy as np, SimPEG as simpeg
|
||||
|
||||
def getEHfields(m1d,sigma,freq,zd):
|
||||
'''Analytic solution for MT 1D layered earth. Returns E and H fields.
|
||||
|
||||
:param SimPEG.mesh, object m1d: Mesh object with the 1D spatial information.
|
||||
:param numpy.array, vector sigma: Physical property of conductivity corresponding with the mesh.
|
||||
:param float, freq: Frequency to calculate data at.
|
||||
:param numpy array, vector zd: location to calculate EH fields at
|
||||
|
||||
Assumes a halfspace with the same conductive as the last cell below.
|
||||
|
||||
'''
|
||||
# Note add an error check for the mesh and sigma are the same size.
|
||||
# Need make the solution e^-iwt dependent
|
||||
|
||||
# Constants: Assume constant
|
||||
mu = 4*np.pi*1e-7*np.ones((m1d.nC+1))
|
||||
eps = 8.85*1e-12*np.ones((m1d.nC+1))
|
||||
# Angular freq
|
||||
w = 2*np.pi*freq
|
||||
# Add the halfspace value to the property
|
||||
sig = np.concatenate((np.array([sigma[0]]),sigma))
|
||||
# Calculate the wave number
|
||||
k = np.sqrt(eps*mu*w**2-1j*mu*sig*w)
|
||||
|
||||
# Initiate the propagation matrix, in the order down up.
|
||||
UDp = np.zeros((2,m1d.nC+1),dtype=complex)
|
||||
UDp[1,0] = 1 # Set the wave amplitude as 1 into the half-space at the bottom of the mesh
|
||||
# Loop over all the layers, starting at the bottom layer
|
||||
for lnr, h in enumerate(m1d.hx): # lnr-number of layer, h-thickness of the layer
|
||||
# Calculate
|
||||
yp1 = k[lnr]/(w*mu[lnr]) # Admittance of the layer below the current layer
|
||||
zp = (w*mu[lnr+1])/k[lnr+1] # Impedance in the current layer
|
||||
# Build the propagation matrix
|
||||
|
||||
# Convert fields to down/up going components in layer below current layer
|
||||
Pj1 = np.array([[1,1],[yp1,-yp1]])
|
||||
# Convert fields to down/up going components in current layer
|
||||
Pjinv = 1./2*np.array([[1,zp],[1,-zp]])
|
||||
# Propagate down and up components through the current layer
|
||||
elamh = np.array([[np.exp(-1j*k[lnr+1]*h),0],[0,np.exp(1j*k[lnr+1]*h)]])
|
||||
|
||||
# The down and up component in current layer.
|
||||
UDp[:,lnr+1] = elamh.dot(Pjinv.dot(Pj1)).dot(UDp[:,lnr])
|
||||
|
||||
# Calculate the fields
|
||||
Ed = np.zeros((zd.size,),dtype=complex)
|
||||
Eu = np.zeros((zd.size,),dtype=complex)
|
||||
Hd = np.zeros((zd.size,),dtype=complex)
|
||||
Hu = np.zeros((zd.size,),dtype=complex)
|
||||
|
||||
# Loop over the layers and calculate the fields
|
||||
for ki,mui,epsi,dlow,dup,Up,Dp in zip(k[1::],mu[1::],eps[1::],m1d.vectorNx[:-1],m1d.vectorNx[1::],UDp[0,1::],UDp[1,1::]):
|
||||
dind = np.logical_and(dup >= zd, zd > dlow)
|
||||
Ed[dind] = Dp*np.exp(-1j*ki*(dup-zd[dind]))
|
||||
Eu[dind] = Up*np.exp(1j*ki*(dup-zd[dind]))
|
||||
Hd[dind] = (ki/(w*mui))*Dp*np.exp(-1j*ki*(dup-zd[dind]))
|
||||
Hu[dind] = -(ki/(w*mui))*Up*np.exp(1j*ki*(dup-zd[dind]))
|
||||
|
||||
# Return return the fields
|
||||
return Ed, Eu, Hd, Hu
|
||||
|
||||
def getImpedance(m1d,sigma,freq):
|
||||
"""Analytic solution for MT 1D layered earth. Returns the impedance at the surface.
|
||||
|
||||
:param SimPEG.mesh, object m1d: Mesh object with the 1D spatial information.
|
||||
:param numpy.array, vector sigma: Physical property corresponding with the mesh.
|
||||
:param numpy.array, vector freq: Frequencies to calculate data at.
|
||||
|
||||
|
||||
"""
|
||||
|
||||
# Define constants
|
||||
mu0 = 4*np.pi*1e-7
|
||||
eps0 = 8.85e-12
|
||||
|
||||
# Initiate the impedances
|
||||
Z1d = np.empty(len(freq) , dtype='complex')
|
||||
h = m1d.hx #vectorNx[:-1]
|
||||
# Start the process
|
||||
for nrFr, fr in enumerate(freq):
|
||||
om = 2*np.pi*fr
|
||||
Zall = np.empty(len(h)+1,dtype='complex')
|
||||
# Calculate the impedance for the bottom layer
|
||||
Zall[0] = (mu0*om)/np.sqrt(mu0*eps0*(om)**2 - 1j*mu0*sigma[0]*om)
|
||||
|
||||
for nr,hi in enumerate(h):
|
||||
# Calculate the wave number
|
||||
# print nr,sigma[nr]
|
||||
k = np.sqrt(mu0*eps0*om**2 - 1j*mu0*sigma[nr]*om)
|
||||
Z = (mu0*om)/k
|
||||
|
||||
Zall[nr+1] = Z *((Zall[nr] + Z*np.tanh(1j*k*hi))/(Z + Zall[nr]*np.tanh(1j*k*hi)))
|
||||
|
||||
#pdb.set_trace()
|
||||
Z1d[nrFr] = Zall[-1]
|
||||
|
||||
return Z1d
|
||||
Binary file not shown.
@@ -0,0 +1,43 @@
|
||||
import numpy as np, SimPEG as simpeg
|
||||
from MT1Danalytic import getEHfields
|
||||
from scipy.constants import mu_0
|
||||
|
||||
def get1DEfields(m1d,sigma,freq,sourceAmp=1.0):
|
||||
"""Function to get 1D electrical fields"""
|
||||
|
||||
# Get the gradient
|
||||
G = m1d.nodalGrad
|
||||
# Mass matrices
|
||||
# Magnetic permeability
|
||||
Mmu = simpeg.Utils.sdiag(m1d.vol*(1.0/mu_0))
|
||||
# Conductivity
|
||||
Msig = m1d.getFaceInnerProduct(sigma)
|
||||
# Set up the solution matrix
|
||||
A = G.T*Mmu*G - 1j*2.*np.pi*freq*Msig
|
||||
# Define the inner part of the solution matrix
|
||||
Aii = A[1:-1,1:-1]
|
||||
# Define the outer part of the solution matrix
|
||||
Aio = A[1:-1,[0,-1]]
|
||||
|
||||
# Set the boundary conditions
|
||||
Ed_low, Eu_low, Hd_low, Hu_low = getEHfields(m1d,sigma,freq,np.array([m1d.vectorNx[0]]))
|
||||
Etot_low = Ed_low + Eu_low
|
||||
## Note: need to use conjugate of the analytic solution. It is derived with e^iwt
|
||||
bc = np.r_[Etot_low.conj(),sourceAmp]
|
||||
# The right hand side
|
||||
rhs = -Aio*bc
|
||||
# Solve the system
|
||||
Aii_inv = simpeg.Solver(Aii)
|
||||
eii = Aii_inv*rhs
|
||||
# Assign the boundary conditions
|
||||
e = np.r_[bc[0],eii,bc[1]]
|
||||
# Return the electrical fields
|
||||
return e
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
hz = [(100.,18)]
|
||||
M = simpeg.Mesh.TensorMesh([hz],'C')
|
||||
sig = np.zeros(M.nC) + 1e-8
|
||||
sig[M.vectorCCx<=0] = sigHalf
|
||||
Binary file not shown.
@@ -0,0 +1 @@
|
||||
from MT1Dsolutions import * # Add the names of the functions
|
||||
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
@@ -0,0 +1,100 @@
|
||||
# Analytic solution of EM fields due to a plane wave
|
||||
|
||||
import numpy as np, SimPEG as simpeg
|
||||
|
||||
def getEHfields(m1d,sigma,freq,zd):
|
||||
'''Analytic solution for MT 1D layered earth. Returns E and H fields.
|
||||
|
||||
:param SimPEG.mesh, object m1d: Mesh object with the 1D spatial information.
|
||||
:param numpy.array, vector sigma: Physical property of conductivity corresponding with the mesh.
|
||||
:param float, freq: Frequency to calculate data at.
|
||||
:param numpy array, vector zd: location to calculate EH fields at
|
||||
|
||||
Assumes a halfspace with the same conductive as the last cell below.
|
||||
|
||||
'''
|
||||
# Note add an error check for the mesh.
|
||||
|
||||
# Constants: Assume constant
|
||||
mu = 4*np.pi*1e-7*np.ones((m1d.nC+1))
|
||||
eps = 8.85*1e-12*np.ones((m1d.nC+1))
|
||||
# Angular freq
|
||||
w = 2*np.pi*freq
|
||||
# Add the halfspace value to the property
|
||||
sig = np.concatenate((np.array([sigma[0]]),sigma))
|
||||
# Calculate the wave number
|
||||
k = np.sqrt(eps*mu*w**2-1j*mu*sig*w)
|
||||
|
||||
# Initiate the propagation matrix, in the order down up.
|
||||
UDp = np.zeros((2,m1d.nC+1),dtype=complex)
|
||||
UDp[1,0] = 1 # Set the wave amplitude as 1 into the half-space at the bottom of the mesh
|
||||
# Loop over all the layers, starting at the bottom layer
|
||||
for lnr, h in enumerate(m1d.hx): # lnr-number of layer, h-thickness of the layer
|
||||
# Calculate
|
||||
yp1 = k[lnr]/(w*mu[lnr]) # Admittance of the layer below the current layer
|
||||
zp = (w*mu[lnr+1])/k[lnr+1] # Impedance in the current layer
|
||||
# Build the propagation matrix
|
||||
|
||||
# Convert fields to down/up going components in layer below current layer
|
||||
Pj1 = np.array([[1,1],[yp1,-yp1]])
|
||||
# Convert fields to down/up going components in current layer
|
||||
Pjinv = 1./2*np.array([[1,zp],[1,-zp]])
|
||||
# Propagate down and up components through the current layer
|
||||
elamh = np.array([[np.exp(-1j*k[lnr+1]*h),0],[0,np.exp(1j*k[lnr+1]*h)]])
|
||||
|
||||
# The down and up component in current layer.
|
||||
UDp[:,lnr+1] = elamh.dot(Pjinv.dot(Pj1)).dot(UDp[:,lnr])
|
||||
|
||||
# Calculate the fields
|
||||
Ed = np.zeros((zd.size,),dtype=complex)
|
||||
Eu = np.zeros((zd.size,),dtype=complex)
|
||||
Hd = np.zeros((zd.size,),dtype=complex)
|
||||
Hu = np.zeros((zd.size,),dtype=complex)
|
||||
|
||||
# Loop over the layers and calculate the fields
|
||||
for ki,mui,epsi,dlow,dup,Up,Dp in zip(k[1::],mu[1::],eps[1::],m1d.vectorNx[:-1],m1d.vectorNx[1::],UDp[0,1::],UDp[1,1::]):
|
||||
dind = np.logical_and(dup >= zd, zd > dlow)
|
||||
Ed[dind] = Dp*np.exp(-1j*ki*(dup-zd[dind]))
|
||||
Eu[dind] = Up*np.exp(1j*ki*(dup-zd[dind]))
|
||||
Hd[dind] = (ki/(w*mui))*Dp*np.exp(-1j*ki*(dup-zd[dind]))
|
||||
Hu[dind] = -(ki/(w*mui))*Up*np.exp(1j*ki*(dup-zd[dind]))
|
||||
|
||||
# Return return the fields
|
||||
return Ed, Eu, Hd, Hu
|
||||
|
||||
def getImpedance(m1d,sigma,freq):
|
||||
"""Analytic solution for MT 1D layered earth. Returns the impedance at the surface.
|
||||
|
||||
:param SimPEG.mesh, object m1d: Mesh object with the 1D spatial information.
|
||||
:param numpy.array, vector sigma: Physical property corresponding with the mesh.
|
||||
:param numpy.array, vector freq: Frequencies to calculate data at.
|
||||
|
||||
|
||||
"""
|
||||
|
||||
# Define constants
|
||||
mu0 = 4*np.pi*1e-7
|
||||
eps0 = 8.85e-12
|
||||
|
||||
# Initiate the impedances
|
||||
Z1d = np.empty(len(freq) , dtype='complex')
|
||||
h = m1d.hx #vectorNx[:-1]
|
||||
# Start the process
|
||||
for nrFr, fr in enumerate(freq):
|
||||
om = 2*np.pi*fr
|
||||
Zall = np.empty(len(h)+1,dtype='complex')
|
||||
# Calculate the impedance for the bottom layer
|
||||
Zall[0] = (mu0*om)/np.sqrt(mu0*eps0*(om)**2 - 1j*mu0*sigma[0]*om)
|
||||
|
||||
for nr,hi in enumerate(h):
|
||||
# Calculate the wave number
|
||||
# print nr,sigma[nr]
|
||||
k = np.sqrt(mu0*eps0*om**2 - 1j*mu0*sigma[nr]*om)
|
||||
Z = (mu0*om)/k
|
||||
|
||||
Zall[nr+1] = Z *((Zall[nr] + Z*np.tanh(1j*k*hi))/(Z + Zall[nr]*np.tanh(1j*k*hi)))
|
||||
|
||||
#pdb.set_trace()
|
||||
Z1d[nrFr] = Zall[-1]
|
||||
|
||||
return Z1d
|
||||
Binary file not shown.
@@ -0,0 +1,42 @@
|
||||
import numpy as np, SimPEG as simpeg
|
||||
from MT1Danalytic import getEHfields
|
||||
from scipy.constants import mu_0
|
||||
|
||||
def get1DEfields(m1d,sigma,freq,sourceAmp=1.0):
|
||||
"""Function to get 1D electrical fields"""
|
||||
|
||||
# Get the gradient
|
||||
G = m1d.nodalGrad
|
||||
# Mass matrices
|
||||
# Magnetic permeability
|
||||
Mmu = simpeg.Utils.sdiag(m1d.vol*(1.0/mu_0))
|
||||
# Conductivity
|
||||
Msig = m1d.getFaceInnerProduct(sigma)
|
||||
# Set up the solution matrix
|
||||
A = G.T*Mmu*G - 1j*2.*np.pi*freq*Msig
|
||||
# Define the inner part of the solution matrix
|
||||
Aii = A[1:-1,1:-1]
|
||||
# Define the outer part of the solution matrix
|
||||
Aio = A[1:-1,[0,-1]]
|
||||
|
||||
# Set the boundary conditions
|
||||
Ed_low, Eu_low, Hd_low, Hu_low = getEHfields(m1d,sigma,freq,np.array([m1d.vectorNx[0]]))
|
||||
Etot_low = Ed_low + Eu_low
|
||||
bc = np.r_[Etot_low,sourceAmp]
|
||||
# The right hand side
|
||||
rhs = -Aio*bc
|
||||
# Solve the system
|
||||
Aii_inv = simpeg.Solver(Aii)
|
||||
eii = Aii_inv*rhs
|
||||
# Assign the boundary conditions
|
||||
e = np.r_[bc[0],eii,bc[1]]
|
||||
# Return the electrical fields
|
||||
return e
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
hz = [(100.,18)]
|
||||
M = simpeg.Mesh.TensorMesh([hz],'C')
|
||||
sig = np.zeros(M.nC) + 1e-8
|
||||
sig[M.vectorCCx<=0] = sigHalf
|
||||
Binary file not shown.
@@ -0,0 +1,2 @@
|
||||
from MT1Dsolutions import * # Add the names of the functions
|
||||
from MT1Danalytic import *
|
||||
Binary file not shown.
@@ -1,4 +1,4 @@
|
||||
# from EM import *
|
||||
import Utils
|
||||
import FDEM
|
||||
# import FDEM
|
||||
import Base
|
||||
Binary file not shown.
Reference in new issue
Block a user