Positivity of advection with backward Euler and different space discretizations.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import nodepy\n",
    "from nodepy import semidisc\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we investigate experimentally the positivity of some discretizations of the advection equation\n",
    "\n",
    "$$ U_t = a U_x.$$\n",
    "\n",
    "If our spatial discretization is $U_x \\approx \\frac{1}{\\Delta x} L u$ and we use the backward Euler method in time, we have\n",
    "\n",
    "$$u^{n+1} = u^n + \\frac{\\Delta t}{\\Delta x} L u^{n+1}$$\n",
    "\n",
    "or\n",
    "\n",
    "$$u^{n+1} = (I-\\nu L)^{-1} u^n,$$\n",
    "\n",
    "where $\\nu = \\Delta t/\\Delta x$ is the CFL number.  The method will be positivity preserving if the matrix\n",
    "$(I-\\nu L)^{-1}$ has non-negative entries."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Advection with centered differences"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After some experimentation, it seems in this case that there is no positivity preservation for even numbers of grid points.  With an odd number of grid points, positivity is obtained for large enough CFL numbers.  It appears that the \"large enough\" value is roughly proportional to the number of points in space.\n",
    "\n",
    "Also, with an even number of grid points, the minimum value approaches zero as the CFL number increases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "m = 57\n",
    "L = semidisc.centered_advection_matrix(m,1.)\n",
    "I = np.eye(m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0068230143449347844"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nu = 100.\n",
    "M = np.linalg.inv(I-nu*L)\n",
    "\n",
    "np.min(M)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Advection with Fourier spectral differencing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, the matrix always has negative values for CFL numbers smaller than 2 and always has only non-negative values for CFL numbers bigger than 3.  The exact cutoff point seems to increase very slowly with m, but is never bigger than about 2.5 for the range of sizes I tested."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def F_matrix(m):\n",
    "    F = np.zeros((m,m),dtype=complex)\n",
    "    for j in range(m):\n",
    "        v = np.zeros(m)\n",
    "        v[j] = 1.\n",
    "        F[:,j] = np.fft.fft(v)\n",
    "    return F"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "m = 21\n",
    "F = F_matrix(m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Finv = np.linalg.inv(F)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "length = 2*np.pi\n",
    "xi = np.fft.fftfreq(m)*m*2*np.pi/length\n",
    "D = np.diag(1.j*xi)\n",
    "L = np.dot(Finv,np.dot(D,F))\n",
    "I = np.eye(m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.0050201735752721842"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nu = 2.02\n",
    "M = np.linalg.inv(I-nu*L)\n",
    "M = np.real(M)\n",
    "\n",
    "np.min(M)"
   ]
  }
 ],
 "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}