Merge pull request #551 from BrandonJLazard/add_checkpoint_tutorial

Added jupyter notebook for reading in checkpoint files

CHANGELOG.md

@@ -27,6 +27,8 @@ This project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.htm
 
 - Added a function called checkpoint_read to the rayleigh_diagnostics.py file to read in checkpoint files    \[Rathish Ratnasingam; 6-21-2024; [#543](https://github.com/geodynamics/Rayleigh/pull/543)\]
 
+-Added a jupyter notebook to read in the checkpoint files, transform it and output the results.         \[Brandon Lazard; 6-21-2024; [#551](https://github.com/geodynamics/Rayleigh/pull/551)]
+
 
 ### Changed
 

post_processing/checkpoint_tutorial.ipynb

@@ -0,0 +1,174 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Tutorial: Reading Checkpoint File\n",
+    "\n",
+    "\n",
+    "This notebook will demonstrate how to read in a checkpoint file, and use the `rayleigh_diagnostics.py` and `spectral_utils.py` to analyze the file. \n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# import the relevant utility\n",
+    "import spectral_utils\n",
+    "\n",
+    "# import checkpoint reading function from Rayleigh\n",
+    "from rayleigh_diagnostics import checkpoint_read\n",
+    "\n",
+    "# import other helpful things\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "\n",
+    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
+    "from matplotlib import cm,colors\n",
+    "plt.rcParams['image.cmap'] = 'seismic'\n",
+    "plt.rcParams['image.origin'] = 'lower'\n",
+    "plt.rcParams['image.interpolation'] = 'none'"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "By default, Rayleigh outputs the checkpoint files in spectral form. Using the `spectral_utils.py` we can easily convert the spectral data into real data. Replace `checkpoint_file` with the path to your checkpoint file"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#Read in the checkpoint data\n",
+    "checkpoint_file = 'sample_checkpoint_T'\n",
+    "\n",
+    "#Specify nr and ntheta of the simulation\n",
+    "nr = 48\n",
+    "ntheta = 48\n",
+    "\n",
+    "#Read in the checkpoint data\n",
+    "checkpoint = checkpoint_read(checkpoint_file, nr, ntheta)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#Specify a few variables based on ntheta\n",
+    "nphi = 2*ntheta\n",
+    "lmax = (nphi) / 3 - 1\n",
+    "dealias = ntheta/(lmax + 1)\n",
+    "\n",
+    "#Initialize the Spectral Transform Class\n",
+    "SHT = spectral_utils.SHT(ntheta, spectral=False, dealias=dealias)\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### First we will perform the spectral transformation in $\\theta$ and $\\phi$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#Create an empty array\n",
+    "spec = np.zeros((ntheta, nphi, nr), dtype = \"float64\")\n",
+    "\n",
+    "#Loop over radius \n",
+    "for i in range(nr):\n",
+    "    \n",
+    "    #We specify .transpose() as ntheta and nphi are swapped in the to_physical function of the SHT class\n",
+    "    spec[:,:,i] = SHT.to_physical(checkpoint[:,:,i].transpose(), th_l_axis=0, phi_m_axis=1)\n",
+    "    "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Next we will perform the Chebyshev transform in $r$. You can change `rmin` and `rmax` as according to your simulation"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#First perform a chebyshev transform in radius to get to real space\n",
+    "rmin=5e10\n",
+    "rmax=6.83177e10\n",
+    "\n",
+    "\n",
+    "#Initialize the Chebyshev Class\n",
+    "cheb = spectral_utils.Chebyshev(nr, rmin=rmin, rmax=rmax)\n",
+    "\n",
+    "#Perform the transformation in radius \n",
+    "phys = cheb.to_physical(spec, axis = 2)\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "You now have real function, $f(\\theta, \\phi, r)$ from your checkpoint file. You can plot and see what it looks like for a particular grid point in radius. `r=0` corresponds to the top of the shell. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#Specify where in the shell to view our function\n",
+    "r = 10\n",
+    "\n",
+    "fig, ax = plt.subplots(1, 1, figsize=(16, 8))\n",
+    "\n",
+    "#Show the image\n",
+    "p0 = ax.imshow(phys[:,:,r], extent=[0, 360, -90, 90])\n",
+    "ax.set_xlabel('Longitude')\n",
+    "ax.set_ylabel('Latitude')\n",
+    "divider = make_axes_locatable(ax)\n",
+    "cax = divider.append_axes('right', size='2%', pad=0.05)\n",
+    "fig.colorbar(p0, cax = cax)\n"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "radev",
+   "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.9.15"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}