{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "be12cf16-dfa2-4a42-8e4e-a5bfafb12389",
   "metadata": {},
   "source": [
    "### THE DANGER OF EXCESSIVELY LARGE SINGLE TIME STEPS IN DIFFUSION\n",
    "\n",
    "->  When the time step is (0.5 / diffusion rate),\n",
    "    a 2-bin system equilibrates in a single step,\n",
    "    and some 3-bin systems can over-shoot equilibrium!\n",
    "\n",
    "->  When the time step is (0.33333 / diffusion rate),\n",
    "    some 3-bin systems equilibrate in a single step\n",
    "\n",
    "So, (0.33333 / diffusion rate) is a - rather lax - upper bound for\n",
    "sensible single time steps!  \n",
    "\n",
    "IMPORTANT: The above is for delta_x = 1; in general, multiply by delta_x**2\n",
    "\n",
    "The **\"Von Neumann stability analysis\"**, which provides a slighly looser max value of time steps, is also discussed.\n",
    "\n",
    "That value of 0.33333 is saved in the Class variable \"time_step_threshold\"\n",
    "\n",
    "LAST REVISED: June 23, 2024 (using v. 1.0 beta34.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "fcdcc8cf-6b07-4e37-9c80-8e261e0e9def",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Added 'D:\\Docs\\- MY CODE\\BioSimulations\\life123-Win7' to sys.path\n"
     ]
    }
   ],
   "source": [
    "import set_path      # Importing this module will add the project's home directory to sys.path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "fe25e919",
   "metadata": {},
   "outputs": [],
   "source": [
    "from life123 import ChemData as chem\n",
    "from life123 import BioSim1D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "61eefdac-4369-45d4-aeb4-513ef0a9c19f",
   "metadata": {},
   "outputs": [],
   "source": [
    "chem_data = chem(diffusion_rates=[10.])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bacef70f-15ac-42f7-b2b1-d74362919af7",
   "metadata": {},
   "source": [
    "## Simulate a 2-bin system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "992cf517-dde9-4341-8640-492aaffcf986",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SYSTEM STATE at Time t = 0:\n",
      "2 bins and 1 species:\n",
      "  Species 0 (Chemical 1). Diff rate: 10.0. Conc: [100.   0.]\n"
     ]
    }
   ],
   "source": [
    "bio = BioSim1D(n_bins=2, chem_data=chem_data)\n",
    "bio.inject_conc_to_bin(bin_address=0, delta_conc=100., species_index=0)\n",
    "bio.describe_state()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c29b3617-7913-4cf2-be84-7e8def8e4476",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[50. 50.]]\n"
     ]
    }
   ],
   "source": [
    "bio.time_step_threshold = 0.51    # To bypass the safety limit (for maximum delta Time) that is typically in place\n",
    "\n",
    "# When the time step is (0.5 / diffusion rate),\n",
    "# a 2-bin system equilibrates in a single step!\n",
    "bio.diffuse(time_step=0.05, n_steps=1)\n",
    "print(bio.system) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d9d808c-a64e-4e75-a173-e2f46918712b",
   "metadata": {},
   "source": [
    "_Note: [[50. 50.]]  : **the two bins have equilibrated!**_"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "104fce42-1ae2-456e-8160-73f48c2040b2",
   "metadata": {},
   "source": [
    "## Start over with a 3-bin system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "5c1bd9a7-9323-4b8a-9130-40dba22dc2ac",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SYSTEM STATE at Time t = 0:\n",
      "3 bins and 1 species:\n",
      "  Species 0 (Chemical 1). Diff rate: 10.0. Conc: [  0. 100.   0.]\n"
     ]
    }
   ],
   "source": [
    "bio = BioSim1D(n_bins=3, chem_data=chem_data)\n",
    "bio.inject_conc_to_bin(bin_address=1, delta_conc=100., species_index=0)\n",
    "bio.describe_state()\n",
    "#3 bins and 1 species:   [[  0. 100.   0.]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "95397550-9261-494d-8c17-b410a2474b61",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[50.  0. 50.]]\n"
     ]
    }
   ],
   "source": [
    "bio.time_step_threshold = 0.51    # To bypass the safety limit (for maximum delta Time) that is typically in place\n",
    "\n",
    "# When the time step is (0.5 / diffusion rate),\n",
    "# a 3-bin system can overshoot equilibrium!\n",
    "bio.diffuse(time_step=0.05, n_steps=1)\n",
    "print(bio.system)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "718a6239-aa63-4574-9ef2-2b06a934844c",
   "metadata": {},
   "source": [
    "_Note: [[50.  0. 50.]] : the diffusion has **over-shot equilibrium!!!**_\n",
    "#### Cleary, a very excessively large time step!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bab05738-3e8d-4c8c-9c5a-cad9dc26ac68",
   "metadata": {},
   "source": [
    "## Start over again with a new 3-bin system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "3ebadd9c-3604-43d0-abe8-f33fdc9b552f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SYSTEM STATE at Time t = 0:\n",
      "3 bins and 1 species:\n",
      "  Species 0 (Chemical 1). Diff rate: 10.0. Conc: [  0. 100.   0.]\n"
     ]
    }
   ],
   "source": [
    "bio = BioSim1D(n_bins=3, chem_data=chem_data)\n",
    "bio.inject_conc_to_bin(bin_address=1, delta_conc=100., species_index=0)\n",
    "bio.describe_state()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "0fd086b3-8b8f-48cf-958f-50799c130bd3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[33.333 33.334 33.333]]\n"
     ]
    }
   ],
   "source": [
    "bio.time_step_threshold = 0.34    # To bypass the safety limit (for maximum delta Time) that is typically in place\n",
    "\n",
    "\n",
    "bio.diffuse(time_step=0.033333, n_steps=1)\n",
    "print(bio.system)  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df3c48e7-7492-499c-a4ac-d3a2f1e9dcaa",
   "metadata": {},
   "source": [
    "#### When the time step is (0.33333 / diffusion rate),\n",
    "#### a 3-bin system, configured as above, equilibrates in *a single step!*\n",
    "#### Again, an excessively large time step!  (note that we're using the default value of 1 for delta_x)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25d75d27-675a-458b-a465-4b94b610a627",
   "metadata": {},
   "source": [
    "## Start over again with a new 3-bin system : same as the last round, but this time use a delta_x = 10 (rather than the default 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "8659eb72-6304-4512-b452-6acd84cc23af",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SYSTEM STATE at Time t = 0:\n",
      "3 bins and 1 species:\n",
      "  Species 0 (Chemical 1). Diff rate: 10.0. Conc: [  0. 100.   0.]\n"
     ]
    }
   ],
   "source": [
    "bio = BioSim1D(n_bins=3, chem_data=chem_data)\n",
    "bio.inject_conc_to_bin(bin_address=1, delta_conc=100., species_index=0)\n",
    "bio.describe_state()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "e55a0df6-0442-44d1-a749-cbb356c53e58",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "bio.time_step_threshold = 1000   # To bypass the safety limit (for maximum delta Time) that is typically in place\n",
    "\n",
    "# When the time step is (0.33333 / diffusion rate),\n",
    "# a 3-bin system, configured as above, equilibrates in a single step!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "a93fc736-b8a7-4442-9ff5-cf9d622e7259",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[33.333 33.334 33.333]]\n"
     ]
    }
   ],
   "source": [
    "bio.diffuse(time_step=3.3333, n_steps=1, delta_x=10)\n",
    "print(bio.system)  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d369711-0001-41b4-bbf0-f0b8964317d1",
   "metadata": {},
   "source": [
    "### This generalizes the previous result to any arbitrary delta_x :\n",
    "#### in this scenario, a time step of (delta_x**2 * 0.33333 / diffusion rate)\n",
    "#### lead to equilibrium in a single step!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78d68b19-4031-421a-b030-186ede007bef",
   "metadata": {},
   "source": [
    "# The physical intuition explored in this experiment suggests enforcing:\n",
    "#### delta_t < delta_x**2 * 0.33333 / diffusion rate\n",
    "\n",
    "#### Interestingly, this is only slightly stricter than the upper bound that emerges from the **\"Von Neumann stability analysis\"** of the diffusion equation in 1-D, which states that solutions may become unstable unless\n",
    "delta_t < delta_x**2 * 0.5 / diffusion rate\n",
    "\n",
    "## NOTE: The `diffuse()` method enforces the stricter upper bound"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91ff8434-597e-4fd6-b76f-7f3572732863",
   "metadata": {},
   "source": [
    "#### To keep in mind that the _\"Von Neumann stability analysis\"_ is ONLY applicable when the diffusion equation is being solved with the **\"explicit Forward-Time Centered Space\"** method, which is the one currently used by the `diffuse()` method.  \n",
    "An explanation can be found at: https://www.youtube.com/watch?v=QUiUGNwNNmo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08fb1f3a",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}