Tidy EoS table format and documentation

doc/RTD/README.md

 SWIFT Documentation
 ===================
 
-This is the main documentation for SWIFT that can be found on ReadTheDocs.
+This is the main documentation for SWIFT that can be found on ReadTheDocs
+or at `swiftsim.com/docs`.
 
 You will need the `sphinx` and `sphinx-autobuild` python packages (pip install
 them!) to build the documentation to html, as well as the `sphinx_rtd_theme`
 package which is used as the theme.
 
-To build the documentation, `make html` and then it is available in
-`build/html`.
+To build the documentation (from this directory), `make html` and the output
+files will be created in `build/html/`.
 
 Please consider adding documentation when you add code!
-
-

doc/RTD/source/Planetary/equations_of_state.rst

 .. Planetary EoS
-    Jacob Kegerreis, 13th March 2020
+    Jacob Kegerreis, 14th July 2022
 
 .. _planetary_eos:
 
 Planetary Equations of State
 ============================
 
-Configuring SWIFT with the ``--with-equation-of-state=planetary`` and 
-``--with-hydro=planetary`` options enables the use of multiple 
+Configuring SWIFT with the ``--with-equation-of-state=planetary`` and
+``--with-hydro=planetary`` options enables the use of multiple
 equations of state (EoS).
-Every SPH particle then requires and carries the additional ``MaterialID`` flag 
-from the initial conditions file. This flag indicates the particle's material 
-and which EoS it should use. 
+Every SPH particle then requires and carries the additional ``MaterialID`` flag
+from the initial conditions file. This flag indicates the particle's material
+and which EoS it should use.
 
-It is important to check that the EoS you use are appropriate 
+It is important to check that the EoS you use are appropriate
 for the conditions in the simulation that you run.
-Please follow the original sources of these EoS for more information and 
+Please follow the original sources of these EoS for more information and
 to check the regions of validity. If an EoS sets particles to have a pressure
 of zero, then particles may end up overlapping, especially if the gravitational
 softening is very small.
 
-So far, we have implemented several Tillotson, ANEOS, SESAME, 
+So far, we have implemented several Tillotson, ANEOS, SESAME,
 and Hubbard \& MacFarlane (1980) materials, with more on the way.
-The material's ID is set by a somewhat arbitrary base type ID 
+The material's ID is set by a somewhat arbitrary base type ID
 (multiplied by 100) plus an individual value:
 
 + Ideal gas: ``0``
@@ -45,32 +45,56 @@ The material's ID is set by a somewhat arbitrary base type ID
     + Forsterite (Stewart et al. 2019): ``400``
     + Iron (Stewart, zenodo.org/record/3866507): ``401``
     + Fe85Si15 (Stewart, zenodo.org/record/3866550): ``402``
-    
-The data files for the tabulated EoS can be downloaded using 
+
+The data files for the tabulated EoS can be downloaded using
 the ``examples/Planetary/EoSTables/get_eos_tables.sh`` script.
 
 To enable one or multiple EoS, the corresponding ``planetary_use_*:``
 flag(s) must be set to ``1`` in the parameter file for a simulation,
-along with the path to any table files, which are set by the 
+along with the path to any table files, which are set by the
 ``planetary_*_table_file:`` parameters,
 as detailed in :ref:`Parameters_eos` and ``examples/parameter_example.yml``.
 
-The data files for the tabulated EoS can be downloaded using 
-the ``examples/EoSTables/get_eos_tables.sh`` script.
-
-Unlike the EoS for an ideal or isothermal gas, these more complicated materials 
-do not always include transformations between the internal energy, 
-temperature, and entropy. At the moment, we have implemented 
-\\(P(\\rho, u)\\) and \\(c_s(\\rho, u)\\), 
-which is sufficient for the :ref:`planetary_sph` hydro scheme, 
-but makes most materials currently incompatible with e.g. entropy-based schemes.
+Unlike the EoS for an ideal or isothermal gas, these more complicated materials
+do not always include transformations between the internal energy,
+temperature, and entropy. At the moment, we have implemented
+\\(P(\\rho, u)\\) and \\(c_s(\\rho, u)\\) (and more in some cases),
+which is sufficient for the :ref:`planetary_sph` hydro scheme,
+but some materials may thus currently be incompatible with
+e.g. entropy-based schemes.
 
-The Tillotson sound speed was derived using 
+The Tillotson sound speed was derived using
 \\(c_s^2 = \\left. ( \\partial P / \\partial \\rho ) \\right|_S \\)
-as described in 
-`Kegerreis et al. (2019)  <https://doi.org/10.1093/mnras/stz1606>`_. 
+as described in
+`Kegerreis et al. (2019)  <https://doi.org/10.1093/mnras/stz1606>`_.
 Note that there is a typo in the sign of
 \\(du = T dS - P dV = T dS + (P / \\rho^2) d\\rho \\) in the appendix;
-the correct version was used in the derivation.
+the correct version was used in the actual derivation.
 
 The ideal gas uses the same equations detailed in :ref:`equation_of_state`.
+
+The data files for the tabulated EoS can be downloaded using
+the ``examples/EoSTables/get_eos_tables.sh`` script.
+
+The format of the data files for SESAME, ANEOS, and similar-EoS tables
+is similar to the SESAME 301 (etc) style. The file contents are:
+
+.. code-block:: python
+
+    # header (12 lines)
+    version_date                                                (YYYYMMDD)
+    num_rho  num_T
+    rho[0]   rho[1]  ...  rho[num_rho]                          (kg/m^3)
+    T[0]     T[1]    ...  T[num_T]                              (K)
+    u[0, 0]                 P[0, 0]     c[0, 0]     s[0, 0]     (J/kg, Pa, m/s, J/K/kg)
+    u[1, 0]                 ...         ...         ...
+    ...                     ...         ...         ...
+    u[num_rho-1, 0]         ...         ...         ...
+    u[0, 1]                 ...         ...         ...
+    ...                     ...         ...         ...
+    u[num_rho-1, num_T-1]   ...         ...         s[num_rho-1, num_T-1]
+
+The ``version_date`` must match the value in the ``sesame.h`` ``SESAME_params``
+objects, so we can ensure that any version updates work with the git repository.
+The header contains a first line that gives the material name, followed by the
+same 11 lines printed here to describe the contents.

doc/RTD/source/Planetary/hydro_scheme.rst

@@ -6,12 +6,15 @@
 Planetary Hydro Scheme
 ======================
 
-This scheme is based on :ref:`minimal` but also allows multiple materials, 
-meaning that different SPH particles can be assigned different 
-`equations of state <equations_of_state.html>`_ (EoS). 
-Every SPH particle then requires and carries the additional ``MaterialID`` flag 
-from the initial conditions file. This flag indicates the particle's material 
-and which EoS it should use. 
+This scheme is based on :ref:`minimal` but also allows multiple materials,
+meaning that different SPH particles can be assigned different
+`equations of state <equations_of_state.html>`_ (EoS).
+Every SPH particle then requires and carries the additional ``MaterialID`` flag
+from the initial conditions file. This flag indicates the particle's material
+and which EoS it should use.
 
-The Balsara viscosity switch is used by default, but can be disabled by 
+The Balsara viscosity switch is used by default, but can be disabled by
 compiling SWIFT with ``make CFLAGS=-DPLANETARY_SPH_NO_BALSARA``.
+
+Note: the boundary-improvement methods presented in Ruiz-Bonilla+2022 are
+in the process of being improved for release with better documentation etc soon.

doc/RTD/source/Planetary/index.rst

 .. Planetary Simulations
-   Jacob Kegerreis, 3rd October 2020
+   Jacob Kegerreis, 14th July 2022
 
 .. _planetary:
-   
+
 Planetary Simulations
 =====================
 
 SWIFT is also designed for running planetary simulations
-with a current focus on giant impacts, as introduced in 
-`Kegerreis et al. (2019) <https://doi.org/10.1093/mnras/stz1606>`_, MNRAS 487:4.
+such as of giant impacts, as introduced in Kegerreis+2019,
+and any other types of simulations with more complicated equations of state
+and/or multiple materials, etc.
 
-Active development is ongoing of more features and examples for planetary 
-simulations, so please let us know if you are interested in using SWIFT 
-or have any implementation requests. 
+Active development is ongoing of more features and examples,
+so please let us know if you are interested in using SWIFT
+or have any implementation requests.
 
 You can find an example simulation in ``swiftsim/examples/Planetary/``
-under ``EarthImpact/``, as well as some hydro tests for comparison with other 
-schemes. The tabulated equation of state files can be downloaded using 
+under ``EarthImpact/``, as well as some hydro tests for comparison with other
+schemes. The tabulated equation of state files can be downloaded using
 ``EoSTables/get_eos_tables.sh``.
 
-Planetary simulations are currently intended to be run with 
+Planetary simulations are currently intended to be run with
 SWIFT configured to use the planetary hydrodynamics scheme and equations of state:
 ``--with-hydro=planetary`` and ``--with-equation-of-state=planetary``.
 These allow for multiple materials to be used,
@@ -28,7 +29,7 @@ chosen from the several available equations of state.
 .. toctree::
    :maxdepth: 1
    :caption: More information:
-   
+
    Hydro Scheme <hydro_scheme>
    Equations of State <equations_of_state>
    Initial Conditions <initial_conditions>

doc/RTD/source/Planetary/initial_conditions.rst

 .. Planetary Initial Conditions
-   Jacob Kegerreis, 13th March 2020
+   Jacob Kegerreis, 14th July 2022
 
 .. _planetary_initial_conditions:
-   
+
 Initial Conditions
 ==================
 
-Tools for the creation of initial conditions are available via 
-the 
-`WoMa <https://github.com/srbonilla/WoMa>`_ and 
-`SEAGen <https://github.com/jkeger/seagen>`_ open-source python packages, 
+Tools for the creation of initial conditions are available via the
+`WoMa <https://github.com/srbonilla/WoMa>`_ and
+`SEAGen <https://github.com/jkeger/seagen>`_ open-source python packages,
 including: creating spherical or spinning planetary (or similar) profiles;
 placing particles to match arbitrary profiles with precise SPH densities;
 and setting the initial target and impactor positions and velocities,
-as presented in 
+as presented in
 `Kegerreis et al. (2019)  <https://doi.org/10.1093/mnras/stz1606>`_ and
-Ruiz-Bonilla et al. (2020).
+`Ruiz-Bonilla et al. (2020)  <https://doi.org/10.1093/mnras/staa3385>`_ .
 
-They are available with documentation and examples at 
+They are available with documentation and examples at
 https://github.com/srbonilla/WoMa and https://github.com/jkeger/seagen,
 or can be installed directly with ``pip``
 (https://pypi.org/project/woma/, https://pypi.org/project/seagen/).
@@ -26,13 +25,13 @@ or can be installed directly with ``pip``
 Settling initial conditions with fixed entropies
 ------------------------------------------------
 
-If the particles' equations of state include specific entropies, 
+If the particles' equations of state include specific entropies,
 and the initial conditions file includes specific entropies for each particle
-(in ``PartType0/Entropies``), 
+(in ``PartType0/Entropies``),
 then configuring SWIFT with ``--enable-planetary-fixed-entropy``
-will override the internal energy of each particle each step such that its 
-specific entropy remains constant. 
+will override the internal energy of each particle each step such that its
+specific entropy remains constant.
 
-This should be used with caution, but may be a convenient way to maintain an 
-entropy profile while initial conditions settle to equilibrium with their 
+This should be used with caution, but may be a convenient way to maintain an
+entropy profile while initial conditions settle to equilibrium with their
 slightly different SPH densities.

examples/Planetary/EarthImpact/earth_impact.yml

@@ -15,7 +15,7 @@ InitialConditions:
 TimeIntegration:
     time_begin:     0                   # The starting time of the simulation (in internal units).
     time_end:       36000               # The end time of the simulation (in internal units).
-    dt_min:         0.0001              # The minimal time-step size of the simulation (in internal units).
+    dt_min:         0.000001            # The minimal time-step size of the simulation (in internal units).
     dt_max:         1000                # The maximal time-step size of the simulation (in internal units).
 
 # Parameters governing the snapshots

src/equation_of_state/planetary/sesame.h

@@ -25,7 +25,6 @@
  *
  * Contains the SESAME and ANEOS-in-SESAME-style EOS functions for
  * equation_of_state/planetary/equation_of_state.h
- *
  */
 
 /* Some standard headers. */
@@ -48,56 +47,98 @@ struct SESAME_params {
   float *table_P_rho_T;
   float *table_c_rho_T;
   float *table_log_s_rho_T;
-  int date, num_rho, num_T;
+  int version_date, num_rho, num_T;
   float u_tiny, P_tiny, c_tiny, s_tiny;
   enum eos_planetary_material_id mat_id;
 };
 
-// Parameter values for each material (cgs units)
+// Parameter values for each material
 INLINE static void set_SESAME_iron(struct SESAME_params *mat,
                                    enum eos_planetary_material_id mat_id) {
   // SESAME 2140
   mat->mat_id = mat_id;
-  mat->date = 20201003;
+  mat->version_date = 20220714;
 }
 INLINE static void set_SESAME_basalt(struct SESAME_params *mat,
                                      enum eos_planetary_material_id mat_id) {
   // SESAME 7530
   mat->mat_id = mat_id;
-  mat->date = 20201003;
+  mat->version_date = 20220714;
 }
 INLINE static void set_SESAME_water(struct SESAME_params *mat,
                                     enum eos_planetary_material_id mat_id) {
   // SESAME 7154
   mat->mat_id = mat_id;
-  mat->date = 20201003;
+  mat->version_date = 20220714;
 }
 INLINE static void set_SS08_water(struct SESAME_params *mat,
                                   enum eos_planetary_material_id mat_id) {
   // Senft & Stewart (2008)
   mat->mat_id = mat_id;
-  mat->date = 20201003;
+  mat->version_date = 20220714;
 }
 INLINE static void set_ANEOS_forsterite(struct SESAME_params *mat,
                                         enum eos_planetary_material_id mat_id) {
   // Stewart et al. (2019)
   mat->mat_id = mat_id;
-  mat->date = 20201104;
+  mat->version_date = 20220714;
 }
 INLINE static void set_ANEOS_iron(struct SESAME_params *mat,
                                   enum eos_planetary_material_id mat_id) {
   // Stewart (2020)
   mat->mat_id = mat_id;
-  mat->date = 20201104;
+  mat->version_date = 20220714;
 }
 INLINE static void set_ANEOS_Fe85Si15(struct SESAME_params *mat,
                                       enum eos_planetary_material_id mat_id) {
   // Stewart (2020)
   mat->mat_id = mat_id;
-  mat->date = 20201104;
+  mat->version_date = 20220714;
+}
+
+/*
+    Skip a line while reading a file.
+*/
+INLINE static int skip_line(FILE* f) {
+    int c;
+
+    // Read each character until reaching the end of the line or file
+    do {
+        c = fgetc(f);
+    } while ((c != '\n') && (c != EOF));
+
+    return c;
+}
+
+/*
+    Skip n lines while reading a file.
+*/
+INLINE static int skip_lines(FILE* f, int n) {
+    int c;
+
+    for (int i = 0; i < n; i++) c = skip_line(f);
+
+    return c;
 }
 
-// Read the tables from file
+/*
+    Read the data from the table file.
+
+    File contents (SESAME-like format plus header info etc)
+    -------------
+    # header (12 lines)
+    version_date                                                (YYYYMMDD)
+    num_rho  num_T
+    rho[0]   rho[1]  ...  rho[num_rho]                          (kg/m^3)
+    T[0]     T[1]    ...  T[num_T]                              (K)
+    u[0, 0]                 P[0, 0]     c[0, 0]     s[0, 0]     (J/kg, Pa, m/s, J/K/kg)
+    u[1, 0]                 ...         ...         ...
+    ...                     ...         ...         ...
+    u[num_rho-1, 0]         ...         ...         ...
+    u[0, 1]                 ...         ...         ...
+    ...                     ...         ...         ...
+    u[num_rho-1, num_T-1]   ...         ...         s[num_rho-1, num_T-1]
+*/
 INLINE static void load_table_SESAME(struct SESAME_params *mat,
                                      char *table_file) {
 
@@ -105,30 +146,26 @@ INLINE static void load_table_SESAME(struct SESAME_params *mat,
   FILE *f = fopen(table_file, "r");
   if (f == NULL) error("Failed to open the SESAME EoS file '%s'", table_file);
 
-  // Ignore header lines
-  char buffer[100];
-  for (int i = 0; i < 6; i++) {
-    if (fgets(buffer, 100, f) == NULL)
-      error("Failed to read the SESAME EoS file header %s", table_file);
-  }
-  float ignore;
+  // Skip header lines
+  skip_lines(f, 12);
 
   // Table properties
-  int date;
-  int c = fscanf(f, "%d", &date);
+  int version_date;
+  int c = fscanf(f, "%d", &version_date);
   if (c != 1) error("Failed to read the SESAME EoS table %s", table_file);
-  if (date != mat->date)
+  if (version_date != mat->version_date)
     error(
-        "EoS file %s date %d does not match expected %d"
+        "EoS file %s version_date %d does not match expected %d (YYYYMMDD)."
         "\nPlease download the file using "
         "examples/Planetary/EoSTables/get_eos_tables.sh",
-        table_file, date, mat->date);
+        table_file, version_date, mat->version_date);
   c = fscanf(f, "%d %d", &mat->num_rho, &mat->num_T);
   if (c != 2) error("Failed to read the SESAME EoS table %s", table_file);
 
   // Ignore the first elements of rho = 0, T = 0
   mat->num_rho--;
   mat->num_T--;
+  float ignore;
 
   // Allocate table memory
   mat->table_log_rho = (float *)malloc(mat->num_rho * sizeof(float));
@@ -206,8 +243,8 @@ INLINE static void prepare_table_SESAME(struct SESAME_params *mat) {
           mat->table_log_u_rho_T[i_rho * mat->num_T + i_T - 1]) {
 
         // Replace it and all elements below it with that value
-        for (int j_u = 0; j_u < i_T; j_u++) {
-          mat->table_log_u_rho_T[i_rho * mat->num_T + j_u] =
+        for (int j_T = 0; j_T < i_T; j_T++) {
+          mat->table_log_u_rho_T[i_rho * mat->num_T + j_T] =
               mat->table_log_u_rho_T[i_rho * mat->num_T + i_T];
         }
         break;
@@ -264,10 +301,10 @@ INLINE static void prepare_table_SESAME(struct SESAME_params *mat) {
 INLINE static void convert_units_SESAME(struct SESAME_params *mat,
                                         const struct unit_system *us) {
 
+  // Convert input table values from all-SI to internal units
   struct unit_system si;
   units_init_si(&si);
 
-  // All table values in SI, apart from sound speeds in km/s
   // Densities (log)
   for (int i_rho = 0; i_rho < mat->num_rho; i_rho++) {
     mat->table_log_rho[i_rho] +=
@@ -285,7 +322,7 @@ INLINE static void convert_units_SESAME(struct SESAME_params *mat,
           units_cgs_conversion_factor(&si, UNIT_CONV_PRESSURE) /
           units_cgs_conversion_factor(us, UNIT_CONV_PRESSURE);
       mat->table_c_rho_T[i_rho * mat->num_T + i_T] *=
-          1e3f * units_cgs_conversion_factor(&si, UNIT_CONV_SPEED) /
+          units_cgs_conversion_factor(&si, UNIT_CONV_SPEED) /
           units_cgs_conversion_factor(us, UNIT_CONV_SPEED);
       mat->table_log_s_rho_T[i_rho * mat->num_T + i_T] +=
           logf(units_cgs_conversion_factor(
@@ -301,7 +338,7 @@ INLINE static void convert_units_SESAME(struct SESAME_params *mat,
       units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
   mat->P_tiny *= units_cgs_conversion_factor(&si, UNIT_CONV_PRESSURE) /
                  units_cgs_conversion_factor(us, UNIT_CONV_PRESSURE);
-  mat->c_tiny *= 1e3f * units_cgs_conversion_factor(&si, UNIT_CONV_SPEED) /
+  mat->c_tiny *= units_cgs_conversion_factor(&si, UNIT_CONV_SPEED) /
                  units_cgs_conversion_factor(us, UNIT_CONV_SPEED);
   mat->s_tiny *=
       units_cgs_conversion_factor(&si,