{ "cells": [ { "cell_type": "markdown", "id": "555197f4", "metadata": {}, "source": [ "# Steady-state advection-diffusion\n", "\n", "Advection-diffusion problems study the *transport of substances* (e.g., chemicals, heat, or particles) influenced by both directional flow and random motion, governed by the advection-diffusion equation. At steady state, the interplay between convection and diffusion shapes the solution based on the velocity field and boundary conditions, with analytic solutions for simple geometries and numerical methods used for more complex cases." ] }, { "cell_type": "markdown", "id": "c6cc27fa", "metadata": {}, "source": [ "## Problem setup\n", "\n", "Consider a one-dimensional steady-state advection-diffusion problem where a field variable $\\phi$ is transported through the advection-diffusion process from $x = 0$ to $x = L$ in a one-dimensional domain. The fluid density is $\\rho = 1.0$ kg/m$^3$, $L = 1.0$ m, and the diffusion coefficient is $\\Gamma = 0.1$ kg/(m$\\cdot$s). As shown in the figure below, determine the distribution of $\\phi$ on a grid discretized into 5 nodes when the velocity $u = 0.1$ m/s.\n", "\n", "![Adiabatic rod](images/advection_diffusion1d.jpg)\n", "\n", "The mathematical model for one-dimensional steady-state advection-diffusion problem is\n", "\n", "$$\\frac{\\mathrm{d}}{\\mathrm{d}x} \\left( \\rho u \\phi \\right) = \\frac{\\mathrm{d}}{\\mathrm{d}x} \\left( k\\frac{\\mathrm{d} \\phi}{\\mathrm{d}x} \\right).$$\n", "\n", "The boundary conditions are\n", "\n", "$$\\phi|_{x=0} = \\phi_A = 1, \\quad \\phi|_{x=L} = \\phi_B = 0.$$\n", "\n", "The analytical solution for the temperature distribution at steady state is\n", "\n", "$$\\frac{\\phi - \\phi_A}{\\phi_B - \\phi_A} = \\frac{\\exp \\left( \\frac{\\rho u x}{\\Gamma} \\right) - 1}{\\exp \\left( \\frac{\\rho u L}{\\Gamma} \\right) - 1}.$$" ] }, { "cell_type": "markdown", "id": "8630aa0e", "metadata": {}, "source": [ "## Solve problem\n", "\n", "### Define grid\n", "\n", "Divide the rod evenly into 5 control volumes, as a result, the length of each control volume becomes $\\delta x=0.2$ m.\n", "\n", "![Grid](images/advection_diffusion1d_grid.jpg)" ] }, { "cell_type": "markdown", "id": "84b40255", "metadata": {}, "source": [ "### Discrete advection-diffusion equation\n", "\n", "Let $A_w = A_e = 1$, then $F = \\rho u$ and $D = \\frac{\\Gamma}{\\delta x}, with $F_e = F_w = F$ and $D_e = D_w = D$ holding for all control volumes. The discretized equations for the internal nodes 2, 3, and 4 are identical, while the boundary nodes 1 and 5 require special treatment.\n", "\n", "**Internal Nodes 2, 3, 4**\n", "\n", "The discrete equation satisfied by the internal nodes is\n", "\n", "$$a_P \\phi_P = a_W \\phi_W + a_E \\phi_E,$$\n", "\n", "where\n", "\n", "$$a_W = D_w + \\frac{F_w}{2} = D + \\frac{F}{2}, \\quad a_E = D_e - \\frac{F_e}{2} = D - \\frac{F}{2}, \\quad a_P = a_W + a_E + (F_e - F_w) = a_W + a_E.$$\n", "\n", "**Boundary Nodes 1**\n", "\n", "The discrete equation for the boundary node 1 is\n", "\n", "$$a_P T_P = a_W T_W + a_E T_E + S_u,$$\n", "\n", "where\n", "\n", "$$a_W = 0, \\quad a_E = D_e - \\frac{F_e}{2} = D - \\frac{F}{2}, \\quad a_P = a_W + a_E + (F_e - F_w) - S_p = a_W + a_E - S_P,$$\n", "\n", "$$S_P = -(2D_A + F_w) = -(2D + F), \\quad S_u = (2D_A + F_A) \\phi_A = (2D + F) \\phi_A$$\n", "\n", "**Boundary Nodes 5**\n", "\n", "The discrete equation for the boundary node 1 is\n", "\n", "$$a_P T_P = a_W T_W + a_E T_E + S_u,$$\n", "\n", "where\n", "\n", "$$a_W = D_w + \\frac{F_w}{2} = D + \\frac{F}{2}, \\quad a_E = 0, \\quad a_P = a_W + a_E + (F_e - F_w) - S_p = a_W + a_E - S_P,$$\n", "\n", "$$S_P = -(2D_B + F_e) = -(2D - F), \\quad S_u = (2D_B - F_B) \\phi_A = (2D - F) \\phi_B$$" ] }, { "cell_type": "markdown", "id": "e9abfb83", "metadata": {}, "source": [ "### Solve algebraic equations" ] }, { "cell_type": "code", "execution_count": 20, "id": "1b8669f8", "metadata": {}, "outputs": [], "source": [ "# Step 0: import required libraries\n", "import numpy as np # for array operation\n", "from matplotlib import pyplot as plt # for plotting figures\n", "import time # for time measurement" ] }, { "cell_type": "code", "execution_count": 21, "id": "ab0f1ec9", "metadata": {}, "outputs": [], "source": [ "# Step 1: parameter declarations\n", "nx = 5 # number of spatial grid points\n", "L = 1.0 # length of the domain\n", "dx = L / nx # spatial grid size\n", "x = np.linspace(0.5*dx, L-0.5*dx, nx) # spatial grid points\n", "rho = 1.0 # fluid density\n", "u = 0.1 # fluid velocity\n", "F = rho*u # advection flux\n", "Gamma = 0.1 # diffusion coefficient\n", "D = Gamma/dx # diffusion conductance per unit area\n", "phiA = 1 # phi value at x = 0\n", "phiB = 0 # phi value at x = L" ] }, { "cell_type": "code", "execution_count": 22, "id": "a354d166", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkYAAAG4CAYAAACgrSiYAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAANWBJREFUeJzt3XtYFfXe///X4qwQoKgcFBFNDTUtIQ3MnZVimpSVaXd3KqVdm06mfK00y9PVvdlZlllqJ9S6LzMrzbrbpNJJ8VAmQrWTToqhxmFDCrIsFJjfHw7r5xJQQA6L5fNxXXNdrM98ZtZ7HGpefD6zZlkMwzAEAAAAubR0AQAAAI6CYAQAAGAiGAEAAJgIRgAAACaCEQAAgIlgBAAAYCIYAQAAmAhGAAAAJoIRAACAiWAEAABgcqhgtG3bNsXFxSkkJEQWi0UbN248Z/8NGzZoxIgR6tixo3x9fRUdHa3Nmzc3T7EAAMDpOFQwslqtGjBggF5++eU69d+2bZtGjBihlJQUpaen67rrrlNcXJwyMjKauFIAAOCMLI76JbIWi0UffPCBxo4dW6/t+vbtqwkTJmju3LlNUxgAAHBabi1dQGOqrKzU8ePH1b59+1r7lJWVqayszG6bP/74QwEBAbJYLM1RJgAAuECGYej48eMKCQmRi0vjTYA5VTBavHixrFarxo8fX2ufpKQkLViwoBmrAgAATeXQoUPq0qVLo+3PaabS1q5dq6lTp+rDDz/U8OHDa+139ohRcXGxunbtqkOHDsnX1/dCywYAAM2gpKREoaGhOnbsmPz8/Bptv04xYrRu3TpNmTJF77333jlDkSR5enrK09OzWruvry/BCACAVqaxb4NxqE+lNcTatWsVHx+vt99+WzfddFNLlwMAAFoxhxoxKi0t1a+//mp7nZ2drczMTLVv315du3bV7NmzdeTIEb311luSToeiSZMm6cUXX9TVV1+tvLw8SVKbNm0adVgNAABcHBxqxGjPnj268sordeWVV0qSEhMTdeWVV9o+ep+bm6ucnBxb/1dffVXl5eV68MEHFRwcbFseeeSRFqkfAAC0bg5783VzKSkpkZ+fn4qLi7nHCACAVqKprt8ONWIEAADQkghGAAAAJoIRAACAiWAEAABgIhgBAACYCEYAAAAmghEAAICJYAQAAGAiGAEAAJgIRgAAACaCEQAAgIlgBAAAYCIYAQAAmAhGAAAAJoIRAACAiWAEAABgIhgBAACYCEYAAAAmghEAAICJYAQAAGAiGAEAAJgIRgAAACaCEQAAgIlgBAAAYCIYAQAAmAhGAAAAJoIRAACAiWAEAABgIhgBAACYCEYAAAAmghEAAICJYAQAAGAiGAEAAJgIRgAAACaCEQAAgIlgBAAAYCIYAQAAmAhGAAAAJoIRAACAiWAEAABgIhgBAACYCEYAAAAmghEAAICJYAQAAGAiGAEAAJgIRgAAACaCEQAAgIlgBAAAYCIYAQAAmAhGAAAAJoIRAACAyaGC0bZt2xQXF6eQkBBZLBZt3LjxvNts3bpVkZGR8vLyUvfu3fXKK680faEAAMApOVQwslqtGjBggF5++eU69c/Oztbo0aM1dOhQZWRk6IknntC0adO0fv36Jq4UAAA4I7eWLuBMo0aN0qhRo+rc/5VXXlHXrl21ZMkSSVJERIT27Nmj5557TrfffnsTVQkAAJyVQ40Y1deuXbsUGxtr1zZy5Ejt2bNHp06dqnGbsrIylZSU2C0AAABSKw9GeXl5CgwMtGsLDAxUeXm5CgsLa9wmKSlJfn5+tiU0NLQ5SgUAAK1Aqw5GkmSxWOxeG4ZRY3uV2bNnq7i42LYcOnSoyWsEAACtg0PdY1RfQUFBysvLs2srKCiQm5ubAgICatzG09NTnp6ezVEeAABoZVr1iFF0dLRSU1Pt2rZs2aKoqCi5u7u3UFUAAKC1cqhgVFpaqszMTGVmZko6/XH8zMxM5eTkSDo9DTZp0iRb/4SEBP32229KTExUVlaWVq5cqeTkZM2cObMlygcAAK2cQ02l7dmzR9ddd53tdWJioiRp8uTJWr16tXJzc20hSZLCw8OVkpKiGTNmaNmyZQoJCdHSpUv5qD4AAGgQi1F1t/JFqqSkRH5+fiouLpavr29LlwMAAOqgqa7fDjWVBgAA0JIIRgAAACaCEQAAgIlgBAAAYCIYAQAAmAhGAAAAJoIRAACAiWDkBKxWqywWiywWi6xWa0uXAwBAq0UwAgAAMBGMAAAATA71XWmouzOnzGr7WZK8vb2brSYAAFo7glEr5ePjU2N7YGCg3euL/KvwAACoF6bSAAAATIwYtVKlpaW2n61Wq22kKD8/n+kzAAAaiGDUStUWfry9vQlGAAA0EFNpAAAAJoIRAACAiak0J+Dt7c2nzwAAaASMGAEAAJgIRgAAACaCEQAAgIlgBAAAYCIYAQAAmAhGAAAAJoIRAACAiWAEAABgIhgBAACYCEYAAAAmghEAAICJYAQAAGAiGAEAAJgIRgAAACaCEQAAgIlgBAAAYCIYAQAAmAhGAAAAJoIRAACAiWAEAABgIhgBAACYCEYAAAAmghEAAICJYAQAAGAiGAEAAJgIRgAAACaCEQAAgIlgBAAAYCIYAQAAmAhGAAAAJoIRAACAiWAEAABgcrhgtHz5coWHh8vLy0uRkZFKS0s7Z/81a9ZowIABatu2rYKDg3XPPfeoqKiomaoF7FmtVlksFlksFlmt1pYuBwBQTw4VjNatW6fp06drzpw5ysjI0NChQzVq1Cjl5OTU2H/79u2aNGmSpkyZoh9++EHvvfeevvnmG02dOrWZKwcAAM7AoYLR888/rylTpmjq1KmKiIjQkiVLFBoaqhUrVtTY/6uvvlK3bt00bdo0hYeH65prrtHf//537dmzp5krBwAAzsBhgtHJkyeVnp6u2NhYu/bY2Fjt3Lmzxm1iYmJ0+PBhpaSkyDAM5efn6/3339dNN91U6/uUlZWppKTEbgEuhNVqtVvO1w4AcFwOE4wKCwtVUVGhwMBAu/bAwEDl5eXVuE1MTIzWrFmjCRMmyMPDQ0FBQfL399dLL71U6/skJSXJz8/PtoSGhjbqceDi4+PjY1vO/P0NDAy0WwcAcHwOE4yqWCwWu9eGYVRrq7Jv3z5NmzZNc+fOVXp6ujZt2qTs7GwlJCTUuv/Zs2eruLjYthw6dKhR6wcAAK2XW0sXUKVDhw5ydXWtNjpUUFBQbRSpSlJSkoYMGaJHH31UktS/f395e3tr6NChevrppxUcHFxtG09PT3l6ejb+AeCiVVpaavvZarXafl/z8/Pl7e3dUmUBABrAYUaMPDw8FBkZqdTUVLv21NRUxcTE1LjNiRMn5OJifwiurq6STo80Ac3B29vbbjlfOwDAcTlMMJKkxMREvfHGG1q5cqWysrI0Y8YM5eTk2KbGZs+erUmTJtn6x8XFacOGDVqxYoUOHDigHTt2aNq0aRo0aJBCQkJa6jAAAEAr5TBTaZI0YcIEFRUVaeHChcrNzVW/fv2UkpKisLAwSVJubq7dM43i4+N1/Phxvfzyy/p//+//yd/fX9dff72eeeaZljoEAADQilmMi3zOqaSkRH5+fiouLpavr29LlwMAAOqgqa7fDjWVBgAA0JIIRgAAACaCEQAAgIlgBAAAYCIYAQAAmAhGAAAAJoIRAACAiWAEAABgIhgBAACYCEYAAAAmghEAAICJYAQAAGAiGAEAAJgIRgAAACaCEQAAgIlgBAAAYCIYAQAAmAhGAAAAJoIRAACAiWAEAABgIhgBAACYCEYAAAAmghEAAICJYAQAAGAiGAEAAJgIRgAAACaCEQAAgIlgBAAAYCIYAQAAmAhGAAAAJoIRAACAiWAEAABgIhgBAACYCEYAAAAmghEAAICJYAQAAGAiGAEAAJgIRgAAACaCEQAAgIlgBAAAYCIYAQAAmAhGAAAAJoIRAACAiWAEwClZrVZZLBZZLBZZrdaWLgdAK0EwAgAAMBGMAAAATG4tXQAANJYzp8xq+1mSvL29m60mAK0LwQiA0/Dx8amxPTAw0O61YRjNUQ6AVoipNAAAABMjRgCcRmlpqe1nq9VqGynKz89n+gxAnRCMADiN2sKPt7c3wQhAnTjcVNry5csVHh4uLy8vRUZGKi0t7Zz9y8rKNGfOHIWFhcnT01M9evTQypUrm6laAADgTC5oxGjv3r1KS0uTh4eHhgwZov79+19QMevWrdP06dO1fPlyDRkyRK+++qpGjRqlffv2qWvXrjVuM378eOXn5ys5OVmXXnqpCgoKVF5efkF1AACAi5PFaODHM5YsWaLExET5+/vLzc1NhYWF6tu3r1avXq3IyMgGFTN48GANHDhQK1assLVFRERo7NixSkpKqtZ/06ZNuvPOO3XgwAG1b9++Qe9ZUlIiPz8/FRcXy9fXt0H7AAAAzauprt/1mkpbuXKl9u7dq7KyMv3jH//QP//5TxUVFamgoEC//fabbrnlFg0bNkzbt2+vdyEnT55Uenq6YmNj7dpjY2O1c+fOGrf56KOPFBUVpUWLFqlz587q1auXZs6cqT///LPW9ykrK1NJSYndAgAAINVzKu3ZZ5/Vr7/+KkmqrKzUN998oxdeeEEDBw7UFVdcoaefflqdO3fWzJkz9dVXX9WrkMLCQlVUVFR73khgYKDy8vJq3ObAgQPavn27vLy89MEHH6iwsFAPPPCA/vjjj1rvM0pKStKCBQvqVRsAALg41GvEKCsrS8ePH9fOnTvl7u4uFxcXvfvuu7rpppsUEBCgsLAwvffee8rIyND//d//KTs7u94FWSwWu9eGYVRrq1JZWSmLxaI1a9Zo0KBBGj16tJ5//nmtXr261lGj2bNnq7i42LYcOnSo3jUCAADnVO9PpXl5eemqq67SkCFDNGDAAH311Vc6fvy4vvvuOyUlJalXr146deqU4uPj1aNHjzrP+3Xo0EGurq7VRocKCgqqjSJVCQ4OVufOneXn52dri4iIkGEYOnz4cI3beHp6ytfX124BAACQLuDj+osXL9aiRYs0depU7d27V7169VJcXJx8fHwUEhKioqIi5eTk6N13363T/jw8PBQZGanU1FS79tTUVMXExNS4zZAhQ/T777/bPdTt559/louLi7p06dLQQwMAABepBgejK664Qunp6frtt9909dVXy8vLS/7+/nrppZf0zDPPSJK6dOmiG2+8sc77TExM1BtvvKGVK1cqKytLM2bMUE5OjhISEiSdngabNGmSrf9dd92lgIAA3XPPPdq3b5+2bdumRx99VPfee6/atGnT0EMDAAAXqQt6jlGPHj2Umpqq/Px8ffXVVzp58qSuvvpqhYaGNmh/EyZMUFFRkRYuXKjc3Fz169dPKSkpCgsLkyTl5uYqJyfH1t/Hx0epqal6+OGHFRUVpYCAAI0fP15PP/30hRwWAAC4SDX4OUbOgucYAQDQ+jjEc4wAAACcGcEIAADARDACAAAwEYwAAABMBCMAAAATwQgAAMBEMAIAADARjAAAAEwEIwAAABPBCAAAwEQwAgAAMBGMAAAATAQjAAAAE8EIAADARDACAAAwEYwAAABMBCMAAAATwQgAAMBEMAIAADARjAAAAEwEIwAAABPBCAAAwEQwAgAAMBGMAAAATAQjAAAAE8EIAADARDACAAAwEYwAAABMBCMAAAATwQgAAMBEMAIAADARjAAAAEwEIwAAABPBCAAAwEQwAgAAMBGMAAAATAQjAAAAE8EIAADARDACAAAwEYwAAABMBCMAAAATwQgAAMBEMAIAADARjAAAAEwEIwAAABPBCAAAwEQwAgA0OavVKovFIovFIqvV2tLlALUiGAEAAJgIRgAAACa3li4AAOCczpwyq+1nSfL29m62moDzcbgRo+XLlys8PFxeXl6KjIxUWlpanbbbsWOH3NzcdMUVVzRtgQCAOvHx8bEtgYGBtvbAwEC7dYAjcahgtG7dOk2fPl1z5sxRRkaGhg4dqlGjRiknJ+ec2xUXF2vSpEm64YYbmqlSAADgjCyGYRgtXUSVwYMHa+DAgVqxYoWtLSIiQmPHjlVSUlKt2915553q2bOnXF1dtXHjRmVmZtb5PUtKSuTn56fi4mL5+vpeSPkAgDOcPX1WNWqUn59vN33GVBoaoqmu3w4zYnTy5Emlp6crNjbWrj02NlY7d+6sdbtVq1Zp//79mjdvXlOXCACoB29vb7vlfO2AI3CYm68LCwtVUVFhNw8tnZ6LzsvLq3GbX375RbNmzVJaWprc3Op2KGVlZSorK7O9LikpaXjRAADAqTjMiFEVi8Vi99owjGptklRRUaG77rpLCxYsUK9eveq8/6SkJPn5+dmW0NDQC64ZAAA4B4e5x+jkyZNq27at3nvvPd1666229kceeUSZmZnaunWrXf9jx46pXbt2cnV1tbVVVlbKMAy5urpqy5Ytuv7666u9T00jRqGhodxjBABAK9JU9xg5zFSah4eHIiMjlZqaaheMUlNTdcstt1Tr7+vrq++//96ubfny5fr888/1/vvvKzw8vMb38fT0lKenZ+MWDwAAnILDBCNJSkxM1MSJExUVFaXo6Gi99tprysnJUUJCgiRp9uzZOnLkiN566y25uLioX79+dtt36tRJXl5e1doBAADqwqGC0YQJE1RUVKSFCxcqNzdX/fr1U0pKisLCwiRJubm5532mEQAAQEM5zD1GLYXnGAEA0Po4/XOMAAAAWhrBCAAAwEQwAgAAMBGMAAAATAQjAAAAE8EIAADARDACAAAwEYwAAABMBCMAAAATwQgAAMBEMAIAADARjAAAAEwEIwAAABPBCAAAwEQwAgAAMBGMAAAATAQjAAAAE8EIAADARDACAAAwEYwAAABMBCMAAAATwQgAAMBEMAIAADARjAAAAEwEIwAAABPBCAAAwEQwAgAAMBGMAAAATAQjAAAAE8EIAADARDACAAAwEYwAAABMBCMAAAATwQgAAMBEMAIAADARjAAAAEwEIwAAABPBCAAAwEQwAgAAMBGMAAAATAQjAAAAE8EIAADARDACAAAwEYwAAABMBCMAAAATwQgAAMBEMAIAADARjAAAAEwEIwAAABPBCAAAwEQwAgAAMDlcMFq+fLnCw8Pl5eWlyMhIpaWl1dp3w4YNGjFihDp27ChfX19FR0dr8+bNzVgtAABwJg4VjNatW6fp06drzpw5ysjI0NChQzVq1Cjl5OTU2H/btm0aMWKEUlJSlJ6eruuuu05xcXHKyMho5soBAIAzsBiGYbR0EVUGDx6sgQMHasWKFba2iIgIjR07VklJSXXaR9++fTVhwgTNnTu3Tv1LSkrk5+en4uJi+fr6NqhuAADQvJrq+u0wI0YnT55Uenq6YmNj7dpjY2O1c+fOOu2jsrJSx48fV/v27WvtU1ZWppKSErsFAABAcqBgVFhYqIqKCgUGBtq1BwYGKi8vr077WLx4saxWq8aPH19rn6SkJPn5+dmW0NDQC6obAAA4D4cJRlUsFovda8MwqrXVZO3atZo/f77WrVunTp061dpv9uzZKi4uti2HDh264JoBAIBzcGvpAqp06NBBrq6u1UaHCgoKqo0inW3dunWaMmWK3nvvPQ0fPvycfT09PeXp6XnB9QIAAOfjMCNGHh4eioyMVGpqql17amqqYmJiat1u7dq1io+P19tvv62bbrqpqcsEAABOzGFGjCQpMTFREydOVFRUlKKjo/Xaa68pJydHCQkJkk5Pgx05ckRvvfWWpNOhaNKkSXrxxRd19dVX20ab2rRpIz8/vxY7DgAA0Do5VDCaMGGCioqKtHDhQuXm5qpfv35KSUlRWFiYJCk3N9fumUavvvqqysvL9eCDD+rBBx+0tU+ePFmrV69u7vIBAEAr51DPMWoJPMcIAIDWx+mfYwQAANDSHGoqzZFVVFTo1KlTLV0GUC/u7u5ydXVt6TIAoNUgGJ2HYRjKy8vTsWPHWroUoEH8/f0VFBRUp+eBAcDFjmB0HlWhqFOnTmrbti0XF7QahmHoxIkTKigokCQFBwe3cEUAHIHVapWPj48kqbS0VN7e3i1ckWMhGJ1DRUWFLRQFBAS0dDlAvbVp00bS6QeldurUiWk1ADgPbr4+h6p7itq2bdvClQANV/X7yz1yAHB+jBjVAdNnaM34/QVgtVrP+7MkptVEMAIAwOlV3VN0trO/i/Qif7ShJKbScAaLxaKNGzees098fLzGjh1b530ePHhQFotFmZmZF1RbUzjzeOta57BhwzR9+vQmrw0A0DIYMXJS8fHxOnbs2HmDzplyc3PVrl07SaeDQnh4uDIyMnTFFVfY+rz44otO+RdFaGiocnNz1aFDB0nSl19+qeuuu05Hjx6Vv7+/rd+GDRvk7u7eQlUCQMOUlpbafrZarbaRovz8fKbPzkIwgk1QUNB5+zjrl/O6urrW6fjbt2/fDNUAQOOqLfx4e3sTjM7CVFo9GYYhq9XaIsuFjNQMGzZM06ZN02OPPab27dsrKChI8+fPt+tz5tRSeHi4JOnKK6+UxWLRsGHDJFWfStu0aZOuueYa+fv7KyAgQGPGjNH+/fvrVVtZWZkee+wxhYaGytPTUz179lRycrJt/datWzVo0CB5enoqODhYs2bNUnl5eb2O7ZdfftHf/vY3eXl5qU+fPkpNTbVbf+ZU2sGDB3XddddJktq1ayeLxaL4+Hjbe505lXb06FFNmjRJ7dq1U9u2bTVq1Cj98ssvtvWrV6+Wv7+/Nm/erIiICPn4+OjGG29Ubm5uvf6NAADNgxGjejpx4kStN7E1tQt9ENebb76pxMREff3119q1a5fi4+M1ZMgQjRgxolrf3bt3a9CgQfr000/Vt29feXh41LhPq9WqxMREXX755bJarZo7d65uvfVWZWZmysWlbrl70qRJ2rVrl5YuXaoBAwYoOztbhYWFkqQjR45o9OjRio+P11tvvaUff/xR9913n7y8vOzCz7mOrbKyUrfddps6dOigr776SiUlJee8Tyg0NFTr16/X7bffrp9++km+vr625wGdLT4+Xr/88os++ugj+fr66vHHH9fo0aO1b98+25TbiRMn9Nxzz+l///d/5eLiorvvvlszZ87UmjVr6vTvAwBoPgSji0j//v01b948SVLPnj318ssv67PPPqsxGHXs2FGSFBAQcM4ppttvv93udXJysjp16qR9+/apX79+563p559/1rvvvqvU1FQNHz5cktS9e3fb+uXLlys0NFQvv/yyLBaLLrvsMv3+++96/PHHNXfuXFv4Otexffrpp8rKytLBgwfVpUsXSdI//vEPjRo1qsaaXF1dbVNmnTp1srvH6ExVgWjHjh2KiYmRJK1Zs0ahoaHauHGj7rjjDkmnnx/0yiuvqEePHpKkhx56SAsXLjzvvw0ANAVvb2+nvFe0sRCM6qlt27Z2N7E193tfiP79+9u9Dg4Otn1dREPt379fTz31lL766isVFhaqsrJSkpSTk1OnYJSZmSlXV1dde+21Na7PyspSdHS03bN4hgwZotLSUh0+fFhdu3aVdO5jy8rKUteuXW2hSJKio6Prd6C11Obm5qbBgwfb2gICAtS7d29lZWXZ2tq2bWsLRWfXBgBwLASjerJYLK32RrWzP01lsVhsQaah4uLiFBoaqtdff10hISGqrKxUv379dPLkyTptX9sUVRXDMKo9oLDqL50z2891bDX9ZdQYDz2s7S+us2uuqTb+WgMAx8TN16hR1T1FFRUVtfYpKipSVlaWnnzySd1www2KiIjQ0aNH6/U+l19+uSorK7V169Ya1/fp00c7d+60CxI7d+7UJZdcos6dO9fpPfr06aOcnBz9/vvvtrZdu3adc5u6HH+fPn1UXl6ur7/+2tZWVFSkn3/+WREREXWqDQDgWAhGqFGnTp3Upk0bbdq0Sfn5+SouLq7Wp127dgoICNBrr72mX3/9VZ9//rkSExPr9T7dunXT5MmTde+992rjxo3Kzs7Wl19+qXfffVeS9MADD+jQoUN6+OGH9eOPP+rDDz/UvHnzlJiYWOebu4cPH67evXtr0qRJ+vbbb5WWlqY5c+acc5uwsDBZLBZ9/PHH+s9//lPj9GnPnj11yy236L777tP27dv17bff6u6771bnzp11yy231OvfAQDgGAhGqJGbm5uWLl2qV199VSEhITVe6F1cXPTOO+8oPT1d/fr104wZM/Tss8/W+71WrFihcePG6YEHHtBll12m++67z/b9PZ07d1ZKSop2796tAQMGKCEhQVOmTNGTTz5Z5/27uLjogw8+UFlZmQYNGqSpU6fqf/7nf865TefOnbVgwQLNmjVLgYGBeuihh2rst2rVKkVGRmrMmDGKjo6WYRhKSUnhIZAA0EpZjIv8ZoeSkhL5+fmpuLhYvr6+duv++usvZWdnKzw8XF5eXi1UIXBh+D0G4IzOdf2+EIwYAQAAmAhGAAAAJoIRAACAiWAEAABgIhgBAACYCEYAAAAmghEAAICJYAQAAGAiGAEAAJgIRrhg3bp105IlSy5oH19++aUsFouOHTvWKDUdPHhQFotFmZmZjbK/szVWvU1dJwCgfghGTm7nzp1ydXXVjTfe2NKl2AwbNkzTp0+3a4uJiVFubq78/PxapqhmEB8fr7Fjx9q1hYaGKjc3V/369WuZogAAdghGzcRqtcpischisdi+ILU5rFy5Ug8//LC2b9+unJycZnvf+vLw8FBQUJAsFktLl9KsXF1dFRQUJDc3t5YuBQAggpFTs1qtevfdd3X//fdrzJgxWr16td36qumgzz77TFFRUWrbtq1iYmL0008/2frs379ft9xyiwIDA+Xj46OrrrpKn376aa3vee+992rMmDF2beXl5QoKCtLKlSsVHx+vrVu36sUXX7QFxYMHD9Y4NbVjxw5de+21atu2rdq1a6eRI0fq6NGjkqRNmzbpmmuukb+/vwICAjRmzBjt37+/Xv8+y5cvV8+ePeXl5aXAwECNGzfOtq6srEzTpk1Tp06d5OXlpWuuuUbffPNNrfuaP3++rrjiCru2JUuWqFu3brb1b775pj788EPbcX/55Zc1TqVt3bpVgwYNkqenp4KDgzVr1iyVl5fb1g8bNkzTpk3TY489pvbt2ysoKEjz58+v17EDAGpGMHJi69atU+/evdW7d2/dfffdWrVqlQzDqNZvzpw5Wrx4sfbs2SM3Nzfde++9tnWlpaUaPXq0Pv30U2VkZGjkyJGKi4urdfRp6tSp2rRpk3Jzc21tKSkpKi0t1fjx4/Xiiy8qOjpa9913n3Jzc5Wbm6vQ0NBq+8nMzNQNN9ygvn37ateuXdq+fbvi4uJUUVEh6XToS0xM1DfffKPPPvtMLi4uuvXWW1VZWVmnf5s9e/Zo2rRpWrhwoX766Sdt2rRJf/vb32zrH3vsMa1fv15vvvmm9u7dq0svvVQjR47UH3/8Uaf9n23mzJkaP368brzxRttxx8TEVOt35MgRjR49WldddZW+/fZbrVixQsnJyXr66aft+r355pvy9vbW119/rUWLFmnhwoVKTU1tUG0AgDMYF7ni4mJDklFcXFxt3Z9//mns27fP+PPPPxu079LSUtuSn59vSDIkGfn5+XbrmkpMTIyxZMkSwzAM49SpU0aHDh2M1NRU2/ovvvjCkGR8+umntrZ//etfhqRzHnOfPn2Ml156yfY6LCzMeOGFF+zWP/PMM7bXY8eONeLj422vr732WuORRx6x22dVLUePHjUMwzD+67/+yxgyZEidj7WgoMCQZHz//feGYRhGdna2IcnIyMiosf/69esNX19fo6SkpNq60tJSw93d3VizZo2t7eTJk0ZISIixaNGiGuudN2+eMWDAALv9vPDCC0ZYWJjt9eTJk41bbrnFrs/ZdT7xxBNG7969jcrKSlufZcuWGT4+PkZFRYVhGKf//a655hq7/Vx11VXG448/XuOxXujvMQA4onNdvy8EI0ZNyMfHx7YEBgba2qumpaqWpvDTTz9p9+7duvPOOyVJbm5umjBhglauXFmtb//+/W0/BwcHS5IKCgoknR6Zeeyxx9SnTx/5+/vLx8dHP/744znvV5o6dapWrVpl28+//vUvu1GouqgaMarN/v37ddddd6l79+7y9fVVeHi4JNX5PqoRI0YoLCxM3bt318SJE7VmzRqdOHHCtu9Tp05pyJAhtv7u7u4aNGiQsrKy6nUc9ZWVlaXo6Gi7e62GDBmi0tJSHT582NZ25jmTTp+3qnMGAGg47vh0UsnJySovL1fnzp1tbYZhyN3dXUePHlW7du1s7e7u7rafqy7IVVNSjz76qDZv3qznnntOl156qdq0aaNx48bp5MmTtb73pEmTNGvWLO3atUu7du1St27dNHTo0HrV36ZNm3Ouj4uLU2hoqF5//XWFhISosrJS/fr1O2ddZ7rkkku0d+9effnll9qyZYvmzp2r+fPn65tvvrFNN559I7hhGLXeHO7i4lJtmvLUqVN1quV871FTPWees6p1dZ1GBADUjhGjJlRaWmpb8vPzbe35+fl26xpbeXm53nrrLS1evFiZmZm25dtvv1VYWJjWrFlT532lpaUpPj5et956qy6//HIFBQXp4MGD59wmICBAY8eO1apVq7Rq1Srdc889dus9PDxs9wrVpn///vrss89qXFdUVKSsrCw9+eSTuuGGGxQREWG7Kbs+3NzcNHz4cC1atEjfffedDh48qM8//1yXXnqpPDw8tH37dlvfU6dOac+ePYqIiKhxXx07dlReXp5dODr72UR1Oe4+ffpo586ddvvZuXOnLrnkEruQCwBoGowYNSFvb+9a22tb1xg+/vhjHT16VFOmTKn2XKBx48YpOTlZDz30UJ32demll2rDhg2Ki4uTxWLRU089VaeRialTp2rMmDGqqKjQ5MmT7dZ169ZNX3/9tQ4ePCgfHx+1b9++2vazZ8/W5ZdfrgceeEAJCQny8PDQF198oTvuuEPt27dXQECAXnvtNQUHBysnJ0ezZs2q0/FU+fjjj3XgwAH97W9/U7t27ZSSkqLKykr17t1b3t7euv/++/Xoo4+qffv26tq1qxYtWqQTJ05oypQpNe5v2LBh+s9//qNFixZp3Lhx2rRpkz755BP5+vraHffmzZv1008/KSAgoMZnNj3wwANasmSJHn74YT300EP66aefNG/ePCUmJsrFhb9jAKCp8X9aJ5ScnKzhw4fXeOG9/fbblZmZqb1799ZpXy+88ILatWunmJgYxcXFaeTIkRo4cOB5txs+fLiCg4M1cuRIhYSE2K2bOXOmXF1d1adPH3Xs2LHG+4J69eqlLVu26Ntvv9WgQYMUHR2tDz/8UG5ubnJxcdE777yj9PR09evXTzNmzNCzzz5bp+Op4u/vrw0bNuj6669XRESEXnnlFa1du1Z9+/aVJP3zn//U7bffrokTJ2rgwIH69ddftXnzZrspyDNFRERo+fLlWrZsmQYMGKDdu3dr5syZdn3uu+8+9e7dW1FRUerYsaN27NhRbT+dO3dWSkqKdu/erQEDBighIUFTpkzRk08+Wa/jAwA0jMU4+8aIi0xJSYn8/PxUXFxs99e9JP3111/Kzs5WeHi4vLy8Luh9rFar7Ubr0tLSJh0xcgQnTpxQSEiIVq5cqdtuu62ly7moNebvMQA4inNdvy8EU2nNxNvbu8ZnCDmbyspK5eXlafHixfLz89PNN9/c0iUBAFBnBCM0qpycHIWHh6tLly5avXo1X3UBAGhVuGqhUXXr1u2iGBkDADgnbr4GAAAwEYzqgBEQtGb8/gJA3RGMzqHq6cJVXxUBtEZVv79nPy0bAFAd9xidg6urq/z9/W3fQdW2bdtavxICcDSGYejEiRMqKCiQv7+/XF1dW7okAHB4BKPzCAoKkiS+oBOtlr+/v+33GABwbgSj87BYLAoODlanTp0a9KWgQEtyd3dnpAgA6oFgVEeurq5cYAAAcHIOd/P18uXLbV9dEBkZqbS0tHP237p1qyIjI+Xl5aXu3bvrlVdeaaZKAQCAs3GoYLRu3TpNnz5dc+bMUUZGhoYOHapRo0bV+CWjkpSdna3Ro0dr6NChysjI0BNPPKFp06Zp/fr1zVw5AABwBg71JbKDBw/WwIEDtWLFCltbRESExo4dq6SkpGr9H3/8cX300UfKysqytSUkJOjbb7/Vrl276vSeTfUldM2p6tNHAADUR2v+tLXTf4nsyZMnlZ6erlmzZtm1x8bGaufOnTVus2vXLsXGxtq1jRw5UsnJyTp16lSNz20pKytTWVmZ7XVxcbGk0//ArZXValVISEhLlwEAaGV+//13eXt7t3QZDVJ13W7s8R2HCUaFhYWqqKhQYGCgXXtgYKDy8vJq3CYvL6/G/uXl5SosLFRwcHC1bZKSkrRgwYJq7aGhoRdQPQAArY8z/FFdVFQkPz+/RtufwwSjKmcP6RmGcc5hvpr619ReZfbs2UpMTLS9PnbsmMLCwpSTk9Oo/7BomJKSEoWGhurQoUOtdmrTWXAuHAfnwnFwLhxHcXGxunbtqvbt2zfqfh0mGHXo0EGurq7VRocKCgqqjQpVCQoKqrG/m5ubAgICatzG09NTnp6e1dr9/Pz4JXcgvr6+nA8HwblwHJwLx8G5cBwuLo37OTKH+VSah4eHIiMjlZqaateempqqmJiYGreJjo6u1n/Lli2Kiorie6EAAEC9OUwwkqTExES98cYbWrlypbKysjRjxgzl5OQoISFB0ulpsEmTJtn6JyQk6LffflNiYqKysrK0cuVKJScna+bMmS11CAAAoBVzmKk0SZowYYKKioq0cOFC5ebmql+/fkpJSVFYWJgkKTc31+6ZRuHh4UpJSdGMGTO0bNkyhYSEaOnSpbr99tvr/J6enp6aN29ejdNraH6cD8fBuXAcnAvHwblwHE11LhzqOUYAAAAtyaGm0gAAAFoSwQgAAMBEMAIAADARjAAAAEwXRTBavny5wsPD5eXlpcjISKWlpZ2z/9atWxUZGSkvLy91795dr7zySjNV6vzqcy42bNigESNGqGPHjvL19VV0dLQ2b97cjNU6v/r+t1Flx44dcnNz0xVXXNG0BV5E6nsuysrKNGfOHIWFhcnT01M9evTQypUrm6la51bfc7FmzRoNGDBAbdu2VXBwsO655x4VFRU1U7XOa9u2bYqLi1NISIgsFos2btx43m0a5fptOLl33nnHcHd3N15//XVj3759xiOPPGJ4e3sbv/32W439Dxw4YLRt29Z45JFHjH379hmvv/664e7ubrz//vvNXLnzqe+5eOSRR4xnnnnG2L17t/Hzzz8bs2fPNtzd3Y29e/c2c+XOqb7no8qxY8eM7t27G7GxscaAAQOap1gn15BzcfPNNxuDBw82UlNTjezsbOPrr782duzY0YxVO6f6nou0tDTDxcXFePHFF40DBw4YaWlpRt++fY2xY8c2c+XOJyUlxZgzZ46xfv16Q5LxwQcfnLN/Y12/nT4YDRo0yEhISLBru+yyy4xZs2bV2P+xxx4zLrvsMru2v//978bVV1/dZDVeLOp7LmrSp08fY8GCBY1d2kWpoedjwoQJxpNPPmnMmzePYNRI6nsuPvnkE8PPz88oKipqjvIuKvU9F88++6zRvXt3u7alS5caXbp0abIaL0Z1CUaNdf126qm0kydPKj09XbGxsXbtsbGx2rlzZ43b7Nq1q1r/kSNHas+ePTp16lST1ersGnIuzlZZWanjx483+hcGXowaej5WrVql/fv3a968eU1d4kWjIefio48+UlRUlBYtWqTOnTurV69emjlzpv7888/mKNlpNeRcxMTE6PDhw0pJSZFhGMrPz9f777+vm266qTlKxhka6/rtUE++bmyFhYWqqKio9iW0gYGB1b58tkpeXl6N/cvLy1VYWKjg4OAmq9eZNeRcnG3x4sWyWq0aP358U5R4UWnI+fjll180a9YspaWlyc3Nqf/X0awaci4OHDig7du3y8vLSx988IEKCwv1wAMP6I8//uA+owvQkHMRExOjNWvWaMKECfrrr79UXl6um2++WS+99FJzlIwzNNb126lHjKpYLBa714ZhVGs7X/+a2lF/9T0XVdauXav58+dr3bp16tSpU1OVd9Gp6/moqKjQXXfdpQULFqhXr17NVd5FpT7/bVRWVspisWjNmjUaNGiQRo8ereeff16rV69m1KgR1Odc7Nu3T9OmTdPcuXOVnp6uTZs2KTs72/Ydn2hejXH9duo/+zp06CBXV9dqSb+goKBaqqwSFBRUY383NzcFBAQ0Wa3OriHnosq6des0ZcoUvffeexo+fHhTlnnRqO/5OH78uPbs2aOMjAw99NBDkk5fnA3DkJubm7Zs2aLrr7++WWp3Ng35byM4OFidO3eWn5+frS0iIkKGYejw4cPq2bNnk9bsrBpyLpKSkjRkyBA9+uijkqT+/fvL29tbQ4cO1dNPP80sQzNqrOu3U48YeXh4KDIyUqmpqXbtqampiomJqXGb6Ojoav23bNmiqKgoubu7N1mtzq4h50I6PVIUHx+vt99+mzn7RlTf8+Hr66vvv/9emZmZtiUhIUG9e/dWZmamBg8e3FylO52G/LcxZMgQ/f777yotLbW1/fzzz3JxcVGXLl2atF5n1pBzceLECbm42F9KXV1dJf3/oxVoHo12/a7XrdqtUNVHL5OTk419+/YZ06dPN7y9vY2DBw8ahmEYs2bNMiZOnGjrX/VxvxkzZhj79u0zkpOT+bh+I6nvuXj77bcNNzc3Y9myZUZubq5tOXbsWEsdglOp7/k4G59Kazz1PRfHjx83unTpYowbN8744YcfjK1btxo9e/Y0pk6d2lKH4DTqey5WrVpluLm5GcuXLzf2799vbN++3YiKijIGDRrUUofgNI4fP25kZGQYGRkZhiTj+eefNzIyMmyPTmiq67fTByPDMIxly5YZYWFhhoeHhzFw4EBj69attnWTJ082rr32Wrv+X375pXHllVcaHh4eRrdu3YwVK1Y0c8XOqz7n4tprrzUkVVsmT57c/IU7qfr+t3EmglHjqu+5yMrKMoYPH260adPG6NKli5GYmGicOHGimat2TvU9F0uXLjX69OljtGnTxggODjb++7//2zh8+HAzV+18vvjii3NeA5rq+m0xDMb6AAAAJCe/xwgAAKA+CEYAAAAmghEAAICJYAQAAGAiGAEAAJgIRgAAACaCEQAAgIlgBAAAYCIYAQAAmAhGAJxKUVGROnXqpIMHD17QfsaNG6fnn3++cYoC0GrwlSAAnMrMmTN19OhRJScnX9B+vvvuO1133XXKzs6Wr69vI1UHwNExYgTAafz5559KTk7W1KlTL3hf/fv3V7du3bRmzZpGqAxAa0EwAuCw1q5dKy8vLx05csTWNnXqVPXv31/FxcXV+n/yySdyc3NTdHS0XfuwYcP08MMPa/r06WrXrp0CAwP12muvyWq16p577tEll1yiHj166JNPPrHb7uabb9batWub5uAAOCSCEQCHdeedd6p3795KSkqSJC1YsECbN2/WJ598Ij8/v2r9t23bpqioqBr39eabb6pDhw7avXu3Hn74Yd1///264447FBMTo71792rkyJGaOHGiTpw4Ydtm0KBB2r17t8rKyprmAAE4HO4xAuDQPv74Y40bN05z587Vc889p7S0NPXt27fGvmPHjlVAQEC1+4uGDRumiooKpaWlSZIqKirk5+en2267TW+99ZYkKS8vT8HBwdq1a5euvvpqSafvMxowYIAOHjyosLCwJjxKAI7CraULAIBzGTNmjPr06aMFCxZoy5YttYYi6fQ9Rl5eXjWu69+/v+1nV1dXBQQE6PLLL7e1BQYGSpIKCgpsbW3atJEku1EkAM6NqTQADm3z5s368ccfVVFRYQsvtenQoYOOHj1a4zp3d3e71xaLxa7NYrFIkiorK21tf/zxhySpY8eODaodQOtDMALgsPbu3as77rhDr776qkaOHKmnnnrqnP2vvPJK7du3r9He/9///re6dOmiDh06NNo+ATg2ghEAh3Tw4EHddNNNmjVrliZOnKiFCxdq/fr1Sk9Pr3WbkSNH6ocffqh11Ki+0tLSFBsb2yj7AtA6EIwAOJw//vhDo0aN0s0336wnnnhCkhQZGam4uDjNmTOn1u0uv/xyRUVF6d13373gGv766y998MEHuu+++y54XwBaDz6VBsCppKSkaObMmfr3v/8tF5eG/+23bNkyffjhh9qyZUsjVgfA0fGpNABOZfTo0frll1905MgRhYaGNng/7u7ueumllxqxMgCtASNGAAAAJu4xAgAAMBGMAAAATAQjAAAAE8EIAADARDACAAAwEYwAAABMBCMAAAATwQgAAMBEMAIAADARjAAAAEz/H0AT3uRVJrysAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Step 2.1: set initial condition\n", "phi = np.zeros(nx) # a numpy array with nx elements all equal to zero\n", "\n", "# Step 2.2: calculate analytical solution\n", "phi_true = (np.exp(F*x/Gamma) - 1) / (np.exp(F*L/Gamma) - 1) * (phiB - phiA) + phiA\n", "\n", "# Step 2.3: visualize initial condition and analytical solution\n", "plt.figure() # create a new figure\n", "plt.plot(x, phi, '-k', label='Initial condition')\n", "plt.scatter(x, phi_true, s=50, c='k', marker='+', label='Analytical solution')\n", "plt.xlabel('$x$ (m)')\n", "plt.ylabel(r'$\\phi$')\n", "ax = plt.gca()\n", "ax.set_xlim(0, 1.0)\n", "ax.set_ylim(0, 1.2)\n", "plt.legend(loc='lower left')\n", "plt.show() # show the figure" ] }, { "cell_type": "code", "execution_count": 23, "id": "a17c56eb", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0.70967742 0. 0. 0. 0. ]\n", "[0.70967742 0.39032258 0. 0. 0. ]\n", "[0.82299688 0.39032258 0.21467742 0. 0. ]\n", "[0.82299688 0.54925312 0.21467742 0.11807258 0. ]\n", "[0.869138 0.54925312 0.35522188 0.11807258 0.04478615]\n", "[0.869138 0.63787575 0.35522188 0.2155258 0.04478615]\n", "[0.89486715 0.63787575 0.44781827 0.2155258 0.08175117]\n", "[0.89486715 0.69369516 0.44781827 0.28308807 0.08175117]\n", "[0.91107279 0.69369516 0.50892197 0.28308807 0.10737823]\n", "[0.91107279 0.73010492 0.50892197 0.32822729 0.10737823]\n", "[0.92164336 0.73010492 0.54925999 0.32822729 0.12450001]\n", "[0.92164336 0.75407084 0.54925999 0.35811799 0.12450001]\n", "[0.92860121 0.75407084 0.57589206 0.35811799 0.13583786]\n", "[0.92860121 0.76988209 0.57589206 0.37786767 0.13583786]\n", "[0.93319158 0.76988209 0.5934756 0.37786767 0.14332912]\n", "[0.93319158 0.78031939 0.5934756 0.39090968 0.14332912]\n", "[0.93622176 0.78031939 0.60508502 0.39090968 0.14827609]\n", "[0.93622176 0.78721023 0.60508502 0.399521 0.14827609]\n", "[0.93822232 0.78721023 0.61275008 0.399521 0.15154245]\n", "[0.93822232 0.79175981 0.61275008 0.40520664 0.15154245]\n", "[0.93954317 0.79175981 0.61781089 0.40520664 0.15369907]\n", "******************************************\n", "Final temperature difference: 0.0085383\n", "Number of iterations: 21\n", "Elapsed time: 0.010 seconds\n" ] } ], "source": [ "# Step 3: finite volume calculations\n", "phi_old = np.zeros_like(phi) # placeholder array to advance the solution\n", "phi_diff = 1 # phi difference for convergence\n", "cnt = 0 # counter for the number of iterations\n", "t_start = time.perf_counter() # start time for performance measurement\n", "\n", "while phi_diff > 1e-2: # loop until convergence\n", " cnt += 1 # increment the counter\n", " phi_old = phi.copy() # copy the current solution to the old solution\n", "\n", " # left boundary condition\n", " a_W, a_E = 0, D - F/2 # coefficients for the left boundary\n", " S_P = -(2*D + F)\n", " a_P = a_W + a_E - S_P\n", " S_u = (2*D + F) * phiA\n", " phi[0] = (a_E * phi_old[1] + S_u) / a_P # update the left boundary condition\n", "\n", " # internal nodes\n", " a_W, a_E = D + F/2, D - F/2 # coefficients for the internal nodes\n", " a_P = 2*D\n", " for i in range(1, nx-1): # loop over the interior grid points\n", " phi[i] = (a_W * phi_old[i-1] + a_E * phi[i+1]) / a_P\n", " \n", " # right boundary condition\n", " a_W, a_E = D + F/2, 0 # coefficients for the right boundary\n", " S_P = -(2*D - F)\n", " a_P = a_W + a_E - S_P\n", " S_u = (2*D - F) * phiB\n", " phi[-1] = (a_W * phi_old[-2] + S_u) / a_P # update the right boundary condition\n", "\n", " # calculate the temperature difference for convergence\n", " phi_diff = np.sum(np.abs(phi - phi_old)) # total difference \n", " print(phi)\n", "\n", "# stop the timer and print the iteration results\n", "t_end = time.perf_counter()\n", "print('******************************************')\n", "print('Final temperature difference: {:.7f}'.format(phi_diff))\n", "print('Number of iterations: {}'.format(cnt))\n", "print('Elapsed time: {:.3f} seconds'.format(t_end - t_start))" ] }, { "cell_type": "code", "execution_count": 24, "id": "2df4435e", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkYAAAG4CAYAAACgrSiYAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAVL9JREFUeJzt3XtcT/fjB/DX6fPpnqIQkRC1LkQhcp+Uy8LGZDbJsLW5J/fNbb7CXDKEUcxGLnPbyCVzS7mVcqu5E1aa0E3S5fz+kM9vnxWrVOfTp9fz8TiPh8/7c875vD4dfF6dcz7nCKIoiiAiIiIiaEgdgIiIiEhVsBgRERERFWAxIiIiIirAYkRERERUgMWIiIiIqACLEREREVEBFiMiIiKiAixGRERERAVYjIiIiIgKsBgRERERFVCpYnTy5El4eHjAzMwMgiBgz549b51/165d6N69O2rVqgVDQ0O0a9cOhw4dqpiwREREpHZUqhhlZmbCwcEBK1euLNb8J0+eRPfu3REaGoro6Gh07doVHh4eiImJKeekREREpI4EVb2JrCAI2L17N/r161ei5ezs7ODp6YmZM2eWTzAiIiJSW3KpA5Sl/Px8pKenw9jY+I3zZGdnIzs7W2mZJ0+ewMTEBIIgVERMIiIiekeiKCI9PR1mZmbQ0Ci7A2BqVYyWLFmCzMxMDBw48I3z+Pv7Y86cORWYioiIiMrL/fv3Ub9+/TJbn9ocSgsJCcGIESOwd+9euLq6vnG+f+8xSk1NRYMGDXD//n0YGhq+a2wiIiKqAGlpaTA3N8ezZ89gZGRUZutViz1G27Ztw/Dhw7Fjx463liIA0NbWhra2dqFxQ0NDFiMiIqJKpqxPg1Gpb6WVRkhICLy9vbFlyxb07t1b6jhERERUianUHqOMjAzcvHlT8fjOnTuIjY2FsbExGjRogGnTpuHhw4fYtGkTgFelyMvLC8uXL0fbtm2RlJQEANDV1S3T3WpERERUNajUHqOoqCi0bNkSLVu2BAD4+vqiZcuWiq/eJyYmIiEhQTH/2rVrkZubi1GjRqFu3bqKady4cZLkJyIiospNZU++rihpaWkwMjJCamoqzzEiIiKqJMrr81ul9hgRERERSYnFiIiIiKgAixERERFRARYjIiIiogIsRkREREQFWIyIiIiICrAYERERERVgMSIiIiIqwGJEREREVIDFiIiIiKgAixERERFRARYjIiIiogIsRkREREQFWIyIiIiICrAYERERERVgMSIiIiIqwGJEREREVIDFiIiIiKgAixERERFRARYjIiIiogIsRkREREQFWIyIiIiICrAYERERERVgMSIiIiIqwGJEREREVIDFiIiIiKgAixERERFRARYjIiIiogIsRkREREQFWIyIiIiICrAYERERERVgMSIiIiIqwGJEREREVIDFiIiIiKgAixERERFRARYjIiIiogIsRkREREQFWIyIiIiICrAYERERERVgMSIiIiIqwGJEREREVIDFiIiIiKgAixERERFRARYjIiIiogIsRkREREQFWIyIiIiICrAYERERERVgMSIiIiIqwGJEREREVIDFiIiIiKgAixERERFRAZUqRidPnoSHhwfMzMwgCAL27Nnzn8ucOHECTk5O0NHRQePGjbFmzZryD0pERERqSaWKUWZmJhwcHLBy5cpizX/nzh306tULHTt2RExMDKZPn46xY8di586d5ZyUiIiI1JFc6gD/1LNnT/Ts2bPY869ZswYNGjRAQEAAAMDGxgZRUVFYvHgx+vfvX04piYiISF2p1B6jkjp9+jTc3NyUxtzd3REVFYWcnJwil8nOzkZaWprSRERERARU8mKUlJQEU1NTpTFTU1Pk5ubi8ePHRS7j7+8PIyMjxWRubl4RUYmIiKgSqNTFCAAEQVB6LIpikeOvTZs2DampqYrp/v375Z6RiIiIKgeVOseopOrUqYOkpCSlseTkZMjlcpiYmBS5jLa2NrS1tSsiHhEREVUylXqPUbt27RAWFqY0dvjwYbRq1QqampoSpSIiIqLKSqWKUUZGBmJjYxEbGwvg1dfxY2NjkZCQAODVYTAvLy/F/D4+Prh37x58fX0RHx+P4OBgBAUFwc/PT4r4REREVMmp1KG0qKgodO3aVfHY19cXADB06FBs3LgRiYmJipIEAI0aNUJoaCgmTJiAVatWwczMDD/88AO/qk9ERESlIoivz1auotLS0mBkZITU1FQYGhpKHYeIiIiKobw+v1XqUBoRERGRlFiMiIiIiAqwGBEREREVYDEiIiIiKsBiRERERFSAxYiIiIioAIsRERERUQEWIzWQmZkJQRAgCAIyMzOljkNERFRpsRgRERERFWAxIiIiIiqgUvdKo+L75yGzf/756tWrsLOzUzzW19ev0FxERESVGYtRJWVgYFDkuLOzs9LjKn4rPCIiohLhoTQ1d+DAAZYjIiKiYmIxqqQyMjIU06NHjxTjR48exaBBgyCXv9oZ2KtXLzg4OGDTpk14+fKlVHGJiIgqBRajSkpfX19peq1NmzYICQnBrVu3MGHCBBgYGODy5csYOnQoLC0tsWTJEqSlpUmYnIiISHWxGKmpBg0aYOnSpUhISIC/vz/q1KmDBw8ewM/PDw0aNMDUqVPx119/SR2TiIhIpbAYqbkaNWpg6tSpuHv3LtavXw9ra2ukpqZi4cKFaNiwIYYPH474+HipYxIREakEFiM1oK+vD1EUIYriG7+er62tjeHDhyMuLg579+5F+/btkZOTg+DgYNja2qJPnz4IDw/nidpERFSlsRhVMRoaGujTpw9OnTqFiIgI9OvXD4Ig4Pfff0enTp3g4uKCXbt2IS8vT+qoREREFY7FqApzcXHB7t27ER8fjy+++ALa2to4c+YM+vfvDxsbG6xduxZZWVlSxyQiIqowLEYEa2trrF27Fvfu3cOMGTNQvXp13LhxAz4+PrCwsMC8efPw5MkTqWMSERGVO0Gs4ieVpKWlwcjICKmpqTA0NJQ6jkrIyMhAUFCQ4lttwKvzmIYPH44JEyagYcOG0gYkIqIqr7w+v7nHiAoxMDDAuHHjcPPmTfzyyy9wcHBAZmYmfvjhBzRp0gSDBw9GTEyM1DGJiIjKHIsRvZGmpiY+/fRTxMTE4NChQ3B1dUVeXh5CQkLg6OiI7t274/Dhw/wmGxERqQ0WI/pPgiDAzc0NYWFhuHDhAj755BPIZDIcOXIE7u7uaNmyJTZv3oycnBypoxIREb0TFiMqkZYtW2LLli24efMmxo4dCz09PVy8eBGfffYZmjRpgoCAAGRkZEgdk4iIqFR48jVPvn4nKSkpWL16NVasWIHk5GQAr662/dVXX2HMmDGoU6eOxAmJiEgdldfnN4sRi1GZyMrKwqZNm7BkyRLcuHEDwKurbXt5eWHixImwtraWOCEREakTfiuNVJquri6+/PJLxMfHY9euXWjbti2ys7Oxbt062NjY4MMPP0RkZKTUMYmIiN6KxYjKlEwmU5Sg8PBweHh4QBRF7NmzB+3bt0eHDh2wd+9e5OfnSx2ViIioEBYjKheCIKBDhw747bffEBcXh88//xxaWlqK+7PZ2tpi/fr1ePHihdRRiYiIFFiMqNzZ2NggKCgId+7cwZQpU2BkZIRr165h5MiRaNSoEfz9/fH06VOpYxIREfHka558XfHS0tKwfv16LFu2DA8ePADw6mrbI0eOxPjx49GgQQOJExIRkarjydekNgwNDeHr64tbt25h06ZNsLe3R0ZGBpYtWwZLS0sMGTIEly5dkjomERFVQSxGJBktLS1FCTpw4AC6du2K3Nxcxf3ZevTogaNHj/KWI0REVGFYjEhygiAoStD58+cxcOBAaGho4NChQ+jWrRtatWqFbdu2ITc3V+qoRESk5liMSKW8LkE3btzAqFGjoKuriwsXLmDQoEFo2rQpVqxYgczMTKljEhGRmuLJ1zz5WqU9fvwYq1atwsqVK/H48WMAgLGxMUaNGoXRo0ejdu3aEickIiIp8JYg5YTFqHJ4/vw5Nm7ciCVLluD27dsAAB0dHXh7e2PixIlo0qSJxAmJiKgi8VtpVKXp6enh66+/xvXr17F9+3a0bt0aL168wJo1a2BlZYUBAwbg7NmzUsckIqJKjsWIKhWZTIaPP/4YZ8+exbFjx9CrVy+IooidO3eibdu26Ny5M/bt28dbjhARUamwGFGlJAgCunTpgv379+Py5csYOnQoNDU1cfLkSXh4eKBZs2bYsGEDsrOzpY5KRESVCIsRVXr29vbYuHEjbt++DT8/P1SrVk1xf7bGjRtj0aJFSE1NlTomERFVAjz5midfq53U1FSsXbsWy5cvx19//QUAqFatGr788kuMGzcO9evXlzghERG9K558TVRMRkZGmDx5Mm7fvo3g4GDY2toiPT0dixcvRuPGjeHt7Y0rV65IHZOIiFQQixGpLW1tbQwbNgyXL1/Gvn370KlTJ+Tk5OCnn35Cs2bN0Lt3bxw/fpy3HCEiIgUWI1J7Ghoa6N27N06cOIEzZ86gf//+EAQBoaGh6Nq1K5ydnbFjxw7k5eVJHZWIiCTGYkRVirOzM3799Vdcu3YNPj4+0NHRUdyfzdraGqtXr0ZWVpbUMYmISCI8+ZonX1dpycnJWLlyJVatWoUnT54AAGrWrIkxY8bg66+/Rs2aNSVOSEREReEtQcoJixEBQGZmJoKDg7F06VLcvXsXAKCrq4vPP/8cvr6+aNy4sbQBiYhISZX5VlpgYCAaNWoEHR0dODk5ITw8/K3zb968GQ4ODtDT00PdunUxbNgwpKSkVFBaUhf6+voYM2YMbty4ga1bt8LR0RFZWVlYtWoVmjZtCk9PT0RFRf3nejIzMyEIAgRBQGZmZgUkJyKisqRSxWjbtm0YP348ZsyYgZiYGHTs2BE9e/ZEQkJCkfOfOnUKXl5eGD58OK5evYodO3bg/PnzGDFiRAUnJ3Uhl8sVJejIkSNwd3dHfn6+4v5s77//Pg4ePMhvshERqSmVKkZLly7F8OHDMWLECNjY2CAgIADm5uZYvXp1kfOfOXMGDRs2xNixY9GoUSN06NABX375ZbF+syd6G0EQ0K1bNxw8eBCxsbH47LPPIJPJcOzYMfTs2RMODg74+eefkZOTI3VUIiIqQypTjF6+fIno6Gi4ubkpjbu5uSEyMrLIZVxcXPDgwQOEhoZCFEU8evQIv/76K3r37v3G18nOzkZaWprSRPQ2r0vQ7du3MWHCBBgYGODy5cvw8vJC48aN4e/vj8TERGRmZiodPnv9+N/jRESkulSmGD1+/Bh5eXkwNTVVGjc1NUVSUlKRy7i4uGDz5s3w9PSElpYW6tSpg+rVq2PFihVvfB1/f38YGRkpJnNz8zJ9H6S+GjRogKVLlyIhIQHz58+HqakpHjx4gOnTp8PMzAwGBgZKf39NTU1hYGCgmIiISPWpTDF6TRAEpceiKBYaey0uLg5jx47FzJkzER0djYMHD+LOnTvw8fF54/qnTZuG1NRUxXT//v0yzU/qr0aNGpg2bRru3r2LdevWSR2HiIjKkMoUo5o1a0ImkxXaO5ScnFxoL9Jr/v7+aN++PSZNmoTmzZvD3d0dgYGBCA4ORmJiYpHLaGtrw9DQUGkiKg0dHR2MGDECaWlp2Lp1K9q1a6f0vKenJ+Lj45GRkYGMjAyJUhIRUUmoTDHS0tKCk5MTwsLClMbDwsLg4uJS5DLPnz+HhobyW5DJZADAbw1RhalWrRo8PT0RGRmJI0eOKMa3bduGFi1aYNasWXjx4oWECYmIqLhUphgBgK+vL9avX4/g4GDEx8djwoQJSEhIUBwamzZtGry8vBTze3h4YNeuXVi9ejVu376NiIgIjB07Fm3atIGZmZlUb4OqsLZt2yr+3LFjR2RnZ2PJkiWwtLTEwoULebsRIiIVp1LFyNPTEwEBAZg7dy5atGiBkydPIjQ0FBYWFgCAxMREpWsaeXt7Y+nSpVi5ciXs7e3x8ccfw9raGrt27ZLqLRAphIaGIjQ0FM2aNUNqaiqmTp2Kpk2bIigoCLm5uVLHIyKiIvCWILwlCJWzvLw8bNmyBd98842i2Nva2sLf3x8eHh5v/HIBERG9WZW5JQiRupHJZBgyZAiuXbuGJUuWwNjYGHFxcejbty86duyIiIgIqSMSEVEBFiOiCqKjowNfX1/cunUL06ZNg66uLiIiItChQwf069cP8fHxUkckIqryWIyIKlj16tUxf/583LhxAyNGjICGhgb27t0Le3t7jBw5Eg8fPpQ6IhFRlcViRCSRevXqYd26dbhy5Qr69euH/Px8rF+/Hk2bNsX06dPx7NkzqSMSEVU5LEZEErOxscHu3bsRERGB9u3bIysrC/7+/rC0tMTSpUuRnZ0tdUQioiqDxYhIRbi4uCA8PBx79+6FjY0Nnjx5gokTJ8La2ho///wz8vLypI5IRKT2WIyIVIggCOjTpw8uXbqEoKAg1KtXD/fu3YOXlxccHR1x8OBBXtWdiKgcsRgRqSC5XI7PP/8c169fx4IFC2BkZIRLly6hZ8+e6NatG86fPy91RCIitcRiRKTC9PT0MGXKFNy+fRsTJ06ElpYWjh07hjZt2mDgwIG4ceOG1BGJiNQKixFRJWBsbIzFixfj+vXrGDp0KARBwI4dO2Bra4tRo0bh0aNHUkckIlILLEZElYiFhQU2btyI2NhY9OrVC7m5uQgMDISlpSVmzZqF9PR0qSMSEVVqLEZElVDz5s2xf/9+xWG1zMxMzJ07F5aWlli5ciVevnwpdUQiokqJxYioEuvSpQvOnDmDHTt2oGnTpvj7778xZswY2NjYYOvWrcjPz5c6IhFRpcJiRFTJCYKAAQMG4OrVq1i9ejVMTU1x+/ZtfPLJJ2jTpg3++OMPqSMSEVUaLEZEakJTUxM+Pj64efMm5s6dCwMDA0RHR8PV1RXu7u6IiYmROiIRkcpjMSJSMwYGBvj2229x69YtjB07Fpqamjh8+DAcHR3x2Wef4c6dO1JHJCJSWSxGRGqqdu3aWL58OeLj4/HJJ58AADZv3gxra2uMHz8ejx8/ljghEZHqYTEiUnOWlpbYsmWL4rBaTk4Oli9fjsaNG+N///sfMjMzpY5IRKQyWIyIqghHR0eEhYXh8OHDaNmyJdLT0/HNN9+gSZMmWLt2LXJzc6WOSEQkORYjoiqme/fuiIqKwpYtW9CoUSMkJSXBx8cHdnZ22LVrF29SS0RVGosRURWkoaGBTz75BPHx8Vi+fDlq1qyJ69evo3///mjXrh1OnjwpdUQiIkmwGBFVYdra2hg7dixu3bqFb7/9Fnp6ejh79iw6d+4MDw8PXLlyReqIREQVisWIiGBoaIi5c+fi5s2b8PHxgUwmw759+9C8eXMMGzYMCQkJUkckIqoQLEZEpFC3bl2sXr0acXFxGDBgAERRxMaNG2FlZYXJkyfjyZMnUkckIipXLEZEVIiVlRV27NiBM2fOoFOnTsjOzsb3338PS0tLLFq0CFlZWVJHJCIqFyxGRPRGzs7OOH78OPbv3w97e3s8e/YMU6ZMgZWVFTZs2IC8vDypIxIRlSkWIyJ6K0EQ0KtXL8TGxmLjxo0wNzfHgwcP8Pnnn8PBwQG///47v+JPRGqDxYiIikUmk2Ho0KG4fv06Fi9ejBo1auDq1avo06cPOnfujNOnT0sdkYjonbEYEVGJ6OjoYOLEibh9+zamTJkCHR0dhIeHw8XFBR999BH+/PNPqSMSEZUaixERlUr16tWxYMEC3LhxA8OHD4eGhgZ2794Ne3t7fPnll/jrr7+kjkhEVGIsRkT0TurXr4/169fj0qVL6NOnD/Ly8vDjjz+iSZMmmDFjBlJTU6WOSERUbCxGRFQm7OzssHfvXsVhtaysLMyfPx+WlpYICAhAdna21BGJiP4TixERlakOHTrg1KlT2L17N9577z2kpKRgwoQJsLa2xi+//IL8/HypIxIRvRGLERGVOUEQ0K9fP1y+fBnr1q2DmZkZ7t27hyFDhsDR0RGHDh3iV/yJSCWxGBFRuZHL5RgxYgRu3LiB+fPnw9DQEBcvXkSPHj3g6uqKqKgoqSMSESlhMSKicqenp4dp06bh9u3b8PX1hZaWFo4ePYrWrVvD09MTN2/elDoiEREAFiMiqkAmJiZYsmQJrl+/jiFDhkAQBGzfvh02NjYYPXo0kpOTpY5IRFUcixERVTgLCwts2rQJMTEx6NmzJ3Jzc7Fq1SpYWlpizpw5SE9PlzoiEVVRLEZEJBkHBweEhoYqDqtlZGRg9uzZaNKkCQIDA5GTkyN1RCKqYliMiEhyXbt2xdmzZ7F9+3Y0adIEycnJGDVqFGxtbbF9+3Z+g42IKgyLERGpBEEQ8PHHHyMuLg6BgYEwNTXFzZs34enpiTZt2uDYsWNSRySiKoDFiIhUiqamJr766ivcvHkTc+bMgYGBAaKiovD++++jZ8+euHjxotQRiUiNsRgRkUoyMDDAzJkzcevWLYwePRpyuRwHDx5Ey5YtMWTIENy9e1fqiESkhliMiEil1a5dGytWrMCff/6JQYMGQRRF/PLLL7C2toavry8eP35c5HKZmZkQBAGCICAzM7OCUxNRZcViRESVgqWlJUJCQhAVFYVu3brh5cuXWLZsGSwtLTF//nw8f/5c6ohEpAZYjIioUnFyckJYWBgOHTqEFi1aIC0tDTNmzECTJk2wbt065ObmSh2RiCoxFiMiqnQEQYCbmxuio6OxefNmNGzYEImJifjiiy9ga2uLkJAQZGRkKObPzMxUmoiI3kQQq/gFQtLS0mBkZITU1FQYGhpKHYeISiE7Oxtr1qzB+PHjizV/Ff9vj0gtlNfnN/cYEVGlp62tjXHjxkkdg4jUAIsREamNjIwMZGRk4ObNm/Dy8lKM6+rqYtGiRUhLS1M6xEZE9G88lMZDaURqKTMzEwYGBkpjHTp0QHBwMJo2bSpRKiIqK1XmUFpgYCAaNWoEHR0dODk5ITw8/K3zZ2dnY8aMGbCwsIC2tjYsLS0RHBxcQWmJqDJYtmwZDAwMcOrUKTRv3hxLly5FXl6e1LGISAW9UzG6cOECli9fjtWrV+PSpUvvHGbbtm0YP348ZsyYgZiYGHTs2BE9e/ZEQkLCG5cZOHAg/vjjDwQFBeHatWsICQnBe++9985ZiEh9jBw5EpcvX4arqytevHiBiRMnokOHDoiPj5c6GhGpmFIfSgsICICvry+qV68OuVyOx48fw87ODhs3boSTk1Opwjg7O8PR0RGrV69WjNnY2KBfv37w9/cvNP/BgwcxaNAg3L59G8bGxqV6TR5KI6o6RFFEUFAQJk6ciLS0NGhra2P27Nnw8/ODXC6XOh4RlYBKHEoLDg7GhQsXkJ2djfnz52PBggVISUlBcnIy7t27h759+6JLly44depUiYO8fPkS0dHRcHNzUxp3c3NDZGRkkcv89ttvaNWqFRYtWoR69erBysoKfn5+yMrKeuPrZGdnIy0tTWkioqpBEASMGDECV65cQc+ePZGdnY1p06ahXbt2uHLlitTxiEgFlKgYff/993B2doaBgQFSUlJw/vx5LFu2DMePH0e1atUwb948LFq0CH5+fiUO8vjxY+Tl5cHU1FRp3NTUFElJSUUuc/v2bZw6dQpXrlzB7t27ERAQgF9//RWjRo164+v4+/vDyMhIMZmbm5c4KxFVbubm5ti/fz82btyI6tWrIyoqCo6Ojvjuu++Qk5MjdTwiklCJilF8fDzS09MRGRkJTU1NaGhoYPv27ejduzdMTExgYWGBHTt2ICYmBr///jvu3LlT4kCCICg9FkWx0Nhr+fn5EAQBmzdvRps2bdCrVy8sXboUGzdufONeo2nTpiE1NVUx3b9/v8QZiajyEwQBQ4cOxdWrV9GnTx/k5ORg5syZaNOmDWJjY6WOR0QSKfHJ1zo6OmjdujXat28PBwcHnDlzBunp6bh06RL8/f1hZWWFnJwceHt7w9LSstjH/WrWrAmZTFZo71BycnKhvUiv1a1bF/Xq1YORkZFizMbGBqIo4sGDB0Uuo62tDUNDQ6WJiKouMzMz7NmzB5s3b4axsTFiY2PRunVrzJw5Ey9fvpQ6HhFVsFJ/K23JkiVYtGgRRowYgQsXLsDKygoeHh4wMDCAmZkZUlJSkJCQgO3btxdrfVpaWoqbQ/5TWFgYXFxcilymffv2+Ouvv5Qu2Hb9+nVoaGigfv36pX1rRFTFCIKAwYMHIy4uDv3790dubi6+++47ODk5ISoqSup4RFSRxHdw8+ZN0dXVVZTJZKKGhoaooaEhamlpib/88kup1rd161ZRU1NTDAoKEuPi4sTx48eL+vr64t27d0VRFMWpU6eKQ4YMUcyfnp4u1q9fXxwwYIB49epV8cSJE2LTpk3FESNGFPs1U1NTRQBiampqqTITkfrZvn27WKtWLRGAKJPJxKlTp4pZWVlSxyKifyivz+8yufL1o0ePcObMGbx8+RJt27Z9pxOaAwMDsWjRIiQmJsLe3h7Lli1Dp06dAADe3t64e/cujh8/rpj/zz//xJgxYxAREQETExMMHDgQ8+bNg66ubrFej1/XJ6Ki/P333xg7diy2bt0K4NVh+uDgYLRt21biZEQElN/nN28JwmJERG+xZ88e+Pj44NGjR9DQ0MCECRPw3XffFfuXLyIqHypxHSMioqqmX79+iIuLw5AhQ5Cfn48lS5bAwcGhVNdrIyLVx2JERPQfjI2NsWnTJvz+++8wMzPDjRs30KlTJ4wbNw6ZmZlSxyOiMsRiRERUTB988AGuXr2Kzz//HKIo4ocffkDz5s2VznskosqNxYiIqASqV6+OoKAgHDx4EObm5rh9+za6du2Kr7/+Gunp6VLHI6J3xGJERFQK7u7uuHLlCr788ksAwOrVq9GsWbNC12IjosqFxYiIqJQMDQ2xZs0a/PHHH2jYsCHu3bsHNzc3jBw5EqmpqVLHI6JSYDEiInpH77//Pi5fvozRo0cDANavXw87OzuEhoZKnIyISorFiIioDBgYGGDFihU4ceIELC0t8fDhQ/Tu3Rve3t54+vSp1PGIqJhYjIiIylCnTp1w6dIlTJgwAYIg4KeffoKdnR1+++03qaMRUTGwGBERlTE9PT0sXboUp06dgrW1NRITE9G3b198+umnSElJkToeEb0FixERUTlxcXFBTEwMJk+eDA0NDWzZsgW2trbYuXOn1NGI6A1YjIiIypGuri4WLlyIM2fOwM7ODsnJyRgwYAA+/vhjJCcnSx2PiP6FxYiIqAK0bt0a0dHR+OabbyCTyfDrr7/C1tYWISEhqOL38iZSKSxGREQVRFtbG9999x3OnTuH5s2bIyUlBYMHD8ZHH32EpKQkqeMREViMiIgqnKOjI86fP4/Zs2dDLpdjz549sLW1xc8//8y9R0QSYzEiIpKAlpYWZs2ahejoaDg6OuLp06fw8vKCh4cHHj58KHU8oiqLxYiISELNmzfH2bNnMX/+fGhpaWH//v2wtbVFUFAQ9x4RSYDFiIhIYnK5HNOmTUNMTAycnZ2RlpaGESNGoEePHkhISJA6HlGVwmJERKQibG1tERERge+//x7a2to4fPgw7OzssGbNGuTn50sdj6hKYDEiIlIhMpkMfn5+uHjxIlxcXJCRkYGvvvoKrq6uuHPnjtTxiNQeixERkQqytrbGyZMnERAQAF1dXRw7dgzNmjXDypUrufeIqByxGBERqSiZTIZx48bh8uXL6Ny5MzIzMzFmzBh06dIFN27ckDoekVpiMSIiUnGWlpY4evQoVq1aBX19fYSHh8PBwQFLly5FXl6e1PGI1AqLERFRJaChoYGvv/4aV65cgaurK7KysjBx4kR07NgRf/75p9TxiNQGixERUSXSsGFDHD58GD/++COqVauG06dPo0WLFli4cCFyc3OljkdU6bEYERFVMoIgYOTIkbh69Sp69OiB7OxsTJ06FS4uLrhy5YrU8YgqNRYjIqJKytzcHKGhodiwYQOMjIxw/vx5ODo6Yt68ecjJyZE6HlGlxGJERFSJCYIAb29vxMXFwcPDAzk5Ofj222/h7OyMixcvSh2PqNJhMSIiUgNmZmbYu3cvfvnlFxgbGyMmJgatWrXCrFmz8PLlS6njEVUaLEZERGpCEAR8+umnuHr1Kj766CPk5uZi7ty5aNWqFaKjo6WOR1QpsBgREamZOnXq4Ndff8W2bdtQs2ZNXL58Gc7Ozpg+fTpevHghdTwilcZiRESkhgRBwMCBAxEXFwdPT0/k5eXB398fjo6OOHv2rNTxiFQWixERkRqrVasWtm7dil27dsHU1BTx8fFwcXHBpEmTkJWVJXU8IpXDYkREVAV8+OGHuHr1Kj777DPk5+dj8eLFaNGiBSIiIqSORqRSWIyIiKoIExMT/Pzzz/j9999hZmaG69evo2PHjhg/fjwyMzOljkekEliMiIiqmA8++ABXr17FsGHDIIoili9fjubNm+P48eNSRyOSHIsREVEVVL16dQQHB+PAgQOoX78+bt++ja5du2LUqFFIT0+XOh6RZFiMiIiqsB49euDq1av44osvAACBgYFo1qwZjhw5InEyImmwGBERVXGGhoZYu3Ytjhw5goYNG+LevXvo3r07vvjiC6Smpkodj6hCsRgREREAoFu3brh8+TJGjRoFAFi3bh3s7e1x8OBBiZMRVRwWIyIiUjAwMMDKlStx/PhxWFpa4sGDB+jZsyeGDRuGp0+fSh2PqNyxGBERUSGdO3fGpUuXMH78eAiCgI0bN8LOzg6///671NGIyhWLERERFUlPTw/Lli1DeHg4rKyskJiYiD59+uDTTz9FSkqK1PGIygWLERERvVX79u0RGxuLSZMmQUNDA1u2bIGtrS127twpdTSiMsdiRERE/0lXVxeLFi1CZGQkbG1tkZycjAEDBmDgwIFITk6WOh5RmWExIiKiYnN2dsaFCxcwffp0yGQy7NixA3Z2dti2bRtEUZQ6HtE7YzEiIqIS0dbWxv/+9z+cO3cOzZs3x+PHjzFo0CD0798fSUlJRS6TmZkJQRAgCALvy0YqjcWIiIhKxdHREefPn8fs2bMhl8uxe/du2Nra4ueff+beI6q0WIyIiKjUtLS0MGvWLERFRaFly5Z4+vQpvLy80KdPHzx8+FDqeEQlxmJERETvzMHBAWfPnsW8efOgpaWFffv2wc7ODqtXr0ZGRobS4bPMzEyliUiVqFwxCgwMRKNGjaCjowMnJyeEh4cXa7mIiAjI5XK0aNGifAMSEVGRNDU1MWPGDFy4cAGtW7dGamoqvv76a1SrVg2mpqaK+UxNTWFgYKCYiFSJShWjbdu2Yfz48ZgxYwZiYmLQsWNH9OzZEwkJCW9dLjU1FV5eXujWrVsFJSUiojexs7NDZGSk1DGISkWlitHSpUsxfPhwjBgxAjY2NggICIC5uTlWr1791uW+/PJLDB48GO3ataugpERE9DZyuRwZGRm4cOEC2rZtq/Tc9evXkZGRoZiIVInKFKOXL18iOjoabm5uSuNubm5v/c1jw4YNuHXrFmbNmlXeEYmIqAT09fXRsmVLnDp1CgsXLlSM9+jRAzdv3oS+vj709fUlTEhUmMoUo8ePHyMvL0/pODTw6lj0m66LcePGDUydOhWbN2+GXC4v1utkZ2cjLS1NaSIiovIjk8kwatQoxePbt2+jXbt2+PnnnyVMRVQ0lSlGrwmCoPRYFMVCYwCQl5eHwYMHY86cObCysir2+v39/WFkZKSYzM3N3zkzEREVX/fu3ZGVlQUvLy+MGjUKL1++lDoSkYLKFKOaNWtCJpMV2juUnJxcaC8SAKSnpyMqKgqjR4+GXC6HXC7H3LlzcfHiRcjlchw9erTI15k2bRpSU1MV0/3798vl/RAR0f/T19eHKIoQRREHDhxQnP4QGBiIzp0748GDBxInJHpFZYqRlpYWnJycEBYWpjQeFhYGFxeXQvMbGhri8uXLiI2NVUw+Pj6wtrZGbGwsnJ2di3wdbW1tGBoaKk1ERFRxZDIZZs+ejX379qF69eo4c+YMHB0dcezYMamjEalOMQIAX19frF+/HsHBwYiPj8eECROQkJAAHx8fAK/29nh5eQEANDQ0YG9vrzTVrl0bOjo6sLe35wl9REQqrnfv3oiOjoaDgwP+/vtvuLq64vvvv+ftREhSKlWMPD09ERAQgLlz56JFixY4efIkQkNDYWFhAQBITEz8z2saERFR5dG4cWNERkbCy8sL+fn5mDx5MgYMGMAvxpBkBLGKV/O0tDQYGRkhNTWVh9WIiCQiiiLWrl2LsWPHIicnB9bW1ti1axdsbW2ljkYqqrw+v1VqjxEREVVNgiDAx8cH4eHhqF+/Pq5du4Y2bdpg+/btUkejKobFiIiIVIazszOio6Px/vvvIzMzE56envD19UVOTo7U0aiKYDEiIiKVUrt2bRw6dAhTp04FACxbtgzdunV748V+icoSixEREakcuVwOf39/7Nq1C9WqVUN4eDgcHR0REREhdTRScyxGRESksj788ENERUXB1tYWiYmJ6NKlC3744Qd+pZ/KDYsRERGpNCsrK5w9exaenp7Izc3FuHHj8OmnnyIzM1PqaKSGWIyIiEjlGRgYICQkBAEBAZDL5QgJCUHbtm1x/fp1qaORmmExIiKiSkEQBIwbNw5Hjx5FnTp1cOXKFbRu3Rp79+6VOhqpERYjIiKqVDp27IgLFy6gQ4cOSEtLQ79+/TB9+nTk5eVJHY3UAIsRERFVOnXr1sXRo0cxYcIEAIC/vz969OiBv//+W+JkVNmxGBERUaWkqamJpUuXIiQkBHp6ejhy5AicnJxw7tw5qaNRJcZiREREldqgQYNw7tw5WFlZ4f79++jYsSPWrl3Lr/RTqbAYERFRpWdnZ4fz58/jww8/xMuXL+Hj44PPP/8cWVlZUkejSobFiIiI1IKhoSF27tyJhQsXQkNDAxs3boSLiwvu3LkjdTSqRFiMiIhIbQiCgMmTJyMsLAy1atVCbGwsnJyccODAAamjUSXBYkRERGrn/fffR3R0NJydnfH06VP07t0bc+bMQX5+vtTRSMWxGBERkVoyNzfHiRMn8NVXX0EURcyePRseHh548uSJ1NFIhbEYERGR2tLW1kZgYCA2btwIHR0dhIaGolWrVoiJiZE6GqkoFiMiIlJ7Q4cOxenTp9GoUSPcuXMHLi4u2Lhxo9SxSAWxGBERUZXQokULREdHo3fv3njx4gWGDRsGHx8fZGdnSx2NVAiLERERVRk1atTAb7/9hrlz50IQBKxduxadOnXC/fv3pY5GKoLFiIiIqhQNDQ18++23CA0NRY0aNXDu3Dk4Ojrijz/+kDoaqQAWIyIiqpJ69OiB6OhotGzZEo8fP4abmxsWLFjAW4lUcSxGRERUZTVq1AgREREYNmwY8vPzMW3aNHz00UdITU2VOhpJhMWIiIiqNF1dXQQFBeHHH3+ElpYW9uzZg9atW+PKlStSRyMJsBgREVGVJwgCRo4ciVOnTsHc3Bw3btyAs7MzQkJCpI5GFYzFiIiIqEDr1q1x4cIFdO/eHc+fP8fgwYMxfvx45OTkSB2NKgiLERER0T/UrFkTBw4cwIwZMwAAy5cvR9euXZGYmChxMqoILEZERET/IpPJMG/ePOzduxeGhoaIiIhAy5YtcfLkSamjUTljMSIiInqDPn36ICoqCvb29nj06BHef/99LFu2jF/pV2MsRkRERG/RtGlTnDlzBoMHD0ZeXh58fX0xaNAgZGRkSB2NygGLERER0X/Q19fHL7/8ghUrVkAul2P79u1o06YNrl27JnU0KmMsRkRERMUgCAJGjx6NEydOwMzMDPHx8WjdujV27doldTQqQyxGREREJeDi4oLo6Gh07twZ6enp6N+/P6ZMmYLc3Fypo1EZYDEiIiIqoTp16uDIkSPw8/MDACxatAhubm5ITk6WOBm9KxYjIiKiUpDL5fj++++xY8cOGBgY4NixY3B0dMSZM2ekjkbvgMWIiIjoHQwYMADnzp3De++9h4cPH6JTp04IDAzkV/orKRYjIiKid2RjY4Nz585hwIAByMnJwahRo+Dt7Y3nz59LHY1KiMWIiIioDFSrVg3bt2/H4sWLIZPJsGnTJrRr1w63bt2SOhqVAIsRERFRGREEARMnTsSRI0dQu3ZtXLp0CU5OTti3b5/U0aiYWIyIiIjKWJcuXXDhwgW0a9cOqamp8PDwwMyZM5GXlyd1NPoPLEZERETloF69ejh+/DhGjx4NAPjuu+/Qu3dvpKSkSJyM3obFiIiIqJxoaWlhxYoV+Pnnn6Grq4tDhw7ByckJ0dHRUkejN2AxIiIiKmefffYZzpw5A0tLS9y7dw/t27dHcHCw1LGoCCxGREREFaB58+aIioqCh4cHsrOzMXz4cHzxxRd48eKF1NHoH1iMiIiIKkj16tWxZ88e/O9//4MgCFi3bh06dOiAe/fuSR2NCrAYERERVSANDQ1Mnz4dBw8ehImJCaKjo+Ho6IjDhw9LHY3AYkRERCQJNzc3REdHo1WrVnjy5Al69OiB//3vf8jPz5c6WpXGYkRERCQRCwsLhIeHY+TIkRBFEd988w369euHZ8+eSR2tymIxIiIikpCOjg5+/PFHrF+/Htra2vj999/RqlUrXLp0SepoVRKLERERkQoYPnw4IiIiYGFhgVu3bqFt27b45ZdfpI5V5ahcMQoMDESjRo2go6MDJycnhIeHv3HeXbt2oXv37qhVqxYMDQ3Rrl07HDp0qALTEhERlZ3XF390d3dHVlYWhgwZgjFjxuDly5dSR6syVKoYbdu2DePHj8eMGTMQExODjh07omfPnkhISChy/pMnT6J79+4IDQ1FdHQ0unbtCg8PD8TExFRwciIiorJhYmKC/fv3Y+bMmQCAlStXokuXLnj48KHEyaoGQRRFUeoQrzk7O8PR0RGrV69WjNnY2KBfv37w9/cv1jrs7Ozg6emp+Av1X9LS0mBkZITU1FQYGhqWKjcREVF52LdvH4YMGYJnz56hdu3a2LZtG7p06SJ1LJVQXp/fKrPH6OXLl4iOjoabm5vSuJubGyIjI4u1jvz8fKSnp8PY2PiN82RnZyMtLU1pIiIiUkUffPABoqKi4ODggOTkZLi6umLx4sVQoX0aakdlitHjx4+Rl5cHU1NTpXFTU1MkJSUVax1LlixBZmYmBg4c+MZ5/P39YWRkpJjMzc3fKTcREVF5srS0RGRkJIYMGYK8vDxMmjQJAwcORHp6utTR1JLKFKPXBEFQeiyKYqGxooSEhGD27NnYtm0bateu/cb5pk2bhtTUVMV0//79d85MRERUnvT09PDTTz8hMDAQmpqa+PXXX9GmTRvEx8dLHU3tqEwxqlmzJmQyWaG9Q8nJyYX2Iv3btm3bMHz4cGzfvh2urq5vnVdbWxuGhoZKExERkaoTBAFfffUVTp48iXr16uHPP/9EmzZtsGPHDqmjqRWVKUZaWlpwcnJCWFiY0nhYWBhcXFzeuFxISAi8vb2xZcsW9O7du7xjEhERSapt27a4cOECunbtioyMDAwcOBB+fn7Izc2VOppaUJliBAC+vr5Yv349goODER8fjwkTJiAhIQE+Pj4AXh0G8/LyUswfEhICLy8vLFmyBG3btkVSUhKSkpKQmpoq1VsgIiIqd7Vr18bhw4cxZcoUAK/OsXV1dcWjR48kTlb5qVQx8vT0REBAAObOnYsWLVrg5MmTCA0NhYWFBQAgMTFR6ZpGa9euRW5uLkaNGoW6desqpnHjxkn1FoiIiCqEXC7HggULsHPnTlSrVg0nTpyAo6Njsb/JTUVTqesYSYHXMSIiosru2rVr+OijjxAXFwe5XI5ly5Zh1KhRxfryUmWl9tcxIiIiotKxtrbG2bNn4enpidzcXIwZMwZDhgxBZmam1NEqHe4xKmbjzMvLQ05OTgUmIypbWlpa0NDg70JE6kwURSxfvhx+fn7Iy8uDvb09du3ahaZNm0odrcyV1x4jFqP/+MGKooikpCQ8e/as4sMRlSENDQ00atQIWlpaUkchonIWHh6OgQMHIikpCYaGhti0aRP69u0rdawyxWJUTv7rB5uYmKi4R42enp5aH68l9ZWfn4+//voLmpqaaNCgAf8eE1UBiYmJGDhwIE6dOgUAmD59OubOnYsXL17AwMAAAJCRkQF9fX0pY5ZaeRUjeZmtSQ3l5eUpSpGJiYnUcYjeSa1atfDXX38hNzcXmpqaUschonJWt25dHD16FJMmTcLy5csxf/58nD9/HuvWrZM6mkrjCQdv8fqcIj09PYmTEL2714fQ8vLyJE5CRBVFU1MTAQEB2LJlC/T09BAWFoaOHTtKHUulsRgVAw87kDrg32OiqqtPnz44duwYmjRponSP0IyMDGRmZiom4qE0IiIitff6nKJ/q1OnjtLjKn7aMQDuMSIV1LBhQwQEBJTZ+rp06YLx48eX2fr+razylndOIiL6byxGasrb2xuCIGDBggVK43v27FH5Qyrnz5/HF198IXWMcnP8+HEIglDoEhC7du3Cd999J00oIlJrGRkZiumf91N79OiR0nPEYqTWdHR0sHDhQjx9+lTqKMXy8uVLAK++PVUVT3g3NjZGtWrVpI5BRGpIX19fafqv8aqMxaiERFFUOlGtIqeSHvt1dXVFnTp14O/v/8Z5Zs+ejRYtWiiNBQQEoGHDhorH3t7e6NevH+bPnw9TU1NUr14dc+bMQW5uLiZNmgRjY2PUr18fwcHBSut5+PAhPD09UaNGDZiYmKBv3764e/duofX6+/vDzMwMVlZWAAofmnr27Bm++OILmJqaQkdHB/b29ti3bx8AICUlBZ988gnq168PPT09NGvWDCEhISX6OV28eBFdu3ZFtWrVYGhoCCcnJ0RFRSme37lzJ+zs7KCtrY2GDRtiyZIlb1zX3bt3IQgCYmNjlfILgoDjx4/j7t276Nq1KwCgRo0aEAQB3t7eAAofSnv69Cm8vLxQo0YN6OnpoWfPnrhx44bi+Y0bN6J69eo4dOgQbGxsYGBggB49eiAxMbFE75+IiP4fT74uoefPn7/xJLbyVtILcclkMsyfPx+DBw/G2LFjUb9+/VK/9tGjR1G/fn2cPHkSERERGD58OE6fPo1OnTrh7Nmz2LZtG3x8fNC9e3eYm5vj+fPn6Nq1Kzp27IiTJ09CLpdj3rx56NGjBy5duqT46vgff/wBQ0NDhIWFFVn88vPz0bNnT6Snp+OXX36BpaUl4uLiIJPJAAAvXryAk5MTpkyZAkNDQ+zfvx9DhgxB48aN4ezsXKz39umnn6Jly5ZYvXo1ZDIZYmNjFdf5iY6OxsCBAzF79mx4enoiMjISX3/9NUxMTBSFpiTMzc2xc+dO9O/fH9euXYOhoSF0dXWLnNfb2xs3btzAb7/9BkNDQ0yZMgW9evVCXFycIt/z58+xePFi/Pzzz9DQ0MBnn30GPz8/bN68ucTZiIiIxUjtffjhh2jRogVmzZqFoKCgUq/H2NgYP/zwAzQ0NGBtbY1Fixbh+fPnmD59OgBg2rRpWLBgASIiIjBo0CBs3boVGhoaWL9+veKcpg0bNqB69eo4fvw43NzcALzajbt+/fo33qbiyJEjOHfuHOLj4xV7lBo3bqx4vl69evDz81M8HjNmDA4ePIgdO3YUuxglJCRg0qRJeO+99wBA6Z5CS5cuRbdu3fDtt98CAKysrBAXF4fvv/++VMVIJpPB2NgYAFC7dm1Ur169yPleF6KIiAi4uLgAADZv3gxzc3Ps2bMHH3/8MYBX19pas2YNLC0tAQCjR4/G3LlzS5yLiKoOfX19fvvsLViMSkhPT0+yE9RKe97NwoUL8f7772PixImlfm07OzulG5CamprC3t5e8Vgmk8HExATJyckAXu1puXnzZqFzZl68eIFbt24pHjdr1uyt9+6KjY1F/fr1FaXo3/Ly8rBgwQJs27YNDx8+RHZ2NrKzs0u0Z83X1xcjRozAzz//DFdXV3z88ceKohEfH1/o/kLt27dHQEAA8vLyFHuuylp8fDzkcrlSuTMxMYG1tTXi4+MVY3p6eoqswKsr3b7eBkREVHIsRiUkCEKlO0GtU6dOcHd3x/Tp0wvt5dDQ0Cj0m8PrK37/079vISEIQpFj+fn5AF4dAnNycirykE6tWrUUf/6vn+WbDjO9tmTJEixbtgwBAQFo1qwZ9PX1MX78eMWJ3MUxe/ZsDB48GPv378eBAwcwa9YsbN26FR9++CFEUSz0Lb63/ab1ujz+c56ifp7/5U2v8e88RW0D/iZIRFR6PPm6iliwYAF+//13REZGKo3XqlULSUlJSh+m/zxxuLQcHR1x48YN1K5dG02aNFGajIyMir2e5s2b48GDB7h+/XqRz4eHh6Nv37747LPP4ODggMaNGyudoFxcVlZWmDBhAg4fPoyPPvoIGzZsAADY2toqbsD4WmRkJKysrIrcW/S69P3zBOh//zyLc2sOW1tb5Obm4uzZs4qxlJQUXL9+HTY2NiV7c0REVGwsRlVEs2bN8Omnn2LFihVK4126dMHff/+NRYsW4datW1i1ahUOHDjwzq/36aefombNmujbty/Cw8Nx584dnDhxAuPGjcODBw+KvZ7OnTujU6dO6N+/P8LCwnDnzh0cOHAABw8eBAA0adIEYWFhiIyMRHx8PL788kskJSUVe/1ZWVkYPXo0jh8/jnv37iEiIgLnz59XlI+JEyfijz/+wHfffYfr16/jp59+wsqVK5XOa/onXV1dtG3bFgsWLEBcXBxOnjyJb775RmkeCwsLCIKAffv24e+//y7y0GzTpk3Rt29fjBw5EqdOncLFixfx2WefoV69eoUO7RERUdlhMapCvvvuu0KHWWxsbBAYGIhVq1bBwcEB586de+OHfkno6enh5MmTaNCgAT766CPY2Njg888/R1ZWFgwNDUu0rp07d6J169b45JNPYGtri8mTJyv2tnz77bdwdHSEu7s7unTpgjp16qBfv37FXrdMJkNKSgq8vLxgZWWFgQMHomfPnpgzZw6AV3u+tm/fjq1bt8Le3h4zZ87E3Llz33ridXBwMHJyctCqVSuMGzcO8+bNU3q+Xr16mDNnDqZOnQpTU1OMHj26yPVs2LABTk5O+OCDD9CuXTuIoojQ0NBCh8+IiKjsCGIVPyEhLS0NRkZGSE1NLfSB/eLFC9y5cweNGjWCjo6ORAmJygb/PhOROnnb5/e74B4jIiIiogIsRkREREQFWIyIiIiICrAYERERERVgMSIiIiIqwGJEREREVIDFiIiIiKgAixERERFRARYjIiIiogIsRvTOGjZsiICAgHdax/HjxyEIAp49e1Ymme7evQtBEMrkhrhFKau85Z2TiIhKhsVIzUVGRkImk6FHjx5SR1Ho0qULxo8frzTm4uKCxMREGBkZSROqAnh7exe6j5u5uTkSExNhb28vTSgiIlLCYlRBMjMzIQgCBEFAZmZmhb1ucHAwxowZg1OnTiEhIaHCXrektLS0UKdOHQiCIHWUCiWTyVCnTh3I5XKpoxAREViM1FpmZia2b9+Or776Ch988AE2btyo9Pzrw0F//PEHWrVqBT09Pbi4uODatWuKeW7duoW+ffvC1NQUBgYGaN26NY4cOfLG1/z888/xwQcfKI3l5uaiTp06CA4Ohre3N06cOIHly5criuLdu3eLPDQVERGBzp07Q09PDzVq1IC7uzuePn0KADh48CA6dOiA6tWrw8TEBB988AFu3bpVop9PYGAgmjZtCh0dHZiammLAgAGK57KzszF27FjUrl0bOjo66NChA86fP//Gdc2ePRstWrRQGgsICEDDhg0Vz//000/Yu3ev4n0fP368yENpJ06cQJs2baCtrY26deti6tSpyM3NVTzfpUsXjB07FpMnT4axsTHq1KmD2bNnl+i9ExFR0ViM1Ni2bdtgbW0Na2trfPbZZ9iwYQNEUSw034wZM7BkyRJERUVBLpfj888/VzyXkZGBXr164ciRI4iJiYG7uzs8PDzeuPdpxIgROHjwIBITExVjoaGhyMjIwMCBA7F8+XK0a9cOI0eORGJiIhITE2Fubl5oPbGxsejWrRvs7Oxw+vRpnDp1Ch4eHsjLywPwqvT5+vri/Pnz+OOPP6ChoYEPP/wQ+fn5xfrZREVFYezYsZg7dy6uXbuGgwcPolOnTornJ0+ejJ07d+Knn37ChQsX0KRJE7i7u+PJkyfFWv+/+fn5YeDAgejRo4fifbu4uBSa7+HDh+jVqxdat26NixcvYvXq1QgKCsK8efOU5vvpp5+gr6+Ps2fPYtGiRZg7dy7CwsJKlY2IiP5BrOJSU1NFAGJqamqh57KyssS4uDgxKyurVOvOyMhQTI8ePRIBiADER48eKT1XXlxcXMSAgABRFEUxJydHrFmzphgWFqZ4/tixYyIA8ciRI4qx/fv3iwDe+p5tbW3FFStWKB5bWFiIy5YtU3p+4cKFisf9+vUTvb29FY87d+4sjhs3Tmmdr7M8ffpUFEVR/OSTT8T27dsX+70mJyeLAMTLly+LoiiKd+7cEQGIMTExRc6/c+dO0dDQUExLSyv0XEZGhqipqSlu3rxZMfby5UvRzMxMXLRoUZF5Z82aJTo4OCitZ9myZaKFhYXi8dChQ8W+ffsqzfPvnNOnTxetra3F/Px8xTyrVq0SDQwMxLy8PFEUX/38OnTooLSe1q1bi1OmTCnyvb72rn+fiYhUyds+v98F9xiVIwMDA8VkamqqGH99WOr1VB6uXbuGc+fOYdCgQQAAuVwOT09PBAcHF5q3efPmij/XrVsXAJCcnAzg1Z6ZyZMnw9bWFtWrV4eBgQH+/PPPt56vNGLECGzYsEGxnv379yvthSqO13uM3uTWrVsYPHgwGjduDENDQzRq1AgAin0eVffu3WFhYYHGjRtjyJAh2Lx5M54/f65Yd05ODtq3b6+YX1NTE23atEF8fHyJ3kdJxcfHo127dkrnWrVv3x4ZGRl48OCBYuyf2wx4td1ebzMiIio9nvGppoKCgpCbm4t69eopxkRRhKamJp4+fYoaNWooxjU1NRV/fv2B/PqQ1KRJk3Do0CEsXrwYTZo0ga6uLgYMGICXL1++8bW9vLwwdepUnD59GqdPn0bDhg3RsWPHEuXX1dV96/MeHh4wNzfHunXrYGZmhvz8fNjb27811z9Vq1YNFy5cwPHjx3H48GHMnDkTs2fPxvnz5xWHG/99Irgoim88OVxDQ6PQYcqcnJxiZfmv1ygqzz+32evninsYkYiI3ox7jMpRRkaGYnr06JFi/NGjR0rPlbXc3Fxs2rQJS5YsQWxsrGK6ePEiLCwssHnz5mKvKzw8HN7e3vjwww/RrFkz1KlTB3fv3n3rMiYmJujXrx82bNiADRs2YNiwYUrPa2lpKc4VepPmzZvjjz/+KPK5lJQUxMfH45tvvkG3bt1gY2OjOCm7JORyOVxdXbFo0SJcunQJd+/exdGjR9GkSRNoaWnh1KlTinlzcnIQFRUFGxubItdVq1YtJCUlKZWjf1+bqDjv29bWFpGRkUrriYyMRLVq1ZRKLhERlQ/uMSpH+vr6bxx/03NlYd++fXj69CmGDx9e6LpAAwYMQFBQEEaPHl2sdTVp0gS7du2Ch4cHBEHAt99+W6w9EyNGjMAHH3yAvLw8DB06VOm5hg0b4uzZs7h79y4MDAxgbGxcaPlp06ahWbNm+Prrr+Hj4wMtLS0cO3YMH3/8MYyNjWFiYoIff/wRdevWRUJCAqZOnVqs9/Pavn37cPv2bXTq1Ak1atRAaGgo8vPzYW1tDX19fXz11VeYNGkSjI2N0aBBAyxatAjPnz/H8OHDi1xfly5d8Pfff2PRokUYMGAADh48iAMHDsDQ0FDpfR86dAjXrl2DiYlJkdds+vrrrxEQEIAxY8Zg9OjRuHbtGmbNmgVfX19oaPD3GCKi8sb/adVQUFAQXF1di/zg7d+/P2JjY3HhwoVirWvZsmWoUaMGXFxc4OHhAXd3dzg6Ov7ncq6urqhbty7c3d1hZmam9Jyfnx9kMhlsbW1Rq1atIs8LsrKywuHDh3Hx4kW0adMG7dq1w969eyGXy6GhoYGtW7ciOjoa9vb2mDBhAr7//vtivZ/Xqlevjl27duH999+HjY0N1qxZg5CQENjZ2QEAFixYgP79+2PIkCFwdHTEzZs3cejQIaVDkP9kY2ODwMBArFq1Cg4ODjh37hz8/PyU5hk5ciSsra3RqlUr1KpVCxEREYXWU69ePYSGhuLcuXNwcHCAj48Phg8fjm+++aZE74+IiEpHEP99YkQVk5aWBiMjI6Smpir9dg8AL168wJ07d9CoUSPo6Oi80+tkZmYqTrTOyMgo1z1GquD58+cwMzNDcHAwPvroI6njEMr27zMRkdTe9vn9LngorYLo6+sXeQ0hdZOfn4+kpCQsWbIERkZG6NOnj9SRiIiIio3FiMpUQkICGjVqhPr162Pjxo281QUREVUq/NSiMtWwYcMqsWeMiIjUE0++JiIiIirAYlQM3ANC6oB/j4mI/huL0Vu8vrrw61tFEFVmr68KLpPJJE5CRKS6eI7RW8hkMlSvXl1xDyo9Pb033hKCSJXl5+fj77//hp6eHk+IJyJ6C/4P+R/q1KkDALxBJ1V6GhoaaNCgAcs9EdFbsBj9B0EQULduXdSuXbtUNwUlUhVaWlq8rQgR0X9gMSommUzGczOIiIjUnMr9+hgYGKi4ZYGTkxPCw8PfOv+JEyfg5OQEHR0dNG7cGGvWrKmgpERERKRuVKoYbdu2DePHj8eMGTMQExODjh07omfPnkXeZBQA7ty5g169eqFjx46IiYnB9OnTMXbsWOzcubOCkxMREZE6UKmbyDo7O8PR0RGrV69WjNnY2KBfv37w9/cvNP+UKVPw22+/IT4+XjHm4+ODixcv4vTp08V6zfK6CR0RERGVH7W/iezLly8RHR2NqVOnKo27ubkhMjKyyGVOnz4NNzc3pTF3d3cEBQUhJydHcR2if8rOzkZ2drbicWpqKoBXP2AiIiKqHF5/bpf1/h2VKUaPHz9GXl4eTE1NlcZNTU2RlJRU5DJJSUlFzp+bm4vHjx+jbt26hZbx9/fHnDlzCo2bm5u/Q3oiIiKSQkpKCoyMjMpsfSpTjF779zVWRFF863VXipq/qPHXpk2bBl9fX8XjZ8+ewcLCAgkJCWX6g6XSSUtLg7m5Oe7fv89DmxLjtlAd3Baqg9tCdaSmpqJBgwYwNjYu0/WqTDGqWbMmZDJZob1DycnJhfYKvVanTp0i55fL5TAxMSlyGW1tbWhraxcaNzIy4l9yFWJoaMjtoSK4LVQHt4Xq4LZQHWV9fTaV+VaalpYWnJycEBYWpjQeFhYGFxeXIpdp165dofkPHz6MVq1aFXl+EREREdHbqEwxAgBfX1+sX78ewcHBiI+Px4QJE5CQkAAfHx8Arw6DeXl5Keb38fHBvXv34Ovri/j4eAQHByMoKAh+fn5SvQUiIiKqxFTmUBoAeHp6IiUlBXPnzkViYiLs7e0RGhoKCwsLAEBiYqLSNY0aNWqE0NBQTJgwAatWrYKZmRl++OEH9O/fv9ivqa2tjVmzZhV5eI0qHreH6uC2UB3cFqqD20J1lNe2UKnrGBERERFJSaUOpRERERFJicWIiIiIqACLEREREVEBFiMiIiKiAlWiGAUGBqJRo0bQ0dGBk5MTwsPD3zr/iRMn4OTkBB0dHTRu3Bhr1qypoKTqryTbYteuXejevTtq1aoFQ0NDtGvXDocOHarAtOqvpP82XouIiIBcLkeLFi3KN2AVUtJtkZ2djRkzZsDCwgLa2tqwtLREcHBwBaVVbyXdFps3b4aDgwP09PRQt25dDBs2DCkpKRWUVn2dPHkSHh4eMDMzgyAI2LNnz38uUyaf36Ka27p1q6ipqSmuW7dOjIuLE8eNGyfq6+uL9+7dK3L+27dvi3p6euK4cePEuLg4cd26daKmpqb466+/VnBy9VPSbTFu3Dhx4cKF4rlz58Tr16+L06ZNEzU1NcULFy5UcHL1VNLt8dqzZ8/Exo0bi25ubqKDg0PFhFVzpdkWffr0EZ2dncWwsDDxzp074tmzZ8WIiIgKTK2eSrotwsPDRQ0NDXH58uXi7du3xfDwcNHOzk7s169fBSdXP6GhoeKMGTPEnTt3igDE3bt3v3X+svr8Vvti1KZNG9HHx0dp7L333hOnTp1a5PyTJ08W33vvPaWxL7/8Umzbtm25ZawqSrotimJrayvOmTOnrKNVSaXdHp6enuI333wjzpo1i8WojJR0Wxw4cEA0MjISU1JSKiJelVLSbfH999+LjRs3Vhr74YcfxPr165dbxqqoOMWorD6/1fpQ2suXLxEdHQ03NzelcTc3N0RGRha5zOnTpwvN7+7ujqioKOTk5JRbVnVXmm3xb/n5+UhPTy/zGwZWRaXdHhs2bMCtW7cwa9as8o5YZZRmW/z2229o1aoVFi1ahHr16sHKygp+fn7IysqqiMhqqzTbwsXFBQ8ePEBoaChEUcSjR4/w66+/onfv3hURmf6hrD6/VerK12Xt8ePHyMvLK3QTWlNT00I3n30tKSmpyPlzc3Px+PFj1K1bt9zyqrPSbIt/W7JkCTIzMzFw4MDyiFillGZ73LhxA1OnTkV4eDjkcrX+r6NClWZb3L59G6dOnYKOjg52796Nx48f4+uvv8aTJ094ntE7KM22cHFxwebNm+Hp6YkXL14gNzcXffr0wYoVKyoiMv1DWX1+q/Ueo9cEQVB6LIpiobH/mr+ocSq5km6L10JCQjB79mxs27YNtWvXLq94VU5xt0deXh4GDx6MOXPmwMrKqqLiVSkl+beRn58PQRCwefNmtGnTBr169cLSpUuxceNG7jUqAyXZFnFxcRg7dixmzpyJ6OhoHDx4EHfu3FHc45MqVll8fqv1r301a9aETCYr1PSTk5MLtcrX6tSpU+T8crkcJiYm5ZZV3ZVmW7y2bds2DB8+HDt27ICrq2t5xqwySro90tPTERUVhZiYGIwePRrAqw9nURQhl8tx+PBhvP/++xWSXd2U5t9G3bp1Ua9ePRgZGSnGbGxsIIoiHjx4gKZNm5ZrZnVVmm3h7++P9u3bY9KkSQCA5s2bQ19fHx07dsS8efN4lKECldXnt1rvMdLS0oKTkxPCwsKUxsPCwuDi4lLkMu3atSs0/+HDh9GqVStoamqWW1Z1V5ptAbzaU+Tt7Y0tW7bwmH0ZKun2MDQ0xOXLlxEbG6uYfHx8YG1tjdjYWDg7O1dUdLVTmn8b7du3x19//YWMjAzF2PXr16GhoYH69euXa151Vppt8fz5c2hoKH+UymQyAP+/t4IqRpl9fpfoVO1K6PVXL4OCgsS4uDhx/Pjxor6+vnj37l1RFEVx6tSp4pAhQxTzv/6634QJE8S4uDgxKCiIX9cvIyXdFlu2bBHlcrm4atUqMTExUTE9e/ZMqregVkq6Pf6N30orOyXdFunp6WL9+vXFAQMGiFevXhVPnDghNm3aVBwxYoRUb0FtlHRbbNiwQZTL5WJgYKB469Yt8dSpU2KrVq3ENm3aSPUW1EZ6eroYExMjxsTEiADEpUuXijExMYpLJ5TX57faFyNRFMVVq1aJFhYWopaWlujo6CieOHFC8dzQoUPFzp07K81//PhxsWXLlqKWlpbYsGFDcfXq1RWcWH2VZFt07txZBFBoGjp0aMUHV1Ml/bfxTyxGZauk2yI+Pl50dXUVdXV1xfr164u+vr7i8+fPKzi1eirptvjhhx9EW1tbUVdXV6xbt6746aefig8ePKjg1Orn2LFjb/0MKK/Pb0EUua+PiIiICFDzc4yIiIiISoLFiIiIiKgAixERERFRARYjIiIiogIsRkREREQFWIyIiIiICrAYERERERVgMSIiIiIqwGJEREREVIDFiIjUSkpKCmrXro27d+++03oGDBiApUuXlk0oIqo0eEsQIlIrfn5+ePr0KYKCgt5pPZcuXULXrl1x584dGBoallE6IlJ13GNERGojKysLQUFBGDFixDuvq3nz5mjYsCE2b95cBsmIqLJgMSIilRUSEgIdHR08fPhQMTZixAg0b94cqampheY/cOAA5HI52rVrpzTepUsXjBkzBuPHj0eNGjVgamqKH3/8EZmZmRg2bBiqVasGS0tLHDhwQGm5Pn36ICQkpHzeHBGpJBYjIlJZgwYNgrW1Nfz9/QEAc+bMwaFDh3DgwAEYGRkVmv/kyZNo1apVkev66aefULNmTZw7dw5jxozBV199hY8//hguLi64cOEC3N3dMWTIEDx//lyxTJs2bXDu3DlkZ2eXzxskIpXDc4yISKXt27cPAwYMwMyZM7F48WKEh4fDzs6uyHn79esHExOTQucXdenSBXl5eQgPDwcA5OXlwcjICB999BE2bdoEAEhKSkLdunVx+vRptG3bFsCr84wcHBxw9+5dWFhYlOO7JCJVIZc6ABHR23zwwQewtbXFnDlzcPjw4TeWIuDVOUY6OjpFPte8eXPFn2UyGUxMTNCsWTPFmKmpKQAgOTlZMaarqwsASnuRiEi98VAaEam0Q4cO4c8//0ReXp6ivLxJzZo18fTp0yKf09TUVHosCILSmCAIAID8/HzF2JMnTwAAtWrVKlV2Iqp8WIyISGVduHABH3/8MdauXQt3d3d8++23b52/ZcuWiIuLK7PXv3LlCurXr4+aNWuW2TqJSLWxGBGRSrp79y569+6NqVOnYsiQIZg7dy527tyJ6OjoNy7j7u6Oq1evvnGvUUmFh4fDzc2tTNZFRJUDixERqZwnT56gZ8+e6NOnD6ZPnw4AcHJygoeHB2bMmPHG5Zo1a4ZWrVph+/bt75zhxYsX2L17N0aOHPnO6yKiyoPfSiMitRIaGgo/Pz9cuXIFGhql/91v1apV2Lt3Lw4fPlyG6YhI1fFbaUSkVnr16oUbN27g4cOHMDc3L/V6NDU1sWLFijJMRkSVAfcYERERERXgOUZEREREBViMiIiIiAqwGBEREREVYDEiIiIiKsBiRERERFSAxYiIiIioAIsRERERUQEWIyIiIqICLEZEREREBViMiIiIiAr8H/hFCKr6dHmaAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Step 4: visualize the final solution\n", "plt.figure() # create a new figure\n", "plt.plot(x, phi, '-k', label='Numerical solution')\n", "plt.scatter(x, phi_true, s=50, c='k', marker='+', label='Analytical solution')\n", "plt.xlabel('$x$ (m)')\n", "plt.ylabel(r'$\\phi$')\n", "ax = plt.gca()\n", "ax.set_xlim(0, 1.0)\n", "ax.set_ylim(0, 1.2)\n", "plt.legend(loc='lower left')\n", "plt.show() # show the figure" ] }, { "cell_type": "markdown", "id": "4e39eea2", "metadata": {}, "source": [ "## Exercise\n", "\n", "1. Determine the distribution of $\\phi$ on a grid discretized into 5 nodes when the velocity $u = 2.5$ m/s.\n", "\n", "2. Determine the distribution of $\\phi$ on a grid discretized into 20 nodes when the velocity $u = 2.5$ m/s." ] } ], "metadata": { "kernelspec": { "display_name": "cfd2025", "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.12.9" } }, "nbformat": 4, "nbformat_minor": 5 }