{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import scipy.sparse as sp\n", "import matplotlib.pyplot as plt\n", "from SimPEG import Mesh, Utils, Solver\n", "from scipy.constants import mu_0, epsilon_0\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Sensitivity computuation for 1D magnetotelluric (MT) problem" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Purpose\n", "\n", "With [SimPEG's](http://simpeg.xyz) mesh class, we discretize the sensitivity function for a 1D magnetotelluric problem. Rather than generating the full sensitivity matrix, we compute the forward (`Jvec`) and adjoint (`Jtvec`) functionals that can evalaute matrix-vector product. There are some milestones to be accomplished:\n", "\n", "- Break apart senstivity function, $J_{\\sigma}$ into to parts then derive each of them (using chain rules):\n", "\n", "$$ \n", "J_{\\sigma} = \\frac{d P(u)}{d \\sigma} = \\frac{d P(u)}{d u} \\frac{d u}{d \\sigma}\n", "$$\n", "\n", "- Compute forward and adjoint sensitivity function: `Jvec` and `Jtvec`\n", "\n", "- Test `Jvec` and `Jtvec`: Order test and Adjoint test" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Discretzation (forward simulation)\n", "\n", "We define physical properties at cell centers, and stagger the electric and magnetic fields\n", "\n", "- $\\sigma$, $\\mu$, $\\epsilon$ : cell centers\n", "- $E_x$: cell centers\n", "- $H_y$: faces\n", "\n", " \n", "\n", "and use a finite difference approach to define the operators, this gives us the discrete system of equations\n", "\n", "$$\n", "\\underbrace{\n", " \\begin{bmatrix}\n", " \\mathbf{Grad} & \\imath \\omega \\mathbf{M}^{f}_{\\mu} \\\\[0.3em]\n", " \\mathbf{M}^{cc}_{\\boldsymbol{\\sigma}} & \\mathbf{Div} \\\\[0.3em]\n", " \\end{bmatrix}\n", "}_{\\mathbf{A}}\n", "\\underbrace{\n", " \\begin{bmatrix}\n", " \\mathbf{e_x} \\\\[0.3em]\n", " \\mathbf{h_y} \\\\[0.3em]\n", " \\end{bmatrix}\n", "}_{\\mathbf{u}}\n", "=\n", "\\underbrace{\n", " \\begin{bmatrix}\n", " - \\mathbf{B}\\mathbf{e_x}^{BC} \\\\[0.3em]\n", " \\boldsymbol{0} \\\\[0.3em]\n", " \\end{bmatrix}\n", "}_{\\mathbf{rhs}}\n", "$$\n", "\n", "with \n", "\n", "- $\\mathbf{e_x}$: Discrete $E_x$, on cell centers $[\\text{nC} \\times 1]$\n", "\n", "- $\\mathbf{h_y}$: Dicrete $H_x$, on cell faces $[(\\text{nC}+1) \\times 1]$\n", "\n", "- $ \\mathbf{Grad}$: Discrete gradient operator $[\\text{nC} \\times (\\text{nC}+1)]$\n", "\n", "- $ \\mathbf{Div}$: Discrete divergence operator $[(\\text{nC}+1) \\times \\text{nC}]$\n", "\n", "- $\\mathbf{M}^{f}_{\\boldsymbol{\\mu}} = \\mathbf{diag}(\\mathbf{Av^{cc2f}} \\boldsymbol{\\mu})$ $[(\\text{nC}+1) \\times (\\text{nC}+1)]$\n", "\n", "- $\\mathbf{M}^{cc}_{\\boldsymbol{\\boldsymbol{\\sigma}}} = \\mathbf{diag}(\\boldsymbol{{\\sigma}})$ $[\\text{nC} \\times \\text{nC}]$. Here we are using the quasi-static assumption for brevity. \n", "\n", "- $\\mathbf{B} \\mathbf{e_x}^{BC}$ handles the boundary conditions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## What are the data?\n", "\n", "The measured data in general can be defined as:\n", "\n", "$$ \\mathbf{d} = P(\\mathbf{u}) $$\n", "\n", "where $P(\\cdot)$ is a evaluation functional which takes a solution vector $\\mathbf{u}$ and ouputs data at a receiver location.\n", "\n", "Here, we use impedace data (one could also consider using apparent resistivity and phase). The impedance is complex, so we treat the real and imaginary components of each as two separate data points\n", "\n", "$$\n", "Z_{xy} = -\\frac{E_x}{H_y} = \\text{Re}[Z_{xy}] + i ~\\text{Im}[Z_{xy}]\n", "$$\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The impedance $Z_{xy}$ can be evaluated from the solution vector $\\mathbf{u}$. We will evaluate data at $z=0m$. The solution vector we obtain from the forward simulation is: \n", "\n", "$$\n", "\\mathbf{u} = \n", " \\begin{bmatrix}\n", " \\mathbf{e_x} \\\\[0.3em]\n", " \\mathbf{h_y} \\\\[0.3em]\n", " \\end{bmatrix}\n", "$$\n", "\n", "At the surface, we specified the boundary condition that $E_x(z=0) = 1$. So what we need $P(\\dot)$ to accomplish is \n", "$$\n", "Z_{xy}\\big|_{z=0} = -\\frac{1}{h_y(z=0)}\n", "$$\n", "\n", "Thus, $P(\\dot)$ can be defined as an interpolation matrix that simply extracts the value of $h_y$ at the surface. We denote this matrix: $\\mathbf{P}_{0}$ (Thinking in terms of matrices is very helpful when we get to the step of taking derivatives!)\n", "\n", "$$\\mathbf{d} = Z_{xy} = - \\mathbf{P}_{0}\\left(\\frac{1}{\\mathbf{u}}\\right) $$\n", "\n", "From complex-values $Z_{xy}$, we can compute real and imagniary part of the $Z_{xy}$ then the data can be defined as:\n", "\n", "$$\n", "\\mathbf{d} = \\begin{bmatrix}\n", " \\text{Re}[Z_{xy}] \\\\[0.3em]\n", " \\text{Im}[Z_{xy}] \\\\[0.3em]\n", "\\end{bmatrix}\n", "$$\n", "\n", "We will set up an example and go through the steps to compute a datum. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Set up Mesh and Model" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "def skin_depth(sigma, f):\n", " \"\"\"\n", " Depth at which the fields propagating through a homogeneous medium \n", " have decayed by a factor of 1/e for a given frequency, f and conductivity, sigma\n", " \"\"\"\n", " return 500./np.sqrt(sigma * f)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The minimum skin depth is 500.00m\n", "The maximum skin depth is 1.58e+05m\n" ] } ], "source": [ "rho_half = 100. # Resistivity of the halfspace in Ohm-m\n", "sigma_half = 1./rho_half # Conductivity is the inverse of conductivity\n", "frequency = np.logspace(-3, 2, 25) # frequencies at which to simulate the MT problem\n", "\n", "skin_depth_min = skin_depth(sigma_half, frequency.max())\n", "skin_depth_max = skin_depth(sigma_half, frequency.min())\n", "\n", "print(\"The minimum skin depth is {:1.2f}m\".format(skin_depth_min))\n", "print(\"The maximum skin depth is {:1.2e}m\".format(skin_depth_max))" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The smallest cell size is 125.00m\n", "The core region of the mesh extends 5.00e+03m\n", "The mesh should extend at least 3.16e+05m\n" ] } ], "source": [ "cs = skin_depth_min / 4.\n", "core_extent = 5000. \n", "domain_extent = 2 * skin_depth_max\n", "\n", "print(\"The smallest cell size is {:1.2f}m\".format(cs))\n", "print(\"The core region of the mesh extends {:1.2e}m\".format(core_extent))\n", "print(\"The mesh should extend at least {:1.2e}m\".format(domain_extent))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Add padding to extend sufficiently far" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "25 padding cells extends 3.82e+05m > 3.16e+05m (2 skin depths)\n" ] } ], "source": [ "npad = 1 # start with 1 cell\n", "padding_fact = 1.3 # the amount by which we will expand each cell of the padding\n", "\n", "def padding_extent(npad):\n", " \"\"\"\n", " given a number of padding cells, this computes how far the padding extends\n", " \"\"\"\n", " padding_widths = cs*padding_fact**(np.arange(npad) + 1)\n", " return padding_widths.sum()\n", "\n", "# keep adding padding until we are beyond the desired extent\n", "padding_z = padding_extent(npad)\n", "while padding_z < domain_extent:\n", " npad+=1\n", " padding_z = padding_extent(npad)\n", " \n", "print(\"{:1.0f} padding cells extends {:1.2e}m > {:1.2e}m (2 skin depths)\".format(\n", " npad, padding_extent(npad), domain_extent\n", "))" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "There are 40 cells in the mesh. The mest extends 3.87e+05m\n" ] } ], "source": [ "ncz = np.ceil(core_extent / cs) # number of cells in the core domain\n", "hz = [(cs, npad, -1.3), (cs, ncz)] # define how to construct the cell widths\n", "mesh = Mesh.TensorMesh([hz], x0='N') # construct a 1D Tensor Mesh\n", "\n", "print(\"There are {:1.0f} cells in the mesh. The mest extends {:1.2e}m\".format(\n", " ncz, mesh.hx.sum()\n", ")) \n", "\n", "sigma = np.ones(mesh.nC) * sigma_half" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Forward simulation function\n", "\n", "Forward simulation function `dpred` takes conductivity model (`nCx1`), and frequency (`1x1`), and outputs real and imaginary part of the impedance, $Z_{xy}$ (`2x1`). By solving $\\mathbf{A} \\mathbf{u}=\\mathbf{rhs}$, we compute solution vector $\\mathbf{u}$, then evaluate real and imaginary part of impedance at $z=0$ ($\\mathbf{d} = P(\\mathbf{u})$). " ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The projection matrix has shape (1, 131) with 1 non-zero entry at (0,130)\n" ] } ], "source": [ "# Projection Matrix\n", "P0 = sp.csr_matrix(\n", " (np.r_[1.], (np.r_[0], np.r_[mesh.nF+mesh.nC-1])), shape=(1, mesh.nF+mesh.nC)\n", " )\n", "\n", "print(\n", " \"The projection matrix has shape {} with {} non-zero entry at ({},{})\".format(\n", " P0.shape, P0.nnz, P0.nonzero()[0][0], P0.nonzero()[1][0]\n", " )\n", ")\n" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "def dpred(sigma, f=100.):\n", " \n", " # angular frequency\n", " omega = f*2.*np.pi\n", " \n", " # physical properties\n", " mu = np.ones(mesh.nC)*mu_0 # magnetic permeability values for all cells\n", " sigmahat = sigma.copy() # here we are ignoring displacement current\n", " \n", " # Grad \n", " mesh.setCellGradBC([['dirichlet', 'dirichlet']]) # Setup boundary conditions\n", " Grad = mesh.cellGrad # Gradient matrix\n", "\n", " # MfMu\n", " Mmu = Utils.sdiag(mesh.aveCC2F * mu) \n", " \n", " # Mccsigma\n", " Msighat = Utils.sdiag(sigmahat) \n", " \n", " # Div\n", " Div = mesh.faceDiv # Divergence matrix\n", " \n", " # Right Hand Side\n", " B = mesh.cellGradBC # a matrix for boundary conditions\n", " Exbc = np.r_[0., 1.] # boundary values for Ex\n", "\n", " # A-matrix\n", " A = sp.vstack([\n", " sp.hstack([Grad, 1j*omega*Mmu]), # Top row of A matrix\n", " sp.hstack((Msighat, Div)) # Bottom row of A matrix\n", " ])\n", " \n", " # Right-hand side\n", " rhs = np.r_[\n", " -B*Exbc, \n", " np.zeros(mesh.nC)\n", " ] \n", "\n", " Ainv = Solver(A) # Factor A matrix\n", " u = Ainv*rhs # Solve A^-1 rhs = u\n", "\n", " # build the projection matrix, a sparse matrix that has 1 entry to grab the value of h_y at the surface\n", " P0 = sp.csr_matrix(\n", " (np.r_[1.], (np.r_[0], np.r_[len(u)-1])), shape=(1, len(u))\n", " )\n", " \n", " Zxy = - P0 * (1./(u))\n", " \n", " return np.r_[Zxy.real, Zxy.imag]" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "At f=100.0Hz, we have two data. \n", "Re[Z] = 0.196, Im[Z] = 0.202\n" ] } ], "source": [ "f = 100.\n", "data = dpred(sigma, f=f)\n", "\n", "print(\"At f={}Hz, we have two data. \\nRe[Z] = {:1.3f}, Im[Z] = {:1.3f}\".format(f, data[0], data[1]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Sensitivity of datum with regard to $\\sigma$:\n", "\n", "The sensitivity function shows how much the data are changed due to changes in model paramters. Understanding how \"sensitive\" our data are to the model is important for survey design. It is essential that we be able to compute the sensitivity when using gradient-based optimization techniques, as the sensitivity gives us the direction in which to take a step and update the model as we search for a solution. \n", "\n", "The sensitivity function can be defined as \n", "\n", "$$ J_{\\sigma} = \\frac{d P(u)}{d \\sigma}$$\n", "\n", "The size of the sensitivity is $[nD \\times n\\sigma]$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To obtain above sensitivity function in discrete space, we first differentiate \n", "\n", "$$\\mathbf{A}\\mathbf{u} = \\mathbf{rhs}$$\n", "\n", "w.r.t ${\\boldsymbol{\\sigma}}$, which yields \n", "\n", "$$ \\frac{d \\mathbf{A}}{d {\\boldsymbol{\\sigma}}}\\mathbf{u} + \\mathbf{A} \\frac{d \\mathbf{u} }{d {\\boldsymbol{\\sigma}}}= 0 $$\n", "\n", "Rearranging and multiplyting by $\\mathbf{A}^{-1}$ on both sides yields\n", "\n", "$$ \\frac{d \\mathbf{u} }{d {\\boldsymbol{\\sigma}}}= -\\mathbf{A}^{-1}\\frac{d \\mathbf{A}}{d {\\boldsymbol{\\sigma}}}\\mathbf{u} $$\n", "\n", "Next, we need to include the evaluation, $\\mathbf{d} = P(\\mathbf{u})$ which is the operation taken on $\\mathbf{u}$ to give us the data. From this, we obtain\n", "\n", "$$ \\mathbf{J}_{{\\boldsymbol{\\sigma}}} = \n", "\\frac{\\partial P(\\mathbf{u})}{\\partial {\\mathbf{u}}}\\Big(\\frac{d \\mathbf{u} }{d {\\boldsymbol{\\sigma}}}\\Big) = \n", "-\\frac{\\partial P(\\mathbf{u})}{\\partial {\\mathbf{u}}} \\Big(\\mathbf{A}^{-1}\\frac{d \\mathbf{A}}{d {\\boldsymbol{\\sigma}}}\\mathbf{u}\\Big) $$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From this, there are two derivatives need to be computed:\n", "\n", "1. $$\\frac{d \\mathbf{A}}{d \\boldsymbol{\\sigma}}\\mathbf{u}=?$$\n", "\n", "2. $$\\frac{\\partial P(\\mathbf{u})}{\\partial \\mathbf{u}}=?$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### First part of the sensitivity, $\\frac{d \\mathbf{A}}{d \\boldsymbol{\\sigma}}\\mathbf{u}=?$\n", "For $\\frac{d \\mathbf{A}}{d \\boldsymbol{\\sigma}}\\mathbf{u}$, keep in mind that we are treating $\\mathbf{u}$ as fixed and that \n", "\n", "$$\n", "\\mathbf{A}\\mathbf{u} = \n", " \\begin{bmatrix}\n", " \\mathbf{Grad} & \\imath \\omega \\mathbf{M}^{f}_{\\mu} \\\\[0.3em]\n", " \\mathbf{M}^{cc}_{\\boldsymbol{\\sigma}} & \\mathbf{Div} \\\\[0.3em]\n", " \\end{bmatrix} \n", " \\begin{bmatrix}\n", " \\mathbf{e_x} \\\\[0.3em]\n", " \\mathbf{h_y} \\\\[0.3em]\n", " \\end{bmatrix}\n", "$$\n", "\n", "Here, we see that the only dependence on $\\boldsymbol{\\sigma}$ is in the matrix $\\mathbf{M}^{cc}_{\\hat{\\boldsymbol{\\sigma}}} = \\mathbf{diag}(\\boldsymbol{\\sigma})$. So lets focus our attention on that block. We are taking the derivative of a Matrix-vector product (which is just a vector) with respect to the vector $\\boldsymbol{\\sigma}$, so the result should be a matrix. We write out the problem: \n", "\n", "$$\n", "\\frac{\\partial}{\\partial \\boldsymbol{\\sigma}} \\mathbf{M}^{cc}_{\\hat{\\boldsymbol{\\sigma}}} \\mathbf{e_x}^{fix} = \\frac{\\partial}{\\partial \\boldsymbol{\\sigma}} \\mathbf{diag}(\\boldsymbol{\\sigma}) \\mathbf{e_x}^{fix}\n", "$$\n", "\n", "and since $\\text{diag}(\\mathbf{x})\\mathbf{y} = \\mathbf{diag}(\\mathbf{y})\\mathbf{x}$, we can interchange the roles of $\\boldsymbol{\\sigma}$ and $\\mathbf{e_x}^{fix}$\n", "\n", "$$\n", "\\frac{\\partial}{\\partial \\boldsymbol{\\sigma}} \\mathbf{M}^{cc}_{\\hat{\\boldsymbol{\\sigma}}} \\mathbf{e_x}^{fix} = \\frac{\\partial}{\\partial \\boldsymbol{\\sigma}} \\mathbf{diag}(\\mathbf{e_x}^{fix}) \\boldsymbol{\\sigma}\n", "$$\n", "\n", "So the derivative is simply: \n", "$$\n", "\\frac{\\partial}{\\partial \\boldsymbol{\\sigma}} \\mathbf{M}^{cc}_{\\hat{\\boldsymbol{\\sigma}}} \\mathbf{e_x}^{fix} =\\text{diag}(\\mathbf{e_x}^{fix})\n", "$$\n", "\n", "Thus the full derivative is \n", "$$ \n", "\\frac{d \\mathbf{A}}{d \\boldsymbol{\\sigma}}\\mathbf{u} =\n", "\\begin{bmatrix}\n", " \\mathbf{0} \\\\[0.3em]\n", " \\mathbf{diag}(\\mathbf{e}_x) \\\\[0.3em]\n", "\\end{bmatrix}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Second part of the Sensitivity: $\\frac{\\partial P(\\mathbf{u})}{\\partial \\mathbf{u}}=?$\n", "For the other one we consider when the data is defined as real and imaginary parts of $Z_{xy} = -\\mathbf{P}_0\\frac{1}{\\mathbf{u}}$:\n", "\n", "Taking derivative of $Z_{xy}$ w.r.t. $\\mathbf{u}$ yields\n", "$$\n", "\\frac{\\partial Z_{\n", "xy}}{\\partial \\mathbf{u}}=\n", "\\mathbf{P}_0\\frac{1}{\\mathbf{u}^2}\n", "$$\n", "\n", "$$ \n", "\\frac{\\partial P(\\mathbf{u})}{\\partial \\mathbf{u}}=\n", "\\begin{bmatrix}\n", " \\frac{\\partial Re[Z_{xy}]}{\\partial \\mathbf{u}} \\\\[0.3em]\n", " \\frac{\\partial Im[Z_{xy}]}{\\partial \\mathbf{u}} \\\\[0.3em]\n", "\\end{bmatrix} \n", "=\n", "\\begin{bmatrix}\n", " Re[\\mathbf{P}_0\\frac{1}{\\mathbf{u}^2}] \\\\[0.3em]\n", " Im[\\mathbf{P}_0\\frac{1}{\\mathbf{u}^2}] \\\\[0.3em]\n", "\\end{bmatrix} \n", "$$\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can form sensitivity matrix $\\mathbf{J}_{\\sigma}$ by combining above equations:\n", "\n", "$$ \\mathbf{J}_{{\\boldsymbol{\\sigma}}} = \n", "-\\frac{\\partial P(\\mathbf{u})}{\\partial {\\mathbf{u}}} \\Big(\\mathbf{A}^{-1}\\frac{d \\mathbf{A}}{d {\\boldsymbol{\\sigma}}}\\mathbf{u}\\Big) $$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Deriving sensitivity for apparent resistivity and phase is possible, but this requires additional details hence, we focus on our attention to real and imaginary parts of impedance. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compute Sensitivity function\n", "\n", "We compute discretized sensitivity matrix, $\\mathbf{J}_{\\sigma}$ shown above. " ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "# angular frequency\n", "omega = f*2*np.pi\n", "\n", "# physical properties\n", "sigma = np.ones(mesh.nC)*sigma_half # conductivity values for all cells\n", "mu = np.ones(mesh.nC)*mu_0 # magnetic permeability values for all cells\n", "sigmahat = sigma.copy()\n", "\n", "# Grad \n", "mesh.setCellGradBC([['dirichlet', 'dirichlet']]) # Setup boundary conditions\n", "Grad = mesh.cellGrad # Gradient matrix\n", "\n", "# MfMu\n", "Mmu = Utils.sdiag(mesh.aveCC2F * mu)\n", "\n", "# Mccsigma\n", "Msighat = Utils.sdiag(sigmahat) \n", "\n", "# Div\n", "Div = mesh.faceDiv # Divergence matrix\n", "\n", "# Right Hand Side\n", "B = mesh.cellGradBC # a matrix for boundary conditions\n", "Exbc = np.r_[0., 1.] # boundary values for Ex\n", "\n", "# A-matrix\n", "A = sp.vstack([\n", " sp.hstack([Grad, 1j*omega*Mmu]), # Top row of A matrix\n", " sp.hstack((Msighat, Div)) # Bottom row of A matrix\n", "])\n", "\n", "# Right-hand side\n", "rhs = np.r_[\n", " -B*Exbc, \n", " np.zeros(mesh.nC)\n", "] " ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "Ainv = Solver(A) # Factorize A matrix\n", "u = Ainv*rhs # Solve A^-1 rhs = u\n", "\n", "Ex = u[:mesh.nC] # Extract Ex from uution vector u\n", "Hy = u[mesh.nC:mesh.nC+mesh.nN] # Extract Hy from solution vector u \n", "\n", "P0 = sp.csr_matrix(\n", " (np.r_[1.], (np.r_[0], np.r_[len(u)-1])), shape=(1, len(u))\n", " )\n", "P0 = P0.tocsr()\n", "\n", "dAdsig_u = sp.vstack((Utils.spzeros(int(mesh.nC+1), mesh.nC), Utils.sdiag(Ex)))\n", "dudsig = - (Ainv * (dAdsig_u.toarray()))\n", "dZdsig = P0 * (Utils.sdiag(1./(u**2)) * dudsig)\n", "\n", "J = np.vstack((dZdsig.real, dZdsig.imag))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Plot the sensitivity " ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0,0.5,'Sensitivity')" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEOCAYAAACetPCkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xd81dX9+PHXOzd7kh1IgISQsEEgTEGGoFhXpYriqFqsq2jV1qq1X9evtlWr1dZV6qqtMsRRwQFqCUOg7L1kBAgrCQTIIPv8/vjcQBJCchNyV/J+Ph6fx733fNb7hst93/M553OOGGNQSimlHOXj7gCUUkp5F00cSimlmkQTh1JKqSbRxKGUUqpJNHEopZRqEk0cSimlmkQTh1JKqSbRxKGUUqpJNHEopZRqEk0cSimlmsTX3QE4Q0xMjElOTnZ3GM1WVFRESEiIu8NQbZR+/tqu1atX5xljYhvbrlUmjuTkZFatWuXuMJotMzOT0aNHuzsM1Ubp56/tEpG9jmynl6qUUko1iccnDhHpIiJvi8hsd8eilFLKTYlDRN4RkRwR2VSnfIKIbBeRnSLyKIAxZrcxZoo74lRKKXU2d7VxvAe8CrxfXSAiNuA1YDyQDawUkc+NMVvcEqFSqs0oLy8nOzubkpISd4fiEoGBgSQlJeHn59es/d2SOIwxi0QkuU7xYGCnMWY3gIjMAK4GNHEopZwqOzubsLAwkpOTERF3h+NUxhiOHj1KdnY2KSkpzTqGJ/WqSgT213idDQwRkWjgWaC/iDxmjPljfTuLyJ3AnQDx8fFkZmY6OVznKSws9Or4lXdri5+/iIgIoqOjKSwsdHcoLuHv78/x48eb/e/sSYmjvjRvjDFHgbsb29kYMw2YBpCRkWGa1Z1w33IoL4bQBAhLgKBIcMOvD+0OqdypLX7+tm7dSnh4uLvDcKnAwED69+/frH09KXFkAx1rvE4CDro0gsw/wu7MM699/CA0HsLirWQSGmcllND4s8ttzbtWqJRS5+u2227jiiuu4Nprr3XJ+TwpcawE0kQkBTgA3ADc6NIIrvobHN8PhUespeAwFOZA4WHIz4L9y6H4aP37BkfbayrxNRKLPalU12BC4yEg1KVvSSnlXYwxGGPw8fHcuyXckjhEZDowGogRkWzgSWPM2yIyFZgH2IB3jDGbXRpYu07W0pDK8jPJpDDHnlzqJJncHVZZVfnZ+/uFNFqD8Ss7AVVV4MEfHKVUy8nKyuKyyy5jzJgxLFu2jAceeIA333yT0tJSUlNTeffddwkNDeWZZ55hzpw5nDp1iuHDh/P3v//dLY357upVNfkc5V8CX7o4nKax+UFEorU0xBg4lX8msZxOLtXPj8CRTbDrv1B6stauFwIs94WQuHpqMHWfx4FvgNPerlJtzdNzNrPl4MnGN2yCnh3CefLKXg1us337dt59912eeeYZJk6cyLfffktISAjPPfccL730Ek888QRTp07liSeeAOCWW25h7ty5XHnllS0aqyM86VJV6yICwVHWEt+z4W3Lis4kk8Ij/LB2CWkJYWdqMCcOwIE1UJQLmLP3D4o8x2WyOpfMAsLd0tivlGpc586dGTp0KHPnzmXLli1ceOGFAJSVlTFs2DAAFixYwPPPP09xcTHHjh2jV69emjjaLP8QiOpiLcCB3Hak1derpbLCSh51L5PVrM3sXWY9Vpaevb9vUAOXyWo8D4kBH5tz37NSHqqxmoGzVI9IbIxh/PjxTJ8+vdb6kpIS7r33XlatWkXHjh156qmn3HbDoiYOb2LzhfD21tIQY6Dk+OkazNmXyQ5D7jbYsxBKTpy9v9ggJLbxGkxoAvgFOue9KtVGDR06lF/84hfs3LmTrl27UlxcTHZ2NnFxcQDExMRQWFjI7NmzXdaLqi5NHK2RiHX5KigS4ro3vG35KXtCqdvQXyPJHNoARTlgqs7avSooiqrUcdj6XoukjtVuyUqdp9jYWN577z0mT55Maal15eD3v/896enp/PznP6dPnz4kJyczaNAgt8UoxtRzzdzLZWRkGJ2Po3HGGIrLKjlZUs6JU+WcPFVhf7S/rlF+srgEivPwL84hoCSP4LI8IiqPkeJzmHE+q2knReSbML71GU6m/yh2B/YiKMCPYH9fgv1tBPvbCLI/D6nxPMjfRkid50H27a3FF39f7V3mSm31BsAePXq4OwyXqu89i8hqY0xGY/tqjcPLlVdWUVBy7i/82q/tS43tK6oa/uEQGuBLRJAfYYG+RARFER4XT1CQH/6BfvgH+XE8wMbM8jLicpaQlvM1V+dncl3pPI5WxrG0cjQLyi9iS0EnisurKC6r5FRZBcXllTTl94qvj5xOJjUTS5C/rz0JnVkXHepP94RwerQPJzZMe5sp5QyaONys+lf/6S/44nLW5lSQtzq7/l/+dRJBUVllg8f39REigvysL/8gPyKC/ekYFUxEkB/h9vLwQL/T24QH+Z4uCwv0xdfm6K/97sAdUFoI278keuNHXLnzY640syC2Owy8FnpfC1EpGGMoKa+iuKyC4rJK+1LBqbJKimo8ry6v3sZaf2bdiVPlHDp+ylpXXklRaQWlFWcup8WGBdCjfTg924fTo30YvTqEkxITis1He5YpdT40cbSA8sqqs37Jn+uX/8la6yvO/at/zfrTT2v/6vejY1QwvWt84df8so8Irl0e5Gdz7Q1CAaHQd5K1FB2FLZ/Bxtnw399bS2IG0uc6gnpdQ1BYPNEtfPr8ojK2Hj7JloMn2XqogC2HTvL2rt2UV1p/4wBfH7onhFkJpYNVM+meEEZYoLbNKOUoTRw15BSUcLSwrNYXe91EcLKeS0DFjfzq97PJ6S/2cPuv/k7RIYTbE0HdX/4/bFnP2BFDm/Gr38OERMOgKdZyfD9s+thKIl8/AvMeg5RR0Oc66HEFBEa0yCkjQ/wZnhrD8NSY02VlFVXsyi20J5OTbDl0knmbDzNj5ZnBmDtFBdOjfRg920dYjx3CSWwX1OqH2FaqOTRx1PCrWetZ/ENevevCAnwJt3/Jhwf60jk6uM4Xvu+Z13USQaCfT5O+gCoO2OgcHdJSb8sztOsIIx6wlpxtsGk2bPwI/nMvzH0Q0i+xkkjapS3exdff14ce7a3aRTVjDIdPlrD1kL1mYk8q87ccOd3+Eh7oe3q/nvYaSs/24fjopS7VxmniqOHuUancMKhTrcs/EUF+hAZ48a9+TxTXHcb+DsY8DgdWWwlk0yewdY51d3uPK6H3T6waic05H1ERoX1EEO0jghjbPf50eXFZBdsOn0kkWw+dZNaq/adrlZ2jg5mU0ZHrBiYRF673sKi2SRNHDRd2jWl8I9VyRCApw1oueRayFluXsrZ+Dus+sG5C7DXRqokkZbhkuJRgf18GdIpkQKfI02VVVYa9x4pZvTef2av388K87bz0zQ7Gdo/jhkEdGZUeqz8s1HkbPnw4S5cudXcYDtHEoTyDzRdSx1jL5S/CD/Oty1mr34MVf4d2naHPtVYSiXNtf3sfHyElJoSUmBCuHZjEnrwiZq7cz+zV2Xyz5Qjx4QFcN7AjkzI60ik62KWxqdbDW5IGgP5MUp7HLxB6XgWT3oeHf4AfvwHRqbDkL/D6UHjjQuv58X1uCS8lJoRHL+vOssfG8vdbBtKrQwSvZ+7kohcWcNNby/l8/UFKKxruMKFUXaGh1lw9mZmZjBo1ikmTJpGens6jjz7KBx98wODBg+nTpw+7du0CYM6cOQwZMoT+/fszbtw4jhw5AkBubi7jx49nwIAB3HXXXXTu3Jm8vPrbbptLaxzKswVGwAU3WkthDmz+1Lqc9e1T1tJxqFUT6XWNNTijC/nZfLi0VwKX9krg4PFTzF6dzcyV+7l/+lraBfsxsX8S1w/qSLeEMJfGpc7TV4/C4Y0te8yEPnDZnxzefP369WzdupWoqCi6dOnCHXfcwYoVK3jllVf429/+xssvv8yIESNYvnw5IsJbb73F888/z4svvsjTTz/N2LFjeeyxx/j666+ZNm1ay74XNHEobxIaB0PuspZje+zdez+CL38NXz0CqWOtBNPzapeP7tuhXRD3X5zG1DFd+X5XHjNW7udfy7N45/s99O/UjhsGdeSKvh0ICdD/cqpxgwYNon17azDT1NRULrnkEgD69OnDggULAMjOzub666/n0KFDlJWVkZKSAsCSJUv49NNPAZgwYQKRkZH1nOH86KdYeaeoFLjo1zDyV3Bks71772yYfTvEpMOoR6xaiIsTiI+PMDItlpFpsRwtLOXTtQeYsXI/j3y8kWfmbOHKfh24dXhyra7BysM0oWbgLAEBZ4bL8fHxOf3ax8eHiooKAO677z4eeughrrrqKjIzM3nqqacAq6u5s2kbh/JuIpDQG8Y9Bb/cANf90xoW/uMp8Powq1ZSdfaovq4QHRrAHSO78M2DF/HxPcO4rE97/rPuIFe9uoR/L9/rkv/gqvU6ceIEiYnWTKT//Oc/T5ePGDGCWbNmATB//nzy8/Nb/NyaOFTr4eMDvX4M9yyFa9+1ksrsn8Ebw6z7RNyUQESEgZ2j+PN1/Vj66Fgu7BrD7z7bxCMfb6CkXBvRVfM89dRTXHfddYwcOZKYmDPte08++STz589nwIABfPXVV7Rv356wsJZtZ9Nh1T1QWxzW2imqKq2xsjKfg7ztENfTuoTV4yorybja8f2w/A1MzhZ2nPBh9ZEqAsOiuLh/NyIioyGwnX2JgCD7Y2A78PV3aZht8fPXmoZVLy0txWaz4evry7Jly7jnnntYt27dWdvpsOpK1cfHZt2B3vPHVm+shc/BR7dCXC8Y/Qh0v9I1CeTIZvj+FasNRgRJ6Es3KSI59Bim+DiBS8sb3t83yEoi0V2thv+eV1uzMypVj3379jFp0iSqqqrw9/fnH//4R4ufQxOHav18bGe67G76xEogs34K8b1h9KPQ7fKWTyDGwN7vYcnLsPMb8AuBIXfD0HuscbuAAGBXbiFT31/G0bwcfn1RAtf1CkVKT1pT+p7Ktx5LTlhTAe9fCV89DF/9BpJHWJflelxl9TZTyi4tLY21a9c69RyaOFTb4WODvtdB74nWr/+Fz8HMm60+9qMehe6Xn/+wJlWVsG2uVcM4sBqCY6xxuTKmQHDUWZunxoby0dQxPPzRen6z8DCLj0fx3E/GEOx/jv+aOVth82ew+RP44lfw5cP2JHKNlURcfC+Laps0cai2x8cG/a63LmNtqk4gN1kJZPRj0O1HTU8g5SWwfjos/Rsc2wWRKXD5S9Z9JX5BDe4aGuDL6zcN4I2Fu/jzvO3sOFzA328ZSHJMPSMkx/WwltGP2pPIp1YSmfuglUhSLrKSSPcrrWHtlcOMMW1mGP3zbdvWXlWq7bL5Qr8b4Bcr4cdvWrMXzrgRpo2C7V/h0Py2p/Jh8Yvwch+Y+wAEhMF178F9q615SBpJGtVEhHtHd+W92wdzpKCEq15dwoJtOQ3tAPE9YezjMHUV3P09jHjIGoZlzi/hz2lWIik56djfoo0LDAzk6NGjbaKLtDGGo0ePEhjY/NGdtVeVB2qLvVo8QmUFbJgJi56H/Cxof4FVA0m/9OwayIkDsPx1axDGskJIvRgu/KX1i/88f7XuP1bMXf9azdbDJ3lwXDpTx3R1fA4QY6zhMla/B6vegbD2cPmfrctwDmqLn7/y8nKys7MpKSlxdyguERgYSFJSEn5+tWe+dLRXlSYOD9QW/+N6lMpyK4EsfB6O74UOA6wEkjYecrfB93+FjbOsL+neE2H4/dC+b4uGcKqskt9+upFP1x5gXI94Xrq+H+FNnd42exV8fj/kbLZ6Yl32PIQlNLqbfv7aLu2Oq1Rz2fyg/83Q93qr3WLRC/DhddbQ7sf3Wt1jM6bAsF9AZGenhBDkb+OlSf3olxTB77/YytWvfs+0WwaSFt+EG7mSMuCuhVZD/cLnYVcmXPL/YMBPXTK3iWq9tI1DqXOx+VlfsvetgSv/ChFJVs3jwc3wo+edljSqiQi3XZjChz8fSkFJBbe+s6Lpd5rb/Kwxve5ZajX+z7kf3rsC8nY6J2jVJmjiUKoxNj8YeCvc/qXVm8nFvZUGp0Txt8n9OXiihPeXZTXvIDFd4dY5VgI8vBHeGG7VpCrKWjJU1UZo4lDKCwxLjWZUeiyvLdjFiVON3Gl+Lj4+VgKcugK6TYD//t7qQZb3Q8sGq1o9TRxKeYlHJnTnZEk5by7cdX4HCkuwZle8Ybo1Oda/J0LBkZYJUrUJmjiU8hI9O4Rzdb8OvLNkD4dPtEC30e4/gps+gqI8+HASlBWd/zFVm6CJQykv8qtLulFlDK98t6NlDpg4wBqC/vAGawj6yoqWOa5q1Tw+cYhIFxF5W0RmuzsWpdytY1QwNw/tzMyV+9mZU3jex8vOL2Zt0BAODHsadnxN/scPcOBkJTtzCtmTV8S+o8UcK9IGdFWbU+/jEJF3gCuAHGNM7xrlE4BXABvwljHmnHM1GmN2A1M0cShlmTqmKx+tyuaFedv4+y2N3qt1TiuzjnHDtOVUVhkglUd9r+TuLf9ib3kFjy89cylMBC7uHscdI7swJCWqzYznpM7N2TcAvge8CrxfXSAiNuA1YDyQDawUkc+xksgf6+z/M2NMAwP2KNX2RIcGcOdFXXjpmx2s2ZfPgE6RTT5GUWkFv5q1ng7tAnnmqt5UGUNlZX8OLa/kt9nTGT8sg+zECVRWQVZeER+u2McN05bTOzGcO0Z04fK+7fGzefwFC+UkTk0cxphFIpJcp3gwsNNek0BEZgBXG2P+iFU7UUo1YsqIFN5ftpc/fbWNmXcObXIt4I9fbWV/fjEzfj6UIV1q3JfS7V2OvzqGQWsfY1CfntB5OABTx3bl07UHeGvxbh6YuY4/fbWN2y5MZvKgTkQEN3EoFOX1nD5WlT1xzK2+VCUi1wITjDF32F/fAgwxxkw9x/7RwLNYNZS37Ammvu3uBO4EiI+PHzhjxowWfieuU1hYSGhoqLvDUB7uu33l/GtLGQ8ODKBfrOO/ATflVfDnVaVcmuzL5O4BZ60vPX6Ykdufwa/8BGsGPMep4KTT66qMYWNeJfOyytlytIoAG4zr5MfEND9sjg7EqDzWmDFjPGOQw3oSx3XApXUSx2BjzH0tdU4d5FC1BeWVVYx7aSFBfja+uH+kQ1/cJ06Vc+lfFhEa6Mvc+0YQ6Gc7a5vMzExG9+0Mb40D/xC449t6ZxncfPAE0xbt5j/rDjIyLYZXbxxARJDWPryZo4McuuMiZTbQscbrJOCgG+JQyqv52Xz49SXd2Ha4gP+sO+DQPk/P2UxuYSkvTepXb9I4LSoFbpxl3SD44fXWiMF19OoQwSs39Oe5n/Rh2a6j/OSNpew7Wtzct6O8iDsSx0ogTURSRMQfuAH43A1xKOX1Lu/Tnj6JEbw4fwelFQ0PgDhv82E+WXOAX4zpSt+kdo0fPGkg/OgFOLjGmgb3HK4f1In3pwwmt6CUH7/+PauyjjX1bSgv49TEISLTgWVANxHJFpEpxpgKYCowD9gKzDLGbHZmHEq1Vj4+wiMTunPg+Cn+vXzfObc7WljKbz/ZSK8O4Uwd09XxE6RdYj1mN3zpd3hqDJ/eO5zwQF9u/Mf/+GytYzUg5Z2cmjiMMZONMe2NMX7GmCRjzNv28i+NMenGmFRjzLPOjEGp1m5EWozVxvDfHzhZcvYlJWMMj3+6iYKSCl6adAH+vk34bx8WDxEdG6xxVOsSG8qn915I/07teGDmOl6av52qqtY3UZzygjvHlVKNe2RCd/KLy/nHot1nrfvPuoN8vfkwD12STreEJkwEVS1xIBxwrLNJZIg//5oyhOsGJvHX/+7k4dkb2sQ83m2NJg6lWoHeiRFc2a8Dby3eQ87JM3d9Hz5RwhP/2cTAzpH8fGSX5h08cSAc3weFuQ5t7u/rw/PX9mXqmK58vCab2auzm3de5bE0cSjVSvxqfDrllVW88p01v4Yxhkc+3kB5peHF6/o1/z6LJHvvTAcuV1UTER4cn86QlCienrOF/ce0t1VroolDqVYiOSaEG4d0YsbK/ezOLWT6iv0s3JHLYz/qTnJMSPMP3L4fiM3hy1XVbD7Ci5P6AfCrWevtY2Kp1kATh1KtyH1j0wjw9eHxTzfx+y+2MKJrDDcPOc+50f1DIK5noz2r6pMUGczTV/ViRdYxptXT/qK8kyYOpVqR2LAA7hjZhWW7j2IT4flr++LTEkOBJA207ueoqmryrhMHJPKjPgm89M12tuw7ArnbYf+KZh1LeQZnj46rlHKxOy/qwqqsY9wytDMd2gW1zEETB8Lq9+DYLohJc2yfvJ2w6WMkP4tXineR7/cDce/UuDkwrieM+S10v8Iau115DU0cSrUyoQG+fPjzoS170MQaDeSOJo6vH4Gd30J4In6RyVR1GcOL24Uu6b25pm8cLPkLzLwZ2l8AYx6HtPGaQLyEXqpSSjUuthv4hzrezlFRClnfw5C74aEtcPuXJPz0HU4MfpAHt6bzfdilcO9y+PEbcCofPrwO3r4EspY4932oFqGJQynVOB8bdOjveM+q/f+DilPQZXSt4scu60GX2BB+/dF6CiuAC26Eqavgir/AyQPw3hWwfmZLR69amCYOpZRjEgfC4U1QXtL4trszrS68nS+sVRzkb+P3V/fm0IkSFmyzT+7p6w8ZP7MSSPII+Oxu2KgzRXsyTRxKKcckZUBVORze2Pi2uzMhaRAEhp+1anBKFBFBfizaUedOdP9guHEmdBwKn9wJmz9rmbhVi9PEoZRyTOJA67Gxy1Wn8uHg2rMuU1XztfkwIi2GhTtyzx7Hyj8EbpplJZ2Pp8DWOecdtmp5mjiUUo4J7wBhHRofemTPYjBV50wcAKPSYskpKGXb4YKzVwaEwU0fWW0qH90O2786r7BVy9PEoZRyXOKAxntW7c60emAlnXsG0ovSYwHOvlxVLTAcbv4YEvrArJ/Cvv81M2DlDJo4lFKOS8qA/D1Q3MAsf7szrUZu27nnH0+ICKR7QhgLz5U4AAIj4JZPIDQB5j4IlRXNj1u1KE0cSinHJTYyUu7xfdbd5V1GN3qoi9JjWZl1jKLSBhJCUCRM+APkbIaVbzU5XOUcmjiUUo7rcAEg575ctTvTeuwyutFDjUqPpbzSsHz30YY37H4FpI6FBc9CYU4TglXOoolDKeW4gDCI63HuGsfuTOvSUmz3Rg+VkRxJkJ+t4ctVYA1DctnzUH4Kvn266TGrFqeJQynVNIkDrMRRtyttVRXsXmjVNhwYcyrA18aw1OjGEwdY42MN+wWs+7c1sq5yK00cSqmmScyAU8esRvKacjZDcZ5Dl6mqjUqPZe/RYrLyihrf+KKHIawDVV/8mhW7ciku08Zyd9HEoZRqmuputtl1LlftWmA9dhnl8KFGVXfL/aHxWsf+Ih8+i78Xn8Pr+eydPzLqhUz+tXwv5ZU6r4eraeJQSjVNbA/wCz77DvLdmVbbRngHhw+VHBNCp6hgFm4/d+JYuCOX299dwUUvLOChzSlsD+zHk8Gz6RNZyf99tonxLy1kzvqDVOnUtC6jiUMp1TQ2X2sOjZoN5BWlsHdpky5TVRuVHsuy3Ucprag8a92c9Qe59Z0VbDp4kvvGpvH9oxfT7bbXCSg/wdv9d/LObRkE+tm4b/parn1zKSdLypv/vpTDNHEopZoucQAc2gAVZdbr/SvqHUbdEaPSYykuq2R1Vn6t8v3HivntJxvp36kd3z8ylofGp9M+IggSekNMOrLzW8Z2j+eL+0fy/LV92ZB9gnv/vUYvXbmAJg6lVNMlZUBlKRyxj5S7e0G9w6g7YlhqNH42qdW7qryyivtnrAXgrzf0x9+3zldVmn3Sp7IibD7CpIyO/GFiH5bszOPxTzeePXiialGaOJRSTXf6DvI11mMDw6g3JiTAl4zOUbUSx8vf7mDtvuP8YWIfOkYFn71T2nioLLMGVLSblNGR+8d2ZdaqbF5bsLPJcSjHOZQ4RORjEblcRDTRKKUgIglC4qw7yBsZRt0Ro7rFsu1wAUdOlrB0Zx6vZ+7i+oyOXNnvHA3tnYaBXwj8ML9W8YPj07mmfyJ/nr+Dz9YeaHY8qmGOJoI3gBuBH0TkTyLS+G2hSqnWS8S6XHVgtXXJqJFh1BtzUZrVLfeztQd4YOY6usSE8ORVPc+9g2+Adb4fvql1I6KI8Kef9GFIShS/mb2BDdnHmx2TOjeHEocx5ltjzE3AACAL+EZElorI7SJy7iEwlVKtV+IAOPqDNVNfI8OoN6ZH+zBiwwL409fbOF5czt8mDyDY37fhndLGw4l9kLu9VnGAr41pt2QQHuTLH77cqu0dTuDwpScRiQZuA+4A1gKvYCWSb5wSmVLKs1W3c2z+pNFh1BsjIoxKj8UY+O2PutOzgwNtJWnjrcedZ38FRQT7MXVMV5bvPsbiH/KaHZeqn6NtHJ8Ai4Fg4EpjzFXGmJnGmPuAUGcGqJTyUIkDrMfzvExV7Z7RqTxxRU9uHZ7s2A4RSRDX86x2jmqTh3QiKTKI5+dt05sDW5ijNY63jDE9jTF/NMYcAhCRAABjTPPrp0op7xUYATHp1vMuo8/7cKmxofxsRAriwACJp6WNh73LoPTsKWgDfG08ND6dTQdO8uWmQ+cdnzrD0cTx+3rKlrVkIEopL5Q8Etp1dmgYdafoOh6qyq1Reetx9QWJdIsP48X5O/TGwBbUYOIQkQQRGQgEiUh/ERlgX0ZjXbZyOhHpISJvishsEbnHFedUSjno0mfhzkyHhlF3ik5DwT/snJerbD7Cw5d2Y09eEbNXZ7s4uNarsRrHpcCfgSTgJeBF+/IQ8NvGDi4i74hIjohsqlM+QUS2i8hOEXm0oWMYY7YaY+4GJgF6WUwpT+IXBMFR7ju/zQ9Sx5zVLbemi3vEMbBzJJ98s4CS4kIXB9g6NZg4jDH/NMaMAW4zxoypsVxljPnEgeO/B0yoWSAiNuA14DKgJzBZRHqKSB8RmVtnibPvcxWwBPiu6W9RKdWqpY2HgoNwZHO9q0WEl5IW81H5/VS90g+WvQZlxS4OsnVpsKO0iNxsjPk3kCwiD9Vdb4x5qaH9jTGLRCS5TvFgYKcxZre5+YAFAAAczElEQVT9HDOAq40xfwSuOMdxPgc+F5EvgA8bOqdSqo3pWqNbbkLv2uuqqmD+7+i8+jW+tw2mnZTTa95vYcnLcOEvYchd59WNuK1q5A4bQuyPLdnlNhHYX+N1NjDkXBvb21MmAgHAlw1sdydwJ0B8fDyZmZktEKp7FBYWenX8yrt54+dvYGgK5n8fsKWgPSWBcSA+SFU5Pba+TFzuErITr+Cf5TezILuKDzL2kLZvJpHzH+fApiX8kH63u8P3Og0mDmPM3+1PXzfGODAxsEPqa0U7ZydrY0wmkNnYQY0x04BpABkZGWb06NHNi84DZGZm4s3xK+/mnZ+/GyDzjwz9313WJFOx3aGyHHI3wvhnSBp+PzfuyGX+uysp7jWJyIlTYf7vSFz6NxIzLocBt7j7DXiVxmoc1ZaKyB5gJvCJMSa/sR0akA10rPE6CTh4HsdTSrV1Fz0MqRdDzhbI2Wo9njwI10yDftcDMLRLNAG+PmRuz2V0tzi4+Ck4vBG+eMi6kTBpoHvfgxdxKHEYY9JEZDBwA/C4iGwBZtjbP5pqJZAmIinAAfsxb2zGcZRSyuJjg46DrOUcAv1sDEuNPjN8u80Xrn0Xpo2CmTfDXQshNM5FAXs3h8eqMsasMMY8hNW4fQz4Z2P7iMh0rBsFu4lItohMMcZUAFOBecBWYJYxpv7uEEop1YJGp8eyJ6+IvUeLrILgKLj+A2to+I9uP2eXXlWbo2NVhYvIrSLyFbAUOISVQBpkjJlsjGlvjPEzxiQZY962l39pjEk3xqQaY549r3eglFIOGt3NqlFkbq/RZNu+L0z4A+xdAjvmuSky7+JojWM9cAHwjP0L/xFjzOrGdlJKKU+SHBNCcnQwC7bn1F7R/xZo1wkWvaC1Dgc4mji6GGMeNMbo+FRKKa82ulscy3YdpaS88kyhzQ9GPAgHVlnzp6sGNTZW1cv2p5+LyFmLC+JTSqkWNapbLKUVVSzffbT2igtugrAOsOjP7gnMizTWq+pf9kf9SyqlWoVhdbvlVvMNsO4m//oRyPoeki90X5AerrGxqqrbMS4wxiysuWC1eSillFcJ9LMxtEuNbrk1DbwVQmKttg51To62cdxaT9ltLRiHUkq5zOhudbrlVvMLguH3We0c2avcE5wXaKyNY7KIzAFS6rRvLACONrSvUkp5qnE94hGB6Sv2n70yY4o1bMn6Ga4PzEs01sZRfc9GDNY8HNUKgA3OCkoppZypY1QwV/TtwL+WZXHXRV2IDPE/vS63zI9tFT3oumYui+Lu54q+HQgJcHR0prahsTaOvcaYTGPMsDptHGvsd4ArpZRX+sWYVIrKKnl3aVat8hfnb+e/5X1oX3mI1z/5hsHPfstrC3a6J0gP1dilqiX2xwIROVljKRCRk64JUSmlWl73hHAu7RXPu9/v4WRJOQCbDpxg5qr9xFzwIwD+NaqAYakxvDBvO19sOOTOcD1KYzWOEfbHMGNMeI0lzBgT7poQlVLKOe4bm0ZBSQXvL83CGMPTczYTFezPzZePhcgUOh1bxhs3D2BAp3Y88vEGdufq1LPg+FhVqSISYH8+WkTuF5F2zg1NKaWcq3diBGO7x/H2kj3MWrWflVn5/PrSbkQE+UHXi2HPYvxMOa/eOAA/m3DvB2tq33HeRjnaHfdjoFJEugJvAynoFK5KqVZg6tiu5BeX89gnG+nVIZxJGfbpgrqOg/Ii2LecDu2C+Mv1F7DtcAFPfa6DeTuaOKrsjeHXAC8bYx4E2jsvLKWUco0BnSIZ0TWGKgNPXtkLm499ktLkkeDjBzu/Bawxru68qAszVu5nx5ECN0bsfo4mjnIRmYx1I+Bce5nO8K6UahVenNSPt36aweCUqDOFAaHQaSjs+u/portHpRLo58M/Fu12Q5Sew9HEcTswDHjWGLPHPntfc2b/U0opjxMfHsi4nvFnr+g6Do5sgpNWj6qoEH8mZXTks3UHyDlZ4uIoPYdDicMYs8UYc78xZrr99R5jzJ+cG5pSSrlZ13HW467vThdNGZFCZZU56/6PtsTRXlUXisg3IrJDRHaLyB4Radt1NaVU6xffC0ITYOeZxNE5OoQJvRP4YPleCkvb5n3Qjl6qeht4CRgBDAIy7I9KKdV6iUDKSNi3vFbxnRelcrKkgpkr6xnrqg1wNHGcMMZ8ZYzJMcYcrV6cGplSSnmCxIFQcPB0OwfABR3bMTglineW7KGqqu1NNeto4lggIi+IyDARGVC9ODUypZTyBB3sX3UH19QqnpTRkQPHT7HtcNvrmuvokI9D7I8ZNcoMMLZlw1FKKQ/Tvi+IDQ6shu6Xny4emRYDwOIfcunZoW2NwORQ4jDGjHF2IEop5ZH8giC+JxyoXeOIL9rOX8JnkLcuFbr9BGJ7gK1tDL/uaK+qeBF5W0S+sr/uKSJTnBuaUkp5iMSB1qUqU6M9Y9ELXFP2OT/P/wu8OcJaStrGoOGOtnG8B8wDOthf7wAecEZASinlcToMgJITcMx+F0JFKexawMGuNzCm9EW2D3oG8nbAnF/WTi6tlKOJI8YYMwuoArCPW6VDRCql2obEgdbjgdXWY9YSKCskesCPOWBLZJYZD2Mfh82fwOp33ReniziaOIpEJBqrQRwRGQqccFpUSinlSWK7g2/QmXaOHV+DbxABaaMZkhLF4h9y4cIHIfVi+OpROLzRvfE6maOJ4yHgcyBVRL4H3gfuc1pUSinlSWy+0L6fVeMwBrZ/DaljwC+IkWkx7DhSyOGCMpg4DYLawZe/cXfETtXY1LGDRCTBGLMGGAX8FigF5gPZLohPKaU8Q+JAOLzBWk7sg/QJAIxMiwWsbrmExMDw+2HfUji03p3ROlVjNY6/A2X258OBx4HXgHxgmhPjUkopz5I4ACpKYPFL1uv0SwHonhBGTGgAi37Is8r73wx+IbD8TTcF6nyNJQ6bMeaY/fn1wDRjzMfGmP8Dujo3NKWU8iCJ9jvIt3wGHfpDWAIAIsKw1GhWZ9m/KoPawQU3wqbZUJjjpmCdq9HEISLVd7RcDPy3xrq2caeLUkoBRKZAUKT1PP2yWquSIoPIKSg9M27VkLugsgxWvePiIF2jscQxHVgoIv8BTgGLAexzj2uvKqVU2yFyZtyqbhNqrYoLC6CiynD8VLlVEJMGXcfDyretez5amQYThzHmWeBXWDcAjjDm9J0tPmivKqVUW9PrGkgZBQl9axXHhgUAkFtQI0kMvRuKcmD7l66M0CUa7Y5rjFlujPnUGFNUo2yHvaeV04nIaBFZLCJvishoV5xTKaXqNeAWuPVzq/ZRQ2yolThyCmpMJ5syGgIiYOe3LgzQNRy9j6NZROQdEckRkU11yieIyHYR2SkijzZyGAMUAoFoF2CllAeKCw8E6tQ4bL7WJFC7MlvdMCROTRxYl7hqXQwUERtWl97LgJ7AZPugiX1EZG6dJQ5YbIy5DHgEeNrJ8SqlVJPVe6kKrJsET2bD0V1uiMp5nNozyhizSESS6xQPBnYaY3YDiMgM4GpjzB+BKxo4XD4Q4Iw4lVLqfIT42wjys5FTJ3FUpozBBhxc8yW2oXcSb6+ZeDt3dKlNBGpO1JvNmYmiziIiE4FLgXbAqw1sdydwJ0B8fDyZmZktEatbFBYWenX8yrvp5695wvyq2LxrH5mZ1r0bJ8sMjy4qYp5PLFsXf8bUhak8MTSIpDBnX+hxPnckDqmn7JwXAI0xnwCfNHZQY8w07HezZ2RkmNGjRzc3PrfLzMzEm+NX3k0/f83TaetSsAmjRw8D4PudeRT/93/kdRzOmPzviDS+vLvDh8+njiAkwLtvg3NH6ssGOtZ4nQQcdEMcSinVYuLCA2q1cWTnFwPQMeNy/CoK+cd4YU9eEb/9dCPGyxvL3ZE4VgJpIpIiIv7ADVgj7yqllNeKDa2dOPYfO4XNR2jX82JA6FuylgfHpfOfdQeZu+GQ+wJtAc7ujjsdWAZ0E5FsEZlinwRqKtaMgluBWcaYzc6MQymlnC02LICTJRWUlFtz3O3PL6Z9RCC+YTHWkOy7F/CLMV1JjQ3hrSV73Bzt+XFq4jDGTDbGtDfG+Bljkowxb9vLvzTGpBtjUu13pyullFeLC6t9L8f+Y8V0jAy2VnYaBofW4+Mj/HRYMuv3H2fd/uPuCvW8eX/zvlJKeYDT93IU2hNH/ik6RgVZK8PbQ3kxlBYycUAiIf423l+W5Z5AW4AmDqWUagHViSPnZCkl5ZXkFpSeqXGExluPhUcIC/Rj4oAk5m44xNFC7xwAUROHUkq1gLgaNY7s/FMAJFXXOELjrMfCIwDcMqwzZRVVzFy1/6zjeANNHEop1QKiQvwRsdo49ld3xa2nxgGQHh/GkJQoZq/2zuH3NHEopVQL8LX5EB0SQG5BCdnH7Ikjqm7iODMj4Pie8ezOLeLg8VOuDvW8aeJQSqkWEhsWYK9xnMLf1+f0cOsERYHYTtc4AC7sGgNYd5h7G00cSinVQuLCAsgpKCU7v5ikdkH4+NhHWPLxsdo5aiSObvFhxIT6s3TXUTdF23zePWCKUkp5kNiwAHYcKcAYSKq+TFUtNK7WpSofH2FYagxLduZhjEGkvmH8PJPWOJRSqoVUX6rad6yYjpFBtVeGxteqcQCM6BpNacFRDqz+wqsme9LEoZRSLSQuLICKKsOJU+UkRTZc46CynEsLP2NhwEMkzb0Jlp1z1giPo4lDKaVaSPVNgMCZu8arhSZYiaOqynr91SO0W/h/7LZ1YUPQYPjmCdj5nQujbT5NHEop1UJO96Kixj0c1ULjwVTCqWPW673fQ9fxfNz7NaYUT8XEdofZP4P8vS6MuHk0cSilVAuJqzE1bMf6GscBCg5DZYU1D3l8LwanRJNb6suecf+A8lPw/SsujLh5NHEopVQLqb5UFexvIzLYr/bKmneP5++BqnKI7X46wWRVxULvibBhJpQWuDLsJtPEoZRSLSQ0wJdgfxsdI4PP7l57eryqHMjdbj2PTT/d+yo7/xQMugPKCq3k4cE0cSilVAtKiAikc3Tw2Stq1jhyt1nPY9KJCQ3A3+bDgfxTkDgQEvrCync8unuu3gColFIt6K839CcssJ6v1oBQ8AuxahzFeRCeCAFh+AAd2gWSffwUiFi1jjn3w77l0HmYy+N3hNY4lFKqBfVOjKBzdEj9K6uHHcndDjHpp4sTI4OsGgdAn2vBPww2zHBBtM2jiUMppVwlNN7qVZX3A8R2O12c2C6IA9Wj5PqHQMpI2L3QTUE2ThOHUkq5SmgcHN4A5UV1EkcwuQXWzIEApIyyel4d3+emQBumiUMppVwlLAFKT1rPY84kjiR7z6rTc3OkXGQ9emitQxOHUkq5SnWXXKhd47AnjtOXq+J6QEgc7NHEoZRSbVt1l9ygKAiJOV2c2M6eOKobyEWsWseeRR7ZLVcTh1JKuUp14qhR2wDr3g8fqVHjACtxFB7hyO71LgzQMZo4lFLKVaovVdXoigvgZ/MhITzwTI0D+CFkIADvf/BPsvKKXBaiIzRxKKWUq4QnAgIJfc5alRQZbA07AhwvLuNnn+dwkFj6VW3lp++sIKegxMXBnpsmDqWUcpXQOJgyHwb89KxViZHWvRzGGH45Yx1HTpQSktSLEVEnyC0o5TezN7gh4Ppp4lBKKVfqOBh8A84qTmwXxOGTJXy27gALd+Ty+OU9iOiQTnDhfu4Z1YXM7bns8ZBLVpo4lFLKAyRGBlFZZXjiP5vpnRjOzUM7Q2QKlBUwuXcwvj7Ch//zjEmeNHEopZQHqO6SW1BSwdNX9cLmIxCVAkBs+SEu6RXPR6uzz9xd7kaaOJRSygN0sk/oNLF/IgM7R1mFUV2sx/w93DykM8eLy/liwyE3RXiGJg6llPIAyTEh/G1yf566uteZwnadAYFjexiWGk1KTAgfr8l2W4zVNHEopZSHuLJfB8IDa0w56xcI4R0gfw8iwtjucazam+/2y1WaOJRSypNFpsCxPQAMT42mrKKKNfvy3RqSxycOERkpIm+KyFsistTd8SillEtFJcOx3QAMTonC5iMs23XUrSE5NXGIyDsikiMim+qUTxCR7SKyU0QebegYxpjFxpi7gbnAP50Zr1JKeZzIFCjKgdJCwgL96JMYwdLWnDiA94AJNQtExAa8BlwG9AQmi0hPEekjInPrLDXGIOZGYLqT41VKKc9i75JLfhYAF3aNZv3+4xSWVrgtJKcmDmPMIuBYneLBwE5jzG5jTBkwA7jaGLPRGHNFnSUHQEQ6ASeMMSedGa9SSnmcyOrEUd3OEUNFlWFlVt2vVtfxdcM5E4H9NV5nA0Ma2WcK8G5DG4jIncCdAPHx8WRmZp5HiO5VWFjo1fEr76afP8/iW17ICGDXym/ZfySMskpDkhxlzneLkUPR7onJDeeUesoanKnEGPNkYwc1xkwDpgFkZGSY0aNHNys4T5CZmYk3x6+8m37+PNDqdqRG+ZA6ejSUl5Cx5HbI8yGq++eQ0Nvl4bijV1U20LHG6yTgoBviUEop7xCTBnuXQmUFbP6EKJOPb1UJvHsZ7F/p8nDckThWAmkikiIi/sANwOduiEMppbzDsKmQuw3WvAfL3yAvuAsTSv+E8Q+B7552eTjO7o47HVgGdBORbBGZYoypAKYC84CtwCxjzGZnxqGUUl6t59WQPBLmPQ6HN5Cd/lMOmmjyet0GWYshZ6tLw3F2r6rJxpj2xhg/Y0ySMeZte/mXxph0Y0yqMeZZZ8aglFJeTwQm/AkqyyCwHb4XTAZgQ+xVYPOHlW+5NByPv3NcKaUUViP4j/4Ml79I5wSrN9WOggDoNRHWz4AS192toIlDKaW8xaAp0OdawgL9iAkNICuvCAbdAWWFsNV1TcWaOJRSygt1iQmxppJNyoCQWNid6bJza+JQSikvlBwTzO68Iqv9I+Ui2LMITIO3xLUYTRxKKeWFUmJCySsspaCkHFJGQeERyN3uknNr4lBKKS+UEhMCQFZeMXQZZRXuWeiSc2viUEopL1SdOHbnFUJksjXN7J5FLjm3Jg6llPJCnaODEcFqIAer1pG1GKqcP62sJg6llPJCgX42OkUFs+NIgVWQMgpKTsCRTQ3v2ALcMTquUkqpFtA9IYxth63EUdZlPHLfevyik51+Xq1xKKWUl+qWEE5WXhEl5ZVkZp2i78vb2HbY+XeQa+JQSikv1SMhjCoDPxwp5H97jlFpzOlGc2fSxKGUUl6qW0IYANsOn2TFnmP079iOAF+b08+riUMppbxU5+gQAv18WJWVz+aDJxjSxTVTyWriUEopL2XzEdLjw/h8/UGqDAxJiXLJeTVxKKWUF+sWH8ap8kp8fYT+ndq55JyaOJRSyot1bx8OQN+kCIL9XXOHhSYOpZTyYt3tDeSDU1zTvgGaOJRSyqv169iO4anRXH1BB5edU+8cV0opLxYa4MuHPx/q0nNqjUMppVSTaOJQSinVJJo4lFJKNYkmDqWUUk2iiUMppVSTaOJQSinVJJo4lFJKNYkmDqWUUk3S6m4AFJErgTwR2duE3SKAEy24bWPbNLY+BshzMB5P1JS/pyee73yP15z9Hd2npbZraL1+/tx/zvM53vl8/jo7tLUxplUtwDRn7uPIto1t48D6Ve7+O7r638CTzne+x3PmZ7CltmtovX7+3H/O8zmes78DjTGt8lLVHCfv48i2jW3TnBi9iavfX0uf73yP58zPYEtt15o/g+54b570GXT2dyBizzbKg4jIKmNMhrvjUG2Tfv5UY1pjjaM1mObuAFSbpp8/1SCtcSillGoSrXEopZRqEk0cSimlmkQTh1JKqSbRxOEFROTHIvIPEfmPiFzi7nhU2yIiPUTkTRGZLSL3uDse5X6aONxERN4RkRwR2VSnfIKIbBeRnSLyKIAx5jNjzM+B24Dr3RCuamWa+Pnbaoy5G5gEaDddpYnDjd4DJtQsEBEb8BpwGdATmCwiPWts8jv7eqXO13s04fMnIlcBS4DvXBum8kSaONzEGLMIOFaneDCw0xiz2xhTBswArhbLc8BXxpg1ro5VtT5N+fzZt//cGDMcuMm1kSpP1OoGOfRyicD+Gq+zgSHAfcA4IEJEuhpj3nRHcKrVq/fzJyKjgYlAAPClG+JSHkYTh2eResqMMeavwF9dHYxqc871+csEMl0bivJkeqnKs2QDHWu8TgIOuikW1fbo5085RBOHZ1kJpIlIioj4AzcAn7s5JtV26OdPOUQTh5uIyHRgGdBNRLJFZIoxpgKYCswDtgKzjDGb3Rmnap3086fOhw5yqJRSqkm0xqGUUqpJNHEopZRqEk0cSimlmkQTh1JKqSbRxKGUUqpJNHEopZRqEk0cyuOISKGLz/eliLSzL/fWKB8tInPP47gdRGR2E7Z/T0T2iMh6EdkhIu+LSOJ5nP/HNUdXFpFMEWl0WHQRad/U9y0ifxaRsc2JU3kfTRyqzbKPOuxjjPmRMeY40A64t7H9HGWMOWiMubaJuz1sjOkHdAPWAgvsd3E3x4+xhkdvqoeAfzRxn78BjzbjXMoLaeJQXkFEOovIdyKywf7YyV6eKiLLRWSliDxTXVsRkVD7dmtEZKOIXG0vTxaRrSLyOrAG6CgiWSISA/wJSBWRdSLygv3UofaZ77aJyAciIvbjZInIH0RkmYisEpEBIjJPRHaJyN01zrXJ/txm/1W+0f4e7mvo/RrLX4DDWPNjICKX2M+3RkQ+EpHQGrE8JyIr7EtXERkOXAW8YH8/qfZDX2ffZoeIjDzH6X8CfG0/9m0i8pmIzLHXhqaKyEMistb+d4+yx7sXiBaRhCb8syovpYlDeYtXgfeNMX2BDzgzWvArwCvGmEHUHpCvBLjGGDMAGAO8WP2lj/Vr/n1jTH/7F161R4FdxpgLjDEP28v6Aw9g/XLvAlxYY/v9xphhwGKsiZGuBYYCz9QT/51ACtC/xntwxBqguz2x/Q4YZ39Pq7BqBtVOGmMGY/2dXjbGLMUaZ+ph+/vZZd/O177dA8CTdU8mIilAvjGmtEZxb+BGrPk6ngWKjTH9sYYs+WmdWGv+fVQrpYlDeYthwIf25/8CRtQo/8j+/MMa2wvwBxHZAHyLNddEvH3dXmPMcgfPu8IYk22MqQLWAck11lUPALgR+J8xpsAYkwuUiEi7OscZB7xpHw8KY0zdSZTOpTrZDcVKXt+LyDrgVqBzje2m13gc1sDxPrE/rq7zXqq1B3LrlC2o8d5OAHPs5RvrHCMH6NDAuVUrofNxKG/V2CBrNwGxwEBjTLmIZAGB9nVFTThPzV/eldT+P1O9rqrOdlWc/X9LHIi5Pv2xpmsV4BtjzORzbGfO8byu6jjrvpdqpzjzd6q7D9R+r3XfZ6B9f9XKaY1DeYulWMN8g5UUltifL8e6Jk+N9QARQI49aYyh9q/zcykAwlog1vrMB+4WEV+A6raBc7E33N+PVQP4Gut9XigiXe3rg0UkvcYu19d4XGZ/3pz3s4P6ayKOSAc2NXNf5UU0cShPFGwf6rt6eQi4H7jdfunpFuCX9m0fAB4SkRVYX7In7OUfABkisgor0Wxr7KTGmKNYl4I21WgcbylvAfuADSKyHqvNoD4v2NfvAAYBY4wxZfbLRLcB0+1/g+VA9xr7BYjI/7D+Lg/ay2YAD9sbslNxgDGmCNhVnaAcJSJ+QFesthfVyumw6sqriUgwcMoYY0TkBmCyMeZqd8flSvbLcBnGmLwWOt41WJf4ftfEfQYYY/6vJWJQnk3bOJS3Gwi8au8xdRz4mZvj8XrGmE9FJLqJu/kCLzojHuV5tMahlFKqSbSNQymlVJNo4lBKKdUkmjiUUko1iSYOpZRSTaKJQymlVJNo4lBKKdUk/x8dPnBeuUPPNwAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(1, 1)\n", "\n", "core_inds = mesh.vectorCCx > -core_extent\n", "\n", "ax.loglog(-mesh.vectorCCx[core_inds], abs(J[0, core_inds]), label=\"real\")\n", "ax.loglog(-mesh.vectorCCx[core_inds], abs(J[1, core_inds]), label=\"imag\")\n", "\n", "ax.grid(True)\n", "ax.legend()\n", "ax.set_xlabel(\"Logarithmic Depth (m)\")\n", "ax.set_ylabel(\"Sensitivity\") " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compute sensitivity-vector products:\n", "\n", "For the 1D MT problem, the sensitivity matrix ($N\\times M$) is small, hence generating sensitivity is not be a big deal. However, for any 3D EM problem, generating the sensitivity matrix will require a huge amount of memory. To minimize that we only compute sensitivity-vector product. In forward case we compute:\n", "\n", "$$ \\mathbf{Jv} = \\mathbf{J}_{\\boldsymbol{\\sigma}} \\mathbf{v} $$\n", "\n", "Similarly, in adjoint case, we compute\n", "\n", "$$ \\mathbf{Jtv} = \\mathbf{J}_{\\boldsymbol{\\sigma}}^{T} \\mathbf{v} $$\n", "\n", "Computing $\\mathbf{Jv}$ and $\\mathbf{Jtv}$ are straight forward from above derivation. \n", "\n", "$$ \\mathbf{J}_{{\\boldsymbol{\\sigma}}}^T \\mathbf{v} \n", "= - \\left(\\frac{d \\mathbf{A}}{d {\\boldsymbol{\\sigma}}}\\mathbf{u} \\right)^T\n", "\\left(\\mathbf{A}^{T}\\right)^{-1} \\frac{\\partial P(\\mathbf{u})}{\\partial {\\mathbf{u}}}^T \\mathbf{v} $$\n", "\n", "One function computes forward sensitivity-vector product as `Jvec` and the other function computes adjoint sensitivity-vector product are named as `Jtvec`. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Jvec\n", "\n", "`Jvec` function takes conductivity ($\\sigma$) and vector ($\\mathbf{v}$), and computes sensitivity-vector product at a given frequency. " ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "def Jvec(sigma, v, f=100.):\n", " mu = np.ones(mesh.nC)*mu_0 # magnetic permeability values for all cells\n", " epsilon = np.ones(mesh.nC)*epsilon_0 # dielectric constant values for all cells\n", " omega = 2*np.pi*f # Angular frequency (rad/s)\n", " sigmahat = sigma # Assume sigmahat = sigma\n", " \n", " Div = mesh.faceDiv # Divergence matrix\n", " mesh.setCellGradBC([['dirichlet', 'dirichlet']]) # Setup boundary conditions\n", " Grad = mesh.cellGrad # Gradient matrix\n", " B = mesh.cellGradBC # a matrix for boundary conditions\n", " Exbc = np.r_[0., 1.] # boundary values for Ex\n", " Msighat = Utils.sdiag(sigmahat) \n", " Mmu = Utils.sdiag(mesh.aveCC2F * mu) \n", "\n", " # A-matrix\n", " A = sp.vstack([\n", " sp.hstack([Grad, 1j*omega*Mmu]), # Top row of A matrix\n", " sp.hstack((Msighat, Div)) # Bottom row of A matrix\n", " ])\n", "\n", " # Right-hand side\n", " rhs = np.r_[\n", " -B*Exbc, \n", " np.zeros(mesh.nC)\n", " ] \n", "\n", " Ainv = Solver(A) # Factorize A matrix\n", " u = Ainv*rhs # Solve A^-1 rhs = u\n", " Ex = u[:mesh.nC] # Extract Ex from uution vector u\n", " Hy = u[mesh.nC:mesh.nC+mesh.nN] # Extract Hy from solution vector u \n", " \n", " P0 = sp.csr_matrix(\n", " (np.r_[1.], (np.r_[0], np.r_[len(u)-1])), shape=(1, len(u))\n", " )\n", " P0 = P0.tocsr()\n", " Zxy = - 1./(P0*u)\n", "\n", " dAdsig_u_v = np.r_[np.zeros_like(Hy), Utils.sdiag(Ex)*v]\n", " dudsig_v = - (Ainv * (dAdsig_u_v))\n", " dZdsig_v = P0 * (Utils.sdiag(1./(u**2)) * dudsig_v)\n", " dZrdsig_v = dZdsig_v.real\n", " dZidsig_v = dZdsig_v.imag\n", " return np.r_[dZrdsig_v, dZidsig_v]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Order test: Jvec\n", "\n", "We have written the `Jvec` function, but how do we make sure this function is working properly? \n", "\n", "Let's consdier a predicted data $ d = \\mathcal{F}[\\sigma + \\triangle \\sigma] $.\n", "Applying Taylor's expansion yields\n", "\n", "$$ \n", "\\mathcal{F}[\\sigma + \\triangle \\sigma] = \\mathcal{F}[\\sigma]\n", "+\\frac{d \\mathcal{F}}{d \\sigma} \\triangle \\sigma \n", "+ \\mathcal{O}(\\triangle \\sigma )^2\n", "$$ \n", "\n", "By rearranging aboe equation, we can consider two misfit functions:\n", "\n", "$$\n", "f^1 = \\|\n", "\\mathcal{F}[\\sigma + \\triangle \\sigma] -\\mathcal{F}[\\sigma] \\|\n", "$$\n", "\n", "$$\n", "f^2 = \\|\n", "\\mathcal{F}[\\sigma + \\triangle \\sigma] -\\mathcal{F}[\\sigma]\n", "-\\frac{d \\mathcal{F}}{d \\sigma} \\triangle \\sigma \\|\n", "$$\n", "\n", "The first misfit function is supposed to have 1st order accuracy, but the other should be 2nd order accuracy. By using `SimPEG`'s `Tests` class we compute this two misfits, and check the accuracy. " ] }, { "cell_type": "code", "execution_count": 14, "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 3.454e-02 7.604e-03 nan\n", " 1 1.00e-02 4.121e-03 9.254e-05 1.915\n", " 2 1.00e-03 4.204e-04 9.461e-07 1.990\n", "========================= PASS! =========================\n", "Testing is important.\n", "\n" ] }, { "data": { "text/plain": [ "True" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from SimPEG import Tests\n", "def derChk(m):\n", " return [dpred(m), lambda mx: Jvec(m, mx)]\n", "Tests.checkDerivative(derChk, sigma, plotIt=False, num=3, eps=1e-20, dx=sigma*3)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Jtvec\n", "\n", "The below function takes conductivity ($\\sigma$) and vector ($\\mathbf{v}$), and computes adjoint sensitivity-vector product (`Jtvec`) at a given frequency. " ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "def misfit(sigma, dobs=None):\n", " r = dpred(sigma) - dobs\n", " return 0.5 * np.linalg.norm(r)**2\n", "\n", "def Jtvec(sigma, v, dtype=\"ri\"):\n", " f = 100.\n", " mu = np.ones(mesh.nC)*mu_0 # magnetic permeability values for all cells\n", " epsilon = np.ones(mesh.nC)*epsilon_0 # dielectric constant values for all cells\n", " omega = 2*np.pi*f # Angular frequency (rad/s)\n", " sigmahat = sigma # Assume sigmahat = sigma\n", "\n", " Div = mesh.faceDiv # Divergence matrix\n", " mesh.setCellGradBC([['dirichlet', 'dirichlet']]) # Setup boundary conditions\n", " Grad = mesh.cellGrad # Gradient matrix\n", " B = mesh.cellGradBC # a matrix for boundary conditions\n", " Exbc = np.r_[0., 1.] # boundary values for Ex\n", " Msighat = Utils.sdiag(sigmahat) \n", " Mmu = Utils.sdiag(mesh.aveCC2F * mu) \n", "\n", " tempUp = sp.hstack((Grad, 1j*omega*Mmu)) # Top row of A matrix\n", " tempDw = sp.hstack((Msighat, Div)) # Bottom row of A matrix\n", " A = sp.vstack((tempUp, tempDw)) # Full A matrix\n", " rhs = np.r_[-B*Exbc, np.zeros(mesh.nC)] # Right-hand side \n", "\n", " Ainv = Solver(A) # Factorize A matrix\n", " u = Ainv*rhs # Solve A^-1 rhs = u\n", " Ex = u[:mesh.nC] # Extract Ex from uution vector u\n", " Hy = u[mesh.nC:mesh.nC+mesh.nN] # Extract Hy from solution vector u \n", " P0 = sp.coo_matrix(\n", " (np.r_[1.], (np.r_[0], np.r_[len(u)-1])), shape=(1, len(u))\n", " )\n", " P0 = P0.tocsr()\n", " Zxy = - 1./(P0*u) \n", " ATinv = Solver(A.T) # Factorize A matrix \n", " \n", " PTvr = (P0.T*np.r_[v[0]]).astype(complex)\n", " PTvi = P0.T*np.r_[v[1]]*-1j\n", " dZrduT_v = Utils.sdiag((1./(u**2)))*PTvr\n", " dZiduT_v = Utils.sdiag((1./(u**2)))*PTvi\n", "\n", " dAdsiguT = sp.hstack((Utils.spzeros(mesh.nC, mesh.nN), Utils.sdiag(Ex)))\n", "\n", " dZrdsigT_v = - (dAdsiguT*(ATinv*dZrduT_v)).real\n", " dZidsigT_v = - (dAdsiguT*(ATinv*dZiduT_v)).real \n", " return dZrdsigT_v + dZidsigT_v" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Order test: Jtvec\n", "\n", "Similarly, `Jtvec` function has to be tested. For this, in this turn we consider a data msifit function:\n", "\n", "$$\n", "\\phi_d = \\frac{1}{2}\\|\n", "\\mathcal{F}[\\sigma] - \\mathbf{d}^{obs}\n", "\\|^2_2\n", "=\\frac{1}{2} \\mathbf{r}^T\\mathbf{r},\n", "$$\n", "\n", "where residual is $\\mathbf{r} = \\mathcal{F}[\\sigma] - \\mathbf{d}^{obs}$. \n", "\n", "By taking derivative w.r.t $\\sigma$, we obtain\n", "\n", "$$\n", "\\frac{d \\phi_d}{d \\sigma} = \\mathbf{J}_{\\sigma}^T \\mathbf{r}\n", "$$\n", "\n", "- Note that this is basically a gradient direction, and for first order optimization method such as steepest descent, we only needs this function that is only `Jvec` is required. \n", "\n", "Then applying taylor expansion to $\\phi_d$ yields\n", "\n", "$$ \n", "\\phi_d[\\sigma + \\triangle \\sigma] = \\phi_d[\\sigma]\n", "+\\frac{d \\phi_d}{d \\sigma} \\triangle \\sigma \n", "+ \\mathcal{O}(\\triangle \\sigma )^2\n", "$$ \n", "\n", "And similarly, we can consider two misfit functions:\n", "\n", "$$\n", "\\phi_d^1 = \\|\n", "\\phi_d[\\sigma + \\triangle \\sigma] -\\phi_d[\\sigma] \\|\n", "$$\n", "\n", "$$\n", "\\phi_d^2 = \\|\n", "\\phi_d[\\sigma + \\triangle \\sigma] -\\phi_d[\\sigma]\n", "-\\frac{d \\phi_d}{d \\sigma} \\triangle \\sigma \\|\n", "$$\n", "\n", "The first data misfit function is supposed to have 1st order accuracy, but the other should be 2nd order accuracy. By using `SimPEG`'s `Tests` class we compute this two misfits, and check the accuracy. " ] }, { "cell_type": "code", "execution_count": 16, "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 4.341e-02 6.094e-02 nan\n", " 1 1.00e-02 1.733e-03 1.981e-05 3.488\n", " 2 1.00e-03 1.754e-04 8.632e-08 2.361\n", " 3 1.00e-04 1.753e-05 1.094e-09 1.897\n", " 4 1.00e-05 1.753e-06 1.116e-11 1.991\n", "========================= PASS! =========================\n", "Testing is important.\n", "\n" ] }, { "data": { "text/plain": [ "True" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sigma0 = sigma*3\n", "dobs_ri = dpred(sigma)\n", "r = dpred(sigma0) - dobs_ri \n", "\n", "Tests.checkDerivative(\n", " lambda m: [misfit(m, dobs=dobs_ri), Jtvec(m, r)],\n", " sigma0,\n", " plotIt=False,\n", " num=5\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Adjoint test\n", "\n", "Both `Jvec` and `Jtvec` functions have passed the order test. These tests are necessary, but not sufficient. \n", "To test that both `Jvec` and `Jtvec` are consistent, we perform adjoint test. Consdier two random vectors: $\\mathbf{w}$ and $\\mathbf{v}$, then we can evalaute\n", "\n", "$$\n", "\\mathbf{w}^T \\mathbf{J}_{\\sigma} \\mathbf{v}, \n", "$$\n", "\n", "which will be a scalar value. Adjoint of above proucts is \n", "\n", "$$\n", "\\mathbf{v}^T \\mathbf{J}_{\\sigma}^T \\mathbf{w}, \n", "$$\n", "\n", "They should have same value: $\\mathbf{w}^T \\mathbf{J}_{\\sigma} \\mathbf{v}=\\mathbf{v}^T \\mathbf{J}_{\\sigma}^T \\mathbf{w}$. We evaluate $\\mathbf{w}^T \\mathbf{J}_{\\sigma} \\mathbf{v}$ and $\\mathbf{v}^T \\mathbf{J}_{\\sigma}^T \\mathbf{w}$ using `Jvec` and `Jtvec`, respectively, and check if they are outputing same values. " ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Adjoint Test 1.1102230246251565e-16 True\n" ] } ], "source": [ "v = np.random.rand(mesh.nC)\n", "w = np.random.rand(dobs_ri.shape[0])\n", "wtJv = w.dot(Jvec(sigma0, v))\n", "vtJtw = v.dot(Jtvec(sigma0, w))\n", "passed = np.abs(wtJv - vtJtw) < 1e-10\n", "print('Adjoint Test', np.abs(wtJv - vtJtw), passed)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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 }