{ "cells": [ { "cell_type": "markdown", "id": "83da039e", "metadata": {}, "source": [ "# 1D transient advection-diffusion\n", "\n", "Compared to the transient diffusion problem, the transient advection-diffusion problem needs to consider the additional *convective term*." ] }, { "cell_type": "markdown", "id": "df14aa27", "metadata": {}, "source": [ "## Problem setup\n", "\n", "Consider a one-dimensional transient advection-diffusion problem, in which a field variable $\\phi$ is transported through the advection-diffusion process from $x = 0$ to $x = L$ in a one-dimensional domain. The velocity $u = 2.0$ m/s, length $L = 1.5$ m, density $\\rho = 1.0$ kg/m$^3$ and the diffusion coefficient $\\Gamma =0.03$ kg/(m$\\cdot$K).\n", "\n", "![domain](images/transient_advection_diffusion1d.jpg)\n", "\n", "The mathematical model for one-dimensional steady-state advection-diffusion problem is\n", "\n", "$$\\frac{\\partial(\\rho\\phi)}{\\partial t} + \\frac{\\partial(\\rho u\\phi)}{\\partial x} = \\frac{\\partial}{\\partial x}\\left(\\Gamma \\frac{\\partial\\phi}{\\partial x}\\right) + S.$$\n", "\n", "The boundary conditions are\n", "\n", "$$\\phi|_{x=0} = 0, \\quad \\frac{\\partial \\phi}{\\partial x}|_{x=L} = 0.$$\n", "\n", "The initial condition is $\\phi = 0$ everywhere at $t = 0$ s.\n", "\n", "The distribution of source term is like the picture below, \n", "\n", "![source](images/transient_advection_diffusion1d_source_term.jpg)\n", "\n", "where $x_1 = 0.6$ m, $x_2 = 0.2$ m, $a = -200$, $b = 100$. Solve the solution of $\\phi$ until the steady state is achieved." ] }, { "cell_type": "markdown", "id": "925dfd3d", "metadata": {}, "source": [ "## Solve problem\n", "\n", "### Define grid\n", "\n", "Divide the rod evenly into 45 control volumes, as a result, the length of each control volume becomes $\\delta x = 0.0333$ m.\n", "\n", "### Discretize 1D transient advection-diffusion equation\n", "\n", "We solve the problem by the ***modified QUICK scheme***. The parameter of the problem is $u = 2.0$ m/s, $\\delta x = 0.0333$ m, $F = \\rho u = 2.0$, $D = \\frac{\\Gamma}{\\delta x} = 0.9$. The modified QUICK scheme uses the following formulas to calculate the field variable values at each interface of the control volume.\n", "\n", "$$\\phi_e = \\phi_P + \\frac{1}{8} \\left( 3\\phi_E - 2\\phi_P - \\phi_W \\right)$$\n", "\n", "$$\\phi_w = \\phi_W + \\frac{1}{8} \\left( 3\\phi_P - 2\\phi_W - \\phi_WW \\right)$$\n", "\n", "The fully implicit discrete equation form at a general internal node is\n", "\n", "$$\\rho(\\phi_P - \\phi_P^0) \\frac{\\delta x}{\\delta t} + F_e \\left[ \\phi_P + \\frac{1}{8} \\left( 3\\phi_E - 2\\phi_P - \\phi_W \\right) \\right] -F_w \\left[ \\phi_W + \\frac{1}{8} (3\\phi_P - 2\\phi_W - \\phi_{WW}) \\right] = D_e(\\phi_E - \\phi_P) - D_W(\\phi_P - \\phi_W) + \\int_{\\delta x} S \\mathrm{d}x$$\n", "\n", "where\n", "\n", "$$\\int_{\\delta x} S \\, \\mathrm{d}x =\n", "\\begin{cases}\n", "\\int_0^{0.6} \\left( ax + b \\right) \\mathrm{d}x = 24, & x\\leq0.6 \\\\\n", "\\int_{0.6}^{0.8} \\left( 100x - 80 \\right) \\mathrm{d}x = -2, & 0.6 < x \\leq0.8 \\\\\n", "0, & x > 0.8\n", "\\end{cases}$$\n", "\n", "Special handling is required for the first and last nodes. A mirrored outer point should be set on the west outer boundary of the control volume where node 1 is located. Since $\\phi_A = 0$ at the boundary $x = 0$, the field variable value at the mirror point of linear extrapolation is\n", "\n", "$$\\phi_0 = -\\phi_P$$\n", "\n", "The flow rate at the left boundary due to diffusion is\n", "\n", "$$\\Gamma\\left. \\frac{\\partial \\phi}{\\partial x} \\right|_A = \\frac{\\phi_A}{3} (9\\phi_P - 8\\phi_A - \\phi_E)$$\n", "\n", "Therefore, the discrete equation of node 1 is\n", "\n", "$$\\rho(\\phi_P - \\phi_P^0) \\frac{\\delta x}{\\delta t} + F_e \\left[ \\phi_P + \\frac{1}{8}(3\\phi_E - \\phi_P) \\right] - F_A\\phi_A = D_e(\\phi_E - \\phi_P) - \\frac{D_A}{3} (9\\phi_P - 8\\phi_A - \\phi_E)$$\n", "\n", "At to the last node, the gradient of $\\phi$ is 0, so\n", "\n", "$$\\phi_B = \\phi_P$$\n", "\n", "Therefore, the discrete equation of the last node is\n", "\n", "$$\\rho(\\phi_P - \\phi_P^0) \\frac{\\delta x}{\\delta t} + F_B\\phi_P - F_W\\left[ \\phi_W + \\frac{1}{8} (3\\phi_P - 2\\phi_W - \\phi_{WW}) \\right] = 0 - D_W(\\phi_P - \\phi_W)$$\n", "\n", "We can summarize the form of the discrete equations and get\n", "\n", "$$a_P\\phi_P = a_W\\phi_W + a_E\\phi_E + a_P^0\\phi_P^0 + S_u$$\n", "\n", "$$a_P = a_W + a_E + a_P^0 + (F_e - F_w)-S_P, \\quad a_P^0 = \\rho\\frac{\\delta x}{\\delta t}$$\n", "\n", "The coefficients are summarized in the following table\n", "\n", "| node | $a_{W}$ | $a_{E}$ | $S_{P}$ | $S_{u}$ |\n", "| :---: | :----: | :----: | :----: | :----: |\n", "| 1 | 0 | $D_e + \\frac{D_A}{3}$ | $- \\left( \\frac{8}{3} D_A + F_A \\right)$ | $\\left( \\frac{8}{3} D_A + F_A \\right) \\phi_A + \\frac{1}{8} F_e (\\phi_P - 3\\phi_E) + S$ |\n", "| 2 | $D_w + F_w$ | $D_e$ | 0 | $\\frac{1}{8} F_w (3\\phi_P - \\phi_W) + \\frac{1}{8} F_e (\\phi_W + 2\\phi_P - 3\\phi_E) + S$ |\n", "| 3-44 | $D_w + F_w$ | $D_e$ | 0 | $\\frac{1}{8} F_w (3\\phi_P - 2\\phi_W - \\phi_{WW}) + \\frac{1}{8} F_e (\\phi_W + 2\\phi_P - 3\\phi_E) + S$ |\n", "| 45 | $D_w + F_w$ | 0 | 0 | $\\frac{1}{8} F_w (3\\phi_P - 2\\phi_W - \\phi_{WW})$ |\n", "\n", "If the difference between the calculation results of two consecutive time steps is small enough (such as $10^{-6}$), it is considered that a stable-state has reached." ] }, { "cell_type": "markdown", "id": "b8105d81", "metadata": {}, "source": [ "### Solve algebraic equations (implicit method)" ] }, { "cell_type": "code", "execution_count": 17, "id": "04a4737e", "metadata": {}, "outputs": [], "source": [ "# Parameter declarations\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "n = 45\n", "ph0 = 0.0\n", "dt = 0.01\n", "L = 1.5\n", "dx = L / n\n", "x = np.linspace(0.5*dx, L-0.5*dx, n)\n", "gamma = 0.03\n", "rou = 1.0\n", "u = 2.0\n", "F = rou * u\n", "D = gamma / dx\n", "ap0 = rou * dx / dt" ] }, { "cell_type": "code", "execution_count": 18, "id": "71ea3c49", "metadata": {}, "outputs": [], "source": [ "# Initialize arrays\n", "a = np.zeros((n, n))\n", "b = np.zeros(n)\n", "phi = np.zeros(n)\n", "phi0 = np.zeros(n)\n", "time = 0.0" ] }, { "cell_type": "code", "execution_count": 19, "id": "bac63e3a", "metadata": {}, "outputs": [], "source": [ "# Renew coefficients\n", "def renew_implicit(a, b, n):\n", " a[0][0] = ap0 + D + D/3.0 + (8.0/3.0*D + F)\n", " a[0][1] = -(D + D/3.0)\n", " a[1][0] = -(D + F)\n", " a[1][1] = 2*D + F + ap0\n", " a[1][2] = -D\n", " \n", " for i in range(2, n - 1):\n", " a[i][i - 1] = -(D + F)\n", " a[i][i] = 2 * D + F + ap0\n", " a[i][i + 1] = -D\n", " \n", " a[n-1][n-2] = -(D + F)\n", " a[n-1][n-1] = ap0 + D + F\n", " \n", " b[0] = ap0*phi0[0] + 1.0/8.0*F*(phi0[0] - 3*phi0[1]) + 24\n", " b[1] = ap0*phi0[1] + 1.0/8.0*F*(3*phi0[1] - phi0[0]) + 1.0/8.0*F*(phi0[0] + 2*phi0[1] - 3*phi0[2]) + 24\n", " \n", " for i in range(2, n-1):\n", " pos = (i-1)*dx + dx/2.0\n", " if pos <= 0.6:\n", " S = 24\n", " elif 0.6 < pos <= 0.8:\n", " S = -2\n", " else:\n", " S = 0\n", " \n", " b[i] = (ap0*phi0[i] + 1.0/8.0*F*(3*phi0[i] - 2*phi0[i-1] - phi0[i-2]) + 1.0/8.0*F*(phi0[i-1] + 2*phi0[i] - 3*phi0[i+1]) + S)\n", " \n", " b[n-1] = ap0*phi0[n-1] + 1.0/8.0*F*(3*phi0[n-1] - 2*phi0[n-2] - phi0[n-3])" ] }, { "cell_type": "code", "execution_count": 20, "id": "4eb166d1", "metadata": {}, "outputs": [], "source": [ "# Output information\n", "def output():\n", " print(\"-----------\")\n", " print(\"time = {}\".format(time))\n", " for i in range(0, n):\n", " print(phi[i])" ] }, { "cell_type": "code", "execution_count": 21, "id": "2344cf81", "metadata": {}, "outputs": [], "source": [ "# Iteration using TDMA method\n", "def TDMA(a, b, T, nx):\n", " C = np.zeros(nx)\n", " phi = np.zeros(nx)\n", " alph = np.zeros(nx)\n", " belt = np.zeros(nx)\n", " D = np.zeros(nx)\n", " A = np.zeros(nx)\n", " Cpi = np.zeros(nx)\n", " \n", " for j in range(0, nx):\n", " belt[j] = -a[j][j-1]\n", " D[j] = a[j][j]\n", " alph[j] = -a[j][j+1] if j < nx-1 else 0\n", " C[j] = b[j]\n", " \n", " for j in range(0, nx):\n", " denom = D[j] - belt[j] * A[j-1]\n", " A[j] = alph[j] / denom if denom != 0 else 0\n", " Cpi[j] = (belt[j] * Cpi[j-1] + C[j]) / denom if denom != 0 else 0\n", " \n", " phi[nx-1] = Cpi[nx-1]\n", " for j in range(nx-2, -1, -1):\n", " phi[j] = A[j]*phi[j+1] + Cpi[j]\n", " \n", " for j in range(0, nx):\n", " T[j] = phi[j]" ] }, { "cell_type": "code", "execution_count": 22, "id": "9750c7a5", "metadata": {}, "outputs": [], "source": [ "# Iteration using Jacobi method\n", "def Jacobi(A, b, phi, k=100):\n", " n = A.shape[1]\n", " D = np.eye(n)\n", " D[np.arange(n), np.arange(n)] = A[np.arange(n), np.arange(n)]\n", " LU = D - A\n", " X = np.zeros(n)\n", "\n", " for i in range(k):\n", " D_inv = np.linalg.inv(D)\n", " X = np.dot(np.dot(D_inv, LU), X) + np.dot(D_inv, b)\n", " \n", " for i in range(0 ,n):\n", " phi[i] = X[i]" ] }, { "cell_type": "code", "execution_count": 23, "id": "61c36ecc", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-----------\n", "time = 1.5100000000000011\n", "6.0\n", "18.000000000000007\n", "30.00000000000001\n", "42.000000000000014\n", "54.00000000000003\n", "66.00000000000003\n", "78.00000000000003\n", "90.0\n", "101.99999999999932\n", "113.99999999998798\n", "125.99999999980719\n", "137.99999999693256\n", "149.99999995123852\n", "161.9999992249244\n", "173.99998768005463\n", "185.99980417266147\n", "197.9968872958205\n", "209.95052311201937\n", "221.21355763260254\n", "221.4993836016016\n", "220.6342082379482\n", "219.64835913782747\n", "218.65006580865105\n", "217.6537877086974\n", "216.71049365806644\n", "216.61158568122656\n", "216.6012148028397\n", "216.60012737671096\n", "216.60001335593998\n", "216.60000140043255\n", "216.60000014685198\n", "216.60000001538333\n", "216.60000000150123\n", "216.59999999981005\n", "216.59999999922226\n", "216.59999999875413\n", "216.5999999991413\n", "216.60000000280468\n", "216.60000001508013\n", "216.60000004536715\n", "216.6000001059163\n", "216.6000002021355\n", "216.6000002982933\n", "216.60000038924136\n", "216.59999858406954\n" ] } ], "source": [ "# Calculate the solution\n", "for i in range(0, n):\n", " phi[i] = ph0\n", " phi0[i] = ph0\n", "\n", "time = 0.0\n", "go = 1\n", "\n", "while go == 1:\n", " go = 0\n", " time += dt\n", " renew_implicit(a, b, n)\n", " # TDMA(a, b, phi, n)\n", " Jacobi(a, b, phi)\n", " \n", " for i in range(0, n):\n", " if abs(phi[i] - phi0[i]) > 1.0e-6:\n", " go = 1\n", " \n", " for i in range(0, n):\n", " phi0[i] = phi[i]\n", " \n", "output()" ] }, { "cell_type": "code", "execution_count": 24, "id": "adde66b6", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj8AAAGyCAYAAAALaqWsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA91ElEQVR4nO3deVxU9eL/8fewCgokoiCBuKFWopZmanm1RU3T8vYt63bbrfSaECq5gPtGbqhg2XIt27Td6paV3hbTrJuZ5tJyXTA1JXfADQTO7w9/zo3EBZyZzyyv5+Mxj+DMZ8b3+QRn3pw5c47NsixLAAAAPsLPdAAAAABXovwAAACfQvkBAAA+hfIDAAB8CuUHAAD4FMoPAADwKZQfAADgUwJMB3BHZWVl2rVrl8LCwmSz2UzHAQAA58GyLBUWFio2NlZ+fmfev0P5qcCuXbsUHx9vOgYAAKiCHTt2KC4u7oz3U34qEBYWJunk5IWHhxtOAwAAzkdBQYHi4+Ptr+NnQvmpwKm3usLDwyk/AAB4mHMdssIBzwAAwKdQfgAAgE+h/AAAAJ9C+QEAAD6F8gMAAHwK5QcAAPgUyg8AAPAplB8AAOBTKD8AAMCnUH4AAIBPofwAAACfQvkBAAA+hQubAnAqy7JUXFyso0eP6ujRowoPDz/nFZcBwJkoPwCqrKioSMOGDdPmzZt15MgRe8H549dHjx5VaWlpucdFRkYqISFBCQkJql+//mlf16xZ85xXZQaAqqL8AKiyyZMna/bs2ec9PiAgQCUlJTpw4IAOHDigNWvWVDiuRo0aql+/vho2bKhGjRqpcePGaty4sRo1aqSEhAQFBLDpAlB1NsuyLNMh3E1BQYEiIiKUn5+v8PBw03EAt/TTTz+pZcuWOnHihEaOHKnmzZsrNDRUoaGhql69ern/nroFBgaqsLBQv/76q7Zt21buv6e+3rNnz1n/3YCAANWvX99ehk4Vo8aNG6tBgwYKDg520QwAcDfn+/pN+akA5Qc4u7KyMnXu3FnLly9Xjx499MEHHzjsbaqjR49q+/bt+vXXX7V161Zt3rzZftu6dauOHz9+xsfabDbVq1evXCH6416jkJAQh2QE4J4oPxeA8gOc3bx58/TQQw8pNDRUGzduVP369V3y75aVlWnXrl3avHmztmzZUq4Ybd68WYcPHz7r4+Pi4s5YjGrUqOGSdQDgPJSfC0D5Ac5sz549atasmQ4ePKhp06YpLS3NdCRJJz9VtnfvXm3atOm0UrRp0ybl5+ef9fExMTFKTEyssBhFRES4aC0AXAjKzwWg/ABndvfdd+vVV19Vq1attGrVKo84+NiyLB04cMBehE7tNTpVlPbv33/Wx0dGRqpRo0Zq2LDhabe4uDiPmAPAF1B+LgDlB6jY0qVL1bVrV/n5+embb77RlVdeaTqSQxw8eLDCt9E2bdp0XgdgJyQkqFGjRqpVq5YCAgLk7+9/1v8GBATIz+9/55g9dbzUuf4LoLwuXbqobdu29u/P9/WbP1cAnJejR4+qf//+kqSBAwd6TfGRpJo1a6pNmzZq06bNafcdPnxYubm52rJli7Zu3Vrulpubq+LiYm3ZskVbtmwxkBzwbTVq1ChXfs4X5QfAeZk4caK2bt2qiy++WBMmTDAdx2Vq1KihpKQkJSUlnXZfaWmpdu3aZS9D+fn5Ki0tVUlJiUpKSir8urS0VCdOnNCpne5/3Ple0dfsnAfOrHnz5lV6HG97VYC3vYDy1q9fryuuuEIlJSVatGiRevfubToSAJzmfF+/ubApgLMqKytTv379VFJSot69e1N8AHg8yg+As3r22Wf19ddfq0aNGsrJyTEdBwAuGOUHwBnt3r1bw4cPlyRNmjRJcXFxhhMBwIWj/AA4o9TUVOXn5+vKK6/Uo48+ajoOADgE5QdAhRYvXqw33nhD/v7+evbZZ+Xv7286EgA4BOUHwGmOHDmiAQMGSDq596dVq1ZmAwGAA1F+AJxm7Nix+vXXX5WQkKBx48aZjgMADkX5AVDOunXrNHPmTEnSU089perVqxtOBACORfkBYGdZlgYOHKjS0lLddttt6tGjh+lIAOBwlB8AdgsXLtTy5csVGhqqGTNmmI4DAE5B+QEgSSosLNTjjz8uSUpPT1e9evUMJwIA56D8AJB08sKlu3btUsOGDTVkyBDTcQDAaSg/APTLL7/YD3KePXu2qlWrZjgRADgP5QfwcZZl6bHHHtOJEyfUo0cP9ezZ03QkAHAqyg/g495//3198sknCgoK0qxZs0zHAQCno/wAPuzYsWMaNGiQJGnIkCFKTEw0nAgAnI/yA/iwadOmKTc3VxdffLHS09NNxwEAl6D8AD5q27ZtyszMlCTNmDFDNWrUMJwIAFyD8gP4qCFDhuj48ePq3Lmz+vTpYzoOALgM5QfwQUuXLtU777wjf39/ZWdny2azmY4EAC7jVuUnMzNTV155pcLCwlSnTh317t1bv/zyS7kxlmVp7Nixio2NVUhIiDp37qyNGzeWG1NUVKTk5GRFRUWpevXquvnmm7Vz505XrgrgtoqLi5WSkiJJevTRR5WUlGQ4EQC4lluVn2XLlunRRx/VN998o6VLl6qkpERdu3bVkSNH7GOmTp2qrKwszZkzR6tWrVJMTIy6dOmiwsJC+5jU1FQtWrRIr732mlasWKHDhw+rZ8+eKi0tNbFagFvJycnRzz//rNq1a2vcuHGm4wCAy9ksy7JMhziTvXv3qk6dOlq2bJn+8pe/yLIsxcbGKjU1VcOGDZN0ci9PdHS0pkyZon79+ik/P1+1a9fWyy+/rDvuuEOStGvXLsXHx2vx4sXq1q3baf9OUVGRioqK7N8XFBQoPj5e+fn5Cg8Pd83KAi6we/duNW3aVIWFhfrnP/+pvn37mo4EAA5TUFCgiIiIc75+u9Wenz/Lz8+XJEVGRkqScnNzlZeXp65du9rHBAcHq1OnTlq5cqUkafXq1Tpx4kS5MbGxsWrevLl9zJ9lZmYqIiLCfouPj3fWKgFGDRs2TIWFhbryyiv1wAMPmI4DAEa4bfmxLEuDBw/WNddco+bNm0uS8vLyJEnR0dHlxkZHR9vvy8vLU1BQkGrWrHnGMX82YsQI5efn2287duxw9OoAxn311Vd6+eWXJUlz5syRn5/b/voDgFMFmA5wJgMHDtS6deu0YsWK0+778ydTLMs656dVzjYmODhYwcHBVQ8LuLnS0lINHDhQkvTggw+qbdu2hhMBgDlu+adfcnKy3n//fX3++eeKi4uzL4+JiZGk0/bg7Nmzx743KCYmRsXFxTp48OAZxwC+5rnnntPatWsVERFhP7EhAPgqtyo/lmVp4MCBeuedd/TZZ5+pQYMG5e5v0KCBYmJitHTpUvuy4uJiLVu2TB06dJAktW7dWoGBgeXG7N69Wxs2bLCPAXzJ/v37lZGRIUkaP3686tSpYzgRAJjlVm97Pfroo1qwYIHee+89hYWF2ffwREREKCQkRDabTampqZo8ebISExOVmJioyZMnKzQ0VHfddZd9bN++fTVkyBDVqlVLkZGRSktLU1JSkm644QaTqwcYMWrUKB04cEDNmzfXgAEDTMcBAOPcqvzMnTtXktS5c+dyy1944QXdf//9kqShQ4fq2LFjGjBggA4ePKirrrpKS5YsUVhYmH38zJkzFRAQoD59+ujYsWO6/vrrNX/+fPn7+7tqVQC3sHbtWj3zzDOSpOzsbAUEuNWvPAAY4dbn+THlfM8TALgzy7LUqVMnLV++XLfffrveeOMN05EAwKm84jw/AKpu4cKFWr58uUJCQjR9+nTTcQDAbVB+AC90+PBhPf7445Kk9PR01atXz3AiAHAflB/AC02cOFG7du1Sw4YNlZaWZjoOALgVyg/gZTZt2qSsrCxJJw/+r1atmuFEAOBeKD+Al0lNTdWJEyd04403qlevXqbjAIDbofwAXuSDDz7Q4sWLFRgYqNmzZ5/zsi8A4IsoP4CXOH78uFJTUyVJgwYNUpMmTcwGAgA3RfkBvERWVpa2bNmiunXrauTIkabjAIDbovwAXmDnzp2aNGmSJGnq1KnlzngOACiP8gN4gccff1xHjx7V1Vdfrb///e+m4wCAW6P8AB5u2bJleu2112Sz2ZSTk8NBzgBwDpQfwIOVlJQoOTlZktSvXz9dfvnlhhMBgPuj/AAe7Omnn9b69esVGRmpiRMnmo4DAB6B8gN4qH379mnUqFGSTl7OolatWoYTAYBnoPwAHiojI0OHDh1Sy5Yt9cgjj5iOAwAeg/IDeKDVq1frueeekyTNmTNH/v7+hhMBgOeg/AAepqysTMnJybIsS3fddZeuueYa05EAwKNQfgAP88orr+jrr79W9erVNXXqVNNxAMDjUH4AD1JQUKBhw4ZJkkaNGqWLL77YcCIA8DyUH8CDTJgwQXl5eUpMTLRfxBQAUDmUH8BD/Pzzz5o1a5Ykafbs2QoODjYbCAA8FOUH8ACWZemxxx5TSUmJevbsqe7du5uOBAAei/IDeID33ntPS5YsUVBQkGbOnGk6DgB4NMoP4OaOHTumQYMGSZLS0tLUuHFjw4kAwLNRfgA3N336dG3btk1xcXFKT083HQcAPB7lB3Bjv/76qzIzMyWdLEHVq1c3nAgAPB/lB3BjaWlpOnbsmDp16qQ+ffqYjgMAXoHyA7ipTz/9VG+99Zb8/PyUnZ0tm81mOhIAeAXKD+CGTpw4oZSUFEnSgAED1KJFC8OJAMB7UH4AN/Tkk0/qxx9/VFRUlMaPH286DgB4FcoP4Gb27NmjMWPGSJImT56smjVrGk4EAN6F8gO4mREjRqigoECtW7fWgw8+aDoOAHgdyg/gRr799ls9//zzkqScnBz5+/sbTgQA3ofyA7iJsrIyDRw4UJJ07733qn379oYTAYB3ovwAbmL+/PlatWqVwsLCNGXKFNNxAMBrUX4AN3Do0CENHz5ckjRmzBjFxMQYTgQA3ovyA7iBcePGae/evWrWrJmSk5NNxwEAr0b5AQzbuHGjcnJyJEnZ2dkKCgoynAgAvBvlBzDIsiylpKSotLRUvXv3VpcuXUxHAgCvR/kBDHr77bf12WefqVq1asrKyjIdBwB8AuUHMOTo0aMaMmSIJGno0KFq0KCB4UQA4BsoP4AhU6ZM0fbt21WvXj0NGzbMdBwA8BmUH8CA3Nxc+7l8srKyFBoaajgRAPgOyg9gwODBg1VUVKTrrrtOt956q+k4AOBTKD+Aiy1ZskTvvvuu/P39lZ2dLZvNZjoSAPgUyg/gQsXFxUpJSZEkJScn67LLLjOcCAB8D+UHcKGcnBz98ssvqlOnjsaOHWs6DgD4JMoP4CK7d+/WuHHjJElPPPGEIiIiDCcCAN9E+QFcZPjw4SosLFTbtm113333mY4DAD6L8gO4wMqVK/XSSy9JOvnWl58fv3oAYApbYMDJSktL7Vdqf/DBB9W2bVvDiQDAt1F+ACebN2+evv/+e0VERCgzM9N0HADweZQfwIkOHjyo9PR0SdLYsWNVp04dw4kAAJQfwIlGjx6t/fv369JLL9Wjjz5qOg4AQJQfwGnWrVunp556StLJg5wDAwMNJwIASJQfwCksy1JycrLKysp022236brrrjMdCQDw/1F+ACd4/fXX9eWXXyokJEQzZswwHQcA8AeUH8DBDh8+rLS0NEnSiBEjVK9ePcOJAAB/RPkBHCwzM1O//fab6tevby9BAAD3QfkBHGjz5s2aPn26JGnmzJkKCQkxnAgA8GeUH8CBBg0apOLiYnXt2lW33HKL6TgAgApQfgAHWbx4sT744AMFBARo9uzZstlspiMBACpA+QEcoKioSI899pgkKTU1Vc2aNTOcCABwJpQfwAFmzZqlzZs3KyYmRqNGjTIdBwBwFpQf4AL99ttvmjBhgiRpypQpCg8PN5wIAHA2lB/gAg0dOlRHjhxR+/btdffdd5uOAwA4B8oPcAGWL1+uBQsWyGazac6cOfLz41cKANydW22pv/zyS/Xq1UuxsbGy2Wx69913y91///33y2azlbu1a9eu3JiioiIlJycrKipK1atX180336ydO3e6cC3gK0pLS5WcnCxJevjhh3XFFVcYTgQAOB9uVX6OHDmili1bas6cOWccc+ONN2r37t322+LFi8vdn5qaqkWLFum1117TihUrdPjwYfXs2VOlpaXOjg8f88wzz+iHH37QRRddpEmTJpmOAwA4TwGmA/xR9+7d1b1797OOCQ4OVkxMTIX35efna968eXr55Zd1ww03SJJeeeUVxcfH69///re6detW4eOKiopUVFRk/76goKCKawBfsX//fo0cOVKSNGHCBEVFRRlOBAA4X2615+d8fPHFF6pTp46aNGmihx9+WHv27LHft3r1ap04cUJdu3a1L4uNjVXz5s21cuXKMz5nZmamIiIi7Lf4+HinrgM838iRI3Xw4EElJSWpf//+puMAACrBo8pP9+7d9eqrr+qzzz7TjBkztGrVKl133XX2vTZ5eXkKCgpSzZo1yz0uOjpaeXl5Z3zeESNGKD8/337bsWOHU9cDnm3NmjV65plnJEk5OTkKCHCrHagAgHPwqK32HXfcYf+6efPmatOmjRISEvThhx/q1ltvPePjLMs666UGgoODFRwc7NCs8E6WZSk5OVmWZenOO+9Up06dTEcCAFSSR+35+bO6desqISFBmzZtkiTFxMSouLhYBw8eLDduz549io6ONhERXubVV1/VV199pdDQUE2bNs10HABAFXh0+dm/f7927NihunXrSpJat26twMBALV261D5m9+7d2rBhgzp06GAqJrxEYWGhhg4dKunkMT9xcXGGEwEAqsKt3vY6fPiwNm/ebP8+NzdXa9euVWRkpCIjIzV27Fj93//9n+rWratt27YpPT1dUVFR+utf/ypJioiIUN++fTVkyBDVqlVLkZGRSktLU1JSkv3TX0BVTZw4Ubt371ajRo00ePBg03EAAFXkVuXnu+++07XXXmv//tQLzH333ae5c+dq/fr1eumll3To0CHVrVtX1157rV5//XWFhYXZHzNz5kwFBASoT58+OnbsmK6//nrNnz9f/v7+Ll8feI9ffvlFM2fOlHTyIqYcIwYAnstmWZZlOoS7KSgoUEREhPLz87lIJWRZlnr06KGPP/5YPXr00Icffmg6EgCgAuf7+u3Rx/wArvCvf/1LH3/8sYKCgjRr1izTcQAAF4jyA5zF8ePHNWjQIEkn34ZNTEw0nAgAcKEoP8BZzJgxQ1u3blVsbKwyMjJMxwEAOADlBziDHTt2aPLkyZKkadOmqUaNGoYTAQAcgfIDnEFaWpqOHj2qjh076m9/+5vpOAAAB6H8ABX4/PPP9cYbb8jPz085OTlnvTwKAMCzUH6APykpKVFKSookqX///mrZsqXhRAAAR6L8AH8yd+5cbdiwQZGRkZowYYLpOAAAB6P8AH+wd+9ejR49WpI0adIkRUZGGk4EAHA0yg/wB+np6Tp06JAuv/xyPfzww6bjAACcgPID/H/fffed5s2bJ0nKycnhenAA4KUoP4CksrIyJScny7Is3X333br66qtNRwIAOAnlB5D08ssv65tvvlGNGjU0ZcoU03EAAE5E+YHPy8/P17BhwyRJo0aNUmxsrOFEAABnovzA540fP16///67mjRpotTUVNNxAABORvmBT/vpp5+UnZ0tSZo9e7aCgoIMJwIAOBvlBz7LsiylpKSopKREN998s2688UbTkQAALkD5gc9699139e9//1vBwcGaOXOm6TgAABeh/MAnHTt2TIMHD5Z08urtDRs2NJwIAOAqlB/4pKlTp2rbtm2Kj4/XiBEjTMcBALgQ5Qc+Z9u2bXriiSckSdOnT1f16tUNJwIAuBLlBz5nyJAhOn78uDp37qzbb7/ddBwAgItRfuBT/v3vf+udd96Rv7+/srOzZbPZTEcCALgY5Qc+48SJE0pJSZEkDRgwQElJSYYTAQBMoPzAZ8yZM0c//fSToqKiNG7cONNxAACGUH7gE37//XeNHTtWkpSZmamaNWuaDQQAMIbyA58wfPhwFRQUqHXr1nrggQdMxwEAGET5gdf7z3/+o/nz50s6+daXv7+/2UAAAKMoP/BqZWVlSk5OliTdd999ateuneFEAADTKD/wai+88IJWrVqlsLAw+4kNAQC+jfIDr3Xo0CH7pSvGjBmjmJgYw4kAAO6A8gOvNWbMGO3du1fNmjWzv/UFAADlB15pw4YNevLJJyVJ2dnZCgoKMpwIAOAuKD/wOpZlKSUlRaWlpfrrX/+qLl26mI4EAHAjlB94nbfeekuff/65qlWrpqysLNNxAABuhvIDr3LkyBENGTJEkjRs2DDVr1/fbCAAgNuh/MCrPPHEE9qxY4cSEhI0bNgw03EAAG6I8gOvsXXrVk2bNk2SlJWVpZCQEMOJAADuiPIDrzF48GAVFRXp+uuv11//+lfTcQAAboryA6/wySef6L333lNAQICys7Nls9lMRwIAuCnKDzxecXGxUlJSJEnJycm69NJLDScCALizSpWfu+++W8eOHZMk7dixwymBgMqaPXu2/vvf/yo6OlpjxowxHQcA4OYCKjO4Ro0aKioqUkhIiBISElSzZk21bNlSLVu2VKtWrdSyZUtddtllCgwMdFZeoJxdu3Zp/Pjxkk5+0isiIsJwIgCAu6tU+Xn66aftX+fm5mrt2rX64YcftHbtWr3//vvatm2bAgIC1KxZM/3www8ODwv82fDhw3X48GG1a9dO9957r+k4AAAPUKny80cJCQlKSEjQLbfcYl9WWFiotWvXat26dQ4JB5zNypUr9fLLL8tmsyknJ0d+fhzCBgA4tyqXn4qEhYWpY8eO6tixoyOfFjhNaWmpBg4cKEnq27ev2rRpYzgRAMBTVLn8rFq1SsOHD9fevXvVuHFjtWrVyn6rV6+eIzMCp/nnP/+pNWvWKCIiQpMmTTIdBwDgQar8PsE999wjf39/9e/fXw0bNtSyZcv0wAMPqH79+qpVq5YjMwLlHDhwQBkZGZKk8ePHq06dOoYTAQA8SZX3/OzYsUMffvihGjVqVG75r7/+qrVr115oLuCMRo8erf3796t58+YaMGCA6TgAAA9T5fJz9dVXa8eOHaeVn1MHQgPO8MMPP2ju3LmSpOzsbAUEOPSwNQCAD6jUK8ctt9xiP69P//79NX78eCUlJfE2F1zCsiwlJyerrKxMffr00bXXXms6EgDAA1Wq/CQmJmrlypWaO3eu9u/fL0lq2rSpbrnlFrVv316XX365kpKSFBQU5JSw8G2vvfaali9frpCQEE2fPt10HACAh7JZlmVV5YE7d+7U2rVry91yc3Pl7++vZs2aefS5fgoKChQREaH8/HyFh4ebjgNJhw8fVtOmTbVr1y5NmDBBI0eONB0JAOBmzvf1u8oHTMTFxSkuLk6XX365JOniiy/W4cOHtWbNGo8uPnBPkyZN0q5du9SwYUOlpaWZjgMA8GBV/qj7V199pQYNGqhevXqqV6+eoqOjNWHCBLVs2VKPPvqoIzPCx23atElZWVmSpJkzZ6patWqGEwEAPFmVy0+/fv102WWXadWqVVq3bp2mTZumTz/9VK1bt9a+ffscmRE+btCgQSouLtaNN96oXr16mY4DAPBwVT7mJyQkROvWrVNiYqJ9mWVZ6tOnjwIDA7VgwQKHhXQ1jvlxHx9++KF69uypwMBArV+/Xk2bNjUdCQDgps739bvKe34uueQS5eXllVtms9k0fvx4/etf/6rq0wJ2RUVFSk1NlSSlpqZSfAAADlHl8nP//ffrkUce0fbt28stz8/PV0RExAUHA7KysrR582bVrVtXo0aNMh0HAOAlqvxpr1N/kTdp0kS33nqrWrVqpdLSUr3yyiuaNm2ao/LBR+3cuVMTJ06UJE2dOlVhYWGGEwEAvEWVj/nZs2eP1qxZox9++MF+np9NmzbJZrPpkksuUVJSklq0aKEWLVroxhtvdHRup+KYH/PuuusuLVy4UB06dNCKFStks9lMRwIAuLnzff2ucvmpyPHjx7V+/XqtXbvWXoo2bNigQ4cOOeqfcAnKj1lffvmlOnXqJJvNpu+++05XXHGF6UgAAA/g9JMcVqRatWq68sordeWVVzryaeFDSkpKlJycLEl65JFHKD4AAIer8gHPgDM888wzWrdunWrWrGk/5gcAAEei/MBt7Nu3z/6prgkTJigqKspwIgCAN6L8wG1kZGTo4MGDatGihfr162c6DgDAS1F+4Ba+//57Pffcc5KknJwcBQQ49HA0AADs3Kr8fPnll+rVq5diY2Nls9n07rvvlrvfsiyNHTtWsbGxCgkJUefOnbVx48ZyY4qKipScnKyoqChVr15dN998s3bu3OnCtUBlWZal5ORkWZalO++8U3/5y19MRwIAeDG3Kj9HjhxRy5YtNWfOnArvnzp1qrKysjRnzhytWrVKMTEx6tKliwoLC+1jUlNTtWjRIr322mtasWKFDh8+rJ49e6q0tNRVq4FKeuWVV7Ry5UqFhoZygkwAgNM59Dw/jmSz2bRo0SL17t1b0sm9A7GxsUpNTdWwYcMkndzLEx0drSlTpqhfv37Kz89X7dq19fLLL+uOO+6QJO3atUvx8fFavHixunXrVuG/VVRUpKKiIvv3BQUFio+P5zw/LlBQUKCmTZsqLy9PkydP1ogRI0xHAgB4KKdf2NTVcnNzlZeXp65du9qXBQcHq1OnTlq5cqUkafXq1Tpx4kS5MbGxsWrevLl9TEUyMzMVERFhv8XHxztvRVDOhAkTlJeXp8aNG2vw4MGm4wAAfIDHlJ9TV5CPjo4utzw6Otp+X15enoKCglSzZs0zjqnIiBEjlJ+fb7/t2LHDwelRkZ9//lmzZs2SJM2aNUvBwcFmAwEAfILHfaTmz9d4sizrnNd9OteY4OBgXnhdzLIspaamqqSkRDfddJNuuukm05EAAD7CY/b8xMTESNJpe3D27Nlj3xsUExOj4uJiHTx48Ixj4B7ef/99ffLJJwoKCrLv/QEAwBU8pvw0aNBAMTExWrp0qX1ZcXGxli1bpg4dOkiSWrdurcDAwHJjdu/erQ0bNtjHwLxjx45p0KBBkqTBgwercePGhhMBAHyJW73tdfjwYW3evNn+fW5urtauXavIyEjVq1dPqampmjx5shITE5WYmKjJkycrNDRUd911lyQpIiJCffv21ZAhQ1SrVi1FRkYqLS1NSUlJuuGGG0ytFv5k+vTpys3N1cUXX6yMjAzTcQAAPsatys93332na6+91v79qU//3HfffZo/f76GDh2qY8eOacCAATp48KCuuuoqLVmyRGFhYfbHzJw5UwEBAerTp4+OHTum66+/XvPnz5e/v7/L1wen2759uzIzMyVJ06ZNU40aNQwnAgD4Grc9z49J53ueAFRenz599Oabb+ovf/mLvvjii3MerA4AwPnyuvP8wPN99tlnevPNN+Xn56fs7GyKDwDACMoPXOLEiRNKSUmRJP3jH/9Qy5YtDScCAPgqyg9cYu7cudq4caNq1aql8ePHm44DAPBhlB843Z49ezR69GhJ0qRJkxQZGWk4EQDAl1F+4HTp6enKz8/XFVdcoYceesh0HACAj6P8wKlWrVql559/XpKUk5PDKQcAAMZRfuA0ZWVlSk5OlmVZuueeezjLNgDALVB+4DQvvfSS/vOf/ygsLExTpkwxHQcAAEmUHzhJfn6+hg0bJkkaPXq06tatazgRAAAnUX7gFOPGjdOePXvUtGlT+/l9AABwB5QfONyPP/6onJwcSdLs2bMVFBRkOBEAAP9D+YFDWZallJQUlZSU6JZbblG3bt1MRwIAoBzKDxzqnXfe0aeffqrg4GBlZWWZjgMAwGkoP3CYo0ePavDgwZKkoUOHqmHDhoYTAQBwOsoPHGbq1Knavn274uPjNXz4cNNxAACoEOUHDrFt2zb7uXxmzJih0NBQw4kAAKgY5QcOMXjwYB0/flzXXnutbrvtNtNxAAA4I8oPLtjSpUu1aNEi+fv7KycnRzabzXQkAADOiPKDC1JcXGw/ieHAgQN12WWXGU4EAMDZUX5wQebMmaOff/5ZtWvX1tixY03HAQDgnCg/qLK8vDx74cnMzNRFF11kNA8AAOeD8oMqGz58uAoLC3XllVfqgQceMB0HAIDzQvlBlXz99dd68cUXJUk5OTny8+NHCQDgGXjFQqWVlpYqOTlZknT//ffrqquuMpwIAIDzR/lBpT3//PNavXq1wsPD9cQTT5iOAwBApVB+UCkHDx5Uenq6JGns2LGKjo42nAgAgMqh/KBSxowZo3379unSSy/VwIEDTccBAKDSKD84b+vXr9dTTz0lSZo9e7YCAwMNJwIAoPIoPzgvlmUpOTlZpaWluvXWW3XDDTeYjgQAQJVQfnBe3njjDS1btkzVqlXTjBkzTMcBAKDKKD84pyNHjigtLU3SyRMb1q9f32wgAAAuAOUH55SZmamdO3eqfv36Gjp0qOk4AABcEMoPzmrLli2aNm2aJCkrK0shISGGEwEAcGEoPzirQYMGqbi4WF26dFHv3r1NxwEA4IJRfnBGH330kf71r38pICBAs2fPls1mMx0JAIALRvlBhYqLi5WamipJSklJ0SWXXGI2EAAADkL5QYVmzZql//73v4qOjtaYMWNMxwEAwGEoPzjNrl27NGHCBEnSlClTFB4ebjgRAACOQ/nBaYYNG6bDhw+rXbt2uueee0zHAQDAoSg/KOerr77SK6+8IpvNppycHPn58SMCAPAuvLLBrrS01H6l9r59+6pNmzaGEwEA4HiUH9g999xzWrt2rS666CJNnjzZdBwAAJyC8gNJ0v79+5WRkSFJGj9+vGrXrm04EQAAzkH5gSRp9OjROnDggJKSkvSPf/zDdBwAAJyG8gP98MMPevrppyVJ2dnZCggIMJwIAADnofz4OMuylJycrLKyMvXp00edO3c2HQkAAKei/Pi4hQsXavny5QoNDdX06dNNxwEAwOkoPz7s8OHDevzxxyVJ6enpio+PN5wIAADno/z4sEmTJmnXrl1q2LChhgwZYjoOAAAuQfnxUZs2bdKMGTMknbyIabVq1QwnAgDANSg/Pio1NVUnTpxQ9+7d1bNnT9NxAABwGcqPD/rggw+0ePFiBQYGatasWbLZbKYjAQDgMpQfH3P8+HGlpqZKkgYNGqQmTZqYDQQAgItRfnzMzJkztWXLFtWtW1cjR440HQcAAJej/PiQnTt3auLEiZKkadOmKSwszHAiAABcj/LjQx5//HEdPXpU11xzje666y7TcQAAMILy4yOWLVum1157TX5+fsrJyeEgZwCAz6L8+ICSkhIlJydLkvr166dWrVqZDQQAgEGUHx/w9NNPa/369YqMjNSECRNMxwEAwCjKj5fbt2+fRo0aJUmaOHGiatWqZTgRAABmUX68XEZGhg4dOqRWrVrpkUceMR0HAADjKD9ebPXq1XruueckSTk5OfL39zecCAAA8yg/XqqsrEzJycmyLEt///vfdc0115iOBACAW6D8eKlXXnlFX3/9tWrUqKGpU6eajgMAgNug/HihgoICDR06VJI0cuRIxcbGGk4EAID7oPx4oQkTJuj3339XYmKi/SKmAADgJMqPl/n55581a9YsSdLs2bMVHBxsNhAAAG7Go8rP2LFjZbPZyt1iYmLs91uWpbFjxyo2NlYhISHq3LmzNm7caDCxa1mWpccee0wlJSXq1auXunfvbjoSAABux6PKjyRddtll2r17t/22fv16+31Tp05VVlaW5syZo1WrVikmJkZdunRRYWGhwcSu895772nJkiUKCgrSzJkzTccBAMAtBZgOUFkBAQHl9vacYlmWZs2apYyMDN16662SpBdffFHR0dFasGCB+vXrd8bnLCoqUlFRkf37goICxwd3smPHjmnQoEGSpLS0NDVq1MhwIgAA3JPH7fnZtGmTYmNj1aBBA915553aunWrJCk3N1d5eXnq2rWrfWxwcLA6deqklStXnvU5MzMzFRERYb/Fx8c7dR2cYdq0adq2bZvi4uKUnp5uOg4AAG7Lo8rPVVddpZdeekmffPKJnnvuOeXl5alDhw7av3+/8vLyJEnR0dHlHhMdHW2/70xGjBih/Px8+23Hjh1OWwdn+PXXX5WZmSlJmj59uqpXr244EQAA7suj3vb64wG8SUlJat++vRo1aqQXX3xR7dq1kyTZbLZyj7Es67RlfxYcHOzRn4pKS0vT8ePH1blzZ/Xp08d0HAAA3JpH7fn5s+rVqyspKUmbNm2yHwf05708e/bsOW1vkDf59NNP9dZbb8nPz0/Z2dnnLHoAAPg6jy4/RUVF+umnn1S3bl01aNBAMTExWrp0qf3+4uJiLVu2TB06dDCY0nlOnDihlJQUSdKAAQOUlJRkOBEAAO7Po972SktLU69evVSvXj3t2bNHEydOVEFBge677z7ZbDalpqZq8uTJSkxMVGJioiZPnqzQ0FDdddddpqM7xZNPPqkff/xRUVFRGj9+vOk4AAB4BI8qPzt37tTf/vY37du3T7Vr11a7du30zTffKCEhQZI0dOhQHTt2TAMGDNDBgwd11VVXacmSJQoLCzOc3PF+//13jRkzRpI0efJk1axZ03AiAAA8g82yLMt0CHdTUFCgiIgI5efnKzw83HScCj344IN64YUX1Lp1a/3nP/+Rv7+/6UgAABh1vq/fHn3Mj6/69ttv9cILL0iScnJyKD4AAFQC5cfDlJWVaeDAgZKke++9V+3btzecCAAAz0L58TDz58/XqlWrFBYWpilTppiOAwCAx6H8eJBDhw5p+PDhkqQxY8ZUeI0zAABwdpQfDzJ27Fjt3btXTZs2VXJysuk4AAB4JMqPh9iwYYPmzJkjScrOzlZQUJDhRAAAeCbKjwewLEuPPfaYSktL1bt373JXrgcAAJVD+fEAb7/9tj777DNVq1ZNWVlZpuMAAODRKD9u7ujRoxoyZIikk2ewbtCggeFEAAB4NsqPm3viiSe0fft21atXT8OGDTMdBwAAj0f5cWO5ubmaOnWqJGnGjBkKDQ01nAgAAM9H+XFjgwcPVlFRka677jr93//9n+k4AAB4BcqPm1qyZIneffdd+fv7Kzs7WzabzXQkAAC8AuXHDRUXFyslJUWSlJycrMsuu8xwIgAAvAflxw3l5OTol19+UZ06dTR27FjTcQAA8CqUHzeze/dujRs3TpKUmZmpiIgIw4kAAPAulB83M3z4cBUWFqpt27a6//77TccBAMDrUH7cyMqVK/XSSy9JOvnWl58f/3sAAHA0Xl3dRGlpqf1K7Q8++KDatm1rOBEAAN6J8uMmnn/+eX3//feKiIhQZmam6TgAAHgtyo8bOHjwoNLT0yVJ48aNU506dQwnAgDAe1F+3MDo0aO1b98+XXrppRowYIDpOAAAeDXKj2Hr1q3TU089JenkQc6BgYGGEwEA4N0oPwZZlqWUlBSVlZXptttu03XXXWc6EgAAXo/yY9Abb7yhZcuWKSQkRNOnTzcdBwAAn0D5MeTIkSNKS0uTJI0YMUIJCQmGEwEA4BsoP4ZMnjxZO3fuVP369e0lCAAAOB/lx4DNmzfb3+aaOXOmQkJCDCcCAMB3UH4MGDRokIqLi9W1a1fdcsstpuMAAOBTKD8utnjxYn3wwQcKCAjQ7NmzZbPZTEcCAMCnUH5cqKioSKmpqZKk1NRUNWvWzGwgAAB8EOXHhWbNmqVNmzYpJiZGo0aNMh0HAACfRPlxkd9++00TJkyQJE2ZMkXh4eGGEwEA4JsoPy4ydOhQHTlyRO3bt9fdd99tOg4AAD6L8uMCK1as0IIFC2Sz2TRnzhz5+THtAACYwquwk5WWlio5OVmS9PDDD+uKK64wnAgAAN9G+XGyZ599VmvXrtVFF12kSZMmmY4DAIDPo/w40f79+zVy5EhJ0oQJExQVFWU4EQAAoPw40ciRI3XgwAElJSWpf//+puMAAABRfpxmzZo1euaZZyRJOTk5CggIMJwIAABIlB+nsCxLycnJsixLd955pzp16mQ6EgAA+P8oP06wYMECffXVVwoNDdW0adNMxwEAAH9A+XGwwsJCPf7445KkjIwMxcXFGU4EAAD+iPLjYBMnTtTu3bvVqFEjDRkyxHQcAADwJ5QfB/rll180c+ZMSScvYhocHGw4EQAA+DPKj4NYlqXU1FSdOHFCPXr0UM+ePU1HAgAAFaD8OMgHH3ygjz/+WEFBQZo1a5bpOAAA4AwoPw5w/PhxpaamSpIGDx6sxMREs4EAAMAZUX4cYMaMGdq6datiY2OVkZFhOg4AADgLys8F2rFjhyZPnixJmjZtmmrUqGE4EQAAOBvKzwVKS0vT0aNH1bFjR/3tb38zHQcAAJwD5ecCfPHFF3rjjTfk5+ennJwc2Ww205EAAMA5UH6qqKSkRCkpKZKk/v37q2XLloYTAQCA80H5qaK5c+dq/fr1ioyM1IQJE0zHAQAA54nyUwV79+7V6NGjJUmTJk1SZGSk4UQAAOB8UX6qID09XYcOHdLll1+uhx9+2HQcAABQCZSfSvruu+80b948SVJOTo78/f0NJwIAAJVB+amEsrIyJScny7Is3X333br66qtNRwIAAJVE+amEl19+Wd98841q1KihqVOnmo4DAACqgPJznvLz8zVs2DBJ0ujRo1W3bl3DiQAAQFVQfs7T+PHj9fvvv6tJkyZ67LHHTMcBAABVRPk5D5Zl6cSJE7LZbJo9e7aCgoJMRwIAAFVksyzLMh3C3RQUFCgiIkL5+fkKDw+3L9+8ebMaN25sMBkAADiTM71+/xl7fiqB4gMAgOej/AAAAJ9C+QEAAD6F8gMAAHyK15afp556Sg0aNFC1atXUunVrLV++3HQkAADgBryy/Lz++utKTU1VRkaG1qxZo44dO6p79+7avn276WgAAMAwr/yo+1VXXaUrrrhCc+fOtS+75JJL1Lt3b2VmZp42vqioSEVFRfbvCwoKFB8ff86PygEAAPfhsx91Ly4u1urVq9W1a9dyy7t27aqVK1dW+JjMzExFRETYb/Hx8a6ICgAADPC68rNv3z6VlpYqOjq63PLo6Gjl5eVV+JgRI0YoPz/fftuxY4crogIAAAMCTAdwFpvNVu57y7JOW3ZKcHCwgoODXRELAAAY5nV7fqKiouTv73/aXp49e/actjcIAAD4Hq8rP0FBQWrdurWWLl1abvnSpUvVoUMHQ6kAAIC78Mq3vQYPHqx77rlHbdq0Ufv27fXss89q+/bt6t+/v+loAADAMK8sP3fccYf279+v8ePHa/fu3WrevLkWL16shIQE09EAAIBhXnmenwuVn5+viy66SDt27OA8PwAAeIhT5+k7dOiQIiIizjjOK/f8XKj9+/dLEuf7AQDAAxUWFlJ+KisyMlKStH379rNOnrc71aB9fQ8Y83AS83AS8/A/zMVJzMNJ7jAPlmWpsLBQsbGxZx1H+amAn9/JD8FFRET49A/yKeHh4cyDmIdTmIeTmIf/YS5OYh5OMj0P57PTwus+6g4AAHA2lB8AAOBTKD8VCA4O1pgxY3z+khfMw0nMw0nMw0nMw/8wFycxDyd50jzwUXcAAOBT2PMDAAB8CuUHAAD4FMoPAADwKZQfAADgU3y2/Dz11FNq0KCBqlWrptatW2v58uVnHb9s2TK1bt1a1apVU8OGDfX000+7KKlzVWYe3nnnHXXp0kW1a9dWeHi42rdvr08++cSFaZ2nsj8Pp3z11VcKCAhQq1atnBvQRSo7D0VFRcrIyFBCQoKCg4PVqFEjPf/88y5K6zyVnYdXX31VLVu2VGhoqOrWrasHHnjAfpkcT/Xll1+qV69eio2Nlc1m07vvvnvOx3jjdrKy8+Ct28mq/Dyc4o7bSZ8sP6+//rpSU1OVkZGhNWvWqGPHjurevbu2b99e4fjc3Fz16NFDHTt21Jo1a5Senq6UlBS9/fbbLk7uWJWdhy+//FJdunTR4sWLtXr1al177bXq1auX1qxZ4+LkjlXZeTglPz9f9957r66//noXJXWuqsxDnz599Omnn2revHn65ZdftHDhQjVr1syFqR2vsvOwYsUK3Xvvverbt682btyoN998U6tWrdJDDz3k4uSOdeTIEbVs2VJz5sw5r/Heup2s7Dx463aysvNwittuJy0f1LZtW6t///7lljVr1swaPnx4heOHDh1qNWvWrNyyfv36We3atXNaRleo7DxU5NJLL7XGjRvn6GguVdV5uOOOO6yRI0daY8aMsVq2bOnEhK5R2Xn46KOPrIiICGv//v2uiOcylZ2HadOmWQ0bNiy3LDs724qLi3NaRleTZC1atOisY7x1O/lH5zMPFfGG7eQfVWYe3HU76XN7foqLi7V69Wp17dq13PKuXbtq5cqVFT7m66+/Pm18t27d9N133+nEiRNOy+pMVZmHPysrK1NhYaH9QrCeqKrz8MILL2jLli0aM2aMsyO6RFXm4f3331ebNm00depUXXzxxWrSpInS0tJ07NgxV0R2iqrMQ4cOHbRz504tXrxYlmXp999/11tvvaWbbrrJFZHdhjduJx3BG7aTVeXO20mfu7Dpvn37VFpaqujo6HLLo6OjlZeXV+Fj8vLyKhxfUlKiffv2qW7duk7L6yxVmYc/mzFjho4cOaI+ffo4I6JLVGUeNm3apOHDh2v58uUKCPCOX6GqzMPWrVu1YsUKVatWTYsWLdK+ffs0YMAAHThwwGOP+6nKPHTo0EGvvvqq7rjjDh0/flwlJSW6+eablZOT44rIbsMbt5OO4A3byapw9+2kz+35OcVms5X73rKs05ada3xFyz1NZefhlIULF2rs2LF6/fXXVadOHWfFc5nznYfS0lLdddddGjdunJo0aeKqeC5TmZ+HsrIy2Ww2vfrqq2rbtq169OihrKwszZ8/36P3/kiVm4cff/xRKSkpGj16tFavXq2PP/5Yubm56t+/vyuiuhVv3U5WlbdtJ8+XJ2wn3a+OOVlUVJT8/f1P+ytuz549p/3VckpMTEyF4wMCAlSrVi2nZXWmqszDKa+//rr69u2rN998UzfccIMzYzpdZeehsLBQ3333ndasWaOBAwdKOlkCLMtSQECAlixZouuuu84l2R2pKj8PdevW1cUXX6yIiAj7sksuuUSWZWnnzp1KTEx0amZnqMo8ZGZm6uqrr9bjjz8uSWrRooWqV6+ujh07auLEiT6zx8Mbt5MXwpu2k5XlCdtJn9vzExQUpNatW2vp0qXlli9dulQdOnSo8DHt27c/bfySJUvUpk0bBQYGOi2rM1VlHqSTf8ncf//9WrBggVcc01DZeQgPD9f69eu1du1a+61///5q2rSp1q5dq6uuuspV0R2qKj8PV199tXbt2qXDhw/bl/33v/+Vn5+f4uLinJrXWaoyD0ePHpWfX/lNqb+/v6T/7fnwBd64nawqb9tOVpZHbCfNHGdt1muvvWYFBgZa8+bNs3788UcrNTXVql69urVt2zbLsixr+PDh1j333GMfv3XrVis0NNQaNGiQ9eOPP1rz5s2zAgMDrbfeesvUKjhEZedhwYIFVkBAgPXkk09au3fvtt8OHTpkahUcorLz8Gfu9imGqqrsPBQWFlpxcXHWbbfdZm3cuNFatmyZlZiYaD300EOmVsEhKjsPL7zwghUQEGA99dRT1pYtW6wVK1ZYbdq0sdq2bWtqFRyisLDQWrNmjbVmzRpLkpWVlWWtWbPG+vXXXy3L8p3tZGXnwVu3k5Wdhz9zt+2kT5Yfy7KsJ5980kpISLCCgoKsK664wlq2bJn9vvvuu8/q1KlTufFffPGFdfnll1tBQUFW/fr1rblz57o4sXNUZh46depkSTrtdt9997k+uINV9ufhj9ztl/pCVHYefvrpJ+uGG26wQkJCrLi4OGvw4MHW0aNHXZza8So7D9nZ2dall15qhYSEWHXr1rX+/ve/Wzt37nRxasf6/PPPz/r77ivbycrOg7duJ6vy8/BH7radtFmWD+2XBQAAPs/njvkBAAC+jfIDAAB8CuUHAAD4FMoPAADwKZQfAADgUyg/AADAp1B+AACAT6H8AAAAn0L5AQAAPoXyAwAAfArlB4DX279/v+rUqaNt27Zd0PPcdtttysrKckwoAMZwbS8AXi8tLU0HDx7UvHnzLuh51q1bp2uvvVa5ubkKDw93UDoArsaeHwBe7dixY5o3b54eeuihC36uFi1aqH79+nr11VcdkAyAKZQfAB5l4cKFqlatmn777Tf7soceekgtWrRQfn7+aeM/+ugjBQQEqH379uWWd+7cWcnJyUpNTVXNmjUVHR2tZ599VkeOHNEDDzygsLAwNWrUSB999FG5x918881auHChc1YOgEtQfgB4lDvvvFNNmzZVZmamJGncuHH65JNP9NFHHykiIuK08V9++aXatGlT4XO9+OKLioqK0rfffqvk5GT94x//0O23364OHTro+++/V7du3XTPPffo6NGj9se0bdtW3377rYqKipyzggCcjmN+AHicDz74QLfddptGjx6t6dOna/ny5brssssqHNu7d2/VqlXrtON9OnfurNLSUi1fvlySVFpaqoiICN1666166aWXJEl5eXmqW7euvv76a7Vr107SyeN+WrZsqW3btikhIcGJawnAWQJMBwCAyurZs6cuvfRSjRs3TkuWLDlj8ZFOHvNTrVq1Cu9r0aKF/Wt/f3/VqlVLSUlJ9mXR0dGSpD179tiXhYSESFK5vUEAPAtvewHwOJ988ol+/vlnlZaW2gvKmURFRengwYMV3hcYGFjue5vNVm6ZzWaTJJWVldmXHThwQJJUu3btKmUHYB7lB4BH+f7773X77bfrmWeeUbdu3TRq1Kizjr/88sv1448/Ouzf37Bhg+Li4hQVFeWw5wTgWpQfAB5j27ZtuummmzR8+HDdc889Gj9+vN5++22tXr36jI/p1q2bNm7ceMa9P5W1fPlyde3a1SHPBcAMyg8Aj3DgwAF1795dN998s9LT0yVJrVu3Vq9evZSRkXHGxyUlJalNmzZ64403LjjD8ePHtWjRIj388MMX/FwAzOHTXgC83uLFi5WWlqYNGzbIz6/qf/M9+eSTeu+997RkyRIHpgPganzaC4DX69GjhzZt2qTffvtN8fHxVX6ewMBA5eTkODAZABPY8wMAAHwKx/wAAACfQvkBAAA+hfIDAAB8CuUHAAD4FMoPAADwKZQfAADgUyg/AADAp1B+AACAT6H8AAAAn/L/AOl/imDlZXLtAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Visualize the results\n", "plt.figure()\n", "plt.plot(x, phi, '-k')\n", "plt.xlabel('$x$ (m)')\n", "plt.ylabel(r'$ phi$')\n", "ax = plt.gca()\n", "ax.set_xlim(0, L)\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "7559eabc", "metadata": {}, "source": [ "## Exercise\n", "\n", "1. Calculate the source term locally and redo the simulation.\n", "\n", "2. Try the Crank-Nicolson scheme to solve the problem by yourself (Hint: change the renew function).\n", "\n", "3. Compare the speed of calculations by TDMA and Jacobi iterations." ] } ], "metadata": { "kernelspec": { "display_name": "cfd-python", "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.10.14" } }, "nbformat": 4, "nbformat_minor": 5 }