{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simple notebook on performing a 1D MT problem.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#from IPython.display import Latex"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Maxwell's equations in 1D are as follows\n",
    "\n",
    "\n",
    "$i \\omega b = - \\partial_z e $ \n",
    "\n",
    "$  s =  \\partial_z(\\mu^{-1} b) - \\sigma(z) e$\n",
    "\n",
    "$b(0) = 1 \\hspace{1cm} ;\\hspace{1cm} b(-\\infty) = 0$\n",
    "\n",
    "\n",
    "\n",
    "where $e = \\widehat{\\overrightarrow{E}}_x $ and $b = \\widehat{\\overrightarrow{B}}_y$\n",
    "\n",
    "In weak form the equations become\n",
    "\n",
    "$i \\omega(b,f) = - (\\partial_ze,f)$ \n",
    "\n",
    "$(s,w) = -(\\mu^{-1} b, \\partial_z w) - (\\sigma(z) e, w)$\n",
    "\n",
    "\n",
    "where f and w are abritrary functions living in same discritizational space as b and e, respectivily.\n",
    "\n",
    "We consider e on nodes and b on cell centers. This way the derivative of any nodal function becomes\n",
    "\n",
    "$ (\\partial_z u)_k \\approx $\n",
    "\n",
    "$h_{k}^{-1} ( u_{k+\\frac{1}{2}} - u_{k-\\frac{1}{2}})  + O(h^2) $\n",
    "\n",
    "Matrix form\n",
    "\n",
    "$ e_z \\approx \\textbf{L}^{-1} \\textbf{G} e = \\begin{bmatrix} h_1^{-1} & & & \\\\\\\\ & h_2^{-1} & & \\\\\\\\ & &  \\ddots & \\\\\\\\ & & & h_n^{-1} \\end{bmatrix}^{(n,n)}\n",
    "\\begin{bmatrix} -1 & 1 & & & & \\\\\\\\ & -1 & 1 & & & \\\\\\\\ & & \\ddots & \\ddots & \\\\\\\\ & & & -1 & 1 \\end{bmatrix}^{(n,n+1)} \n",
    "\\begin{bmatrix} e_1 \\\\\\\\ \\\\\\\\ \\vdots \\\\\\\\ \\\\\\\\ e_{n+1} \\end{bmatrix}^{(n+1,1)} $\n",
    "\n",
    "where $ \\textbf{L} = diag(h) $ is the cell size and $ \\textbf{G}$ is the gradient operator with -1,1 representing the topology of the mesh, taking the difference between adjoint cells.\n",
    "\n",
    "We need to compute 2 inner products, on cell centers and from nodes to cell centers.\n",
    "\n",
    "Cell centers inner product is\n",
    "\n",
    "$ (b,f) \\approx \\sum\\limits_k h_k \\textbf{b}_k \\textbf{f}_k + O(h^2) $\n",
    "\n",
    "and in matrix from\n",
    "\n",
    "$ (b,) \\approx \\textbf{b}^T \\textbf{M}^f \\textbf{f}$  and $ (\\mu^{-1} b,f) \\approx \\textbf{b}^T \\textbf{M}_{\\mu}^f \\textbf{f}$ \n",
    "\n",
    "where $ \\textbf{M}_{\\mu}^f = diag(\\textbf{h} \\odot \\mu^{-1}) $ and $ \\textbf{M}^f = diag(\\textbf{h}) $ are the matrices.\n",
    "Nodes to cell centers inner product is\n",
    "\n",
    "$ (\\sigma e, w) \\approx \\sum\\limits_k \\frac{h_k \\sigma_k}{4} ( e_{k+\\frac{1}{2}} w_{k+\\frac{1}{2}} + e_{k+\\frac{1}{2}} w_{k+\\frac{1}{2}} )$\n",
    "\n",
    "and in matrix from\n",
    "\n",
    "$ (\\sigma e, w ) \\approx (\\textbf{h} \\odot \\sigma )^T ( \\textbf{A}_v (\\textbf{e} \\odot \\textbf{w} = \\textbf{w}^T diag(\\textbf{A}_v^T (\\textbf{h} \\odot \\sigma)) \\textbf{e} $\n",
    "\n",
    "Here $\\odot$ is a point wise Hadamard product and $\\textbf{A}_v$ is the averaging operator/matrix from nodes to cell centers\n",
    "\n",
    "$ \\textbf{A}_v = \\begin{bmatrix} \\frac{1}{2} & \\frac{1}{2} &  & \\\\\\\\ & \\ddots & \\ddots & \\\\\\\\ &  &  \\frac{1}{2} & \\frac{1}{2} \\end{bmatrix} ^{(n+1 , n)} $\n",
    "\n",
    "The sigma mass matrix is defined as \n",
    "\n",
    "$ \\textbf{M}_{\\sigma}^{e} = diag(\\textbf{A}_v^T (\\textbf{h} \\odot \\sigma)) $\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the MT problem there is no source in the domain, so $ s = 0 $. How ever the boundary conditions provide the right hand side where \n",
    "\n",
    "$ (\\partial_z \\mu^{-1} b, w ) = - (\\mu^{-1} b, \\partial_z w ) + (\\mu^{-1} b w )|_0^{end} $\n",
    "\n",
    "where \n",
    "\n",
    "$ (\\mu^{-1} b w )|_0^{end} = \\textbf{bc}^T (\\textbf{BC w}) $\n",
    "\n",
    "here $\\textbf{BC}$ is an matrix operator that extracts the boundary elements from $ \\textbf{w}$ and $\\textbf{bc} $ are the known boundary condintions. For the 1D case with homogenous boundary conditions we have \n",
    "\n",
    "$ \\textbf{B} = \\begin{bmatrix} -1 & 0 \\\\\\\\ \\vdots & \\vdots \\\\\\\\ 0 & 1 \\end{bmatrix} $ \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The weak form is \n",
    "\n",
    "$ (\\mu^{-1} b, \\partial_z w) + (\\sigma e, w) = (\\mu^{-1} b w )|_0^{end} $\n",
    "\n",
    "$ (i \\omega b,f) + (\\partial_z e , f) = 0 $\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the above matrix represntation we get the Maxwells equations in following form\n",
    "\n",
    "$ i \\omega \\textbf{f}^T \\textbf{M}^f \\textbf{b} + \\textbf{f}^T \\textbf{M}^f \\textbf{L}^{-1} \\textbf{G} \\textbf{e} = 0 $ \n",
    "\n",
    "$ \\textbf{w}^T \\textbf{G}^T \\textbf{L}^{-1}  \\textbf{M}^f_{\\mu} \\textbf{b} + \\textbf{w}^T \\textbf{M}_{\\sigma}^e \\textbf{e} = \\textbf{w}^T \\textbf{bc}^T \\textbf{BC} $ \n",
    "\n",
    "Here we use that \n",
    "\n",
    "$ (\\textbf{b},\\textbf{f}) \\approx \\textbf{b}^T \\textbf{M}^f \\textbf{f} = \\textbf{f}^T  \\textbf{M}^f \\textbf{b} $ \n",
    "\n",
    "since $\\textbf{M}^f$ is a symmetric diaoganal matrix of size n by n and $\\textbf{b}$ and $\\textbf{f}$ are vectors of length n."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We eliminate testing vectors and get system of equations to solve\n",
    "\n",
    "$ i \\omega \\textbf{b} + \\textbf{L}^{-1} \\textbf{G} \\textbf{e} = 0 $ \n",
    "\n",
    "$ \\textbf{G}^T \\textbf{L}^{-1}  \\textbf{M}^f_{\\mu} \\textbf{b} +  \\textbf{M}_{\\sigma}^e \\textbf{e} = \\textbf{bc}^T \\textbf{BC} $ \n",
    "\n",
    "and as $ \\textbf{A} \\textbf{x} = \\textbf{bc} $ system that we will solve, where \n",
    "\n",
    "$ \\textbf{A} = \\begin{bmatrix} \\textbf{G}^T \\textbf{L}^{-1}  \\textbf{M}^f_{\\mu}  & \\textbf{M}_{\\sigma}^e  \\\\\\\\ i \\omega  &  \\textbf{L}^{-1} \\textbf{G}  \\end{bmatrix} $\n",
    "\n",
    "$ \\textbf{x} = \\begin{bmatrix} \\textbf{b} \\\\\\\\ \\textbf{e} \\end{bmatrix} $\n",
    "\n",
    "$ \\textbf{bc} = \\begin{bmatrix} \\textbf{bc}^T \\textbf{BC} \\\\\\\\ \\textbf{0} \\end{bmatrix} $ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Efficiency Warning: Interpolation will be slow, use setup.py!\n",
      "\n",
      "            python setup.py build_ext --inplace\n",
      "    \n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "sys.path.append('C:/GudniWork/Codes/python/simpeg')\n",
    "import SimPEG as simpeg, numpy as np, scipy, scipy.sparse as sp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have\n",
    "\n",
    "$ i \\omega b = - \\partial_z e \\hspace{1cm} ; \\hspace{1cm} \\partial_z(\\mu^{-1} b) - \\sigma(z) e \\hspace{1cm} ;\\hspace{1cm} b(0) = 1 \\hspace{1cm} ;\\hspace{1cm} b(-\\infty) = 0 $\n",
    "\n",
    " \n",
    "To deal with boundary: we assume that below depth L both $ \\sigma $ and $ \\mu $ are constants ($ z < - L $). At the boundary we have that \n",
    "\n",
    "$ e = c \\exp(ikz) \\hspace{0.2cm} where \\hspace{0.2cm} k = \\sqrt{i\\omega\\mu\\sigma} $. \n",
    "\n",
    "Therefore for $ z < - L $ we have that\n",
    "\n",
    "$\\omega b - k e = 0 $.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We discretize the e field on the nodes and b field at the cell centers. The system we want to solve is \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\begin{bmatrix} i  \\omega & \\frac{\\partial}{\\partial z} \\\\\\\\ \\frac{1}{\\mu} \\frac{\\partial}{\\partial z} & -\\sigma \\end{bmatrix} \n",
    "\\begin{bmatrix} b \\\\\\\\ e \\end{bmatrix}  = \\begin{bmatrix} s1 \\\\\\\\ s2 \\end{bmatrix}$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    " # Set up  the problem\n",
    "mu = 4*np.pi*1e-7\n",
    "eps0 = 8.85e-12\n",
    "# Frequency\n",
    "fr = np.array([1e1]) #np.logspace(0,5,200) #np.array([2000]) #np.logspace(-4,5,82)\n",
    "omega = 2*np.pi*fr\n",
    "# Mesh\n",
    "sig0 = 1e-2\n",
    "#L = 3*np.sqrt(2/(mu*omega[0]*sig0))\n",
    "#nn=np.ceil(np.log(0.3*L + 1)/np.log(1.3))\n",
    "#h = 5*(1.3**(np.arange(nn+1)))\n",
    "\n",
    "h = np.ones(18)\n",
    "x0 = np.array([0])\n",
    "#sig = sig0*np.ones((len(h),1)) \n",
    "#sig[0:50] = 0.1\n",
    "#sig[50:100] = 1\n",
    "# Make the mesh\n",
    "mesh = simpeg.Mesh.TensorMesh([h],x0)\n",
    "sig = np.zeros(mesh.nC) + 1e-8\n",
    "sig[mesh.vectorCCx<=0] = 1e-2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 10.]), array([ 62.83185307]))"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fr,omega"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Make the operators\n",
    "G = mesh.nodalGrad\n",
    "Av = mesh.aveN2CC\n",
    "Li = scipy.sparse.spdiags(1/mesh.hx,0,mesh.nNx,mesh.nNx)\n",
    "Mmu = scipy.sparse.spdiags(mesh.hx/mu,0,mesh.nCx,mesh.nCx)\n",
    "Msig = scipy.sparse.spdiags(Av.T.dot(mesh.hx*sig.ravel()),0,mesh.nNx,mesh.nNx)\n",
    "# The boundaries\n",
    "bc_b = np.zeros((mesh.nCx,1))\n",
    "bc_b[0] = -1 # Set the top b field to 1\n",
    "bc_e = np.zeros((mesh.nNx,1))\n",
    "# Make the sparse matrix\n",
    "bc = sp.vstack((bc_b,bc_e))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/Gudni/anaconda/lib/python2.7/site-packages/scipy/sparse/linalg/dsolve/linsolve.py:90: SparseEfficiencyWarning: spsolve requires A be CSC or CSR matrix format\n",
      "  SparseEfficiencyWarning)\n"
     ]
    }
   ],
   "source": [
    "b = np.empty((mesh.nCx,len(omega)),dtype=np.complex64)\n",
    "e = np.empty((mesh.nNx,len(omega)),dtype=np.complex64)\n",
    "# Loop all the frequencies\n",
    "for nrOm, om in enumerate(omega):\n",
    "    # Left hand side\n",
    "    A = sp.vstack((sp.hstack(( -G.conj().T.dot(Mmu), - Msig)), sp.hstack((1j*om*scipy.sparse.identity(mesh.nCx) , G))))\n",
    "    #A = A.tocsr\n",
    "    # Solve the system\n",
    "    bef = scipy.sparse.linalg.spsolve(A,bc)\n",
    "    # Sort the output\n",
    "    b[:,nrOm] = bef[0:mesh.nCx]\n",
    "    e[:,nrOm] = bef[mesh.nCx::]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fa8671e1fd0>]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "# Plot the solution\n",
    "z=e[0,:]/(b[0,:]/mu)\n",
    "app_res = ((1./(8e-7*np.pi**2))/fr)*np.abs(z)**2\n",
    "app_phs = np.arctan(z.imag/z.real)*(180/np.pi)\n",
    "ax_res = plt.subplot(2,1,1)\n",
    "ax_res.loglog(fr,app_res)\n",
    "ax_phs = plt.subplot(2,1,2)\n",
    "ax_phs.semilogx(fr,app_phs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-5714285.5+0.00050081j], dtype=complex64)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Calculate the impedance\n",
    "z = e[0,:]/(b[0,:]/mu)\n",
    "z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  4.13555827e+17])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "app_res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}