heat_transfer_v1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 91,
   "id": "961ee96f-e37f-4fea-bdea-b2f225c043bf",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%load_ext memory_profiler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "id": "e5afa264-4c2e-4af3-b5b3-054fb1a7fffa",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# %load_ext autotime"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69a43ae6-9bfd-4504-b25f-9edd0c1d416d",
   "metadata": {},
   "source": [
    "##### v1, 129 seconds to completion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "id": "47e587c4-ca34-46b7-9da9-40641da057fe",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "peak memory: 423.68 MiB, increment: 45.80 MiB\n"
     ]
    }
   ],
   "source": [
    "%%memit\n",
    "import numpy as np\n",
    "\n",
    "# Grid parameters\n",
    "nx, ny, nz = 100, 100, 100  # Number of grid points\n",
    "dx = dy = dz = 1.0          # Grid spacing\n",
    "dt = 1                    # Time step\n",
    "alpha = 0.15             # Thermal diffusivity\n",
    "nt = 1000                   # Number of time steps\n",
    "\n",
    "def second_derivative_z(u_z):\n",
    "    \"\"\"\n",
    "    Compute the second derivative along the z-axis using central differences.\n",
    "    Assumes Neumann boundary conditions (zero derivative at boundaries).\n",
    "    \"\"\"\n",
    "    d2u_dz2 = np.zeros_like(u_z)\n",
    "    d2u_dz2[1:-1] = (u_z[2:] - 2 * u_z[1:-1] + u_z[:-2]) / dz**2\n",
    "    return d2u_dz2\n",
    "\n",
    "def main():\n",
    "    # Initialize the temperature field with initial conditions\n",
    "    u = np.zeros((nx, ny, nz))\n",
    "    u[nx//2,ny//2, 0] = 1.0  # place heat source here \n",
    "\n",
    "    for _ in range(nt):\n",
    "        # Compute second derivative along the z-axis using apply_along_axis\n",
    "        # d2u_dz2 = np.apply_along_axis(second_derivative_z, axis=2, arr=u)\n",
    "        # Compute second derivative along the x-axis\n",
    "        d2u_dz2 = np.zeros_like(u)\n",
    "        # d2u_dx2[1:-1, :, :] = 0\n",
    "        d2u_dz2[:, :, 1:-1] = (u[:, :, 2:] - 2 * u[:, :, 1:-1] + u[:, :, :-2]) / dz**2\n",
    "        \n",
    "        # Compute second derivative along the x-axis\n",
    "        d2u_dx2 = np.zeros_like(u)\n",
    "        # d2u_dx2[1:-1, :, :] = 0\n",
    "        d2u_dx2[1:-1, :, :] = (u[2:, :, :] - 2 * u[1:-1, :, :] + u[:-2, :, :]) / dx**2\n",
    "\n",
    "        # Compute second derivative along the y-axis\n",
    "        d2u_dy2 = np.zeros_like(u)\n",
    "        # d2u_dy2[:, 1:-1, :] = 0\n",
    "\n",
    "        d2u_dy2[:, 1:-1, :] = (u[:, 2:, :] - 2 * u[:, 1:-1, :] + u[:, :-2, :]) / dy**2\n",
    "\n",
    "        # Update the temperature field\n",
    "        u += alpha * dt * (d2u_dx2 + d2u_dy2 + d2u_dz2)\n",
    "\n",
    "    # Optionally, return or save the final temperature field\n",
    "    return u\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    final_temperature = main()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c48f86b4-12eb-426f-8e48-395b90e2f538",
   "metadata": {},
   "source": [
    "##### vectorized using np.roll 122 seconds to completion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "id": "1ba9abba-fe84-4590-9a75-e917de647b08",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# %%memit\n",
    "# %%time\n",
    "\n",
    "# import numpy as np\n",
    "\n",
    "# Grid parameters\n",
    "nx, ny, nz = 100, 100, 100  # Number of grid points\n",
    "dx = dy = dz = 1.0          # Grid spacing\n",
    "dt = 1                      # Time step\n",
    "alpha = 0.15                # Thermal diffusivity\n",
    "nt = 1000                    # Number of time steps\n",
    "\n",
    "def main():\n",
    "    # Initialize the temperature field with initial conditions\n",
    "    u = np.zeros((nx, ny, nz))\n",
    "    u[nx // 2, ny // 2, 0] = 1.0  # place heat source here \n",
    "\n",
    "    for _ in range(nt):\n",
    "        # Compute second derivative along the z-axis\n",
    "        d2u_dz2 = (np.roll(u, -1, axis=2) - 2 * u + np.roll(u, 1, axis=2)) / dz**2\n",
    "        d2u_dz2[:, :, 0] = (u[:, :, 1] - u[:, :, 0]) / dz**2  # Neumann BC at z=0\n",
    "        d2u_dz2[:, :, -1] = (u[:, :, -2] - u[:, :, -1]) / dz**2  # Neumann BC at z=nz-1\n",
    "\n",
    "        # Compute second derivative along the x-axis\n",
    "        d2u_dx2 = (np.roll(u, -1, axis=0) - 2 * u + np.roll(u, 1, axis=0)) / dx**2\n",
    "        d2u_dx2[0, :, :] = (u[1, :, :] - u[0, :, :]) / dx**2  # Neumann BC at x=0\n",
    "        d2u_dx2[-1, :, :] = (u[-2, :, :] - u[-1, :, :]) / dx**2  # Neumann BC at x=nx-1\n",
    "\n",
    "        # Compute second derivative along the y-axis\n",
    "        d2u_dy2 = (np.roll(u, -1, axis=1) - 2 * u + np.roll(u, 1, axis=1)) / dy**2\n",
    "        d2u_dy2[:, 0, :] = (u[:, 1, :] - u[:, 0, :]) / dy**2  # Neumann BC at y=0\n",
    "        d2u_dy2[:, -1, :] = (u[:, -2, :] - u[:, -1, :]) / dy**2  # Neumann BC at y=ny-1\n",
    "\n",
    "        # Update the temperature field\n",
    "        u += alpha * dt * (d2u_dx2 + d2u_dy2 + d2u_dz2)\n",
    "\n",
    "    # Optionally, return or save the final temperature field\n",
    "    return u\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    final_temperature = main()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "id": "9f4f75aa-a8b5-4936-aee1-b958454ef11f",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAewAAAGsCAYAAAAfcpQMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAxGElEQVR4nO3df3RU9Z3/8dcYYBI0MxZsMkkJkGwjYCKKiYtBQLqUsMDhyMq6/kCg2u3XLCCEHBQCbqWsZBApmyIaxCLioQhnT7DFFVmyCkEWXElIhEVEXSNJadIsVmf4ZUKS+/2DZuqYAJnM5MfHeT7O+ZzT+czn3vuZq/Wd9/t+7r02y7IsAQCAbu2arp4AAAC4OgI2AAAGIGADAGAAAjYAAAYgYAMAYAACNgAABiBgAwBgAAI2AAAGIGADAGAAAjYAAAYgYAMAQmbfvn2aPHmy4uPjZbPZ9Nvf/rZDj9fQ0KAnn3xSiYmJioqKUlJSkpYtW6ampibfmO3bt2v8+PG64YYbZLPZVF5e3qFz6igEbABAyJw7d0633HKL1q5d2ynHe+aZZ7Ru3TqtXbtWx48f18qVK/Xss8/queee85vTnXfeqRUrVnTKnDpKj66eAADgu2PChAmaMGHCZb+vr6/Xk08+qd/85jf66quvlJqaqmeeeUZjxoxp1/EOHjyou+++W5MmTZIkDRw4UK+99ppKSkp8Y6ZPny5J+vzzz9t1jO6CDBsA0Gkefvhh/dd//Ze2bt2qI0eO6N5779Xf/u3f6pNPPmnX/kaOHKm3335bH3/8sSTpgw8+0P79+zVx4sRQTrtbIMMGAHSK//3f/9Vrr72m3//+94qPj5ckLViwQLt27dLGjRuVl5cX8D4XLlwoj8ejwYMHKyIiQo2NjVq+fLkeeOCBUE+/yxGwAQCd4vDhw7IsSzfeeKNff11dnfr27SvpUtk6MTHxivuZPXu27xr5tm3btHnzZm3ZskUpKSkqLy9Xdna24uPjNXPmzI75IV2EgA0A6BRNTU2KiIhQaWmpIiIi/L677rrrJEk/+MEPdPz48Svu53vf+57vfz/++ONatGiR7r//fknSzTffrJMnT8rtdhOwAQBoj2HDhqmxsVG1tbUaNWpUq2N69uypwYMHt3mf58+f1zXX+C/HioiI8Lut67uCgA0ACJmzZ8/q008/9X2uqKhQeXm5+vTpoxtvvFHTpk3TjBkz9Mtf/lLDhg3T6dOn9c477+jmm29u10KxyZMna/ny5erfv79SUlJUVlam1atX65FHHvGN+dOf/qTKykr94Q9/kCSdOHFCkuRyueRyuYL8xZ3IAgAgRPbs2WNJatFmzpxpWZZl1dfXWz//+c+tgQMHWj179rRcLpf1d3/3d9aRI0fadTyv12vNmzfP6t+/vxUZGWklJSVZS5Ysserq6nxjNm7c2OqcnnrqqRD84s5jsyzL6qo/FgAAQNtwHzYAAAbosID9wgsvKDExUZGRkUpLS9O7777bUYcCAOA7r0MWnW3btk3Z2dl64YUXdOedd+rFF1/UhAkT9OGHH6p///5X3LapqUl/+MMfFB0dLZvN1hHTAwB0IMuydObMGcXHx7dYwR1KX3/9terr64PeT69evRQZGRmCGXWsDrmGPXz4cN12220qKCjw9Q0ZMkRTpkyR2+2+4ra///3vlZCQEOopAQA6WVVVlfr169ch+/7666+VmJiompqaoPflcrlUUVHR7YN2yDPs+vp6lZaWatGiRX79mZmZOnDgQIvxdXV1qqur831u/vshUhL5NQCYx5L0taTo6OgOO0Z9fb1qampUVVUhh8PR7v14vV4lJCSqvr6+2wfskNcqTp8+rcbGRsXGxvr1x8bGtvqXkNvtltPp9LXmkrmNRqPRaMY2SZ1yWdPhcATdguF2u2Wz2ZSdnX3FccXFxUpLS1NkZKSSkpK0bt26gI/VYRcXvv0PyrKsVv/h5ebmyuPx+FpVVVVHTQkA8J3TEILWPocOHdL69es1dOjQK46rqKjQxIkTNWrUKJWVlWnx4sWaO3euCgsLAzpeyEviN9xwgyIiIlpk07W1tS2ybkmy2+2y2+2hngYAICwEF3Tbu+3Zs2c1bdo0vfTSS3r66aevOHbdunXq37+/8vPzJV1a01VSUqJVq1Zp6tSpbT5myDPsXr16KS0tTUVFRX79RUVFGjFiRKgPBwAIa6HJsL1er1/75tqq1syePVuTJk3Sj3/846vO8ODBg8rMzPTrGz9+vEpKSnTx4sU2/9IOua0rJydH06dPV3p6ujIyMrR+/XpVVlYqKyurIw4HAEBQvn130lNPPaWlS5e2Onbr1q06fPiwDh061KZ919TUtLquq6GhQadPn1ZcXFyb9tMhAfu+++7TF198oWXLlqm6ulqpqanauXOnBgwY0BGHAwCErUYFVxJvlHTpFrRvLkC73KXaqqoqzZs3T7t37w5oVXlr67pa67+SDntb16xZszRr1qyO2j0AAArVNey2rhgvLS1VbW2t0tLSfH2NjY3at2+f1q5dq7q6uhbv+na5XK2u6+rRo4f69u3b5pnyek0AANpo7NixOnr0qF/fww8/rMGDB2vhwoUtgrUkZWRk6I033vDr2717t9LT09WzZ882H5uADQAwWOeuEo+OjlZqaqpf37XXXqu+ffv6+nNzc3Xq1Cm9+uqrkqSsrCytXbtWOTk5+tnPfqaDBw9qw4YNeu211wI6NgEbAGCwrrmt60qqq6tVWVnp+5yYmKidO3dq/vz5ev755xUfH681a9YEdEuXpI55lngwvF6vnE6novSXp+UAAMxhSbogyePxBP0ksctpjhUezxE5HO1/BKrXe0ZO59AOnWuokGEDAAzWqOaV3u3f3gwEbACAwUJzW5cJOu5FpQAAIGTIsAEABut+i846CgEbAGAwAjYAAAYIn4DNNWwAAAxAhg0AMFj4rBInYAMADEZJHAAAdCNk2AAAg4VPhk3ABgAYLHwCNiVxAAAMQIYNADBY+GTYBGwAgMHC57YuSuIAABiADBsAYDBK4gAAGICADQCAAcInYHMNGwAAA5BhAwAMFj4ZNgEbAGAwbusCAADdCBk2AMBgjQouSzYnwyZgAwAMFj7XsCmJAwBgADJsAIDBwifDJmADAAzGKnEAANCNkGEDAAxGSRwAAAMQsAEAMED4BGyuYQMAYAAybACAwciwAQAwQPNtXe1tgd3WVVBQoKFDh8rhcMjhcCgjI0NvvfXWZcfv3btXNputRfvoo48C/aFk2AAAtFW/fv20YsUK/fCHP5Qkbdq0SXfffbfKysqUkpJy2e1OnDghh8Ph+/z9738/4GMTsAEABmuQFBHk9m03efJkv8/Lly9XQUGB3nvvvSsG7JiYGF1//fXtmaAPJXEAgMGCKYf/5fq31+v1a3V1dVc9cmNjo7Zu3apz584pIyPjimOHDRumuLg4jR07Vnv27GnXLyVgAwDCXkJCgpxOp6+53e7Ljj169Kiuu+462e12ZWVl6fXXX9dNN93U6ti4uDitX79ehYWF2r59uwYNGqSxY8dq3759Ac/RZlmWFfBWHcjr9crpdCpKkq2rJwMACJgl6YIkj8fjd902lJpjhcczSw6HPYj91MnpfEFVVVV+c7Xb7bLbW99vfX29Kisr9dVXX6mwsFC//vWvVVxcfNmg/W2TJ0+WzWbTjh07Apor17ABAAYLzcs/mld9t0WvXr18i87S09N16NAh/epXv9KLL77Ypu3vuOMObd68OeCZUhIHACAIlmW16Zp3s7KyMsXFxQV8HDJsAIDBGhRc7hlYdr548WJNmDBBCQkJOnPmjLZu3aq9e/dq165dkqTc3FydOnVKr776qiQpPz9fAwcOVEpKiurr67V582YVFhaqsLAw4JkSsAEABuvcgP3HP/5R06dPV3V1tZxOp4YOHapdu3Zp3LhxkqTq6mpVVlb6xtfX12vBggU6deqUoqKilJKSojfffFMTJ04MeKYsOgMAhFTnLjqbJoejVxD7qZfT+ZsOnWuocA0bAAADUBIHABisUYE+D7zl9mYgYAMADBaa27pMQEkcAAADkGEDAAzWoOCWKJvzPmwCNgDAYOETsCmJAwBgADJsAIDBwifDJmADAAwWPgGbkjgAAAYgwwYAGKxRwWXY5tyHTcAGABgs2JK2OSVxAjYAwGDhE7ADuobtdrt1++23Kzo6WjExMZoyZYpOnDjhN8ayLC1dulTx8fGKiorSmDFjdOzYsZBOGgCAcBNQwC4uLtbs2bP13nvvqaioSA0NDcrMzNS5c+d8Y1auXKnVq1dr7dq1OnTokFwul8aNG6czZ86EfPIAgHDXEIJmhqDeh/1///d/iomJUXFxsUaPHi3LshQfH6/s7GwtXLhQklRXV6fY2Fg988wzevTRR6+6T96HDQBm69z3YafJ4YgIYj+NcjpLv/vvw/Z4PJKkPn36SJIqKipUU1OjzMxM3xi73a677rpLBw4caHUfdXV18nq9fg0AAPhrd8C2LEs5OTkaOXKkUlNTJUk1NTWSpNjYWL+xsbGxvu++ze12y+l0+lpCQkJ7pwQACDvhUxJvd8CeM2eOjhw5otdee63FdzabfzHbsqwWfc1yc3Pl8Xh8raqqqr1TAgCEnfAJ2O26reuxxx7Tjh07tG/fPvXr18/X73K5JF3KtOPi4nz9tbW1LbLuZna7XXa7vT3TAAAgbASUYVuWpTlz5mj79u165513lJiY6Pd9YmKiXC6XioqKfH319fUqLi7WiBEjQjNjAAB8yLBbNXv2bG3ZskW/+93vFB0d7bsu7XQ6FRUVJZvNpuzsbOXl5Sk5OVnJycnKy8tT79699eCDD3bIDwAAhLMGXVqX3l7f0UeTFhQUSJLGjBnj179x40b95Cc/kSQ98cQTunDhgmbNmqUvv/xSw4cP1+7duxUdHR2SCQMAEI6Cug+7I3AfNgCYrXPvwx4oh6P9dyh7vU1yOj834j5sniUOADBYo4IriTeFaiIdjoANADBYg4J7Bpg5ATuoJ50BAIDOQYYNADBY+GTYBGwAgMHCJ2BTEgcAwABk2AAAgzUquCy5W93ZfEUEbACAwRoU3FM7zAnYlMQBADAAGTYAwGDhk2ETsAEABgufgE1JHACANiooKNDQoUPlcDjkcDiUkZGht95664rbFBcXKy0tTZGRkUpKStK6devadWwCNgDAXFaTZDUG0QJbYd6vXz+tWLFCJSUlKikp0d/8zd/o7rvv1rFjx1odX1FRoYkTJ2rUqFEqKyvT4sWLNXfuXBUWFgb8U3lbFwAgpDr1bV1/koI5hNcrOfsEN9c+ffro2Wef1U9/+tMW3y1cuFA7duzQ8ePHfX1ZWVn64IMPdPDgwYCOQ4YNADBXYwiaLv0B8M1WV1d39UM3Nmrr1q06d+6cMjIyWh1z8OBBZWZm+vWNHz9eJSUlunjxYkA/lYANAAh7CQkJcjqdvuZ2uy879ujRo7ruuutkt9uVlZWl119/XTfddFOrY2tqahQbG+vXFxsbq4aGBp0+fTqgObJKHABgrm9kye3eXlJVVZVfSdxut192k0GDBqm8vFxfffWVCgsLNXPmTBUXF182aNts/hd4m69Ef7v/agjYAABzNSm4J5P+edvmVd9t0atXL/3whz+UJKWnp+vQoUP61a9+pRdffLHFWJfLpZqaGr++2tpa9ejRQ3379g1oqpTEAQAIgmVZl73mnZGRoaKiIr++3bt3Kz09XT179gzoOARsAIC5QrTorK0WL16sd999V59//rmOHj2qJUuWaO/evZo2bZokKTc3VzNmzPCNz8rK0smTJ5WTk6Pjx4/r5Zdf1oYNG7RgwYKAfyolcQCAuUJUEm+rP/7xj5o+fbqqq6vldDo1dOhQ7dq1S+PGjZMkVVdXq7Ky0jc+MTFRO3fu1Pz58/X8888rPj5ea9as0dSpUwOeKvdhAwBCqlPvwz4ZgvuwB3TsXEOFDBsAYK4mBbdKPJjsvJMRsAEA5grRbV0mYNEZAAAGIMMGAJirkxeddSUCNgDAXGFUEidgAwDMFUYBm2vYAAAYgAwbAGAurmEDAGAASuIAAKA7IcMGAJjLUnBl7W71cO4rI2ADAMxFSRwAAHQnZNgAAHOFUYZNwAYAmCuMbuuiJA4AgAHIsAEA5qIkDgCAAQjYAAAYgGvYAACgOyHDBgCYq0nBlbUNyrAJ2AAAc1ESBwAA3QkZNgDAXKwSBwDAAGEUsCmJAwBgADJsAIC5wmjRGQEbAGAuSuIAAKA7IcMGAJgrjDJsAjYAwFyWgrsObYVqIh2PgA0AMFcYZdhcwwYAwABk2AAAc4XRbV1k2AAAczWGoAXA7Xbr9ttvV3R0tGJiYjRlyhSdOHHiitvs3btXNputRfvoo48COjYBGwCANiouLtbs2bP13nvvqaioSA0NDcrMzNS5c+euuu2JEydUXV3ta8nJyQEdm5I4AMBcnbzobNeuXX6fN27cqJiYGJWWlmr06NFX3DYmJkbXX399gBP8CzJsAIC5mkLQJHm9Xr9WV1fXpsN7PB5JUp8+fa46dtiwYYqLi9PYsWO1Z8+eNv/EZgRsAEDYS0hIkNPp9DW3233VbSzLUk5OjkaOHKnU1NTLjouLi9P69etVWFio7du3a9CgQRo7dqz27dsX0BwpiQMAzBWiknhVVZUcDoev2263X3XTOXPm6MiRI9q/f/8Vxw0aNEiDBg3yfc7IyFBVVZVWrVp11TL6N5FhAwDM1aTgVoj/uSTucDj82tUC9mOPPaYdO3Zoz5496tevX8DTvuOOO/TJJ58EtE1QAdvtdstmsyk7O9vXZ1mWli5dqvj4eEVFRWnMmDE6duxYMIcBAKB1IbqG3VaWZWnOnDnavn273nnnHSUmJrZr2mVlZYqLiwtom3aXxA8dOqT169dr6NChfv0rV67U6tWr9corr+jGG2/U008/rXHjxunEiROKjo5u7+EAAOhys2fP1pYtW/S73/1O0dHRqqmpkSQ5nU5FRUVJknJzc3Xq1Cm9+uqrkqT8/HwNHDhQKSkpqq+v1+bNm1VYWKjCwsKAjt2uDPvs2bOaNm2aXnrpJX3ve9/z9VuWpfz8fC1ZskT33HOPUlNTtWnTJp0/f15btmxpz6EAALi8Tn5wSkFBgTwej8aMGaO4uDhf27Ztm29MdXW1KisrfZ/r6+u1YMECDR06VKNGjdL+/fv15ptv6p577gno2DbLsgJ+V8nMmTPVp08f/eu//qvGjBmjW2+9Vfn5+frss8/0V3/1Vzp8+LCGDRvmG3/33Xfr+uuv16ZNm1rsq66uzm/5vNfrVUJCgqIk2QKdGACgy1mSLujSLU/fXMgVSl6vV06nU541kiMqiP1ckJxzO3auoRJwSXzr1q06fPiwDh061OK75tJAbGysX39sbKxOnjzZ6v7cbrd+8YtfBDoNAADCSkAl8aqqKs2bN0+bN29WZGTkZcfZbP65sWVZLfqa5ebmyuPx+FpVVVUgUwIAhLNOLol3pYAy7NLSUtXW1iotLc3X19jYqH379mnt2rW+B6DX1NT4rX6rra1tkXU3s9vtbbrfDQCAFngfduvGjh2ro0ePqry83NfS09M1bdo0lZeXKykpSS6XS0VFRb5t6uvrVVxcrBEjRoR88gAAhIuAMuzo6OgWj1+79tpr1bdvX19/dna28vLylJycrOTkZOXl5al379568MEHQzdrAACksHofdsgfTfrEE0/owoULmjVrlr788ksNHz5cu3fv5h5sAEDoNT/pLJjtDdGu27o6UvNSfW7rAgAzdeptXSskx+XXQF99P19LzkXf0du6AADoNiiJAwBggDBaJU7ABgCYK4wCNq/XBADAAGTYAABzcQ0bAAADUBIHAADdCRk2AMBcYZRhE7ABAOayFNx16G716LAroyQOAIAByLABAOaiJA4AgAHC6LYuSuIAABiADBsAYC5K4gAAGICADQCAAbiGDQAAuhMybACAuSiJAwBggCYFF3QpiQMAgFAiwwYAmCuMFp0RsAEA5gqja9iUxAEAMAAZNgDAXJTEAXSkc5b/S3ivtdm6aCaA4SiJAwCA7oQMGwBgLjJsAAAM0BSCFgC3263bb79d0dHRiomJ0ZQpU3TixImrbldcXKy0tDRFRkYqKSlJ69atC+zAImADneKcZfk1Hbf5Nb/vALRd85PO2tsCDNjFxcWaPXu23nvvPRUVFamhoUGZmZk6d+7cZbepqKjQxIkTNWrUKJWVlWnx4sWaO3euCgsLAzo2JXEAANpo165dfp83btyomJgYlZaWavTo0a1us27dOvXv31/5+fmSpCFDhqikpESrVq3S1KlT23xsMmwAgLmCya6/cf3b6/X6tbq6ujYd3uPxSJL69Olz2TEHDx5UZmamX9/48eNVUlKiixcvtu13ioANADBZiK5hJyQkyOl0+prb7b7qoS3LUk5OjkaOHKnU1NTLjqupqVFsbKxfX2xsrBoaGnT69Ok2/1RK4kAn+PZ91tyHDXQvVVVVcjgcvs92u/2q28yZM0dHjhzR/v37rzrW9q3/j1t//m/At/uvhIANADBXo4KrFf+5JO5wOPwC9tU89thj2rFjh/bt26d+/fpdcazL5VJNTY1fX21trXr06KG+ffu2+ZgEbACAuTr50aSWZemxxx7T66+/rr179yoxMfGq22RkZOiNN97w69u9e7fS09PVs2fPNh+ba9hAF7jWZvNrAMwwe/Zsbd68WVu2bFF0dLRqampUU1OjCxcu+Mbk5uZqxowZvs9ZWVk6efKkcnJydPz4cb388svasGGDFixYENCxCdgAAHOFaJV4WxUUFMjj8WjMmDGKi4vztW3btvnGVFdXq7Ky0vc5MTFRO3fu1N69e3XrrbfqX/7lX7RmzZqAbumSJJtlda8nNXi9XjmdTkVJIu8AAPNYki7o0i1PgVwXDkRzrPD8WHK0varccj8XJed/duxcQ4UMGwAAA7DoDABgLkvBLTrrVjXmKyNgAwDM1ajgrp8a9LYuAjYAwFxhFLC5hg0AgAHIsAEA5urkB6d0JQI2AMBclMQBAEB3QoYNADAXJXEAAAxASRwAAHQnZNgAAHM1KbgsmZI4AACdoEnBlcQNCtiUxAEAMAAZNgDAXMEuGjNo0RkBGwBgLgI2AAAG4Bo2AADoTsiwAQDmoiQOAIABKIkDAIDuhAwbAGCuYDNkgzJsAjYAwFyNkqwgtjcoYFMSBwDAAGTYAABzhVFJPOAM+9SpU3rooYfUt29f9e7dW7feeqtKS0t931uWpaVLlyo+Pl5RUVEaM2aMjh07FtJJAwAg6VJJPNhmiIAC9pdffqk777xTPXv21FtvvaUPP/xQv/zlL3X99df7xqxcuVKrV6/W2rVrdejQIblcLo0bN05nzpwJ9dwBAAgbAZXEn3nmGSUkJGjjxo2+voEDB/r+t2VZys/P15IlS3TPPfdIkjZt2qTY2Fht2bJFjz76aGhmDQCAxKKzy9mxY4fS09N17733KiYmRsOGDdNLL73k+76iokI1NTXKzMz09dntdt111106cOBAq/usq6uT1+v1awAAtElTCJohAgrYn332mQoKCpScnKz/+I//UFZWlubOnatXX31VklRTUyNJio2N9dsuNjbW9923ud1uOZ1OX0tISGjP7wAAhKMmBXf9+rsasJuamnTbbbcpLy9Pw4YN06OPPqqf/exnKigo8Btns/k/J86yrBZ9zXJzc+XxeHytqqoqwJ8AAMB3X0ABOy4uTjfddJNf35AhQ1RZWSlJcrlcktQim66trW2RdTez2+1yOBx+DQCANqEk3ro777xTJ06c8Ov7+OOPNWDAAElSYmKiXC6XioqKfN/X19eruLhYI0aMCMF0AQD4hjC6rSugVeLz58/XiBEjlJeXp3/4h3/Q+++/r/Xr12v9+vWSLpXCs7OzlZeXp+TkZCUnJysvL0+9e/fWgw8+2CE/AACAcBBQwL799tv1+uuvKzc3V8uWLVNiYqLy8/M1bdo035gnnnhCFy5c0KxZs/Tll19q+PDh2r17t6Kjo0M+eQBAmGtUcK/XDOaWsE5msyyrW03X6/XK6XQqSsH9MwAAdA1L0gVJHo+nw9YlNccKj11yBBEsvJbkrAtsrvv27dOzzz6r0tJSVVdX6/XXX9eUKVMuO37v3r360Y9+1KL/+PHjGjx4cJvnyrPEAQAIwLlz53TLLbfo4Ycf1tSpU9u83YkTJ/z+KPj+978f0HEJ2AAAc3VBSXzChAmaMGFCwNvFxMT4Pco7ULxeEwBgrhCtEv/2Ezfr6upCPtVhw4YpLi5OY8eO1Z49ewLenoANAAh7CQkJfk/ddLvdIdt3XFyc1q9fr8LCQm3fvl2DBg3S2LFjtW/fvoD2Q0kcAGAuSyFZ6V1VVeV3fdlutwe/0z8bNGiQBg0a5PuckZGhqqoqrVq1SqNHj27zfsiwAQDGCtVzU779xM1QBuzW3HHHHfrkk08C2oYMGwBgrGAfVtZVDzorKytTXFxcQNsQsAEACMDZs2f16aef+j5XVFSovLxcffr0Uf/+/ZWbm6tTp0753mSZn5+vgQMHKiUlRfX19dq8ebMKCwtVWFgY0HEJ2AAAYwX7/o72bFtSUuL3IJScnBxJ0syZM/XKK6+ourra91Is6dI7NRYsWKBTp04pKipKKSkpevPNNzVx4sSAjsuTzgAAIdWZTzr7g6RgjuCVFK+OnWuosOgMAAADUBIHABirK0riXYWADQAwlqmrxNuDkjgAAAYgwwYAGKtJwWXJlMQBAOgE4XQNm5I4AAAGIMMGABgrnBadEbABAMYiYAMAYACuYQMAgG6FDBsAYCxK4gAAGICSOAAA6FbIsAEAxuJJZwAAGCCcrmFTEgcAwABk2AAAY4XTojMCNgDAWJTEAQBAt0KGDQAwVjhl2ARsAICxuIYNAIABwinD5ho2AAAGIMMGABjLUnBlbStUE+kEBGwAgLEoiQMAgG6FDBsAYKxwyrAJ2AAAY4XTbV2UxAEAMAAZNgDAWJTEAQAwQDgFbEriAAAYgAwbAGAsFp0BAGCAJv2lLN6e1p6AvW/fPk2ePFnx8fGy2Wz67W9/e9VtiouLlZaWpsjISCUlJWndunUBH5eADQAwVlMIWqDOnTunW265RWvXrm3T+IqKCk2cOFGjRo1SWVmZFi9erLlz56qwsDCg41ISBwAgABMmTNCECRPaPH7dunXq37+/8vPzJUlDhgxRSUmJVq1apalTp7Z5P2TYAABjBVMO/+YKc6/X69fq6upCNseDBw8qMzPTr2/8+PEqKSnRxYsX27wfAjYAwFihCtgJCQlyOp2+5na7QzbHmpoaxcbG+vXFxsaqoaFBp0+fbvN+KIkDAMJeVVWVHA6H77Pdbg/p/m02m99ny7Ja7b8SAjYAwFihuq3L4XD4BexQcrlcqqmp8eurra1Vjx491Ldv3zbvh4ANADCWCU86y8jI0BtvvOHXt3v3bqWnp6tnz55t3g/XsAEACMDZs2dVXl6u8vJySZdu2yovL1dlZaUkKTc3VzNmzPCNz8rK0smTJ5WTk6Pjx4/r5Zdf1oYNG7RgwYKAjkuGDQAwVldk2CUlJfrRj37k+5yTkyNJmjlzpl555RVVV1f7grckJSYmaufOnZo/f76ef/55xcfHa82aNQHd0iVJNqv5ync34fV65XQ6FSWp7ZfiAQDdhSXpgiSPx9Nh14WbY8V6SVFB7OeCpP+njp1rqFASBwDAAJTEAQDGMmHRWagQsAEAxgqnt3URsAEAxgqnDJtr2AAAGIAMGwBgrHDKsAnYAABjhdM1bEriAAAYIKCA3dDQoCeffFKJiYmKiopSUlKSli1bpqamv/yNYlmWli5dqvj4eEVFRWnMmDE6duxYyCcOAECoXq9pgoAC9jPPPKN169Zp7dq1On78uFauXKlnn31Wzz33nG/MypUrtXr1aq1du1aHDh2Sy+XSuHHjdObMmZBPHgAQ3poUXLD+zpbEDx48qLvvvluTJk3SwIED9fd///fKzMxUSUmJpEvZdX5+vpYsWaJ77rlHqamp2rRpk86fP68tW7Z0yA8AACAcBBSwR44cqbffflsff/yxJOmDDz7Q/v37NXHiREmX3lhSU1OjzMxM3zZ2u1133XWXDhw40Oo+6+rq5PV6/RoAAG3RFIJmioBWiS9cuFAej0eDBw9WRESEGhsbtXz5cj3wwAOS5HtBd2xsrN92sbGxOnnyZKv7dLvd+sUvftGeuQMAwlw43dYVUIa9bds2bd68WVu2bNHhw4e1adMmrVq1Sps2bfIbZ7P5v2fLsqwWfc1yc3Pl8Xh8raqqKsCfAADAd19AGfbjjz+uRYsW6f7775ck3XzzzTp58qTcbrdmzpwpl8sl6VKmHRcX59uutra2RdbdzG63y263t3f+AIAwxn3Yl3H+/Hldc43/JhEREb7buhITE+VyuVRUVOT7vr6+XsXFxRoxYkQIpgsAwF+E021dAWXYkydP1vLly9W/f3+lpKSorKxMq1ev1iOPPCLpUik8OztbeXl5Sk5OVnJysvLy8tS7d289+OCDHfIDAADhK5yuYQcUsJ977jn98z//s2bNmqXa2lrFx8fr0Ucf1c9//nPfmCeeeEIXLlzQrFmz9OWXX2r48OHavXu3oqOjQz55AADChc2yLKurJ/FNXq9XTqdTUZJaX6YGAOjOLEkXJHk8Hjkcjg45RnOsyJUUGcR+vpbkVsfONVR4+QcAwFjNTzoLZntT8PIPAAAMQIYNADAWi84AADAA92EDAIBuhQwbAGAsSuIAABiAkjgAAOhWyLABAMaiJA4AgAEI2AAAGMBScNehu9Wzua+Ca9gAABiADBsAYCxK4gAAGCCcAjYlcQAADECGDQAwFg9OAQDAAI0haO3xwgsvKDExUZGRkUpLS9O777572bF79+6VzWZr0T766KOAjknABgAgANu2bVN2draWLFmisrIyjRo1ShMmTFBlZeUVtztx4oSqq6t9LTk5OaDjErABAMZqCkEL1OrVq/XTn/5U//iP/6ghQ4YoPz9fCQkJKigouOJ2MTExcrlcvhYRERHQcQnYAABjhaok7vV6/VpdXV2rx6uvr1dpaakyMzP9+jMzM3XgwIErznXYsGGKi4vT2LFjtWfPnoB/KwEbABD2EhIS5HQ6fc3tdrc67vTp02psbFRsbKxff2xsrGpqalrdJi4uTuvXr1dhYaG2b9+uQYMGaezYsdq3b19Ac2SVOADAWE0K7l7q5pJ4VVWVHA6Hr99ut19xO5vN5vfZsqwWfc0GDRqkQYMG+T5nZGSoqqpKq1at0ujRo9s8VzJsAICxQnUN2+Fw+LXLBewbbrhBERERLbLp2traFln3ldxxxx365JNP2jxeImADAAzW2bd19erVS2lpaSoqKvLrLyoq0ogRI9q8n7KyMsXFxQV0bEriAAAEICcnR9OnT1d6eroyMjK0fv16VVZWKisrS5KUm5urU6dO6dVXX5Uk5efna+DAgUpJSVF9fb02b96swsJCFRYWBnRcAjYAwFiNCq5U3J7r3/fdd5+++OILLVu2TNXV1UpNTdXOnTs1YMAASVJ1dbXfPdn19fVasGCBTp06paioKKWkpOjNN9/UxIkTAzquzbKsbvU6UK/XK6fTqShJrV++BwB0Z5akC5I8Ho/fQq5Qao4VEyX1DGI/FyXtVMfONVS4hg0AgAEoiQMAjNUVJfGuQsAGABiLt3UBAIBuhQwbAGCsUD3pzAQEbACAsRoV3B1FJl3DpiQOAIAByLABAMYKp0VnBGwAgLHCqSROwAYAGCucAjbXsAEAMAAZNgDAWFzDBgDAAJTEAQBAt0KGDQAwlqXgytrd6v3SV0HABgAYK9iSNiVxAAAQUmTYAABjhVOGTcAGABirScGtEjfpti5K4gAAGIAMGwBgLEriAAAYgIANAIABuIYNAAC6FTJsAICxgs2QTcqwCdgAAGOFU8CmJA4AgAHIsAEAxmpUcC/wMCnDJmADAIwVTgGbkjgAAAYgwwYAGCucFp0RsAEAxqIkDgAAuhUybACAsZoUXIYdzLadjQwbAGCsphC09njhhReUmJioyMhIpaWl6d13373i+OLiYqWlpSkyMlJJSUlat25dwMckYAMAjNUYghaobdu2KTs7W0uWLFFZWZlGjRqlCRMmqLKystXxFRUVmjhxokaNGqWysjItXrxYc+fOVWFhYUDHtVmW1a0qAl6vV06nU1EK7g0sAICuYUm6IMnj8cjhcHTIMZpjxXUKLlZYks4qsLkOHz5ct912mwoKCnx9Q4YM0ZQpU+R2u1uMX7hwoXbs2KHjx4/7+rKysvTBBx/o4MGDbZ5rt7uG3fz3Q7f6KwIA0GbN//3ujHywUcEHbOnSHwDfZLfbZbfbW4yvr69XaWmpFi1a5NefmZmpAwcOtHqMgwcPKjMz069v/Pjx2rBhgy5evKiePXu2aa7dLmCfOXNGkvR1F88DABCcM2fOyOl0dsi+e/XqJZfLpZqamqD3dd111ykhIcGv76mnntLSpUtbjD19+rQaGxsVGxvr1x8bG3vZudTU1LQ6vqGhQadPn1ZcXFyb5tntAnZ8fLyqqqpkWZb69++vqqqqDiupfBd4vV4lJCRwnq6C89Q2nKe24TxdmWVZOnPmjOLj4zvsGJGRkaqoqFB9fX3Q+7IsSzabf57eWnb9Td8e39o+rja+tf4r6XYB+5prrlG/fv185QmHw8H/IdqA89Q2nKe24Ty1Defp8joqs/6myMhIRUZGdvhxvumGG25QREREi2y6tra2RRbdrLVKQG1trXr06KG+ffu2+disEgcAoI169eqltLQ0FRUV+fUXFRVpxIgRrW6TkZHRYvzu3buVnp7e5uvXEgEbAICA5OTk6Ne//rVefvllHT9+XPPnz1dlZaWysrIkSbm5uZoxY4ZvfFZWlk6ePKmcnBwdP35cL7/8sjZs2KAFCxYEdNxuVxJvZrfb9dRTT131OkK44zy1DeepbThPbcN5Cm/33XefvvjiCy1btkzV1dVKTU3Vzp07NWDAAElSdXW13z3ZiYmJ2rlzp+bPn6/nn39e8fHxWrNmjaZOnRrQcbvdfdgAAKAlSuIAABiAgA0AgAEI2AAAGICADQCAAbptwA701WXfZW63W7fffruio6MVExOjKVOm6MSJE35jLMvS0qVLFR8fr6ioKI0ZM0bHjh3rohl3D263WzabTdnZ2b4+ztMlp06d0kMPPaS+ffuqd+/euvXWW1VaWur7nvMkNTQ06Mknn1RiYqKioqKUlJSkZcuWqanpLy9k5DyhU1nd0NatW62ePXtaL730kvXhhx9a8+bNs6699lrr5MmTXT21LjF+/Hhr48aN1v/8z/9Y5eXl1qRJk6z+/ftbZ8+e9Y1ZsWKFFR0dbRUWFlpHjx617rvvPisuLs7yer1dOPOu8/7771sDBw60hg4das2bN8/Xz3myrD/96U/WgAEDrJ/85CfWf//3f1sVFRXWf/7nf1qffvqpbwznybKefvppq2/fvta///u/WxUVFda//du/Wdddd52Vn5/vG8N5QmfqlgH7r//6r62srCy/vsGDB1uLFi3qohl1L7W1tZYkq7i42LIsy2pqarJcLpe1YsUK35ivv/7acjqd1rp167pqml3mzJkzVnJyslVUVGTdddddvoDNebpk4cKF1siRIy/7PefpkkmTJlmPPPKIX98999xjPfTQQ5ZlcZ7Q+bpdSbz51WXffhXZlV5dFm48Ho8kqU+fPpIuvRy9pqbG75zZ7XbdddddYXnOZs+erUmTJunHP/6xXz/n6ZIdO3YoPT1d9957r2JiYjRs2DC99NJLvu85T5eMHDlSb7/9tj7++GNJ0gcffKD9+/dr4sSJkjhP6Hzd7kln7Xl1WTixLEs5OTkaOXKkUlNTJcl3Xlo7ZydPnuz0OXalrVu36vDhwzp06FCL7zhPl3z22WcqKChQTk6OFi9erPfff19z586V3W7XjBkzOE9/tnDhQnk8Hg0ePFgRERFqbGzU8uXL9cADD0ji3yd0vm4XsJsF+uqycDFnzhwdOXJE+/fvb/FduJ+zqqoqzZs3T7t3777iG3zC/Tw1NTUpPT1deXl5kqRhw4bp2LFjKigo8Hv+cbifp23btmnz5s3asmWLUlJSVF5eruzsbMXHx2vmzJm+ceF+ntB5ul1JvD2vLgsXjz32mHbs2KE9e/aoX79+vn6XyyVJYX/OSktLVVtbq7S0NPXo0UM9evRQcXGx1qxZox49evjORbifp7i4ON10001+fUOGDPE9+5h/ny55/PHHtWjRIt1///26+eabNX36dM2fP19ut1sS5wmdr9sF7Pa8uuy7zrIszZkzR9u3b9c777yjxMREv+8TExPlcrn8zll9fb2Ki4vD6pyNHTtWR48eVXl5ua+lp6dr2rRpKi8vV1JSEudJ0p133tnitsCPP/7Y9+IC/n265Pz587rmGv//REZERPhu6+I8odN14YK3y2q+rWvDhg3Whx9+aGVnZ1vXXnut9fnnn3f11LrEP/3TP1lOp9Pau3evVV1d7Wvnz5/3jVmxYoXldDqt7du3W0ePHrUeeOABbi+xLL9V4pbFebKsS7e89ejRw1q+fLn1ySefWL/5zW+s3r17W5s3b/aN4TxZ1syZM60f/OAHvtu6tm/fbt1www3WE0884RvDeUJn6pYB27Is6/nnn7cGDBhg9erVy7rtttt8tzCFI0mtto0bN/rGNDU1WU899ZTlcrksu91ujR492jp69GjXTbqb+HbA5jxd8sYbb1ipqamW3W63Bg8ebK1fv97ve86TZXm9XmvevHlW//79rcjISCspKclasmSJVVdX5xvDeUJn4vWaAAAYoNtdwwYAAC0RsAEAMAABGwAAAxCwAQAwAAEbAAADELABADAAARsAAAMQsAEAMAABGwAAAxCwAQAwAAEbAAADELABADDA/wfdRY+WaGDTKQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.colors import LogNorm\n",
    "# Plot a slice at z = nz // 2\n",
    "plt.imshow(final_temperature[:, :, 99], cmap='hot')#, vmin = 0, vmax = 1)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "id": "729e7696-8f02-4686-b4b9-02b0ea81e485",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0007617260670990519"
      ]
     },
     "execution_count": 90,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "final_temperature[nx//2,ny//2,0]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a8d2a50-0260-48c9-97fa-8b0a0c5b0c93",
   "metadata": {},
   "source": [
    "#### Using CUPY 4 seconds to completion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "id": "3e56fa89-e428-4068-aefc-9359e3f2e98f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# %%memit\n",
    "# %%time\n",
    "import cupy as cp\n",
    "\n",
    "# Grid parameters\n",
    "nx, ny, nz = 100, 100, 100  # Number of grid points\n",
    "dx = dy = dz = 1.0          # Grid spacing\n",
    "dt = 1                      # Time step\n",
    "alpha = 0.15                # Thermal diffusivity\n",
    "nt = 1000                    # Number of time steps\n",
    "\n",
    "def main():\n",
    "    # Initialize the temperature field with initial conditions (using CuPy)\n",
    "    u = cp.zeros((nx, ny, nz))\n",
    "    u[nx // 2, ny // 2, 0] = 1.0  # place heat source here \n",
    "\n",
    "    for _ in range(nt):\n",
    "        # Compute second derivative along the z-axis\n",
    "        d2u_dz2 = (cp.roll(u, -1, axis=2) - 2 * u + cp.roll(u, 1, axis=2)) / dz**2\n",
    "        d2u_dz2[:, :, 0] = (u[:, :, 1] - u[:, :, 0]) / dz**2  # Neumann BC at z=0\n",
    "        d2u_dz2[:, :, -1] = (u[:, :, -2] - u[:, :, -1]) / dz**2  # Neumann BC at z=nz-1\n",
    "\n",
    "        # Compute second derivative along the x-axis\n",
    "        d2u_dx2 = (cp.roll(u, -1, axis=0) - 2 * u + cp.roll(u, 1, axis=0)) / dx**2\n",
    "        d2u_dx2[0, :, :] = (u[1, :, :] - u[0, :, :]) / dx**2  # Neumann BC at x=0\n",
    "        d2u_dx2[-1, :, :] = (u[-2, :, :] - u[-1, :, :]) / dx**2  # Neumann BC at x=nx-1\n",
    "\n",
    "        # Compute second derivative along the y-axis\n",
    "        d2u_dy2 = (cp.roll(u, -1, axis=1) - 2 * u + cp.roll(u, 1, axis=1)) / dy**2\n",
    "        d2u_dy2[:, 0, :] = (u[:, 1, :] - u[:, 0, :]) / dy**2  # Neumann BC at y=0\n",
    "        d2u_dy2[:, -1, :] = (u[:, -2, :] - u[:, -1, :]) / dy**2  # Neumann BC at y=ny-1\n",
    "\n",
    "        # Update the temperature field\n",
    "        u += alpha * dt * (d2u_dx2 + d2u_dy2 + d2u_dz2)\n",
    "\n",
    "    # Optionally, return or save the final temperature field\n",
    "    return u\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    final_temperature = main()\n",
    "    # If you need to bring the result back to the CPU:\n",
    "    final_temperature_cpu = cp.asnumpy(final_temperature)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "id": "619bc945-fe75-41c0-b1a1-52e5ab4eb876",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAewAAAGsCAYAAAAfcpQMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABKt0lEQVR4nO3dfZBU1Zk/8G8zMC/gzLCizDARdKxCVDA6AaPgC6RQXDTWWlqu75Jkt0oW3yZTBkX9RVRk1BiKMq4YLCNahGhtgVldzS5sFNDFXRUlsqiocQITwuwsic7wOgMz9/cHS3v7zPRz+ulzu/ue6e+nqqvm9jn3pW93z+l7nvuckwiCIAARERHF2qBCHwARERHZscEmIiLyABtsIiIiD7DBJiIi8gAbbCIiIg+wwSYiIvIAG2wiIiIPsMEmIiLyABtsIiIiD7DBJiIi8gAbbCIiyqv169fj0ksvRV1dHRKJBH7961/ndH/z589HIpFIedTW1uZ0n7nABpuIiPJq7969OP300/HEE0/kbZ/jx4/Hzp07k4/Nmzfnbd9RGVzoAyAiouIyc+ZMzJw5M215d3c37r33Xvzyl7/EV199hQkTJuCRRx7BtGnTst7n4MGDvbyqDuMVNhERxcr3v/99/Md//AdeeOEFfPjhh7jyyivx13/91/jss8+y3uZnn32Guro61NfX4+qrr8YXX3wR4RHnR4LTaxIRUaEkEgm89NJLuOyyywAAv//97zF27Fj88Y9/RF1dXbLeBRdcgG9/+9tYuHCheh+/+c1vsG/fPpx00kn4n//5HyxYsACffPIJtmzZghEjRkT1UnKOV9hERBQb77//PoIgwEknnYSjjjoq+Vi3bh1+//vfAwD+8Ic/9LmJzHzccsstyW3OnDkTV1xxBU477TRccMEFePXVVwEAzz33XEFeY7YYwyYiotjo7e1FSUkJNm7ciJKSkpSyo446CgDwjW98Ax9//LG4nb/6q79KWzZs2DCcdtppTl3shcAGm4iIYqOhoQE9PT1ob2/Heeed12+dIUOG4OSTT856H11dXfj444/Tbj+u2GATEVFe7dmzB59//nlyuaWlBZs2bcLRRx+Nk046Cddddx1uvPFG/PSnP0VDQwN27dqF119/Haeddhouvvhi9f7uuOMOXHrppRgzZgza29uxYMECdHZ2YtasWVG+rJzjTWdERJRXa9euxXe+850+z8+aNQvLli3DwYMHsWDBAjz//PPYsWMHRowYgcmTJ+P+++/Haaedpt7f1VdfjfXr12PXrl049thjcfbZZ+PBBx/EqaeeGsXLyRs22ERERB7gXeJEREQeyFmD/eSTT6K+vh7l5eWYOHEi3nzzzVztioiIaMDLyU1nL774IhobG/Hkk0/inHPOwc9//nPMnDkTH330EcaMGSOu29vbiz/96U+orKxEIpHIxeEREVEOBUGA3bt3o66uDoMG5a4j98CBA+ju7nbeTmlpKcrLyyM4otzKSQz7rLPOwre+9S0sWbIk+dwpp5yCyy67DM3NzeK6f/zjHzF69OioD4mIiPKstbUVxx13XE62feDAAdTX16Otrc15W7W1tWhpaYl9ox35FXZ3dzc2btyIu+66K+X5GTNmYMOGDX3qd3V1oaurK7l85PdDOQBeXwMl9ipZcfnNm+0xueyzEMerla8bQnrztJ+eLNdzOT6XdbM9Xtf9SlyOyWcBgAMAKisrc7aP7u5utLW1obW1BVVVVVlvp7OzE6NH16O7u7v4Guxdu3ahp6cHNTU1Kc/X1NT0+0uoubkZ999/f5/nE2CDDeTuHLhsN9t1C9Vg56shHWh3cBYifcRln4X4TBdqu77IR1izqqrKqcH2Sc7+x5hvVBAE/b558+bNQ0dHR/LR2tqaq0PyQonxiMog40F65jks1DmNwzEMNLk6j7n6PlPYoQgefoj8CvuYY45BSUlJn6vp9vb2PlfdAFBWVoaysrKoD4OIiIqCa6PrT4Md+Y/y0tJSTJw4EWvWrEl5fs2aNZgyZUrUuyMioqLGK2wnTU1NuOGGGzBp0iRMnjwZS5cuxfbt2zF79uxc7I5yzKUrL183qEXV3eh7t7Lt+DU3V0nnVLqZyjwGl31qbtpyWZfIBzlpsK+66ir8+c9/xgMPPICdO3diwoQJeO2113D88cfnYndERFS0euB2lezPT7vYjSXe2dmJ6upqVKA477AcSGlcLuvyCjs3okpf0vyLK1SqFtO8CicAsB9AR0dHzu7gPtJWdHR8iqqq7NPHOjt3o7r6pJwea1QG+v8nIiKiAYHzYRdYMV9RF+oYo/qVGoc0He0Vm/TacxVrdomrRxXTdjlPUV5tuxwTpVM8d4mzwSYiIo8VT4PNLnEiIiIP8AqbiIg81oPC3Z6YX2ywB5Bsu0viGLPWHlMhXnuu9uOSe6xh7kdzDs24brY52+Z+bfHibGPALrFwl7xyyofiSetilzgREZEHeIVNREQeK56bzthg51nUM3BlK1fd4Lnq9s5ld3pU60Yll8cQ7vxz6abXpIfZ9hPetiYFrFDDmEbVRc6hVKPCBpuIiMgDxdNgM4ZNRETkAV5hExGRx4rnLnE22HkQh4kpNMcQZWy5EHHqOAxbCkT3vvvw7yTbODSgSwmTUsBcUsuiimm7pHwxpp0tdokTERFRjPAKm4iIPFY8V9hssImIyGNssMlBHHKto4xZRxVbjjLeXYhc8ELlaEc5rGncp3e0xaXzlbOd7XmKchhTxrTJxAabiIg8xitsIiIiDxRPWhfvEiciIvIAr7AjEIeYNRBdbrJLrFla13YMvsW/Xfaj4RJDNWUb17XVjSoX2UVUOdvmtuIw7rh5HP5cE+YDu8SJiIg8wAabiIjIA2ywyaIQw43GJVUr225vl2My181XWldUXeS5msrSRuoKjnLKzDh20UY1XKr52lyGQOUwpuSCDTYREXmMV9hEREQeYFoXERERxQivsDNUqCky85GuVKjYsibeHcdUM5tsfw0PMZZdYp3StUOUKWDhbdnqRhV/dYn129bNVfpbpvvUKu6Ydg/cXrE/Z4sNNhEReax4YtjsEiciIvIAr7CJiMhjxXOFzQY7jUINN5qvITbD62rjw/mIU2vj6tnmf7vExgs13WaYJu/XpBme07ZdKearidXmKl5s4zKsadxj2v5EaLPFu8SJiIgoRniFTUREHmOXeFGK43Cj+UjV0h5DtuvaurFd0sfyMayp9vMRVfeV1DVqpoCZnXuaFDFNmpTU1atJmTLLbd3AvnX1Silt2nBGtl3kAz/liw02ERGRB4qnwWYMm4iIyAO8wiYiIo8VzxV2UTfYAylmrV23JM3frutKdW0xbJd0Mc1+NNvNtCyTbWdLE7u0xSfD5eZ2pXi4WVczNKl0DOa6Uny7v/JcyXZYVk0KmDa2nKupOTXHEE9M6yIiIqIYKeorbCIi8t0huPWXskuciIgoD9hgD0hxHG40yqkts81F1saWpXhlVPFuTV3bMeUrhp2v+JIml1oqt8VBpdxqKaatiXcDcm6yFMO2vXfS8UcpvG1bzF16rYUaxjTbY6D8K6oGm4iIBhpeYRMREXmgeO4SH9ANdrF1gUc1M1aUqVpDhDKXrvZSxbouoQOXoUqj+vxJ/05s/2ps3dOZ1rV1iXcLZZqua9/TVjRd+rbu5zh2kZv8aeoGhgHdYBMR0UB3CG4/9dglTkRElAdssImIiDzABtsbUcapw+I43GiuhvY0h6R0Sb+StuUS73aJYUeZPpZpWX/bylZUaVzmtswyTV3NFJnmugeFdW2f6fC6Lilg+RLVMKZmeb6GMbVJ9x0IcrjPYuZ9g01ERMWsB263v/lz6xwbbCIi8ljxpHX5nkVBRERUFLy4ws5VnNqUj1xr7RSZ2eYXFyrnWYph22LlUgzbXNclJq+Jq2c73GsuaYbcdBlCNLxsxp3NuuFzbqsrxWql+LbJ9t51I3P5GtY022FMzXXjkKMdH4cAJBzXz1xzczNWrVqFTz75BBUVFZgyZQoeeeQRjBs3Tlxv3bp1aGpqwpYtW1BXV4e5c+di9uzZqn3zCpuIiDx2KIJH5tatW4ebb74Z//mf/4k1a9bg0KFDmDFjBvbu3Zt2nZaWFlx88cU477zz8MEHH+Duu+/GbbfdhpUrV6r27cUVNhERURz867/+a8rys88+i5EjR2Ljxo04//zz+13nqaeewpgxY7B48WIAwCmnnIL33nsPjz32GK644oqM980Gm4iIPBZNl3hnZ2fKs2VlZSgrK7Ou3dHRAQA4+uij09Z5++23MWPGjJTnLrroIjzzzDM4ePAghgwxA3n9i22DXQK3tyATvo0Pnqu6mtiybV0p/m2LjYe3pYmj2/YjHX+U59SUj7HENXFos76tbjiebJ5/Mz4s5TzbYrPdGZaZ5bZYeaZlhaLJIwfiPzVn4UTTYI8ePTrl2fvuuw/z588X1wyCAE1NTTj33HMxYcKEtPXa2tpQU1OT8lxNTQ0OHTqEXbt2YdSoURkdaWwbbCIionxpbW1FVVVVcjmTq+tbbrkFH374Id566y1r3UQi9UdFEAT9Pi9hg01ERB7rgdsV9uH+iKqqqpQG2+bWW2/Fyy+/jPXr1+O4444T69bW1qKtrS3lufb2dgwePBgjRozIeJ9F1WBHOdxotvvRDoWZ7dCkmm7hKLvEzWVNV3WpUKY5Rpf0sWLuEpe6mM0y8xz3CGW2VC1pGE2zrrktqW6mZdr95EqupuaMyzCmueM6Frhu/SAIcOutt+Kll17C2rVrUV9fb11n8uTJeOWVV1KeW716NSZNmpRx/BpgWhcREXktv2ldN998M5YvX44VK1agsrISbW1taGtrw/79+5N15s2bhxtvvDG5PHv2bGzbtg1NTU34+OOP8Ytf/ALPPPMM7rjjDtW+VQ12c3MzzjzzTFRWVmLkyJG47LLLsHXr1pQ6QRBg/vz5qKurQ0VFBaZNm4YtW7aoDoqIiCiOlixZgo6ODkybNg2jRo1KPl588cVknZ07d2L79u3J5fr6erz22mtYu3YtzjjjDDz44IN4/PHHVSldgLJL/EjC+JlnnolDhw7hnnvuwYwZM/DRRx9h2LBhAIBHH30UixYtwrJly3DSSSdhwYIFuPDCC7F161ZUVlaqDo6IiEiW/y5xm2XLlvV5burUqXj//fdV+zIlgkz2nsb//u//YuTIkVi3bh3OP/98BEGAuro6NDY24s477wQAdHV1oaamBo888ghuuukm6zY7OztRXV2No+Ce1uXS36+ZMjOXU2RmO9yoJo5rS5lyiTVL+5HqSvvsr1wTV892atFcTsUZppky0xbDlsptaVEH0/xtq2umYpnrSsuadc0yaV3zPNj2I50n6Zxqzr+mrklT1+SS4pZpjDsAsAeHc5Q1N3JpHGkrOjomoqoq+ztGOjt7UF29MafHGhWnGLaZMN7S0oK2traUBPGysjJMnToVGzZs6HcbXV1d6OzsTHkQERFRqqwb7P4Sxo/ctt5fgrh5S/sRzc3NqK6uTj7M5HUiIqL08nvTWSFl3WAfSRj/1a9+1aesvwTxdMnh8+bNQ0dHR/LR2tqa7SEREVHRKZ4GO6s87HQJ47W1tQAOX2mHh1prb2/vc9V9RKbjtWYqV3FrTYQkqji0WW6rq5meUopDa4YQLbfsRxPDLhfKNPnemnU1OecuMewoZTtFJpAaj7XlYWcbL7blYR8wll3u8chUlPnFhcjRBuTYs2ZoUk1d08Ab1tQvqvYtCALccsstWLVqFV5//fU+CeP19fWora3FmjVrks91d3dj3bp1mDJlSjRHTERElMQr7H7dfPPNWLFiBf75n/85mTAOANXV1aioqEAikUBjYyMWLlyIsWPHYuzYsVi4cCGGDh2Ka6+9NicvgIiIitkhHL4vPVtxnBqmf6oGe8mSJQCAadOmpTz/7LPP4nvf+x4AYO7cudi/fz/mzJmDL7/8EmeddRZWr16dsxzsQnWBS93cLnI13KjUHa3pbjbLNV3itrpSl7h5DNK2XLrPNbOAabprNZ9TWzejJjVImlVL081tS80Kd3PbUqakYUBdzrEUsshlN7Zm2y5dyNJr18zsFaV0n2uX5pPSUzXYmaRsJxIJzJ8/3zotGRERkbseuP1E8CcSX1STfxAR0UDDBpuIiMgDh+AWHGWDHamophTTDDfqcgwuqVm5Gm7UZRhQM+kuXK6JS+eqrlnukmqmSZXLV5qXZopMW6pWeF3NMKBmKpY03aYmjcusb5bZ1g2Lw9SDtnh2VMOC2j5buZqKkwrLiwabiIiof7zCJiIi8kDxNNhx6EUiIiIii9heYQ+C+/SaLjHEXE2ZqYlZm+Uuw41GlR9t7seMb5t1pbi0tGw7Bpd1s53G02VqVBcuUzRKudaa3GrzHJqx5SEZlvW3HK5vnjNz2dy2VDfTsijZrtWiuq/Btp/w65VytIHUz5Dt+OIZ4+6B21WyP1njsW2wiYiI7A7B7fLOnwabXeJEREQeGHBX2C6zamnqZjsDl2Z2KLM8yuFGo+p+dqkrdaebdYdZ9iN18UeVAqZN45JCIRLNDFy2ui6pWgfS/A3IXeR7hTJAF/aRlnMVghgIctV17TLTV+4UzxX2gGuwiYiomBRPg80fpURERB7gFTYREfkr6HW7SPbnAtv/BjvKmHWuuhs0MTqXYU2lYUw1U2S6xKWHCnVtcehwTNtWV9pvlEOgStNrJmxzc+ZhbNLgYNoiAPLwo2Zs2VweKpTtM5al+yfMmLb5WdSkdeWKlBRknn7pbbbdU5DpPjMpl+RrKs50+8lrG9gLt5Plz7gp/jfYRERUxHrgdvdbPO6cywhj2ERERB7gFTYREfmriK6wvWuwowoLAnL3QpSxZSkf11x2ya0uTfM3oMt5dokXS7FnM75t1tVs13w94XLN6ynRjOFquzGgADHshBGkLjFi2hXC8jCjTIpp24YbDZ9TM2Zt1jXj35rvYabrRSnKEGf47bIdf1RTcZrM/Wqm4oylIophs0uciIjIA95dYRMRESWxS5yIiMgDRdQlHtsGuwT5n15TM/ZztlNmRpmHrYl3m3XDy5qYtVkuxaHNck0MW7Nd85jMugmXAH5UMewoBxMPl0uJ1v0th4LRJUZgepixPLTr67/NuLN5msJxa1setub7kauYnXmKbcuZbsv21kk526Yo79cJk/LKNVNxUv7FtsEmIiKy6kXu7tCLGTbYRETkL8aw/ZWrNBBpu5opATVTMgLZT5lp6+ktE8psy5pu7sos62q7xEs045pK+WNF1CVuG280ESofZvRrl3elLodPjRSqAXRd4i40w41qSNOd2ur2CmUm6fijHLZUcy6k98ejds9bA67BJiKiIsKbzoiIiDzALnEiIiIPsMH2h0vMOsr0K6lM2q4tLi1NmakZRVMT73ZJv6o0lsN1qxR1zbIK86SaFTSBdU28W5qH1PZmZhuclQKfQGpc2qxrxqWlGLaZb2Wei3C58dpLjHWrdn/99xDjeLX3beSD7RT3pPnbtmyWmR8JiXk7gnneeoQyjWyHLbXRpKxRdrxvsImIqIgxhk1EROSBIuoS5+QfREREHvDuCtsl3hVlrEwaMlRatuWkSsOP2vJbpRi2lG7skoctxaGB1Li1WVdaLjWHE5Vi1ma5SxBeevGamwgA+UMisQVND6b5u79lac5M8zyZ449KyfrCa6/YnVpUYuRs5ypmLcWhNantZrnt7ZCGJpWWpWFLM9lWpmUm6SpNGrZUu5+8CeB2YEFUB5J73jXYRERESewSJyIiojiJ7RX2IGT3a0KzjmZmrGzTuLTblTKDzLouM3BJXeK2YUGHCmVS6pbZqz3cWC6R8ro0y5pcM02XuK1b2PYGSXWlLj0z3ydcV9MFbi7busQ1rz28bHzIS40u8uHmfrKkGSLUpUvc1jUdfgvM0yIdo60L3NY9LdWVSHVtw5bGsou8iK6wY9tgExERWRVRWhe7xImIiDzAK2wiIvIXu8T9IXUR5CqNy1zWDGOqSQEDUsOguRqa1MygssW0hwllUmh5uFFWYsZQwxXMDVUby5pgubkfKd4txbDNE6XJszNp5jiUUrfM+LaRQqWKYRuxZjGGLb12y4faLNbEtDWxZWkWUilTDpC/S9Lwo+YxSKPX2mLWUozbth9pu6bw/ydt73AsumjZYBMREXmAMWwiIiKKE15hExGRv3rh1q3t0RV2UTXYmmn+bF0PUn60ZruauLQtZJrtcKNmmRmqlfKwNeFiMWZtVjZj1mZdaV1bzrY0tqo016gtF9l8g6SbF6SEVluScDgga8vDNuPD4XLbhyQc09bMA6sUPjWVxvFqprI0l8OnyXxpUi61uS1brDn8ttvqSlNkavKypak3TbZ7eVzysmOBXeJEREQUJ0V1hU1ERAMM7xKPL1uXQK5mArKlX2W6rjatK/x6XdK6pN5cWxqXVG5L60qZdcvWVR3uBh9ulNmWw+tKKV/msq1PXxNnkGbr0vRl2fpVw/23mqFIgdQuctvrkXIKo/yihV5fqfFaK400NWlUVnM5/Nm0dYGbLz3cna753pnpY5p0UE2aV75m8rKlj8WirSuiBptd4kRERB7w7gqbiIgoqYhuOmODTURE/iqiLvHYNtjZTq9pbiPb+prpNW11pRQwW5qXFEaU0rw0KWCamDUgT69ZYb4gKV4spW4NN8qOFuqay7Z5PMPlUh4aIJ8oTXBTE/O1jZt5UCjbayybMexwqpYmJm+LYbt8WYVcrQrj9R0MXQ2Z8WLNaTJfqllemmGZWW6WSRl5tiFPzVNsHkeYlOaleWu0bVe6j3VCuR3KTGwbbCIiIiteYRMREXkggFscOojqQHKPDTYREfmLV9j+yDYf2raubbvZTplpxrik0Sxt+5FCpi5Dk9qGKpXysFXDgJr50sPT/J3JsiaHW4qrSwF7W7DffEMS4XfP/KqZy4fS/A0gMC4fpDxsMyZvDk0aPmZzOk1puFGXPHJbuRRsNpaHdaSvap6K8EvVxrClGUylcL5ZZq4rDWNqi1mHy23Dskp1NWx52ZRf3jfYRERUxIoorYsDpxARkb96IngorV+/Hpdeeinq6uqQSCTw61//Wqy/du1aJBKJPo9PPvlEtV9eYRMRESns3bsXp59+Or7//e/jiiuuyHi9rVu3oqrq61jgscceq9qvdw22JrYc5bY1U2ZKKaq2PGxzXc1wzlJdKQ/bLDPDoFLotkQKcAO6aS8102ua5SMUdcP7sY07Hj4Zg80zc5SxLN0NoPmqGTHshBGdLQ0tl+5JLRtqBK2lGLbtQyINfC2xBVjN4Kw0XaixXBJ66eXGOOPmR08K9ZvvlFkujWlg7DalrhR3Btz+F/QIZdLY4i7DvptvXSxj2gW46WzmzJmYOXOmer2RI0di+PDh+h3+H3aJExGRv3ojeADo7OxMeXR1mT/N3DU0NGDUqFGYPn063njjDfX6bLCJiKjojR49GtXV1clHc3NzZNseNWoUli5dipUrV2LVqlUYN24cpk+fjvXr16u2412XuIZm6kqTyy8ZqWvLdkxSZk2+pte0Lad0DEtpUOaybS5OTVf1cGHZNoxppVBWanaAHhP62+wCty2H312HtK4+fXZ70vwNYLCxXLUrdbk81P9sG+tW6ks1D0kz76WUj6WYLnRoV9oiAHLqoqaLXPO9sw01HH7p2v9P4XLzFErd67lMMw7vp2Dd4xF1ibe2tqbEl8vKzHhf9saNG4dx48YllydPnozW1lY89thjOP/88zPeDq+wiYjIX71wu0P8/35pVFVVpTyibLD7c/bZZ+Ozzz5TrePUYDc3NyORSKCxsTH5XBAEmD9/Purq6lBRUYFp06Zhy5YtLrshIiLqX0Qx7Hz74IMPMGrUKNU6WXeJv/vuu1i6dCm++c1vpjz/6KOPYtGiRVi2bBlOOukkLFiwABdeeCG2bt2KykqzX5OIiMgve/bsweeff55cbmlpwaZNm3D00UdjzJgxmDdvHnbs2IHnn38eALB48WKccMIJGD9+PLq7u7F8+XKsXLkSK1euVO03qwZ7z549uO666/D0009jwYIFyeeDIMDixYtxzz334PLLLwcAPPfcc6ipqcGKFStw0003ZbyPEsSrv94lLi3VtaVySDM0amLYUkzbNsKmWZ4oEwrNOLU0NKk0taUmZg3IQ5NKKWKDzWD5McKy+eJsMeyI0rr6RFyFGHaf+TWNN6g0FNM+ujO1TPrC2VK1pBi2OT6n+XIqhTLzRt1QmlrC6LE007w092lIaV2a7DcpZg1kP6SxuS0p5cukqeuybsH+Xxcgreu9997Dd77zneRyU1MTAGDWrFlYtmwZdu7cie3btyfLu7u7cccdd2DHjh2oqKjA+PHj8eqrr+Liiy9W7TerBvvmm2/GJZdcggsuuCClwW5paUFbWxtmzJiRfK6srAxTp07Fhg0b+m2wu7q6Um6f7+zs7FOHiIioXwUYmnTatGkIgvTTfC1btixlee7cuZg7d65+RwZ1g/3CCy/g/fffx7vvvtunrK2tDQBQU1OT8nxNTQ22bdvW7/aam5tx//33aw+DiIioqKh6MVpbW3H77bdj+fLlKC83O5W+lkgkUpaDIOjz3BHz5s1DR0dH8tHa2qo5JCIiKmYFGEu8UFRX2Bs3bkR7ezsmTpyYfK6npwfr16/HE088ga1btwI4fKUdvvutvb29z1X3EWVlZU63z0cZN8l2CD/NeraYtbRsG5pUGk7RJWe7z0+z8BNmHNp8K4em+RvQTcWpWbbFvweHE7VrjUIzhi0FxzUxbJexPTUx7K+EYzCWB7elFg3/S+pyuKvQjENLudWKXOo+y+ZnxAzJD01fZsaww6/cHKHVFpfWDAks3WeiGdLYJP2vkKbTBFLfOts0nhIv2rIimg9b1d5Nnz4dmzdvxqZNm5KPSZMm4brrrsOmTZtw4oknora2FmvWrEmu093djXXr1mHKlCmRHzwREVGxUF1hV1ZWYsKECSnPDRs2DCNGjEg+39jYiIULF2Ls2LEYO3YsFi5ciKFDh+Laa6+N7qiJiIiAopoPO/KhSefOnYv9+/djzpw5+PLLL3HWWWdh9erVecvB1nT3ZDs7F6Ab1lSqq0nzsqV1Sd3nUhe5rQu8RJMTJqV1SWWA3K1t604PZ2cNN8r6pG6Fu8HNLnAzdDM8zd9A33wxzWxdmqFJzT7k8H47FMdgMdjY7/BQxoZmeFHzvbF1ie9L87dt2XhpJUbd8oNpq6q+H5rvnTZtU6qrLc90PZceYJcUsZw5MtKZy/qecG6w165dm7KcSCQwf/58zJ8/33XTRERE9H8G9OQfREQ0wLFLnIiIyANFdJe4dw22Jo6jiTtrty2tF+XQpIPS/N1fXWmaP3NZSlsxl1XjmkrLmrqaFDCzXJwi01w2Y9ZSXTNmHYehSW0nVbNfYz+loeXKbrGqOEVmrj4jlsD0kINpi/rErKXviy2tS/MdzdXQpNIQqB61R9kpogY7TsN1ExERURreXWETERElMYZNRETkgSLqEmeDnYZ2CNGotqvJw5bqugxNaq6bs3k8zeWhQplUt8+yFIcG5OFGzbojhLq2GHZ4nFaXPGxzjknppNq+0uFtK8YMHfqn1CJpCFHbe+UQl9Z89jT3aUjfgVx9R7XxyPC6Zsw6V9Nt2upKw6VS9NhgExGRv3iFTURE5IEAbnHo9NNax44XDXaubmXPtlvbtp2ous/NdV1mAtKkl/Q53y5TfUljoJoze2U75CkADA73u5pd02bl4Wn+Bvqmbmnqavp6Jbauaqmr3SSliEnpYkDKeRxs9GsPM8YBDS/uNjajed81nydLv7bmMx5V+pXmu++S1hUl6Rg17WB4O/1PpkyuvGiwiYiI+sUucSIiIg8UUVoXB04hIiLygPdX2FHFofO1X23cShr2UFp2iXcnbHlepUKZNF6qy5yftthmStzalm4VXtYMNxrl0KRSWpf54myx5nTb6W/d8LIQs+6zbJSVGzFszXsnve+2cUClz56xHP4cDzICwC7fD2nZdgWUq/tmpPi2WdejC8rMsEuciIjIA2ywiYiIPMAYNhEREcVJbK+wB2Hg/JpwmRJUUzfbIRL77NMlaCflxtrquowlmRI41cSwzYBrVHXN+i5Dk0pfU7OueQxS/NvhPLmMdWuWS3NZSnUVn1PNvSJmucu4Cpq6mq+dbd3wfj26gMwOu8SJiIg80Au3RtejXzQD5SKWiIhoQOMVNhER+auIbjorqgbbFk+SyjXxI6muddxuoX5UY5Kb++2zXZfBk13m/JQOyqybMHcUjrGaK2vG+JYGvrbVlZZdYtgmTc62y2svSV9mnv8hof96ts9EruayFJa1Y4lraL6j4vfOA5qpOfOmiGLY7BInIiLyQFFdYRMR0QDDLnF/5WrovzjKtgvfpW6fFWx9i9nmwFj7LKUuZlv3c4lQJi2bXcYu65rCdc1u7Vwdk3lSHc5pSffXf7v0N7uM1SvQfsY1YS0qMHaJExERUZwMuCtsIiIqIkV0hc0Gm4iI/MUYNlEO2GKbqgCN9NHVfKxtdTX7sS1nymW7Lq9HU9cs6+63FoCI3/f8i/nhEUc6IyIiojjhFTYREfmrB26XnoxhExER5QFj2EQ5YH4xzF+2qi+ONHynbWhPTV3NfqRlzVdNs13butryqPYTEun7nn8xPzwqImywiYjIX+wSJyIi8gC7xP0V/rHkMoSg+aPLnGAoDqQfhprPoPrz2pvmb6DvQWl+vYbr2rYrdhPbupB7hDJp2TYTlmZdzWxdmu1q1jVPqsM5Vb13AtvnR/rsCbSfcam+RxdkNMAMuAabiIiKCLvEiYiIPFBEDTYHTiEiIvJAUV1hmz+kzBh3uNz8JSPFtMwyTfzLFlvrSfO3lrTfPtu1xRGlg5Lil+bolQcVB2XWDYwdJcKxXHNlM857IMMyAOhS1HUZxlSKF2uOv8tYdnntPenLzPMffn9snwnzvQx/LjSfJ8XnVLNZLc13VPzeeSCWxxzA7caxIKoDyb2iarCJiGiA6QGQcFzfE2ywiYjIX0XUYDOGTURE5IHYXmGHQxK+/6rQhFeiSll1CQVag/JSLqwZnwwva2KZ0nb6Wy4Nx1j3GIXScrWlbjjXusxS1yQNTeoSw96TYVl/y9meJ2M/0vshva/9lfcKZVJdxedUk95tlkc1tICN5mtnW9eji0Z3HDiFiIjIA+wSJyIiojjx/grblqqVj/1qhi3VZEwBWfcAqvZjbicwuiETZpdmd5q/+9tRtl2lZk+vua5ZXhruvtV09XYYZeZwo+Fl29fF7MouT/O3jS3dKrxsHr/ttXcIZdKyUSa9P7b3Tnrfzbrm50n67BnL4c+xdgRd6fuhiRCZorqYy1W3vJfYJU5EROQBdokTERFRnPAKm4iI/NULt6tkdolHK1cpXuZ7nO22tXHpbLftEnfTDE1qrluiSeHRxDalUTTNunuN5aHm8r6v/x5sxmLNlb8K/S3FrAH5K2LGrI8ylsMvwiWtyzxRUrz+K2PZjHF/leZvoO95Cm370L7UIrNqeFkz2qtZrvk8WVLANJ9xzfdDk04pffdz+X9DQ/q/UYjtqPXCrUvcowabXeJEREQe8OIKm4iIqF+uXREe3XTGBpuIiPzFBptsMa5wLEGTh62JnZnlUc5kqUmPrtCsLMUgbbHNfUKZVNdcrtplFEpxapeYte0gNTncmv1IMWwzZm2ei68UdUPL0vkG5PdZqmtb1+Gz5zLSbVQzfkp1taFTTTw82+1o68airWMMm4iIiOKEV9hEROQvdokTERF5oIi6xL1rsM1zK/Xp2344meOO9wplEs04xS75oLa60kyW5rImvtcnhq2JS2db14x77jaWxbC0EYUvtcW0JeF4su0FSFNzRjm9pjSW+FeW5f8J/S3ErAGgO3QezfNv5mGHyzUxa3PZpe7B9Iu2+zSk74tteHPNdzTb/G6T7X+bR22Ql9avX4+f/OQn2LhxI3bu3ImXXnoJl112mbjOunXr0NTUhC1btqCurg5z587F7NmzVftlDJuIiPx1ZKSzbB9Z/LrZu3cvTj/9dDzxxBMZ1W9pacHFF1+M8847Dx988AHuvvtu3HbbbVi5cqVqv95dYRMRESX1AAgc1s+iwZ45cyZmzpyZcf2nnnoKY8aMweLFiwEAp5xyCt577z089thjuOKKKzLezoBrsMNdRbZubWlqTtuwpdK2bd3R0jFoUkakrjqXmSzN5WHGxkqyHUK0z3CixrKUbWUum120Q9L8DQBHd6YuD25D5qQXa3aBm0OThg9aE2Qx32lNWtdXxrKUumWch0PGeQpvyijqc/73CWW27vS9GZYB4tvRY3xONdli0vdD873ThrykutryqNfrj0f3Z6l1dqZ+yMvKylBWVhbJtt9++23MmDEj5bmLLroIzzzzDA4ePIghQzJLDmaXOBER+as3ggeA0aNHo7q6Ovlobm6O7BDb2tpQU1OT8lxNTQ0OHTqEXbvMe0nSUzfYO3bswPXXX48RI0Zg6NChOOOMM7Bx48ZkeRAEmD9/Purq6lBRUYFp06Zhy5Yt2t0QERHZucSvjzwAtLa2oqOjI/mYN29epIeZSKTeyh4EQb/PS1QN9pdffolzzjkHQ4YMwW9+8xt89NFH+OlPf4rhw4cn6zz66KNYtGgRnnjiCbz77ruora3FhRdeiN27zX4xIiKieKiqqkp5RNUdDgC1tbVoa0sNQbW3t2Pw4MEYMWJExttRxbAfeeQRjB49Gs8++2zyuRNOOCH5dxAEWLx4Me655x5cfvnlAIDnnnsONTU1WLFiBW666SbN7jKiSfOyCcdnNNsx4zpSNEKTAmYu29JLwuUuQ5NaY9pS+tUwY3lfmr8BOVXLFsM2T3J4udQoM9/M4X/5+u/BmhQqW8xaimHna2hSM+grpG5JMWsg9f2xxaWlulIKGCB/RsypOIXha6V3zmVoUs33zmUKXJNmWFOXlDAp1cwLBbjpTGvy5Ml45ZVXUp5bvXo1Jk2alHH8GlC2by+//DImTZqEK6+8EiNHjkRDQwOefvrpZHlLSwva2tpSgutlZWWYOnUqNmzY0O82u7q60NnZmfIgIiLKSEQxbI09e/Zg06ZN2LRpE4DDbd+mTZuwfft2AMC8efNw4403JuvPnj0b27ZtQ1NTEz7++GP84he/wDPPPIM77rhDtV9Vg/3FF19gyZIlGDt2LP7t3/4Ns2fPxm233Ybnn38eAJKX/P0F183ugCOam5tTAv2jR49WvQAiIipiBcjDfu+999DQ0ICGhgYAQFNTExoaGvDjH/8YALBz585k4w0A9fX1eO2117B27VqcccYZePDBB/H444+rUroAZZd4b28vJk2ahIULFwIAGhoasGXLFixZsiTl10R/wfV0gfV58+ahqakpudzZ2clGm4iIYmvatGnJm8b6s2zZsj7PTZ06Fe+//77TflUN9qhRo3DqqaemPHfKKackR2upra0FcPhKe9SoUck67e3tfa66j0iX65YuLKHJZo2SFDMyuymiHJpUkw8a1ayXZnp0n/JQXDFhy8MOv7Wa3GozDm2GeaTEeNuHJHzShxshmFIpEmqLWZvL4QNxGZrUfOelGLa5LAw3+pVR1Ywth8s1I6Da4t3SvQyKPOygK20RAN1YAy7fJc19JvkamjTbui7rSkO05pTrWOIu8e88U3WJn3POOdi6dWvKc59++imOP/54AIcv+2tra7FmzZpkeXd3N9atW4cpU6ZEcLhEREQhEaV1+UB1hf3DH/4QU6ZMwcKFC/G3f/u3eOedd7B06VIsXboUwOGu8MbGRixcuBBjx47F2LFjsXDhQgwdOhTXXnttTl4AERFRMVA12GeeeSZeeuklzJs3Dw888ADq6+uxePFiXHfddck6c+fOxf79+zFnzhx8+eWXOOuss7B69WpUVlZGfvA22h9O4Q5Ms0snqtm7bN1eZpdaqVDmkqpVLpTZlsM9msPMLktxFi1LXWl4UU2XuMk8yVLOTqWREDf0T1//PdgMFkhpXOZylGldQqrZIaO/Wep+tg03Gu4G/8ooM5c125XSvszjNZf3pi9ymNhLLNd878zvvmYY0yiHNc1Xl3QsZgXrQdF0iavHEv/ud7+L7373u2nLE4kE5s+fj/nz57scFxERkR1j2ERERBQnA262LiIiKiLsEo8vaUpMIDWmou0+kKbmlKbbtNWV0h1scapwRNUM42pCs1LMzoy8mrFBc78pYWkjtabEjGlLcWkpdcs8qS6zU0o3CtiCmeGw9TDjzJQby+brS4Q/JQ5pXYHxKZHylczzL8SArbHlr9L8DchpXmZZhMOa9oQ+b5qXrolZm8vmML8uQ5O6/C/QpIRluh2z3HbfTyxi1qYiarDZJU5EROQB766wiYiIkgJ4dZXsgg02ERF5y3XsE4/GTfG/wbbFtKW6pnB8QIpZA3I8SRqq1IyHmXXNMKi0HynuZsbopLRmW10zxh2OFZp1q8wYpEtudaZlJltwUDNmpZQjbHs9Q0I7LjHfeXM5xGUMWmEoTwCpMWLbdKdfZVhmLktlgJynbZmKc2+av4G+M3FKediaZZehSc2PXvhd18ShzXJbLDmqRiiWMWtDMTXYjGETERF5wPsrbCIiKl5ZTmmdsr4vYttgpzuJLj2jtnWlN07KMtKkmml6a4HUYzZ7XM1O1SFCmZTWImVi2cr7ZGYZL6gi3KXpkqplknJTNDlt5oi50hiuZrewNLQqkPr6XD64Uhe5y7iymnQrsxtbSuv6Sijrb10hrWu/cS72pfkbkKMB2tm6pBRJ6btl+971ZlgGZD9Tlq1M2xUvSbduPhtBdokTERFRrMT2CpuIiMiGXeJEREQeKKYuce8abFtcWhpe1IWUumWLK4TXlY63v2Up5uUyvaaU1qWJaVtnvQzl2pSaMdJcvUG2oKMU3DTj0kOFMttQq9L4tZo8QSk3yPZGS3NQWoYBFYcMdRnGVNhWd1faIgByWpf0OTbr5mt6Ten77TK9pm3dbBuhfKWLUXa8a7CJiIiO6IXbDwl2iRMREeVBMcWweZc4ERGRB4rqCtvll5Qmt1ozFacZD5NSlW3DmIbjcH1iycK6Zl0zFGuGQUvS/G1bHr4vfZmalMBuC1BqYtjh5TKjTJoetL/lMNucrGGauR+l8TnNZcswoKqc7XCeti1m/VXqYk9oP9IhmOVmmbkcPhUuQ5O65GxLw5hq86Ol/zkSl/i2DzFr3nRGRETkATbYREREHmAMm4iIiGLF+ytsKS9bM/WmjdRtYv7qkeLHmqk4gdT4mFlXin9rhtM2z4sZs5Neny1WLhFj2rbAmxRINF+ANF64GfiUYti2PGzbzQBSXelnvhmnlpLzNcFZM4YtBYxt8W5FznaPsZ+vMl9VzK2W8rJtoX3NmAbSsuZjq8m77m9ZKtNcNebi6pRjieeG9w02EREVL3aJExERUazE9gq7F0DQz/O2LtfwryWX6TRtwsdhmzJTGppUk+YljVBpbtu2H82omSYprctFZairtFTbJa6ZcjI83KimS9yWtiVNr6k5UZrXrol9mMu2LvG9QpnUV22U2YYblXrTzVk9pek1pcN3Sesyv2eaCIU0haZ2aFLN6LXSMKYSW904Xo1ypDMiIiIPFFMMm13iREREHuAVNhEReauYbjrzrsF2SdWKS5qXdAxmXSmtS4pTa1K1bMegjXFnK3xOK424Z4Umt8Z8sWZOjxSblVK3zDJpvFdAPskSW4BSyg2y5StJMWzNeRLSvPb3pi3qd1Phcilmbda1Ta8pDS9qW9bcJtCd5m9AjlNr07iyjUubbPuV6sYRu8SJiIgoVry7wiYiIjqimK6w2WATEZG3GMP2iBSXloYt7W/dbElTb/a3HGbGw6Rta6bitE17KQ1Nmi+a+N5B44lhxhSOJVLA0gyEavKwh6T5G+g79KgU/NcEn2xBxnC5GTTV5GXbzpMQ7+4x4t1Zhrut62rqSmnlmpi1uay5fUK63cBclvK5+1uWyjTfJYlL4xXeZ39jaORKMV1hM4ZNRETkAe+vsImIqHgFcOsZyGdvgKsB12CHuzdsXdWa7gWp693WpSLtxzZ7l9RlHlVXtm07UXXDaIZTlFJn+isvD3XRDjW6axNlRuVw36k0FClQXF3ixnIQOo9SbzmQfVe1y7ouaV222bvCy5oucbMsV2ldLgZa9zG7xImIiChWBtwVNhERFY9iusJmg01ERN5iWlcM9ABI9PN8Loci1Qw3qvlVJu3XNgyoNDWnGb6UtiMNa5rLuEi2UwJqZ40cJpSVd6VfLjGDptKUmS4xbBdRxrBDyz1GmSYGLI1qaosta2LY0n41WWlmXduUmdJIt9Lwoy7TXmpuXdBMzekyZaZPV5/FILYNNhERkQ27xImIiDxQTA027xInIiLygHdX2FFOkSnlZUf5q0szfZ0Ul9bI1y8x2w0b2b72YUZZn7xroVxKpTbLy40NDzGWw2Fq85wmpOk0gbzEsIODaYsA6GKz0nCdLvFi2yye0oyf0n41U2ba6prHJE2ZqcnDloYf1cSszXJbXclAi1nzpjMiIiIP9CLa+cHjjA02ERF5i1fYHtF0kdvqht+4KLuUNeli0rq2Htfw8dtmAStEl7mtuzbcVa2Z9QhI7eK0jTa6TyiTsrr6pMoZByHNkBblyKSa4S2l82absUrqUtbUtXVzdwllLt3cmuOXur1ts2hFldaVqzCcS9hKs58wn8bn9on3DTYRERWvYrpLnA02ERF5q5gabKZ1EREReWDAXWFL02tKdW31ffoVlmuaIRI1MTspBm+LS4eX+wxNaixrRhsNl9vuIbANB5stl5ipFKvV3BdgGxo223i3bV1pWUrFMuvajsElVUszveZBoSzKexdyJY7/B3nTGRERkQfYJU5ERESxwitsIiLyVjFdYce2we7F19NrZtsN4DKMqTRsaX/bjopmKs5cscV0pHIpZ1WTWy0NPdrfspSHbcYrw3Fp2wyZUh62LWadj9k1bfcFmOcp25k5NXFpTczath9NXfP1dAllillIretmm4etiVn3V1+S7fCjPjVeRwRwi0P7lDPOLnEiIiIPxPYKm4iIyIZd4kRERB5gWlfMSCdU06dvi2lLOdy2mLZUV0P6tWcbHzxfNPExKb5n5kCHy824oRlrNs9FeFtSzNrclm2GTGl6TVtMW9quRDP2sy3OKU3v6JJfLMWpbTFfTVzaJa6uyTnXrCvdJ6CZMtMlZu1SV1sed8V0hc0YNhERkQfYYBMRkbd6Inhk48knn0R9fT3Ky8sxceJEvPnmm2nrrl27FolEos/jk08+Ue3Tiy5xiaar2iR1kWum4jS5pPNoPjy56iJ3SSfRDKdodkf3CmW27lupm1szZaZZtyTN34BbGpe5brbnWJsaJA2NmavhOTVd4pr0K9t2w9vSDNlqlmtS5QbClJka6Y4jn6lShYhhv/jii2hsbMSTTz6Jc845Bz//+c8xc+ZMfPTRRxgzZkza9bZu3Yqqqqrk8rHHHqvaL6+wiYiIFBYtWoS/+7u/w9///d/jlFNOweLFizF69GgsWbJEXG/kyJGora1NPkpKdJd3qgb70KFDuPfee1FfX4+KigqceOKJeOCBB9Db+/VvlCAIMH/+fNTV1aGiogLTpk3Dli1bVAdFRESUiai6xDs7O1MeXV3m9DKHdXd3Y+PGjZgxY0bK8zNmzMCGDRvEY21oaMCoUaMwffp0vPHGG+rXqmqwH3nkETz11FN44okn8PHHH+PRRx/FT37yE/zsZz9L1nn00UexaNEiPPHEE3j33XdRW1uLCy+8ELt371YfHBERkaQXbo31kcvN0aNHo7q6Ovlobm7ud3+7du1CT08PampqUp6vqalBW1tbv+uMGjUKS5cuxcqVK7Fq1SqMGzcO06dPx/r161WvVRXDfvvtt/E3f/M3uOSSSwAAJ5xwAn71q1/hvffeA3D46nrx4sW45557cPnllwMAnnvuOdTU1GDFihW46aabVAeXjXA8QtvfL6V1+TAVZ1QxbZf7AFymBAwva6a9BFKHobTFsEsyLANSz4VLDFvT8aVJw9GkEZnlmnQlTd1cDgOabVzdlu4mvT4pZm3W1Zx/7fSaklzVNfmUs6zV2tqaEl8uKysT6ycSiZTlIAj6PHfEuHHjMG7cuOTy5MmT0draisceewznn39+xseoatPOPfdc/Pa3v8Wnn34KAPjd736Ht956CxdffDEAoKWlBW1tbSldBWVlZZg6dWraroKurq4+XRFERESZ6I3gAQBVVVUpj3QN9jHHHIOSkpI+V9Pt7e19rrolZ599Nj777LOM6wPKBvvOO+/ENddcg5NPPhlDhgxBQ0MDGhsbcc011wBA8gVougqam5tTuiFGjx6tegFERFS88p3WVVpaiokTJ2LNmjUpz69ZswZTpkzJeDsffPABRo0apdq3qkv8xRdfxPLly7FixQqMHz8emzZtQmNjI+rq6jBr1qxkPU1Xwbx589DU1JRc7uzsZKNNRESx1dTUhBtuuAGTJk3C5MmTsXTpUmzfvh2zZ88GcLhd27FjB55//nkAwOLFi3HCCSdg/Pjx6O7uxvLly7Fy5UqsXLlStV9Vg/2jH/0Id911F66++moAwGmnnYZt27ahubkZs2bNQm1tLYDDV9rhXw5SV0FZWZk1VpAtlxxtG+lXWb6m4pS4xLNt8dZsh1eUhhMF5Fxqc10pf9oWl5bysKW4tK2uKVfTa4a5TNFoyy+W4rgudaOKNWvi6prXajumqO7T0Mass40fa/7/+BijLkQe9lVXXYU///nPeOCBB7Bz505MmDABr732Go4//ngAwM6dO7F9+/Zk/e7ubtxxxx3YsWMHKioqMH78eLz66qvJcHKmVA32vn37MGhQ6r+pkpKSZFpXfX09amtrsWbNGjQ0NCQPdN26dXjkkUdUB0ZERGRTqLHE58yZgzlz5vRbtmzZspTluXPnYu7cuVnu6WuqBvvSSy/FQw89hDFjxmD8+PH44IMPsGjRIvzgBz8AcLgrvLGxEQsXLsTYsWMxduxYLFy4EEOHDsW1117rfLBERERhhWqwC0HVYP/sZz/D//t//w9z5sxBe3s76urqcNNNN+HHP/5xss7cuXOxf/9+zJkzB19++SXOOussrF69GpWVlZEfvJami1yTxqXZT1xSvqRuIFsXUbazCElDkZrLttm6pC5z832NarjRgd4l7tItLKVm2YYFlbrTNbONRdklnm3Yx6WuKV8zcPnYDV6sEkEQ5HPYV6vOzk5UV1fjKAD936YWHU1M2yXHVtpPlA1GtnU1MV/butL0lJr9SNuxbYsNdvplNtiZrcsG200AYA+Ajo6OlNzmKB1pK+YBKHfYzgEAzcjtsUbF+8k/iIioeB0Z6cxlfV9w8g8iIiIPxPYKuwdfd4lH1bVoiiqmbYt3S8OlxvGGB013m1lXilPbzpM0NKwt3UfqppdSwmzDi0pDk5py9Tk1Se+PZrhLTXe6pgtZ01Wt3Y/UzS19RlyGcM3XcKMuU2YWugu8kHjTGRERkQcKkYddKOwSJyIi8gCvsImIyFvsEo8Z6YRGGTfUTM0pxVulWK0tbh7Vh8cWO9OcNylerMnZtsW7pXNqy8sOn0czDhpV+pvJJZ1PI8r0nqhSkFzixVGlj+VqGNb+jkOqm+1wo4xZR4Nd4kRERBQrXlxhExER9Ydd4h7JVXd5lMOYatKV4pD2ZTsmzWhN2b5222hrmnU1I5S5pHXZtpUtzWcgyrSibNOVtN3Nmv1km6ql6T7vrzzb/Zh8Hr2sP+mOKZ/DZ7LBJiIi8kAAtx8osRqb24IxbCIiIg/wCpuIiLzFLvEBwiWVyVSIYUzN/fjwwdKkgGnS32znTRPDTreerW6+0rhsooqDatKKcjXkqVk/ylQzlyFDXWZIk+pmul4m5Znux4UP/3OKqcFmlzgREZEHBvQVNhERDWzFNHAKG2wiIvJWMXWJF1WDnauYdq6GMbXtJ1cfNJdYsyZnW3qttrpRTZlpKqYYtm272Q656RIvzmVudabHoF030/Vs6zJmTTZF1WATEdHAwi5xIiIiDxRTlzjvEiciIvJAUV9hRxXTztW442a5ZtxxW353rkixQCm+ba6riVm7rGsek1TXLDPFIYZtiiqfOFdjb5vlcTmmTLfjsq72yi6q77BPV5SZ6IXba2KXOBERUR4whk1EROSBHrj1bvnU48AGO0TqVtUo1NScklx1kWteq2ZdbXqM1HVtCztIdaVu8Kimz3Sh/WeTry7yqOpGuW6m2zHXjSoEkUl5GLvAycQGm4iIvMUrbCIiIg8UUwybaV1EREQe4BV2GnEYxtTcr63rRkoBKxTNcKNS3NA2LKtmCFfpvbTFv6W6vqV1aeq6TCmZqyk/bWWFSNVyqRvld9anbl5X7BInIiLyALvEiYiIKFZ4hU1ERN7iSGfURyGGMTX3a9unJmc7vF/NtJeF4jINqUl6PS5TmOYrLztX/5wKNW2kdJ+DZluMWRenHgAJx/V9wS5xIiIiD/AKm4iIvFVMN52xwc5SHIcx1dSNw4c0qpnJzHVdus9dusALcU5zGbvLtps7X93nmnW1703cu8F96sbNtWLqEmeDTURE3iqmBpsxbCIiIg/wCpuIiLzFGDap5GoYUyB/MW3NMUSVPpYvUQ33avLhi16IWGy+pph0iXdrt5VtXaZu5R67xImIiChWeIVNRETeCuDWkxFEdSB5wAabiIi85dql7VOXOBvsmHOZmjNMmsrSVjdXcWmXHOh8iSrfPl/yNaaySy61Zr/5ine71s/kGIhcscEmIiJv8QqbiIjIA71wu0vcp14RNtg5EGWaV5h2pq+wfA1r6pIClul2KD2Xfz5xSG3yIVXLlKt/+PzMk4kNNhEReYtd4kRERB5gg01EROQBxrApUr6lBpmiGtZUMwyo7UukiXFrUtqkdX36Yh/hcvVQqDh1PrYbR74fP+UeG2wiIvKW6w9pn36Is8EmIiJvFVODzck/iIiIPMArbI9FlZetjatr1s12aFXNdm3b1sS3NevaRHW/Qq5im3EYxjRf241yXZ+uyIpBD9wm8PDp/WSDTURE3iqmBptd4kRERB7gFTYREXmrmG46Y4OdZ7kaZxzQxYvDcnlMYS5x5yhj5dJ+TC7nIg55tVEdQy5jy/n4hxnHmHUcPh8DAbvEiYiIKFZ4hU1ERN7qhdsVtsu6+cYGm5zkqztdwyXdzeTbsLJRdrPGsauQ3chkch1LnA02ERFRHvSgeBpsxrCJiIg8ELsr7CA4/HvHp189LnL1OgsV08lVN6rLL+ioxOEYbKI8/9luS7uepn62n804xjgH8v+4I6/tyP/zXCqmK+zYNdi7d+8GABwo8HEQEZGb3bt3o7q6OifbLi0tRW1tLdra2py3VVtbi9LS0giOKrcSQT5+Ain09vbiT3/6E4IgwJgxY9Da2oqqqqpCH1ZsdXZ2YvTo0TxPFjxPmeF5ygzPkywIAuzevRt1dXUYNCh3kdcDBw6gu7vbeTulpaUoLy+P4IhyK3ZX2IMGDcJxxx2Hzs5OAEBVVRW/EBngecoMz1NmeJ4yw/OUXq6urMPKy8u9aGijwpvOiIiIPMAGm4iIyAOxbbDLyspw3333oaysrNCHEms8T5nhecoMz1NmeJ6oEGJ30xkRERH1FdsrbCIiIvoaG2wiIiIPsMEmIiLyABtsIiIiD8S2wX7yySdRX1+P8vJyTJw4EW+++WahD6lgmpubceaZZ6KyshIjR47EZZddhq1bt6bUCYIA8+fPR11dHSoqKjBt2jRs2bKlQEccD83NzUgkEmhsbEw+x/N02I4dO3D99ddjxIgRGDp0KM444wxs3LgxWc7zBBw6dAj33nsv6uvrUVFRgRNPPBEPPPAAenu/Hv2c54nyKoihF154IRgyZEjw9NNPBx999FFw++23B8OGDQu2bdtW6EMriIsuuih49tlng//+7/8ONm3aFFxyySXBmDFjgj179iTrPPzww0FlZWWwcuXKYPPmzcFVV10VjBo1Kujs7CzgkRfOO++8E5xwwgnBN7/5zeD2229PPs/zFAR/+ctfguOPPz743ve+F/zXf/1X0NLSEvz7v/978Pnnnyfr8DwFwYIFC4IRI0YE//Iv/xK0tLQE//RP/xQcddRRweLFi5N1eJ4on2LZYH/7298OZs+enfLcySefHNx1110FOqJ4aW9vDwAE69atC4IgCHp7e4Pa2trg4YcfTtY5cOBAUF1dHTz11FOFOsyC2b17dzB27NhgzZo1wdSpU5MNNs/TYXfeeWdw7rnnpi3neTrskksuCX7wgx+kPHf55ZcH119/fRAEPE+Uf7HrEu/u7sbGjRsxY8aMlOdnzJiBDRs2FOio4qWjowMAcPTRRwMAWlpa0NbWlnLOysrKMHXq1KI8ZzfffDMuueQSXHDBBSnP8zwd9vLLL2PSpEm48sorMXLkSDQ0NODpp59OlvM8HXbuuefit7/9LT799FMAwO9+9zu89dZbuPjiiwHwPFH+xW7yj127dqGnpwc1NTUpz9fU1EQyjZrvgiBAU1MTzj33XEyYMAEAkuelv3O2bdu2vB9jIb3wwgt4//338e677/Yp43k67IsvvsCSJUvQ1NSEu+++G++88w5uu+02lJWV4cYbb+R5+j933nknOjo6cPLJJ6OkpAQ9PT146KGHcM011wDg54nyL3YN9hGJROqU5EEQ9HmuGN1yyy348MMP8dZbb/UpK/Zz1traittvvx2rV68WZ/Ap9vPU29uLSZMmYeHChQCAhoYGbNmyBUuWLMGNN96YrFfs5+nFF1/E8uXLsWLFCowfPx6bNm1CY2Mj6urqMGvWrGS9Yj9PlD+x6xI/5phjUFJS0udqur29vc8v2WJz66234uWXX8Ybb7yB4447Lvl8bW0tABT9Odu4cSPa29sxceJEDB48GIMHD8a6devw+OOPY/DgwclzUeznadSoUTj11FNTnjvllFOwfft2APw8HfGjH/0Id911F66++mqcdtppuOGGG/DDH/4Qzc3NAHieKP9i12CXlpZi4sSJWLNmTcrza9aswZQpUwp0VIUVBAFuueUWrFq1Cq+//jrq6+tTyuvr61FbW5tyzrq7u7Fu3bqiOmfTp0/H5s2bsWnTpuRj0qRJuO6667Bp0yaceOKJPE8AzjnnnD5pgZ9++imOP/54APw8HbFv3z4MGpT6L7KkpCSZ1sXzRHlXwBve0jqS1vXMM88EH330UdDY2BgMGzYs+MMf/lDoQyuIf/iHfwiqq6uDtWvXBjt37kw+9u3bl6zz8MMPB9XV1cGqVauCzZs3B9dccw3TS4Ig5S7xIOB5CoLDKW+DBw8OHnrooeCzzz4LfvnLXwZDhw4Nli9fnqzD8xQEs2bNCr7xjW8k07pWrVoVHHPMMcHcuXOTdXieKJ9i2WAHQRD84z/+Y3D88ccHpaWlwbe+9a1kClMxAtDv49lnn03W6e3tDe67776gtrY2KCsrC84///xg8+bNhTvomDAbbJ6nw1555ZVgwoQJQVlZWXDyyScHS5cuTSnneQqCzs7O4Pbbbw/GjBkTlJeXByeeeGJwzz33BF1dXck6PE+UT5xek4iIyAOxi2ETERFRX2ywiYiIPMAGm4iIyANssImIiDzABpuIiMgDbLCJiIg8wAabiIjIA2ywiYiIPMAGm4iIyANssImIiDzABpuIiMgDbLCJiIg88P8Bnt02AtQfWusAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.colors import LogNorm\n",
    "# Plot a slice at z = nz // 2\n",
    "plt.imshow(final_temperature_cpu[:, :, 0], cmap='hot')#, vmin = 0, vmax = 1)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "b36e9e4b-011e-4faf-8622-9f5ee52fbbc8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# Grid parameters\n",
    "nr, ntheta, nz = 100, 100, 100  # Number of grid points in r, theta, z\n",
    "r_max = 10.0                    # Maximum radius\n",
    "z_max = 10.0                    # Maximum height\n",
    "\n",
    "dr = r_max / nr                 # Radial grid spacing\n",
    "dtheta = 2 * np.pi / ntheta     # Angular grid spacing\n",
    "dz = z_max / nz                 # Axial grid spacing\n",
    "\n",
    "dt = 0.01                       # Time step\n",
    "alpha = 0.01                    # Thermal diffusivity\n",
    "nt = 10                       # Number of time steps\n",
    "\n",
    "# Create radial and angular coordinate arrays\n",
    "r = np.linspace(0 + dr/2, r_max - dr/2, nr)  # Avoid r=0 to prevent division by zero\n",
    "theta = np.linspace(0, 2 * np.pi, ntheta, endpoint=False)\n",
    "z = np.linspace(0, z_max, nz)\n",
    "\n",
    "def second_derivative_z(u_z):\n",
    "    \"\"\"\n",
    "    Compute the second derivative along the z-axis using central differences.\n",
    "    Assumes Neumann boundary conditions (zero gradient at boundaries).\n",
    "    \"\"\"\n",
    "    d2u_dz2 = np.zeros_like(u_z)\n",
    "    d2u_dz2[1:-1] = (u_z[2:] - 2 * u_z[1:-1] + u_z[:-2]) / dz**2\n",
    "    # Neumann boundary conditions\n",
    "    d2u_dz2[0] = (u_z[1] - u_z[0]) / dz**2\n",
    "    d2u_dz2[-1] = (u_z[-2] - u_z[-1]) / dz**2\n",
    "    return d2u_dz2\n",
    "\n",
    "def main():\n",
    "    # Initialize the temperature field\n",
    "    u = np.zeros((nr, ntheta, nz))\n",
    "    T0 = 100.0  # Initial temperature at z=0\n",
    "    u[:, :, 0] = T0  # Uniform temperature at z=0 across r and theta\n",
    "\n",
    "    for _ in range(nt):\n",
    "        # Compute second derivative along the z-axis using apply_along_axis\n",
    "        d2u_dz2 = np.apply_along_axis(second_derivative_z, axis=2, arr=u)\n",
    "\n",
    "        # Compute second derivative along the r-axis\n",
    "        d2u_dr2 = np.zeros_like(u)\n",
    "        for i in range(1, nr - 1):\n",
    "            u_plus = u[i + 1, :, :]\n",
    "            u_minus = u[i - 1, :, :]\n",
    "            u_i = u[i, :, :]\n",
    "            d2u_dr2[i, :, :] = ((u_plus - 2 * u_i + u_minus) / dr**2 +\n",
    "                                (u_plus - u_minus) / (2 * dr * r[i]))\n",
    "        # Boundary at r=0 (special case)\n",
    "        i = 0\n",
    "        u_plus = u[i + 1, :, :]\n",
    "        u_i = u[i, :, :]\n",
    "        d2u_dr2[i, :, :] = (2 * (u_plus - u_i) / dr**2)\n",
    "\n",
    "        # Boundary at r=r_max (Neumann boundary condition)\n",
    "        i = nr - 1\n",
    "        u_minus = u[i - 1, :, :]\n",
    "        u_i = u[i, :, :]\n",
    "        d2u_dr2[i, :, :] = (2 * (u_minus - u_i) / dr**2)\n",
    "\n",
    "        # Compute second derivative along the theta-axis\n",
    "        d2u_dtheta2 = np.zeros_like(u)\n",
    "        r_squared = r[:, None] ** 2  # Shape: (nr, 1)\n",
    "        for j in range(ntheta):\n",
    "            jp = (j + 1) % ntheta  # Periodic boundary condition\n",
    "            jm = (j - 1) % ntheta\n",
    "            u_jp = u[:, jp, :]\n",
    "            u_jm = u[:, jm, :]\n",
    "            u_j = u[:, j, :]\n",
    "            numerator = (u_jp - 2 * u_j + u_jm) / dtheta**2  # Shape: (nr, nz)\n",
    "            d2u_dtheta2[:, j, :] = numerator / r_squared  # Broadcast over nz\n",
    "\n",
    "        # Update the temperature field\n",
    "        u += alpha * dt * (d2u_dr2 + d2u_dtheta2 + d2u_dz2)\n",
    "\n",
    "    # Optionally, return or save the final temperature field\n",
    "    return u\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    final_temperature = main()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c0bf4a09-c0e8-4571-8ee7-e5d331bd19b3",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'u' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[1;32mIn[5], line 4\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mmatplotlib\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mpyplot\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mplt\u001b[39;00m\n\u001b[0;32m      3\u001b[0m z_index \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m99\u001b[39m  \u001b[38;5;66;03m# Middle of the z-axis\u001b[39;00m\n\u001b[1;32m----> 4\u001b[0m plt\u001b[38;5;241m.\u001b[39mimshow(u[:, :, z_index], extent\u001b[38;5;241m=\u001b[39m(\u001b[38;5;241m0\u001b[39m, r_max, \u001b[38;5;241m0\u001b[39m, \u001b[38;5;241m2\u001b[39m \u001b[38;5;241m*\u001b[39m np\u001b[38;5;241m.\u001b[39mpi), aspect\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mauto\u001b[39m\u001b[38;5;124m'\u001b[39m, cmap\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mhot\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[0;32m      5\u001b[0m plt\u001b[38;5;241m.\u001b[39mcolorbar(label\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mTemperature\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[0;32m      6\u001b[0m plt\u001b[38;5;241m.\u001b[39mxlabel(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mRadius (r)\u001b[39m\u001b[38;5;124m'\u001b[39m)\n",
      "\u001b[1;31mNameError\u001b[0m: name 'u' is not defined"
     ]
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "z_index = 99  # Middle of the z-axis\n",
    "plt.imshow(u[:, :, z_index], extent=(0, r_max, 0, 2 * np.pi), aspect='auto', cmap='hot')\n",
    "plt.colorbar(label='Temperature')\n",
    "plt.xlabel('Radius (r)')\n",
    "plt.ylabel('Angle (θ)')\n",
    "plt.title(f'Temperature Distribution at z={z[z_index]:.2f}')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "104c38f9-abf8-4eff-9b41-60d40365d794",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The memory_profiler extension is already loaded. To reload it, use:\n",
      "  %reload_ext memory_profiler\n"
     ]
    }
   ],
   "source": [
    "%load_ext memory_profiler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c74ca72e-0043-48bc-b726-84ec1832cea3",
   "metadata": {},
   "outputs": [
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "Cell \u001b[1;32mIn[4], line 1\u001b[0m\n\u001b[1;32m----> 1\u001b[0m get_ipython()\u001b[38;5;241m.\u001b[39mrun_cell_magic(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mmemit\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mimport cupy as cp\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m# Grid parameters\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124mnr, ntheta, nz = 100, 100, 100  # Number of grid points in r, theta, z\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124mr_max = 10.0                    # Maximum radius\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124mz_max = 10.0                    # Maximum height\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124mdr = r_max / nr                 # Radial grid spacing\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124mdtheta = 2 * cp.pi / ntheta     # Angular grid spacing\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124mdz = z_max / nz                 # Axial grid spacing\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124mdt = 0.01                       # Time step\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124malpha = 0.01                    # Thermal diffusivity\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124mnt = 100                        # Number of time steps\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m# Create radial and angular coordinate arrays\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124mr = cp.linspace(0 + dr/2, r_max - dr/2, nr)  # Avoid r=0 to prevent division by zero\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124mtheta = cp.linspace(0, 2 * cp.pi, ntheta, endpoint=False)\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124mz = cp.linspace(0, z_max, nz)\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124mdef second_derivative_z(u_z):\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m    \u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m    Compute the second derivative along the z-axis using central differences.\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m    Assumes Neumann boundary conditions (zero gradient at boundaries).\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m    \u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m    d2u_dz2 = cp.zeros_like(u_z)\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m    d2u_dz2[1:-1] = (u_z[2:] - 2 * u_z[1:-1] + u_z[:-2]) / dz**2\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m    # Neumann boundary conditions\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m    d2u_dz2[0] = (u_z[1] - u_z[0]) / dz**2\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m    d2u_dz2[-1] = (u_z[-2] - u_z[-1]) / dz**2\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m    return d2u_dz2\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124mdef main():\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m    # Initialize the temperature field\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m    u = cp.zeros((nr, ntheta, nz))\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m    T0 = 100.0  # Initial temperature at z=0\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m    u[:, :, 0] = T0  # Uniform temperature at z=0 across r and theta\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m    r_squared = r[:, cp.newaxis] ** 2  # Shape: (nr, 1)\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m    for _ in range(nt):\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m        # Compute second derivative along the z-axis using apply_along_axis\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m        d2u_dz2 = cp.apply_along_axis(second_derivative_z, axis=2, arr=u)\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m        # Compute second derivative along the r-axis\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m        d2u_dr2 = cp.zeros_like(u)\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m        u_plus = cp.roll(u, -1, axis=0)\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m        u_minus = cp.roll(u, 1, axis=0)\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m        u_i = u\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m        # Interior points\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m        d2u_dr2[1:-1, :, :] = ((u_plus[1:-1] - 2 * u_i[1:-1] + u_minus[1:-1]) / dr**2 +\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m                               (u_plus[1:-1] - u_minus[1:-1]) / (2 * dr * r[1:-1, cp.newaxis, cp.newaxis]))\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m        # Boundary at r=0\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m        d2u_dr2[0, :, :] = (2 * (u[1, :, :] - u[0, :, :]) / dr**2)\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m        # Boundary at r=r_max\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m        d2u_dr2[-1, :, :] = (2 * (u[-2, :, :] - u[-1, :, :]) / dr**2)\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m        # Compute second derivative along the theta-axis\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m        d2u_dtheta2 = cp.zeros_like(u)\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m        u_jp = cp.roll(u, -1, axis=1)\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m        u_jm = cp.roll(u, 1, axis=1)\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m        numerator = (u_jp - 2 * u + u_jm) / dtheta**2  # Shape: (nr, ntheta, nz)\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m        d2u_dtheta2 = numerator / r_squared[:, :, cp.newaxis]  # Broadcast over nz\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m        # Update the temperature field\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m        u += alpha * dt * (d2u_dr2 + d2u_dtheta2 + d2u_dz2)\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m    # Optionally, return or save the final temperature field\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m    return u\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124mif __name__ == \u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m__main__\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m:\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m    final_temperature = main()\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m    # Transfer data back to CPU for further processing or visualization\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m    final_temperature_cpu = cp.asnumpy(final_temperature)\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m'\u001b[39m)\n",
      "File \u001b[1;32m~\\anaconda3\\Lib\\site-packages\\IPython\\core\\interactiveshell.py:2493\u001b[0m, in \u001b[0;36mInteractiveShell.run_cell_magic\u001b[1;34m(self, magic_name, line, cell)\u001b[0m\n\u001b[0;32m   2491\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mbuiltin_trap:\n\u001b[0;32m   2492\u001b[0m     args \u001b[38;5;241m=\u001b[39m (magic_arg_s, cell)\n\u001b[1;32m-> 2493\u001b[0m     result \u001b[38;5;241m=\u001b[39m fn(\u001b[38;5;241m*\u001b[39margs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs)\n\u001b[0;32m   2495\u001b[0m \u001b[38;5;66;03m# The code below prevents the output from being displayed\u001b[39;00m\n\u001b[0;32m   2496\u001b[0m \u001b[38;5;66;03m# when using magics with decorator @output_can_be_silenced\u001b[39;00m\n\u001b[0;32m   2497\u001b[0m \u001b[38;5;66;03m# when the last Python token in the expression is a ';'.\u001b[39;00m\n\u001b[0;32m   2498\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mgetattr\u001b[39m(fn, magic\u001b[38;5;241m.\u001b[39mMAGIC_OUTPUT_CAN_BE_SILENCED, \u001b[38;5;28;01mFalse\u001b[39;00m):\n",
      "File \u001b[1;32m~\\anaconda3\\Lib\\site-packages\\memory_profiler.py:1113\u001b[0m, in \u001b[0;36mMemoryProfilerMagics.memit\u001b[1;34m(self, line, cell)\u001b[0m\n\u001b[0;32m   1111\u001b[0m \u001b[38;5;28;01mwhile\u001b[39;00m counter \u001b[38;5;241m<\u001b[39m repeat:\n\u001b[0;32m   1112\u001b[0m     counter \u001b[38;5;241m+\u001b[39m\u001b[38;5;241m=\u001b[39m \u001b[38;5;241m1\u001b[39m\n\u001b[1;32m-> 1113\u001b[0m     tmp \u001b[38;5;241m=\u001b[39m memory_usage((_func_exec, (stmt, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mshell\u001b[38;5;241m.\u001b[39muser_ns)),\n\u001b[0;32m   1114\u001b[0m                        timeout\u001b[38;5;241m=\u001b[39mtimeout, interval\u001b[38;5;241m=\u001b[39minterval,\n\u001b[0;32m   1115\u001b[0m                        max_usage\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mTrue\u001b[39;00m, max_iterations\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m1\u001b[39m,\n\u001b[0;32m   1116\u001b[0m                        include_children\u001b[38;5;241m=\u001b[39minclude_children)\n\u001b[0;32m   1117\u001b[0m     mem_usage\u001b[38;5;241m.\u001b[39mappend(tmp)\n\u001b[0;32m   1119\u001b[0m result \u001b[38;5;241m=\u001b[39m MemitResult(mem_usage, baseline, repeat, timeout, interval,\n\u001b[0;32m   1120\u001b[0m                      include_children)\n",
      "File \u001b[1;32m~\\anaconda3\\Lib\\site-packages\\memory_profiler.py:379\u001b[0m, in \u001b[0;36mmemory_usage\u001b[1;34m(proc, interval, timeout, timestamps, include_children, multiprocess, max_usage, retval, stream, backend, max_iterations)\u001b[0m\n\u001b[0;32m    376\u001b[0m \u001b[38;5;66;03m# When there is an exception in the \"proc\" - the (spawned) monitoring processes don't get killed.\u001b[39;00m\n\u001b[0;32m    377\u001b[0m \u001b[38;5;66;03m# Therefore, the whole process hangs indefinitely. Here, we are ensuring that the process gets killed!\u001b[39;00m\n\u001b[0;32m    378\u001b[0m \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[1;32m--> 379\u001b[0m     returned \u001b[38;5;241m=\u001b[39m f(\u001b[38;5;241m*\u001b[39margs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkw)\n\u001b[0;32m    380\u001b[0m     parent_conn\u001b[38;5;241m.\u001b[39msend(\u001b[38;5;241m0\u001b[39m)  \u001b[38;5;66;03m# finish timing\u001b[39;00m\n\u001b[0;32m    381\u001b[0m     ret \u001b[38;5;241m=\u001b[39m parent_conn\u001b[38;5;241m.\u001b[39mrecv()\n",
      "File \u001b[1;32m~\\anaconda3\\Lib\\site-packages\\memory_profiler.py:889\u001b[0m, in \u001b[0;36m_func_exec\u001b[1;34m(stmt, ns)\u001b[0m\n\u001b[0;32m    886\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m_func_exec\u001b[39m(stmt, ns):\n\u001b[0;32m    887\u001b[0m     \u001b[38;5;66;03m# helper for magic_memit, just a function proxy for the exec\u001b[39;00m\n\u001b[0;32m    888\u001b[0m     \u001b[38;5;66;03m# statement\u001b[39;00m\n\u001b[1;32m--> 889\u001b[0m     exec(stmt, ns)\n",
      "File \u001b[1;32m<string>:73\u001b[0m\n",
      "File \u001b[1;32m<string>:43\u001b[0m, in \u001b[0;36mmain\u001b[1;34m()\u001b[0m\n",
      "File \u001b[1;32m~\\anaconda3\\Lib\\site-packages\\cupy\\lib\\_shape_base.py:60\u001b[0m, in \u001b[0;36mapply_along_axis\u001b[1;34m(func1d, axis, arr, *args, **kwargs)\u001b[0m\n\u001b[0;32m     58\u001b[0m buff[ind0] \u001b[38;5;241m=\u001b[39m res\n\u001b[0;32m     59\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m ind \u001b[38;5;129;01min\u001b[39;00m inds:\n\u001b[1;32m---> 60\u001b[0m     out \u001b[38;5;241m=\u001b[39m func1d(inarr_view[ind], \u001b[38;5;241m*\u001b[39margs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs)\n\u001b[0;32m     61\u001b[0m     buff[ind] \u001b[38;5;241m=\u001b[39m cupy\u001b[38;5;241m.\u001b[39masarray(out)\n\u001b[0;32m     63\u001b[0m \u001b[38;5;66;03m# restore the inserted axes back to where they belong\u001b[39;00m\n",
      "File \u001b[1;32m<string>:30\u001b[0m, in \u001b[0;36msecond_derivative_z\u001b[1;34m(u_z)\u001b[0m\n",
      "File \u001b[1;32mcupy\\_core\\core.pyx:1594\u001b[0m, in \u001b[0;36mcupy._core.core._ndarray_base.__setitem__\u001b[1;34m()\u001b[0m\n",
      "File \u001b[1;32mcupy\\_core\\_routines_indexing.pyx:51\u001b[0m, in \u001b[0;36mcupy._core._routines_indexing._ndarray_setitem\u001b[1;34m()\u001b[0m\n",
      "File \u001b[1;32mcupy\\_core\\_routines_indexing.pyx:996\u001b[0m, in \u001b[0;36mcupy._core._routines_indexing._scatter_op\u001b[1;34m()\u001b[0m\n",
      "File \u001b[1;32mcupy\\_core\\_routines_indexing.pyx:263\u001b[0m, in \u001b[0;36mcupy._core._routines_indexing._prepare_slice_list\u001b[1;34m()\u001b[0m\n",
      "File \u001b[1;32m~\\anaconda3\\Lib\\site-packages\\numpy\\core\\numeric.py:1878\u001b[0m, in \u001b[0;36misscalar\u001b[1;34m(element)\u001b[0m\n\u001b[0;32m   1874\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m_frombuffer\u001b[39m(buf, dtype, shape, order):\n\u001b[0;32m   1875\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m frombuffer(buf, dtype\u001b[38;5;241m=\u001b[39mdtype)\u001b[38;5;241m.\u001b[39mreshape(shape, order\u001b[38;5;241m=\u001b[39morder)\n\u001b[1;32m-> 1878\u001b[0m \u001b[38;5;129m@set_module\u001b[39m(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mnumpy\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[0;32m   1879\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21misscalar\u001b[39m(element):\n\u001b[0;32m   1880\u001b[0m \u001b[38;5;250m    \u001b[39m\u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[0;32m   1881\u001b[0m \u001b[38;5;124;03m    Returns True if the type of `element` is a scalar type.\u001b[39;00m\n\u001b[0;32m   1882\u001b[0m \n\u001b[1;32m   (...)\u001b[0m\n\u001b[0;32m   1952\u001b[0m \n\u001b[0;32m   1953\u001b[0m \u001b[38;5;124;03m    \"\"\"\u001b[39;00m\n\u001b[0;32m   1954\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m (\u001b[38;5;28misinstance\u001b[39m(element, generic)\n\u001b[0;32m   1955\u001b[0m             \u001b[38;5;129;01mor\u001b[39;00m \u001b[38;5;28mtype\u001b[39m(element) \u001b[38;5;129;01min\u001b[39;00m ScalarType\n\u001b[0;32m   1956\u001b[0m             \u001b[38;5;129;01mor\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(element, numbers\u001b[38;5;241m.\u001b[39mNumber))\n",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "%%memit\n",
    "import cupy as cp\n",
    "\n",
    "# Grid parameters\n",
    "nr, ntheta, nz = 100, 100, 100  # Number of grid points in r, theta, z\n",
    "r_max = 10.0                    # Maximum radius\n",
    "z_max = 10.0                    # Maximum height\n",
    "\n",
    "dr = r_max / nr                 # Radial grid spacing\n",
    "dtheta = 2 * cp.pi / ntheta     # Angular grid spacing\n",
    "dz = z_max / nz                 # Axial grid spacing\n",
    "\n",
    "dt = 0.01                       # Time step\n",
    "alpha = 0.01                    # Thermal diffusivity\n",
    "nt = 100                        # Number of time steps\n",
    "\n",
    "# Create radial and angular coordinate arrays\n",
    "r = cp.linspace(0 + dr/2, r_max - dr/2, nr)  # Avoid r=0 to prevent division by zero\n",
    "theta = cp.linspace(0, 2 * cp.pi, ntheta, endpoint=False)\n",
    "z = cp.linspace(0, z_max, nz)\n",
    "\n",
    "def second_derivative_z(u_z):\n",
    "    \"\"\"\n",
    "    Compute the second derivative along the z-axis using central differences.\n",
    "    Assumes Neumann boundary conditions (zero gradient at boundaries).\n",
    "    \"\"\"\n",
    "    d2u_dz2 = cp.zeros_like(u_z)\n",
    "    d2u_dz2[1:-1] = (u_z[2:] - 2 * u_z[1:-1] + u_z[:-2]) / dz**2\n",
    "    # Neumann boundary conditions\n",
    "    d2u_dz2[0] = (u_z[1] - u_z[0]) / dz**2\n",
    "    d2u_dz2[-1] = (u_z[-2] - u_z[-1]) / dz**2\n",
    "    return d2u_dz2\n",
    "\n",
    "def main():\n",
    "    # Initialize the temperature field\n",
    "    u = cp.zeros((nr, ntheta, nz))\n",
    "    T0 = 100.0  # Initial temperature at z=0\n",
    "    u[:, :, 0] = T0  # Uniform temperature at z=0 across r and theta\n",
    "\n",
    "    r_squared = r[:, cp.newaxis] ** 2  # Shape: (nr, 1)\n",
    "\n",
    "    for _ in range(nt):\n",
    "        # Compute second derivative along the z-axis using apply_along_axis\n",
    "        d2u_dz2 = cp.apply_along_axis(second_derivative_z, axis=2, arr=u)\n",
    "\n",
    "        # Compute second derivative along the r-axis\n",
    "        d2u_dr2 = cp.zeros_like(u)\n",
    "        u_plus = cp.roll(u, -1, axis=0)\n",
    "        u_minus = cp.roll(u, 1, axis=0)\n",
    "        u_i = u\n",
    "\n",
    "        # Interior points\n",
    "        d2u_dr2[1:-1, :, :] = ((u_plus[1:-1] - 2 * u_i[1:-1] + u_minus[1:-1]) / dr**2 +\n",
    "                               (u_plus[1:-1] - u_minus[1:-1]) / (2 * dr * r[1:-1, cp.newaxis, cp.newaxis]))\n",
    "        # Boundary at r=0\n",
    "        d2u_dr2[0, :, :] = (2 * (u[1, :, :] - u[0, :, :]) / dr**2)\n",
    "        # Boundary at r=r_max\n",
    "        d2u_dr2[-1, :, :] = (2 * (u[-2, :, :] - u[-1, :, :]) / dr**2)\n",
    "\n",
    "        # Compute second derivative along the theta-axis\n",
    "        d2u_dtheta2 = cp.zeros_like(u)\n",
    "        u_jp = cp.roll(u, -1, axis=1)\n",
    "        u_jm = cp.roll(u, 1, axis=1)\n",
    "        numerator = (u_jp - 2 * u + u_jm) / dtheta**2  # Shape: (nr, ntheta, nz)\n",
    "        d2u_dtheta2 = numerator / r_squared[:, :, cp.newaxis]  # Broadcast over nz\n",
    "\n",
    "        # Update the temperature field\n",
    "        u += alpha * dt * (d2u_dr2 + d2u_dtheta2 + d2u_dz2)\n",
    "\n",
    "    # Optionally, return or save the final temperature field\n",
    "    return u\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    final_temperature = main()\n",
    "\n",
    "    # Transfer data back to CPU for further processing or visualization\n",
    "    final_temperature_cpu = cp.asnumpy(final_temperature)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "621496b8-17fc-48be-a179-cf76160dc0a4",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "Cell \u001b[1;32mIn[47], line 85\u001b[0m\n\u001b[0;32m     82\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m T\n\u001b[0;32m     84\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;18m__name__\u001b[39m \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m__main__\u001b[39m\u001b[38;5;124m\"\u001b[39m:\n\u001b[1;32m---> 85\u001b[0m     final_temperature \u001b[38;5;241m=\u001b[39m main()\n",
      "Cell \u001b[1;32mIn[47], line 80\u001b[0m, in \u001b[0;36mmain\u001b[1;34m()\u001b[0m\n\u001b[0;32m     78\u001b[0m \u001b[38;5;28;01mglobal\u001b[39;00m T\n\u001b[0;32m     79\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m n \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mrange\u001b[39m(nt):\n\u001b[1;32m---> 80\u001b[0m     T \u001b[38;5;241m=\u001b[39m heat_diffusion(T, alpha, dr, dtheta, dphi, dt, r, theta)\n\u001b[0;32m     81\u001b[0m     \u001b[38;5;66;03m# Optional: Add code to monitor progress or visualize intermediate results\u001b[39;00m\n\u001b[0;32m     82\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m T\n",
      "Cell \u001b[1;32mIn[47], line 64\u001b[0m, in \u001b[0;36mheat_diffusion\u001b[1;34m(T, alpha, dr, dtheta, dphi, dt, r, theta)\u001b[0m\n\u001b[0;32m     61\u001b[0m             term_theta \u001b[38;5;241m=\u001b[39m (\u001b[38;5;241m1\u001b[39m \u001b[38;5;241m/\u001b[39m (r[i]\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m2\u001b[39m \u001b[38;5;241m*\u001b[39m sin_theta[j])) \u001b[38;5;241m*\u001b[39m (sin_theta[j] \u001b[38;5;241m*\u001b[39m d2T_dtheta2 \u001b[38;5;241m+\u001b[39m cos_theta[j] \u001b[38;5;241m*\u001b[39m dT_dtheta)\n\u001b[0;32m     62\u001b[0m             term_phi \u001b[38;5;241m=\u001b[39m d2T_dphi2 \u001b[38;5;241m/\u001b[39m (r[i]\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m2\u001b[39m \u001b[38;5;241m*\u001b[39m sin_theta[j]\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m2\u001b[39m)\n\u001b[1;32m---> 64\u001b[0m             T_new[i, j, k] \u001b[38;5;241m=\u001b[39m T[i, j, k] \u001b[38;5;241m+\u001b[39m alpha \u001b[38;5;241m*\u001b[39m dt \u001b[38;5;241m*\u001b[39m (term_r \u001b[38;5;241m+\u001b[39m term_theta \u001b[38;5;241m+\u001b[39m term_phi)\n\u001b[0;32m     66\u001b[0m \u001b[38;5;66;03m# Boundary conditions\u001b[39;00m\n\u001b[0;32m     67\u001b[0m \u001b[38;5;66;03m# At r = 0 (center), assume symmetry (no flux across center)\u001b[39;00m\n\u001b[0;32m     68\u001b[0m T_new[\u001b[38;5;241m0\u001b[39m, :, :] \u001b[38;5;241m=\u001b[39m T_new[\u001b[38;5;241m1\u001b[39m, :, :]\n",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "# Physical constants and parameters\n",
    "R_moon = 1737.4e3  # Radius of the Moon in meters\n",
    "alpha = 1e-6       # Thermal diffusivity in m^2/s (typical for rock)\n",
    "T0 = 300.0         # Maximum initial temperature in Kelvin\n",
    "\n",
    "# Grid parameters\n",
    "nr = 50           # Number of radial grid points\n",
    "ntheta = 50       # Number of polar angle grid points (θ from 0 to π/2)\n",
    "nphi = 100        # Number of azimuthal angle grid points (φ from 0 to 2π)\n",
    "nt = 1000         # Number of time steps\n",
    "dt = 10.0         # Time step in seconds\n",
    "\n",
    "# Spatial discretization\n",
    "r = np.linspace(0, R_moon, nr)\n",
    "theta = np.linspace(0, np.pi / 2, ntheta)\n",
    "phi = np.linspace(0, 2 * np.pi, nphi, endpoint=False)\n",
    "\n",
    "dr = r[1] - r[0]\n",
    "dtheta = theta[1] - theta[0]\n",
    "dphi = phi[1] - phi[0]\n",
    "\n",
    "# Create 3D mesh grids\n",
    "R, Theta, Phi = np.meshgrid(r, theta, phi, indexing='ij')\n",
    "\n",
    "# Initial temperature distribution based on incident angle\n",
    "T = np.zeros((nr, ntheta, nphi))\n",
    "T[:, :, :] = T0 * np.cos(Theta)\n",
    "\n",
    "def heat_diffusion(T, alpha, dr, dtheta, dphi, dt, r, theta):\n",
    "    \"\"\"\n",
    "    Perform one time step of the heat diffusion equation in spherical coordinates.\n",
    "    \"\"\"\n",
    "    T_new = np.copy(T)\n",
    "    \n",
    "    # Precompute sines and cosines\n",
    "    sin_theta = np.sin(theta)\n",
    "    sin_theta[sin_theta == 0] = 1e-10  # Avoid division by zero\n",
    "    cos_theta = np.cos(theta)\n",
    "    \n",
    "    # Compute derivatives\n",
    "    for i in range(1, nr - 1):\n",
    "        for j in range(1, ntheta - 1):\n",
    "            for k in range(nphi):\n",
    "                # Radial derivatives\n",
    "                dT_dr = (T[i + 1, j, k] - T[i - 1, j, k]) / (2 * dr)\n",
    "                d2T_dr2 = (T[i + 1, j, k] - 2 * T[i, j, k] + T[i - 1, j, k]) / dr**2\n",
    "                \n",
    "                # Theta derivatives\n",
    "                dT_dtheta = (T[i, j + 1, k] - T[i, j - 1, k]) / (2 * dtheta)\n",
    "                d2T_dtheta2 = (T[i, j + 1, k] - 2 * T[i, j, k] + T[i, j - 1, k]) / dtheta**2\n",
    "                \n",
    "                # Phi derivatives (periodic boundary conditions)\n",
    "                k_plus = (k + 1) % nphi\n",
    "                k_minus = (k - 1) % nphi\n",
    "                d2T_dphi2 = (T[i, j, k_plus] - 2 * T[i, j, k] + T[i, j, k_minus]) / dphi**2\n",
    "                \n",
    "                # Heat equation in spherical coordinates\n",
    "                term_r = (1 / r[i]**2) * (r[i]**2 * d2T_dr2 + 2 * r[i] * dT_dr)\n",
    "                term_theta = (1 / (r[i]**2 * sin_theta[j])) * (sin_theta[j] * d2T_dtheta2 + cos_theta[j] * dT_dtheta)\n",
    "                term_phi = d2T_dphi2 / (r[i]**2 * sin_theta[j]**2)\n",
    "                \n",
    "                T_new[i, j, k] = T[i, j, k] + alpha * dt * (term_r + term_theta + term_phi)\n",
    "    \n",
    "    # Boundary conditions\n",
    "    # At r = 0 (center), assume symmetry (no flux across center)\n",
    "    T_new[0, :, :] = T_new[1, :, :]\n",
    "    # At r = R_moon, insulated boundary (Neumann BC)\n",
    "    T_new[-1, :, :] = T_new[-2, :, :]\n",
    "    # At θ = 0 and θ = π/2, symmetry boundary conditions\n",
    "    T_new[:, 0, :] = T_new[:, 1, :]\n",
    "    T_new[:, -1, :] = T_new[:, -2, :]\n",
    "    \n",
    "    return T_new\n",
    "\n",
    "def main():\n",
    "    global T\n",
    "    for n in range(nt):\n",
    "        T = heat_diffusion(T, alpha, dr, dtheta, dphi, dt, r, theta)\n",
    "        # Optional: Add code to monitor progress or visualize intermediate results\n",
    "    return T\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    final_temperature = main()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1293f4b-54a9-4e88-b3ac-02e6c5b0c8e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Extract a radial slice at θ = π/4 and φ = 0\n",
    "theta_index = ntheta // 2\n",
    "phi_index = 0\n",
    "temperature_slice = final_temperature[:, theta_index, phi_index]\n",
    "\n",
    "plt.plot(r, temperature_slice)\n",
    "plt.xlabel('Radius (m)')\n",
    "plt.ylabel('Temperature (K)')\n",
    "plt.title('Temperature Profile along Radius at θ = π/4, φ = 0')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "3b057707-c6a2-4100-b718-af4045776d9d",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "ename": "MemoryError",
     "evalue": "Unable to allocate 18.3 GiB for an array with shape (50, 50, 100, 100, 98) and data type float64",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mMemoryError\u001b[0m                               Traceback (most recent call last)",
      "Cell \u001b[1;32mIn[49], line 57\u001b[0m\n\u001b[0;32m     54\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m T\n\u001b[0;32m     56\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;18m__name__\u001b[39m \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m__main__\u001b[39m\u001b[38;5;124m\"\u001b[39m:\n\u001b[1;32m---> 57\u001b[0m     final_temperature \u001b[38;5;241m=\u001b[39m main()\n",
      "Cell \u001b[1;32mIn[49], line 52\u001b[0m, in \u001b[0;36mmain\u001b[1;34m()\u001b[0m\n\u001b[0;32m     50\u001b[0m \u001b[38;5;28;01mglobal\u001b[39;00m T\n\u001b[0;32m     51\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m n \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mrange\u001b[39m(nt):\n\u001b[1;32m---> 52\u001b[0m     T \u001b[38;5;241m=\u001b[39m heat_diffusion_depth(T, alpha, dz, dt)\n\u001b[0;32m     53\u001b[0m     \u001b[38;5;66;03m# Optional: Monitor progress or visualize intermediate results\u001b[39;00m\n\u001b[0;32m     54\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m T\n",
      "Cell \u001b[1;32mIn[49], line 36\u001b[0m, in \u001b[0;36mheat_diffusion_depth\u001b[1;34m(T, alpha, dz, dt)\u001b[0m\n\u001b[0;32m     34\u001b[0m \u001b[38;5;66;03m# Compute second derivative along the depth axis using central differences\u001b[39;00m\n\u001b[0;32m     35\u001b[0m d2T_dz2 \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mzeros_like(T)\n\u001b[1;32m---> 36\u001b[0m d2T_dz2[\u001b[38;5;241m.\u001b[39m\u001b[38;5;241m.\u001b[39m\u001b[38;5;241m.\u001b[39m, \u001b[38;5;241m1\u001b[39m:\u001b[38;5;241m-\u001b[39m\u001b[38;5;241m1\u001b[39m] \u001b[38;5;241m=\u001b[39m (T[\u001b[38;5;241m.\u001b[39m\u001b[38;5;241m.\u001b[39m\u001b[38;5;241m.\u001b[39m, \u001b[38;5;241m2\u001b[39m:] \u001b[38;5;241m-\u001b[39m \u001b[38;5;241m2\u001b[39m \u001b[38;5;241m*\u001b[39m T[\u001b[38;5;241m.\u001b[39m\u001b[38;5;241m.\u001b[39m\u001b[38;5;241m.\u001b[39m, \u001b[38;5;241m1\u001b[39m:\u001b[38;5;241m-\u001b[39m\u001b[38;5;241m1\u001b[39m] \u001b[38;5;241m+\u001b[39m T[\u001b[38;5;241m.\u001b[39m\u001b[38;5;241m.\u001b[39m\u001b[38;5;241m.\u001b[39m, :\u001b[38;5;241m-\u001b[39m\u001b[38;5;241m2\u001b[39m]) \u001b[38;5;241m/\u001b[39m dz\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m2\u001b[39m\n\u001b[0;32m     38\u001b[0m \u001b[38;5;66;03m# Neumann boundary conditions at the surface (z=0) and bottom (z=depth)\u001b[39;00m\n\u001b[0;32m     39\u001b[0m \u001b[38;5;66;03m# Surface boundary (z=0): assume constant temperature (due to incident radiation)\u001b[39;00m\n\u001b[0;32m     40\u001b[0m d2T_dz2[\u001b[38;5;241m.\u001b[39m\u001b[38;5;241m.\u001b[39m\u001b[38;5;241m.\u001b[39m, \u001b[38;5;241m0\u001b[39m] \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m0\u001b[39m  \u001b[38;5;66;03m# No change at the surface\u001b[39;00m\n",
      "\u001b[1;31mMemoryError\u001b[0m: Unable to allocate 18.3 GiB for an array with shape (50, 50, 100, 100, 98) and data type float64"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "# Physical constants and parameters\n",
    "R_moon = 1737.4e3   # Radius of the Moon in meters\n",
    "alpha = 1e-6        # Thermal diffusivity in m^2/s (typical for lunar rock)\n",
    "T0 = 300.0          # Maximum initial temperature in Kelvin\n",
    "\n",
    "# Grid parameters\n",
    "nr = 50             # Number of radial grid points\n",
    "ntheta = 50         # Number of polar angle grid points (θ from 0 to π/2)\n",
    "nphi = 100          # Number of azimuthal angle grid points (φ from 0 to 2π)\n",
    "nz = 100            # Number of depth grid points\n",
    "nt = 1000           # Number of time steps\n",
    "dt = 10.0           # Time step in seconds\n",
    "\n",
    "# Spatial discretization\n",
    "r = np.linspace(0, R_moon, nr)\n",
    "theta = np.linspace(0, np.pi / 2, ntheta)\n",
    "phi = np.linspace(0, 2 * np.pi, nphi, endpoint=False)\n",
    "z = np.linspace(0, 1000, nz)  # Depth from the surface into the moon (in meters)\n",
    "\n",
    "dz = z[1] - z[0]\n",
    "\n",
    "# Create 4D mesh grids\n",
    "R, Theta, Phi, Z = np.meshgrid(r, theta, phi, z, indexing='ij')\n",
    "\n",
    "# Initial temperature distribution based on incident angle, uniform along depth\n",
    "T = T0 * np.cos(Theta)[..., np.newaxis] * np.ones(nz)\n",
    "\n",
    "def heat_diffusion_depth(T, alpha, dz, dt):\n",
    "    \"\"\"\n",
    "    Perform one time step of the heat diffusion equation along the depth axis.\n",
    "    \"\"\"\n",
    "    # Compute second derivative along the depth axis using central differences\n",
    "    d2T_dz2 = np.zeros_like(T)\n",
    "    d2T_dz2[..., 1:-1] = (T[..., 2:] - 2 * T[..., 1:-1] + T[..., :-2]) / dz**2\n",
    "\n",
    "    # Neumann boundary conditions at the surface (z=0) and bottom (z=depth)\n",
    "    # Surface boundary (z=0): assume constant temperature (due to incident radiation)\n",
    "    d2T_dz2[..., 0] = 0  # No change at the surface\n",
    "    # Bottom boundary (z=max): insulated boundary (zero gradient)\n",
    "    d2T_dz2[..., -1] = (T[..., -2] - T[..., -1]) / dz**2\n",
    "\n",
    "    # Update temperature field\n",
    "    T_new = T + alpha * dt * d2T_dz2\n",
    "\n",
    "    return T_new\n",
    "\n",
    "def main():\n",
    "    global T\n",
    "    for n in range(nt):\n",
    "        T = heat_diffusion_depth(T, alpha, dz, dt)\n",
    "        # Optional: Monitor progress or visualize intermediate results\n",
    "    return T\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    final_temperature = main()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21583c78-4e51-494a-abb4-ce89c28c66f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Example: Temperature profile along depth at a specific location\n",
    "r_index = nr // 2      # Midpoint in radial direction\n",
    "theta_index = ntheta // 4  # Some angle\n",
    "phi_index = nphi // 2  # Midpoint in azimuthal direction\n",
    "\n",
    "depth_profile = final_temperature[r_index, theta_index, phi_index, :]\n",
    "\n",
    "plt.figure(figsize=(8, 6))\n",
    "plt.plot(z, depth_profile)\n",
    "plt.xlabel('Depth (m)')\n",
    "plt.ylabel('Temperature (K)')\n",
    "plt.title('Temperature Profile Along Depth')\n",
    "plt.gca().invert_xaxis()  # Depth increases into the moon\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "773c62fc-19fc-4a65-bbd7-251ca8c389aa",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\seane\\anaconda3\\Lib\\site-packages\\cupy\\_environment.py:217: UserWarning: CUDA path could not be detected. Set CUDA_PATH environment variable if CuPy fails to load.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "ename": "RuntimeError",
     "evalue": "CuPy failed to load nvrtc64_120_0.dll: FileNotFoundError: Could not find module 'nvrtc64_120_0.dll' (or one of its dependencies). Try using the full path with constructor syntax.",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mFileNotFoundError\u001b[0m                         Traceback (most recent call last)",
      "File \u001b[1;32mcupy_backends\\cuda\\_softlink.pyx:25\u001b[0m, in \u001b[0;36mcupy_backends.cuda._softlink.SoftLink.__init__\u001b[1;34m()\u001b[0m\n",
      "File \u001b[1;32m~\\anaconda3\\Lib\\ctypes\\__init__.py:376\u001b[0m, in \u001b[0;36mCDLL.__init__\u001b[1;34m(self, name, mode, handle, use_errno, use_last_error, winmode)\u001b[0m\n\u001b[0;32m    375\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m handle \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m--> 376\u001b[0m     \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_handle \u001b[38;5;241m=\u001b[39m _dlopen(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_name, mode)\n\u001b[0;32m    377\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n",
      "\u001b[1;31mFileNotFoundError\u001b[0m: Could not find module 'nvrtc64_120_0.dll' (or one of its dependencies). Try using the full path with constructor syntax.",
      "\nThe above exception was the direct cause of the following exception:\n",
      "\u001b[1;31mRuntimeError\u001b[0m                              Traceback (most recent call last)",
      "Cell \u001b[1;32mIn[48], line 17\u001b[0m\n\u001b[0;32m     14\u001b[0m dt \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m10.0\u001b[39m           \u001b[38;5;66;03m# Time step in seconds\u001b[39;00m\n\u001b[0;32m     16\u001b[0m \u001b[38;5;66;03m# Spatial discretization\u001b[39;00m\n\u001b[1;32m---> 17\u001b[0m r \u001b[38;5;241m=\u001b[39m cp\u001b[38;5;241m.\u001b[39mlinspace(\u001b[38;5;241m0\u001b[39m, R_moon, nr)\n\u001b[0;32m     18\u001b[0m theta \u001b[38;5;241m=\u001b[39m cp\u001b[38;5;241m.\u001b[39mlinspace(\u001b[38;5;241m0\u001b[39m, cp\u001b[38;5;241m.\u001b[39mpi \u001b[38;5;241m/\u001b[39m \u001b[38;5;241m2\u001b[39m, ntheta)\n\u001b[0;32m     19\u001b[0m phi \u001b[38;5;241m=\u001b[39m cp\u001b[38;5;241m.\u001b[39mlinspace(\u001b[38;5;241m0\u001b[39m, \u001b[38;5;241m2\u001b[39m \u001b[38;5;241m*\u001b[39m cp\u001b[38;5;241m.\u001b[39mpi, nphi, endpoint\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mFalse\u001b[39;00m)\n",
      "File \u001b[1;32m~\\anaconda3\\Lib\\site-packages\\cupy\\_creation\\ranges.py:161\u001b[0m, in \u001b[0;36mlinspace\u001b[1;34m(start, stop, num, endpoint, retstep, dtype, axis)\u001b[0m\n\u001b[0;32m    159\u001b[0m scalar_stop \u001b[38;5;241m=\u001b[39m cupy\u001b[38;5;241m.\u001b[39misscalar(stop)\n\u001b[0;32m    160\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m scalar_start \u001b[38;5;129;01mand\u001b[39;00m scalar_stop:\n\u001b[1;32m--> 161\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m _linspace_scalar(start, stop, num, endpoint, retstep, dtype)\n\u001b[0;32m    163\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m scalar_start:\n\u001b[0;32m    164\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m (\u001b[38;5;28misinstance\u001b[39m(start, cupy\u001b[38;5;241m.\u001b[39mndarray) \u001b[38;5;129;01mand\u001b[39;00m start\u001b[38;5;241m.\u001b[39mdtype\u001b[38;5;241m.\u001b[39mkind \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m'\u001b[39m):\n",
      "File \u001b[1;32m~\\anaconda3\\Lib\\site-packages\\cupy\\_creation\\ranges.py:105\u001b[0m, in \u001b[0;36m_linspace_scalar\u001b[1;34m(start, stop, num, endpoint, retstep, dtype)\u001b[0m\n\u001b[0;32m    103\u001b[0m     _linspace_ufunc_underflow(start, stop \u001b[38;5;241m-\u001b[39m start, div, ret)\n\u001b[0;32m    104\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m--> 105\u001b[0m     _linspace_ufunc(start, step, ret)\n\u001b[0;32m    107\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m endpoint:\n\u001b[0;32m    108\u001b[0m     \u001b[38;5;66;03m# Here num == div + 1 > 1 is ensured.\u001b[39;00m\n\u001b[0;32m    109\u001b[0m     ret[\u001b[38;5;241m-\u001b[39m\u001b[38;5;241m1\u001b[39m] \u001b[38;5;241m=\u001b[39m stop\n",
      "File \u001b[1;32mcupy\\_core\\_kernel.pyx:1375\u001b[0m, in \u001b[0;36mcupy._core._kernel.ufunc.__call__\u001b[1;34m()\u001b[0m\n",
      "File \u001b[1;32mcupy\\_core\\_kernel.pyx:1402\u001b[0m, in \u001b[0;36mcupy._core._kernel.ufunc._get_ufunc_kernel\u001b[1;34m()\u001b[0m\n",
      "File \u001b[1;32mcupy\\_core\\_kernel.pyx:1082\u001b[0m, in \u001b[0;36mcupy._core._kernel._get_ufunc_kernel\u001b[1;34m()\u001b[0m\n",
      "File \u001b[1;32mcupy\\_core\\_kernel.pyx:94\u001b[0m, in \u001b[0;36mcupy._core._kernel._get_simple_elementwise_kernel\u001b[1;34m()\u001b[0m\n",
      "File \u001b[1;32mcupy\\_core\\_kernel.pyx:82\u001b[0m, in \u001b[0;36mcupy._core._kernel._get_simple_elementwise_kernel_from_code\u001b[1;34m()\u001b[0m\n",
      "File \u001b[1;32mcupy\\_core\\core.pyx:2281\u001b[0m, in \u001b[0;36mcupy._core.core.compile_with_cache\u001b[1;34m()\u001b[0m\n",
      "File \u001b[1;32mcupy\\_core\\core.pyx:2219\u001b[0m, in \u001b[0;36mcupy._core.core.assemble_cupy_compiler_options\u001b[1;34m()\u001b[0m\n",
      "File \u001b[1;32mcupy_backends\\cuda\\libs\\nvrtc.pyx:57\u001b[0m, in \u001b[0;36mcupy_backends.cuda.libs.nvrtc.getVersion\u001b[1;34m()\u001b[0m\n",
      "File \u001b[1;32mcupy_backends\\cuda\\libs\\_cnvrtc.pxi:72\u001b[0m, in \u001b[0;36mcupy_backends.cuda.libs.nvrtc.initialize\u001b[1;34m()\u001b[0m\n",
      "File \u001b[1;32mcupy_backends\\cuda\\libs\\_cnvrtc.pxi:76\u001b[0m, in \u001b[0;36mcupy_backends.cuda.libs.nvrtc._initialize\u001b[1;34m()\u001b[0m\n",
      "File \u001b[1;32mcupy_backends\\cuda\\libs\\_cnvrtc.pxi:143\u001b[0m, in \u001b[0;36mcupy_backends.cuda.libs.nvrtc._get_softlink\u001b[1;34m()\u001b[0m\n",
      "File \u001b[1;32mcupy_backends\\cuda\\_softlink.pyx:32\u001b[0m, in \u001b[0;36mcupy_backends.cuda._softlink.SoftLink.__init__\u001b[1;34m()\u001b[0m\n",
      "\u001b[1;31mRuntimeError\u001b[0m: CuPy failed to load nvrtc64_120_0.dll: FileNotFoundError: Could not find module 'nvrtc64_120_0.dll' (or one of its dependencies). Try using the full path with constructor syntax."
     ]
    }
   ],
   "source": [
    "import cupy as cp\n",
    "\n",
    "# Physical constants and parameters\n",
    "R_moon = 1737.4e3   # Radius of the Moon in meters\n",
    "alpha = 1e-6        # Thermal diffusivity in m^2/s (typical for lunar rock)\n",
    "T0 = 300.0          # Maximum initial temperature in Kelvin\n",
    "\n",
    "# Grid parameters\n",
    "nr = 50             # Number of radial grid points\n",
    "ntheta = 50         # Number of polar angle grid points\n",
    "nphi = 100          # Number of azimuthal angle grid points\n",
    "nz = 100            # Number of depth grid points\n",
    "nt = 1000           # Number of time steps\n",
    "dt = 10.0           # Time step in seconds\n",
    "\n",
    "# Spatial discretization\n",
    "r = cp.linspace(0, R_moon, nr)\n",
    "theta = cp.linspace(0, cp.pi / 2, ntheta)\n",
    "phi = cp.linspace(0, 2 * cp.pi, nphi, endpoint=False)\n",
    "z = cp.linspace(0, 1000, nz)  # Depth from the surface into the moon (in meters)\n",
    "\n",
    "dz = z[1] - z[0]\n",
    "\n",
    "# Create 4D mesh grids\n",
    "R, Theta, Phi, Z = cp.meshgrid(r, theta, phi, z, indexing='ij')\n",
    "\n",
    "# Initial temperature distribution based on incident angle, uniform along depth\n",
    "T = T0 * cp.cos(Theta)[..., cp.newaxis] * cp.ones(nz)\n",
    "\n",
    "def heat_diffusion_depth(T, alpha, dz, dt):\n",
    "    \"\"\"\n",
    "    Perform one time step of the heat diffusion equation along the depth axis.\n",
    "    \"\"\"\n",
    "    # Compute second derivative along the depth axis using central differences\n",
    "    d2T_dz2 = cp.zeros_like(T)\n",
    "    d2T_dz2[..., 1:-1] = (T[..., 2:] - 2 * T[..., 1:-1] + T[..., :-2]) / dz**2\n",
    "\n",
    "    # Neumann boundary conditions at the surface (z=0) and bottom (z=depth)\n",
    "    # Surface boundary (z=0): assume constant temperature (due to incident radiation)\n",
    "    d2T_dz2[..., 0] = 0  # No change at the surface\n",
    "    # Bottom boundary (z=max): insulated boundary (zero gradient)\n",
    "    d2T_dz2[..., -1] = (T[..., -2] - T[..., -1]) / dz**2\n",
    "\n",
    "    # Update temperature field\n",
    "    T_new = T + alpha * dt * d2T_dz2\n",
    "\n",
    "    return T_new\n",
    "\n",
    "def main():\n",
    "    global T\n",
    "    for n in range(nt):\n",
    "        T = heat_diffusion_depth(T, alpha, dz, dt)\n",
    "        # Optional: Monitor progress or visualize intermediate results\n",
    "    return T\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    final_temperature = main()\n",
    "\n",
    "    # Transfer data back to CPU for further processing or visualization\n",
    "    final_temperature_cpu = cp.asnumpy(final_temperature)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a727d04-db05-4123-8796-adc16611cac4",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}