add cooling documentation

CHANGELOG.md

@@ -6,6 +6,7 @@
 - [[PR 1]](https://github.com/parthenon-hpc-lab/athenapk/pull/1) Add isotropic thermal conduction and RKL2 supertimestepping
 
 ### Changed (changing behavior/API/variables/...)
+- [[PR 97]](https://github.com/parthenon-hpc-lab/athenapk/pull/97) Fixed Schure cooling curve. Removed SD one. Added description of cooling function conventions.
 - [[PR 84]](https://github.com/parthenon-hpc-lab/athenapk/pull/84) Bump Parthenon to latest develop (2024-02-15)
 
 ### Fixed (not changing behavior/API/variables/...)
@@ -18,6 +19,11 @@
 ### Removed (removing behavior/API/varaibles/...)
 
 ### Incompatibilities (i.e. breaking changes)
+- [[PR 97]](https://github.com/parthenon-hpc-lab/athenapk/pull/97)
+  - Removes original `schure.cooling` cooling curve as it had unknown origin.
+  - To avoid confusion, only cooling table for a single solar metallicity are supported
+    from now on (i.e., the parameters to specify temperature and lambda columns have been removed).
+  - Added `schure.cooling_#Z` curves (and associated notebook to calculate it from the paper tables).
 - [[PR 84]](https://github.com/parthenon-hpc-lab/athenapk/pull/84) Bump Parthenon to latest develop (2024-02-15)
   - Updated access to block dimension: `pmb->block_size.nx1` -> `pmb->block_size.nx(X1DIR)` (and similarly x2 and x3)
   - Update access to mesh size: `pmesh->mesh_size.x1max` -> `pmesh->mesh_size.xmax(X1DIR)` (and similarly x2, x3, and min)

docs/cooling/cooling.ipynb

@@ -0,0 +1,369 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "33ca3da4-e1e2-4dff-bf33-88072e57508e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import unyt\n",
+    "\n",
+    "import scipy as sp\n",
+    "import scipy.interpolate\n",
+    "\n",
+    "import matplotlib.pyplot as plt\n",
+    "import matplotlib as mpl"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "eadb3f5c-e2aa-44b9-8828-555cac0c8c42",
+   "metadata": {},
+   "source": [
+    "## Load data from [Schure 2009](https://doi.org/10.1051/0004-6361/200912495) Tables 2 and 4"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6e9bfaf9-9e06-48bf-933d-b3e201e4c5f5",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "elemental_rates_data = np.genfromtxt(\"schure_table_4.txt\",\n",
+    "                                skip_header=4,\n",
+    "                                dtype=float, delimiter='\\t', names=True) \n",
+    "elements = elemental_rates_data.dtype.names[1:]\n",
+    "\n",
+    "log_temperature = elemental_rates_data[\"T\"]\n",
+    "\n",
+    "cooling_rates = unyt.unyt_array(np.empty( (log_temperature.size,len(elements)) ),\"cm**3*erg/s\")\n",
+    "\n",
+    "for j,element in enumerate(elements):\n",
+    "    cooling_rates[:,j] = unyt.unyt_array(10**elemental_rates_data[element],\"cm**3*erg/s\")\n",
+    "\n",
+    "\n",
+    "schure_table2_data = np.loadtxt(\"schure_table_2.txt\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "0dda7ed5-24d0-4ab5-a703-7baa71d1e1b8",
+   "metadata": {},
+   "source": [
+    "## Compute Cooling rates ($\\Lambda_N$) from Schure Table 4 and eq. 3"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "dd992955-c664-4dc5-9a40-0c25ebafdc88",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "z1_log_cr = np.log10( cooling_rates.sum(axis=1) )\n",
+    "\n",
+    "def compute_cr(n_over_n_sun):\n",
+    "    return (cooling_rates*n_over_n_sun).sum(axis=1)\n",
+    "\n",
+    "z05_log_cr = np.log10( compute_cr(np.array( [1,1.0] + [1./2.,]*(len(elements)-2) ) ))\n",
+    "z03_log_cr = np.log10( compute_cr(np.array( [1,1] + [0.3,]*(len(elements)-2) ) ))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "59b73548-47cf-438d-bc8b-41a691b7f5e8",
+   "metadata": {},
+   "source": [
+    "## Compute $\\Lambda_{hd}$ from $n_e/n_H$ in Table 2 and eq. 1"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "17526d49-46cd-41ce-b12f-f682de1cdba3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#ne_over_nh = sp.interpolate.CubicSpline(schure_table2_data[:,0],schure_table2_data[:,3])\n",
+    "\n",
+    "start = np.where(schure_table2_data[:,0] == log_temperature[0])[0][0]\n",
+    "end = np.where(schure_table2_data[:,0] == log_temperature[-1])[0][0]\n",
+    "ne_over_nh = schure_table2_data[start:end+1,3]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "12c2e08a-721a-4535-87b6-8d0688627d5d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "log_x_hd = np.log10(ne_over_nh)\n",
+    "schure_data = np.vstack((log_temperature,log_x_hd+z1_log_cr,log_x_hd+z05_log_cr,log_x_hd+z03_log_cr))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "abb949f6-e0eb-43d6-80db-9242e6b11ded",
+   "metadata": {},
+   "source": [
+    "## Save the schure data with a thorough note"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9783bc26-4023-41f4-b7ad-f5dbbab8fc7d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "header=\"\"\"Cooling rates generated from Schure 2009 (doi.org/10.1051/0004-6361/200912495)\n",
+    "containing temperatures in the first column (in log10 K) and collisional ionisation\n",
+    "equilibrium (CIE) cooling rates (in log10 erg cm^3/s).\n",
+    "Cooling rates are in the convention Lambda_hd from eq. 1,\n",
+    "where the proton ratio n_e/n_H is taken from table 2.\n",
+    "Lambda_N is computed from eq. 3. The cooling rate Lambda_N(X_i,T) from eq. 3 is contained\n",
+    "in table 4. n_i/n_i(solar) is taken to be 1.0 for all elements for Z=1 while for Z=0.5\n",
+    "and Z=0.3 n_i/n_i(solar) is set to 1 for H and He and set to 0.5 and 0.3 respectively for\n",
+    "all other elements. Made by Forrest Glines ([email protected])\n",
+    "-----------------------------------------------------------------------------------------\\n\"\"\""
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b82b388b-b1db-42c0-b491-ec0acf998dd3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "for i, z in enumerate([1.0, 0.5, 0.3]):\n",
+    "    np.savetxt(f\"schure.cooling_{z:.1f}Z\",\n",
+    "               schure_data.T[:,[0, i + 1]],\n",
+    "               header=header + f\"log10 T [K] Z={z:.1f} log10 Lambda_N [erg cm^3/s]\",\n",
+    "               fmt=(\"%1.2f\",\"%2.4f\"))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ebd552ea-937d-4d26-9517-cffaf98dbeb4",
+   "metadata": {},
+   "source": [
+    "### Process Gnat Sternberg Cooling Table"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1a9b0fd3-a27e-46eb-b356-4d46e53bcc55",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "gnat_sternberg_data = np.loadtxt(\"gnat_sternberg_cie_table.txt\", skiprows=23)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fdf5d5db-a0eb-43e7-9283-29339d1d6d5a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "header=\"\"\"Adapted from: http://wise-obs.tau.ac.il/~orlyg/cooling/CIEcool/tab13.txt\n",
+    "Title: Time-Dependent Ionization in Radiatively Cooling Gas \n",
+    "Authors: Orly Gnat and Amiel Sternberg\n",
+    "Table: CIE Cooling Efficiencies\n",
+    "-----------------------------------------------------------------------------------------\n",
+    "Our assumed Z=1 solar abundances are listed in Table 1.\n",
+    "-----------------------------------------------------------------------------------------\\n\"\"\""
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0abce35b-6775-4c43-b6a8-2ed0cd6f0e1c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "for i, z in enumerate([1e-3, 1e-2, 1e-1, 1, 2]):\n",
+    "    np.savetxt(f\"gnat-sternberg.cooling_{z:.1g}Z\",\n",
+    "               np.log10(gnat_sternberg_data[:,[0, i + 1]]),\n",
+    "               header=header + f\"log10 T [K] Z={z:.1g} log10 Lambda_N [erg cm^3/s]\",\n",
+    "               fmt=(\"%1.2f\",\"%2.4f\"))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "6c7c385d-21ed-42e5-ba60-8ae76cea7185",
+   "metadata": {},
+   "source": [
+    "## Load other cooling tables (See [Sutherland and Dopita 1993](https://ui.adsabs.harvard.edu/link_gateway/1993ApJS...88..253S/doi:10.1086/191823) )"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9177f991-65fc-46a9-a8e1-8c121da5f68f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "enzo_schure_data = np.loadtxt(\"enzo_schure.cooling\") #Schure cooling table with enzo roots\n",
+    "#gnat_sternberg_data = np.loadtxt(\"gnat-sternberg.cooling\")\n",
+    "#sutherland_dopita_data = np.loadtxt(\"cooling_data/sutherland_dopita.cooling\") #SD table with PLUTO roots\n",
+    "sutherland_dopita_table6_data = np.loadtxt(\"sutherland_dopita_table_6.txt\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "702382fd-8067-4a7f-bb85-45838f236ff1",
+   "metadata": {},
+   "source": [
+    "## Make some plots"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ef31148e-36d4-4ff7-b950-712e84224f12",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig,ax = plt.subplots()\n",
+    "\n",
+    "ax.plot(log_temperature,z1_log_cr,label=\"$\\\\Lambda_N$ Schure $Z_\\\\odot$\")\n",
+    "ax.plot(log_temperature,z05_log_cr,label=\"$\\\\Lambda_N$ Schure $0.5 Z_\\\\odot$\")\n",
+    "ax.plot(log_temperature,z03_log_cr,label=\"$\\\\Lambda_N$ Schure $0.3 Z_\\\\odot$\")\n",
+    "\n",
+    "\n",
+    "ax.plot(log_temperature,log_x_hd+z1_log_cr,label=\"$\\\\Lambda_{hd}$ Schure $Z_\\\\odot$\",linestyle=\"--\")\n",
+    "ax.plot(log_temperature,log_x_hd+z05_log_cr,label=\"$\\\\Lambda_{hd}$ Schure $0.5 Z_\\\\odot$\",linestyle=\"--\")\n",
+    "ax.plot(log_temperature,log_x_hd+z03_log_cr,label=\"$\\\\Lambda_{hd}$ Schure $0.3 Z_\\\\odot$\",linestyle=\"--\")\n",
+    "\n",
+    "#ax.plot(enzo_schure_data[:,0],enzo_schure_data[:,1],label=\"Legacy Schure $Z_\\\\odot$\",linestyle=\"--\")\n",
+    "#ax.plot(enzo_schure_data[:,0],enzo_schure_data[:,2],label=\"Legacy Schure $0.5 Z_\\\\odot$\",linestyle=\"--\")\n",
+    "#ax.plot(schure_table2_data[:,0],schure_table2_data[:,1],label=\"Schure Table 2 $Z=1$\",linestyle=\"--\")\n",
+    "\n",
+    "#ax.plot(gnat_sternberg_data[:,0],gnat_sternberg_data[:,4],label=\"GS $Z=1$\",linestyle=\":\")\n",
+    "#ax.plot(sutherland_dopita_data[:,0],sutherland_dopita_data[:,2],label=\"SD $Z=1$\",linestyle=\"-\")\n",
+    "#ax.plot(sutherland_dopita_data[:,0],sutherland_dopita_data[:,1],label=\"SD $Z=1/3$\",linestyle=\":\")\n",
+    "\n",
+    "#ax.plot(sutherland_dopita_table6_data[:,0],sutherland_dopita_table6_data[:,5],label=\"SD Table 6 $Z=1$\",linestyle=\"--\")\n",
+    "\n",
+    "ax.legend(ncols=2)\n",
+    "\n",
+    "#ax.set_xlim(log_temperature.min(),log_temperature.max())\n",
+    "\n",
+    "ax.set_xlim(3.8,8.2)\n",
+    "ax.set_ylim(-23,-20)\n",
+    "#ax.set_ylim(-26,-20)\n",
+    "\n",
+    "ax.set_ylabel(\"Cooling Rate -- $ \\\\log \\\\Lambda$ [ cm${}^{3}$ erg/s ]\")\n",
+    "ax.set_xlabel(\"Temperature -- $ \\\\log T$ [ K ]\")\n",
+    "\n",
+    "ax.grid()\n",
+    "\n",
+    "plt.show()\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0261a25f-2c1a-49b7-987e-3c5bf1deedf5",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig,ax = plt.subplots()\n",
+    "\n",
+    "ax.plot(log_temperature,log_x_hd+z1_log_cr,label=\"My $\\\\Lambda_{hd}$\",linestyle=\"-\")\n",
+    "ax.plot(log_temperature,schure_table2_data[start:end+1,2],label=\"Schure Table 2 $\\\\Lambda_{hd}$\",linestyle=\":\")\n",
+    "\n",
+    "ax.plot(log_temperature,z1_log_cr,label=\"My $\\\\Lambda_{N}$\",linestyle=\"-\")\n",
+    "ax.plot(log_temperature,enzo_schure_data[:,1],label=\"Enzo $\\\\Lambda_{N}$\",linestyle=\"--\")\n",
+    "ax.plot(log_temperature,schure_table2_data[start:end+1,2],label=\"Schure Table 1 $\\\\Lambda_{N}$\",linestyle=\":\")\n",
+    "\n",
+    "ax.legend()\n",
+    "\n",
+    "#ax.set_xlim(log_temperature.min(),log_temperature.max())\n",
+    "\n",
+    "ax.set_xlim(3.8,8.2)\n",
+    "ax.set_ylim(-23,-20)\n",
+    "#ax.set_ylim(-26,-20)\n",
+    "\n",
+    "ax.set_ylabel(\"Cooling Rate -- $ \\\\log \\\\Lambda$ [ cm${}^{3}$ erg/s ]\")\n",
+    "ax.set_xlabel(\"Temperature -- $ \\\\log T$ [ K ]\")\n",
+    "\n",
+    "ax.grid()\n",
+    "\n",
+    "plt.show()\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "693b185a-a4c4-4602-b841-b43bfa909601",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig,ax = plt.subplots()\n",
+    "\n",
+    "ax.plot(log_temperature,\n",
+    "        (10**(log_x_hd+z1_log_cr) - 10**schure_table2_data[start:end+1,2])/10**(log_x_hd+z1_log_cr),\n",
+    "        label=\"$\\\\Lambda_{hd}$: Mine - Schure Table 2\",linestyle=\"-\")\n",
+    "\n",
+    "ax.plot(log_temperature,\n",
+    "        (10**(z1_log_cr) - 10**schure_table2_data[start:end+1,1])/10**(z1_log_cr),\n",
+    "        label=\"$\\\\Lambda_{N}$: Mine - Schure Table 2\",linestyle=\"--\")\n",
+    "ax.plot(log_temperature,\n",
+    "        (10**(z1_log_cr) - 10**enzo_schure_data[:,1])/10**(z1_log_cr),\n",
+    "        label=\"$\\\\Lambda_{N}$: Mine - Enzo\",linestyle=\"-\")\n",
+    "\n",
+    "\n",
+    "ax.legend()\n",
+    "\n",
+    "#ax.set_xlim(log_temperature.min(),log_temperature.max())\n",
+    "\n",
+    "ax.set_xlim(3.8,8.2)\n",
+    "#ax.set_ylim(-23,-20)\n",
+    "#ax.set_ylim(-26,-20)\n",
+    "\n",
+    "ax.set_ylabel(\"Relative difference in $\\\\Lambda$\")\n",
+    "ax.set_xlabel(\"Temperature -- $ \\\\log T$ [ K ]\")\n",
+    "\n",
+    "ax.grid()\n",
+    "\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fb9b70d5-32e4-43db-9fcd-f1c7c481b5a4",
+   "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.12.5"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}