mirror of
https://git.gfz-potsdam.de/naaice/tug.git
synced 2025-12-13 09:28:23 +01:00
267 lines
20 KiB
Plaintext
267 lines
20 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 2,
|
|
"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": "iVBORw0KGgoAAAANSUhEUgAAAa4AAAGdCAYAAABKG5eZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAkE0lEQVR4nO3df3RU9Z3/8dfFhAlyklGUJDMQQuQglB+HQkACKj+kBENFWKmg7oGwtrZsqRVTTiG2HnH/aLCtLgdBWbv8kLWLnG4IZDdsJRzzQ0tgQRLXWsS4RpKVpBw4JQO4DAE+3z/8ZuqYmcGBmSSf+Hycc8/x3vv+fPKeTyKv3Jk7E8cYYwQAgCV6dXUDAABEg+ACAFiF4AIAWIXgAgBYheACAFiF4AIAWIXgAgBYheACAFgloasbiJUrV67oxIkTSk5OluM4Xd0OACAKxhidPXtWXq9XvXpFvqbqMcF14sQJZWRkdHUbAIDr0NTUpIEDB0as6THBlZycLEm6S7OVoMQu7gYAEI1LatPb2hP4tzySHhNc7U8PJihRCQ7BBQBW+f+fmvtVXurh5gwAgFUILgCAVeIWXC+99JKysrKUlJSk7OxsvfXWWxHrq6qqlJ2draSkJN12223auHFjvFoDAFgsLsG1Y8cOLV++XD/72c9UW1uru+++W3l5eWpsbAxZ39DQoNmzZ+vuu+9WbW2tnnrqKf34xz9WcXFxPNoDAFjMiccfkpw4caLGjRunl19+OXDsG9/4hubNm6eioqIO9StXrlRpaamOHj0aOLZ06VK9++67qqmp+Upf0+fzye12a5rmcnMGAFjmkmlTpXartbVVKSkpEWtjfsV18eJFvfPOO8rNzQ06npubq/3794ccU1NT06F+1qxZOnz4sNra2kKO8fv98vl8QRsAoOeLeXCdOnVKly9fVlpaWtDxtLQ0tbS0hBzT0tISsv7SpUs6depUyDFFRUVyu92BjTcfA8DXQ9xuzvjyvfjGmIj354eqD3W8XWFhoVpbWwNbU1PTdXYMALBBzN+AfOutt+qGG27ocHV18uTJDldV7dLT00PWJyQk6JZbbgk5xuVyyeVyxaZpAIA1Yn7F1bt3b2VnZ6u8vDzoeHl5uSZPnhxyzKRJkzrU7927V+PHj1diIjdaAAD+Ki5PFRYUFOif//mftXnzZh09elRPPvmkGhsbtXTpUkmfP823ePHiQP3SpUt1/PhxFRQU6OjRo9q8ebM2bdqkFStWxKM9AIDF4vJZhQsXLtTp06f1D//wD2pubtaoUaO0Z88eZWZmSpKam5uD3tOVlZWlPXv26Mknn9SGDRvk9Xq1bt06zZ8/Px7tAQAsFpf3cXUF3scFAPbq0vdxAQAQTwQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqMQ+uoqIiTZgwQcnJyUpNTdW8efN07NixiGMqKyvlOE6H7YMPPoh1ewAAy8U8uKqqqrRs2TIdOHBA5eXlunTpknJzc3X+/Pmrjj127Jiam5sD29ChQ2PdHgDAcgmxnvD3v/990P6WLVuUmpqqd955R1OmTIk4NjU1VTfddFOsWwIA9CBxf42rtbVVktSvX7+r1o4dO1Yej0czZsxQRUVFxFq/3y+fzxe0AQB6vrgGlzFGBQUFuuuuuzRq1KiwdR6PR6+88oqKi4u1c+dODRs2TDNmzFB1dXXYMUVFRXK73YEtIyMjHg8BANDNOMYYE6/Jly1bprKyMr399tsaOHBgVGPnzJkjx3FUWloa8rzf75ff7w/s+3w+ZWRkaJrmKsFJvK6+AQCd65JpU6V2q7W1VSkpKRFr43bF9fjjj6u0tFQVFRVRh5Yk5eTkqL6+Pux5l8ullJSUoA0A0PPF/OYMY4wef/xxlZSUqLKyUllZWdc0T21trTweT4y7AwDYLubBtWzZMv3rv/6rdu/ereTkZLW0tEiS3G63+vTpI0kqLCzUp59+qm3btkmS1q5dq8GDB2vkyJG6ePGiXnvtNRUXF6u4uDjW7QEALBfz4Hr55ZclSdOmTQs6vmXLFi1ZskSS1NzcrMbGxsC5ixcvasWKFfr000/Vp08fjRw5UmVlZZo9e3as2wMAWC6uN2d0Jp/PJ7fbzc0ZAGChbnFzBgAA8UBwAQCsQnABAKxCcAEArEJwAQCsQnABAKxCcAEArEJwAQCsQnABAKxCcAEArEJwAQCsQnABAKxCcAEArEJwAQCsQnABAKxCcAEArEJwAQCsQnABAKxCcAEArEJwAQCsQnABAKxCcAEArEJwAQCsQnABAKxCcAEArEJwAQCsQnABAKxCcAEArEJwAQCsQnABAKxCcAEArEJwAQCsQnABAKxCcAEArEJwAQCsEvPgWr16tRzHCdrS09MjjqmqqlJ2draSkpJ02223aePGjbFuCwDQQyTEY9KRI0dq3759gf0bbrghbG1DQ4Nmz56txx57TK+99pr+8Ic/6Ic//KH69++v+fPnx6M9AIDF4hJcCQkJV73Kardx40YNGjRIa9eulSR94xvf0OHDh/XrX/+a4AIAdBCX17jq6+vl9XqVlZWlhx56SB9//HHY2pqaGuXm5gYdmzVrlg4fPqy2traw4/x+v3w+X9AGAOj5Yh5cEydO1LZt2/TGG2/oN7/5jVpaWjR58mSdPn06ZH1LS4vS0tKCjqWlpenSpUs6depU2K9TVFQkt9sd2DIyMmL6OAAA3VPMgysvL0/z58/X6NGj9a1vfUtlZWWSpFdffTXsGMdxgvaNMSGPf1FhYaFaW1sDW1NTUwy6BwB0d3F5jeuL+vbtq9GjR6u+vj7k+fT0dLW0tAQdO3nypBISEnTLLbeEndflcsnlcsW0VwBA9xf393H5/X4dPXpUHo8n5PlJkyapvLw86NjevXs1fvx4JSYmxrs9AIBlYh5cK1asUFVVlRoaGnTw4EF95zvfkc/nU35+vqTPn+JbvHhxoH7p0qU6fvy4CgoKdPToUW3evFmbNm3SihUrYt0aAKAHiPlThf/7v/+rhx9+WKdOnVL//v2Vk5OjAwcOKDMzU5LU3NysxsbGQH1WVpb27NmjJ598Uhs2bJDX69W6deu4FR4AEJJj2u+EsJzP55Pb7dY0zVWCw1OMAGCTS6ZNldqt1tZWpaSkRKzlswoBAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAViG4AABWIbgAAFYhuAAAVol5cA0ePFiO43TYli1bFrK+srIyZP0HH3wQ69YAAD1AQqwnPHTokC5fvhzY/+Mf/6iZM2fqwQcfjDju2LFjSklJCez3798/1q0BAHqAmAfXlwNnzZo1GjJkiKZOnRpxXGpqqm666aZYtwMA6GHi+hrXxYsX9dprr+nRRx+V4zgRa8eOHSuPx6MZM2aooqIinm0BACwW8yuuL9q1a5fOnDmjJUuWhK3xeDx65ZVXlJ2dLb/fr3/5l3/RjBkzVFlZqSlTpoQd5/f75ff7A/s+ny+WrQMAuinHGGPiNfmsWbPUu3dv/fu//3tU4+bMmSPHcVRaWhq2ZvXq1Xr22Wc7HJ+muUpwEqPuFQDQdS6ZNlVqt1pbW4Pudwglbk8VHj9+XPv27dP3vve9qMfm5OSovr4+Yk1hYaFaW1sDW1NT07W2CgCwSNyeKtyyZYtSU1P17W9/O+qxtbW18ng8EWtcLpdcLte1tgcAsFRcguvKlSvasmWL8vPzlZAQ/CUKCwv16aefatu2bZKktWvXavDgwRo5cmTgZo7i4mIVFxfHozUAgOXiElz79u1TY2OjHn300Q7nmpub1djYGNi/ePGiVqxYoU8//VR9+vTRyJEjVVZWptmzZ8ejNQCA5eJ6c0Zn8vl8crvd3JwBABbqFjdnAAAQDwQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCoEFwDAKgQXAMAqBBcAwCpRB1d1dbXmzJkjr9crx3G0a9euoPPGGK1evVper1d9+vTRtGnT9P7771913uLiYo0YMUIul0sjRoxQSUlJtK0BAL4Gog6u8+fPa8yYMVq/fn3I87/85S/1wgsvaP369Tp06JDS09M1c+ZMnT17NuycNTU1WrhwoRYtWqR3331XixYt0oIFC3Tw4MFo2wMA9HCOMcZc82DHUUlJiebNmyfp86str9er5cuXa+XKlZIkv9+vtLQ0Pffcc/rBD34Qcp6FCxfK5/PpP//zPwPH7r33Xt18883avn37V+rF5/PJ7XZrmuYqwUm81ocEAOgCl0ybKrVbra2tSklJiVgb09e4Ghoa1NLSotzc3MAxl8ulqVOnav/+/WHH1dTUBI2RpFmzZkUc4/f75fP5gjYAQM8X0+BqaWmRJKWlpQUdT0tLC5wLNy7aMUVFRXK73YEtIyPjOjoHANgiLncVOo4TtG+M6XDsescUFhaqtbU1sDU1NV17wwAAayTEcrL09HRJn19BeTyewPGTJ092uKL68rgvX11dbYzL5ZLL5brOjgEAtonpFVdWVpbS09NVXl4eOHbx4kVVVVVp8uTJYcdNmjQpaIwk7d27N+IYAMDXU9RXXOfOndNHH30U2G9oaFBdXZ369eunQYMGafny5frFL36hoUOHaujQofrFL36hG2+8UY888khgzOLFizVgwAAVFRVJkp544glNmTJFzz33nObOnavdu3dr3759evvtt2PwEAEAPUnUwXX48GFNnz49sF9QUCBJys/P19atW/XTn/5U//d//6cf/vCH+stf/qKJEydq7969Sk5ODoxpbGxUr15/vdibPHmyXn/9df385z/X008/rSFDhmjHjh2aOHHi9Tw2AEAPdF3v4+pOeB8XANiry97HBQBAvBFcAACrEFwAAKsQXAAAqxBcAACrEFwAAKsQXAAAqxBcAACrEFwAAKsQXAAAqxBcAACrEFwAAKsQXAAAqxBcAACrEFwAAKsQXAAAqxBcAACrEFwAAKsQXAAAqxBcAACrEFwAAKsQXAAAqxBcAACrEFwAAKsQXAAAqxBcAACrEFwAAKsQXAAAqxBcAACrEFwAAKsQXAAAqxBcAACrEFwAAKsQXAAAqxBcAACrRB1c1dXVmjNnjrxerxzH0a5duwLn2tratHLlSo0ePVp9+/aV1+vV4sWLdeLEiYhzbt26VY7jdNguXLgQ9QMCAPRsUQfX+fPnNWbMGK1fv77Duc8++0xHjhzR008/rSNHjmjnzp368MMPdf/991913pSUFDU3NwdtSUlJ0bYHAOjhEqIdkJeXp7y8vJDn3G63ysvLg469+OKLuuOOO9TY2KhBgwaFnddxHKWnp0fbDgDgaybur3G1trbKcRzddNNNEevOnTunzMxMDRw4UPfdd59qa2sj1vv9fvl8vqANANDzxTW4Lly4oFWrVumRRx5RSkpK2Lrhw4dr69atKi0t1fbt25WUlKQ777xT9fX1YccUFRXJ7XYHtoyMjHg8BABAN+MYY8w1D3YclZSUaN68eR3OtbW16cEHH1RjY6MqKysjBteXXblyRePGjdOUKVO0bt26kDV+v19+vz+w7/P5lJGRoWmaqwQnMerHAgDoOpdMmyq1W62trVfNi6hf4/oq2tratGDBAjU0NOjNN9+MKrQkqVevXpowYULEKy6XyyWXy3W9rQIALBPzpwrbQ6u+vl779u3TLbfcEvUcxhjV1dXJ4/HEuj0AgOWivuI6d+6cPvroo8B+Q0OD6urq1K9fP3m9Xn3nO9/RkSNH9B//8R+6fPmyWlpaJEn9+vVT7969JUmLFy/WgAEDVFRUJEl69tlnlZOTo6FDh8rn82ndunWqq6vThg0bYvEYAQA9SNTBdfjwYU2fPj2wX1BQIEnKz8/X6tWrVVpaKkn65je/GTSuoqJC06ZNkyQ1NjaqV6+/XuydOXNG3//+99XS0iK3262xY8equrpad9xxR7TtAQB6uOu6OaM78fl8crvd3JwBABaK5uYMPqsQAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYJWog6u6ulpz5syR1+uV4zjatWtX0PklS5bIcZygLScn56rzFhcXa8SIEXK5XBoxYoRKSkqibQ0A8DUQdXCdP39eY8aM0fr168PW3HvvvWpubg5se/bsiThnTU2NFi5cqEWLFundd9/VokWLtGDBAh08eDDa9gAAPVxCtAPy8vKUl5cXscblcik9Pf0rz7l27VrNnDlThYWFkqTCwkJVVVVp7dq12r59e7QtAgB6sLi8xlVZWanU1FTdfvvteuyxx3Ty5MmI9TU1NcrNzQ06NmvWLO3fvz/sGL/fL5/PF7QBAHq+mAdXXl6efvvb3+rNN9/U888/r0OHDumee+6R3+8PO6alpUVpaWlBx9LS0tTS0hJ2TFFRkdxud2DLyMiI2WMAAHRfUT9VeDULFy4M/PeoUaM0fvx4ZWZmqqysTA888EDYcY7jBO0bYzoc+6LCwkIVFBQE9n0+H+EFAF8DMQ+uL/N4PMrMzFR9fX3YmvT09A5XVydPnuxwFfZFLpdLLpcrZn0CAOwQ9/dxnT59Wk1NTfJ4PGFrJk2apPLy8qBje/fu1eTJk+PdHgDAMlFfcZ07d04fffRRYL+hoUF1dXXq16+f+vXrp9WrV2v+/PnyeDz65JNP9NRTT+nWW2/V3/zN3wTGLF68WAMGDFBRUZEk6YknntCUKVP03HPPae7cudq9e7f27dunt99+OwYPEQDQk0QdXIcPH9b06dMD++2vM+Xn5+vll1/We++9p23btunMmTPyeDyaPn26duzYoeTk5MCYxsZG9er114u9yZMn6/XXX9fPf/5zPf300xoyZIh27NihiRMnXs9jAwD0QI4xxnR1E7Hg8/nkdrs1TXOV4CR2dTsAgChcMm2q1G61trYqJSUlYi2fVQgAsArBBQCwCsEFALAKwQUAsArBBQCwCsEFALAKwQUAsArBBQCwCsEFALAKwQUAsArBBQCwCsEFALAKwQUAsArBBQCwCsEFALAKwQUAsArBBQCwCsEFALAKwQUAsArBBQCwCsEFALAKwQUAsArBBQCwCsEFALAKwQUAsArBBQCwCsEFALAKwQUAsArBBQCwCsEFALAKwQUAsArBBQCwCsEFALAKwQUAsArBBQCwStTBVV1drTlz5sjr9cpxHO3atSvovOM4Ibdf/epXYefcunVryDEXLlyI+gEBAHq2qIPr/PnzGjNmjNavXx/yfHNzc9C2efNmOY6j+fPnR5w3JSWlw9ikpKRo2wMA9HAJ0Q7Iy8tTXl5e2PPp6elB+7t379b06dN12223RZzXcZwOYwEA+LK4vsb15z//WWVlZfrud7971dpz584pMzNTAwcO1H333afa2tqI9X6/Xz6fL2gDAPR8cQ2uV199VcnJyXrggQci1g0fPlxbt25VaWmptm/frqSkJN15552qr68PO6aoqEhutzuwZWRkxLp9AEA35BhjzDUPdhyVlJRo3rx5Ic8PHz5cM2fO1IsvvhjVvFeuXNG4ceM0ZcoUrVu3LmSN3++X3+8P7Pt8PmVkZGia5irBSYzq6wEAutYl06ZK7VZra6tSUlIi1kb9GtdX9dZbb+nYsWPasWNH1GN79eqlCRMmRLzicrlccrlc19MiAMBCcXuqcNOmTcrOztaYMWOiHmuMUV1dnTweTxw6AwDYLOorrnPnzumjjz4K7Dc0NKiurk79+vXToEGDJH3+tN3vfvc7Pf/88yHnWLx4sQYMGKCioiJJ0rPPPqucnBwNHTpUPp9P69atU11dnTZs2HAtjwkA0INFHVyHDx/W9OnTA/sFBQWSpPz8fG3dulWS9Prrr8sYo4cffjjkHI2NjerV668Xe2fOnNH3v/99tbS0yO12a+zYsaqurtYdd9wRbXsAgB7uum7O6E58Pp/cbjc3ZwCAhaK5OYPPKgQAWIXgAgBYheACAFiF4AIAWIXgAgBYheACAFiF4AIAWIXgAgBYheACAFiF4AIAWIXgAgBYJW5/j6urlHz4nlKSry+PZ3m/GZtmAOBr4I0Tddc9h+/sFd18+1er5YoLAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYBWCCwBgFYILAGAVggsAYJUe8xeQjTGSJN+5K9c91yXTdt1zAMDXhe/s9f+72/5vd/u/5ZH0mOA6e/asJClz3CcxmO3jGMwBAF8PN98eu7nOnj0rt9sdscYxXyXeLHDlyhWdOHFCycnJchwnZI3P51NGRoaampqUkpLSyR1eO/rufLb2Tt+di75jxxijs2fPyuv1qlevyK9i9Zgrrl69emngwIFfqTYlJaXbfLOiQd+dz9be6btz0XdsXO1Kqx03ZwAArEJwAQCs8rUKLpfLpWeeeUYul6urW4kKfXc+W3un785F312jx9ycAQD4evhaXXEBAOxHcAEArEJwAQCsQnABAKzS44LrpZdeUlZWlpKSkpSdna233norYn1VVZWys7OVlJSk2267TRs3buykTj9XVFSkCRMmKDk5WampqZo3b56OHTsWcUxlZaUcx+mwffDBB53UtbR69eoOXz89PT3imK5e63aDBw8OuX7Lli0LWd9V611dXa05c+bI6/XKcRzt2rUr6LwxRqtXr5bX61WfPn00bdo0vf/++1edt7i4WCNGjJDL5dKIESNUUlLSaX23tbVp5cqVGj16tPr27Suv16vFixfrxIkTEefcunVryO/BhQsXOqVvSVqyZEmHr5+Tk3PVebtyvSWFXDfHcfSrX/0q7Jydsd7Xo0cF144dO7R8+XL97Gc/U21tre6++27l5eWpsbExZH1DQ4Nmz56tu+++W7W1tXrqqaf04x//WMXFxZ3Wc1VVlZYtW6YDBw6ovLxcly5dUm5urs6fP3/VsceOHVNzc3NgGzp0aCd0/FcjR44M+vrvvfde2NrusNbtDh06FNR3eXm5JOnBBx+MOK6z1/v8+fMaM2aM1q9fH/L8L3/5S73wwgtav369Dh06pPT0dM2cOTPwuZ2h1NTUaOHChVq0aJHeffddLVq0SAsWLNDBgwc7pe/PPvtMR44c0dNPP60jR45o586d+vDDD3X//fdfdd6UlJSg9W9ublZSUlKn9N3u3nvvDfr6e/bsiThnV6+3pA5rtnnzZjmOo/nz50ecN97rfV1MD3LHHXeYpUuXBh0bPny4WbVqVcj6n/70p2b48OFBx37wgx+YnJycuPV4NSdPnjSSTFVVVdiaiooKI8n85S9/6bzGvuSZZ54xY8aM+cr13XGt2z3xxBNmyJAh5sqVKyHPd4f1lmRKSkoC+1euXDHp6elmzZo1gWMXLlwwbrfbbNy4Mew8CxYsMPfee2/QsVmzZpmHHnoo5j0b07HvUP7rv/7LSDLHjx8PW7Nlyxbjdrtj21wEofrOz883c+fOjWqe7rjec+fONffcc0/Ems5e72j1mCuuixcv6p133lFubm7Q8dzcXO3fvz/kmJqamg71s2bN0uHDh9XW1jV/2qS1tVWS1K9fv6vWjh07Vh6PRzNmzFBFRUW8W+ugvr5eXq9XWVlZeuihh/Txx+E/Vb87rrX0+c/Na6+9pkcffTTshzO36+r1/qKGhga1tLQEranL5dLUqVPD/rxL4b8PkcbEW2trqxzH0U033RSx7ty5c8rMzNTAgQN13333qba2tnMa/ILKykqlpqbq9ttv12OPPaaTJ09GrO9u6/3nP/9ZZWVl+u53v3vV2u6w3uH0mOA6deqULl++rLS0tKDjaWlpamlpCTmmpaUlZP2lS5d06tSpuPUajjFGBQUFuuuuuzRq1KiwdR6PR6+88oqKi4u1c+dODRs2TDNmzFB1dXWn9Tpx4kRt27ZNb7zxhn7zm9+opaVFkydP1unTp0PWd7e1brdr1y6dOXNGS5YsCVvTHdb7y9p/pqP5eW8fF+2YeLpw4YJWrVqlRx55JOKHvQ4fPlxbt25VaWmptm/frqSkJN15552qr6/vtF7z8vL029/+Vm+++aaef/55HTp0SPfcc4/8fn/YMd1tvV999VUlJyfrgQceiFjXHdY7kh7z6fDtvvxbszEm4m/SoepDHe8MP/rRj/Tf//3fevvttyPWDRs2TMOGDQvsT5o0SU1NTfr1r3+tKVOmxLtNSZ//T9xu9OjRmjRpkoYMGaJXX31VBQUFIcd0p7Vut2nTJuXl5cnr9Yat6Q7rHU60P+/XOiYe2tra9NBDD+nKlSt66aWXItbm5OQE3Qhx5513aty4cXrxxRe1bt26eLcqSVq4cGHgv0eNGqXx48crMzNTZWVlEYOgu6y3JG3evFl/+7d/e9XXqrrDekfSY664br31Vt1www0dfpM5efJkh9942qWnp4esT0hI0C233BK3XkN5/PHHVVpaqoqKiq/851m+KCcnp0t/G+rbt69Gjx4dtofutNbtjh8/rn379ul73/te1GO7er3b7+CM5ue9fVy0Y+Khra1NCxYsUENDg8rLy6P+0xq9evXShAkTuvR74PF4lJmZGbGH7rLekvTWW2/p2LFj1/Tz3h3W+4t6THD17t1b2dnZgTvE2pWXl2vy5Mkhx0yaNKlD/d69ezV+/HglJibGrdcvMsboRz/6kXbu3Kk333xTWVlZ1zRPbW2tPB5PjLv76vx+v44ePRq2h+6w1l+2ZcsWpaam6tvf/nbUY7t6vbOyspSenh60phcvXlRVVVXYn3cp/Pch0phYaw+t+vp67du375p+cTHGqK6urku/B6dPn1ZTU1PEHrrDerfbtGmTsrOzNWbMmKjHdof1DtJVd4XEw+uvv24SExPNpk2bzJ/+9CezfPly07dvX/PJJ58YY4xZtWqVWbRoUaD+448/NjfeeKN58sknzZ/+9CezadMmk5iYaP7t3/6t03r++7//e+N2u01lZaVpbm4ObJ999lmg5st9/+M//qMpKSkxH374ofnjH/9oVq1aZSSZ4uLiTuv7Jz/5iamsrDQff/yxOXDggLnvvvtMcnJyt17rL7p8+bIZNGiQWblyZYdz3WW9z549a2pra01tba2RZF544QVTW1sbuPtuzZo1xu12m507d5r33nvPPPzww8bj8RifzxeYY9GiRUF31f7hD38wN9xwg1mzZo05evSoWbNmjUlISDAHDhzolL7b2trM/fffbwYOHGjq6uqCfub9fn/YvlevXm1+//vfm//5n/8xtbW15u/+7u9MQkKCOXjwYKf0ffbsWfOTn/zE7N+/3zQ0NJiKigozadIkM2DAgG693u1aW1vNjTfeaF5++eWQc3TFel+PHhVcxhizYcMGk5mZaXr37m3GjRsXdFt5fn6+mTp1alB9ZWWlGTt2rOndu7cZPHhw2G9svEgKuW3ZsiVs388995wZMmSISUpKMjfffLO56667TFlZWaf2vXDhQuPxeExiYqLxer3mgQceMO+//37Yno3p+rX+ojfeeMNIMseOHetwrrusd/tt+F/e8vPzjTGf3xL/zDPPmPT0dONyucyUKVPMe++9FzTH1KlTA/Xtfve735lhw4aZxMREM3z48JgHcKS+Gxoawv7MV1RUhO17+fLlZtCgQaZ3796mf//+Jjc31+zfv7/T+v7ss89Mbm6u6d+/v0lMTDSDBg0y+fn5prGxMWiO7rbe7f7pn/7J9OnTx5w5cybkHF2x3teDP2sCALBKj3mNCwDw9UBwAQCsQnABAKxCcAEArEJwAQCsQnABAKxCcAEArEJwAQCsQnABAKxCcAEArEJwAQCsQnABAKzy/wBVnTUJwLaO3QAAAABJRU5ErkJggg==",
|
|
"text/plain": [
|
|
"<Figure size 640x480 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {},
|
|
"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": 9,
|
|
"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": 10,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def sub(delta_t, delta_x):\n",
|
|
" return delta_t / (2 * delta_x**2)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 38,
|
|
"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",
|
|
" for i_sv in range(0, sv.shape[0]):\n",
|
|
" i_alpha = i_sv-1\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": 39,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Test A:\n",
|
|
"[[ 2.70902018 -0.77802779 0. 0. 0. ]\n",
|
|
" [-0.34610519 2.12413297 -0.77802779 0. 0. ]\n",
|
|
" [ 0. -0.77802779 3.74324914 -1.96522135 0. ]\n",
|
|
" [ 0. 0. -1.96522135 4.21016215 -1.2449408 ]\n",
|
|
" [ 0. 0. 0. 0.43089455 5.13559523]]\n",
|
|
"\n",
|
|
"TestB:\n",
|
|
"[2.45811954 1.70379428 2.70189936 0.92506681 4.0009449 ]\n",
|
|
"\n",
|
|
"Result\n",
|
|
"[1.38151174 1.65087116 1.70267339 1.21472714 0.67714168]\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"row = 5\n",
|
|
"col = 5\n",
|
|
"\n",
|
|
"testConcentrations = np.random.random((row,col))\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": 17,
|
|
"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.11.4"
|
|
},
|
|
"orig_nbformat": 4
|
|
},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 2
|
|
}
|