465 lines
40 KiB
Plaintext
465 lines
40 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 12,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# imports\n",
|
|
"import numpy as np\n",
|
|
"import matplotlib.pyplot as plt\n",
|
|
"\n",
|
|
"%matplotlib notebook\n",
|
|
"%matplotlib inline\n",
|
|
"from ipywidgets import interact, IntSlider"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 13,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# set up environment\n",
|
|
"dim = 2\n",
|
|
"grid_size = {'x': 20, 'y': 20} # to start out\n",
|
|
"domain_size = {'x': 20, 'y': 20} # to start out\n",
|
|
"\n",
|
|
"# boundary_condition\n",
|
|
"# constant boundary condition\n",
|
|
"bc_left = 0\n",
|
|
"bc_right = 0\n",
|
|
"bc_top = 0\n",
|
|
"bc_bottom = 0\n",
|
|
"\n",
|
|
"delta_x = domain_size['x'] / grid_size['x']\n",
|
|
"delta_y = domain_size['y'] / grid_size['y'] \n",
|
|
"\n",
|
|
"alpha_x = np.ones((grid_size['x'],grid_size['y']))\n",
|
|
"alpha_y = np.ones((grid_size['x'],grid_size['y']))\n",
|
|
"\n",
|
|
"# time dimension\n",
|
|
"max_stable_time_step_x = delta_x**2 / (2 * np.max(alpha_x))\n",
|
|
"max_stable_time_step_y = delta_y**2 / (2 * np.max(alpha_y))\n",
|
|
"max_stable_time_step = max_stable_time_step_x if max_stable_time_step_x <= max_stable_time_step_y else max_stable_time_step_y\n",
|
|
"\n",
|
|
"time_step = 0.1\n",
|
|
"if time_step > max_stable_time_step:\n",
|
|
" print('Warning: Time step more than maximum stable time step!')"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 29,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAQUAAAD4CAYAAADl7fPiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAANR0lEQVR4nO3df6hf9X3H8edrSbSQ6dQ5U39krXRBiGVkJcSVuaGztVFkaUfpImPLNkFXJqwwGG6DWrp/OoaTDcWSdkE7WrVsyxpo/BGygRVaa5T4q9WZSYq5TZO16Uxdu7roe3/ck3I/N9+b3H1/3O/33jwfcPme8/l8vud8Dl945Zzz/ea8U1VI0nE/Ne4JSJoshoKkhqEgqWEoSGoYCpIay8c9gV7OyJn1NlaOexrSkvU//Ddv1I/Tq28iQ+FtrOSKXDPuaUhL1hO1e84+Lx8kNQYKhSQbk7yUZF+S23r0n5nkwa7/iSTvHGR/kkav71BIsgy4G7gOWAvcmGTtrGE3Ad+vql8A7gT+qt/9SVoYg5wpbAD2VdUrVfUG8ACwadaYTcB93fI/Atck6XlzQ9JkGCQULgZenbF+oGvrOaaqjgGvAT/ba2NJbk6yJ8me/+XHA0xL0iAm5kZjVW2tqvVVtX4FZ457OtJpa5BQmAJWz1i/pGvrOSbJcuBngO8NsE9JIzZIKDwJrElyaZIzgM3AjlljdgBbuuUPA/9a/l9taaL1/eOlqjqW5FbgEWAZsK2qXkjySWBPVe0A/h74hyT7gCNMB4ekCZZJ/If77JxX/qJRGp0najdH60jPbwIn5kajpMlgKEhqGAqSGoaCpIahIKlhKEhqGAqSGoaCpIahIKlhKEhqGAqSGoaCpIahIKlhKEhqGAqSGoaCpIahIKlhKEhqDFIhanWSf0vyjSQvJPnjHmOuSvJakr3d38cHm66kURuk6vQx4E+q6ukkZwFPJdlVVd+YNe4rVXXDAPuRtID6PlOoqoNV9XS3/APgm5xYIUrSIjOUewpdNelfAp7o0f3eJM8keSjJ5SfZhmXjpAkwyOUDAEl+Gvgn4GNVdXRW99PAO6rq9STXA/8CrOm1naraCmyF6Ue8DzovSf0Z6EwhyQqmA+HzVfXPs/ur6mhVvd4t7wRWJDl/kH1KGq1Bvn0I0xWgvllVfzPHmLcfLz2fZEO3P2tJShNskMuHXwF+B3guyd6u7c+Bnweoqk8zXT/yo0mOAT8CNltLUppsg9SSfBzoWXZqxpi7gLv63YekhecvGiU1DAVJDUNBUsNQkNQwFCQ1DAVJDUNBUsNQkNQwFCQ1DAVJDUNBUsNQkNQwFCQ1DAVJDUNBUmPgZzQKHvn23nmN+8BF60Y6D2kYPFOQ1DAUJDUGDoUk+5M815WF29OjP0n+Lsm+JM8mec+g+5Q0OsO6p3B1VX13jr7rmK71sAa4Arine5U0gRbi8mET8Lma9jXgnCQXLsB+JfVhGKFQwKNJnkpyc4/+i4FXZ6wfoEfNScvGSZNhGJcPV1bVVJILgF1JXqyqx/6/G7FsnDQZBj5TqKqp7vUwsB3YMGvIFLB6xvolXZukCTRoLcmVSc46vgxcCzw/a9gO4He7byF+GXitqg4Osl9JozPo5cMqYHtXLnI58IWqejjJH8JPSsftBK4H9gE/BH5/wH1KGqFMYmnHs3NeXZFrxj0Nacl6onZztI70LPvoLxolNQwFSQ1DQVLDUJDUMBQkNQwFSQ1DQVLDUJDUMBQkNQwFSQ1DQVLDUJDUMBQkNQwFSQ1DQVLDUJDUMBQkNQwFSY2+QyHJZV2puON/R5N8bNaYq5K8NmPMxweesaSR6vvBrVX1ErAOIMkyph/bvr3H0K9U1Q397kfSwhrW5cM1wH9U1beGtD1JYzKsUNgM3D9H33uTPJPkoSSXz7UBy8ZJk2HgR7wnOQP4NnB5VR2a1Xc28FZVvZ7keuBvq2rNqbbpI96l0Rr1I96vA56eHQgAVXW0ql7vlncCK5KcP4R9ShqRYYTCjcxx6ZDk7enKRyXZ0O3ve0PYp6QRGahsXFc/8v3ALTPaZpaM+zDw0STHgB8Bm2sSS1JJ+gnLxkmnIcvGSZo3Q0FSw1CQ1DAUJDUMBUkNQ0FSw1CQ1DAUJDUMBUkNQ0FSw1CQ1DAUJDUMBUkNQ0FSw1CQ1DAUJDUMBUkNQ0FSY16hkGRbksNJnp/Rdl6SXUle7l7PneO9W7oxLyfZMqyJSxqN+Z4p3AtsnNV2G7C7q+Owu1tvJDkPuB24AtgA3D5XeEiaDPMKhap6DDgyq3kTcF+3fB/wwR5v/QCwq6qOVNX3gV2cGC6SJsgg9xRWVdXBbvk7wKoeYy4GXp2xfqBrkzShhnKjsavlMNCz4q0lKU2GQULhUJILAbrXwz3GTAGrZ6xf0rWdoKq2VtX6qlq/gjMHmJakQQwSCjuA498mbAG+1GPMI8C1Sc7tbjBe27VJmlDz/UryfuCrwGVJDiS5CfgU8P4kLwPv69ZJsj7JZwGq6gjwl8CT3d8nuzZJE8qycdJpyLJxkubNUJDUMBQkNQwFSQ1DQVLDUJDUMBQkNQwFSQ1DQVLDUJDUMBQkNQwFSQ1DQVLDUJDUMBQkNQwFSQ1DQVLDUJDUMBQkNU4ZCnPUkfzrJC8meTbJ9iTnzPHe/UmeS7I3yZ4hzlvSiMznTOFeTiz1tgt4d1X9IvDvwJ+d5P1XV9W6qlrf3xQlLaRThkKvOpJV9WhVHetWv8Z0kRdJS8Aw7in8AfDQHH0FPJrkqSQ3n2wjlo2TJsPyQd6c5C+AY8Dn5xhyZVVNJbkA2JXkxe7M4wRVtRXYCtN1HwaZl6T+9X2mkOT3gBuA3645KspU1VT3ehjYDmzod3+SFkZfoZBkI/CnwG9U1Q/nGLMyyVnHl5muI/l8r7GSJsd8vpLsVUfyLuAspi8J9ib5dDf2oiQ7u7euAh5P8gzwdeDLVfXwSI5C0tBYS1I6DVlLUtK8GQqSGoaCpIahIKlhKEhqGAqSGoaCpIahIKlhKEhqGAqSGoaCpIahIKlhKEhqGAqSGoaCpIahIKlhKEhqGAqSGv2WjftEkqnu+Yx7k1w/x3s3Jnkpyb4ktw1z4pJGo9+ycQB3duXg1lXVztmdSZYBdwPXAWuBG5OsHWSykkavr7Jx87QB2FdVr1TVG8ADwKY+tiNpAQ1yT+HWrur0tiTn9ui/GHh1xvqBrq0ny8ZJk6HfULgHeBewDjgI3DHoRKpqa1Wtr6r1Kzhz0M1J6lNfoVBVh6rqzap6C/gMvcvBTQGrZ6xf0rVJmmD9lo27cMbqh+hdDu5JYE2SS5OcAWwGdvSzP0kL55RVp7uycVcB5yc5ANwOXJVkHdOl5vcDt3RjLwI+W1XXV9WxJLcCjwDLgG1V9cIoDkLS8Fg2TjoNWTZO0rwZCpIahoKkhqEgqWEoSGoYCpIahoKkhqEgqWEoSGoYCpIahoKkhqEgqWEoSGoYCpIahoKkhqEgqWEoSGoYCpIa83lG4zbgBuBwVb27a3sQuKwbcg7wX1W1rsd79wM/AN4EjlXV+qHMWtLInDIUmC4bdxfwueMNVfVbx5eT3AG8dpL3X11V3+13gpIW1ilDoaoeS/LOXn1JAnwE+PUhz0vSmAx6T+FXgUNV9fIc/QU8muSpJDefbEOWjZMmw3wuH07mRuD+k/RfWVVTSS4AdiV5sStYe4Kq2gpshelHvA84L0l96vtMIcly4DeBB+caU1VT3ethYDu9y8tJmiCDXD68D3ixqg706kyyMslZx5eBa+ldXk7SBDllKHRl474KXJbkQJKbuq7NzLp0SHJRkp3d6irg8STPAF8HvlxVDw9v6pJGwbJx0mnIsnGS5s1QkNQwFCQ1DAVJDUNBUsNQkNQwFCQ1DAVJDUNBUsNQkNQwFCQ1DAVJDUNBUsNQkNQwFCQ1DAVJjYl8yEqS/wS+Nav5fGAp1o9YqscFS/fYlsJxvaOqfq5Xx0SGQi9J9izFClNL9bhg6R7bUj2u47x8kNQwFCQ1FlMobB33BEZkqR4XLN1jW6rHBSyiewqSFsZiOlOQtAAMBUmNRREKSTYmeSnJviS3jXs+w5Jkf5LnkuxNsmfc8xlEkm1JDid5fkbbeUl2JXm5ez13nHPsxxzH9YkkU93ntjfJ9eOc47BNfCgkWQbcDVwHrAVuTLJ2vLMaqqurat0S+N77XmDjrLbbgN1VtQbY3a0vNvdy4nEB3Nl9buuqameP/kVr4kOB6UrV+6rqlap6A3gA2DTmOWmWqnoMODKreRNwX7d8H/DBhZzTMMxxXEvaYgiFi4FXZ6wf6NqWggIeTfJUkpvHPZkRWFVVB7vl7zBddHipuDXJs93lxaK7LDqZxRAKS9mVVfUepi+N/ijJr417QqNS0999L5Xvv+8B3gWsAw4Cd4x1NkO2GEJhClg9Y/2Srm3Rq6qp7vUwsJ3pS6Wl5FCSCwG618Njns9QVNWhqnqzqt4CPsMS+9wWQyg8CaxJcmmSM4DNwI4xz2lgSVYmOev4MnAt8PzJ37Xo7AC2dMtbgC+NcS5DczzoOh9iiX1uy8c9gVOpqmNJbgUeAZYB26rqhTFPaxhWAduTwPTn8IWqeni8U+pfkvuBq4DzkxwAbgc+BXwxyU1M/1f4j4xvhv2Z47iuSrKO6cuh/cAt45rfKPgzZ0mNxXD5IGkBGQqSGoaCpIahIKlhKEhqGAqSGoaCpMb/AZlpw/ulk6QYAAAAAElFTkSuQmCC",
|
|
"text/plain": [
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"# initialize environment with start values\n",
|
|
"C_t = np.zeros((grid_size['x'],grid_size['y']))\n",
|
|
"C_t[5,4] = 2000\n",
|
|
"\n",
|
|
"\n",
|
|
"plt.imshow(C_t, vmin=0, vmax=2000)\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 30,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# calculate mean between two cells\n",
|
|
"def alpha_interblock(alpha1, alpha2, harmonic=False):\n",
|
|
" if not harmonic:\n",
|
|
" return 0.5 * (alpha1 + alpha2)\n",
|
|
" else:\n",
|
|
" return 2 / ((1/alpha1) + (1/alpha2))\n"
|
|
]
|
|
},
|
|
{
|
|
"attachments": {},
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"### Implementation for constant boundary conditions"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 31,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# simulate one time step\n",
|
|
"def simulate_constant_boundary(C_t): \n",
|
|
" C_t1 = np.zeros((grid_size['x'],grid_size['y']))\n",
|
|
"\n",
|
|
" # inner cells\n",
|
|
" for i in range(1, grid_size['y']-1): #rows\n",
|
|
" for j in range(1, grid_size['x']-1): #columns\n",
|
|
" C_t1[i,j] = C_t[i,j] \\\n",
|
|
" + time_step/delta_y**2 * (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) * C_t[i+1,j]\n",
|
|
" - (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j])) * C_t[i,j] \n",
|
|
" + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j]) * C_t[i-1,j]) \\\n",
|
|
" + time_step/delta_x**2 * (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) * C_t[i,j+1]\n",
|
|
" - (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j])) * C_t[i,j] \n",
|
|
" + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j]) * C_t[i,j-1])\n",
|
|
" \n",
|
|
" # boundary conditions\n",
|
|
" # left without corners / looping over rows\n",
|
|
" for i in range(1, grid_size['y']-1):\n",
|
|
" j = 0\n",
|
|
" C_t1[i,0] = C_t[i,0] \\\n",
|
|
" + time_step/delta_x**2 * (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) * C_t[i,j+1] \n",
|
|
" - (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) + 2 * alpha_x[i,j]) * C_t[i,j]\n",
|
|
" + 2 * alpha_x[i,j] * bc_left) \\\n",
|
|
" + time_step/delta_y**2 * (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) * C_t[i+1,j]\n",
|
|
" - (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j])) * C_t[i,j] \n",
|
|
" + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j]) * C_t[i-1,j]) \n",
|
|
"\n",
|
|
" # right without corners / looping over rows\n",
|
|
" n = grid_size['x']-1 # maximum index in x-direction (columns)\n",
|
|
" for i in range(1, grid_size['y']-1):\n",
|
|
" j = n\n",
|
|
" C_t1[i,n] = C_t[i,n] \\\n",
|
|
" + time_step/delta_x**2 * (2 * alpha_x[i,n] * bc_right \n",
|
|
" - (alpha_interblock(alpha_x[i,n-1], alpha_x[i,n]) + 2 * alpha_x[i,n]) * C_t[i,n]\n",
|
|
" + alpha_interblock(alpha_x[i,n-1], alpha_x[i,n]) * C_t[i,n-1]) \\\n",
|
|
" + time_step/delta_y**2 * (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) * C_t[i+1,j]\n",
|
|
" - (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j])) * C_t[i,j] \n",
|
|
" + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j]) * C_t[i-1,j])\n",
|
|
"\n",
|
|
" # top without corners / looping over columns\n",
|
|
" for j in range(1, grid_size['x']-1):\n",
|
|
" i = 0\n",
|
|
" C_t1[0,j] = C_t[0,j] \\\n",
|
|
" + time_step/delta_y**2 * (alpha_interblock(alpha_y[1,j], alpha_y[0,j]) * C_t[1,j] \n",
|
|
" - (alpha_interblock(alpha_y[1,j], alpha_y[0,j]) + 2 * alpha_y[0,j]) * C_t[0,j]\n",
|
|
" + 2 * alpha_y[0,j] * bc_top) \\\n",
|
|
" + time_step/delta_x**2 * (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) * C_t[i,j+1]\n",
|
|
" - (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j])) * C_t[i,j] \n",
|
|
" + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j]) * C_t[i,j-1])\n",
|
|
"\n",
|
|
" # bottom without corners / looping over columns\n",
|
|
" m = grid_size['y']-1 # maximum index in y-direction (rows)\n",
|
|
" for j in range(1, grid_size['x']-1):\n",
|
|
" i = m\n",
|
|
" C_t1[m,j] = C_t[m,j] \\\n",
|
|
" + time_step/delta_y**2 * (2 * alpha_y[m,j] * bc_bottom \n",
|
|
" - (alpha_interblock(alpha_y[m,j], alpha_y[m-1,j]) + 2 * alpha_y[m,j]) * C_t[m,j]\n",
|
|
" + alpha_interblock(alpha_y[m,j], alpha_y[m-1,j]) * C_t[m-1,j]) \\\n",
|
|
" + time_step/delta_x**2 * (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) * C_t[i,j+1]\n",
|
|
" - (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j])) * C_t[i,j] \n",
|
|
" + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j]) * C_t[i,j-1])\n",
|
|
"\n",
|
|
" # corner top left\n",
|
|
" i = 0\n",
|
|
" j = i\n",
|
|
" C_t1[i,j] = C_t[i,j] \\\n",
|
|
" + time_step/delta_x**2 * (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) * C_t[i,j+1] \n",
|
|
" - (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) + 2 * alpha_x[i,j]) * C_t[i,j]\n",
|
|
" + 2 * alpha_x[i,j] * bc_left) \\\n",
|
|
" + time_step/delta_y**2 * (alpha_interblock(alpha_y[1,j], alpha_y[0,j]) * C_t[1,j] \n",
|
|
" - (alpha_interblock(alpha_y[1,j], alpha_y[0,j]) + 2 * alpha_y[0,j]) * C_t[0,j]\n",
|
|
" + 2 * alpha_y[0,j] * bc_top)\n",
|
|
" \n",
|
|
"\n",
|
|
" # corner top right\n",
|
|
" i = 0\n",
|
|
" j = grid_size['x']-1\n",
|
|
" n = j\n",
|
|
" C_t1[i,j] = C_t[i,j] \\\n",
|
|
" + time_step/delta_x**2 * (2 * alpha_x[i,n] * bc_right \n",
|
|
" - (alpha_interblock(alpha_x[i,n-1], alpha_x[i,n]) + 2 * alpha_x[i,n]) * C_t[i,n]\n",
|
|
" + alpha_interblock(alpha_x[i,n-1], alpha_x[i,n]) * C_t[i,n-1]) \\\n",
|
|
" + time_step/delta_y**2 * (alpha_interblock(alpha_y[1,j], alpha_y[0,j]) * C_t[1,j] \n",
|
|
" - (alpha_interblock(alpha_y[1,j], alpha_y[0,j]) + 2 * alpha_y[0,j]) * C_t[0,j]\n",
|
|
" + 2 * alpha_y[0,j] * bc_top)\n",
|
|
"\n",
|
|
" # corner bottom left\n",
|
|
" i = grid_size['y']-1\n",
|
|
" m = i\n",
|
|
" j = 0\n",
|
|
" C_t1[i,j] = C_t[i,j] \\\n",
|
|
" + time_step/delta_x**2 * (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) * C_t[i,j+1] \n",
|
|
" - (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) + 2 * alpha_x[i,j]) * C_t[i,j]\n",
|
|
" + 2 * alpha_x[i,j] * bc_left) \\\n",
|
|
" + time_step/delta_y**2 * (2 * alpha_y[m,j] * bc_bottom \n",
|
|
" - (alpha_interblock(alpha_y[m,j], alpha_y[m-1,j]) + 2 * alpha_y[m,j]) * C_t[m,j]\n",
|
|
" + alpha_interblock(alpha_y[m,j], alpha_y[m-1,j]) * C_t[m-1,j])\n",
|
|
"\n",
|
|
" # corner bottom right\n",
|
|
" i = grid_size['y']-1\n",
|
|
" j = grid_size['x']-1\n",
|
|
" m = i \n",
|
|
" n = j\n",
|
|
" C_t1[i,j] = C_t[i,j] \\\n",
|
|
" + time_step/delta_x**2 * (2 * alpha_x[i,n] * bc_right \n",
|
|
" - (alpha_interblock(alpha_x[i,n-1], alpha_x[i,n]) + 2 * alpha_x[i,n]) * C_t[i,n]\n",
|
|
" + alpha_interblock(alpha_x[i,n-1], alpha_x[i,n]) * C_t[i,n-1]) \\\n",
|
|
" + time_step/delta_y**2 * (2 * alpha_y[m,j] * bc_bottom \n",
|
|
" - (alpha_interblock(alpha_y[m,j], alpha_y[m-1,j]) + 2 * alpha_y[m,j]) * C_t[m,j]\n",
|
|
" + alpha_interblock(alpha_y[m,j], alpha_y[m-1,j]) * C_t[m-1,j])\n",
|
|
"\n",
|
|
" return C_t1"
|
|
]
|
|
},
|
|
{
|
|
"attachments": {},
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"### Implementation for closed boundary conditions"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 32,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# simulate one time step\n",
|
|
"def simulate_closed_boundary(C_t): \n",
|
|
" C_t1 = np.zeros((grid_size['x'],grid_size['y']))\n",
|
|
"\n",
|
|
" # inner cells\n",
|
|
" for i in range(1, grid_size['y']-1): #rows\n",
|
|
" for j in range(1, grid_size['x']-1): #columns\n",
|
|
" C_t1[i,j] = C_t[i,j] \\\n",
|
|
" + time_step/delta_y**2 * (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) * C_t[i+1,j]\n",
|
|
" - (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j])) * C_t[i,j] \n",
|
|
" + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j]) * C_t[i-1,j]) \\\n",
|
|
" + time_step/delta_x**2 * (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) * C_t[i,j+1]\n",
|
|
" - (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j])) * C_t[i,j] \n",
|
|
" + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j]) * C_t[i,j-1])\n",
|
|
" \n",
|
|
" # boundary conditions\n",
|
|
" # left without corners / looping over rows\n",
|
|
" for i in range(1, grid_size['y']-1):\n",
|
|
" j = 0\n",
|
|
" C_t1[i,0] = C_t[i,0] \\\n",
|
|
" + time_step/delta_x**2 * (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) * (C_t[i,j+1] - C_t[i,j]))\\\n",
|
|
" + time_step/delta_y**2 * (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) * C_t[i+1,j]\n",
|
|
" - (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j])) * C_t[i,j] \n",
|
|
" + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j]) * C_t[i-1,j]) \n",
|
|
"\n",
|
|
" # right without corners / looping over rows\n",
|
|
" n = grid_size['x']-1 # maximum index in x-direction (columns)\n",
|
|
" for i in range(1, grid_size['y']-1):\n",
|
|
" j = n\n",
|
|
" C_t1[i,n] = C_t[i,n] \\\n",
|
|
" + time_step/delta_x**2 * ( \n",
|
|
" - (alpha_interblock(alpha_x[i,n-1], alpha_x[i,n]) * (C_t[i,n] - C_t[i,n-1])))\\\n",
|
|
" + time_step/delta_y**2 * (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) * C_t[i+1,j]\n",
|
|
" - (alpha_interblock(alpha_y[i+1,j], alpha_y[i,j]) + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j])) * C_t[i,j] \n",
|
|
" + alpha_interblock(alpha_y[i-1,j], alpha_y[i,j]) * C_t[i-1,j])\n",
|
|
"\n",
|
|
" # top without corners / looping over columns\n",
|
|
" for j in range(1, grid_size['x']-1):\n",
|
|
" i = 0\n",
|
|
" C_t1[0,j] = C_t[0,j] \\\n",
|
|
" + time_step/delta_y**2 * (alpha_interblock(alpha_y[1,j], alpha_y[0,j]) * (C_t[1,j] - C_t[0,j])) \\\n",
|
|
" + time_step/delta_x**2 * (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) * C_t[i,j+1]\n",
|
|
" - (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j])) * C_t[i,j] \n",
|
|
" + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j]) * C_t[i,j-1])\n",
|
|
"\n",
|
|
" # bottom without corners / looping over columns\n",
|
|
" m = grid_size['y']-1 # maximum index in y-direction (rows)\n",
|
|
" for j in range(1, grid_size['x']-1):\n",
|
|
" i = m\n",
|
|
" C_t1[m,j] = C_t[m,j] \\\n",
|
|
" + time_step/delta_y**2 * ( \n",
|
|
" - (alpha_interblock(alpha_y[m,j], alpha_y[m-1,j]) * (C_t[m,j] - C_t[m-1,j]))) \\\n",
|
|
" + time_step/delta_x**2 * (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) * C_t[i,j+1]\n",
|
|
" - (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j])) * C_t[i,j] \n",
|
|
" + alpha_interblock(alpha_x[i,j-1], alpha_x[i,j]) * C_t[i,j-1])\n",
|
|
"\n",
|
|
" # corner top left\n",
|
|
" i = 0\n",
|
|
" j = i\n",
|
|
" C_t1[i,j] = C_t[i,j] \\\n",
|
|
" + time_step/delta_x**2 * (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) * (C_t[i,j+1] - C_t[i,j]))\\\n",
|
|
" + time_step/delta_y**2 * (alpha_interblock(alpha_y[1,j], alpha_y[0,j]) * (C_t[1,j]-C_t[0,j]))\n",
|
|
" \n",
|
|
"\n",
|
|
" # corner top right\n",
|
|
" i = 0\n",
|
|
" j = grid_size['x']-1\n",
|
|
" n = j\n",
|
|
" C_t1[i,j] = C_t[i,j] \\\n",
|
|
" + time_step/delta_x**2 * (- (alpha_interblock(alpha_x[i,n-1], alpha_x[i,n]) * (C_t[i,n] - C_t[i,n-1])))\\\n",
|
|
" + time_step/delta_y**2 * (alpha_interblock(alpha_y[1,j], alpha_y[0,j]) * (C_t[1,j] - C_t[0,j])) \n",
|
|
" \n",
|
|
"\n",
|
|
" # corner bottom left\n",
|
|
" i = grid_size['y']-1\n",
|
|
" m = i\n",
|
|
" j = 0\n",
|
|
" C_t1[i,j] = C_t[i,j] \\\n",
|
|
" + time_step/delta_x**2 * (alpha_interblock(alpha_x[i,j+1], alpha_x[i,j]) * (C_t[i,j+1] - C_t[i,j])) \\\n",
|
|
" + time_step/delta_y**2 * (- (alpha_interblock(alpha_y[m,j], alpha_y[m-1,j] * (C_t[m,j] - C_t[m-1,j]))))\n",
|
|
"\n",
|
|
" # corner bottom right\n",
|
|
" i = grid_size['y']-1\n",
|
|
" j = grid_size['x']-1\n",
|
|
" m = i \n",
|
|
" n = j\n",
|
|
" C_t1[i,j] = C_t[i,j] \\\n",
|
|
" + time_step/delta_x**2 * (- alpha_interblock(alpha_x[i,n-1], alpha_x[i,n]) * (C_t[i,n] - C_t[i,n-1]))\\\n",
|
|
" + time_step/delta_y**2 * (- alpha_interblock(alpha_y[m,j], alpha_y[m-1,j] * (C_t[m,j] - C_t[m-1,j])))\n",
|
|
"\n",
|
|
" return C_t1"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 33,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAUEAAAD8CAYAAADpLRYuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAAaO0lEQVR4nO3dfaxddZ3v8ffHysOIcAGrHZ4UxltN1IxVGmDiw4WrQG2M4MRx2pkIKHeqI9zozUzmopOI0ZB474hmHAzcOjTAjRYZEW1mKlC5ZjomghQGgfIwFITYWtuUGsAnpO3n/rF+G7anZ5+zztr7dD+szytZOWv/9trrgYZvfr/fd631lW0iItrqRcM+gYiIYUoQjIhWSxCMiFZLEIyIVksQjIhWSxCMiFZLEIyIoZJ0gqTvSXpA0mZJHyvtR0vaIOmR8veo0i5JX5K0RdK9kt7cta/zy/aPSDq/1vFzn2BEDJOkY4BjbN8t6XDgLuBc4AJgt+3PSboEOMr2/5S0HPjvwHLgVODvbZ8q6WhgE7AUcNnPybZ/PtPx0xOMiKGyvd323WX9GeBB4DjgHODastm1VIGR0n6dK7cDR5ZAejawwfbuEvg2AMtmO/6LB3kxg3KwDvGhHDbs04iYWL/hl/zWz6qffZx9xmF+cvfeWtvede+zm4HfdDWttr166naSTgTeBNwBLLK9vXz1M2BRWT8O+EnXz7aWtl7tMxrJIHgoh3Gq3jHs04iYWHf4tr73sWv3Xu645fha2x50zKO/sb10pm0kvRS4Efi47aelF2K0bUual7m7DIcjoiGz1/tqLbORdBBVAPyq7W+W5h1lmNuZN9xZ2rcBJ3T9/PjS1qt9Rn0FQUnLJD1csjSXTPP9IZK+Xr6/o3R1I2ICGNiHay0zUdXluxp40PYXur5aB3QyvOcD3+5qP69kiU8DnirD5luAsyQdVTLJZ5W2GTUeDktaAHwZOJNq7H2npHW2H+ja7ELg57b/s6QVwP8C/rTpMSNitOxj9l5eDW8BPgDcJ+me0vZJ4HPADZIuBJ4A3l++W0+VGd4C/Ar4IIDt3ZI+C9xZtvuM7d2zHbyfOcFTgC22HwOQdD1V1qY7CJ4DfLqsfwO4QpKc+3Iixp4xz9UY6s66H/v7QK8kzX7JgRI/LuqxrzXAmrkcv5/hcJ1MzPPb2N4DPAW8bLqdSVolaZOkTc/xbB+nFREHgoG9uNYyykYmO1zS5asBjtDRo/1fLSIAZp3vGwf9BME6mZjONlslvRj4T8CTfRwzIkaEgb0TMLPVz3D4TmCxpJMkHQysoMradOvO7rwP+H+ZD4yYHPtqLqOscU/Q9h5JF1OloBcAa2xvlvQZYJPtdVRp7/8raQuwmypQRsQE8BjM99XR15yg7fVU6erutk91rf8G+JN+jhERo8mG58Y/Bo5OYiQixo3Y2/POlvGRIBgRjRjYl55gRLRZeoIR0VrVzdIJghHRUgae8/i/iCpBMCIaMWLvBLyNL0EwIhrb5wyHI6KlMicYES0n9mZOMCLaqnqzdIJgRLSULX7rBcM+jb4lCEZEY/smYE5w/PuyETEUVWLkRbWW2UhaI2mnpPu72r4u6Z6yPN6pPyLpREm/7vruqq7fnCzpvlLc7UvqrtvZQ3qCEdHQQBMj1wBXANd1Gmw/X5RN0uVU5Tk6HrW9ZJr9XAn8BVXx9vXAMuA7Mx04PcGIaKSTGKmzzLoveyPVO0f3U3pz7wfWzrSPUpv4CNu3l5c3XwecO9uxEwQjorG9Vq2lT28Ddth+pKvtJEn/LulfJb2ttB1HVfCtY7rib/vJcDgiGjHiOdcOIQslber6vLoUV6tjJb/bC9wOvNL2k5JOBr4l6fV1T2SqBMGIaKSTGKlpl+2lcz1GKdD2x8DJzx/Xfhaqury275L0KPAaqsJux3f9fLrib/tpPByWdIKk70l6QNJmSR+bZpvTJT3VlcX51HT7iojxY+oNhfscDr8TeMj288NcSS+XtKCs/wGwGHjM9nbgaUmnlXnE84Bvz3aAfnqCe4C/sn23pMOBuyRtsP3AlO3+zfa7+zhORIyoQT0xImktcDrVsHkrcKntq6mKs01NiLwd+Iyk56iK2X3Ediep8lGqTPPvUWWFZ8wMQ3/V5rZTjc2x/YykB6kmIacGwYiYQDYDu0XG9soe7RdM03YjcGOP7TcBb5jLsQcyJyjpROBNVPfmTPVHkn4E/BT4a9ube+xjFbAK4FBeMojTioh5VCVG8tgckl5KFZU/bvvpKV/fDbzK9i8kLQe+RTV+30/JFK0GOEJHT0D5lojJNwkvVe3rCiQdRBUAv2r7m1O/t/207V+U9fXAQZIW9nPMiBgNRuxzvWWUNe4JluzL1cCDtr/QY5vfp7rJ0ZJOoQq6TzY9ZkSMlknoCfYzHH4L8AHgvs6DzcAngVcC2L4KeB/wl5L2AL8GVpTHWSJizFV1h1scBG1/H2Z+j47tK6geio6IiaO8Xj8i2qsquZnscES0lK12D4cjIlJoKSJaq3qfYOYEI6K1UnIzIlqsukUmPcGIaKk8OxwRrZfi6xHRWtWrtDIcjogWy5xgRLRW9RaZDIcjoqWqx+YSBCOitSajJzj+VxARQ7MP1VpmI2mNpJ2S7u9q+7SkbV3VKpd3ffcJSVskPSzp7K72ZaVti6RL6lxDeoIR0ciAs8PXUL1277op7V+0/fnuBkmvo6pC93rgWOC7kl5Tvv4ycCawFbhT0rppKmD+jgTBiGhsUMNh2xtLwbY6zgGuL0XYfyxpC3BK+W6L7ccAJF1ftp0xCGY4HBGNzLHGyEJJm7qWVTUPc7Gke8tw+ajSdhzwk65ttpa2Xu0zSk8wIhoxsKd+T3CX7aVzPMSVwGfLoT4LXA58aI77mNUgSm4+DjwD7AX2TL3QUpDp74HlwK+AC2zf3e9xI2L45jM7bHtHZ13SV4B/Lh+3ASd0bXp8aWOG9p4G1RM8w/auHt+9i6rW8GLgVKrofuqAjhsRwzLP5TQlHWN7e/n4XqCTOV4HfE3SF6gSI4uBH1LVPFos6SSq4LcC+LPZjnMghsPnANeVKnO3SzpyysVFxBga5EtVJa0FTqeaO9wKXAqcLmlJOdTjwIcBbG+WdANVwmMPcJHtvWU/FwO3AAuANbY3z3bsQQRBA7dKMvB/bK+e8n2vycrfCYJlonQVwKG8ZACnFRHzbVA9Qdsrp2m+eobtLwMum6Z9PbB+LsceRBB8q+1tkl4BbJD0kO2Nc91JCZ6rAY7Q0alNHDHi8lLVwva28nenpJuo7tfpDoIzTWJGxJgyYs++8b/Lrq8rkHSYpMM768BZvDB52bEOOE+V04CnMh8YMRkG9djcMPXbE1wE3FTdBcOLga/ZvlnSRwBsX0U1Pl8ObKG6ReaDfR4zIkaBMxymPJ7yxmnar+paN3BRP8eJiNGTOcGIaL0EwYhoLSP2TkBiJEEwIhob9aRHHQmCEdGIkxiJiLZzgmBEtNf8vkDhQEkQjIjG0hOMiNayYe++BMGIaLFkhyOitUyGwxHRakmMRETLeQLe/JkgGBGNTcJwePwf/IuIoaiywy+qtcym1BXeKen+rra/k/RQqTt8k6QjS/uJkn4t6Z6yXNX1m5Ml3Sdpi6QvlWqXM0oQjIjG7HpLDdcAy6a0bQDeYPsPgf8APtH13aO2l5TlI13tVwJ/wQsVLqfucz8JghHRmK1ay+z78UZg95S2W23vKR9vpyrN0ZOkY4AjbN9e3mN6HXDubMdOEIyIRky9AFiC4EJJm7qWVXM83IeA73R9PknSv0v6V0lvK23HUVWz7OhUtpxREiMR0dgcksO7bC9tcgxJf0tVX/irpWk78ErbT0o6GfiWpNc32Tf00ROU9Nquicl7JD0t6eNTtjld0lNd23yq6fEiYsQYvE+1lqYkXQC8G/jzMsTF9rO2nyzrdwGPAq+hqmLZPWSuVdmycU/Q9sPAknKiC8rBbppm03+z/e6mx4mI0TWft8hIWgb8DfBfbP+qq/3lwG7beyX9AVUC5DHbu0tn7DTgDuA84B9mO86ghsPvoMrWPDGg/UXEGBjUzdKS1gKnU80dbgUupcoGHwJsKHe63F4ywW8HPiPpOWAf8BHbnaTKR6kyzb9HNYfYPY84rUEFwRXA2h7f/ZGkHwE/Bf7a9ubpNioTpasADuUlAzqtiJgvg3x22PbKaZqv7rHtjcCNPb7bBLxhLsfuOzss6WDgPcA/TfP13cCrbL+Rqlv6rV77sb3a9lLbSw/ikH5PKyLmmwGr3jLCBnGLzLuAu23vmPqF7adt/6KsrwcOkrRwAMeMiBEwwJulh2YQw+GV9BgKS/p9YIdtSzqFKug+OYBjRsTQ9Zf5HRV9BUFJhwFnAh/uavsIgO2rgPcBfylpD/BrYEUnzR0RE2AC/m/uKwja/iXwsiltV3WtXwFc0c8xImJEeTLeIpMnRiKiubb3BCOi7dITjIg22zfsE+hfgmBENNO5T3DMJQhGRGOTcK9HgmBENJcgGBGtluFwRLSZ0hOMiNayoO2PzUVEy6UnGBGtliAYEa2WIBgRrTUhN0un7nBENCbXW2bdj7RG0k5J93e1HS1pg6RHyt+jSrskfUnSFkn3Snpz12/OL9s/Iun8OteQIBgRzbnmMrtrgGVT2i4BbrO9GLitfIbqbfaLy7IKuBKqoElVoOlU4BTg0k7gnEmCYEQ0NqieoO2NwO4pzecA15b1a4Fzu9qvc+V24EhJxwBnAxts77b9c2AD+wfW/WROMCKaqz8nuFDSpq7Pq22vnuU3i2xvL+s/AxaV9eOAn3Rtt7W09WqfUYJgRDRTf6gLsMv20saHquoUzUsuutZweC6TltP8ds4TlRExJgY3JzidHWWYS/m7s7RvA07o2u740tarfUZ15wSvof6k5fOaTlRGxHjQvnpLQ+uATsfpfODbXe3nlSzxacBTZdh8C3CWpKNKnDmrtM2oVhCc46Rlt0YTlRExJgbUE5S0FvgB8FpJWyVdCHwOOFPSI8A7y2eA9cBjwBbgK8BHAWzvBj4L3FmWz5S2GfUzJ9hr0rJbo4nKiBh9dTO/ddhe2eOrd0yzrYGLeuxnDbBmLsceSGJkEJOWklZR3fPDobxkEKcVEfOt5U+M9Jq07FZ7otL2attLbS89iEP6OK2IOGDmNzFyQPQTBHtNWnZrNFEZEeNhUDdLD1PdW2RqT1pKWirpH6H5RGVEjAHPe3b4gKg1JzjHSctNwH/r+jznicqIGBMj3surI0+MRERzCYIR0WajPt9XR94iExGtlp5gRDQ3AT3BBMGIaMajn/mtI0EwIppLTzAi2kpMRmIkQTAimksQjIjWGoNH4upIEIyI5pIYiYg2S08wItotQTAiWmsM3hVYRx6bi4jGBvE+QUmvlXRP1/K0pI9L+rSkbV3ty7t+8wlJWyQ9LOnsfq4hPcGIaG4APUHbDwNLACQtoHr7/E3AB4Ev2v589/aSXgesAF4PHAt8V9JrbO9tcvz0BCOisXl4qeo7gEdtPzHDNucA19t+1vaPqarOndL0GhIEI6KZuvVFqt7iQkmbupZVPfa6Aljb9fliSfdKWtNVs3ygVSwTBCOiEc1hAXZ1CqmVZfV++5MOBt4D/FNpuhJ4NdVQeTtw+XxcR4JgRDQ32Gpz7wLutr0DwPYO23tt76Mqst4Z8tauYlnHrEGwdEN3Srq/q+3vJD1Uuqk3STqyx28fl3RfyexsanqSETGaBlxtbiVdQ+FOSd/ivUAnBq0DVkg6RNJJwGLgh02voU5P8Bpg2ZS2DcAbbP8h8B/AJ2b4/Rm2l9he2uwUI2JkDagnKOkw4Ezgm13N/7t0ou4FzgD+B4DtzcANwAPAzcBFTTPDUOMWGdsbJZ04pe3Wro+3A+9regIRMaYG+FJV278EXjal7QMzbH8ZcNkgjj2IOcEPAd/p8Z2BWyXdNUM2CABJqzqZo+d4dgCnFRHzbrBzgkPR183Skv4W2AN8tccmb7W9TdIrgA2SHrK9cboNS7ZoNcAROnrE/7NFBEzGCxQa9wQlXQC8G/hz29P+p7C9rfzdSXUHeOMbGiNiBE1AT7BREJS0DPgb4D22f9Vjm8MkHd5ZB87ihexOREyAAWeHh6LOLTJrgR8Ar5W0VdKFwBXA4VRD3HskXVW2PVbS+vLTRcD3Jf2IKn39L7ZvnperiIgDz1QvVa2zjLA62eGV0zRf3WPbnwLLy/pjwBv7OruIGFkptBQRkSAYEW2m6XOiYyVBMCKaGYPMbx0JghHRWOYEI6LVBvXY3DAlCEZEc+kJRkRrjcGN0HUkCEZEcwmCEdFWuVk6IlpP+8Y/CiYIRkQzuU8wItpuEm6RSbW5iGhucDVG9ivKJuloSRskPVL+HlXaJelLkraUYm9v7ucSEgQjorEBv09walG2S4DbbC8GbiufoSrNubgsq6jqEzeWIBgRzRiw6y3NnANcW9avBc7tar/OlduBI6eU55yTBMGIaEz76i3Awk4htbJMLbw2XVG2Rba3l/WfUb2oGeA44Cddv91a2hpJYiQiGpnjfYK7Zqk9vl9Rtu4vbVuan7sS0xOMiGbqDoVrDId7FGXb0Rnmlr87y+bbgBO6fn58aWukTo2RNZJ2Srq/q+3TkraVTM49kpb3+O0ySQ+XLM4l020TEeNrEImRGYqyrQPOL5udD3y7rK8DzitZ4tOAp7qGzXNWZzh8DVVhpeumtH/R9ud7/UjSAuDLwJlUY/Y7Ja2z/UDDc42IUTOYAeoi4CZJUMWkr9m+WdKdwA2luNsTwPvL9uupahltAX4FfLCfg9cptLRR0okN9n0KsKUUXELS9VRZnQTBiAkxiFm6XkXZbD8JvGOadgMX9X/kSj9zgheXGxXXdG5inGJOGRxJqzqZo+d4to/TiogDwsBe11tGWNMgeCXwamAJsB24vN8Tsb3a9lLbSw/ikH53FxEHwCQUX290i4ztHZ11SV8B/nmazQaawYmIETQB1eYa9QSn3J39XqpMzlR3AoslnSTpYGAFVVYnIiZEK3qCktYCp1Pd8b0VuBQ4XdISqlmBx4EPl22PBf7R9nLbeyRdDNwCLADW2N48HxcREUPQlldp2V45TfPVPbb9KVXquvN5PVU6OyImjACNeNKjjjw2FxGNaQLmBBMEI6KZtgyHIyKm19drskZGgmBENDbqmd86EgQjorn0BCOitZzscES03fjHwATBiGgut8hERLslCEZEaxmYgOLrCYIR0YhwhsMR0XL7xr8rmGpzEdFMZzhcZ5mBpBMkfU/SA5I2S/pYae9Z0E3SJ0oBt4clnd3PZaQnGBGNDWg4vAf4K9t3l6pzd0naUL7br6CbpNdRvZ/09cCxwHclvcb23iYHT08wIpobQN1h29tt313WnwEeZIZ6RFQF2663/aztH1NVnTul6SUkCEZEQ4Mrvt5RKlu+CbijNE1X0G1ORdxmkyAYEc3Mrdrcwk41ybKsmro7SS8FbgQ+bvtp5qGg23QyJxgRjc1hTnCX7aU99yMdRBUAv2r7mzBjQbeBFnGbtSdYuqE7Jd3f1fb1rozN45Lu6fHbxyXdV7bb1PQkI2JEDWA4LElUJTsetP2FrvZeBd3WASskHSLpJGAx8MOml1CnJ3gNcAVwXafB9p92nejlwFMz/P4M27uanmBEjCgD+waSHX4L8AHgvq4O1SeBldMVdLO9WdINwANUmeWLmmaGoV6hpY1lsnI/JYK/H/ivTU8gIsbVYN4sbfv7VHWbpupZpM32ZcBlfR+c/hMjbwN22H6kx/cGbpV013QTod0krepMmj7Hs32eVkQcEAPODg9Dv4mRlcDaGb5/q+1tkl4BbJD0kO2N021oezWwGuAIHT3a/9UiomSHW/zYnKQXA38MfL3XNra3lb87gZvo44bGiBg1Bu+rt4ywfobD7wQesr11ui8lHVYegUHSYcBZvJDdiYhJMAHD4Tq3yKwFfgC8VtJWSReWr1YwZSgs6VhJncnMRcD3Jf2IKn39L7ZvHtypR8RQdbLDdZYRVic7vLJH+wXTtP0UWF7WHwPe2Of5RcQoG/FeXh15YiQimksQjIjWsmFv43uUR0aCYEQ0l55gRLRagmBEtNfoZ37rSBCMiGYMHvEboetIEIyI5ibgsbkEwYhoxp6IkpsJghHRXBIjEdFmTk8wItpr9F+OUEeCYEQ0M7jX6w9VgmBENGLAE/DYXOoOR0QzHtxLVSUtk/SwpC2SLjkAZ/+89AQjojEPYDgsaQHwZeBMYCtwp6R1th/oe+c1pCcYEc0Npid4CrDF9mO2fwtcD5wz7+dejGRP8Bl+vuu7/sYTU5oXApNYv3hSrwsm99om4bpe1e8OnuHnt3zX31hYc/NDJW3q+ry6FFcDOA74Sdd3W4FT+z2/ukYyCNp++dQ2SZtsLx3G+cynSb0umNxrm9Trmivby4Z9DoOQ4XBEDNs24ISuz8eXtgMiQTAihu1OYLGkkyQdTFXEbd2BOvhIDod7WD37JmNpUq8LJvfaJvW6hsL2HkkXA7cAC4A1tjcfqOPLE/DYS0REUxkOR0SrJQhGRKuNRRAc5iM180nS45Luk3TPlHuoxo6kNZJ2Srq/q+1oSRskPVL+HjXMc2yix3V9WtK28u92j6TlwzzH6M/IB8GuR2reBbwOWCnpdcM9q4E6w/aSCbjv7Bpg6n1jlwC32V4M3FY+j5tr2P+6AL5Y/t2W2F5/gM8pBmjkgyBDfqQm6rG9Edg9pfkc4Nqyfi1w7oE8p0HocV0xQcYhCE73SM1xQzqXQTNwq6S7JK0a9snMg0W2t5f1nwGLhnkyA3axpHvLcHnshvnxgnEIgpPsrbbfTDXUv0jS24d9QvPF1b1Yk3I/1pXAq4ElwHbg8qGeTfRlHILgUB+pmU+2t5W/O4GbqIb+k2SHpGMAyt+dQz6fgbC9w/ZeV0V3v8Lk/bu1yjgEwaE+UjNfJB0m6fDOOnAWcP/Mvxo764Dzy/r5wLeHeC4D0wnsxXuZvH+3Vhn5x+aG/UjNPFoE3CQJqn+Hr9m+ebin1JyktcDpwEJJW4FLgc8BN0i6EHgCeP/wzrCZHtd1uqQlVMP7x4EPD+v8on95bC4iWm0chsMREfMmQTAiWi1BMCJaLUEwIlotQTAiWi1BMCJaLUEwIlrt/wPp94uV7Lx+gQAAAABJRU5ErkJggg==",
|
|
"text/plain": [
|
|
"<Figure size 432x288 with 2 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"records = []\n",
|
|
"\n",
|
|
"for i in range(1000):\n",
|
|
" C_t = simulate_closed_boundary(C_t)\n",
|
|
" records.append(C_t)\n",
|
|
"\n",
|
|
"# with np.printoptions(precision=2,floatmode='fixed',linewidth=200):\n",
|
|
"# print(C_t)\n",
|
|
"\n",
|
|
"plt.imshow(C_t, vmin=0, vmax=2000)\n",
|
|
"plt.colorbar()\n",
|
|
"plt.show()\n"
|
|
]
|
|
},
|
|
{
|
|
"attachments": {},
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"### Widget"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 37,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"application/vnd.jupyter.widget-view+json": {
|
|
"model_id": "324063b36c0147808a548646753cb24a",
|
|
"version_major": 2,
|
|
"version_minor": 0
|
|
},
|
|
"text/plain": [
|
|
"interactive(children=(IntSlider(value=0, description='w', max=999), Output()), _dom_classes=('widget-interact'…"
|
|
]
|
|
},
|
|
"metadata": {},
|
|
"output_type": "display_data"
|
|
},
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"<function __main__.update(w=1)>"
|
|
]
|
|
},
|
|
"execution_count": 37,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"def update(w = 1):\n",
|
|
" fig = plt.figure(figsize = (10,7))\n",
|
|
" y = records[w]\n",
|
|
" plt.imshow(y, vmin=0, vmax=50)\n",
|
|
" \n",
|
|
"interact(update, w = IntSlider(min=0, max = 999, step = 1, value = 0))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 35,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAQUAAAD4CAYAAADl7fPiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAAREUlEQVR4nO3dfYxc1X3G8e/jtaGVMbGJwQHjEppaSFZU3MhyGpVWpiTURihOqjS1VbVOS7UkClIjNapoK0GU/kNVUaTECOIkFqRKgL45sRQLsGglgpQXFmRenEBxLUfs4thNTHFwiGDtX//Yu9Gc9cz63rn37NyZPh9pNXfuPXPuuTvLw7z8fI4iAjOzWYsGPQAzaxeHgpklHApmlnAomFnCoWBmicWDHkA3YxcujSUXLx/0MErL8wWOKgygQre5vmyq1G+ma6siR7+Zxlrht1XaW6+e4PSpU127bmUoLLl4OWvu+Hjj/Vb6jzfKPxVnSraNM+X7rNT2dIU/m1xtK/y+VGkM5ZtW6VdnSjas8DyU7hMqBUiOUJj83F09j/ntg5klaoWCpM2SXpR0SNKtXY6fL+mh4vh3Jb2zzvnMLL++Q0HSGHA3sAVYB2yXtG5Os5uAVyPi14C7gL/v93xmtjDqvFLYCByKiMMR8SbwILB1TputwP3F9r8C10nK8RbJzBpSJxRWAy933J8s9nVtExHTwGvA27t1Jmlc0oSkidMnT9UYlpnV0ZoPGiNiV0RsiIgNYxcuHfRwzP7fqhMKU8CajvuXF/u6tpG0GHgb8JMa5zSzzOqEwpPAWklXSjoP2AbsndNmL7Cj2P4I8B/hf6tt1mp9Fy9FxLSkW4BHgDFgd0QclPRZYCIi9gJfBv5J0iHgBDPBYWYtVquiMSL2Afvm7LutY/vnwB/UOce5x1ClcYUqwQptS/fbhnLkXEb52oZJA7/b1nzQaGbt4FAws4RDwcwSDgUzSzgUzCzhUDCzhEPBzBIOBTNLOBTMLOFQMLNEKyduhQrly5lKlwdeiZujzDpr2/JNK8k1hpJtlal8u9JMQwv8x+hXCmaWcCiYWcKhYGYJh4KZJRwKZpZwKJhZwqFgZok6K0StkfSfkr4v6aCkv+jSZpOk1yQdKH5u69aXmbVHneKlaeAvI+JpScuApyTtj4jvz2n3rYi4scZ5zGwB9f1KISKORsTTxfZPgR9w9gpRZjZkGilzLlaT/g3gu10Ov0/SM8ArwKcj4mCPPsaBcYDFK99GnGn+445K1aJVSqLPlG1XpSS7dFMoef6Zjiu0rUBVypErjLdKmXG1kuTmZ+DOVrqc63fQQ+3/8iRdAPwb8KmIODnn8NPAFRFxNfB54Ou9+ulcNm7RMi8bZzYotUJB0hJmAuGrEfHvc49HxMmIeL3Y3gcskbSyzjnNLK863z6ImRWgfhAR/9ijzTtml56XtLE4n9eSNGuxOp8p/Bbwx8Bzkg4U+/4G+BWAiLiXmfUjPyFpGngD2Oa1JM3arc5akk9wjs9WImInsLPfc5jZwnNFo5klHApmlnAomFnCoWBmCYeCmSXaO5tz6YaZSocrtC1dvlypFDjTrMtV+s1VPp2pJLpSOXDJfrPN5jzo0vR5+vQrBTNLOBTMLOFQMLOEQ8HMEg4FM0s4FMws4VAws4RDwcwSDgUzS7S3ovF0pWkwy/VZpZIuR0VjpWrCNlT9ZZqMNVvbDM9vG6oUh23iVjMbLQ4FM0s0McX7EUnPFcvCTXQ5Lkmfk3RI0rOS3lP3nGaWT1OfKVwbET/ucWwLsLb4eS9wT3FrZi20EG8ftgJfiRnfAZZLunQBzmtmfWgiFAJ4VNJTxdJvc60GXu64P0mXNScljUuakDRx+uSpBoZlZv1o4u3DNRExJekSYL+kFyLi8aqdRMQuYBfA+b+62mtDmA1I7VcKETFV3B4H9gAb5zSZAtZ03L+82GdmLVR3LcmlkpbNbgPXA8/PabYX+JPiW4jfBF6LiKN1zmtm+dR9+7AK2FMsF7kY+FpEPCzp4/CLpeP2ATcAh4CfAX9a85xmllGtUIiIw8DVXfbf27EdwCerdaw8k6FWGUKOktlcpcsV+s1WulyhLD1b6XKG8ulcpcuVrqsClzmbWeMcCmaWcCiYWcKhYGYJh4KZJRwKZpZwKJhZwqFgZgmHgpklHApmlmjxbM4Z8irTDLqly4xzlE5TsWT2dJ5+K5UuV5mpe8CzRA98NumqbRvo068UzCzhUDCzhEPBzBIOBTNLOBTMLOFQMLOEQ8HMEn2HgqSriqXiZn9OSvrUnDabJL3W0ea22iM2s6z6Ll6KiBeB9QCSxpiZtn1Pl6bfiogb+z2PmS2spt4+XAf8d0T8sKH+zGxAmipz3gY80OPY+yQ9A7wCfDoiDnZrVCw5Nw4wdtFyqFIKW1au0tKyM0q3YhbjCqXLVUqih6h0uUrbbM9Dtr/FCm17aGIp+vOADwL/0uXw08AVEXE18Hng6736iYhdEbEhIjaMLVtad1hm1qcm3j5sAZ6OiGNzD0TEyYh4vdjeByyRtLKBc5pZJk2EwnZ6vHWQ9A4Vy0dJ2lic7ycNnNPMMqn1mUKxfuQHgJs79nUuGfcR4BOSpoE3gG3FilFm1lJ1l407Bbx9zr7OJeN2AjvrnMPMFpYrGs0s4VAws4RDwcwSDgUzSzgUzCzRztmcgzxlzlXHUFLpUtgKZbBVZn7OVgo8bKXLFcabYzbnVpQ5N3B+v1Iws4RDwcwSDgUzSzgUzCzhUDCzhEPBzBIOBTNLOBTMLOFQMLOEQ8HMEu0sc0bVSmzLGvAMuvlmaK7SNk/5dKXxVpklOtfs0znKnKuUxg+4zHk+fqVgZolSoSBpt6Tjkp7v2HeRpP2SXipuV/R47I6izUuSdjQ1cDPLo+wrhfuAzXP23Qo8FhFrgceK+wlJFwG3A+8FNgK39woPM2uHUqEQEY8DJ+bs3grcX2zfD3yoy0N/D9gfESci4lVgP2eHi5m1SJ3PFFZFxNFi+0fAqi5tVgMvd9yfLPaZWUs18kFjsZZDrc9IJY1LmpA0cfr115sYlpn1oU4oHJN0KUBxe7xLmylgTcf9y4t9Z0nWkrzgghrDMrM66oTCXmD224QdwDe6tHkEuF7SiuIDxuuLfWbWUmW/knwA+DZwlaRJSTcBdwAfkPQS8P7iPpI2SPoSQEScAP4OeLL4+Wyxz8xaqlRFY0Rs73Houi5tJ4A/77i/G9jd1+jMbMG1s8w5gLIlq5Fn1uccZajVyoYzzeZcpRS3QtlwlTLnRblmic5R5pypHDlX+XQTfbrM2cwSDgUzSzgUzCzhUDCzhEPBzBIOBTNLOBTMLOFQMLOEQ8HMEg4FM0u0tsxZ082XL2ebQbdkqXWl8+cqXc40Q/MwlS5XaZttNudMz28T/ErBzBIOBTNLOBTMLOFQMLOEQ8HMEg4FM0s4FMwscc5Q6LGO5D9IekHSs5L2SFre47FHJD0n6YCkiQbHbWaZlHmlcB9nL/W2H3h3RPw68F/AX8/z+GsjYn1EbOhviGa2kM4ZCt3WkYyIRyNiurj7HWYWeTGzEdBEmfOfAQ/1OBbAo5IC+EJE7OrViaRxYBxg8fIVLMpQ5lxJjtmcq5w/Q8lu9bajWbpcpW2lsbagjL303+I87WqFgqS/BaaBr/Zock1ETEm6BNgv6YXilcdZisDYBXD+mjULXO1tZrP6/vZB0seAG4E/KhaYPUtETBW3x4E9wMZ+z2dmC6OvUJC0Gfgr4IMR8bMebZZKWja7zcw6ks93a2tm7VHmK8lu60juBJYx85bggKR7i7aXSdpXPHQV8ISkZ4DvAd+MiIezXIWZNeacnyn0WEfyyz3avgLcUGwfBq6uNTozW3CuaDSzhEPBzBIOBTNLOBTMLOFQMLNEK2dzVlQrLx24BkpLz5JtNufypcuVnoM2lGVnKJ8edJk1kKXkfr52fqVgZgmHgpklHApmlnAomFnCoWBmCYeCmSUcCmaWcCiYWcKhYGaJVlY00oaKxgyzRFabgLNC5WGuSroWVPMNum2+SWbL/zFUes4a4FcKZpZwKJhZot9l4z4jaaqYn/GApBt6PHazpBclHZJ0a5MDN7M8+l02DuCuYjm49RGxb+5BSWPA3cAWYB2wXdK6OoM1s/z6WjaupI3AoYg4HBFvAg8CW/vox8wWUJ3PFG4pVp3eLWlFl+OrgZc77k8W+7qSNC5pQtLE6VOnagzLzOroNxTuAd4FrAeOAnfWHUhE7IqIDRGxYWzp0rrdmVmf+gqFiDgWEacj4gzwRbovBzcFrOm4f3mxz8xarN9l4y7tuPthui8H9ySwVtKVks4DtgF7+zmfmS2cc1Y0FsvGbQJWSpoEbgc2SVrPTN3fEeDmou1lwJci4oaImJZ0C/AIMAbsjoiDOS7CzJqTbdm44v4+4KyvK88pYNFbJct8y1cDD5dKk7FW6DdX2XCm8Q56ktdqfeYpXc4yyasnbjWzshwKZpZwKJhZwqFgZgmHgpklHApmlnAomFnCoWBmCYeCmSUcCmaWaOdszoCmyzas0mn5phUmU86i2szPmfqtIEspbs5+S7bNNetylVmiF2WYUXq+6/crBTNLOBTMLOFQMLOEQ8HMEg4FM0s4FMws4VAws0SZORp3AzcCxyPi3cW+h4CriibLgf+NiPVdHnsE+ClwGpiOiA2NjNrMsilTvHQfsBP4yuyOiPjD2W1JdwKvzfP4ayPix/0O0MwWVpmJWx+X9M5uxyQJ+Cjwuw2Py8wGpG6Z828DxyLipR7HA3hUUgBfiIhdvTqSNA6MAyy5cAVjb5UbwJmxCqOtUro8RJ+2DLoUuDX9ZpjVOtesy1VKlxdNl/8lLCr5zwPmew7qhsJ24IF5jl8TEVOSLgH2S3qhWLD2LEVg7AL45UvXZKrQN7Nz6fv/h5IWA78PPNSrTURMFbfHgT10X17OzFqkzovk9wMvRMRkt4OSlkpaNrsNXE/35eXMrEXOGQrFsnHfBq6SNCnppuLQNua8dZB0maTZFaFWAU9Iegb4HvDNiHi4uaGbWQ79LhtHRHysy75fLBsXEYeBq2uOz8wW2BB9xm5mC8GhYGYJh4KZJRwKZpZwKJhZopWzOesMjL1RsvEvle/3TIWrrVKyWnrm5ypl1m2o6cw1hjZcW1ktKMkuW7oMsPjnJTueZ5Zqv1Iws4RDwcwSDgUzSzgUzCzhUDCzhEPBzBIOBTNLOBTMLOFQMLOEQ8HMEopoX82ppP8Bfjhn90pgFNePGNXrgtG9tlG4risi4uJuB1oZCt1ImhjFFaZG9bpgdK9tVK9rlt8+mFnCoWBmiWEKhZ6rSw25Ub0uGN1rG9XrAoboMwUzWxjD9ErBzBaAQ8HMEkMRCpI2S3pR0iFJtw56PE2RdETSc5IOSJoY9HjqkLRb0nFJz3fsu0jSfkkvFbcrBjnGfvS4rs9ImiqetwOSbhjkGJvW+lCQNAbcDWwB1gHbJa0b7KgadW1ErB+B773vAzbP2Xcr8FhErAUeK+4Pm/s4+7oA7iqet/URsa/L8aHV+lBgZqXqQxFxOCLeBB4Etg54TDZHRDwOnJizeytwf7F9P/ChhRxTE3pc10gbhlBYDbzccX+y2DcKAnhU0lOSxgc9mAxWRcTRYvtHzCw6PCpukfRs8fZi6N4WzWcYQmGUXRMR72HmrdEnJf3OoAeUS8x89z0q33/fA7wLWA8cBe4c6GgaNgyhMAWs6bh/ebFv6EXEVHF7HNjDzFulUXJM0qUAxe3xAY+nERFxLCJOR8QZ4IuM2PM2DKHwJLBW0pWSzgO2AXsHPKbaJC2VtGx2G7geeH7+Rw2dvcCOYnsH8I0BjqUxs0FX+DAj9ry1coWoThExLekW4BFgDNgdEQcHPKwmrAL2SIKZ5+FrEfHwYIfUP0kPAJuAlZImgduBO4B/lnQTM/8U/qODG2F/elzXJknrmXk7dAS4eVDjy8FlzmaWGIa3D2a2gBwKZpZwKJhZwqFgZgmHgpklHApmlnAomFni/wCOWVT0CjszYQAAAABJRU5ErkJggg==",
|
|
"text/plain": [
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"import matplotlib.animation as animation\n",
|
|
"\n",
|
|
"\n",
|
|
"frames = []\n",
|
|
"fig = plt.figure()\n",
|
|
"for i in records:\n",
|
|
" frames.append([plt.imshow(i, vmin=0, vmax=10, animated=True)])\n",
|
|
"\n",
|
|
"ani = animation.ArtistAnimation(fig, frames, interval=50, blit=True, repeat_delay=1000)\n",
|
|
"\n",
|
|
"ani.save(\"/Users/hannessigner/Desktop/test.gif\")"
|
|
]
|
|
}
|
|
],
|
|
"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.9.7"
|
|
},
|
|
"orig_nbformat": 4
|
|
},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 2
|
|
}
|