Use single metallicity cooling tables

parthenon-hpc-lab · Aug 30, 2024 · a8cbd3a · a8cbd3a
1 parent 85895c4
commit a8cbd3a
Show file tree

Hide file tree

Showing 19 changed files with 1,473 additions and 351 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -6,10 +6,10 @@
 - [[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/...)
-- [[PR 97]](https://github.com/parthenon-hpc-lab/athenapk/pull/97) Fixed Schure cooling curve. Removed SD one. Added description of cooling function conventions.
 
 ### Infrastructure
 - [[PR 112]](https://github.com/parthenon-hpc-lab/athenapk/pull/112) Add dev container configuration
@@ -21,7 +21,9 @@
 ### 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.
-  - Added `updated_schure.cooling` curve (and associated notebook to calculate it from the paper tables).
+  - 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)

diff --git a/docs/cooling/cooling.ipynb b/docs/cooling/cooling.ipynb
@@ -122,15 +122,14 @@
    "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) for Z=1 solar metalicity, Z=0.5,\n",
-    "and Z=0.3 across the columns. Cooling rates are in the convention Lambda_hd from eq. 1,\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.5 n_i/n_i(solar) is set to 1 for H and He and set to 0.5 and 0.3 respectively for\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",
-    "log10 T [K] Z=1 log10 Lambda_N [erg cm^3/s] Z=0.5 log10 Lambda_N [erg cm^3/s] Z=0.3 log10 Lambda_N [erg cm^3/s]\"\"\""
+    "-----------------------------------------------------------------------------------------\\n\"\"\""
    ]
   },
   {
@@ -140,7 +139,59 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "np.savetxt(\"updated_schure.cooling\",schure_data.T,header=header,fmt=(\"%1.2f\",\"%2.4f\",\"%2.4f\",\"%2.4f\"))"
+    "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\"))"
    ]
   },
   {

diff --git a/docs/cooling/plot_lambda.py b/docs/cooling/plot_lambda.py
@@ -25,7 +25,14 @@
 
     ## Gnat and Sternberg (2007) [NOTE: n_H n_e convention for Lambda]
     # Columns: Temperature Lambda(Z=1e-3) Lambda(Z=1e-2) Lambda(Z=1e-1) Lambda(Z=1) Lambda(Z=2)
-    gs_T, gs_LambdaZem3, gs_LambdaZem2, gs_LambdaZem1, gs_LambdaZ1, gs_LambdaZ2 = (
+    (
+        gs_T,
+        gs_LambdaZem3,
+        gs_LambdaZem2,
+        gs_LambdaZem1,
+        gs_LambdaZ1,
+        gs_LambdaZ2
+    ) = (
         np.loadtxt("gnat_sternberg_cie_table.txt", unpack=True, skiprows=23)
     )
     # NOTE: There is not enough information in this table alone to compute n_e!

diff --git a/docs/input.md b/docs/input.md
@@ -216,15 +216,17 @@ than any other timestep constraint (e.g., the hyperbolic one).
 ### Cooling
 
 Tabular cooling (e.g., for optically thin cooling) is enabled through the `cooling` block in the input file.
+The tabulated table itself is a text file containing (apart from comments) only two columns (two floating
+point numbers per line with the first one being the log10 temperature and the second one the
+log10 cooling rate scaled to a source function with $S_{cool} = n_H^2 \Lambda(T)$).
+
 A possible block might look like:
 
 ```
 <cooling>
-enable_cooling = tabular                   # To disable, set to `none`
-table_filename = updated_schure.cooling    # Path to the cooling table (in a text file)
-log_temp_col = 0                           # Column in the file that contains the log10 temperatures
-log_lambda_col = 1                         # Column in the file that contains the cooling rates
-lambda_units_cgs = 1                       # Conversion factor of the cooling rate relative to CGS units
+enable_cooling = tabular                # To disable, set to `none`
+table_filename = schure.cooling_1.0Z    # Path to the cooling table (in a text file)
+lambda_units_cgs = 1                    # Conversion factor of the cooling rate relative to CGS units
 
 integrator = townsend              # Other possible options are `rk12` and `rk45` for error bound subcycling
 #max_iter = 100                    # Max number of iteration for subcycling. Unsued for Townsend integrator

diff --git a/inputs/cloud.in b/inputs/cloud.in
@@ -62,9 +62,7 @@ He_mass_fraction = 0.24
 
 <cooling>
 enable_cooling = none                     # turn on with "tabular"
-table_filename = updated_schure.cooling
-log_temp_col = 0
-log_lambda_col = 1
+table_filename = schure.cooling_1.0Z
 lambda_units_cgs = 1
 
 integrator = rk12  

diff --git a/inputs/cluster/cluster.in b/inputs/cluster/cluster.in
@@ -86,9 +86,7 @@ code_time_cgs = 3.15576e+16             # 1 Gyr in s
 
 <cooling>
 enable_cooling=tabular
-table_filename=updated_schure.cooling
-log_temp_col=0
-log_lambda_col=1
+table_filename=schure.cooling_1.0Z
 lambda_units_cgs=1.0
 
 integrator=rk45

diff --git a/inputs/cluster/cooling.in b/inputs/cluster/cooling.in
@@ -70,9 +70,7 @@ code_time_cgs = 3.15576e+16             # 1 Gyr in s
 
 <cooling>
 enable_cooling = tabular
-table_filename = updated_schure.cooling
-log_temp_col = 0
-log_lambda_col = 1
+table_filename = schure.cooling_1.0Z
 lambda_units_cgs = 1 #erg cm^3/s in cgs, as used in schure.cooling
 
 integrator = rk12

diff --git a/inputs/cooling_tables/gnat-sternberg.cooling b/inputs/cooling_tables/gnat-sternberg.cooling