{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from MT1D import MT1DProblem, MT1DSurvey, MT1DSrc, ZxyRx, Survey, AppResPhaRx\n", "from SimPEG import Maps, DataMisfit\n", "from scipy.constants import mu_0\n", "import numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Put all things together to SimPEG Problem\n", "\n", "In the [1_MT1D_NumericalSetup](./1_MT1D_NumericalSetup.ipynb) notebook, we coded up a 1D MT simulation function that computes complex impedances at multiple frequencies for a given conductivity model. To invert MT data in 1D using gradient-based optimization, sensitivity functions were derived, coded up, and tested in [Appendix_A_MT1D_Sensitivity](./Appendix_A_MT1D_Sensitivity.ipynb) notebook. Basically those are three funcionalities:\n", "\n", "- `dpred(m)`: takes model and predict MT data for a given model, `m`\n", "\n", "- `Jvec(m,v)`: takes model and a vector `v`, then computes sensitivity-vector product. \n", "\n", "- `Jtvec(m,v)`: takes model and a vector `v`, then computes adjoint sensitivity-vector product. \n", "\n", "Now we put all these things to `SimPEG`'s framework. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## `SimPEG` framework\n", "\n", "Our setup of the simluation follows the [SimPEG framework](http://simpeg.xyz). \n", "\n", "\"plane_wave\"\n", "\n", "\n", "For more details, see the [SimPEG docs](http://docs.simpeg.xyz)\n", "\n", "- `Mesh`: mesh geometry and differential operators\n", "- `Problem`: physics engine. contains the machinery to construct and solve the PDE\n", "- `Survey`: sources and receivers\n", "- `Src`: sources\n", "- `Rx`: receivers\n", "\n", "`Problem` class in `SimPEG` is a physics. It discretizes the earth and computes fields in a discretized domain, and here we use `Mesh` class. For an MT problem, Maxwell's equations are used, and electromagnetic fields are computed. To do this, we need to know details about MT survey (e.g. frequencies and location of source and receivers), and they are defined in `Survey` class. `Survey` class takes `Src` class, and `Src` class takes `Rx` class. Therefore, both information in sources and receivers can be aware in the `Survey`.\n", "\n", "### Mapping\n", "\n", "For simulation, we do need to input discretized physical property. For instance in 1D MT problem values of conductivity at each cell is required. However, when consdiering inversion, inversion model can be different from physical property (e.g. conductivity). For instance, often in electromagnetic (EM) inversions logarithmic conductivity was used as an inversion model: $\\mathbf{m} = log(\\boldsymbol{\\sigma})$. While inversion process, still simulation is required to predict data hence updated model, $\\mathbf{m}$ has to be converted to $\\boldsymbol{\\sigma}$, which is a discretized physical property. Here we define a mapping:\n", "\n", "$$ \\sigma = \\mathcal{M}(\\mathbf{m}) $$\n", "\n", "mapping is a function which takes an inversion model, and outputs physical property used in simulation. For more details see: http://docs.simpeg.xyz/content/api_core/api_Maps.html. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Codes:\n", "\n", "Taking snippets of codes written in [1_MT1D_NumericalSetup](./1_MT1D_NumericalSetup.ipynb) and [Appendix_A_MT1D_Sensitivity](./Appendix_A_MT1D_Sensitivity.ipynb), we generated the `Problem`, `Survey`, `Src`, and `Rx` classes for 1D MT problem. And here are lists of classes:\n", "\n", "- `MT1DProblem`: This includes functions: `dpred(m)`, `Jvec(m,v)`, `Jtvec(m,v)`\n", "- `MT1DSurvey`: Need to input `MT1DSrc`\n", "- `MT1DSrc`: Need to input `ZxyRx`\n", "- `ZxyRx`: Receiver class / Real and Imaginary part of $Z_{xy}$" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ ">> Smallest cell size = 50 m\n", ">> Padding distance = 316227 m\n", ">> # of padding cells 17\n", ">> # of core cells cells 16\n" ] } ], "source": [ "rxloc = np.r_[0.]\n", "srcloc = np.r_[0.]\n", "frequency = np.logspace(-3, 2, 25)\n", "\n", "# Rx\n", "rx = ZxyRx(rxloc, component=\"both\", frequency=frequency)\n", "rxList = [rx]\n", "\n", "# Src\n", "src = MT1DSrc(rxList, loc=srcloc)\n", "\n", "# Survey\n", "survey = MT1DSurvey([src])\n", "mesh = survey.setMesh(\n", " sigma=0.01, max_depth_core=5000., \n", " ncell_per_skind=10, n_skind=2, \n", " core_meshType=\"log\", max_hz_core=1000.\n", ")\n", "\n", "# Physical property: Conductivity\n", "sigma = np.ones(mesh.nC) * 0.01\n", "\n", "# Problem\n", "prob = MT1DProblem(mesh, sigmaMap=Maps.ExpMap(mesh), verbose=False)\n", "prob.pair(survey)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Order and Adjoint tests \n", "\n", "We perform both order and adjoint tests shown in [Appendix_A_MT1D_Sensitivity](./Appendix_A_MT1D_Sensitivity.ipynb) using 1D MT problem in `SimPEG`'s framework. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Order test: `Jvec`" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "==================== checkDerivative ====================\n", "iter h |ft-f0| |ft-f0-h*J0*dx| Order\n", "---------------------------------------------------------\n", " 0 1.00e-01 2.669e-01 5.676e-02 nan\n", " 1 1.00e-02 2.151e-02 4.914e-04 2.063\n", " 2 1.00e-03 2.106e-03 4.846e-06 2.006\n", "========================= PASS! =========================\n", "You are awesome.\n", "\n" ] }, { "data": { "text/plain": [ "True" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from SimPEG import Tests\n", "def derChk(m):\n", " return [survey.dpred(m), lambda mx: prob.Jvec(m, mx)]\n", "Tests.checkDerivative(derChk, np.log(sigma), plotIt=False, num=3, eps=1e-20, dx=np.log(sigma)*2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Order test: `Jtvec`" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "SimPEG.DataMisfit.l2_DataMisfit assigning default eps of 1e-5 * ||dobs||\n" ] } ], "source": [ "# Evaluate predicted data\n", "pred = survey.dpred(np.log(sigma))\n", "survey.dobs = pred\n", "dmis = DataMisfit.l2_DataMisfit(survey)\n", "m0 = np.log(sigma)*2." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "==================== checkDerivative ====================\n", "iter h |ft-f0| |ft-f0-h*J0*dx| Order\n", "---------------------------------------------------------\n", " 0 1.00e-01 1.074e+06 4.010e+05 nan\n", " 1 1.00e-02 7.038e+04 3.059e+03 2.118\n", " 2 1.00e-03 6.762e+03 2.982e+01 2.011\n", " 3 1.00e-04 6.735e+02 2.975e-01 2.001\n", "========================= PASS! =========================\n", "You are awesome.\n", "\n" ] }, { "data": { "text/plain": [ "True" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Tests.checkDerivative(\n", " lambda m: [dmis(m), dmis.deriv(m)],\n", " m0,\n", " plotIt=False,\n", " num=4,\n", " dx = m0*1.\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Adjoint test" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Adjoint Test 8.881784197001252e-16 True\n" ] } ], "source": [ "v = np.random.rand(mesh.nC)\n", "w = np.random.rand(pred.shape[0])\n", "wtJv = w.dot(prob.Jvec(m0, v))\n", "vtJtw = v.dot(prob.Jtvec(m0, w))\n", "passed = np.abs(wtJv - vtJtw) < 1e-10\n", "print('Adjoint Test', np.abs(wtJv - vtJtw), passed)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Ready to run an inversion\n", "\n", "Once you have created the `Survey` and `Problem` classes so that you can compute predicted data (`dpred`), forward and adjoint sensitivity-vector products (`Jvec` and `Jtvec`), you are all set to run inversion using `SimPEG`'s inverersion framework. Check out the examples in [3_MT1D_5layer_inversion](./3_MT1D_5layer_inversion.ipynb) and [1_MT1D_NumericalSetup](./1_MT1D_NumericalSetup.ipynb).\n" ] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.6" } }, "nbformat": 4, "nbformat_minor": 2 }