{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c37ac15a-66a3-452b-8d2f-1e02daafe7cc",
   "metadata": {},
   "source": [
    "# HSE Convergence"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4ee377c-72a0-4afb-98bc-87a01a8fd382",
   "metadata": {},
   "source": [
    "Here we explore the convergence of the HSE well-balanced method.\n",
    "\n",
    "We use reflecting boundary conditions and the `HSEPPMInterpolant` reconstruction of the pressure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e9274186-dea7-4254-a5e1-1c10e4ede36c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ppmpy.euler import Euler\n",
    "from ppmpy.gravity import constant_gravity\n",
    "from ppmpy.initial_conditions import hse\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf604c93-4360-4431-a3f0-64a098aa8d4a",
   "metadata": {},
   "source": [
    "## Convergence testing with HSE reconstruction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "80932eef-a2f0-4c5b-b267-06255b543d56",
   "metadata": {},
   "outputs": [],
   "source": [
    "params = {\"base_density\": 1.0, \"base_pressure\": 1.0, \"g_const\": -1.0}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "898eb0ed-1bec-41f0-9431-45b2ad4c5787",
   "metadata": {},
   "outputs": [],
   "source": [
    "simulations = []\n",
    "for nx in [32, 64, 128, 256, 512]:\n",
    "    dt = 0.015625 * (32 / nx)\n",
    "    e = Euler(nx, 0.5, fixed_dt=dt, init_cond=hse, grav_func=constant_gravity,\n",
    "              use_limiting=True, use_flattening=True,\n",
    "              use_hse_reconstruction=True,\n",
    "              bc_left_type=\"reflect\", bc_right_type=\"reflect\", params=params)\n",
    "    e.evolve(0.5, verbose=False)\n",
    "    simulations.append(e)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2aa51975-6eff-4d44-8bd6-765fa150ceff",
   "metadata": {},
   "source": [
    "Richardson convergence testing--compare adjacent resolution runs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "313edca0-aa4c-4dac-ab14-7a6cff0455bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "from itertools import pairwise\n",
    "ivar = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "6935138d-852d-4d56-ac43-4960cd4f4b6b",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 64 ->  32 : 3.6537626526206627e-05\n",
      "128 ->  64 : 9.208712352368843e-06\n",
      "256 -> 128 : 2.3117807006939062e-06\n",
      "512 -> 256 : 5.794725187474326e-07\n"
     ]
    }
   ],
   "source": [
    "for coarse, fine in pairwise(simulations):\n",
    "    _, cd = fine.grid.coarsen(fine.U[:, ivar])\n",
    "    err = coarse.grid.norm(coarse.U[:, ivar] - cd)\n",
    "    print(f\"{fine.grid.nx:3d} -> {coarse.grid.nx:3d} : {err}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33af9113-8f88-4570-9d84-e6ff78229ba3",
   "metadata": {},
   "source": [
    "Compare to the initial conditions.  If we were in HSE, then the density should be the same as it was originally."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "66950c33-44e4-4975-9e96-c569396b8f5c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 32 : 1.5642584519965996e-05\n",
      " 64 : 3.955429302849995e-06\n",
      "128 : 9.94413687738086e-07\n",
      "256 : 2.5011318512955466e-07\n",
      "512 : 6.259701007232347e-08\n"
     ]
    }
   ],
   "source": [
    "for s in simulations:\n",
    "    print(f\"{s.grid.nx:3d} : {s.grid.norm(s.U_init[:, ivar] - s.U[:, ivar]) }\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4ae9e92-a5b2-41cf-b38b-314b83b676b3",
   "metadata": {},
   "source": [
    "Visualize the atmosphere at the start and end."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "d4e9df43-d976-4e00-b99b-213f52b60839",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fbde5466cc0>]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABLsklEQVR4nO3dd3hUZfr/8feZkoQACUsLoYcYkLJSEvqiYEHRRRELioIFXbEhoK6wdnTFBmLQYEF0UVSkWX5iyXcVBSkhFFsQYoiGEggJkISSqef3B2Q0S4KZSDKZyed1XXNd5MxzMvecZTm3z3me+zZM0zQRERERCRBLoAMQERGRuk3JiIiIiASUkhEREREJKCUjIiIiElBKRkRERCSglIyIiIhIQCkZERERkYBSMiIiIiIBZQt0AJXh9XrZvXs3DRs2xDCMQIcjIiIilWCaJsXFxbRs2RKLpeL5j6BIRnbv3k2bNm0CHYaIiIhUwY4dO2jdunWF7wdFMtKwYUM4/mWioqICHY6IiIhUQlFREW3atPHdxysSFMlI6aOZqKgoJSMiIiJB5o+WWGgBq4iIiASU38nI119/zfDhw2nZsiWGYfD+++//4TlfffUViYmJRERE0KFDB1566aWqxisiIiIhxu9k5PDhw3Tv3p0XXnihUuOzs7O58MILGTRoEJs2beJf//oXEyZMYMmSJVWJV0REREKM32tGhg0bxrBhwyo9/qWXXqJt27bMmjULgM6dO5Oens6zzz7LZZdd5u/Hi4iISIip9jUja9asYejQoWWOnX/++aSnp+Nyuco9x+FwUFRUVOYlIiIioanak5E9e/YQExNT5lhMTAxut5v8/Pxyz5k+fTrR0dG+l2qMiIiIhK4a2U3zv1t6TNMs93ipqVOnUlhY6Hvt2LGjJsIUERGRAKj2OiMtWrRgz549ZY7l5eVhs9lo0qRJueeEh4cTHh5e3aGJiIhILVDtyUj//v356KOPyhz7/PPPSUpKwm63V/fHV8jjNUnL3k9ecQnNG0bQJ64xVov63oiIiNQ0v5ORQ4cO8fPPP/t+zs7OZvPmzTRu3Ji2bdsydepUdu3axfz58wEYP348L7zwApMnT+bmm29mzZo1vPbaa7zzzjun9ptU1pfTydx3hLFZg8ktLPEdjo2OYH78ChKaRcKQqYGJTUREpA7ye81Ieno6PXv2pGfPngBMnjyZnj178tBDDwGQm5tLTk6Ob3xcXBzLly9nxYoV9OjRg8cee4zk5OSAbevN3HeEhIxkLj/0dpnjVxx6m4SMZDL3HQlIXCIiInWV3zMjgwcP9i1ALc8bb7xxwrGzzjqLjRs3+h/dKebxmozNGszlrt3cbV8MwGzPSO60LmWyfTEzXZezKGswq7ymHtmIiIjUkKBolHeqpGXvJ7ewhNmMBOBu+2Im2ZZgMUxmuC5ntmckFJaQlr2f/vHlL64VERGRU6tONcrLK/5tjcgLnhF4TbAYJl7T4BXP38sdJyIiItWrTiUjzRtG+P58h/V9LAaYxxOSD8MewMB7wjgRERGpXnUqGekT15jY6AgmWJdyt30xM1yXM9p1P27ToJNlJwvtjxEbfWybr4iIiNSMOpWMWC0G8+NX+BarzvaMZI23K/e6xgPQx7qV18Ke1eJVERGRGlSnkhGAhGaRZHaZwKIGo33HlnkH8YJ5JQCnF61m8/+9G8AIRURE6pY6tZsGgCFTSQBW/U8F1t7tLmD9C4fpffBjOq6cQGbTViT0GBToaEVEREKeYZ6saEgtUVRURHR0NIWFhURFRVXb57icDn6aeQF/LdlIPo1wXv85Ldt3qrbPExERCWWVvX/Xucc0J2MPC6f9rUvYbm1PUw7imj+Swv37Ah2WiIhISFMy8j8aRjemwQ3LyKMx7bw72fHSSBwlKhEvIiJSXZSMlKN56w4cvvwdDpn16Ob8ju9SxuL1eAMdloiISEhSMlKBuG79+OWcObhNC72LUlk7b3KgQxIREQlJSkZOotuZl7Kp+yMADNj1OusWPxfokEREREKOkpE/0HvkXaxtfSMAid9P49svl+DxmqzJKuCDzbtYk1WAx1vrNySJiIjUWnWvzkgV9L1xBunP7yKp8DNOX/EPrl21nzWHW/rej42OYH78ChKaRcKQqQGNVUREJNhoZqQSDIuFM26bT64lhnDDzSvuB2lBge/9Kw69TUJGMpn7tOtGRETEX0pGKslqD+da6zPkexvS0DjKR2H305Aj3Gld6ut1MzZrsB7ZiIiI+EnJSCWlZe8nq9jGJc7HOWRG0MxSxLfhN/u6/yZ7RpJbWEJa9v5AhyoiIhJUlIxUUl5xCQC7aMZVzgcwTbAYJh7T4AXPiBPGiYiISOUoGamk5g0jfH8eYtmMYYBpgtUwec8+rdxxIiIi8seUjFRSn7jGxEZHMMG61PdoZrLrVgB6W7fxhv1JYqMj6BPXONChioiIBBVt7a0kq8U4tn0349hi1dmekQA0cx3kX/Z3GGz9jrkRyVgt5wQ6VBERkaCimRE/JDSLJLPLBBY1GO079orn7/zHvAiA0w9+xeYVSwIYoYiISPAxTNOs9XtRi4qKiI6OprCwkKioqECHg8drkpa9n7ziEpo3jCCpbTTfzb6SxKL/ctgMZ8clizi911mBDlNERCSgKnv/VjJyiricJWydOYxuJRvZTxSHrvmYtglnBDosERGRgKns/VuPaU4Re1gEcbctIcsaT2OKsL59OftycwIdloiISK2nZOQUqh/VmEb/+IBdRgytzL0Uzh1BUaGKoImIiJyMkpFTrElMG4xrl7GfKE7zZPHri5fiKFHPGhERkYooGakGLeO7sn/EAo6Y4fzVuZnvX7gGj8cT6LBERERqJSUj1eS0HmeSfc7LuEwrSYe+IO2l8bjdHtZkFfDB5l2sySpQUz0REREVPateXc+8lA1FeSSm/5P++95jxmNhzHZc6Hs/NjriWCG1ZpEwZGpAYxUREQkUzYxUs8S/38JP9fsAcLfxFpdaVvreu+LQ2yRkJJO5T2tKRESk7lIyUs08XpMb3FNY7+kIwLP2lzjL8i13Wpcy2X6stPzYrMF6ZCMiInWWkpFqlpa9n9zCEq50PUSGty1Ww+QN+1O+ZnvJnpHkFpaQlq0twCIiUjcpGalmecUlAJhYuMT5OF7TwDDANOETb58TxomIiNQ1SkaqWfOGEb4/j7d+iMUw8ZpgGLAs7CFaG/tOGCciIlKXKBmpZn3iGhMbHcEE61Lfo5lejpfJ90bR0CjhI/u/6NzwCH3iGgc6VBERkYBQMlLNrBaD+fErfItVZ3tGcpCGXOicTqE3kr9YDvOW624OHcwPdKgiIiIBoWSkBiQ0iySzywQWNRjtO5bHX7gh7GkOE04TCtk9ZzhHDhUGNE4REZFAMEzTrPV7Sivbgri283hN0rL3k1dcQvOGEfSJa8yOLWk0XnQpURzm24gkTp/0/wgPrxfoUEVERP60yt6/NTNSg6wWg/7xTbikRyv6xzfBajFo37UvuRfN56gZRveSdL6ffTUetzvQoYqIiNQYJSO1QKfe55J19ss4TStJh75kfcqNmF5voMMSERGpEUpGaoluZ40ko9+zeE2Dfvs/YM3ciYEOSUREpEYoGalFegy7kY1nPATAgN3/YfX8hwMdkoiISLWrUjKSkpJCXFwcERERJCYmsnLlypOOf/HFF+ncuTP16tWjU6dOzJ8/v6rxhrykyyaTFj8BgAHbZ7FuyaxAhyQiIlKt/E5GFi5cyMSJE7n//vvZtGkTgwYNYtiwYeTk5JQ7fs6cOUydOpVHHnmEH3/8kUcffZTbb7+djz766FTEH5L6jHmMdS3HAJD03SOkf/IGHq/JmqwCPti8izVZBWqsJyIiIcPvrb19+/alV69ezJkzx3esc+fOjBgxgunTp58wfsCAAQwcOJBnnnnGd2zixImkp6ezatWqSn1mqGzt9Yfp9ZL+wlh67/8It2nhDsv9fHq0s+/92OgI5sevIKFZJAyZGtBYRUREylMtW3udTicbNmxg6NChZY4PHTqU1atXl3uOw+EgIqJs35V69eqRlpaGy+Wq8JyioqIyr7rGsFjoddsb5FjbYjO8JHufoKeR6Xv/ikNvk5CRTOa+IwGNU0RE5M/yKxnJz8/H4/EQExNT5nhMTAx79uwp95zzzz+fuXPnsmHDBkzTJD09nXnz5uFyucjPL78E+vTp04mOjva92rRp40+YocNi5RrLM/zibU6Y4eHtsMfpaOzgTutSX3n5sVmD9chGRESCWpUWsBqGUeZn0zRPOFbqwQcfZNiwYfTr1w+73c4ll1zC9ddfD4DVai33nKlTp1JYWOh77dixoyphBr207P3sKPYwzPkku7yNqWe4+CzsPl/DvWTPSHILS0jL3h/oUEVERKrMr2SkadOmWK3WE2ZB8vLyTpgtKVWvXj3mzZvHkSNH+OWXX8jJyaF9+/Y0bNiQpk2blntOeHg4UVFRZV51UV5xCQBHiWCY80m8JhgGmCYs8Zx5wjgREZFg5FcyEhYWRmJiIqmpqWWOp6amMmDAgJOea7fbad26NVarlXfffZe///3vWCwqc3IyzRv+ttbmOuvnWI4nIoYBH4dNpRkHThgnIiISbPzOBiZPnszcuXOZN28eW7ZsYdKkSeTk5DB+/Hg4/ohl7NixvvHbtm3jrbfeIjMzk7S0NK666ip++OEHnnjiiVP7TUJQn7jGxEZHMMG61Pdopr9jNge9kfzFcphPwqbQqWEJfeIaBzpUERGRKrP5e8KoUaMoKChg2rRp5Obm0q1bN5YvX067du0AyM3NLVNzxOPxMGPGDLZu3YrdbmfIkCGsXr2a9u3bn9pvEoKsFuPY9t2MY4tVZ3tGAvB31xN8EjaVppZiFrgmc+jgIKIbNw90uCIiIlXid52RQKiLdUZ8vpxO5r4jjM0aTG7hb2tDkhrkM999L5E42GrrRMs7P6VhtGZIRESk9qjs/VvJSJDweE3SsveTV1xC84YR9IlrzI6f1tPovUtpxCEy7N1oP/ETIuvXzesjIiK1j5KROuLnzV8T8/6VNOQo34Un0nHi/yOiXmSgwxIREameCqxS+5zW40x2X/QmR8xwznBsYEvySJwOR6DDEhERqTQlIyGgU+/z+GXoazhMOz2PruG72aNwV1BqX0REpLZRMhIiugwcztbBKThNK0mHvmTjC9fi9XgCHZaIiMgfUjISQs4YciUZA57DbVroU/gpaSnjML1ePF6TNVkFfLB5F2uyCtTLRkREahW/64xI7dbj/OvY4DxKz/Qp9CtYRuZT2xhjTmNP0W/rSGKjI47VL2kWCUOmBjReERERzYyEoMTh49l4xsMAJDh+ZObRB8q8f8Wht0nISCZz35EARSgiIvIbJSMhquelE3mSY92RB1i3MN8+HYA7rUuZbD9W0XVs1mA9shERkYDTY5oQlZa9n5dKhmKxHuWf9oWcaf2eny3XYjO8zCgtLV9YQlr2fvrHNwl0uCIiUodpZiRE5RUfKx2f4rmEZPcIAGyGF7dp8fW4+f04ERGRQFEyEqKaN4zw/dll/jYBZjO8vGp/ttxxIiIigaBkJET1iWtMbHQEE6xLudu+mBmuy3nZfREA51k38op9BrHRx3rciIiIBJLWjIQoq8U4tn0349hi1WOPZkyseLnJ9glDrRtoFf4cVss5gQ5VRETqOM2MhLCEZpFkdpnAogajjx8xeNx9LW+awwDoUriSdUuTAxqjiIiIuvbWAR6vSVr2fvKKS2jeMILe7Rqx8eV/0GffYrymQXrPx+kz4o5AhykiIiGmsvdvPaapA6wW44Ttu71vfZW0FC998peStOkB0gwrfS65NWAxiohI3aXHNHWUYbHQ+7a5pDW5BIthkrhxKus/einQYYmISB2kZKQOMyxWkm57nfWNh2M1THqlT2H9/5sb6LBERKSOUTJSx1msVhJv/w/r/3IRVsOk5/p7SV8+L9BhiYhIHaJkRI4lJHe8SXqjYdgMLz3W3c2GT97A4zVZk1XAB5t3sSarQH1sRESkWmg3jfh43W42zh5NUuFneEyDeyx3s+xoL9/7sdERx2qXNIuEIVMDGquIiNR+lb1/a2ZEfCw2Gz3vfJtfbHFYDZNnvTM417LB9/4Vh94mISOZzH1HAhqniIiEFiUjUpbFyjXG02zxtsFqmLxkn8nZlo3caV3KZPuxaq5jswbrkY2IiJwySkakjLTs/ewqdvF35xP85G2NzTB5zf6sr79NsmckuYUlpGXvD3SoIiISIpSMSBl5xSUAeLBykXM6HtPAMMA0YZvZ+oRxIiIif5aSESmjecMI359vs36A1TB9CUmK/Xkusqw9YZyIiMifoWREyugT15jY6AgmWJf6Hs0kON7kR29brIbJbHsyV0espU9c40CHKiIiIULJiJRhtRjMj1/hW6w62zMSLxaGO5/ge097LAb820wm/cMXAx2qiIiECCUjcoKEZpFkdpnAogajfce8WLil3rNsD+uIxYDemx5g3ZJZAY1TRERCg4qeSYU8XpO07P3kFZfQvGEEfeIaY8Fk/Zyb6LNvCQBru9xPvyv/GehQRUSkFqrs/VvJiPjN9HpJe3k8ffcuBGBNpyn0v1oVWUVEpCxVYJVqY1gs9LnlJdbFXgtA/61PsuatRwMdloiIBCklI1IlhsVCn5tns7b1DQD0/3kmq//zYKDDEhGRIKRkRKrMsFjoN+451rb9BwADspNZPe8+OL7eRB1/RUSkMrRmRE6JtW9Mpd8vKQBsjejBdd4H2VPk8L2vjr8iInWP1oxIjep3/XTWxd8FQKeSzcw6ej/wW56rjr8iIlIRJSNyyiRd8yhPMxaAftafWGR/FDDV8VdERE7KFugAJHSkZe8npeQCDlkNptn/Q2/rNrIs12I1TGYcr+bK8Y6//eObBDpcERGpJTQzIqdMaSff+Z7z+ZdrHICv0d6LnhEnjBMREUHJiJxKv+/k24RCAEzzWEKyPGwKNtwnjBMREVEyIqdMeR1/b3Pdhds0ON2yk8/C7qNtlEUdf0VEpAwlI3LKlNfx9xNvX25y3YvLtBBvyeVd110cPVwY6FBFRKQWUTIip1R5HX9XeHtwm/VhnKaNluZedj5/AYX78wMap4iI1B5VSkZSUlKIi4sjIiKCxMREVq5cedLxCxYsoHv37kRGRhIbG8sNN9xAQUFBVWOW2mzIVBKufIxV953NOzf34/mrevDOzf146YG7+PXixRRSn9PdW8h/8TwK8nYFOloREakF/E5GFi5cyMSJE7n//vvZtGkTgwYNYtiwYeTk5JQ7ftWqVYwdO5Zx48bx448/smjRItavX89NN910KuKXWspqMegf34RLerSif3wTrBaDhMQh7L9iGQVEE+/ZzqGXhrJ31/ZAhyoiIgHmdzIyc+ZMxo0bx0033UTnzp2ZNWsWbdq0Yc6cOeWOX7t2Le3bt2fChAnExcXxt7/9jVtuuYX09PRTEb8EmbiufTl6zUfspQntvDvxzL2AXdt/CnRYIiISQH4lI06nkw0bNjB06NAyx4cOHcrq1avLPWfAgAHs3LmT5cuXY5ome/fuZfHixVx00UUVfo7D4aCoqKjMS0JH64TueG/4hF1GC1qae7HPH8avWzcFOiwREQkQv5KR/Px8PB4PMTExZY7HxMSwZ8+ecs8ZMGAACxYsYNSoUYSFhdGiRQsaNWrE7NmzK/yc6dOnEx0d7Xu1adPGnzAlCMS260T4zZ/zi6UNzdlPw3cu5ufvVqvbr4hIHVSlBayGYZT52TTNE46VysjIYMKECTz00ENs2LCBTz/9lOzsbMaPH1/h7586dSqFhYW+144dO6oSptRyTVu2o9GtqfxsjacxRbRecjE3Pp7C1a+u5a53N3P1q2v521NfkPneg/Dl9ECHKyIi1cSvZKRp06ZYrdYTZkHy8vJOmC0pNX36dAYOHMi9997LGWecwfnnn09KSgrz5s0jNze33HPCw8OJiooq85LQ1KhZLDF3prLXaEqE4eJVz4P0t/zoe1/dfkVEQp9fyUhYWBiJiYmkpqaWOZ6amsqAAQPKPefIkSNYLGU/xmq1wvEZFZHIqMaMsj3PL97mhBke3rRPZ4hlk7r9iojUEX4/ppk8eTJz585l3rx5bNmyhUmTJpGTk+N77DJ16lTGjh3rGz98+HCWLl3KnDlz2L59O9988w0TJkygT58+tGzZ8tR+GwlKadn7+aXY4Hzn0/zsbYnN8DLP/oyvpHyyZyS5x7v9iohI6LH5e8KoUaMoKChg2rRp5Obm0q1bN5YvX067du0AyM3NLVNz5Prrr6e4uJgXXniBu+++m0aNGnH22Wfz1FNPndpvIkGrtIuvgzAucD7J1vDrsBompgmF1D9hnIiIhBbDDIJnJUVFRURHR1NYWKj1IyFoTVYBV7+6FoA7jzfZ85gGVuPYX81Z7pHMcl/GOzf3p398kwBHKyIilVXZ+7d600jAldftN97xFqs9XQCYaFvKv8Pnk9g2OtChiohINVAyIgFXXrdfMBjteoD/8/QE4BrjM75LvgKnQ49qRERCjZIRqRXK6/YL8GDkg2T85Ww8pkFS8Rdsee4iDhcXBixOERE59bRmRGoVj9ckLXs/ecUlNG8YQZ+4xlgtBt+vWEz8l7cRaTj4ydaZFuM/oFHT8mvbiIhI7VDZ+7eSEQkaW9P/jxb/byzRHCbb0pbIcR8Q06pDoMMSEZEKaAGrhJxOSedycNSH5NGYOG8O3leHkpP5XaDDEhGRP0nJiASVdp2T8NzwKTuMlsSyjwYLLiJz8zdw/BGPmuyJiAQfPaaRoLR/704OvHIx8Z4snKaNzJhh3FR4A7mFv+22iY2OYH78ChKaRcKQqQGNV0SkLtJjGglpjWNa03zC/5ERdgZhhpuueR9x3+Fny4xRkz0RkeCgZESCVsPoxrS76xNSvb0BGGFbzcv2mXC8kqua7ImIBAclIxLUvtvjYLxzAu+5zwLgfGs6P4dfqyZ7IiJBRMmIBLW84hI8WPmn+x+85P47ADbDi8c0eMEzosw4ERGpnZSMSFBr3jDi+J8MDpsRvuNWw+TTsCmE4/yfcSIiUtsoGZGgVl6TvTudd+A2DTpZdvJ/YfcQ39BNn7jGgQ5VREQqYAt0ACJ/RmmTvYSM3zfZg3xXNK/bn6KNJZ9FrjvZt6sbLdrEBzpcEREph2ZGJOiV12RvjbcrN9qe4hARNKYQ47Xz2J6xPqBxiohI+VT0TEJGeU329u3MxPHGSNp5d1BEfXKGvka3AcMCHaqISJ2gRnkixxUW7CX3pRGc7srAYdr5od8zJA67IdBhiYiEPFVgFTkuukkM7Selsqn+3wg3XPRcO4m17zwR6LBEROQ4JSNSJ0RENuCMSR+wrulILIZJv61PseblO/B6PGqwJyISYHpMI3WK6fWy7s0H6Jf9IgC/2OK41niSncUe3xg12BMROTX0mEakHIbFQr/rnmB9j3/jMQ3au7N5y3En9TnqG6MGeyIiNUvJiNRJvS6+nQmWf+E0rbS35PFF2N0046Aa7ImIBICSEamT0rL38/HRrlzmfJTDZjgxloOkhd+mBnsiIgGgZETqpNLGed+bHRjmfBLTBMMA04Q0b+cTxomISPVRMiJ10u8b511i+QbDAK9pYBiwIOxxRlhWnTBORESqh5IRqZPKa7DX2fE6W72tsRkms8JSuCf8fXq3axToUEVEQp6SEamTShvslS5Wne0ZiYMwLnA+SZqnEwB3GO+xcfZonA49qhERqU5KRqTOKq/BnomFuyKnsyVqIF4T+hR+yrYZQyncvy+gsYqIhDIVPZM6r7wGe1aLwXcrFhP/5e3UN0r41dIa27VLaNXh9ECHKyISNNQoT+QU2P7DWhosvprm7KeAaPKH/4dOiUMCHZaISFBQBVaRU6BDt34YN39BlrUDTSik7YdXsvHT+YEOS0QkpCgZEfkDzVrF0eKuL/iuXh/qGU56rJnAmremYXq9arInInIK6DGNSCW5XU42vnQzfQreB+DnsE5cZz7OrmKXb4ya7ImI/EaPaUROMZs9jN63v8660yYBcJpzK/MdE4jkt62/arInIuI/JSMifjAsFpJGP8w9xmRcpoV4Sy5fhk2mOQfUZE9EpIqUjIj4KS17P4uPJnGl82Ffk7114beryZ6ISBUpGRHxU2nzvE1mAhf8T5O9bWabE8aJiMjJKRkR8dPvm+eN+J8mey+HPcdt1g8AU032REQqScmIiJ/Ka7KX4JjPJk88AP+0L2RW2Muc0ULJiIhIZSgZEfFTeU32PFi51PUY//X0BGCE5WtyZp1Lwd6dgQ5XRKTWUzIiUgXlNdkDeCDyQX5sdhEObHR2ZeB4aQjbM9ICFqeISDBQ0TORP6GiJns52zZjfWcUrcw9HDYjyDzzeXqcc1WgwxURqVFqlCcSYIUFe9j58hV0dX6HxzRY33Eyfa9+AMOiCUkRqRuqtQJrSkoKcXFxREREkJiYyMqVKysce/3112MYxgmvrl27VuWjRYJGdJMWdLwnlbTGw7EaJv0yZ5A2eyxOR4l62oiI/I7fMyMLFy5kzJgxpKSkMHDgQF5++WXmzp1LRkYGbdu2PWF8YWEhR48e9f3sdrvp3r07d955J4888kilPlMzIxLMTK+XtHceJ2nbTKyGSa7RnLHWp8k8FOYbo542IhKKqm1mZObMmYwbN46bbrqJzp07M2vWLNq0acOcOXPKHR8dHU2LFi18r/T0dA4cOMANN9zg70eLBCXDYqHvNQ/x41kv4zBtxJp5vOe8nXhjl2+MetqISF3mVzLidDrZsGEDQ4cOLXN86NChrF69ulK/47XXXuPcc8+lXbt2FY5xOBwUFRWVeYkEu66Dr2Ss7WkKvZH8xXKYj8P+xd8s36unjYjUeX4lI/n5+Xg8HmJiYsocj4mJYc+ePX94fm5uLp988gk33XTTScdNnz6d6Oho36tNmzYnHS8SDNKy97PucAvOds5gp7cJEYaLN+3T1dNGROq8Ki1gNQyjzM+maZ5wrDxvvPEGjRo1YsSIEScdN3XqVAoLC32vHTt2VCVMkVqltFdNAdGc7ZyJ53gJeYDWRj5huMqMExGpK2z+DG7atClWq/WEWZC8vLwTZkv+l2mazJs3jzFjxhAWFnbSseHh4YSHh/sTmkit9/teNbdYP8JqmLhNCzbDyyjbCk6z7GK8c6J62ohInePXzEhYWBiJiYmkpqaWOZ6amsqAAQNOeu5XX33Fzz//zLhx46oWqUiQK6+nzWmOt1jsHgRAoiWTj8IfoNGB7wIdqohIjfL7Mc3kyZOZO3cu8+bNY8uWLUyaNImcnBzGjx8Pxx+xjB079oTzXnvtNfr27Uu3bt1OTeQiQaa8njYA97hv5TX3BQC0MA7Q4aMrWP/+iwGOVkSk5vj1mAZg1KhRFBQUMG3aNHJzc+nWrRvLly/37Y7Jzc0lJyenzDmFhYUsWbKE559//tRFLhKEfD1tsgZD4W9rQ+bW/wcD27YjOvtjYt076b35X6zJ/Z7eNyVjs5/8saaISLBTOXiRAKiop43X4yHt9Xvpt/M1AL4P70Xbf7xLdJOTr8kSEamN1JtGJIht+uR1Oq29j0jDwU6jBa4rFhDXJSnQYYmI+KVae9OISPXqOewG9lzxEbuN5rQ299B84UVs/PwtOD6ror42IhJKNDMiUosd2JfL7ldH0dX5LQBb6/Xkevf95BY7fWPU10ZEaivNjIiEgL80i6XjPamsbXYFAJ2ObuI/jruI5LfFr+prIyLBTsmISC1nDwun962v8jC34jYtdLTs4quwibQx9qqvjYiEBL+39opIzUvL3s9/SgbxvRHD/LAnaWYp4uuwSRgGzCitWXK8r03/+CaBDldExC+aGREJAqX9ajaaHTnH8SxeEwwDTBO8WDDwlhknIhJMlIyIBIHf96u50roCi4Gv0d699vd42f4cDTmivjYiEpSUjIgEgfL62sQ7FvC5JxGAodYNfBD+IM0d2YEOVUTEb0pGRIJARX1t/uG6mwXuswHoYOQSu/BCNiyfF+BoRUT8o2REJEj4+to0GF3m+Av172BL/DjyrDFEGg4S0yax9qXbcLucFf4uEZHaREXPRIJMRX1t3C4n6+dNon/usUqtP4Z1p8W4t2kS0zrQIYtIHaXeNCJ11MZPXqfT2inUN0rYSxMKL55Hx16DK0xiRESqi5IRkTrs1y0bsLx3LW3M3XhMg63RAxnnmExu4W9bf1VGXkSqm8rBi9Rh7Ton0mjiN2yKHIjVMOlStIpXjkwmnN/WkaiMvIjUFkpGREJUw+jGdJv0IbPMq/Ga8FfrL6wIm0RL8lVGXkRqFZWDFwlh6TmFzHIMZ4OlHa/YZxJrOcA34RNURl5EahXNjIiEsNLy8Cu9Z3Ce85kyZeQNUBl5EakVlIyIhLDfl4e/1LKyTBn5yfbF/Mf+FI0pUhl5EQkoJSMiIayiMvKfuHsDcKb1ez4O/xfR+RsCHaqI1GFKRkRCWEVl5G91T+I/7vMAiDX2k7D8KtYumIbp9QY4YhGpi5SMiIS4isrIv1T/Vn5KuIXdYe2xGx76Zc5g04yLKTpYELBYRaRuUtEzkTqiogqsptdL2qKn6ZnxNGGGh51GLI5L5xF/xoBAhywiQU4VWEXEL9s2riDqw5towT4cpp1vz3iA3pdOwLBYVEpeRKpEyYiI+O1g/h5+fW0M3Y+mAbArPJ5D7c7j+l+HqpS8iPhN5eBFxG+Nmrbgr/d8ypq42/GYBq0cWXTa9hI3H365zDiVkheRU0nJiIiUYbFa6X/dE2Sc9xb7zGgAbrR9xmxbMoBKyYvIKadkRETKdSi2Pxc6nmCttzMAw21ryQq/xlevJNkzktzjpeRFRP4MJSMiUq684hL28Reucf6LF90XA2A1TLwmfOTtX2aciMifoWRERMpVWiLeg5USMww41tPGYsBnYfdxseWbMuNERKpKyYiIlKu8UvL9HC+ww9uUcMNNctiLzAh7lb82U/NvEflzlIyISLnKKyW/l8YMdj7HGs+xdSSXWb4kf9ZAfslYH+hwRSSIKRkRkQqVV0reg5XJkf/mx+bDOUwE7b07aLFwGGmLZ6q3jYhUiYqeicgfqqgCa8Henex6/TrOKEkHYEPDs+l402s0jG4c6JBFpBZQBVYRqRFej4e0BY+QmPUidsPDTqMFJRe/ymk9z4STJDIiEvqUjIhIjfpp/X+J/vgWYtmH07SS16QPR2KSuG77EJWSF6mjVA5eRGrU6b3PIXLCGjbV/xthhofW+9fQcctsrj30RplxKiUvIv9LyYiInDLRjZvR4+6PWHP6FBzmsS2/t9s/5EnbK6BS8iJSASUjInJKGRYL9P4HI53T2O5tAcBVthX8HH6tSsmLSLmUjIjIKZdXXMKPZnuGO//NUs/fALAZXrymwRLPmWXGiYgoGRGRU660RPxh6pF9fHbkWCl5ky/C7+Yiy9oy40SkblMyIiKnXHml5M9yPsdub2MiDBcvhiUzI+wVujbVP0EiomRERKpBeaXkc8wYznTOYq3ndAAus6ygcNYAMjd9HehwRSTAlIyISLUor5S8GxuTIp/gx+YXc4h6tDF30/79EayZ/yBejyeg8YpI4FSp6FlKSgrPPPMMubm5dO3alVmzZjFo0KAKxzscDqZNm8Zbb73Fnj17aN26Nffffz833nhjpT5PRc9EgldFFVgL9+eRNW8cvQ4dmxn5IbwHzce+QfNWcYEOWUROkWqrwLpw4ULGjBlDSkoKAwcO5OWXX2bu3LlkZGTQtm3bcs+55JJL2Lt3L48//jinnXYaeXl5uN1uBgwYcEq/jIgEF9PrJf392XT99t9EGg4O0oDt/Z+k1/ljVEZeJARUWzLSt29fevXqxZw5c3zHOnfuzIgRI5g+ffoJ4z/99FOuuuoqtm/fTuPGVWuepWREJLTtyPyWkndvJMHzMwDb7QncwCP8Wvxb8qEy8iLBp1rKwTudTjZs2MDQoUPLHB86dCirV68u95wPP/yQpKQknn76aVq1akXHjh255557OHr0aIWf43A4KCoqKvMSkdDVJqE77f75DWtir8VrQgdXJssct9DV+MU3RmXkRUKXX8lIfn4+Ho+HmJiYMsdjYmLYs2dPueds376dVatW8cMPP7Bs2TJmzZrF4sWLuf322yv8nOnTpxMdHe17tWnTxp8wRSQIhYVH0OfmF7jV+jDFZgSNLYf4IOx+brJ+zATrEpWRFwlhVdpNYxhln9uapnnCsVJerxfDMFiwYAF9+vThwgsvZObMmbzxxhsVzo5MnTqVwsJC32vHjh1VCVNEgkxa9n4+O9KJQY7nyfS2xGaYPGBfwGT7El5xXagy8iIhyq9kpGnTplit1hNmQfLy8k6YLSkVGxtLq1atiI6O9h3r3Lkzpmmyc+fOcs8JDw8nKiqqzEtEQl9pefiDNOQ85zO4zd/+ibrK9iUjLKsAU2XkRUKMX8lIWFgYiYmJpKamljmemppa4c6YgQMHsnv3bg4dOuQ7tm3bNiwWC61bt65q3CISgn5fHv5O6zJshhenaQUgyjjKrLAUXrAnE+3VOjKRUOL3Y5rJkyczd+5c5s2bx5YtW5g0aRI5OTmMHz8ejj9iGTt2rG/86NGjadKkCTfccAMZGRl8/fXX3Hvvvdx4443Uq1fv1H4bEQlq5ZWR7+h4k+dclwHgMQ3+bl1H1w+H8d2KxYEOV0ROEZu/J4waNYqCggKmTZtGbm4u3bp1Y/ny5bRr1w6A3NxccnJyfOMbNGhAamoqd955J0lJSTRp0oQrr7ySxx9//NR+ExEJeqVl5BMyfisjD/C85zJMDCbbF3OAhjTjAM1WjGPdDx9zxg3J1KvfMNChi8ifUKUKrDVNdUZE6pAvp5O57whjswaTW/jb2pDSOiNxjWykZ+6g375FAOwwWnJ0+Bw69hocwKBFpDzVVvQsEJSMiNQ9f1SB9fuvlxHzxWSasx+3aSG97TgSx/wbe1i4qreK1BJKRkQk5BUW7OXnN8aTWPwFAAVGY/a2OZ9xe68od1ZF1VtFala1VGAVEalNopvEkHj3MtKTnqGI+jQx99Ml5x0eO/IY8Nt/Z6l6q0jtpmRERIJe0t//weEbv2a1txsA51o38WXYZJpzgDutS1W9VaSWUzIiIiHhF9dfuMY5hUdcY3GZFuIse1kXfrtvi7Cqt4rUXkpGRCQk5BWXYGLhDc8FDHM+ideE0i4V3Sy/0JRC3zgRqV2UjIhISPh99dZhljQsBr5y8udb0/k8/F4usqwtM05EagclIyISEsqr3nqa4y3mu88FoLFxiBfDkglfdgMH9uUGOlwR+R0lIyISEkqrt5YuVi2t3vqQ+0ZmuY792WMa9Dr0Fd4X+7Lp8zcDHLGIlPK7HLyISG2V0CySzC4TWJQ1GH5XZ2Rhg2u5KL41jT37KM78hvbeHJqsvoP0Hz6g43UpRDVpDpUotCYi1UNFz0Qk5JwsqXCUHGHj/Cn02TUfq2GSTyPcbf9GcYMOFZagV7E0kaqp7P1bMyMiEnKsFoP+8U3KfS88IpL+/0hm64ZLifj4Dtp5d0LO/6MFcI3rV55llG/ssWJpi8nsMoGEGoxfpK7RmhERqZM6JQ4h5p51rI4Zjdc8Nmtyh/0DZthSAFQsTaQGKRkRkTorIrIBxtDHucL5ENneGAAus60iK/waFUsTqUFKRkSkTssrLmGD2YlhzieZ574AAKthYprwvRlXZpyIVA8lIyJSp5UWQSshnANmAwDM49Vb3wh7hpn2FBpRrGJpItVIyYiI1GnlFUvr7Hid9Z6OAIy0ruL/wu8lbOsHx7IUETnllIyISJ1WXrG0EsK5wvUIC9xnA9DUKCIxbRKbnv07+bm/BjpkkZCjZERE6jxfsbQGo8scf6H+HWztdBs7GvbEZVrpeXgVYS/3Z/37szG93oDFKxJqVPRMROS4kxVLy/p+Hd73byPB8zMA30Uk0Xz0HFq07fiH54rUVZW9fysZERGpJLfLyfp3HqNX1hzCDReHzQgONEviSLOeXLd9iKq3ivyPyt6/9ZhGRKSSbPYw+o99jL3X/pct9q7UN0ponb+Kjltmc/PhV8qMPVa9NZnMfUcCFq9IsFAyIiLip7YJ3ek0ZSVrOt3HYTMcgBttnzLfPh0rHlVvFfGTkhERkSqwWK3Q5xbOdz7N156/AnCm9Xsyw8eoequIn5SMiIhUUV5xCTvNZox1TeFe1z8wTbAYx8qRNDCOUo8S3zgRqZiSERGRKvqtKqtBC/ZjGOAxDQwDbrF9zOdh93GW5VtVbxX5A0pGRESqqLzqrfGOBSx1DwSgjWUf/wl7irD3b6Jg745AhytSaykZERGpovKqtwJMdt/ObNcIALwmJBZ/gW1OX9YvfV7F0kTKoWRERORPqKh669sNxpLZZQIHTr+an63xRHOY3t89RMaTZ7Ej81vfOI/XZE1WAR9s3sWarALtvJE6SUXPREROgZNVYHW5nKS/+296/JxCPcOJw7Szr9EZHI7ty/XZ56pYmoQsFT0TEalBVotB//gmXNKjFf3jm5QpBW+3h9F/zKPsv24l30UkEW64aF24gU4/pTDh8Owyv0fF0qQuUjIiIlJDWnU4nb/+M5W0nk+Rbx77r8SrbV/yvv0BGnJExdKkzlIyIiJSgwyLBU+3KzjX8Qzvuc8CoId1O9+F36RiaVJnKRkREalhecUlHKQh/3TfwtXO+zFNMI4/1elt2Uo7Y49vnEhdoGRERKSG/b4IWpKxFcMAt3nsn+Mzrd/zedh93GVdQuNwbQOWukHJiIhIDSuvWNppjrd4zX0BAOGGi0n2JcS9dx4/fP1+oMMVqXZKRkREalhFxdIec49lputyAA4RQWszl25fXMeGGSPI3/1LgKMWqT5KRkREAqCiYmmLGowms8sE7P1uYW2zK/CYBonFXxLxcj/WvfM4HrfLN1YF0yRUqOiZiEgAnaxYGsDP367C+9EkOrq3HfvZGs9f4pPYb4thbNZgFUyTWq2y929bjUYlIiJllBZLq8hp3f+Gp+sa1i6dRZeMmZzmycK7NYsmBlzj+oVnuco39ljBtMVkdplAQg3FL3Iq6DGNiEgtZ7XZ6HflPThvTSMt6nxKJ07usH/IHNtzgKmCaRLUlIyIiASJpjGt8Vwyh1GOB8n0tgJgmG0928OvVcE0CWpKRkREgkhecQnrzM5c6JzOU66rME2wGCamCY2NYqI47BsnEiyUjIiIBJHSgmkubNhwYxjgMQ0MA26wfcYX4XdzhXUFzerbAx2qSKVVKRlJSUkhLi6OiIgIEhMTWblyZYVjV6xYgWEYJ7x++umnPxO3iEidVF7BtHjHAha7BwHQ1CjiGfsrNFk4nJ83V/xvs0ht4ncysnDhQiZOnMj999/Ppk2bGDRoEMOGDSMnJ+ek523dupXc3FzfKyFBa71FRPxVUcG0e9y3Mst17M8O00Yn9090WDactNljOZi/J8BRi5yc33VG+vbtS69evZgzZ47vWOfOnRkxYgTTp08/YfyKFSsYMmQIBw4coFGjRlUKUnVGRER+58vpZO47UmGdkZYRDrZkZpFU/H8AHKQB27pOJPHSSVhttj+sbSJyqlRLnRGn08mGDRuYMmVKmeNDhw5l9erVJz23Z8+elJSU0KVLFx544AGGDBlS4ViHw4HD4SjzZURE5LghU0kAVpWbVJwDQBLw45pPqJc6hQ7eX+jz4+Psz3iR3JjB3HRgjIqlSa3i12Oa/Px8PB4PMTExZY7HxMSwZ0/504CxsbG88sorLFmyhKVLl9KpUyfOOeccvv766wo/Z/r06URHR/tebdq08SdMEZE6obRg2iU9WtE/vskJsxtd+w+j7dT1rO10H0VmJI3NA3Tds4xXj0yiKYW+cceKpSWTue9IAL6FiJ+PaXbv3k2rVq1YvXo1/fv39x3/97//zZtvvlnpRanDhw/HMAw+/PDDct8vb2akTZs2ekwjIlJFebk7WDnnTi6zfAlAiWnjKffVRHGYSfalzHRdzqIGo1l139l6ZCOnTGUf0/g1M9K0aVOsVusJsyB5eXknzJacTL9+/cjMzKzw/fDwcKKiosq8RESk6rKORHK382ZGOKaxx9uICMPNw/Y3mWRfyhL331QsTQLKr2QkLCyMxMREUlNTyxxPTU1lwIABlf49mzZtIjY21p+PFhGRP6G0CNpm8zQGOF/Abf72z/9ltlXMtT9DnJGrYmkSEH43yps8eTJjxowhKSmJ/v3788orr5CTk8P48eMBmDp1Krt27WL+/PkAzJo1i/bt29O1a1ecTidvvfUWS5YsYcmSJaf+24iISLlKi6UB3G59H5vhxWnaCDPceEyDc62bONPyHSu/WkdR+yeJalRx8z6RU83vZGTUqFEUFBQwbdo0cnNz6datG8uXL6ddu3YA5Obmlqk54nQ6ueeee9i1axf16tWja9eufPzxx1x44YWn9puIiEiFSoulXXHobSYfL5Y22zOSO48XT9vubUEHyx7OOfAe+2d9yvpuE+k14i6sNjV3l+rnd52RQFCdERGRPy/zvQdJyEhm5vGGeqUmHO/4+0PTYTTc/z3tvDsB2G6No+ScJ+gy4Lf/eFSNEvFHtdQZERGR4JXQLJLMLhNYlDUYfldnZFGD0QyPb0m3ZpE4B/yHNYufoWtmCh082fD51WxcfRYJpyWwxxlZYaE11SiRP0MzIyIidUxlZjf25+1m28Kp9M7/AKth4jYt2Awvya4RzPRc6RtXOquS2WUCCVc+FoBvI7VZZe/fSkZERKRC239Yx5EP/0k352bfsU/cvbnNfRd3WN/n7uM9clSjRMpTLXVGRESkbunQrS/FVy7mH85J/OptDsAw23qywq/1dQ1WjRL5s5SMiIjISeUdcvK5tzfnOZ9huutqTBNKJ0C6WH6lnXGsEKZqlEhVaQGriIicVGmNEid2wnBhGOAxDayGyTDres6xbORNz1CizdMCHaoEKc2MiIjISZXWKJlwvCbJDNflxDsW8B/3eQCEGR7G2T6h5wdDWLdgGk6HZkjEP0pGRETkpKwWg/nxK5h8fLHq7OM1Sh5238BM1+UAFBBNNIfpmzmDfU92Z9On/8H0en2/w+M1WZNVwAebd7EmqwCPt9bvnZAapN00IiLyx76cTua+IxXWGYlvEs76A/WJ/2EWTTkIwE/2LjSL+yv7bTGqT1JHaWuviIiccn9Uo+RQ8UG+f3caPXa+ST3D6Tv+qvtC/u2+1vez6pPUDUpGREQkYPbuzOLXRVNJOvg5FuPYbSbN04lxrnu53vqp6pPUEaozIiIiARPTOh7PxXMY7vw3qz1dAOhj3cp34Tdxt30xs1wjVZ9EfJSMiIhItcgrLuFHsz2jXfczznk3pgnG8QmQS62rGG5ZjYFX9UlEyYiIiFSP0vokYNDF+BXDALd57LbTzpLH7LAX+CDsQey/fh3QOCXwlIyIiEi1KK8+yWmOt0h2jQDAYdo4w5LNhZvG8/2TZ7P9+zWBDlkCRMmIiIhUi4rqk8z0XMlM1+WEG25+Djsdp2nlryUb6LDkAjbMvJzcX7eW+T2qURL6VA5eRESqTUKzSDK7TGBR1mD4XZ2RRQ1GMzy+JQnNItnZ9mL2vP8AScVfkFiUinPel6yLuYzu7Zqy47BVNUrqAG3tFRGRavdH9UkAtm1aieOTB/irczMcf4wTbrh53nUpz3mu8I1TjZLgUdn7t2ZGRESk2lktBv3jm5x0TMeegzC7f8m3Xy+j/tePcZo3G4C77MvobMnhVtdEbrN+4HvssyhrMKu8pmqUhACtGRERkVrDsFjoPvgy9l79OXc5b2OHtxkAQ60b+Dl8jG8hrGqUhBYlIyIiUuvkH3bxgfdvnON8lmmuMWVqlJxn3cCZlm8BUzVKQoSSERERqXVKa5Q4sVOfo2VqlJxhyWZ+2FMsDHsM+660AEcqp4KSERERqXUqqlGS4hoOHEtM+lp+4sL11/PdU+epRkmQUzIiIiK1TkU1Sp72XM1M1+XYDC/b7Qm4TQtnHE2jw5IL2DhjBDszvyvze1SjJDhoN42IiNRKlalRktP6IvZ++DC9i7+gV/GXuN86i7QmF9G1XSy7HRGqURIkVGdERERqtcrUKPn5uzUc+uQRehxdC8cf49gML3Ncw3nKc7VvnGqU1KzK3r+VjIiISMjYkpaKN/VRurq+9x1b4+nMLa7JXGf9jLtLa5Q0GM2q+85WjZJqVtn7t9aMiIhIyOjc5zyKRi3jWudUvvV2AKC/dQvfht/M3fbFzHZdoholtZCSERERCSl5h5ys8v6VS5yPcYtzEt7f1Si51vZfxls/JJIS1SipRZSMiIhISCmtUQIGHY0dWAxwHa9R8hfjEFPs7/J1+ETC0lIoOVIc0FjlGCUjIiISUsqrUZLgeIuZrssAOOCtT1OjiGG7X+Dw091Ie+fflBw9XOZ3aEtwzdLWXhERCSmlNUoSMsrWKEn2XAYYTLYvJiNqEI2Kt9LSzKPJ1qfZ99Rcvus8nsTGDrYfcGlLcA1TMiIiIiHnj2qUdGkWiXPAYtZ++CLtf0yhBfk02/IExUSSwBFGuXYyi8t9511x6G0SMo5vCQ7Qdwpl2torIiIhqzI1SkqOHmHzhy8Qt2UOMfy2w+YzTxK3ue7iNusH2hJcRaozIiIi4oevt+zgi7ee5jbbhzQ3DgL4ugU/5xrJ855jMyXv3NyP/vFNAhxtcFCdERERET8ccFh4w3MBZzqe4zHXNb5EBOBy60qusn6BHbe2BFcDJSMiIiK/2xJcQjiRODCMY2XlAdpY9vGkfS4rwidh3/g6jpIjAY42tCgZERERqWBL8GmOt3jedSkAh8wIWhkFXJjzDIVPdiXt3ScoOXLohN+jbcH+024aERGRk2wJfs5zBR6sx7YEN+hPk0PbiKGA5j89Rf5PL/Ntx3F0HzGJiHUvkLnviLYFV4FmRkRERI7zbQluMLrM8UUNRpPZZQJdEs8i+r4fWNvlAXJpSlMO0nfbDA4/3ZWt6f8lISOZyw+9XebcY9uCk8ncp0c7FdFuGhERkf9RmS3BDsdRNn30Mm1+TKGVuReAI2YYkYaT2a5LmOEZxZ3HH/nU1W3B2torIiJSA5wOB5+9m0y3rFeJs+z1HXebFmyGlxm/e+RT17YFa2uviIhIDQgLD8fbfTTnOp9lovM2sryxANgML6YJ9Y0SmnGsbom2BZdPyYiIiMif1LxhBB6svO/9Gx94BgDgPV6nZLzt/7Eq/C6m2V6n/tHdgQ61VqpSMpKSkkJcXBwREREkJiaycuXKSp33zTffYLPZ6NGjR1U+VkREpFb6/bbgyfYlzHBdTgfHApa5BwIQbrgYa0vlrE/PJ33WKHZs23zC76jLW4L9XjOycOFCxowZQ0pKCgMHDuTll19m7ty5ZGRk0LZt2wrPKywspFevXpx22mns3buXzZtP/B+iIlozIiIitV3mew+SkJHMTNflJB9fIwIcT1AWs8vSglbePQB4TYNvGw4iauh9xBesDNktwdW2gLVv37706tWLOXPm+I517tyZESNGMH369ArPu+qqq0hISMBqtfL+++8rGRERkdDy5fQ/TCp+iurHkf8+Q68j3/je32VpSSvv7jILXfldEpPZZQIJVz5W41/nVKjs/duvomdOp5MNGzYwZcqUMseHDh3K6tWrKzzv9ddfJysri7feeovHH3/8Dz/H4XDgcDh8PxcVFfkTpoiISM0bMpUEYFW524LPAeB0gMQhZP24nv2fPUXPwv/SyntsHcnd9sW0N/Zwt/tW7rQuY3LpluCswazymiG9JdivZCQ/Px+Px0NMTEyZ4zExMezZs6fcczIzM5kyZQorV67EZqvcx02fPp1HH33Un9BERERqBavF+MPtu/FdexPfdTGfrVxD/mfPcLn1K8INN5fZVnGpdRUWA55zXXbscU9hCWnZ+0N6S3CVFrAaRtnszDTNE44BeDweRo8ezaOPPkrHjh0r/funTp1KYWGh77Vjx46qhCkiIlKrlTRsy/3ucQxyPM+r7gsxTSidABlpXcm11lQicIT8lmC/ZkaaNm2K1Wo9YRYkLy/vhNkSgOLiYtLT09m0aRN33HEHAF6vF9M0sdlsfP7555x99tknnBceHk54eLj/30ZERCSIlHYKzuMvFJmRvk7BNsNLO0sej1teZ5JtMemrr6SwzT+JbtIi0CFXC79mRsLCwkhMTCQ1NbXM8dTUVAYMGHDC+KioKL7//ns2b97se40fP55OnTqxefNm+vbt++e/gYiISJD6o07Bhd5ImhjFnJ/3GvbkM1iXcjN7c7ad8HuCfVuw3117J0+ezJgxY0hKSqJ///688sor5OTkMH78eDj+iGXXrl3Mnz8fi8VCt27dypzfvHlzIiIiTjguIiJS11SmU/CPjYYQXvwrp3m20zfvPdyvLSa90Tk0HXov7fO+CIltwX4nI6NGjaKgoIBp06aRm5tLt27dWL58Oe3atQMgNzeXnJyc6ohVREQk5Pg6BWcNht8lFIsajGZ4fEu6NovEPGsp3369DMvqZP7q3ExSYSosSmWXpSUJ3t1c7trNbH7bFnysU/DxbcEB+l7+UKM8ERGRWqAynYIBtm38muIvZtCj+Cusxm+38A/d/Znovp3bre/Xmk7B6torIiISwj7/+hv2fj6TK6xfEWG4ADCP98NJdo1gpudKCHCnYHXtFRERCWFHo9rzoPtGBjqSSXaP8CUiAGNtqdxre5dmHAiKbcFKRkRERIJQ6bbgAqJxmTbftmCARsZhbrd9yDfhE2icehe//Liu3N9RW3bh+L2AVURERAKvdFvwFYfeZvLxbcGzPSOZYF3CZPsSdnqb0NpSwKDDxxa7/vBRL+h/O10HjcT46qlatQtHMyMiIiJBqHRbcGkPm9Jtwcmey5jpupzWlgJ+jL2U9AaD8ZgG3Uo20u3LceQ8fgZbNn9DQkYylx96u8zvPLYLJ5nMfUdq9LtoZkRERCRIVWZbMEOmsjP7J3Z8MpMz9n5AO+8OKNzBYTOcu+2LqYeDpz1Xc+fxLsGBaM6n3TQiIiJBrrLbggsPFLBq4bP0zH2Xlsb+3843DayG6XvUwynahVPZ+7dmRkRERIJcZToFA0T/pQnufndw5rv9uNCSxk22jznDko3VMHGbFl8iAtToLhytGREREalDmjeMwI2ND70DSPUkwvGZEZvh5U7r0jLjaopmRkREROqQinbh3Hm8WZ9xfM1Jn7jGNRaTkhEREZE6pKLmfLM9IzGAyfbFDI9vidVyTo3FpGRERESkjvmjXTgJzSJrNB4lIyIiInXNkKkkAKvK3YVTczMipZSMiIiI1FGV3YVT3bSbRkRERAJKyYiIiIgElJIRERERCSglIyIiIhJQSkZEREQkoJSMiIiISEApGREREZGAUjIiIiIiAaVkRERERAIqKCqwmqYJQFFRUaBDERERkUoqvW+X3scrEhTJSHFxMQBt2rQJdCgiIiLip+LiYqKjoyt83zD/KF2pBbxeL7t376Zhw4YYhuHXuUVFRbRp04YdO3YQFRVVbTHKb3TNa56uec3TNa95uuaB8Weuu2maFBcX07JlSyyWileGBMXMiMVioXXr1n/qd0RFRekvbw3TNa95uuY1T9e85umaB0ZVr/vJZkRKaQGriIiIBJSSEREREQmokE9GwsPDefjhhwkPDw90KHWGrnnN0zWvebrmNU/XPDBq4roHxQJWERERCV0hPzMiIiIitZuSEREREQkoJSMiIiISUEpGREREJKBCIhlJSUkhLi6OiIgIEhMTWbly5UnHf/XVVyQmJhIREUGHDh146aWXaizWUOHPNV+6dCnnnXcezZo1Iyoqiv79+/PZZ5/VaLyhwN+/56W++eYbbDYbPXr0qPYYQ42/19zhcHD//ffTrl07wsPDiY+PZ968eTUWbyjw95ovWLCA7t27ExkZSWxsLDfccAMFBQU1Fm+w+/rrrxk+fDgtW7bEMAzef//9PzynWu6hZpB79913Tbvdbr766qtmRkaGedddd5n169c3f/3113LHb9++3YyMjDTvuusuMyMjw3z11VdNu91uLl68uMZjD1b+XvO77rrLfOqpp8y0tDRz27Zt5tSpU0273W5u3LixxmMPVv5e81IHDx40O3ToYA4dOtTs3r17jcUbCqpyzS+++GKzb9++ZmpqqpmdnW2uW7fO/Oabb2o07mDm7zVfuXKlabFYzOeff97cvn27uXLlSrNr167miBEjajz2YLV8+XLz/vvvN5csWWIC5rJly046vrruoUGfjPTp08ccP358mWOnn366OWXKlHLH//Of/zRPP/30MsduueUWs1+/ftUaZyjx95qXp0uXLuajjz5aDdGFpqpe81GjRpkPPPCA+fDDDysZ8ZO/1/yTTz4xo6OjzYKCghqKMPT4e82feeYZs0OHDmWOJScnm61bt67WOENVZZKR6rqHBvVjGqfTyYYNGxg6dGiZ40OHDmX16tXlnrNmzZoTxp9//vmkp6fjcrmqNd5QUJVr/r+8Xi/FxcU0bty4mqIMLVW95q+//jpZWVk8/PDDNRBlaKnKNf/www9JSkri6aefplWrVnTs2JF77rmHo0eP1lDUwa0q13zAgAHs3LmT5cuXY5ome/fuZfHixVx00UU1FHXdU1330KBolFeR/Px8PB4PMTExZY7HxMSwZ8+ecs/Zs2dPuePdbjf5+fnExsZWa8zBrirX/H/NmDGDw4cPc+WVV1ZTlKGlKtc8MzOTKVOmsHLlSmy2oP6/eUBU5Zpv376dVatWERERwbJly8jPz+e2225j//79WjdSCVW55gMGDGDBggWMGjWKkpIS3G43F198MbNnz66hqOue6rqHBvXMSCnDMMr8bJrmCcf+aHx5x6Vi/l7zUu+88w6PPPIICxcupHnz5tUYYeip7DX3eDyMHj2aRx99lI4dO9ZghKHHn7/nXq8XwzBYsGABffr04cILL2TmzJm88cYbmh3xgz/XPCMjgwkTJvDQQw+xYcMGPv30U7Kzsxk/fnwNRVs3Vcc9NKj/k6lp06ZYrdYTsua8vLwTMrdSLVq0KHe8zWajSZMm1RpvKKjKNS+1cOFCxo0bx6JFizj33HOrOdLQ4e81Ly4uJj09nU2bNnHHHXfA8RulaZrYbDY+//xzzj777BqLPxhV5e95bGwsrVq1KtMuvXPnzpimyc6dO0lISKj2uINZVa759OnTGThwIPfeey8AZ5xxBvXr12fQoEE8/vjjmumuBtV1Dw3qmZGwsDASExNJTU0tczw1NZUBAwaUe07//v1PGP/555+TlJSE3W6v1nhDQVWuOcdnRK6//nrefvttPc/1k7/XPCoqiu+//57Nmzf7XuPHj6dTp05s3ryZvn371mD0wakqf88HDhzI7t27OXTokO/Ytm3bsFgstG7dutpjDnZVueZHjhzBYil7G7NarfC7/1qXU6va7qF/avlrLVC6Fey1114zMzIyzIkTJ5r169c3f/nlF9M0TXPKlCnmmDFjfONLtyVNmjTJzMjIMF977TVt7fWTv9f87bffNm02m/niiy+aubm5vtfBgwcD+C2Ci7/X/H9pN43//L3mxcXFZuvWrc3LL7/c/PHHH82vvvrKTEhIMG+66aYAfovg4u81f/31102bzWampKSYWVlZ5qpVq8ykpCSzT58+AfwWwaW4uNjctGmTuWnTJhMwZ86caW7atMm3nbqm7qFBn4yYpmm++OKLZrt27cywsDCzV69e5ldffeV777rrrjPPOuusMuNXrFhh9uzZ0wwLCzPbt29vzpkzJwBRBzd/rvlZZ51lAie8rrvuugBFH5z8/Xv+e0pGqsbfa75lyxbz3HPPNevVq2e2bt3anDx5snnkyJEARB68/L3mycnJZpcuXcx69eqZsbGx5jXXXGPu3LkzAJEHpy+//PKk/z7X1D3UMDWXJSIiIgEU1GtGREREJPgpGREREZGAUjIiIiIiAaVkRERERAJKyYiIiIgElJIRERERCSglIyIiIhJQSkZEREQkoJSMiIiISEApGREREZGAUjIiIiIiAaVkRERERALq/wPpAt+SmxYplQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "s = simulations[0]\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(s.grid.x[s.grid.lo:s.grid.hi+1], s.U_init[s.grid.lo:s.grid.hi+1, ivar], marker=\"o\")\n",
    "ax.plot(s.grid.x[s.grid.lo:s.grid.hi+1], s.U[s.grid.lo:s.grid.hi+1, ivar], marker=\"x\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a519e449-59f3-45f0-af79-8beaac6acdaa",
   "metadata": {},
   "source": [
    "## Convergence of the well-balanced method"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48b3f3e3-334c-4b69-b5d7-0506319e2cab",
   "metadata": {},
   "source": [
    "Now we do the same for the reconstruction that does the characteristic tracing only on\n",
    "the pressure perturbation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "404610cf-5d8e-4e43-a637-29c817d55983",
   "metadata": {},
   "outputs": [],
   "source": [
    "simulations = []\n",
    "for nx in [32, 64, 128, 256, 512]:\n",
    "    dt = 0.015625 * (32 / nx)\n",
    "    e = Euler(nx, 0.5, fixed_dt=dt, init_cond=hse, grav_func=constant_gravity,\n",
    "              use_limiting=True, use_flattening=True,\n",
    "              use_hse_reconstruction=True, hse_as_perturbation=True,\n",
    "              bc_left_type=\"reflect\", bc_right_type=\"reflect\", params=params)\n",
    "    e.evolve(0.5, verbose=False)\n",
    "    simulations.append(e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "f23c236b-e6b8-4980-bfa1-d261a19b4900",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 64 ->  32 : 3.4760439589959424e-05\n",
      "128 ->  64 : 8.776428605358436e-06\n",
      "256 -> 128 : 2.204958615561335e-06\n",
      "512 -> 256 : 5.525999756839299e-07\n"
     ]
    }
   ],
   "source": [
    "for coarse, fine in pairwise(simulations):\n",
    "    _, cd = fine.grid.coarsen(fine.U[:, ivar])\n",
    "    err = coarse.grid.norm(coarse.U[:, ivar] - cd)\n",
    "    print(f\"{fine.grid.nx:3d} -> {coarse.grid.nx:3d} : {err}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6243e975-5d5a-4ff2-98eb-3305a4fdecb8",
   "metadata": {},
   "source": [
    "As expected, this also converges 2nd-order."
   ]
  }
 ],
 "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.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}