Skip to content
This repository has been archived by the owner on Jul 12, 2023. It is now read-only.

Commit

Permalink
Merge 81983ec into 22fd00a
Browse files Browse the repository at this point in the history
  • Loading branch information
sgkang committed Aug 9, 2018
2 parents 22fd00a + 81983ec commit 2d5d01e
Show file tree
Hide file tree
Showing 10 changed files with 1,135 additions and 1,712 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -934,9 +934,9 @@
],
"metadata": {
"kernelspec": {
"display_name": "Python [conda env:simpeg]",
"display_name": "Python 3",
"language": "python",
"name": "conda-env-simpeg-py"
"name": "python3"
},
"language_info": {
"codemirror_mode": {
Expand Down
158 changes: 122 additions & 36 deletions notebooks/develop/Test_global_em1d_inversion_td_new.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,36 @@
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"from SimPEG import Mesh, Maps, Utils\n",
"import numpy as np\n",
"from matplotlib.colors import LogNorm\n",
"from simpegEM1D import (\n",
" GlobalEM1DProblemTD, GlobalEM1DSurveyTD, get_vertical_discretization_time, EM1DSurveyTD\n",
")\n",
"from pymatsolver import PardisoSolver\n",
"from simpegEM1D import skytem_HM_2015\n",
"\n",
"def get_xyz(\n",
" dx=100, \n",
" dy=100, \n",
" nx=5, \n",
" ny=5, \n",
" x0=np.r_[0, 0, 0.]\n",
"):\n",
" x = np.arange(nx)*dx + x0[0]\n",
" y = np.arange(ny)*dy + x0[1]\n",
" z = np.r_[0.]\n",
" return Utils.ndgrid(x, y, z)\n",
"xyz_topo = get_xyz()\n",
"xyz = get_xyz(x0=np.r_[0, 0, 30.])"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
Expand All @@ -19,57 +49,32 @@
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-5-b93fa3ba971d>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m 63\u001b[0m \u001b[0ma\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0ma_global\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 64\u001b[0m \u001b[0minput_currents\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0minput_currents_global\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 65\u001b[0;31m \u001b[0mtime_input_currents\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mtime_input_currents_global\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 66\u001b[0m )\n\u001b[1;32m 67\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/anaconda2/envs/simpeg/lib/python3.6/site-packages/properties/base/base.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(cls, *args, **kwargs)\u001b[0m\n\u001b[1;32m 250\u001b[0m \u001b[0;32mwith\u001b[0m \u001b[0mhandlers\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlisteners_disabled\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 251\u001b[0m \u001b[0mobj\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_reset\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 252\u001b[0;31m \u001b[0mobj\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__init__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 253\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mobj\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 254\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m<ipython-input-6-8835dccb0d03>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m 38\u001b[0m \u001b[0ma\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0ma_global\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 39\u001b[0m \u001b[0minput_currents\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0minput_currents_global\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 40\u001b[0;31m \u001b[0mtime_input_currents\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mtime_input_currents_global\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 41\u001b[0m )\n\u001b[1;32m 42\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/anaconda2/envs/simpeg/lib/python3.6/site-packages/properties/base/base.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(cls, *args, **kwargs)\u001b[0m\n\u001b[1;32m 251\u001b[0m \u001b[0;32mwith\u001b[0m \u001b[0mhandlers\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlisteners_disabled\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 252\u001b[0m \u001b[0mobj\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_reset\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 253\u001b[0;31m \u001b[0mobj\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__init__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 254\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mobj\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 255\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/Projects/simpegEM1D/simpegEM1D/GlobalEM1D.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, **kwargs)\u001b[0m\n\u001b[1;32m 517\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m__init__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 518\u001b[0m \u001b[0mGlobalEM1DSurvey\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__init__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 519\u001b[0;31m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mset_parameters\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 520\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 521\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mset_parameters\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/Projects/simpegEM1D/simpegEM1D/GlobalEM1D.py\u001b[0m in \u001b[0;36mset_parameters\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 554\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 555\u001b[0m \u001b[0;31m# List\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 556\u001b[0;31m \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtime_input_currents\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 557\u001b[0m self.time_input_currents = [\n\u001b[1;32m 558\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mempty\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mfloat\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mn_sounding\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mValueError\u001b[0m: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()"
]
}
],
"source": [
"from SimPEG import Mesh, Maps\n",
"import numpy as np\n",
"from matplotlib.colors import LogNorm\n",
"from simpegEM1D import (\n",
" GlobalEM1DProblemTD, GlobalEM1DSurveyTD, get_vertical_discretization_time, EM1DSurveyTD\n",
")\n",
"from pymatsolver import PardisoSolver\n",
"from simpegEM1D import skytem_HM_2015\n",
"wave = skytem_HM_2015()\n",
"time = np.logspace(-5, -2, 21)\n",
"hz = get_vertical_discretization_time(time, facter_tmax=0.5, factor_tmin=10.)\n",
"\n",
"time_input_currents = wave.current_times[-7:]\n",
"input_currents = wave.currents[-7:]\n",
"\n",
"from scipy.spatial import Delaunay\n",
"def PolygonInd(mesh, pts):\n",
" hull = Delaunay(pts)\n",
" inds = hull.find_simplex(mesh.gridCC)>=0\n",
" return inds\n",
"\n",
"n_sounding = 5\n",
"dx = 20.\n",
"hx = np.ones(n_sounding) * dx\n",
"n_sounding = xyz.shape[0]\n",
"hz = get_vertical_discretization_time(time, facter_tmax=0.5, factor_tmin=10.)\n",
"hx = np.ones(n_sounding)\n",
"mesh = Mesh.TensorMesh([hx, hz], x0='00')\n",
"inds = np.logical_and(mesh.gridCC[:,1]>25, mesh.gridCC[:,1]<75.) \n",
"sigma = np.ones(mesh.nC) * 1./20.\n",
"sigma[inds] = 1./100.\n",
"x0 = np.r_[0., 75.]\n",
"x1 = np.r_[dx*n_sounding, 75.]\n",
"x2 = np.r_[dx*n_sounding, 90.]\n",
"x3 = np.r_[0., 200.]\n",
"pts = np.vstack((x0, x1, x2, x3, x0))\n",
"poly_inds = PolygonInd(mesh, pts)\n",
"sigma[poly_inds] = 0.1\n",
"sigma_em1d = sigma.reshape(mesh.vnC, order='F').flatten()\n",
"\n",
"x = mesh.vectorCCx\n",
"y = np.zeros_like(x)\n",
"z = np.ones_like(x) * 30.\n",
"rx_locations = np.c_[x, y, z]\n",
"src_locations = np.c_[x, y, z]\n",
"topo = np.c_[x, y, z-30.].astype(float)\n",
"rx_locations = xyz.copy()\n",
"src_locations = xyz.copy()\n",
"topo = xyz_topo.copy()\n",
"mapping = Maps.ExpMap(mesh)\n",
"\n",
"rx_type_global = np.array([\"dBzdt\"], dtype=str).repeat(n_sounding, axis=0)\n",
Expand All @@ -96,12 +101,93 @@
")\n",
"\n",
"prob = GlobalEM1DProblemTD(\n",
" [], sigmaMap=mapping, hz=hz, parallel=True,\n",
" [], sigmaMap=mapping, hz=hz, parallel=True, ncpu=2,\n",
" Solver=PardisoSolver\n",
")\n",
"prob.pair(survey)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"> \u001b[0;32m/Users/sgkang/Projects/simpegEM1D/simpegEM1D/GlobalEM1D.py\u001b[0m(556)\u001b[0;36mset_parameters\u001b[0;34m()\u001b[0m\n",
"\u001b[0;32m 554 \u001b[0;31m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0m\u001b[0;32m 555 \u001b[0;31m \u001b[0;31m# List\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0m\u001b[0;32m--> 556 \u001b[0;31m \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtime_input_currents\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0m\u001b[0;32m 557 \u001b[0;31m self.time_input_currents = [\n",
"\u001b[0m\u001b[0;32m 558 \u001b[0;31m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mempty\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mfloat\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mn_sounding\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0m\n",
"ipdb> self.time_input_currents\n",
"array([[-0.004 , -0.00391 , -0.0038 , 0. , 0.00019437,\n",
" 0.00019504, 0.00019637],\n",
" [-0.004 , -0.00391 , -0.0038 , 0. , 0.00019437,\n",
" 0.00019504, 0.00019637],\n",
" [-0.004 , -0.00391 , -0.0038 , 0. , 0.00019437,\n",
" 0.00019504, 0.00019637],\n",
" [-0.004 , -0.00391 , -0.0038 , 0. , 0.00019437,\n",
" 0.00019504, 0.00019637],\n",
" [-0.004 , -0.00391 , -0.0038 , 0. , 0.00019437,\n",
" 0.00019504, 0.00019637],\n",
" [-0.004 , -0.00391 , -0.0038 , 0. , 0.00019437,\n",
" 0.00019504, 0.00019637],\n",
" [-0.004 , -0.00391 , -0.0038 , 0. , 0.00019437,\n",
" 0.00019504, 0.00019637],\n",
" [-0.004 , -0.00391 , -0.0038 , 0. , 0.00019437,\n",
" 0.00019504, 0.00019637],\n",
" [-0.004 , -0.00391 , -0.0038 , 0. , 0.00019437,\n",
" 0.00019504, 0.00019637],\n",
" [-0.004 , -0.00391 , -0.0038 , 0. , 0.00019437,\n",
" 0.00019504, 0.00019637],\n",
" [-0.004 , -0.00391 , -0.0038 , 0. , 0.00019437,\n",
" 0.00019504, 0.00019637],\n",
" [-0.004 , -0.00391 , -0.0038 , 0. , 0.00019437,\n",
" 0.00019504, 0.00019637],\n",
" [-0.004 , -0.00391 , -0.0038 , 0. , 0.00019437,\n",
" 0.00019504, 0.00019637],\n",
" [-0.004 , -0.00391 , -0.0038 , 0. , 0.00019437,\n",
" 0.00019504, 0.00019637],\n",
" [-0.004 , -0.00391 , -0.0038 , 0. , 0.00019437,\n",
" 0.00019504, 0.00019637],\n",
" [-0.004 , -0.00391 , -0.0038 , 0. , 0.00019437,\n",
" 0.00019504, 0.00019637],\n",
" [-0.004 , -0.00391 , -0.0038 , 0. , 0.00019437,\n",
" 0.00019504, 0.00019637],\n",
" [-0.004 , -0.00391 , -0.0038 , 0. , 0.00019437,\n",
" 0.00019504, 0.00019637],\n",
" [-0.004 , -0.00391 , -0.0038 , 0. , 0.00019437,\n",
" 0.00019504, 0.00019637],\n",
" [-0.004 , -0.00391 , -0.0038 , 0. , 0.00019437,\n",
" 0.00019504, 0.00019637],\n",
" [-0.004 , -0.00391 , -0.0038 , 0. , 0.00019437,\n",
" 0.00019504, 0.00019637],\n",
" [-0.004 , -0.00391 , -0.0038 , 0. , 0.00019437,\n",
" 0.00019504, 0.00019637],\n",
" [-0.004 , -0.00391 , -0.0038 , 0. , 0.00019437,\n",
" 0.00019504, 0.00019637],\n",
" [-0.004 , -0.00391 , -0.0038 , 0. , 0.00019437,\n",
" 0.00019504, 0.00019637],\n",
" [-0.004 , -0.00391 , -0.0038 , 0. , 0.00019437,\n",
" 0.00019504, 0.00019637]])\n",
"ipdb> if not self.time_input_currents\n",
"*** SyntaxError: invalid syntax\n",
"ipdb> if not self.time_input_currents\n",
"*** SyntaxError: invalid syntax\n",
"ipdb> if not self.time_input_currents: \n",
"*** SyntaxError: unexpected EOF while parsing\n",
"ipdb> quit\n"
]
}
],
"source": [
"debug"
]
},
{
"cell_type": "markdown",
"metadata": {},
Expand Down Expand Up @@ -512,7 +598,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.3"
"version": "3.6.4"
}
},
"nbformat": 4,
Expand Down

0 comments on commit 2d5d01e

Please sign in to comment.