diff --git a/.gitignore b/.gitignore
index 5284b8c7208812c41b9044cd482e1047d3b13fd8..9657e827707dcd56690f67781e1c745c557eafe9 100644
--- a/.gitignore
+++ b/.gitignore
@@ -27,6 +27,7 @@ examples/swift_mpi
examples/*/*.xmf
examples/*/*.h5
examples/*/*.png
+examples/*/*.pdf
examples/*/*.mp4
examples/*/*.txt
examples/*/dependency_graph.csv
diff --git a/doc/RTD/source/EquationOfState/index.rst b/doc/RTD/source/EquationOfState/index.rst
index 3558041e9513b967a2530165acec5e5f4f11a364..11e6069130a84ba54dec01e9892464574d9c2c6b 100644
--- a/doc/RTD/source/EquationOfState/index.rst
+++ b/doc/RTD/source/EquationOfState/index.rst
@@ -1,19 +1,21 @@
-.. Equation of State
+.. Equations of State
Loic Hausammann, 6th April 2018
+ Jacob Kegerreis, 3rd February 2019
.. _equation_of_state:
-Equation of State
-=================
+Equations of State
+==================
-Currently (if the documentation was well updated), we have two different
-equation of states implemented: ideal gas and isothermal. They describe the
-relations between our main thermodynamical variables: the internal energy
-(\\(u\\)), the density (\\(\\rho\\)), the entropy (\\(A\\)) and the pressure
-(\\(P\\)).
+Currently (if the documentation was well updated), we have two different gas
+equations of state (EoS) implemented: ideal and isothermal; as well as a variety
+of EoS for "planetary" materials.
+The EoS describe the relations between our main thermodynamical variables:
+the internal energy (\\(u\\)), the density (\\(\\rho\\)), the entropy (\\(A\\))
+and the pressure (\\(P\\)).
-Equations
----------
+Gas EoS
+-------
In the following section, the variables not yet defined are: \\(\\gamma\\) for
the adiabatic index and \\( c_s \\) for the speed of sound.
@@ -37,12 +39,55 @@ the adiabatic index and \\( c_s \\) for the speed of sound.
"\\( c_s\\)", "", "\\(\\sqrt{ u \\gamma \\left( \\gamma - 1 \\right) } \\)", ""
+
+Planetary EoS
+-------------
+Configuring SWIFT with the ``--with-equation-of-state=planetary`` and
+``--with-hydro=planetary`` options enables the use of multiple 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.
+
+So far, we have implemented several Tillotson, SESAME, and Hubbard \& MacFarlane
+(1980) materials, with more on their way.
+The material's ID is set by a base type ID (multiplied by 100), plus a minor
+type:
+
++ Tillotson (Melosh, 2007): ``1``
+ + Iron: ``100``
+ + Granite: ``101``
+ + Water: ``102``
++ Hubbard \& MacFarlane (1980): ``2``
+ + Hydrogen-helium atmosphere: ``200``
+ + Ice H20-CH4-NH3 mix: ``201``
+ + Rock SiO2-MgO-FeS-FeO mix: ``202``
++ SESAME (and similar): ``3``
+ + Iron (2140): ``300``
+ + Basalt (7530): ``301``
+ + Water (7154): ``302``
+ + Senft \& Stewart (2008) water (in a SESAME-style table): ``303``
+
+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 only implemented
+\\(P(\\rho, u)\\) and \\(c_s(\\rho, u)\\).
+This is sufficient for the simple :ref:`planetary_sph` hydrodynamics scheme,
+but makes these materials currently incompatible with other entropy-based
+schemes.
+
+The Tillotson sound speed was derived using
+\\(c_s^2 = \\left. \\dfrac{\\partial P}{\\partial \\rho} \\right|_S \\)
+as described in Kegerreis et al. (2019).
+The table files for the HM80 and SESAME-style EoS can be downloaded using
+the ``examples/EoSTables/get_eos_tables.sh`` script.
+
+
How to Implement a New Equation of State
----------------------------------------
See :ref:`new_option` for a full list of required changes.
-You will need to provide a ``equation_of_state.h`` file containing: the
+You will need to provide an ``equation_of_state.h`` file containing: the
definition of ``eos_parameters``, IO functions and transformations between the
different variables: \\(u(\\rho, A)\\), \\(u(\\rho, P)\\), \\(P(\\rho,A)\\),
\\(P(\\rho, u)\\), \\(A(\\rho, P)\\), \\(A(\\rho, u)\\), \\(c_s(\\rho, A)\\),
diff --git a/doc/RTD/source/HydroSchemes/index.rst b/doc/RTD/source/HydroSchemes/index.rst
index cd6c169245e83440a1258d216991763488586c0c..462bb7378162ff1addab3212a6901412195a3377 100644
--- a/doc/RTD/source/HydroSchemes/index.rst
+++ b/doc/RTD/source/HydroSchemes/index.rst
@@ -15,6 +15,7 @@ schemes available in SWIFT, as well as how to implement your own.
traditional_sph
minimal_sph
+ planetary
hopkins_sph
gizmo
adding_your_own
diff --git a/doc/RTD/source/HydroSchemes/planetary.rst b/doc/RTD/source/HydroSchemes/planetary.rst
new file mode 100755
index 0000000000000000000000000000000000000000..20f41758baadba2cddb99e79d3435bb3301065e0
--- /dev/null
+++ b/doc/RTD/source/HydroSchemes/planetary.rst
@@ -0,0 +1,26 @@
+.. Planetary SPH
+ Jacob Kegerreis, 3rd February 2019
+
+.. _planetary_sph:
+
+Planetary (Density-Energy, Multi-Material) SPH
+==============================================
+
+.. toctree::
+ :maxdepth: 2
+ :hidden:
+ :caption: Contents:
+
+This scheme is the same as the Minimal SPH scheme but also allows multiple
+materials, meaning that different SPH particles can be assigned different
+:ref:`equation_of_state` (EoS).
+
+To use the planetary scheme and the corresponding planetary EoS, use
+
+.. code-block:: bash
+
+ ./configure --with-hydro=planetary --with-equation-of-state=planetary
+
+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.
\ No newline at end of file
diff --git a/doc/RTD/source/conf.py b/doc/RTD/source/conf.py
index 2249faa2851846c28e743400b2c826bfa6780c0a..fac755bbb4ee9cd25bf3526bc435c69be3a9d5b5 100644
--- a/doc/RTD/source/conf.py
+++ b/doc/RTD/source/conf.py
@@ -18,7 +18,7 @@
# -- Project information -----------------------------------------------------
-project = 'SWIFT: SPH WIth Fine-grained inter-dependent Tasking'
+project = 'SWIFT: SPH With Inter-dependent Fine-grained Tasking'
copyright = '2018, SWIFT Collaboration'
author = 'SWIFT Team'
diff --git a/doc/RTD/source/index.rst b/doc/RTD/source/index.rst
index e04efe8c889fb8a005c88f691f1e01a387f19ebb..83422b4e5caf05bacb3824d06426b9cdeba3921e 100644
--- a/doc/RTD/source/index.rst
+++ b/doc/RTD/source/index.rst
@@ -3,7 +3,7 @@
You can adapt this file completely to your liking, but it should at least
contain the root `toctree` directive.
-Welcome to SWIFT: SPH WIth Fine-grained inter-dependent Tasking's documentation!
+Welcome to SWIFT: SPH With Inter-dependent Fine-grained Tasking's documentation!
================================================================================
Want to get started using SWIFT? Check out the on-boarding guide available
diff --git a/examples/EoSTables/get_eos_tables.sh b/examples/EoSTables/get_eos_tables.sh
new file mode 100755
index 0000000000000000000000000000000000000000..979431777fceee421468138aa0180295f45adef2
--- /dev/null
+++ b/examples/EoSTables/get_eos_tables.sh
@@ -0,0 +1,9 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/HM80_HHe.txt
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/HM80_ice.txt
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/HM80_rock.txt
+
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/SESAME_basalt_7530.txt
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/SESAME_iron_2140.txt
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/SESAME_water_7154.txt
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/EoS/SS08_water.txt
diff --git a/examples/GiantImpactUranus/README.md b/examples/GiantImpactUranus/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..c7d8687886ee9b737ecb492d20d73e780233e9df
--- /dev/null
+++ b/examples/GiantImpactUranus/README.md
@@ -0,0 +1,2 @@
+An example planetary simulation of a giant impact onto the young Uranus with
+~10^6 SPH particles, as described in Kegerreis et al. (2018), ApJ, 861, 52.
diff --git a/examples/GiantImpactUranus/configuration.yml b/examples/GiantImpactUranus/configuration.yml
new file mode 100644
index 0000000000000000000000000000000000000000..ccce852862bec6d1eeba2c132457678564979b8a
--- /dev/null
+++ b/examples/GiantImpactUranus/configuration.yml
@@ -0,0 +1,2 @@
+with-hydro: planetary
+with-equation-of-state: planetary
diff --git a/examples/GiantImpactUranus/get_init_cond.sh b/examples/GiantImpactUranus/get_init_cond.sh
new file mode 100755
index 0000000000000000000000000000000000000000..a6fc4a06fc49f68506783c78026ee467411b44a5
--- /dev/null
+++ b/examples/GiantImpactUranus/get_init_cond.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/uranus_1e6.hdf5
diff --git a/examples/GiantImpactUranus/output_list.txt b/examples/GiantImpactUranus/output_list.txt
new file mode 100644
index 0000000000000000000000000000000000000000..0d233fb11ceea02e40c42b8d512e9c1afbe6e835
--- /dev/null
+++ b/examples/GiantImpactUranus/output_list.txt
@@ -0,0 +1,7 @@
+# Time
+4000
+9000
+14000
+20000
+30000
+40000
\ No newline at end of file
diff --git a/examples/GiantImpactUranus/plot_solution.py b/examples/GiantImpactUranus/plot_solution.py
new file mode 100644
index 0000000000000000000000000000000000000000..faeb071487adc04f0ba37d20967f693ed84652ec
--- /dev/null
+++ b/examples/GiantImpactUranus/plot_solution.py
@@ -0,0 +1,144 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2019 Jacob Kegerreis (jacob.kegerreis@durham.ac.uk)
+ #
+ # This program is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU Lesser General Public License as published
+ # by the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ #
+ # This program is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ # GNU General Public License for more details.
+ #
+ # You should have received a copy of the GNU Lesser General Public License
+ # along with this program. If not, see <http://www.gnu.org/licenses/>.
+ #
+ ##############################################################################
+
+# Plot the snapshots from the example giant impact on Uranus, showing the
+# particles in a thin slice near z=0, coloured by their material, similarish
+# (but not identical) to Fig. 2 in Kegerreis et al. (2018).
+
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import numpy as np
+# import swiftsimio as sw
+import h5py
+
+font_size = 20
+params = {
+ 'axes.labelsize' : font_size,
+ 'font.size' : font_size,
+ 'xtick.labelsize' : font_size,
+ 'ytick.labelsize' : font_size,
+ 'text.usetex' : True,
+ 'font.family' : 'serif',
+ }
+matplotlib.rcParams.update(params)
+
+# Snapshot output times
+output_list = [4000, 9000, 14000, 20000, 30000, 40000]
+
+# Material IDs ( = type_id * type_factor + unit_id )
+type_factor = 100
+type_HM80 = 2
+id_body = 10000
+# Name and ID
+Di_mat_id = {
+ 'HM80_HHe' : type_HM80 * type_factor, # Hydrogen-helium atmosphere
+ 'HM80_ice' : type_HM80 * type_factor + 1, # H20-CH4-NH3 ice mix
+ 'HM80_ice_2' : type_HM80 * type_factor + 1 + id_body,
+ 'HM80_rock' : type_HM80 * type_factor + 2, # SiO2-MgO-FeS-FeO rock mix
+ 'HM80_rock_2' : type_HM80 * type_factor + 2 + id_body,
+ }
+# ID and colour
+Di_id_colour = {
+ Di_mat_id['HM80_HHe'] : '#33DDFF',
+ Di_mat_id['HM80_ice'] : 'lightsteelblue',
+ Di_mat_id['HM80_ice_2'] : '#A080D0',
+ Di_mat_id['HM80_rock'] : 'slategrey',
+ Di_mat_id['HM80_rock_2'] : '#706050',
+ }
+
+def get_snapshot_slice(snapshot):
+ """ Load and select the particles to plot. """
+ # Load particle data
+ # data = load("uranus_1e6_%06d.hdf5" % snapshot)
+ # id = data.gas.particle_ids
+ # pos = data.gas.coordinates
+ # mat_id = data.gas.material
+ with h5py.File("uranus_1e6_%06d.hdf5" % snapshot, 'r') as f:
+ id = f['PartType0/ParticleIDs'].value
+ pos = (f['PartType0/Coordinates'].value
+ - 0.5 * f['Header'].attrs['BoxSize'])
+ mat_id = f['PartType0/MaterialID'].value
+
+ # Edit the material ID of particles in the impactor
+ num_in_target = 869104
+ sel_id = np.where(num_in_target < id)[0]
+ mat_id[sel_id] += id_body
+
+ # Select particles in a thin slice around z=0
+ z_min = -0.1
+ z_max = 0.1
+ sel_z = np.where((z_min < pos[:, 2]) & (pos[:, 2] < z_max))[0]
+ pos = pos[sel_z]
+ mat_id = mat_id[sel_z]
+
+ return pos, mat_id
+
+def plot_snapshot_slice(pos, mat_id):
+ """ Plot the particles, coloured by their material. """
+ colour = np.empty(len(pos), dtype=object)
+ for id, c in Di_id_colour.items():
+ sel_c = np.where(mat_id == id)[0]
+ colour[sel_c] = c
+
+ ax.scatter(pos[:, 0], pos[:, 1], c=colour, edgecolors='none', marker='.',
+ s=10, alpha=0.5, zorder=0)
+
+# Set up the figure
+fig = plt.figure(figsize=(12, 8))
+gs = matplotlib.gridspec.GridSpec(2, 3)
+axes = [plt.subplot(gs[i_y, i_x]) for i_y in range(2) for i_x in range(3)]
+
+# Plot each snapshot
+for i_ax, ax in enumerate(axes):
+ plt.sca(ax)
+ ax.set_rasterization_zorder(1)
+
+ # Load and select the particles to plot
+ pos, mat_id = get_snapshot_slice(output_list[i_ax])
+
+ # Plot the particles, coloured by their material
+ plot_snapshot_slice(pos, mat_id)
+
+ # Axes etc.
+ ax.set_aspect('equal')
+ ax.set_facecolor('k')
+
+ ax.set_xlim(-13, 13)
+ ax.set_ylim(-13, 13)
+
+ if i_ax in [0, 3]:
+ ax.set_ylabel(r"y Postion $(R_\oplus)$")
+ else:
+ ax.set_yticklabels([])
+ if 2 < i_ax:
+ ax.set_xlabel(r"x Postion $(R_\oplus)$")
+ else:
+ ax.set_xticklabels([])
+
+ # Corner time labels
+ x = ax.get_xlim()[0] + 0.04 * (ax.get_xlim()[1] - ax.get_xlim()[0])
+ y = ax.get_ylim()[0] + 0.89 * (ax.get_ylim()[1] - ax.get_ylim()[0])
+ ax.text(x, y, "%.1f h" % (output_list[i_ax] / 60**2), color='w')
+
+plt.subplots_adjust(wspace=0, hspace=0)
+plt.tight_layout()
+
+# Save
+plt.savefig("uranus_1e6.pdf", dpi=200)
\ No newline at end of file
diff --git a/examples/GiantImpactUranus/run.sh b/examples/GiantImpactUranus/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..8c61d42f4e75d73595d387bed32eb7ce34a7d2a3
--- /dev/null
+++ b/examples/GiantImpactUranus/run.sh
@@ -0,0 +1,23 @@
+#!/bin/bash
+
+# Get the initial conditions if they are not present.
+if [ ! -e uranus_1e6.hdf5 ]
+then
+ echo "Fetching initial conditions file for the Uranus impact example..."
+ ./get_init_cond.sh
+fi
+
+# Get the EoS tables if they are not present.
+cd ../EoSTables
+if [ ! -e HM80_HHe.txt ] || [ ! -e HM80_ice.txt ] || [ ! -e HM80_rock.txt ]
+then
+ echo "Fetching equations of state tables for the Uranus impact example..."
+ ./get_eos_tables.sh
+fi
+cd ../GiantImpactUranus
+
+# Run SWIFT
+../swift -s -G -t 8 uranus_1e6.yml 2>&1 | tee output.log
+
+# Plot the solution
+python3 plot_solution.py
diff --git a/examples/GiantImpactUranus/system.yml b/examples/GiantImpactUranus/system.yml
new file mode 100644
index 0000000000000000000000000000000000000000..c99fc7158854ec538e68c44ff74bbb0ba1adfb48
--- /dev/null
+++ b/examples/GiantImpactUranus/system.yml
@@ -0,0 +1,9 @@
+input:
+ - uranus_1e6.yml
+ - output_list.txt
+
+output:
+ - uranus_1e6.pdf
+
+swift_parameters: -s -G
+swift_threads: 8
diff --git a/examples/GiantImpactUranus/uranus_1e6.yml b/examples/GiantImpactUranus/uranus_1e6.yml
new file mode 100644
index 0000000000000000000000000000000000000000..355748d847097623f171078c2ca8372e06a06efa
--- /dev/null
+++ b/examples/GiantImpactUranus/uranus_1e6.yml
@@ -0,0 +1,63 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 5.9724e27 # Grams
+ UnitLength_in_cgs: 6.371e8 # Centimeters
+ UnitVelocity_in_cgs: 6.371e8 # Centimeters per second
+ UnitCurrent_in_cgs: 1 # Amperes
+ UnitTemp_in_cgs: 1 # Kelvin
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: uranus_1e6.hdf5 # The initial conditions file to read
+ periodic: 0 # Are we running with periodic ICs?
+
+# Parameters governing the time integration
+TimeIntegration:
+ time_begin: 0 # The starting time of the simulation (in internal units).
+ time_end: 40000 # 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_max: 100 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: uranus_1e6 # Common part of the name of output files
+ time_first: 0 # Time of the first output (in internal units)
+ delta_time: 1000 # Time difference between consecutive outputs (in internal units)
+ int_time_label_on: 1 # Enable to label the snapshots using the time rounded to an integer (in internal units)
+ output_list_on: 1 # Enable the output list
+ output_list: output_list.txt # File containing the output times (see documentation in "Parameter File" section)
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ time_first: 0 # Time of the first output (in internal units)
+ delta_time: 1000 # Time between statistics output
+
+# Parameters controlling restarts
+Restarts:
+ enable: 0 # Whether to enable dumping restarts at fixed intervals.
+
+# Parameters for the hydrodynamics scheme
+SPH:
+ resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+ delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
+ CFL_condition: 0.2 # Courant-Friedrich-Levy condition for time integration.
+ h_max: 0.5 # Maximal allowed smoothing length (in internal units).
+ viscosity_alpha: 1.5 # Override for the initial value of the artificial viscosity.
+
+# Parameters for the self-gravity scheme
+Gravity:
+ eta: 0.025 # Constant dimensionless multiplier for time integration.
+ theta: 0.7 # Opening angle (Multipole acceptance criterion)
+ comoving_softening: 0.003 # Comoving softening length (in internal units).
+ max_physical_softening: 0.003 # Physical softening length (in internal units).
+
+# Parameters for the task scheduling
+Scheduler:
+ max_top_level_cells: 64 # Maximal number of top-level cells in any dimension. The nu
+
+# Parameters related to the equation of state
+EoS:
+ planetary_use_HM80: 1 # Whether to initialise the Hubbard & MacFarlane (1980) EOS
+ planetary_HM80_HHe_table_file: ../EoSTables/HM80_HHe.txt
+ planetary_HM80_ice_table_file: ../EoSTables/HM80_ice.txt
+ planetary_HM80_rock_table_file: ../EoSTables/HM80_rock.txt
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index eae7d537c831d8a61058ed06d4b2632263a2036e..ca0d47b5c2cc8c05168fa1bd16ac62e64d3068ce 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -184,13 +184,13 @@ EoS:
planetary_use_ANEOS: 0 # (Optional) Whether to prepare the ANEOS EOS
planetary_use_SESAME: 0 # (Optional) Whether to prepare the SESAME EOS
# (Optional) Table file paths
- planetary_HM80_HHe_table_file: ./equation_of_state/planetary_HM80_HHe.txt
- planetary_HM80_ice_table_file: ./equation_of_state/planetary_HM80_ice.txt
- planetary_HM80_rock_table_file: ./equation_of_state/planetary_HM80_rock.txt
- planetary_SESAME_iron_table_file: ./equation_of_state/planetary_SESAME_iron_2140.txt
- planetary_SESAME_basalt_table_file: ./equation_of_state/planetary_SESAME_basalt_7530.txt
- planetary_SESAME_water_table_file: ./equation_of_state/planetary_SESAME_water_7154.txt
- planetary_SS08_water_table_file: ./equation_of_state/planetary_SS08_water.txt
+ planetary_HM80_HHe_table_file: ./EoSTables/planetary_HM80_HHe.txt
+ planetary_HM80_ice_table_file: ./EoSTables/planetary_HM80_ice.txt
+ planetary_HM80_rock_table_file: ./EoSTables/planetary_HM80_rock.txt
+ planetary_SESAME_iron_table_file: ./EoSTables/planetary_SESAME_iron_2140.txt
+ planetary_SESAME_basalt_table_file: ./EoSTables/planetary_SESAME_basalt_7530.txt
+ planetary_SESAME_water_table_file: ./EoSTables/planetary_SESAME_water_7154.txt
+ planetary_SS08_water_table_file: ./EoSTables/planetary_SS08_water.txt
# Parameters related to external potentials --------------------------------------------