mirror of
https://git.gfz-potsdam.de/naaice/tug.git
synced 2025-12-13 09:28:23 +01:00
295 lines
13 KiB
Plaintext
295 lines
13 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 1,
|
|
"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": 4,
|
|
"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": 5,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def sub(delta_t, delta_x):\n",
|
|
" return delta_t / (2 * delta_x**2)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 15,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# CONSTANT\n",
|
|
"\n",
|
|
"# creates the coeffiecient matrix for a given row\n",
|
|
"def create_coeff_matrix(alpha_x, row_index, s_x, transpose=False):\n",
|
|
" \n",
|
|
"\n",
|
|
" r = row_index\n",
|
|
" if (transpose):\n",
|
|
" alpha_x = alpha_x.transpose()\n",
|
|
"\n",
|
|
" # this is always a square matrix -> column by column, so \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]) + 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]) + 2 * alpha_x[n,n]) # alpha_x[n,n]?\n",
|
|
"\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_y, alpha_x, col_index, row_index, s_y, s_x, bc_left, bc_right, transpose=False):\n",
|
|
"\n",
|
|
" # \n",
|
|
" sv = np.zeros(alpha_y.shape[1])\n",
|
|
"\n",
|
|
" # inner\n",
|
|
" # TODO: iterate over columns of alpha_y alpha_y[i_sv, row_index+ -1 .. 0 .. 1]\n",
|
|
" # TODO: error with current indices: i_sv starts from 0 --> i_alpha = -1 and i_alpha-1 = -2\n",
|
|
" for i_sv in range(0, sv.shape[0]):\n",
|
|
" print(i_sv)\n",
|
|
" i_alpha = i_sv-1\n",
|
|
" print(i_alpha)\n",
|
|
" print(concentrations[i_alpha-1,col_index])\n",
|
|
" sv[i_sv] = s_y * alpha_interblock(alpha_y[i_alpha+1,col_index],alpha_y[i_alpha,col_index])*concentrations[i_alpha+1,col_index] \\\n",
|
|
" + (1 - s_y *(alpha_interblock(alpha_y[i_alpha+1,col_index], alpha_y[i_alpha,col_index])) \\\n",
|
|
" + alpha_interblock(alpha_y[i_alpha,col_index], alpha_y[i_alpha-1,col_index]))*concentrations[i_alpha,col_index] \\\n",
|
|
" + s_y * alpha_interblock(alpha_y[i_alpha,col_index], alpha_y[i_alpha-1,col_index])*concentrations[i_alpha-1,col_index]\n",
|
|
"\n",
|
|
" # first\n",
|
|
" sv[0] += 2 * s_x * alpha_x[row_index,0] * bc_left # alpha_x[row_index,0]?\n",
|
|
"\n",
|
|
" # last\n",
|
|
" n = sv.shape[0]-1\n",
|
|
" sv[n] += 2 * s_x * alpha_x[row_index,n] * bc_right\n",
|
|
"\n",
|
|
" return sv\n",
|
|
"\n",
|
|
"\n",
|
|
"def solve_row(A, b):\n",
|
|
" return np.linalg.solve(A,b)\n",
|
|
"\n",
|
|
"def calc_sweep(concentrations, alpha_x, alpha_y, s_y, s_x, bc_left, bc_right, transpose=False):\n",
|
|
" for row in range(concentrations.shape[0]):\n",
|
|
" A = create_coeff_matrix(alpha_x, row, s_x, transpose)\n",
|
|
" b = create_solution_vector(concentrations, alpha_y, alpha_x, )\n",
|
|
"\n",
|
|
"def calc_round():\n",
|
|
" calc_sweep()\n",
|
|
" calc_sweep()\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 17,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"[[0.56726545 0.41188268 0.11335637 0.29221913 0.5365537 ]\n",
|
|
" [0.8120275 0.96249867 0.18743648 0.30176834 0.45374814]\n",
|
|
" [0.17014599 0.26796332 0.4616269 0.23477139 0.71348435]\n",
|
|
" [0.00654629 0.61879048 0.22742861 0.25794781 0.91561469]\n",
|
|
" [0.04155862 0.46501565 0.8463777 0.81575064 0.1985415 ]]\n",
|
|
"Test A:\n",
|
|
"[[ 7.77257437 -1.49738434 0. 0. 0. ]\n",
|
|
" [-2.78866756 5.2860519 -1.49738434 0. 0. ]\n",
|
|
" [ 0. -1.49738434 2.94508723 -0.44770289 0. ]\n",
|
|
" [ 0. 0. -0.44770289 1.90753379 -0.4598309 ]\n",
|
|
" [ 0. 0. 0. 1.69958064 4.50889084]]\n",
|
|
"0\n",
|
|
"-1\n",
|
|
"0.618790482597486\n",
|
|
"1\n",
|
|
"0\n",
|
|
"0.4650156531117613\n",
|
|
"2\n",
|
|
"1\n",
|
|
"0.41188268496900426\n",
|
|
"3\n",
|
|
"2\n",
|
|
"0.9624986706480608\n",
|
|
"4\n",
|
|
"3\n",
|
|
"0.26796331918576755\n",
|
|
"\n",
|
|
"TestB:\n",
|
|
"[6.33240747 1.20531698 0.62203484 1.97054078 6.92562211]\n",
|
|
"\n",
|
|
"Result\n",
|
|
"[1.01442248 1.03665212 0.96541206 1.49411867 0.97279954]\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"row = 5\n",
|
|
"col = 5\n",
|
|
"\n",
|
|
"testConcentrations = np.random.random((row,col))\n",
|
|
"print(testConcentrations)\n",
|
|
"testAlpha = np.random.random((row,col))\n",
|
|
"testA = create_coeff_matrix(testAlpha, 1, 3)\n",
|
|
"print(\"Test A:\")\n",
|
|
"print(testA)\n",
|
|
"# print(test.shape[0])\n",
|
|
"testb = create_solution_vector(testConcentrations,testAlpha,testAlpha,1,1,3,3,1,1)\n",
|
|
"print(\"\\nTestB:\")\n",
|
|
"print(testb)\n",
|
|
"\n",
|
|
"# needed to solve: TODO calc boundary in functions\n",
|
|
"# testA[0,0] = 1\n",
|
|
"# testA[4,4] = 1\n",
|
|
"# testb[0] = 1\n",
|
|
"# test[4] = 1\n",
|
|
"\n",
|
|
"result = solve_row(testA, testb)\n",
|
|
"print(\"\\nResult\")\n",
|
|
"print(result)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 8,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"3\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"test2 = np.random.random_sample(3)\n",
|
|
"print(test2.shape[0])"
|
|
]
|
|
},
|
|
{
|
|
"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
|
|
}
|