Lecture-11-Scalars-Vectors.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Lecture 11:  Scalar and Vector Fields, Accumulation "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "## Background\n",
    "\n",
    "Engineers and scientists are typically interested in how certain quantities vary within space.  A scalar field is a mapping of values to positions within some coordinate system.  In physical problems it is common to map values to three dimensions or fewer.  For example, scalar fields may be generated from temperature measurements within a chemical reactor, pressure data over the Earth's surface, or x-ray intensity from telescope observations.  Any mapping of a numerical quantity to a spatial coordinate is a scalar field.\n",
    "\n",
    "Materials engineers interested in phase transformations typically deal with scalar fields corresponding to concentration as a function of position.  This could be indicated using the symbols $c(x,y)$ where the $c$ corresponds to the field of interest and the $(x,y)$ indicating the dimensions of and position within the field."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "\n",
    "## What Skills Will I Learn?\n",
    "\n",
    "* Define a scalar and vector field\n",
    "* Define and use the gradient and divergence operators\n",
    "* Develop a geometric intuition for the divergence and gradient\n",
    "* Learn or review the definition of the chemical potential and ideal solution models\n",
    "\n",
    "## What Steps Should I Take?\n",
    "\n",
    "1. Read through all of the activities below.  The computer code below is only for visualization and you are not required to understand it all completely before you finish this lecture's assignment.\n",
    "1. Understand how to use the $\\vec{\\nabla}$ operator to compute the gradients and divergence of scalar and vector fields respectively.\n",
    "1. Use the visualizations to build your geometric understanding for the concepts of a scalar and vector field.\n",
    "1. Compute the chemical potential from the Gibbs free energy and then apply the gradient and divergence operators to that expression and analyze what you get.\n",
    "\n",
    "## Reading and Reference\n",
    "* Essential Mathematical Methods for Physicists, H. Weber and G. Arfken, Academic Press, 2003\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quick Review:  Product Rule and Chain Rule\n",
    "\n",
    "There are three ways that functions can be combined:  addition, multiplication and compostion (one function placed as an argument into another function).  The product and chain rule are the result of computing differentials for products and compositions.  The sum rule is stated as:\n",
    "\n",
    "$$\n",
    "{\\frac {d(af+bg)}{dx}}=a{\\frac {df}{dx}}+b{\\frac {dg}{dx}}\n",
    "$$\n",
    "\n",
    "The product rule is:\n",
    "\n",
    "$$\n",
    "{\\dfrac{d}{dx}}(u\\cdot v)={\\dfrac {du}{dx}}\\cdot v+u\\cdot {\\dfrac {dv}{dx}}\n",
    "$$\n",
    "\n",
    "Assuming that $z(y)$ and $y(x)$ then the chain rule is stated as:\n",
    "\n",
    "$$\n",
    "\\frac{dz}{dx}={\\frac{dz}{dy}}\\cdot {\\frac{dy}{dx}}\n",
    "$$\n",
    "\n",
    "Reviewing these rules in the context of geometry will help build a better intuition for why these are the results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Scalar Fields\n",
    "\n",
    "A scalar field is just a number linked to an independent variable or variables.  The Python code here is just to help us visualize the field."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "$$\n",
    "f(x,y) = \\sin(2x) \\cdot \\sin(y)\n",
    "$$\n",
    "\n",
    "Note that this function takes an $(x,y)$ position as input and returns a single scalar value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "\n",
    "import sympy as sp\n",
    "from sympy.plotting import plot3d\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "x, y = sp.symbols('x y')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We use `plot3d()` to visualize this field so that we use the height of the surface as a proxy for the magnitude of the scalar field.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Plot our scalar function over a specified range.\n",
    "plot3d(sp.sin(2*x)*sp.sin(y), (x, -3, 3), (y, -2, 2));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Projections Onto Two Dimensions\n",
    "----\n",
    "\n",
    "* A field $f(x,y)$ contains three pieces of information:  f, x, y\n",
    "* Colors and other glyphs can help access the additional information\n",
    "* Contours and heatmaps are two such methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "An alternative to the three dimensional plot is to project the scalar values onto a two dimensional surface.  Two common options for this are _contour plots_ and _heat maps_.  Rather than using a third dimension to represent the scalar values a contour plot traces single-valued lines through the domain whereas a heat-map uses colors to represent the value of the scalar field.\n",
    "\n",
    "\n",
    "### Visualizing the Scalar Field by Contours\n",
    "----\n",
    "\n",
    "As previously stated, _iso-lines_ are plotted within the domain and the color, annotations and position of the lines quantify the scalar field.  Topographic maps use this representation to help the reader understand the locations and incline of mountains and valleys.  `matplotlib` has a `.contour()` method that will generate contour plots.  We demonstrate this next."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "delta = 0.025\n",
    "xnp = np.arange(-3.0, 3.0, delta)\n",
    "ynp = np.arange(-2.0, 2.0, delta)\n",
    "X, Y = np.meshgrid(xnp, ynp)\n",
    "Z = np.sin(2*X)*np.sin(Y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "contours = 10\n",
    "plt.figure()\n",
    "CS = plt.contour(X, Y, Z, contours)\n",
    "plt.clabel(CS, inline=1, fontsize=10)\n",
    "plt.title(r'A Simple Contour Plot of Your Scalar Field')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "### Visualizing the Scalar Field by Value\n",
    "----\n",
    "\n",
    "Color choice requires consideration of the medium in which you distribute your data as well as the capabilities of the reader.  An [essay](http://blogs.nature.com/onyourwavelength/2018/07/18/heuristics-for-better-figures/) has been written on the topic and provides some guidance for preparation of figures.  `matplotlib` has a module called `cm` that provides some standard color maps.  User specified maps are possible."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### DIY:  Using Surface and Contour Plots\n",
    "----\n",
    "\n",
    "Plot the electric field around a point charge if the potential is given by:\n",
    "\n",
    "$$\n",
    "V(x,y) = \\frac{1}{x^2+y^2}\n",
    "$$\n",
    "\n",
    "Use three dimensional surface plots and contour plots to help visualize this potential."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "### Gradients\n",
    "----\n",
    "\n",
    "Scalar fields have associated vector fields.  One such vector field is known as the _gradient_ and is indicated with the symbol $\\nabla$ and is called [nabla](https://en.wikipedia.org/wiki/Nabla_symbol).  The gradient is a type of directional derivative and is a vector quantity.  In a physical problem of three dimensions the gradient can be thought of as pointing in the direction of the maximum spatial rate of change.  Each basis vector magnitude is multiplied by the partial derivative of the field with respect to the basis vector's coordinate.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "You can visualize the gradient operator as a vector with components:\n",
    "\n",
    "$$\\overrightarrow{\\nabla} = \\frac{\\partial}{\\partial x} \\hat{i}  + \\frac{\\partial}{\\partial y} \\hat{j} + \\frac{\\partial}{\\partial z} \\hat{k}  $$\n",
    "\n",
    "When applied to a scalar field the result is a vector field - the gradient.  Geometrically, the gradient is a vector field that \"points uphill\".  The following illustrations are meant to convey some of that geometric intuition."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "### Quiver Plots\n",
    "----\n",
    "\n",
    "One way to visualize a vector field in a physical problem is to use small oriented arrow to indicate direction and length to indicate magnitude.  These are known as quiver plots and they can be used in conjunction with heat maps to help visualize vector fields.  In the example below the scalar field is sampled with two different point densities.  A higher point density (100 points) for the filled contour plot and a lower point density (30 points) for the quivers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "x0, x1 = (-3,3)\n",
    "y0, y1 = (-2,2)\n",
    "\n",
    "# Read the docstrings to understand why the numbers are given\n",
    "# as complex quantities.  use: ?np.mgrid\n",
    "Y, X =  np.mgrid[y0:y1:100j, x0:x1:100j]\n",
    "Y1, X1 =  np.mgrid[y0:y1:25j, x0:x1:25j]\n",
    "\n",
    "Z = np.sin(2*X)*np.sin(Y)\n",
    "\n",
    "# u and v here are the results of applying the gradient operation\n",
    "# to our scalar field.  Probably wise to check this in a seperate\n",
    "# code block.\n",
    "u = (2*np.sin(Y1)*np.cos(2*X1))\n",
    "v = (np.sin(2*X1)*np.cos(Y1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "filled_contour = plt.contourf(X,Y,Z,5)\n",
    "fig.colorbar(filled_contour, ax=ax)\n",
    "plt.quiver(X1,Y1,u,v, color='white')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Divergence\n",
    "----\n",
    "\n",
    "If we interpret the gradient field as a measurement of the flow of some quantity through space, then the divergence of that vector field is a measurement of the accumulation of that quantity.  To compute the divergence we compute the dot product of _nabla_ and the vector field.  The resulting quantity is a scalar quantity:\n",
    "\n",
    "$$\n",
    "\\overrightarrow{\\nabla} \\cdot \\overrightarrow{F} = \\frac{\\partial F_x}{\\partial x} + \\frac{\\partial F_y}{\\partial y} + \\frac{\\partial F_z}{\\partial z}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "### Accumulation and Diffusion\n",
    "----\n",
    "\n",
    "One of the key problems in science and engineering is the diffusion of an extensive quantity such as energy, charge, or mass.  The gradient of a potential provides a driving force for this diffusion.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Using the diffusion of mass as an example:\n",
    "\n",
    "$$\\overrightarrow{J} = -M \\overrightarrow{\\nabla} \\cdot {\\mu} $$\n",
    "\n",
    "and then computing the accumulation based on the vector field:\n",
    "\n",
    "$$\\frac{\\partial X(x,t)}{\\partial t} = - \\overrightarrow{\\nabla} \\cdot \\overrightarrow{J} $$\n",
    "\n",
    "where the minus sign indicates that accumulation occurs antiparallel to the gradient.  You may have seen "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Plotting the Flux (or Accumulation)\n",
    "----\n",
    "\n",
    "In this section we keep track of the vector components by hand using the mathematical definitions.  Writing this out we have a scalar field:\n",
    "\n",
    "$$\n",
    "F(x,y)\n",
    "$$\n",
    "\n",
    "the flux written as minus the gradient of the field (with $M=1$):\n",
    "\n",
    "$$\n",
    "\\overrightarrow{J} = -\\nabla F(x,y) = -\\left(\\frac{\\partial F}{\\partial x} \\hat{i} + \\frac{\\partial F}{\\partial y}\\hat{j} \\right)\n",
    "$$\n",
    "\n",
    "and the accumulation as minus the divergence of the flux:\n",
    "\n",
    "$$\n",
    "A = -\\nabla \\cdot \\overrightarrow{J} = \\left(\\frac{\\partial^2 F}{\\partial x^2} + \\frac{\\partial^2 F}{\\partial y^2}\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "x, y = sp.symbols('x y')\n",
    "\n",
    "# Note NEW field function:\n",
    "concentrationFunction = sp.sin(sp.pi*x)*sp.cos(sp.pi*y)\n",
    "\n",
    "\n",
    "fluxX = -sp.diff(concentrationFunction,x)\n",
    "fluxY = -sp.diff(concentrationFunction,y)\n",
    "accumulationFunction = sp.diff(concentrationFunction,x,2) + sp.diff(concentrationFunction,y,2)\n",
    "\n",
    "# We use lambdify to permit the functions to take arguments and vectorize the computations.\n",
    "myConcentrationFunction = sp.lambdify((x,y), concentrationFunction, 'numpy')\n",
    "myFluxX = sp.lambdify((x,y), fluxX, 'numpy')\n",
    "myFluxY = sp.lambdify((x,y), fluxY, 'numpy')\n",
    "myAccumulationFunction = sp.lambdify((x,y), accumulationFunction, 'numpy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "concentrationFunction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fluxX"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fluxY"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We use our quiver plotting capability to plot:\n",
    "\n",
    "* the scalar value for the concentration as `z`\n",
    "* the $\\hat{i}$ component of the flux as `u`\n",
    "* the $\\hat{j}$ component of the flux as `v`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "x0, x1 = (-1,1)\n",
    "y0, y1 = (-1,1)\n",
    "plotResolution = 200\n",
    "\n",
    "Y, X =  np.mgrid[y0:y1:200j, x0:x1:200j]\n",
    "# Quivers are on a seperate grid since they clutter things up.\n",
    "Y1, X1 =  np.mgrid[y0:y1:20j, x0:x1:20j]\n",
    "\n",
    "Z = myConcentrationFunction(X,Y)\n",
    "u = myFluxX(X1,Y1)\n",
    "v = myFluxY(X1,Y1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "from matplotlib import cm\n",
    "\n",
    "plt.contourf(X,Y,Z,20, cmap=cm.coolwarm)\n",
    "plt.colorbar()\n",
    "plt.quiver(X1,Y1,u,v, color='white')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "accumulationFunction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "In this example we use the contour plot to show locations of high (red) and low (blue) concentrations.  From this concentration field we compute the flux and show that as a quiver plot overlaid on the contour plot.  Examination of this figure confirms our physical understanding that mass flows from high to low concentration areas when the chemical potential is proportional to the concentration of species.  \n",
    "\n",
    "The actual RATE of accumulation at this particular INSTANT in time is given by the `accumulationFunction`.  This needs to be recomputed for every increment of time to be practical.  I'll show you how to do that in the upcoming lectures."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Homework:  Using the ideal solution model for the Gibbs free energy, *derive* the accumulation of species in a diffusive problem.\n",
    "\n",
    "The Gibbs free energy is an auxiliary function re-cast in variables that are more appropriate for experimental observation.\n",
    "\n",
    "$$\n",
    "G = G(T,P,n_i,...)\n",
    "$$\n",
    "\n",
    "The flux is the divergence of the chemical potential.  The chemical potential is the derivative of the Gibbs free energy (per mole, hence the switch to mole fraction) with respect to composition.  Review the chain rule before you start this problem.  The ideal solution model in terms of mole fraction for a two component system is:\n",
    "\n",
    "$$\n",
    "G(X_B) = (1-X_B)G_A + X_B G_B + RT(X_B \\ln X_B + (1-X_B) \\ln (1-X_B))\n",
    "$$\n",
    "\n",
    "The ideal solution model requires that the activity of a species be proportional to the mole fraction of that species in solution and that the heat of mixing be zero.\n",
    "\n",
    "Compute the divergence of the gradient of the derivative (w.r.t. concentration) of this function.  Note that the concentration is a function of spatial coordinates.  In this problem assume that the mole fraction $X$ is only a function of $x$, i.e. $X(x)$.  $G_A$ and $G_B$ are constants that depend on the melting temperature and heat of fusion, $R$ is the gas constant and $T$ is the absolute temperature."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#  Advanced and Optional Activities"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "There are multiple ways to interact with Python and get at the gradient of a function.  In the first instance we can use the coordinate system capabilities of `sympy` so that we can access the built in method `.gradient()`.  We start by defining a coordinate system and then calling gradient on our scalar function.  Scalars and vectors are objects of the coordinate system.  [See this page](http://docs.sympy.org/latest/modules/vector/intro.html) for more information on the vector module.\n",
    "\n",
    "For the purposes of numeric computing, `NumPy` has functions and operators that work as you would expect for vector operations.  What is presented below is for convenience, only."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import sympy as sp\n",
    "import sympy.vector as spv\n",
    "\n",
    "C = spv.CoordSys3D('C')\n",
    "spv.gradient(sp.sin(2*C.x)*sp.sin(C.y))\n",
    "# The gradient function should return something that looks like u+v\n",
    "# where u is a vector in the x direction and v is a vector in the\n",
    "# y direction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Define your example scalar field (a concentration like \n",
    "# variable C(x,y).\n",
    "exampleField = sp.sin(sp.pi*C.x)*sp.cos(sp.pi*C.y)\n",
    "exampleField"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We then use the built in `sympy` function `gradient` to compute the gradient:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Compute the gradient.\n",
    "gradientOfField = -spv.gradient(exampleField)\n",
    "gradientOfField"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We then compute the divergence of the gradient.  Note the absence of `C.i` and `C.j` in the answer indicating that these are not components of a vector.  (Compare this to the last slide.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Compute the divergence.\n",
    "accumulation = -spv.divergence(gradientOfField)\n",
    "accumulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Optional:  Read [this](http://blogs.nature.com/onyourwavelength/2018/07/18/heuristics-for-better-figures/) article and write a 200 word essay on choosing colors for heat maps."
   ]
  }
 ],
 "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}