{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-cell"
    ]
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.linalg import solve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solution strategy for channel flow\n",
    "\n",
    "## Approach\n",
    "\n",
    "Rearrange the discretized equation so that all unknowns are on the left-hand side, and only known values are on the right-hand side.\n",
    "\n",
    "\\begin{equation}\n",
    "  v_{i-1} - 2v_i + v_{i+1} = \\frac{dP/dx}{\\eta}\\Delta y^2\n",
    "\\end{equation}\n",
    "\n",
    "There is one such equation for each internal grid point. At the boundaries, boundary conditions have to be used. For example, $v_0 = v_\\mathrm{surf}$ and $v_{\\mathrm{ny}-1} = v_\\mathrm{bott}$.\n",
    "\n",
    "For five grid points we thus have five equations in total:\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{matrix}\n",
    "i=0:\\qquad & v_0 &    &      &    &     &    &     &   &     & = & v_\\mathrm{surf} \\newline\n",
    "i=1:\\qquad & v_0 & -2 & v_1  & +  & v_2 &    &     &   &     & = & \\frac{dP/dx}{\\eta}\\Delta y^2 \\newline\n",
    "i=2:\\qquad &     &    & v_1  & -2 & v_2 & +  & v_3 &   &     & = & \\frac{dP/dx}{\\eta}\\Delta y^2 \\newline\n",
    "i=3:\\qquad &     &    &      &    & v_2 & -2 & v_3 & + & v_4 & = & \\frac{dP/dx}{\\eta}\\Delta y^2 \\newline\n",
    "i=4:\\qquad &     &    &      &    &     &    &     &   & v_4 & = & v_\\mathrm{bott}\n",
    "\\end{matrix}\n",
    "\\end{equation}\n",
    "\n",
    "This system of equations can be written in matrix format as\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{pmatrix}\n",
    "1 &    &    &    &   \\newline\n",
    "1 & -2 &  1 &    &   \\newline \n",
    "  &  1 & -2 &  1 &   \\newline\n",
    "  &    &  1 & -2 & 1 \\newline\n",
    "  &    &    &    & 1\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "v_0 \\newline\n",
    "v_1 \\newline\n",
    "v_2 \\newline\n",
    "v_3 \\newline\n",
    "v_4\n",
    "\\end{pmatrix} =\n",
    "\\begin{pmatrix}\n",
    "v_\\mathrm{surf} \\newline\n",
    "\\frac{dP/dx}{\\eta}\\Delta y^2 \\newline\n",
    "\\frac{dP/dx}{\\eta}\\Delta y^2 \\newline\n",
    "\\frac{dP/dx}{\\eta}\\Delta y^2 \\newline\n",
    "v_\\mathrm{bott}\n",
    "\\end{pmatrix}\n",
    "\\end{equation}\n",
    "\n",
    "and is typically written with symbols\n",
    "\n",
    "\\begin{equation}\n",
    "  \\mathbf{Ax}=\\mathbf{b}\n",
    "\\end{equation}\n",
    "\n",
    "Python comes with a module, `scipy.linalg` that includes function `solve`, which can be used to solve this type of systems of linear equations:\n",
    "\n",
    "```python\n",
    "from scipy.linalg import solve\n",
    "x = solve(A, b)\n",
    "```\n",
    "  \n",
    "where `b` has to be a 1-D NumPy array of size `ny`, and includes the numerical values of the right hand side expressions; and `A` a 2-D NumPy array of size (`ny`,`ny`), and includes the numerical coefficients, as above. `x` will then contain the wanted velocity values. The following script utilizes this function to solve the channel flow problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise - Exploring a channel flow\n",
    "\n",
    "- Complete the following script\n",
    "- Once complete, answer the following:\n",
    "    - Which direction is the flow mostly going?\n",
    "    - How can you make it go more to the left or the right?\n",
    "    - What happens if you decrease the viscosity? Why?\n",
    "    - What happens if you remove (or comment out) the lines where B.C. are set?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "raises-exception"
    ]
   },
   "outputs": [],
   "source": [
    "# Define physical parameters\n",
    "visc = 1e19                    # Viscosity of the rock in the channel, Pa s\n",
    "dpdx = -2800*9.81*4000/500e3   # Pressure gradient within the channel, Pa m^-1\n",
    "\n",
    "# Define problem geometry\n",
    "ny = 5      # number of grid points\n",
    "L = 10e3    # size of the problem (width of the channel)\n",
    "\n",
    "# Define values for boundary conditions, m s^-1\n",
    "vx_surf = 0\n",
    "vx_bott = -0.01/(60*60*24*365.25)\n",
    "\n",
    "# Calculate grid spacing and y coordinates\n",
    "dy = L / (ny-1)\n",
    "y = np.linspace(0, L, ny)\n",
    "\n",
    "# Create the empty (zero) coefficient and right hand side arrays\n",
    "A = np.zeros((ny,ny))  # 2-dimensional array, ny rows, ny columns\n",
    "b = np.zeros(ny)\n",
    "\n",
    "# Set B.C. values in the coefficient array and in the r.h.s. array\n",
    "A[0, 0] = 1\n",
    "b[0] = vx_surf\n",
    "A[ny-1, ny-1] = 1\n",
    "b[ny-1] = vx_bott\n",
    "\n",
    "# Set remaining (internal grid point) coefficients in the array `A`.\n",
    "# TODO: Complete the for loop so that it will write the coefficient \n",
    "# values in the array `A`. The for loop loops over all the rows\n",
    "# of the matrix. On each row, you need to fill in three coefficients.\n",
    "# In python, elements of 2D arrays are referenced to like 'A[row, col]'\n",
    "for iy in range(1, ny-1):\n",
    "    .\n",
    "    .\n",
    "    .\n",
    "    \n",
    "\n",
    "# Set remaining (internal grid point) coefficients in the r.h.s. array b\n",
    "# TODO: Create a for loop that loops over the internal grid points (i.e.\n",
    "# from 1 to ny-1), and fills in the `b` array with r.h.s. values\n",
    ".\n",
    ".\n",
    ".\n",
    "\n",
    "# For debugging purposes, plot A (if it is small enough) and b\n",
    "if ny < 15:\n",
    "    print(\"A is:\\n\", A, \"\\nb is:\\n\", b.T)\n",
    "    \n",
    "# Solve it!\n",
    "vx = solve(A, b)\n",
    "\n",
    "# Plot it!\n",
    "plt.plot(vx * (60*60*25*365.25), -y)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating strain rates\n",
    "\n",
    "The strain rate in 2D is defined as \n",
    "\n",
    "\\begin{equation}\n",
    "  \\underline{\\underline{\\dot\\epsilon}} = \\begin{pmatrix}\n",
    "  \\frac{\\partial v_x}{\\partial x} & \\frac{\\partial v_x}{\\partial y} \\newline\n",
    "  \\frac{\\partial v_y}{\\partial x} & \\frac{\\partial v_y}{\\partial y} \\newline\n",
    "  \\end{pmatrix}\n",
    "\\end{equation}\n",
    "\n",
    "In our case the only relevant component is the shear strain rate $\\dot\\epsilon_{xy} = \\frac{\\partial v_x}{\\partial y}$ . Why?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise - Finite difference strain rates\n",
    " \n",
    "- Discretize the shear strain rate expression\n",
    "- Complete the script below so that it will loop over all the grid points, calculate the shear strain, and store it in the array `e`\n",
    "- How does the resulting plot look like? Note that the *shear stress* is defined as $\\sigma_{xy} = \\eta\\dot\\epsilon_{xy}$. Where is the shear stress smallest?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "raises-exception"
    ]
   },
   "outputs": [],
   "source": [
    "e = np.zeros(ny-1)\n",
    "\n",
    "# TODO: Write here a for-loop that loops over the grid points\n",
    "# and calculate the shear strain rate on every step\n",
    ".\n",
    ".\n",
    ".\n",
    " \n",
    "plt.plot(e, -(y[1:] + y[:-1]) / 2.0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Velocity boundary conditions\n",
    "\n",
    "Like in the case of heat equation, the FD formulation for velocity field may have different boundary conditions:\n",
    "\n",
    "1. No-slip:\n",
    "    - Velocity parallel to the surface is fixed at zero\n",
    "    - $v_x = 0$\n",
    "2. Free-slip: \n",
    "    - Shear stress at the surface is zero\n",
    "    - $0 = \\sigma_{xy} = \\eta\\dot\\epsilon_{xy} = \\eta\\frac{\\partial v_x}{\\partial y} \\Rightarrow \\frac{\\partial v_x}{\\partial y} = 0$\n",
    "3. Other:\n",
    "    - Velocity is fixed at non-zero value\n",
    "    - For example, $v_x = a \\neq 0$, also incoming/outgoing flow ($v_y \\neq 0$) possible (in 2D models)\n",
    " \n",
    "We used cases 1 and 3 previously."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise - Exploring channel flow boundary conditions\n",
    " \n",
    "- Discretize the free-slip boundary condition $\\frac{\\partial v_x}{\\partial y} = 0$\n",
    "- If you want to apply this at the upper surface (i.e. at $x=x_0$), how would your system of equations change? How would the coefficients in the matrix $A$ change?\n",
    "- Modify your script for the channel flow model, and change the upper surface velocity boundary condition from $v_x=0$ to a free-slip boundary condition\n",
    "    - How does the resulting velocity plot change?\n",
    "    - How does the plot for shear strain rate change?"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "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.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}