diff --git a/doc/RTD/README.md b/doc/RTD/README.md
index 3394ce7b8b97a71a7f14fb235e49e2efb51b9f9f..1d0582386390ab20d48f835805cdaa43e8da496e 100644
--- a/doc/RTD/README.md
+++ b/doc/RTD/README.md
@@ -1,15 +1,14 @@
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!
-
-
diff --git a/doc/RTD/source/Planetary/equations_of_state.rst b/doc/RTD/source/Planetary/equations_of_state.rst
index 58a374817f52ecee55fabdd0bccd3e7f1636c303..4b0d2dfd988234048350fccb22b7e7608e5c35bd 100644
--- a/doc/RTD/source/Planetary/equations_of_state.rst
+++ b/doc/RTD/source/Planetary/equations_of_state.rst
@@ -1,28 +1,28 @@
.. 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.
diff --git a/doc/RTD/source/Planetary/hydro_scheme.rst b/doc/RTD/source/Planetary/hydro_scheme.rst
index ae920722c1fff88d23bf26ed85d534b05d1e6a42..f0c787dce100142ecbe9b650fbee4d2110f34676 100644
--- a/doc/RTD/source/Planetary/hydro_scheme.rst
+++ b/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.
diff --git a/doc/RTD/source/Planetary/index.rst b/doc/RTD/source/Planetary/index.rst
index f3d2134ea3537166cfe9167e6e6a6f31b9cca861..b851bc51144702511c31aea9f89c58e3e15ea3fa 100644
--- a/doc/RTD/source/Planetary/index.rst
+++ b/doc/RTD/source/Planetary/index.rst
@@ -1,25 +1,26 @@
.. 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>
diff --git a/doc/RTD/source/Planetary/initial_conditions.rst b/doc/RTD/source/Planetary/initial_conditions.rst
index fd348a9fd134fca3a35be4a7eeb16669ea56534a..aad89d5060f0ea72864f5fc91c2b69d7ef5cff52 100755
--- a/doc/RTD/source/Planetary/initial_conditions.rst
+++ b/doc/RTD/source/Planetary/initial_conditions.rst
@@ -1,23 +1,22 @@
.. 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.
diff --git a/examples/Planetary/EarthImpact/earth_impact.yml b/examples/Planetary/EarthImpact/earth_impact.yml
index 4833f9ffbf90b9478a65463122ea7d0dadcfb369..39f561770efce4fbd41183a7462371f79bd8b95e 100644
--- a/examples/Planetary/EarthImpact/earth_impact.yml
+++ b/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
diff --git a/src/equation_of_state/planetary/sesame.h b/src/equation_of_state/planetary/sesame.h
index 51b923c28b85296050186e9bc7c7a3feadbb4265..aff3052399780b0f678bd26bcc76b2c02083326a 100644
--- a/src/equation_of_state/planetary/sesame.h
+++ b/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,