mirror of
https://git.gfz-potsdam.de/naaice/tug.git
synced 2025-12-13 09:28:23 +01:00
611 lines
52 KiB
Plaintext
611 lines
52 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 40,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# imports\n",
|
|
"import numpy as np\n",
|
|
"import matplotlib.pyplot as plt\n",
|
|
"\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 2,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# set up environment\n",
|
|
"dim = 2\n",
|
|
"grid_size = {'x': 20, 'y': 20}\n",
|
|
"domain_size = {'x': 20, 'y': 20}\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.random.random((grid_size['x'],grid_size['y']))\n",
|
|
"alpha_y = np.ones((grid_size['x'],grid_size['y']))\n",
|
|
"\n",
|
|
"time_step = 0.1"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 3,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAQUAAAD4CAYAAADl7fPiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAANaklEQVR4nO3df6xf9V3H8efLUliCRUCk45cbmZUEllmXpriIBmTjR0PsZpZZYhSVBFxG4hITg5qMZf4zY5BoIJBuNjCzAUata7JCaaoJI9kYhZRfG0glXegda906YcgcFt7+cU/N/dx+b3v9nu/33u+9PB/Jzfecz+dzzvmcfJNXzznfb7/vVBWSdMRPLPYEJE0WQ0FSw1CQ1DAUJDUMBUmNExZ7AoOcmJPqHZy82NOQlq3/5r94o36cQX0TGQrv4GQuzuWLPQ1p2Xq0ds3Z5+2DpEavUEhyVZLnk+xNcvOA/pOS3N/1P5rk3X2OJ2n8hg6FJCuAO4CrgQuBa5NcOGvY9cAPqurngNuAvxj2eJIWRp8rhfXA3qp6sareAO4DNs4asxG4p1v+B+DyJAMfbkiaDH1C4RzgpRnr+7u2gWOq6jDwCvDTg3aW5IYku5Ps/h9+3GNakvqYmAeNVbW5qtZV1bqVnLTY05HetvqEwhRw3oz1c7u2gWOSnAD8FPD9HseUNGZ9QuExYE2S85OcCGwCts0asw24rlv+KPAv5f/Vliba0F9eqqrDSW4CdgArgC1V9WySzwC7q2ob8LfA3yXZCxxiOjgkTbBM4j/cp+T08huN0vg8Wrt4tQ4N/CRwYh40SpoMhoKkhqEgqWEoSGoYCpIahoKkhqEgqWEoSGoYCpIahoKkhqEgqWEoSGoYCpIahoKkhqEgqWEoSGoYCpIahoKkRp8KUecl+dck30zybJI/HDDm0iSvJNnT/X2q33QljVufqtOHgT+qqieSrAIeT7Kzqr45a9xXq+qaHseRtICGvlKoqper6olu+YfAtzi6QpSkJWYkzxS6atK/CDw6oPsDSZ5M8kCSi46xD8vGSROgz+0DAEl+EvhH4JNV9eqs7ieAd1XVa0k2AP8MrBm0n6raDGyG6Z947zsvScPpdaWQZCXTgfDFqvqn2f1V9WpVvdYtbwdWJjmjzzEljVefTx/CdAWob1XVX80x5p1HSs8nWd8dz1qS0gTrc/vwy8BvA08n2dO1/SnwswBVdRfT9SM/nuQw8CNgk7UkpcnWp5bkI8DAslMzxtwO3D7sMSQtPL/RKKlhKEhqGAqSGoaCpIahIKlhKEhqGAqSGoaCpIahIKlhKEhqGAqSGoaCpIahIKlhKEhqGAqSGoaCpIahIKlhKEhq9A6FJPuSPN2Vhds9oD9J/ibJ3iRPJXl/32NKGp/edR86l1XV9+bou5rpWg9rgIuBO7tXSRNoIW4fNgJfqGlfB05NctYCHFfSEEYRCgU8lOTxJDcM6D8HeGnG+n4G1Jy0bJw0GUZx+3BJVU0lORPYmeS5qnr4/7sTy8ZJk6H3lUJVTXWvB4GtwPpZQ6aA82asn9u1SZpAfWtJnpxk1ZFl4ArgmVnDtgG/030K8UvAK1X1cp/jShqfvrcPq4GtXbnIE4AvVdWDSf4A/q903HZgA7AXeB34vZ7HlDRGvUKhql4EfmFA+10zlgv4RJ/jSFo4fqNRUsNQkNQwFCQ1DAVJDUNBUsNQkNQwFCQ1DAVJDUNBUsNQkNQwFCQ1DAVJDUNBUsNQkNQwFCQ1DAVJDUNBUsNQkNQYOhSSXNCVijvy92qST84ac2mSV2aM+VTvGUsaq6F/o7GqngfWAiRZwfTPtm8dMPSrVXXNsMeRtLBGdftwOfDvVfXtEe1P0iIZVShsAu6do+8DSZ5M8kCSi+bagWXjpMmQ6V9g77GD5ETgO8BFVXVgVt8pwFtV9VqSDcBfV9Wa4+3zlJxeF+fyXvOSNLdHaxev1qEM6hvFlcLVwBOzAwGgql6tqte65e3AyiRnjOCYksZkFKFwLXPcOiR5Z7ryUUnWd8f7/giOKWlMelWI6upHfgi4cUbbzJJxHwU+nuQw8CNgU/W9X5E0Vr2fKYyDzxSk8Rr3MwVJy4ihIKlhKEhqGAqSGoaCpIahIKlhKEhqGAqSGoaCpIahIKlhKEhqGAqSGoaCpIahIKlhKEhqGAqSGoaCpIahIKkxr1BIsiXJwSTPzGg7PcnOJC90r6fNse113ZgXklw3qolLGo/5XincDVw1q+1mYFdXx2FXt95IcjpwC3AxsB64Za7wkDQZ5hUKVfUwcGhW80bgnm75HuDDAza9EthZVYeq6gfATo4OF0kTpM8zhdVV9XK3/F1g9YAx5wAvzVjf37VJmlAjedDY1XLo9Vvx1pKUJkOfUDiQ5CyA7vXggDFTwHkz1s/t2o5SVZural1VrVvJST2mJamPPqGwDTjyacJ1wJcHjNkBXJHktO4B4xVdm6QJNd+PJO8FvgZckGR/kuuBzwIfSvIC8MFunSTrknweoKoOAX8OPNb9faZrkzShLBsnvQ1ZNk7SvBkKkhqGgqSGoSCpYShIahgKkhqGgqSGoSCpYShIahgKkhqGgqSGoSCpYShIahgKkhqGgqSGoSCpYShIahgKkhqGgqTGcUNhjjqSf5nkuSRPJdma5NQ5tt2X5Okke5LsHuG8JY3JfK4U7uboUm87gfdW1fuAfwP+5BjbX1ZVa6tq3XBTlLSQjhsKg+pIVtVDVXW4W/0600VeJC0Do3im8PvAA3P0FfBQkseT3HCsnVg2TpoMJ/TZOMmfAYeBL84x5JKqmkpyJrAzyXPdlcdRqmozsBmm6z70mZek4Q19pZDkd4FrgN+qOSrKVNVU93oQ2AqsH/Z4khbGUKGQ5Crgj4Ffr6rX5xhzcpJVR5aZriP5zKCxkibHfD6SHFRH8nZgFdO3BHuS3NWNPTvJ9m7T1cAjSZ4EvgF8paoeHMtZSBoZa0lKb0PWkpQ0b4aCpIahIKlhKEhqGAqSGoaCpIahIKlhKEhqGAqSGoaCpIahIKlhKEhqGAqSGoaCpIahIKlhKEhqGAqSGoaCpMawZeM+nWSq+33GPUk2zLHtVUmeT7I3yc2jnLik8Ri2bBzAbV05uLVVtX12Z5IVwB3A1cCFwLVJLuwzWUnjN1TZuHlaD+ytqher6g3gPmDjEPuRtID6PFO4qas6vSXJaQP6zwFemrG+v2sbyLJx0mQYNhTuBN4DrAVeBm7tO5Gq2lxV66pq3UpO6rs7SUMaKhSq6kBVvVlVbwGfY3A5uCngvBnr53ZtkibYsGXjzpqx+hEGl4N7DFiT5PwkJwKbgG3DHE/Swjlu1emubNylwBlJ9gO3AJcmWct0qfl9wI3d2LOBz1fVhqo6nOQmYAewAthSVc+O4yQkjY5l46S3IcvGSZo3Q0FSw1CQ1DAUJDUMBUkNQ0FSw1CQ1DAUJDUMBUkNQ0FSw1CQ1DAUJDUMBUkNQ0FSw1CQ1DAUJDUMBUkNQ0FSYz6/0bgFuAY4WFXv7druBy7ohpwK/GdVrR2w7T7gh8CbwOGqWjeSWUsam+OGAtNl424HvnCkoap+88hykluBV46x/WVV9b1hJyhpYR03FKrq4STvHtSXJMDHgF8b8bwkLZK+zxR+BThQVS/M0V/AQ0keT3LDsXZk2ThpMszn9uFYrgXuPUb/JVU1leRMYGeS57qCtUepqs3AZpj+ifee85I0pKGvFJKcAPwGcP9cY6pqqns9CGxlcHk5SROkz+3DB4Hnqmr/oM4kJydZdWQZuILB5eUkTZDjhkJXNu5rwAVJ9ie5vuvaxKxbhyRnJ9nera4GHknyJPAN4CtV9eDopi5pHCwbJ70NWTZO0rwZCpIahoKkhqEgqWEoSGr0/UbjWPz8+15nx4498xp75dlrxzoXaSnZ8Z098xq3/srX5+zzSkFSw1CQ1DAUJDUMBUkNQ0FSw1CQ1DAUJDUMBUkNQ0FSw1CQ1JjIH1lJ8h/At2c1nwEsx/oRy/W8YPme23I4r3dV1c8M6pjIUBgkye7lWGFquZ4XLN9zW67ndYS3D5IahoKkxlIKhc2LPYExWa7nBcv33JbreQFL6JmCpIWxlK4UJC0AQ0FSY0mEQpKrkjyfZG+Smxd7PqOSZF+Sp5PsSbJ7sefTR5ItSQ4meWZG2+lJdiZ5oXs9bTHnOIw5zuvTSaa6921Pkg2LOcdRm/hQSLICuAO4GrgQuDbJhYs7q5G6rKrWLoPPve8GrprVdjOwq6rWALu69aXmbo4+L4DbuvdtbVVtH9C/ZE18KDBdqXpvVb1YVW8A9wEbF3lOmqWqHgYOzWreCNzTLd8DfHgh5zQKc5zXsrYUQuEc4KUZ6/u7tuWggIeSPJ7khsWezBisrqqXu+XvMl10eLm4KclT3e3FkrstOpalEArL2SVV9X6mb40+keRXF3tC41LTn30vl8+/7wTeA6wFXgZuXdTZjNhSCIUp4LwZ6+d2bUteVU11rweBrUzfKi0nB5KcBdC9Hlzk+YxEVR2oqjer6i3gcyyz920phMJjwJok5yc5EdgEbFvkOfWW5OQkq44sA1cAzxx7qyVnG3Bdt3wd8OVFnMvIHAm6zkdYZu/bRFaImqmqDie5CdgBrAC2VNWzizytUVgNbE0C0+/Dl6rqwcWd0vCS3AtcCpyRZD9wC/BZ4O+TXM/0f4X/2OLNcDhznNelSdYyfTu0D7hxseY3Dn7NWVJjKdw+SFpAhoKkhqEgqWEoSGoYCpIahoKkhqEgqfG/uHzLBjDxX7wAAAAASUVORK5CYII=",
|
|
"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[19,0] = 2000\n",
|
|
"C_t[19,19] = 2000\n",
|
|
"\n",
|
|
"\n",
|
|
"plt.imshow(C_t, vmin=0, vmax=2000)\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 41,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# calculate mean between two cells\n",
|
|
"def alpha_interblock(alpha1, alpha2, harmonic=True):\n",
|
|
" if not harmonic:\n",
|
|
" return 0.5 * (alpha1 + alpha2)\n",
|
|
" else:\n",
|
|
" return 2 / ((1/alpha1) + (1/alpha2))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 217,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def sub(delta_t, delta_x):\n",
|
|
" return delta_t / (2 * delta_x**2)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 260,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"##############\n",
|
|
"### CLOSED ###\n",
|
|
"##############\n",
|
|
"\n",
|
|
"# creates the coeffiecient matrix for a given row\n",
|
|
"def create_coeff_matrix(alpha_x, row_index, s_x):\n",
|
|
" # alias\n",
|
|
" r = row_index\n",
|
|
"\n",
|
|
" # a square matrix of column^2 dimension for the coefficients\n",
|
|
" cm = np.zeros((alpha_x.shape[1],alpha_x.shape[1]))\n",
|
|
"\n",
|
|
" # top left\n",
|
|
" cm[0,0] = 1 + s_x * (alpha_interblock(alpha_x[r,0], alpha_x[r,1]) + 2 * alpha_x[r,0]) # alpha_x[0,0]\n",
|
|
" cm[0,1] = -s_x * alpha_interblock(alpha_x[r, 0], alpha_x[r, 1])\n",
|
|
"\n",
|
|
" # inner cells\n",
|
|
" for i_cm in range(1, cm.shape[0]-1):\n",
|
|
" j_alpha = i_cm #- 1 # ???????\n",
|
|
" cm[i_cm,i_cm-1] = -s_x * alpha_interblock(alpha_x[r, j_alpha-1], alpha_x[r, j_alpha])\n",
|
|
" cm[i_cm,i_cm] = 1 + s_x * (alpha_interblock(alpha_x[r,j_alpha], alpha_x[r, j_alpha+1]) \\\n",
|
|
" + alpha_interblock(alpha_x[r, j_alpha], alpha_x[r, j_alpha-1]))\n",
|
|
" cm[i_cm,i_cm+1] = -s_x * alpha_interblock(alpha_x[r, j_alpha], alpha_x[r, j_alpha+1])\n",
|
|
"\n",
|
|
"\n",
|
|
" # bottom right\n",
|
|
" n = cm.shape[0]-1\n",
|
|
" cm[n,n-1] = -s_x * alpha_interblock(alpha_x[r, n-1], alpha_x[r, n]) ## changed to - !!!\n",
|
|
" cm[n,n] = 1 + s_x * (alpha_interblock(alpha_x[r, n], alpha_x[r, n-1]) + 2 * alpha_x[r,n]) # alpha_x[n,n]? -> changed to [r,n]!!! and changed to +alpha_interblock\n",
|
|
"\n",
|
|
" # print(cm)\n",
|
|
" return cm\n",
|
|
"\n",
|
|
"# creates the solution vector for a given row\n",
|
|
"# s_x, alpha_x, row_index?\n",
|
|
"def create_solution_vector(concentrations, alpha_x, alpha_y, row_index, s_x, s_y, bc_left, bc_right, bc_top, bc_bottom):\n",
|
|
"\n",
|
|
" sv = np.zeros(concentrations.shape[1])\n",
|
|
" # print(concentrations.shape[1])\n",
|
|
" # print(sv.shape[0])\n",
|
|
" # print(row_index)\n",
|
|
"\n",
|
|
" # inner rows\n",
|
|
" if row_index > 0 and row_index < concentrations.shape[0]-1:\n",
|
|
" for i in range(0, sv.shape[0]):\n",
|
|
" sv[i] = s_y * alpha_interblock(alpha_y[row_index+1,i],alpha_y[row_index,i]) * concentrations[row_index+1,i] \\\n",
|
|
" + (\n",
|
|
" 1 - s_y * (\n",
|
|
" alpha_interblock(alpha_y[row_index+1,i], alpha_y[row_index,i]) \\\n",
|
|
" + alpha_interblock(alpha_y[row_index,i], alpha_y[row_index-1,i])\n",
|
|
" )\n",
|
|
" ) * concentrations[row_index,i] \\\n",
|
|
" + s_y * alpha_interblock(alpha_y[row_index,i], alpha_y[row_index-1,i]) * concentrations[row_index-1,i]\n",
|
|
"\n",
|
|
" # first row\n",
|
|
" # TODO check formula for alpha boundary cases\n",
|
|
" # added factor of 2 in front of alpha_y!!!\n",
|
|
" if row_index == 0:\n",
|
|
" for i in range(0, sv.shape[0]):\n",
|
|
" sv[i] = s_y * alpha_interblock(alpha_y[row_index+1,i],alpha_y[row_index,i]) * concentrations[row_index+1,i] \\\n",
|
|
" + (\n",
|
|
" 1 - s_y * (\n",
|
|
" alpha_interblock(alpha_y[row_index+1,i], alpha_y[row_index,i]) \\\n",
|
|
" + 2 * alpha_y[row_index,i]\n",
|
|
" )\n",
|
|
" ) * concentrations[row_index,i] \\\n",
|
|
" + s_y * alpha_y[row_index,i] * bc_top\n",
|
|
"\n",
|
|
"\n",
|
|
" # last row\n",
|
|
" # TODO check formula for alpha boundary cases\n",
|
|
" # added factor of 2 in front of alpha_y!!!\n",
|
|
" if row_index == sv.shape[0]-1:\n",
|
|
" for i in range(0, sv.shape[0]):\n",
|
|
" i = i\n",
|
|
" sv[i] = s_y * alpha_y[row_index,i]*bc_bottom \\\n",
|
|
" + (\n",
|
|
" 1 - s_y * (\n",
|
|
" 2 * alpha_y[row_index,i] \\\n",
|
|
" + alpha_interblock(alpha_y[row_index,i], alpha_y[row_index-1,i])\n",
|
|
" )\n",
|
|
" ) * concentrations[row_index,i] \\\n",
|
|
" + s_y * alpha_interblock(alpha_y[row_index,i], alpha_y[row_index-1,i]) * concentrations[row_index-1,i]\n",
|
|
"\n",
|
|
"\n",
|
|
" # first column (additional fixed concentration change from perpendicular dimension)\n",
|
|
" # TODO check correctness of additional summand\n",
|
|
" i = 0\n",
|
|
" sv[0] += 2 * s_x * alpha_x[row_index,0] * bc_left # alpha_x[row_index,0]?\n",
|
|
"\n",
|
|
" # last column (additional fixed concentration change from perpendicular direction)\n",
|
|
" # TODO check correctness of additional summand\n",
|
|
" n = sv.shape[0]-1\n",
|
|
" sv[n] += 2 * s_x * alpha_x[row_index,n] * bc_right\n",
|
|
"\n",
|
|
" # print(sv)\n",
|
|
" return sv\n",
|
|
"\n",
|
|
"def calc_sweep(concentrations, alpha_x, alpha_y, s_x, s_y, bc_left, bc_right, bc_top, bc_bottom):\n",
|
|
"\n",
|
|
" concentrations_t1 = np.zeros((concentrations.shape[0], concentrations.shape[1]))\n",
|
|
"\n",
|
|
" for i in range(concentrations.shape[0]):\n",
|
|
" A = create_coeff_matrix(alpha_x, i, s_x)\n",
|
|
" # print(i)\n",
|
|
" b = create_solution_vector(concentrations, alpha_x, alpha_y, i, s_x, s_y, bc_left, bc_right, bc_top, bc_bottom)\n",
|
|
" concentrations_t1[i,:] = np.linalg.solve(A,b)\n",
|
|
"\n",
|
|
" # print(concentrations_t1)\n",
|
|
" return concentrations_t1\n",
|
|
"\n",
|
|
"def calc_iteration(concentrations, alpha_x, alpha_y, s_x, s_y, bc_left, bc_right, bc_top, bc_bottom):\n",
|
|
" concentrations = calc_sweep(concentrations, alpha_x, alpha_y, s_x, s_y, bc_left, bc_right, bc_top, bc_bottom)\n",
|
|
" # concentrations = calc_sweep(concentrations, alpha_x, alpha_y, s_x, s_y, bc_left, bc_right, bc_top, bc_bottom, True)\n",
|
|
" concentrations = calc_sweep(concentrations.transpose(), alpha_y.transpose(), alpha_x.transpose(), s_y, s_x, bc_top, bc_bottom, bc_left, bc_right).transpose()\n",
|
|
" \n",
|
|
" # print(concentrations)\n",
|
|
" return concentrations\n",
|
|
"\n",
|
|
"\n",
|
|
"def sim_constant(concentrations, alpha_x, alpha_y, iterations, delta_t, delta_x, delta_y, bc_left, bc_right, bc_top, bc_bottom):\n",
|
|
" s_x = sub(delta_t, delta_x)\n",
|
|
" s_y = sub(delta_t, delta_y)\n",
|
|
" for i in range(iterations):\n",
|
|
" concentrations = calc_iteration(concentrations, alpha_x, alpha_y, s_x, s_y, bc_left, bc_right, bc_top, bc_bottom)\n",
|
|
"\n",
|
|
" return concentrations\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 279,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"################\n",
|
|
"### CONSTANT ###\n",
|
|
"################\n",
|
|
"\n",
|
|
"# creates the coeffiecient matrix for a given row\n",
|
|
"def create_coeff_matrix(alpha_x, row_index, s_x):\n",
|
|
" # alias\n",
|
|
" r = row_index\n",
|
|
"\n",
|
|
" # a square matrix of column^2 dimension for the coefficients\n",
|
|
" cm = np.zeros((alpha_x.shape[1],alpha_x.shape[1]))\n",
|
|
"\n",
|
|
" # top left\n",
|
|
" cm[0,0] = 1 + s_x * (alpha_interblock(alpha_x[r,0], alpha_x[r,1])) # alpha_x[0,0]\n",
|
|
" cm[0,1] = -s_x * alpha_interblock(alpha_x[r, 0], alpha_x[r, 1])\n",
|
|
"\n",
|
|
" # inner cells\n",
|
|
" for i_cm in range(1, cm.shape[0]-1):\n",
|
|
" j_alpha = i_cm #- 1 # ???????\n",
|
|
" cm[i_cm,i_cm-1] = -s_x * alpha_interblock(alpha_x[r, j_alpha-1], alpha_x[r, j_alpha])\n",
|
|
" cm[i_cm,i_cm] = 1 + s_x * (alpha_interblock(alpha_x[r,j_alpha], alpha_x[r, j_alpha+1]) \\\n",
|
|
" + alpha_interblock(alpha_x[r, j_alpha], alpha_x[r, j_alpha-1]))\n",
|
|
" cm[i_cm,i_cm+1] = -s_x * alpha_interblock(alpha_x[r, j_alpha], alpha_x[r, j_alpha+1])\n",
|
|
"\n",
|
|
"\n",
|
|
" # bottom right\n",
|
|
" n = cm.shape[0]-1\n",
|
|
" cm[n,n-1] = -s_x * alpha_interblock(alpha_x[r, n-1], alpha_x[r, n])\n",
|
|
" cm[n,n] = 1 + s_x * (alpha_interblock(alpha_x[r, n], alpha_x[r, n-1])) # alpha_x[n,n]? \n",
|
|
"\n",
|
|
" # print(cm)\n",
|
|
" return cm\n",
|
|
"\n",
|
|
"# creates the solution vector for a given row\n",
|
|
"# s_x, alpha_x, row_index?\n",
|
|
"def create_solution_vector(concentrations, alpha_x, alpha_y, row_index, s_x, s_y, bc_left, bc_right, bc_top, bc_bottom):\n",
|
|
"\n",
|
|
" sv = np.zeros(concentrations.shape[1])\n",
|
|
" # print(concentrations.shape[1])\n",
|
|
" # print(sv.shape[0])\n",
|
|
" # print(row_index)\n",
|
|
"\n",
|
|
" # inner rows\n",
|
|
" if row_index > 0 and row_index < concentrations.shape[0]-1:\n",
|
|
" for i in range(0, sv.shape[0]):\n",
|
|
" sv[i] = s_y * alpha_interblock(alpha_y[row_index+1,i],alpha_y[row_index,i]) * concentrations[row_index+1,i] \\\n",
|
|
" + (\n",
|
|
" 1 - s_y * (\n",
|
|
" alpha_interblock(alpha_y[row_index+1,i], alpha_y[row_index,i]) \\\n",
|
|
" + alpha_interblock(alpha_y[row_index,i], alpha_y[row_index-1,i])\n",
|
|
" )\n",
|
|
" ) * concentrations[row_index,i] \\\n",
|
|
" + s_y * alpha_interblock(alpha_y[row_index,i], alpha_y[row_index-1,i]) * concentrations[row_index-1,i]\n",
|
|
"\n",
|
|
" # first row\n",
|
|
" # TODO check formula for alpha boundary cases\n",
|
|
" # added factor of 2 in front of alpha_y!!!\n",
|
|
" if row_index == 0:\n",
|
|
" for i in range(0, sv.shape[0]):\n",
|
|
" sv[i] = s_y * alpha_interblock(alpha_y[row_index+1,i],alpha_y[row_index,i]) * concentrations[row_index+1,i] \\\n",
|
|
" + (\n",
|
|
" 1 - s_y * (\n",
|
|
" alpha_interblock(alpha_y[row_index+1,i], alpha_y[row_index,i]) \\\n",
|
|
" )\n",
|
|
" ) * concentrations[row_index,i]\n",
|
|
"\n",
|
|
"\n",
|
|
" # last row\n",
|
|
" # TODO check formula for alpha boundary cases\n",
|
|
" # added factor of 2 in front of alpha_y!!!\n",
|
|
" if row_index == sv.shape[0]-1:\n",
|
|
" for i in range(0, sv.shape[0]):\n",
|
|
" i = i\n",
|
|
" sv[i] = (\n",
|
|
" 1 - s_y * (\n",
|
|
" alpha_interblock(alpha_y[row_index,i], alpha_y[row_index-1,i])\n",
|
|
" )\n",
|
|
" ) * concentrations[row_index,i] \\\n",
|
|
" + s_y * alpha_interblock(alpha_y[row_index,i], alpha_y[row_index-1,i]) * concentrations[row_index-1,i]\n",
|
|
"\n",
|
|
" # print(sv)\n",
|
|
" return sv\n",
|
|
"\n",
|
|
"def calc_sweep(concentrations, alpha_x, alpha_y, s_x, s_y, bc_left, bc_right, bc_top, bc_bottom):\n",
|
|
"\n",
|
|
" concentrations_t1 = np.zeros((concentrations.shape[0], concentrations.shape[1]))\n",
|
|
"\n",
|
|
" for i in range(concentrations.shape[0]):\n",
|
|
" A = create_coeff_matrix(alpha_x, i, s_x)\n",
|
|
" # print(i)\n",
|
|
" b = create_solution_vector(concentrations, alpha_x, alpha_y, i, s_x, s_y, bc_left, bc_right, bc_top, bc_bottom)\n",
|
|
" concentrations_t1[i,:] = np.linalg.solve(A,b)\n",
|
|
"\n",
|
|
" # print(concentrations_t1)\n",
|
|
" return concentrations_t1\n",
|
|
"\n",
|
|
"def calc_iteration(concentrations, alpha_x, alpha_y, s_x, s_y, bc_left, bc_right, bc_top, bc_bottom):\n",
|
|
" concentrations = calc_sweep(concentrations, alpha_x, alpha_y, s_x, s_y, bc_left, bc_right, bc_top, bc_bottom)\n",
|
|
" # concentrations = calc_sweep(concentrations, alpha_x, alpha_y, s_x, s_y, bc_left, bc_right, bc_top, bc_bottom, True)\n",
|
|
" concentrations = calc_sweep(concentrations.transpose(), alpha_y.transpose(), alpha_x.transpose(), s_y, s_x, bc_top, bc_bottom, bc_left, bc_right).transpose()\n",
|
|
" \n",
|
|
" # print(concentrations)\n",
|
|
" return concentrations\n",
|
|
"\n",
|
|
"\n",
|
|
"def sim_closed(concentrations, alpha_x, alpha_y, iterations, delta_t, delta_x, delta_y, bc_left, bc_right, bc_top, bc_bottom):\n",
|
|
" s_x = sub(delta_t, delta_x)\n",
|
|
" s_y = sub(delta_t, delta_y)\n",
|
|
" for i in range(iterations):\n",
|
|
" concentrations = calc_iteration(concentrations, alpha_x, alpha_y, s_x, s_y, bc_left, bc_right, bc_top, bc_bottom)\n",
|
|
"\n",
|
|
" return concentrations\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 280,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Test Concentrations t=0:\n",
|
|
"[[2000. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
|
|
" 0. 0. 0. 0. 0. 0. 0. 0.]\n",
|
|
" [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
|
|
" 0. 0. 0. 0. 0. 0. 0. 0.]\n",
|
|
" [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
|
|
" 0. 0. 0. 0. 0. 0. 0. 0.]\n",
|
|
" [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
|
|
" 0. 0. 0. 0. 0. 0. 0. 0.]\n",
|
|
" [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
|
|
" 0. 0. 0. 0. 0. 0. 0. 0.]\n",
|
|
" [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
|
|
" 0. 0. 0. 0. 0. 0. 0. 0.]\n",
|
|
" [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
|
|
" 0. 0. 0. 0. 0. 0. 0. 0.]\n",
|
|
" [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
|
|
" 0. 0. 0. 0. 0. 0. 0. 0.]\n",
|
|
" [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
|
|
" 0. 0. 0. 0. 0. 0. 0. 0.]\n",
|
|
" [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
|
|
" 0. 0. 0. 0. 0. 0. 0. 0.]\n",
|
|
" [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
|
|
" 0. 0. 0. 0. 0. 0. 0. 0.]\n",
|
|
" [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
|
|
" 0. 0. 0. 0. 0. 0. 0. 0.]\n",
|
|
" [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
|
|
" 0. 0. 0. 0. 0. 0. 0. 0.]\n",
|
|
" [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
|
|
" 0. 0. 0. 0. 0. 0. 0. 0.]\n",
|
|
" [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
|
|
" 0. 0. 0. 0. 0. 0. 0. 0.]\n",
|
|
" [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
|
|
" 0. 0. 0. 0. 0. 0. 0. 0.]\n",
|
|
" [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
|
|
" 0. 0. 0. 0. 0. 0. 0. 0.]\n",
|
|
" [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
|
|
" 0. 0. 0. 0. 0. 0. 0. 0.]\n",
|
|
" [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
|
|
" 0. 0. 0. 0. 0. 0. 0. 0.]\n",
|
|
" [ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
|
|
" 0. 0. 0. 0. 0. 0. 0. 0.]]\n",
|
|
"\n",
|
|
"Test Alpha:\n",
|
|
"[[1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n",
|
|
" [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n",
|
|
" [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n",
|
|
" [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n",
|
|
" [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n",
|
|
" [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n",
|
|
" [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n",
|
|
" [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n",
|
|
" [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n",
|
|
" [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n",
|
|
" [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n",
|
|
" [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n",
|
|
" [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n",
|
|
" [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n",
|
|
" [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n",
|
|
" [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n",
|
|
" [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n",
|
|
" [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n",
|
|
" [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n",
|
|
" [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]]\n",
|
|
"\n",
|
|
"Test Concentrations t=t+100:\n",
|
|
"[[6.28591428e+01 5.97567801e+01 5.40111841e+01 4.64274274e+01\n",
|
|
" 3.79691932e+01 2.95580781e+01 2.19168076e+01 1.54897130e+01\n",
|
|
" 1.04428201e+01 6.72160139e+00 4.13431473e+00 2.43233389e+00\n",
|
|
" 1.37011777e+00 7.39680705e-01 3.83120860e-01 1.90615726e-01\n",
|
|
" 9.12968457e-02 4.24222400e-02 1.98821745e-02 1.12119045e-02]\n",
|
|
" [5.97567801e+01 5.68075320e+01 5.13455053e+01 4.41360389e+01\n",
|
|
" 3.60952541e+01 2.80992628e+01 2.08351211e+01 1.47252306e+01\n",
|
|
" 9.92742306e+00 6.38986213e+00 3.93026893e+00 2.31228800e+00\n",
|
|
" 1.30249670e+00 7.03174355e-01 3.64212236e-01 1.81208039e-01\n",
|
|
" 8.67909629e-02 4.03285243e-02 1.89009056e-02 1.06585500e-02]\n",
|
|
" [5.40111841e+01 5.13455053e+01 4.64086509e+01 3.98923724e+01\n",
|
|
" 3.26247065e+01 2.53975274e+01 1.88318306e+01 1.33094043e+01\n",
|
|
" 8.97290440e+00 5.77547886e+00 3.55237479e+00 2.08996222e+00\n",
|
|
" 1.17726205e+00 6.35564358e-01 3.29193342e-01 1.63784942e-01\n",
|
|
" 7.84460386e-02 3.64509491e-02 1.70835894e-02 9.63373367e-03]\n",
|
|
" [4.64274274e+01 4.41360389e+01 3.98923724e+01 3.42910501e+01\n",
|
|
" 2.80438435e+01 2.18314387e+01 1.61876371e+01 1.14406194e+01\n",
|
|
" 7.71301119e+00 4.96453891e+00 3.05358279e+00 1.79650883e+00\n",
|
|
" 1.01196167e+00 5.46324221e-01 2.82971023e-01 1.40787758e-01\n",
|
|
" 6.74313630e-02 3.13328401e-02 1.46848680e-02 8.28105286e-03]\n",
|
|
" [3.79691932e+01 3.60952541e+01 3.26247065e+01 2.80438435e+01\n",
|
|
" 2.29347645e+01 1.78541471e+01 1.32385436e+01 9.35634629e+00\n",
|
|
" 6.30784062e+00 4.06009008e+00 2.49727546e+00 1.46921755e+00\n",
|
|
" 8.27600633e-01 4.46793869e-01 2.31418841e-01 1.15138785e-01\n",
|
|
" 5.51465932e-02 2.56245656e-02 1.20095517e-02 6.77239541e-03]\n",
|
|
" [2.95580781e+01 2.80992628e+01 2.53975274e+01 2.18314387e+01\n",
|
|
" 1.78541471e+01 1.38990120e+01 1.03058789e+01 7.28368424e+00\n",
|
|
" 4.91049795e+00 3.16067973e+00 1.94406720e+00 1.14374953e+00\n",
|
|
" 6.44266629e-01 3.47817980e-01 1.80153846e-01 8.96326973e-02\n",
|
|
" 4.29302593e-02 1.99480907e-02 9.34913903e-03 5.27214237e-03]\n",
|
|
" [2.19168076e+01 2.08351211e+01 1.88318306e+01 1.61876371e+01\n",
|
|
" 1.32385436e+01 1.03058789e+01 7.64163230e+00 5.40072685e+00\n",
|
|
" 3.64104995e+00 2.34358978e+00 1.44149246e+00 8.48070648e-01\n",
|
|
" 4.77712648e-01 2.57901062e-01 1.33580985e-01 6.64611067e-02\n",
|
|
" 3.18320504e-02 1.47911669e-02 6.93222614e-03 3.90920308e-03]\n",
|
|
" [1.54897130e+01 1.47252306e+01 1.33094043e+01 1.14406194e+01\n",
|
|
" 9.35634629e+00 7.28368424e+00 5.40072685e+00 3.81696598e+00\n",
|
|
" 2.57331359e+00 1.65633306e+00 1.01877541e+00 5.99374289e-01\n",
|
|
" 3.37623616e-01 1.82271685e-01 9.44084172e-02 4.69714152e-02\n",
|
|
" 2.24973151e-02 1.04536635e-02 4.89935375e-03 2.76283093e-03]\n",
|
|
" [1.04428201e+01 9.92742306e+00 8.97290440e+00 7.71301119e+00\n",
|
|
" 6.30784062e+00 4.91049795e+00 3.64104995e+00 2.57331359e+00\n",
|
|
" 1.73487080e+00 1.11666292e+00 6.86835725e-01 4.04084817e-01\n",
|
|
" 2.27618334e-01 1.22883517e-01 6.36480556e-02 3.16670837e-02\n",
|
|
" 1.51671896e-02 7.04762745e-03 3.30303535e-03 1.86263918e-03]\n",
|
|
" [6.72160139e+00 6.38986213e+00 5.77547886e+00 4.96453891e+00\n",
|
|
" 4.06009008e+00 3.16067973e+00 2.34358978e+00 1.65633306e+00\n",
|
|
" 1.11666292e+00 7.18748669e-01 4.42087091e-01 2.60092297e-01\n",
|
|
" 1.46508290e-01 7.90949197e-02 4.09675600e-02 2.03827617e-02\n",
|
|
" 9.76247810e-03 4.53625956e-03 2.12602408e-03 1.19890202e-03]\n",
|
|
" [4.13431473e+00 3.93026893e+00 3.55237479e+00 3.05358279e+00\n",
|
|
" 2.49727546e+00 1.94406720e+00 1.44149246e+00 1.01877541e+00\n",
|
|
" 6.86835725e-01 4.42087091e-01 2.71918412e-01 1.59977266e-01\n",
|
|
" 9.01141476e-02 4.86496108e-02 2.51982789e-02 1.25370053e-02\n",
|
|
" 6.00469362e-03 2.79015723e-03 1.30767241e-03 7.37419252e-04]\n",
|
|
" [2.43233389e+00 2.31228800e+00 2.08996222e+00 1.79650883e+00\n",
|
|
" 1.46921755e+00 1.14374953e+00 8.48070648e-01 5.99374289e-01\n",
|
|
" 4.04084817e-01 2.60092297e-01 1.59977266e-01 9.41191353e-02\n",
|
|
" 5.30166932e-02 2.86219374e-02 1.48248578e-02 7.37587358e-03\n",
|
|
" 3.53273052e-03 1.64152815e-03 7.69340539e-04 4.33844533e-04]\n",
|
|
" [1.37011777e+00 1.30249670e+00 1.17726205e+00 1.01196167e+00\n",
|
|
" 8.27600633e-01 6.44266629e-01 4.77712648e-01 3.37623616e-01\n",
|
|
" 2.27618334e-01 1.46508290e-01 9.01141476e-02 5.30166932e-02\n",
|
|
" 2.98639564e-02 1.61225501e-02 8.35074541e-03 4.15478134e-03\n",
|
|
" 1.98996399e-03 9.24662067e-04 4.33364493e-04 2.44381787e-04]\n",
|
|
" [7.39680705e-01 7.03174355e-01 6.35564358e-01 5.46324221e-01\n",
|
|
" 4.46793869e-01 3.47817980e-01 2.57901062e-01 1.82271685e-01\n",
|
|
" 1.22883517e-01 7.90949197e-02 4.86496108e-02 2.86219374e-02\n",
|
|
" 1.61225501e-02 8.70402492e-03 4.50828782e-03 2.24302732e-03\n",
|
|
" 1.07431492e-03 4.99194087e-04 2.33958978e-04 1.31933544e-04]\n",
|
|
" [3.83120860e-01 3.64212236e-01 3.29193342e-01 2.82971023e-01\n",
|
|
" 2.31418841e-01 1.80153846e-01 1.33580985e-01 9.44084172e-02\n",
|
|
" 6.36480556e-02 4.09675600e-02 2.51982789e-02 1.48248578e-02\n",
|
|
" 8.35074541e-03 4.50828782e-03 2.33508742e-03 1.16178582e-03\n",
|
|
" 5.56446119e-04 2.58559763e-04 1.21180077e-04 6.83355566e-05]\n",
|
|
" [1.90615726e-01 1.81208039e-01 1.63784942e-01 1.40787758e-01\n",
|
|
" 1.15138785e-01 8.96326973e-02 6.64611067e-02 4.69714152e-02\n",
|
|
" 3.16670837e-02 2.03827617e-02 1.25370053e-02 7.37587358e-03\n",
|
|
" 4.15478134e-03 2.24302732e-03 1.16178582e-03 5.78028166e-04\n",
|
|
" 2.76850968e-04 1.28642322e-04 6.02912313e-05 3.39992757e-05]\n",
|
|
" [9.12968457e-02 8.67909629e-02 7.84460386e-02 6.74313630e-02\n",
|
|
" 5.51465932e-02 4.29302593e-02 3.18320504e-02 2.24973151e-02\n",
|
|
" 1.51671896e-02 9.76247810e-03 6.00469362e-03 3.53273052e-03\n",
|
|
" 1.98996399e-03 1.07431492e-03 5.56446119e-04 2.76850968e-04\n",
|
|
" 1.32599868e-04 6.16142144e-05 2.88769419e-05 1.62842106e-05]\n",
|
|
" [4.24222400e-02 4.03285243e-02 3.64509491e-02 3.13328401e-02\n",
|
|
" 2.56245656e-02 1.99480907e-02 1.47911669e-02 1.04536635e-02\n",
|
|
" 7.04762745e-03 4.53625956e-03 2.79015723e-03 1.64152815e-03\n",
|
|
" 9.24662067e-04 4.99194087e-04 2.58559763e-04 1.28642322e-04\n",
|
|
" 6.16142144e-05 2.86298280e-05 1.34180382e-05 7.56666548e-06]\n",
|
|
" [1.98821745e-02 1.89009056e-02 1.70835894e-02 1.46848680e-02\n",
|
|
" 1.20095517e-02 9.34913903e-03 6.93222614e-03 4.89935375e-03\n",
|
|
" 3.30303535e-03 2.12602408e-03 1.30767241e-03 7.69340539e-04\n",
|
|
" 4.33364493e-04 2.33958978e-04 1.21180077e-04 6.02912313e-05\n",
|
|
" 2.88769419e-05 1.34180382e-05 6.28867725e-06 3.54629466e-06]\n",
|
|
" [1.12119045e-02 1.06585500e-02 9.63373367e-03 8.28105286e-03\n",
|
|
" 6.77239541e-03 5.27214237e-03 3.90920308e-03 2.76283093e-03\n",
|
|
" 1.86263918e-03 1.19890202e-03 7.37419252e-04 4.33844533e-04\n",
|
|
" 2.44381787e-04 1.31933544e-04 6.83355566e-05 3.39992757e-05\n",
|
|
" 1.62842106e-05 7.56666548e-06 3.54629466e-06 1.99981734e-06]]\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"row = 20\n",
|
|
"col = 20\n",
|
|
"domain_row = row\n",
|
|
"domain_col = col\n",
|
|
"iterations = 100\n",
|
|
"delta_t = 0.1\n",
|
|
"delta_x = domain_col / col\n",
|
|
"delta_y = domain_row / row\n",
|
|
"bc_left = 0\n",
|
|
"bc_right = 0\n",
|
|
"bc_top = 0\n",
|
|
"bc_bottom = 0\n",
|
|
"\n",
|
|
"\n",
|
|
"testConcentrations = np.zeros((row,col))\n",
|
|
"testConcentrations[0,0] = 2000\n",
|
|
"print(\"Test Concentrations t=0:\")\n",
|
|
"print(testConcentrations)\n",
|
|
"\n",
|
|
"# seed = np.random.seed(42)\n",
|
|
"# testAlpha = np.random.random((row,col)) * 10\n",
|
|
"testAlpha= np.ones((row,col))\n",
|
|
"print(\"\\nTest Alpha:\")\n",
|
|
"print(testAlpha)\n",
|
|
"\n",
|
|
"testConcentrations = sim_closed(testConcentrations, testAlpha, testAlpha, iterations, delta_t, delta_x, delta_y, bc_left, bc_right, bc_top, bc_bottom)\n",
|
|
"print(\"\\nTest Concentrations t=t+\" + str(iterations) + \":\")\n",
|
|
"print(testConcentrations)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 281,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAfUAAAGdCAYAAAAL7+omAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA3wUlEQVR4nO3de3xU5Z3H8e8JkAlgEouFTCIBIxtBLqUUkBCVi5TQWKldtKBYLvVSXdCK1Bc1Ul+N3W6iVGlWUVxd5LIWsV0u8lqsEl6SIAUsSGIpUsQ1hVQz5gULCbdc59k/aEbHXGAyZ8J4zufd1/Oqc87zPPObk9HfPM95zjmWMcYIAAB85cVc7AAAAIA9SOoAADgESR0AAIcgqQMA4BAkdQAAHIKkDgCAQ5DUAQBwCJI6AAAO0fliB2AXv9+vTz/9VPHx8bIs62KHAwAIgTFGJ0+eVEpKimJiIjferKmpUV1dnS19xcbGKi4uzpa+7OKYpP7pp58qNTX1YocBAAhDeXm5evfuHZG+a2pqlNb3EvkqG23pz+v1qqysLKoSu2OSenx8vCTp8N4rlHBJeL/yyhtO2RGSJGn72T629LPl/wbZ0o8klZTb9C9MeTd7+pHU7e/2zK5cUtFgSz+S1PWzGlv66fR/p23pR5J00p6+zNmztvQjSaau3p6OGu35D60kGbv64i7aHaZB9dquNwL/LY+Euro6+SobVfZeXyXEh5cnqk/6lTb8sOrq6kjqkdA05Z5wSUzYf6z4Bvumfrp2sucQd6mNtaUfSYrpZtMX0MYvcqdYe5J65y72JfXONv3b0amTfTEpxp4EaiwbE6hdZ7ss+/69M7b1RVLvMP841B1x+jQhPvw8Ea0ck9QBALgQjcavxjB/rzUavz3B2IykDgBwFb+M/GHOwoTbPlIiNv/w/PPPKy0tTXFxcRo+fLjeeeedNusXFxdr+PDhiouL05VXXqkXXnghUqEBAFzMb9P/olFEkvprr72mefPmaeHChSopKdH111+v7OxsHTlypMX6ZWVluvHGG3X99derpKREjz76qH7yk59o7dq1kQgPAABHikhSX7x4se666y7dfffduvrqq1VQUKDU1FQtXbq0xfovvPCC+vTpo4KCAl199dW6++67deedd+qpp56KRHgAABdrNMaWEo1sT+p1dXV67733lJWVFbQ9KytLO3bsaLHNzp07m9WfNGmS9uzZo/r6llf71tbWqrq6OqgAAHA+TefUwy3RyPakfvToUTU2NiopKSloe1JSknw+X4ttfD5fi/UbGhp09OjRFtvk5+crMTExULjxDADA7SK2UO7L1xoaY9q8/rCl+i1tb5KTk6OqqqpAKS8vDzNiAIAb+GXUGGaJ1pG67Ze0ff3rX1enTp2ajcorKyubjcabeL3eFut37txZl112WYttPB6PPB6PPUEDAFyDS9pCEBsbq+HDh6uwsDBoe2FhoTIzM1tsM3r06Gb1N2/erBEjRqhLly52hwgAQIf65JNP9MMf/lCXXXaZunXrpm9+85t67733AvuNMcrNzVVKSoq6du2qcePGaf/+/SG/T0Sm3+fPn6///M//1Msvv6wDBw7ooYce0pEjR3TfffdJOjd1PnPmzED9++67T4cPH9b8+fN14MABvfzyy1q2bJkefvjhSIQHAHCxjl79fvz4cV177bXq0qWL/vCHP+iDDz7Q008/rUsvvTRQZ9GiRVq8eLGWLFmi3bt3y+v1auLEiTp58mRIny0id5SbNm2ajh07pl/+8peqqKjQ4MGD9cYbb6hv376SpIqKiqBr1tPS0vTGG2/ooYce0nPPPaeUlBQ988wzuuWWWyIRHgDAxfz/KOH2caGefPJJpaamavny5YFtV1xxReCfjTEqKCjQwoULNWXKFEnSypUrlZSUpNWrV+vee++94PeK2G1i58yZozlz5rS4b8WKFc22jR07Vnv37o1UOAAA2O7Ll1O3tN5r48aNmjRpkn7wgx+ouLhYl19+uebMmaN77rlH0rkbsPl8vqBLuz0ej8aOHasdO3aElNSd+ZgaAABaEe7K96YiSampqUGXV+fn5zd7v48//lhLly5Venq63nrrLd133336yU9+olWrVklSYKF4KJeCt8ZxD3QpbzgV9qNT/9ZwiU3RSBX1X7Oln2M13W3pR5Ia6jrZ0k9snS3dSJI61duzkjTGpn4kyQr3MU5NovTOU9HG+DlO6BiNRjY8pe3c/5eXlyshISGwvaWrsvx+v0aMGKG8vDxJ0rBhw7R//34tXbo0aH1ZqJeCt4SROgDAVfw2FUlKSEgIKi0l9eTkZA0cODBo29VXXx1YW+b1eiUppEvBW0NSBwAggq699lodPHgwaNuHH34YWDyelpYmr9cbdGl3XV2diouLW70UvDWOm34HAKAtfllqVGjT2i31caEeeughZWZmKi8vT1OnTtWf/vQnvfjii3rxxRclnZt2nzdvnvLy8pSenq709HTl5eWpW7dumj59ekhxkdQBAK7iN+dKuH1cqJEjR2r9+vXKycnRL3/5S6WlpamgoEB33HFHoM6CBQt09uxZzZkzR8ePH9eoUaO0efNmxcfHhxQXSR0AgAi76aabdNNNN7W637Is5ebmKjc3N6z3IakDAFyl0Ybp93DbRwpJHQDgKk5O6qx+BwDAIRipAwBcxW8s+U2Yq9/DbB8pJHUAgKsw/Q4AAKIeI3UAgKs0KkaNYY5pG22KxW4kdQCAqxgbzqkbzqkDAHDxcU4dAABEPUbqAABXaTQxajRhnlMP897xkUJSBwC4il+W/GFOVPsVnVmd6XcAABzCcSP17Wf7qGun8D5WRf3XbIpGer+6ty39fFqdYEs/kqTqLrZ00/mMfQtFOtXY86s3pt6+X89WvT0XrVgN9l38Yhpt6svvt6cfKfxnWAIdzMkL5RyX1AEAaIs959Sj88cs0+8AADgEI3UAgKucWygX5gNdmH4HAODi89twm1hWvwMAgIhipA4AcBUnL5QjqQMAXMWvGMfefIakDgBwlUZjqTHMp6yF2z5SOKcOAIBDMFIHALhKow2r3xuZfgcA4OLzmxj5w1wo54/ShXJMvwMA4BCM1AEArsL0OwAADuFX+KvXbXzOoa1sn37Pz8/XyJEjFR8fr169eun73/++Dh482GaboqIiWZbVrPz1r3+1OzwAABzL9pF6cXGx5s6dq5EjR6qhoUELFy5UVlaWPvjgA3Xv3r3NtgcPHlRCwufPDe/Zs6fd4QEAXM6em89E55I025P6m2++GfR6+fLl6tWrl9577z2NGTOmzba9evXSpZdeandIAAAE2HObWJck9S+rqqqSJPXo0eO8dYcNG6aamhoNHDhQP//5zzV+/PhW69bW1qq2tjbwurq6WpK05f8GqUttbFgxH6tpe0YhFJ9WJ5y/0gU4VWlfTJ7j9nwZY6vsWygSe9qeM1SdzjTY0o8kWXU29dXQaE8/ktRo05m8KL0cB0B4IvpTwxij+fPn67rrrtPgwYNbrZecnKwXX3xRa9eu1bp169S/f39NmDBB27Zta7VNfn6+EhMTAyU1NTUSHwEA4DBNz1MPt0SjiI7U77//fv35z3/W9u3b26zXv39/9e/fP/B69OjRKi8v11NPPdXqlH1OTo7mz58feF1dXU1iBwCcF9Pv7fDAAw9o48aN2rZtm3r37h1y+4yMDL3yyiut7vd4PPJ4POGECABwIXuuU3dJUjfG6IEHHtD69etVVFSktLS0dvVTUlKi5ORkm6MDAMC5bE/qc+fO1erVq/X6668rPj5ePp9PkpSYmKiuXbtKOjd1/sknn2jVqlWSpIKCAl1xxRUaNGiQ6urq9Morr2jt2rVau3at3eEBAFzObyz5w735TJQ+etX2pL506VJJ0rhx44K2L1++XLNnz5YkVVRU6MiRI4F9dXV1evjhh/XJJ5+oa9euGjRokDZt2qQbb7zR7vAAAC7nt2H63TXXqZsLuFRmxYoVQa8XLFigBQsW2B0KAACuwr3fAQCuYs+jV10yUgcAIJo1ylJjmNeZh9s+UqLzpwYAAAgZI3UAgKsw/Q4AgEM0Kvzpcxuf6GCr6PypAQAAQsZIHQDgKky/AwDgEDzQBQAAhzA2PDrVcEkbAADuk5ubK8uygorX6w3sN8YoNzdXKSkp6tq1q8aNG6f9+/e3671I6gAAV2mafg+3hGLQoEGqqKgIlH379gX2LVq0SIsXL9aSJUu0e/dueb1eTZw4USdPngz5szlu+r2kvLdiusWF1UdDXSebopFU3cWWbjzH7fv9FXfUnmkjzwm/Lf1IUpfqBlv66XS23pZ+JMmqqbOlH1NvX0xqtOdCGmNTP+c6s+l7YFc/knQBz6CAe12Mp7R17tw5aHTexBijgoICLVy4UFOmTJEkrVy5UklJSVq9erXuvffekN6HkToAAO1UXV0dVGpra1usd+jQIaWkpCgtLU233XabPv74Y0lSWVmZfD6fsrKyAnU9Ho/Gjh2rHTt2hBwPSR0A4CqN/3j0arhFklJTU5WYmBgo+fn5zd5v1KhRWrVqld566y299NJL8vl8yszM1LFjx+Tz+SRJSUlJQW2SkpIC+0LhuOl3AADaYuf0e3l5uRISEgLbPR5Ps7rZ2dmBfx4yZIhGjx6tfv36aeXKlcrIyJAkWVZwPMaYZtsuBCN1AADaKSEhIai0lNS/rHv37hoyZIgOHToUOM/+5VF5ZWVls9H7hSCpAwBcxa8YW0p71dbW6sCBA0pOTlZaWpq8Xq8KCwsD++vq6lRcXKzMzMyQ+2b6HQDgKo3GUmOY0++htH/44Yc1efJk9enTR5WVlfrVr36l6upqzZo1S5Zlad68ecrLy1N6errS09OVl5enbt26afr06SHHRVIHACCC/v73v+v222/X0aNH1bNnT2VkZGjXrl3q27evJGnBggU6e/as5syZo+PHj2vUqFHavHmz4uPjQ34vkjoAwFU6+jr1NWvWtLnfsizl5uYqNzc3rJgkkjoAwGWMDU9pMzzQBQCAi69RlhrDfCBLuO0jJTp/agAAgJAxUgcAuIrfhH7v9pb6iEYkdQCAq/htOKcebvtIic6oAABAyBipAwBcxS9L/jAXuoXbPlJI6gAAV+noO8p1JKbfAQBwCOeN1Mu7SXFxYXURW2dTLJI6n7Hn11xslX1LLT0n/Lb0E3e80ZZ+JKlLtT0HPeZUrS39SJJqbfoi1NXb048k09BgT0c2Lt010boMGGiFkxfKOS+pAwDQBr9suE1slJ5Tj86fGgAAIGSM1AEArmJsWP1uonSkTlIHALhKRz+lrSOR1AEAruLkhXK2R5WbmyvLsoKK1+tts01xcbGGDx+uuLg4XXnllXrhhRfsDgsAAMeLyEh90KBB2rJlS+B1p06dWq1bVlamG2+8Uffcc49eeeUV/fGPf9ScOXPUs2dP3XLLLZEIDwDgYky/h9pp587nHZ03eeGFF9SnTx8VFBRIkq6++mrt2bNHTz31FEkdAGA7J98mNiInBQ4dOqSUlBSlpaXptttu08cff9xq3Z07dyorKyto26RJk7Rnzx7V17d+047a2lpVV1cHFQAA3Mz2pD5q1CitWrVKb731ll566SX5fD5lZmbq2LFjLdb3+XxKSkoK2paUlKSGhgYdPXq01ffJz89XYmJioKSmptr6OQAAztQ0/R5uiUa2J/Xs7GzdcsstGjJkiL797W9r06ZNkqSVK1e22saygg+OMabF7V+Uk5OjqqqqQCkvL7chegCA0zk5qUf8krbu3btryJAhOnToUIv7vV6vfD5f0LbKykp17txZl112Wav9ejweeTweW2MFAOCrLOIX2tXW1urAgQNKTk5ucf/o0aNVWFgYtG3z5s0aMWKEunTpEunwAAAu4+SRuu1J/eGHH1ZxcbHKysr07rvv6tZbb1V1dbVmzZol6dy0+cyZMwP177vvPh0+fFjz58/XgQMH9PLLL2vZsmV6+OGH7Q4NAABHJ3Xbp9///ve/6/bbb9fRo0fVs2dPZWRkaNeuXerbt68kqaKiQkeOHAnUT0tL0xtvvKGHHnpIzz33nFJSUvTMM89wORsAACGyPamvWbOmzf0rVqxotm3s2LHau3ev3aEAANCMUfjXmRt7QrEd934HALgKd5QDAMAhSOpfId3+bqlTbHgHu1O9fRMrnWrs6Sv2tN+WfiSpS3WDTf3U2dKPJHU6VWtLP9aZGlv6kSRTa09MpsGe4y1JprHRpo7s+z7Z2heAsDguqQMA0BZG6gAAOISTk3p0PuUdAACEjJE6AMBVjLFkwhxph9s+UkjqAABX4XnqAAAg6jFSBwC4ipMXypHUAQCu4uRz6ky/AwDgEIzUAQCuwvQ7AAAO4eTpd5I6AMBVjA0j9WhN6pxTBwDAIRipAwBcxUgyYT5A075nedqLpA4AcBW/LFncUQ4AAEQzRuoAAFdh9ftXyCUVDercpSGsPmLq7TtbYldfnc6E95mC+jpbb0s/MadqbelHkqwzNbb0Y2rsi0l19hwn1dvUjyT57fk+GZv6sVW4JzmBC+Q3liyHXqfO9DsAAA5BUgcAuIox9pT2ys/Pl2VZmjdv3hdiMsrNzVVKSoq6du2qcePGaf/+/SH3TVIHALhK0zn1cEt77N69Wy+++KK+8Y1vBG1ftGiRFi9erCVLlmj37t3yer2aOHGiTp48GVL/JHUAADrAqVOndMcdd+ill17S1772tcB2Y4wKCgq0cOFCTZkyRYMHD9bKlSt15swZrV69OqT3IKkDAFzFzpF6dXV1UKmtbX2x7ty5c/Xd735X3/72t4O2l5WVyefzKSsrK7DN4/Fo7Nix2rFjR0ifjaQOAHCVpqe0hVskKTU1VYmJiYGSn5/f4nuuWbNGe/fubXG/z+eTJCUlJQVtT0pKCuy7UI67pA0AgLaEu9CtqQ9JKi8vV0JCQmC7x+NpVre8vFwPPvigNm/erLi4uFb7tKzg8/TGmGbbzoekDgBAOyUkJAQl9Za89957qqys1PDhwwPbGhsbtW3bNi1ZskQHDx6UdG7EnpycHKhTWVnZbPR+Pky/AwBc5dxIPdxz6hf+fhMmTNC+fftUWloaKCNGjNAdd9yh0tJSXXnllfJ6vSosLAy0qaurU3FxsTIzM0P6bIzUAQCu0tG3iY2Pj9fgwYODtnXv3l2XXXZZYPu8efOUl5en9PR0paenKy8vT926ddP06dNDioukDgDARbZgwQKdPXtWc+bM0fHjxzVq1Cht3rxZ8fHxIfVDUgcAuIpR+M9DD7d9UVFR0GvLspSbm6vc3Nyw+iWpAwBcxclPaWOhHAAADmF7Ur/iiitkWVazMnfu3BbrFxUVtVj/r3/9q92hAQDw+fx7uCUK2T79vnv3bjU2NgZe/+Uvf9HEiRP1gx/8oM12Bw8eDLrWr2fPnnaHBgCAZMP0u6J0+t32pP7lZPzEE0+oX79+Gjt2bJvtevXqpUsvvdTucAAACGLnHeWiTUTPqdfV1emVV17RnXfeed5b3Q0bNkzJycmaMGGCtm7dGsmwAABwpIiuft+wYYNOnDih2bNnt1onOTlZL774ooYPH67a2lr913/9lyZMmKCioiKNGTOm1Xa1tbVBT8Oprq6WJHX9rEadw/xUVqN9P8Gs+sbzV7qQfuoabOlHkqyaOns6qrWpH0mmjScbhaSu3p5+JJkGe465afTb0s+5vuz5PsnYF1PUDlmAVjh59XtEk/qyZcuUnZ2tlJSUVuv0799f/fv3D7wePXq0ysvL9dRTT7WZ1PPz8/X444/bGi8AwAWMFf458ShN6hGbfj98+LC2bNmiu+++O+S2GRkZOnToUJt1cnJyVFVVFSjl5eXtDRUAAEeI2Eh9+fLl6tWrl7773e+G3LakpCToSTUt8Xg8LT7iDgCAtjh5oVxEkrrf79fy5cs1a9Ysdf7SCe6cnBx98sknWrVqlSSpoKBAV1xxhQYNGhRYWLd27VqtXbs2EqEBANwuGu4TGyERSepbtmzRkSNHdOeddzbbV1FRoSNHjgRe19XV6eGHH9Ynn3yirl27atCgQdq0aZNuvPHGSIQGAIBjRSSpZ2VlybQyN7FixYqg1wsWLNCCBQsiEQYAAM2w+h0AACeJ0unzcPFAFwAAHIKROgDAVZh+BwDAKVj9DgCAU1j/KOH2EX04pw4AgEMwUgcAuAvT7wAAOISDkzrT7wAAOAQjdQCAuzj40askdQCAq/CUtq+QTv93Wp06NYTXiY1/Lauh0Z6O7OpHkqmvt6ejOpv6kWQawvybNbHrs0kyjX6b+rHvbydjT0wAnMlxSR0AgDY5eKEcSR0A4C4OPqfO6ncAAByCkToAwFUsc66E20c0IqkDANyFc+oAADgE59QBAEC0Y6QOAHAXpt8BAHAIByd1pt8BAHAIRuoAAHdx8EidpA4AcBdWvwMAgGjHSB0A4CrcUQ4AAKdw8Dl1pt8BAHAIkjoAAA7B9DsAwFUs2XBO3ZZI7Oe8pH7ytBRTf7GjCDCNjfZ01Oi3px9Jsikm09BgSz+SjcfJb9+JLttiMjb+7exiovSEINARuKQNAABEO+eN1AEAaIuDV7+T1AEA7uLgpM70OwAADkFSBwC4StMd5cItF2rp0qX6xje+oYSEBCUkJGj06NH6wx/+ENhvjFFubq5SUlLUtWtXjRs3Tvv372/XZws5qW/btk2TJ09WSkqKLMvShg0bgva3N7i1a9dq4MCB8ng8GjhwoNavXx9qaAAAnJ+xqVyg3r1764knntCePXu0Z88e3XDDDbr55psDuXHRokVavHixlixZot27d8vr9WrixIk6efJkyB8t5KR++vRpDR06VEuWLGlxf3uC27lzp6ZNm6YZM2bo/fff14wZMzR16lS9++67oYYHAEBUmTx5sm688UZdddVVuuqqq/Rv//ZvuuSSS7Rr1y4ZY1RQUKCFCxdqypQpGjx4sFauXKkzZ85o9erVIb9XyEk9Oztbv/rVrzRlypRm+9obXEFBgSZOnKicnBwNGDBAOTk5mjBhggoKCkINDwCAttk4Uq+urg4qtbW1bb51Y2Oj1qxZo9OnT2v06NEqKyuTz+dTVlZWoI7H49HYsWO1Y8eOkD+arefU2xvczp07g9pI0qRJk9psU1tb2+xgAgBwPnaeU09NTVViYmKg5Ofnt/ie+/bt0yWXXCKPx6P77rtP69ev18CBA+Xz+SRJSUlJQfWTkpIC+0Jh6yVtbQV3+PDhNtuF+oHy8/P1+OOPhxEtAADhKS8vV0JCQuC1x+NpsV7//v1VWlqqEydOaO3atZo1a5aKi4sD+y0r+A51xphm2y5ERFa/tye4UNvk5OSoqqoqUMrLy9sfMADAPZpuExtukQIr2ptKa0k9NjZW//RP/6QRI0YoPz9fQ4cO1b//+7/L6/VKUrNBbGVlZbPB7oWwNam3Nziv1xtyG4/H0+xgAgBwXh28+r3FEIxRbW2t0tLS5PV6VVhYGNhXV1en4uJiZWZmhtyvrUm9vcGNHj06qI0kbd68uV0fCACAtnT0deqPPvqo3nnnHf3tb3/Tvn37tHDhQhUVFemOO+6QZVmaN2+e8vLytH79ev3lL3/R7Nmz1a1bN02fPj3kzxbyOfVTp07po48+CrwuKytTaWmpevTooT59+gSCS09PV3p6uvLy8poFN3PmTF1++eWBBQUPPvigxowZoyeffFI333yzXn/9dW3ZskXbt28P+QMBABBNPvvsM82YMUMVFRVKTEzUN77xDb355puaOHGiJGnBggU6e/as5syZo+PHj2vUqFHavHmz4uPjQ36vkJP6nj17NH78+MDr+fPnS5JmzZqlFStWXFBwR44cUUzM55MEmZmZWrNmjX7+85/rscceU79+/fTaa69p1KhRIX8gAADa1MH3fl+2bFmb+y3LUm5urnJzc8OLSZJljDMerFxdXa3ExERN6HW3OsfEXuxwPsfz1C+sL56n3rGc8a89HKTB1KtIr6uqqipia6Sa8sSVj+WpU1xcWH011tTo4399NKLxtgf3fgcAwCEc9+hVc/asjGXTCMsOfptGaTaOrKJxVGzXaNZEYUy2YoQNhM/Bj151XFIHAKBNDk7qTL8DAOAQjNQBAK4S6nXmrfURjRipAwDgECR1AAAcgul3AIC7OHihHEkdAOAqTj6nTlIHALhPlCblcHFOHQAAh2CkDgBwF86pAwDgDE4+p870OwAADsFIHQDgLky/AwDgDEy/AwCAqMdIHQDgLky/AwDgEA5O6ky/AwDgEI4bqZu6ehnrYkfxBf4o/Dln/DZ1Y+NnsykmW5ko/NsBCJuTF8o5LqkDANAmB0+/k9QBAO7i4KTOOXUAAByCkToAwFU4pw4AgFMw/Q4AAKIdI3UAgKsw/Q4AgFMw/Q4AAKIdI3UAgLs4eKROUgcAuIr1jxJuH9GI6XcAAByCkToAwF2YfgcAwBmcfElbyNPv27Zt0+TJk5WSkiLLsrRhw4bAvvr6ev3sZz/TkCFD1L17d6WkpGjmzJn69NNP2+xzxYoVsiyrWampqQn5AwEA0CZjU4lCISf106dPa+jQoVqyZEmzfWfOnNHevXv12GOPae/evVq3bp0+/PBDfe973ztvvwkJCaqoqAgqcXFxoYYHAIBrhTz9np2drezs7Bb3JSYmqrCwMGjbs88+q2uuuUZHjhxRnz59Wu3Xsix5vd5QwwEAIHRROtIOV8TPqVdVVcmyLF166aVt1jt16pT69u2rxsZGffOb39S//uu/atiwYa3Wr62tVW1tbeB1dXX1uX9obJSs8Bb1G79D/9pNjP9iRxA5xuF/OwBh45x6O9XU1OiRRx7R9OnTlZCQ0Gq9AQMGaMWKFdq4caNeffVVxcXF6dprr9WhQ4dabZOfn6/ExMRASU1NjcRHAADgKyNiSb2+vl633Xab/H6/nn/++TbrZmRk6Ic//KGGDh2q66+/Xr/73e901VVX6dlnn221TU5OjqqqqgKlvLzc7o8AAHAiBy+Ui8j0e319vaZOnaqysjK9/fbbbY7SWxITE6ORI0e2OVL3eDzyeDzhhgoAcBmm30PQlNAPHTqkLVu26LLLLgu5D2OMSktLlZycbHd4AAA4Vsgj9VOnTumjjz4KvC4rK1Npaal69OihlJQU3Xrrrdq7d6/+53/+R42NjfL5fJKkHj16KDY2VpI0c+ZMXX755crPz5ckPf7448rIyFB6erqqq6v1zDPPqLS0VM8995wdnxEAgM9xR7nP7dmzR+PHjw+8nj9/viRp1qxZys3N1caNGyVJ3/zmN4Pabd26VePGjZMkHTlyRDExn08SnDhxQj/+8Y/l8/mUmJioYcOGadu2bbrmmmtCDQ8AgDY5efo95KQ+btw4mTYuG2prX5OioqKg17/5zW/0m9/8JtRQAADAF3DvdwCAuzh4+p1HrwIA3KWDL2nLz8/XyJEjFR8fr169eun73/++Dh48GBySMcrNzVVKSoq6du2qcePGaf/+/SF/NJI6AMBVms6ph1suVHFxsebOnatdu3apsLBQDQ0NysrK0unTpwN1Fi1apMWLF2vJkiXavXu3vF6vJk6cqJMnT4b02Zh+BwAggt58882g18uXL1evXr303nvvacyYMTLGqKCgQAsXLtSUKVMkSStXrlRSUpJWr16te++994Lfi5E6AMBdbJx+r66uDipffCZJa6qqqiSdu9RbOndpuM/nU1ZWVqCOx+PR2LFjtWPHjpA+GkkdAOAqljG2FElKTU0Neg5J0/1XWmOM0fz583Xddddp8ODBkhS4n0tSUlJQ3aSkpMC+C8X0OwAA7VReXh50K/Tz3b78/vvv15///Gdt37692T7LsoJeG2OabTsfkjoAwF1svKQtISHhgp9v8sADD2jjxo3atm2bevfuHdju9XolnRuxf/H26JWVlc1G7+fD9DsAwFU6evW7MUb333+/1q1bp7fffltpaWlB+9PS0uT1elVYWBjYVldXp+LiYmVmZob02RipAwAQQXPnztXq1av1+uuvKz4+PnCePDExUV27dpVlWZo3b57y8vKUnp6u9PR05eXlqVu3bpo+fXpI70VSBwC4SwffUW7p0qWSFHj+SZPly5dr9uzZkqQFCxbo7NmzmjNnjo4fP65Ro0Zp8+bNio+PDyksxyV109goY3FWoU0XcH9+AHCqjn6gy4U8E8WyLOXm5io3N7f9QYlz6gAAOIbjRuoAALTJwQ90IakDAFyF56kDAOAUDh6pc04dAACHYKQOAHCdaJ0+DxdJHQDgLsaEf2lvlF4azPQ7AAAOwUgdAOAqrH4HAMApWP0OAACiHSN1AICrWP5zJdw+ohFJHQDgLky/AwCAaMdIHQDgKqx+BwDAKRx88xmSOgDAVRipf5UYO1ZAAADw1eO8pA4AQFscvPqdpA4AcBUnT79zSRsAAA7BSB0A4C6sfgcAwBmYfv+Cbdu2afLkyUpJSZFlWdqwYUPQ/tmzZ8uyrKCSkZFx3n7Xrl2rgQMHyuPxaODAgVq/fn2ooQEA4GohJ/XTp09r6NChWrJkSat1vvOd76iioiJQ3njjjTb73Llzp6ZNm6YZM2bo/fff14wZMzR16lS9++67oYYHAEDbjE0lCoU8/Z6dna3s7Ow263g8Hnm93gvus6CgQBMnTlROTo4kKScnR8XFxSooKNCrr74aaogAALSK6fcQFRUVqVevXrrqqqt0zz33qLKyss36O3fuVFZWVtC2SZMmaceOHa22qa2tVXV1dVABAMDNbE/q2dnZ+u1vf6u3335bTz/9tHbv3q0bbrhBtbW1rbbx+XxKSkoK2paUlCSfz9dqm/z8fCUmJgZKamqqbZ8BAOBgfmNPiUK2r36fNm1a4J8HDx6sESNGqG/fvtq0aZOmTJnSajvLsoJeG2OabfuinJwczZ8/P/C6urqaxA4AOD/uKNd+ycnJ6tu3rw4dOtRqHa/X22xUXllZ2Wz0/kUej0cej8e2OAEA7mDJhnPqtkRiv4jfUe7YsWMqLy9XcnJyq3VGjx6twsLCoG2bN29WZmZmpMMDAMAxQh6pnzp1Sh999FHgdVlZmUpLS9WjRw/16NFDubm5uuWWW5ScnKy//e1vevTRR/X1r39d//zP/xxoM3PmTF1++eXKz8+XJD344IMaM2aMnnzySd188816/fXXtWXLFm3fvt2GjwgAwBdwR7nP7dmzR+PHjw+8bjqvPWvWLC1dulT79u3TqlWrdOLECSUnJ2v8+PF67bXXFB8fH2hz5MgRxcR8PkmQmZmpNWvW6Oc//7kee+wx9evXT6+99ppGjRoVzmcDAKAZJ1/SFnJSHzdunEwbv1Deeuut8/ZRVFTUbNutt96qW2+9NdRwAADAP3DvdwCAu7D6HQAAZ7CMkRXmOfFw20cKz1MHAMAhGKkDANzF/48Sbh9RiKQOAHAVpt8BAEDUY6QOAHAXVr8DAOAQ3FEOAABncPId5TinDgCAQzBSBwC4C9PvAAA4g+U/V8LtIxox/Q4AgEOQ1AEA7tI0/R5uCcG2bds0efJkpaSkyLIsbdiw4UshGeXm5iolJUVdu3bVuHHjtH///pA/GkkdAOAuxqYSgtOnT2vo0KFasmRJi/sXLVqkxYsXa8mSJdq9e7e8Xq8mTpyokydPhvQ+nFMHACDCsrOzlZ2d3eI+Y4wKCgq0cOFCTZkyRZK0cuVKJSUlafXq1br33nsv+H0YqQMAXKXp3u/hFkmqrq4OKrW1tSHHU1ZWJp/Pp6ysrMA2j8ejsWPHaseOHSH1RVIHALiLjefUU1NTlZiYGCj5+fkhh+Pz+SRJSUlJQduTkpIC+y4U0+8AALRTeXm5EhISAq89Hk+7+7IsK+i1MabZtvMhqQMA3MUo/Oeh/2OhXEJCQlBSbw+v1yvp3Ig9OTk5sL2ysrLZ6P18mH4HALiKnefU7ZCWliav16vCwsLAtrq6OhUXFyszMzOkvhipAwDcxciG28SGVv3UqVP66KOPAq/LyspUWlqqHj16qE+fPpo3b57y8vKUnp6u9PR05eXlqVu3bpo+fXpI70NSBwAgwvbs2aPx48cHXs+fP1+SNGvWLK1YsUILFizQ2bNnNWfOHB0/flyjRo3S5s2bFR8fH9L7kNQBAO5yER7oMm7cOJk22liWpdzcXOXm5oYVFkkdAOAufkmhLSpvuY8oxEI5AAAcgpE6AMBV7Fi9bufqdzuR1AEA7nIRzql3FKbfAQBwCEbqAAB3cfBInaQOAHAXByd1pt8BAHAIRuoAAHdx8HXqJHUAgKtwSRsAAE7BOfXPbdu2TZMnT1ZKSoosy9KGDRuC9luW1WL59a9/3WqfK1asaLFNTU1NyB8IAAC3Cnmkfvr0aQ0dOlQ/+tGPdMsttzTbX1FREfT6D3/4g+66664W635RQkKCDh48GLQtLi4u1PAAAGib30hWmCNtf3SO1ENO6tnZ2crOzm51v9frDXr9+uuva/z48bryyivb7NeyrGZtAQCwHdPv7fPZZ59p06ZNuuuuu85b99SpU+rbt6969+6tm266SSUlJW3Wr62tVXV1dVABAMDNIprUV65cqfj4eE2ZMqXNegMGDNCKFSu0ceNGvfrqq4qLi9O1116rQ4cOtdomPz9fiYmJgZKammp3+AAARzKfj9bbW+TCkfrLL7+sO+6447znxjMyMvTDH/5QQ4cO1fXXX6/f/e53uuqqq/Tss8+22iYnJ0dVVVWBUl5ebnf4AAAnCjeh2zF9HyERu6TtnXfe0cGDB/Xaa6+F3DYmJkYjR45sc6Tu8Xjk8XjCCREAAEeJ2Eh92bJlGj58uIYOHRpyW2OMSktLlZycHIHIAACu5jf2lCgU8kj91KlT+uijjwKvy8rKVFpaqh49eqhPnz6SpOrqav3+97/X008/3WIfM2fO1OWXX678/HxJ0uOPP66MjAylp6erurpazzzzjEpLS/Xcc8+15zMBANA64z9Xwu0jCoWc1Pfs2aPx48cHXs+fP1+SNGvWLK1YsUKStGbNGhljdPvtt7fYx5EjRxQT8/kkwYkTJ/TjH/9YPp9PiYmJGjZsmLZt26Zrrrkm1PAAAHAty5goPdsfourqaiUmJmqcblZnq8vFDgcAEIIGU68iva6qqiolJCRE5D2a8sS3U/9FnWPCW5PV4K/VlvKlEY23Pbj3OwDAXfw2XJLmlHPqAAB8pXFHOQAAEO0YqQMA3MXIhpG6LZHYjqQOAHAXpt8BAEC0Y6QOAHAXv19SmDeP8Tvk5jMAAHylMf0OAACiHSN1AIC7OHikTlIHALiLg+8ox/Q7AAAOwUgdAOAqxvhlwnx0arjtI4WkDgBwF2PCnz7nnDoAAFHA2HBOPUqTOufUAQBwCEbqAAB38fslK8xz4pxTBwAgCjD9DgAAoh0jdQCAqxi/XybM6XcuaQMAIBow/Q4AAKIdI3UAgLv4jWQ5c6ROUgcAuIsxksK9pC06kzrT7wAAOAQjdQCAqxi/kQlz+t1E6UidpA4AcBfjV/jT71zSBgDARefkkTrn1AEAcAjHjNSbfjU1qD7sewoAADpWg+oldcwIuMHUhj193hRvtHFMUj958qQkabveuMiRAADa6+TJk0pMTIxI37GxsfJ6vdrusydPeL1excbG2tKXXSwTrScGQuT3+/Xpp58qPj5elmW1WKe6ulqpqakqLy9XQkJCB0fYfsTd8b6qsRN3xyJu+xhjdPLkSaWkpCgmJnJnhmtqalRXV2dLX7GxsYqLi7OlL7s4ZqQeExOj3r17X1DdhISEqPkih4K4O95XNXbi7ljEbY9IjdC/KC4uLuoSsZ1YKAcAgEOQ1AEAcAhXJXWPx6Nf/OIX8ng8FzuUkBB3x/uqxk7cHYu4EW0cs1AOAAC3c9VIHQAAJyOpAwDgECR1AAAcgqQOAIBDOC6pP//880pLS1NcXJyGDx+ud955p836xcXFGj58uOLi4nTllVfqhRde6KBIz8nPz9fIkSMVHx+vXr166fvf/74OHjzYZpuioiJZltWs/PWvf+2gqKXc3Nxm7+/1ettsc7GPdZMrrriixeM3d+7cFutfrOO9bds2TZ48WSkpKbIsSxs2bAjab4xRbm6uUlJS1LVrV40bN0779+8/b79r167VwIED5fF4NHDgQK1fv77D4q6vr9fPfvYzDRkyRN27d1dKSopmzpypTz/9tM0+V6xY0eLfoKampkPilqTZs2c3e/+MjIzz9nsxj7ekFo+bZVn69a9/3WqfHXG8ERmOSuqvvfaa5s2bp4ULF6qkpETXX3+9srOzdeTIkRbrl5WV6cYbb9T111+vkpISPfroo/rJT36itWvXdljMxcXFmjt3rnbt2qXCwkI1NDQoKytLp0+fPm/bgwcPqqKiIlDS09M7IOLPDRo0KOj99+3b12rdaDjWTXbv3h0Ud2FhoSTpBz/4QZvtOvp4nz59WkOHDtWSJUta3L9o0SItXrxYS5Ys0e7du+X1ejVx4sTAcxBasnPnTk2bNk0zZszQ+++/rxkzZmjq1Kl69913OyTuM2fOaO/evXrssce0d+9erVu3Th9++KG+973vnbffhISEoONfUVFh653Bzne8Jek73/lO0Pu/8Ubb9xC/2MdbUrNj9vLLL8uyLN1yyy1t9hvp440IMQ5yzTXXmPvuuy9o24ABA8wjjzzSYv0FCxaYAQMGBG279957TUZGRsRiPJ/KykojyRQXF7daZ+vWrUaSOX78eMcF9iW/+MUvzNChQy+4fjQe6yYPPvig6devn/H7/S3uj4bjLcmsX78+8Nrv9xuv12ueeOKJwLaamhqTmJhoXnjhhVb7mTp1qvnOd74TtG3SpEnmtttusz1mY5rH3ZI//elPRpI5fPhwq3WWL19uEhMT7Q2uDS3FPWvWLHPzzTeH1E80Hu+bb77Z3HDDDW3W6ejjDfs4ZqReV1en9957T1lZWUHbs7KytGPHjhbb7Ny5s1n9SZMmac+ePaqvvziP1auqqpIk9ejR47x1hw0bpuTkZE2YMEFbt26NdGjNHDp0SCkpKUpLS9Ntt92mjz/+uNW60XispXPfm1deeUV33nlnqw8CanKxj/cXlZWVyefzBR1Tj8ejsWPHtvp9l1r/O7TVJtKqqqpkWZYuvfTSNuudOnVKffv2Ve/evXXTTTeppKSkYwL8gqKiIvXq1UtXXXWV7rnnHlVWVrZZP9qO92effaZNmzbprrvuOm/daDjeCJ1jkvrRo0fV2NiopKSkoO1JSUny+XwttvH5fC3Wb2ho0NGjRyMWa2uMMZo/f76uu+46DR48uNV6ycnJevHFF7V27VqtW7dO/fv314QJE7Rt27YOi3XUqFFatWqV3nrrLb300kvy+XzKzMzUsWPHWqwfbce6yYYNG3TixAnNnj271TrRcLy/rOk7Hcr3valdqG0iqaamRo888oimT5/e5oNFBgwYoBUrVmjjxo169dVXFRcXp2uvvVaHDh3qsFizs7P129/+Vm+//baefvpp7d69WzfccINqa2tbbRNtx3vlypWKj4/XlClT2qwXDccb7eOYp7Q1+fJoyxjT5gispfotbe8I999/v/785z9r+/btbdbr37+/+vfvH3g9evRolZeX66mnntKYMWMiHaakc/+BazJkyBCNHj1a/fr108qVKzV//vwW20TTsW6ybNkyZWdnKyUlpdU60XC8WxPq9729bSKhvr5et912m/x+v55//vk262ZkZAQtSrv22mv1rW99S88++6yeeeaZSIcqSZo2bVrgnwcPHqwRI0aob9++2rRpU5tJMlqOtyS9/PLLuuOOO857bjwajjfaxzEj9a9//evq1KlTs1/AlZWVzX4pN/F6vS3W79y5sy677LKIxdqSBx54QBs3btTWrVsv+BGyX5SRkXFRf0V3795dQ4YMaTWGaDrWTQ4fPqwtW7bo7rvvDrntxT7eTVcahPJ9b2oXaptIqK+v19SpU1VWVqbCwsKQH/8ZExOjkSNHXtS/QXJysvr27dtmDNFyvCXpnXfe0cGDB9v1fY+G440L45ikHhsbq+HDhwdWMjcpLCxUZmZmi21Gjx7drP7mzZs1YsQIdenSJWKxfpExRvfff7/WrVunt99+W2lpae3qp6SkRMnJyTZHd+Fqa2t14MCBVmOIhmP9ZcuXL1evXr303e9+N+S2F/t4p6Wlyev1Bh3Turo6FRcXt/p9l1r/O7TVxm5NCf3QoUPasmVLu37UGWNUWlp6Uf8Gx44dU3l5eZsxRMPxbrJs2TINHz5cQ4cODbltNBxvXKCLtUIvEtasWWO6dOlili1bZj744AMzb9480717d/O3v/3NGGPMI488YmbMmBGo//HHH5tu3bqZhx56yHzwwQdm2bJlpkuXLua///u/Oyzmf/mXfzGJiYmmqKjIVFRUBMqZM2cCdb4c929+8xuzfv168+GHH5q//OUv5pFHHjGSzNq1azss7p/+9KemqKjIfPzxx2bXrl3mpptuMvHx8VF9rL+osbHR9OnTx/zsZz9rti9ajvfJkydNSUmJKSkpMZLM4sWLTUlJSWCV+BNPPGESExPNunXrzL59+8ztt99ukpOTTXV1daCPGTNmBF398cc//tF06tTJPPHEE+bAgQPmiSeeMJ07dza7du3qkLjr6+vN9773PdO7d29TWloa9J2vra1tNe7c3Fzz5ptvmv/93/81JSUl5kc/+pHp3Lmzeffddzsk7pMnT5qf/vSnZseOHaasrMxs3brVjB492lx++eVRfbybVFVVmW7dupmlS5e22MfFON6IDEcldWOMee6550zfvn1NbGys+da3vhV0adisWbPM2LFjg+oXFRWZYcOGmdjYWHPFFVe0+qWPFEktluXLl7ca95NPPmn69etn4uLizNe+9jVz3XXXmU2bNnVo3NOmTTPJycmmS5cuJiUlxUyZMsXs37+/1ZiNufjH+oveeustI8kcPHiw2b5oOd5Nl9J9ucyaNcsYc+6ytl/84hfG6/Uaj8djxowZY/bt2xfUx9ixYwP1m/z+9783/fv3N126dDEDBgyw/cdJW3GXlZW1+p3funVrq3HPmzfP9OnTx8TGxpqePXuarKwss2PHjg6L+8yZMyYrK8v07NnTdOnSxfTp08fMmjXLHDlyJKiPaDveTf7jP/7DdO3a1Zw4caLFPi7G8UZk8OhVAAAcwjHn1AEAcDuSOgAADkFSBwDAIUjqAAA4BEkdAACHIKkDAOAQJHUAAByCpA4AgEOQ1AEAcAiSOgAADkFSBwDAIUjqAAA4xP8D7z4vfwrZU84AAAAASUVORK5CYII=",
|
|
"text/plain": [
|
|
"<Figure size 640x480 with 2 Axes>"
|
|
]
|
|
},
|
|
"metadata": {},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"plt.imshow(testConcentrations)\n",
|
|
"plt.colorbar()\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": []
|
|
}
|
|
],
|
|
"metadata": {
|
|
"kernelspec": {
|
|
"display_name": "forschungsprojekt",
|
|
"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
|
|
}
|