OP_low.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sympy import *\n",
    "from sympy.abc import *\n",
    "from sympy.galgebra.ga import *\n",
    "import numpy as np\n",
    "from numpy import linalg as LA\n",
    "from __future__ import print_function\n",
    "init_printing()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Operational intensity of  differential operation\n",
    "\n",
    "We consder  differential operation on a vector $u$ at a given point in 3D with a 1D stencil size $k$ (number of points in the stencil) for every order, the subindex $i$ represent the dimension number $1$ for z, $2$ for $x$ and 3 for $y$, \n",
    "\n",
    "First order :\n",
    "$\n",
    "\\frac{d u}{dx_i}\n",
    "$\n",
    "\n",
    "Second order :\n",
    "$\n",
    "\\frac{d^2 u}{dx_i^2}\n",
    "$\n",
    "\n",
    "Second order cross derivative\n",
    "$\n",
    "\\frac{d^2 u}{dx_i dx_j}\n",
    "$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Arithmetic operations\n",
    "k = symbols('k')\n",
    "s = symbols('s')\n",
    "# 1D stencil\n",
    "#                multiplication         addition\n",
    "AI_dxi   =              k            +    k - 1\n",
    "AI_dxxi  =           k + 1           +    k - 1\n",
    "AI_dxxij =            2*k            +   2*k-1\n",
    "\n",
    "# square stencil (all uses the same stencil mask)\n",
    "#                multiplication         addition\n",
    "AI_dxis   =          k**2            +    k**2 - 1\n",
    "AI_dxxis  =          k**2            +    k**2 - 1\n",
    "AI_dxxijs =          k**2            +    k**2 - 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# I/O operations\n",
    "#                    load          \n",
    "IO_dxi   =           k             \n",
    "IO_dxxi  =           k             \n",
    "IO_dxxij =          2*k   \n",
    "IO_square =          k**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2*k - 1)/(4*k)\n",
      "1/2\n",
      "(4*k - 1)/(8*k)\n",
      "(2*k**2 - 1)/(4*k**2)\n",
      "(2*k**2 - 1)/(4*k**2)\n",
      "(2*k**2 - 1)/(4*k**2)\n"
     ]
    }
   ],
   "source": [
    "# Operational intensity in single precision\n",
    "\n",
    "print(AI_dxi/(4*IO_dxi))\n",
    "print(AI_dxxi/(4*IO_dxxi))\n",
    "print(AI_dxxij/(4*IO_dxxij))\n",
    "\n",
    "print(AI_dxis/(4*IO_square))\n",
    "print(AI_dxxis/(4*IO_square))\n",
    "print(AI_dxxijs/(4*IO_square))\n",
    "\n",
    "OI_dxi   = lambdify(k,AI_dxi/(4*IO_dxi))\n",
    "OI_dxxi  = lambdify(k,AI_dxxi/(4*IO_dxxi))\n",
    "OI_dxxij = lambdify(k,AI_dxxij/(4*IO_dxxij))\n",
    "\n",
    "\n",
    "OI_dxis   = lambdify(k,AI_dxis/(4*IO_dxxij))\n",
    "OI_dxxis  = lambdify(k,AI_dxxis/(4*IO_dxxij))\n",
    "OI_dxxijs = lambdify(k,AI_dxxijs/(4*IO_dxxij))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Operational intensity of  wave equations\n",
    "\n",
    "We now consider geophysical wave equations to obtain the theoretical expression of the operational intensity. We write directly the expression of a single time step as a function of differential operators. An operation on a wavefield is counted only once as we consider the minimum of arithmetic operations required.\n",
    "\n",
    "## Acoustic isotropic\n",
    "\n",
    "$ u(x,y,z,t+dt) = dt^2 v^2(x,y,z) ( 2 u(x,y,z,t) + u(x,y,z,t-dt) + \\nabla^2 u(x,y,z,t) +q )  $\n",
    "\n",
    "## VTI\n",
    "$ p(x,y,z,t+dt) = dt^2 v^2(x,y,z) \\left( 2 p(x,y,z,t) + p(x,y,z,t-dt) +(1+2\\epsilon)(\\frac{d^2 p(x,t)}{dx^2}+\\frac{d^2 p(x,t)}{dyx^2}) + \\sqrt{(1+2\\delta)} \\frac{d^2 r(x,t)}{dz^2} + q \\right)  $\n",
    "$ r(x,y,z,t+dt) = dt^2 v^2(x,y,z) \\left( 2 r(x,y,z,t) + r(x,y,z,t-dt) +\\sqrt{(1+2\\delta)}(\\frac{d^2 p(x,t)}{dx^2}+ \\frac{d^2 p(x,t)}{dy^2}) +  \\frac{d^2 r(x,t)}{dz^2} + q \\right)   $\n",
    "\n",
    "\n",
    "## TTI\n",
    "$  p(x,y,z,t+dt) = dt^2 v^2(x,y,z) \\left( 2 p(x,y,z,t) + p(x,y,z,t-dt) + (1+2\\epsilon) (G_{\\bar{x}\\bar{x}} + G_{\\bar{y}\\bar{y}}) p(x,y,z,t) + \\sqrt{(1+2\\delta)} G_{\\bar{z}\\bar{z}} r(x,y,z,t) + q \\right)  $\n",
    "$  r(x,y,z,t+dt) = dt^2 v^2(x,y,z) \\left( 2 r(x,y,z,t) + r(x,y,z,t-dt) + \\sqrt{(1+2\\delta)}(G_{\\bar{x}\\bar{x}} + G_{\\bar{y}\\bar{y}}) p(x,y,z,t) + G_{\\bar{z}\\bar{z}} r(x,y,z) +q  \\right)  $\n",
    "\n",
    "where\n",
    "$\n",
    "\\begin{cases} \n",
    "G_{\\bar{x}\\bar{x}} & = cos(\\phi)^2 cos(\\theta)^2 \\frac{d^2}{dx^2} +sin(\\phi)^2 cos(\\theta)^2 \\frac{d^2}{dy^2}+ sin(\\theta)^2 \\frac{d^2}{dz^2} + sin(2\\phi) cos(\\theta)^2 \\frac{d^2}{dx dy} - sin(\\phi) sin(2\\theta) \\frac{d^2}{dy dz} -cos(\\phi) sin(2\\theta) \\frac{d^2}{dx dz} \\\\\n",
    "G_{\\bar{y}\\bar{y}} & =  sin(\\phi)^2 \\frac{d^2}{dx^2} +cos(\\phi)^2  \\frac{d^2}{dy^2} - sin(2\\phi)^2 \\frac{d^2}{dx dy}\\\\\n",
    "G_{\\bar{z}\\bar{z}} & = cos(\\phi)^2 sin(\\theta)^2 \\frac{d^2}{dx^2} +sin(\\phi)^2 sin(\\theta)^2 \\frac{d^2}{dy^2}+ cos(\\theta)^2 \\frac{d^2}{dz^2} + sin(2\\phi) sin(\\theta)^2 \\frac{d^2}{dx dy} + sin(\\phi) sin(2\\theta) \\frac{d^2}{dy dz} +cos(\\phi) sin(2\\theta) \\frac{d^2}{dx dz} \\\\ \n",
    "\\end{cases}\n",
    "$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Arithmetic\n",
    "#                    dxi             dxxi             dxxij        multiplications     additions   duplicates \n",
    "AI_acou =        0*AI_dxi    +    3*AI_dxxi   +    0*AI_dxxij    +        3          +     5     -  2 * 2\n",
    "AI_vti  = 2 * (  0*AI_dxi    +    3*AI_dxxi   +    0*AI_dxxij    +        5          +     5     -  2    )\n",
    "AI_tti  = 2 * (  0*AI_dxi    +    3*AI_dxxi   +    3*AI_dxxij    +        44         +     17    -  8    )\n",
    "\n",
    "AI_acoums =        0*AI_dxi  +  3*s*AI_dxxi   +    0*AI_dxxij    +        3*s        +     5*s   -  2 * 2 *s\n",
    "AI_vtims  = 2 * (  0*AI_dxi  +  3*s*AI_dxxi   +    0*AI_dxxij    +        5*s        +     5*s   -  2*s  )\n",
    "AI_ttims  = 2 * (  0*AI_dxi  +  3*s*AI_dxxi   +  3*s*AI_dxxij    +        44*s       +     17*s  -  8*s  )\n",
    "\n",
    "AI_acous =        0*AI_dxis  +    3*AI_dxxis  +    0*AI_dxxijs    +        3          +     5     -  2 * 2\n",
    "AI_vtis  = 2 * (  0*AI_dxis  +    3*AI_dxxis  +    0*AI_dxxijs    +        5          +     5     -  2 * 2  )\n",
    "AI_ttis  = 2 * (  0*AI_dxis  +    3*AI_dxxis  +    3*AI_dxxijs    +        44         +     17    -  3*k**2 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# I/O operations (we load a point once only)\n",
    "#                     dxi              dxxi             dxxij         duplicate   other load/write\n",
    "IO_acou   =         0*IO_dxi    +    3*IO_dxxi   +    0*IO_dxxij  -     2       +     3\n",
    "IO_vti    =  2 * (  0*IO_dxi    +    3*IO_dxxi   +    0*IO_dxxij  -     2       +     2      ) + 3\n",
    "IO_tti    =  2 * (  0*IO_dxi    +    3*IO_dxxi   +    3*IO_dxxij  -  3*k +2     +     4     )  + 7\n",
    "\n",
    "IO_acoums =         0*IO_dxi    +    3*s*IO_dxxi +    0*IO_dxxij  -     2*s     +     3*s+1\n",
    "IO_vtims  =  2 * (  0*IO_dxi    +    3*s*IO_dxxi +    0*IO_dxxij  -     2*s     +     2*s   ) + 3\n",
    "IO_ttims  =  2 * (  0*IO_dxi    +    3*s*IO_dxxi +  3*s*IO_dxxij  -  s*(3*k +2) +     4*s   )  + 7\n",
    "\n",
    "IO_acous  =         0*IO_square +    3*IO_square +    0*IO_square -     2       +     3\n",
    "IO_vtis   =  2 * (  0*IO_square +    3*IO_square +    0*IO_square -     2       +     2      ) + 3\n",
    "IO_ttis   =  2 * (  0*IO_square +    3*IO_square +    3*IO_square -  3*IO_square+     4     )  + 7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(3*k + 2)/(2*(3*k + 1))\n",
      "(3*k + 4)/(3*(2*k + 1))\n",
      "(9*k + 25)/(12*k + 19)\n",
      "s*(3*k + 2)/(2*(3*k*s + s + 1))\n",
      "s*(3*k + 4)/(3*(2*k*s + 1))\n",
      "s*(9*k + 25)/(12*k*s + 4*s + 7)\n",
      "(6*k**2 + 1)/(4*(3*k**2 + 1))\n",
      "1/2\n",
      "(9*k**2 + 55)/(6*(2*k**2 + 5))\n"
     ]
    }
   ],
   "source": [
    "print(simplify(AI_acou/(4*IO_acou)))\n",
    "print(simplify(AI_vti/(4*IO_vti)))\n",
    "print(simplify(AI_tti/(4*IO_tti)))\n",
    "\n",
    "print(simplify(AI_acoums/(4*IO_acoums)))\n",
    "print(simplify(AI_vtims/(4*IO_vtims)))\n",
    "print(simplify(AI_ttims/(4*IO_ttims)))\n",
    "\n",
    "print(simplify(AI_acous/(4*IO_acous)))\n",
    "print(simplify(AI_vtis/(4*IO_vtis)))\n",
    "print(simplify(AI_ttis/(4*IO_ttis)))\n",
    "\n",
    "OI_acou  = lambdify(k,AI_acou/(4*IO_acou))\n",
    "OI_vti   = lambdify(k,AI_vti/(4*IO_vti))\n",
    "OI_tti   = lambdify(k,AI_tti/(4*IO_tti))\n",
    "\n",
    "OI_acoums  = lambdify((k,s),AI_acoums/(4*IO_acoums))\n",
    "OI_vtims   = lambdify((k,s),AI_vtims/(4*IO_vtims))\n",
    "OI_ttims   = lambdify((k,s),AI_ttims/(4*IO_ttims))\n",
    "\n",
    "\n",
    "OI_acous  = lambdify(k,AI_acous/(4*IO_acous))\n",
    "OI_vtis   = lambdify(k,AI_vtis/(4*IO_vtis))\n",
    "OI_ttis   = lambdify(k,AI_ttis/(4*IO_ttis))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1/2\n",
      "1/2\n",
      "3/4\n",
      "11/20\n",
      "13/18\n",
      "13/10\n",
      "1/2\n",
      "1/2\n",
      "3/4\n"
     ]
    }
   ],
   "source": [
    "print(limit(OI_acou(k),k,oo))\n",
    "print(limit(OI_vti(k),k,oo))\n",
    "print(limit(OI_tti(k),k,oo))\n",
    "\n",
    "print(limit(OI_acoums(3,s),s,oo))\n",
    "print(limit(OI_vtims(3,s),s,oo))\n",
    "print(limit(OI_ttims(3,s),s,oo))\n",
    "\n",
    "print(limit(OI_acous(k),k,oo))\n",
    "print(limit(OI_vtis(k),k,oo))\n",
    "print(limit(OI_ttis(k),k,oo))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "kk=[3,5,7,9,11,13,15,17,19,21,23,25,27,29,31]\n",
    "ss=[2,4,8,16,32,64]\n",
    "OI_wave=np.zeros((15,6))\n",
    "OI_wavems=np.zeros((15,6,3))\n",
    "OI=np.zeros((15,3))\n",
    "for i in range(0,15):\n",
    "    OI_wave[i,0]=OI_acou(kk[i])\n",
    "    OI_wave[i,1]=OI_vti(kk[i])\n",
    "    OI_wave[i,2]=OI_tti(kk[i])\n",
    "    OI_wave[i,3]=OI_acous(kk[i])\n",
    "    OI_wave[i,4]=OI_vtis(kk[i])\n",
    "    OI_wave[i,5]=OI_ttis(kk[i])\n",
    "    OI[i,0]=OI_dxi(kk[i])\n",
    "    OI[i,1]=OI_dxxi(kk[i])\n",
    "    OI[i,2]=OI_dxxij(kk[i])\n",
    "    for j in range(0,6):\n",
    "        OI_wavems[i,j,0]=OI_acoums(kk[i],ss[j])\n",
    "        OI_wavems[i,j,1]=OI_vtims(kk[i],ss[j])\n",
    "        OI_wavems[i,j,2]=OI_ttims(kk[i],ss[j])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "plt.hold(\"off\")\n",
    "acou = plt.plot(OI_wave[:,0],label='acou')   # this is how you'd plot a single line...\n",
    "vti  = plt.plot(OI_wave[:,1],label='vti')   # this is how you'd plot a single line...\n",
    "tti  = plt.plot(OI_wave[:,2],label='tti')   # this is how you'd plot a single line...\n",
    "fig = plt.figure()\n",
    "plt.hold(\"off\")\n",
    "acou = plt.plot(OI_wave[:,3],label='acous')   # this is how you'd plot a single line...\n",
    "vti  = plt.plot(OI_wave[:,4],label='vtis')   # this is how you'd plot a single line...\n",
    "tti  = plt.plot(OI_wave[:,5],label='ttis')   # this is how you'd plot a single line...\n",
    "fig = plt.figure()\n",
    "plt.hold(\"off\")\n",
    "acou = plt.plot(OI_wavems[:,2,0],label='acous')   # this is how you'd plot a single line...\n",
    "vti  = plt.plot(OI_wavems[:,2,1],label='vtis')   # this is how you'd plot a single line...\n",
    "tti  = plt.plot(OI_wavems[:,2,2],label='ttis')   # this is how you'd plot a single line...\n",
    "fig = plt.figure()\n",
    "plt.hold(\"off\")\n",
    "acou = plt.plot(OI_wavems[:,5,0],label='acous')   # this is how you'd plot a single line...\n",
    "vti  = plt.plot(OI_wavems[:,5,1],label='vtis')   # this is how you'd plot a single line...\n",
    "tti  = plt.plot(OI_wavems[:,5,2],label='ttis')   # this is how you'd plot a single line...\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}