diff --git a/examples/GiantImpacts/GiantImpact/README.md b/examples/GiantImpacts/GiantImpact/README.md
deleted file mode 100644
index 5c72f1bb37717e4e08b90cdd8c2ccc2da8fe3c59..0000000000000000000000000000000000000000
--- a/examples/GiantImpacts/GiantImpact/README.md
+++ /dev/null
@@ -1,6 +0,0 @@
-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.
-
-This example requires the code to be configured to use the Planetary
-hydrodynamics and equation of state:
- --with-hydro=planetary --with-equation-of-state=planetary
diff --git a/examples/GiantImpacts/GiantImpact/get_init_cond.sh b/examples/GiantImpacts/GiantImpact/get_init_cond.sh
deleted file mode 100755
index e81035ead7e50ef204bb6bb476e94c7fe0eae0f6..0000000000000000000000000000000000000000
--- a/examples/GiantImpacts/GiantImpact/get_init_cond.sh
+++ /dev/null
@@ -1,2 +0,0 @@
-#!/bin/bash
-wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/GiantImpacts/uranus_1e6.hdf5
diff --git a/examples/GiantImpacts/GiantImpact/output_list.txt b/examples/GiantImpacts/GiantImpact/output_list.txt
deleted file mode 100644
index 0d233fb11ceea02e40c42b8d512e9c1afbe6e835..0000000000000000000000000000000000000000
--- a/examples/GiantImpacts/GiantImpact/output_list.txt
+++ /dev/null
@@ -1,7 +0,0 @@
-# Time
-4000
-9000
-14000
-20000
-30000
-40000
\ No newline at end of file
diff --git a/examples/GiantImpacts/GiantImpact/plot_solution.py b/examples/GiantImpacts/GiantImpact/plot_solution.py
deleted file mode 100644
index faeb071487adc04f0ba37d20967f693ed84652ec..0000000000000000000000000000000000000000
--- a/examples/GiantImpacts/GiantImpact/plot_solution.py
+++ /dev/null
@@ -1,144 +0,0 @@
-###############################################################################
- # 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/GiantImpacts/GiantImpact/run.sh b/examples/GiantImpacts/GiantImpact/run.sh
deleted file mode 100755
index 5a4f0a74dd098b1fff4659a7d72be3845ad47fc6..0000000000000000000000000000000000000000
--- a/examples/GiantImpacts/GiantImpact/run.sh
+++ /dev/null
@@ -1,23 +0,0 @@
-#!/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 ../GiantImpact
-
-# 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/GiantImpacts/GiantImpact/system.yml b/examples/GiantImpacts/GiantImpact/system.yml
deleted file mode 100644
index c99fc7158854ec538e68c44ff74bbb0ba1adfb48..0000000000000000000000000000000000000000
--- a/examples/GiantImpacts/GiantImpact/system.yml
+++ /dev/null
@@ -1,9 +0,0 @@
-input:
- - uranus_1e6.yml
- - output_list.txt
-
-output:
- - uranus_1e6.pdf
-
-swift_parameters: -s -G
-swift_threads: 8
diff --git a/examples/Planetary/EarthImpact/README.md b/examples/Planetary/EarthImpact/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..b49289f105332fcc621b2b2e95058e4e332d2bb8
--- /dev/null
+++ b/examples/Planetary/EarthImpact/README.md
@@ -0,0 +1,12 @@
+An example planetary simulation of the first 10 hours of a grazing giant impact
+onto the proto-Earth (a roughly canonical Moon-forming impact) with the simple
+Tillotson equation of state and 10^5 SPH particles.
+
+See Kegerreis et al. (2019), Mon. Not. R. Astron. Soc., 487:4 for more about
+using SWIFT for planetary simulations.
+
+More examples and planetary-specific documentation coming soon!
+
+This example requires the code to be configured to use the planetary
+hydrodynamics scheme and equations of state:
+ --with-hydro=planetary --with-equation-of-state=planetary
diff --git a/examples/GiantImpacts/GiantImpact/configuration.yml b/examples/Planetary/EarthImpact/configuration.yml
similarity index 100%
rename from examples/GiantImpacts/GiantImpact/configuration.yml
rename to examples/Planetary/EarthImpact/configuration.yml
diff --git a/examples/GiantImpacts/GiantImpact/uranus_1e6.yml b/examples/Planetary/EarthImpact/earth_impact.yml
similarity index 75%
rename from examples/GiantImpacts/GiantImpact/uranus_1e6.yml
rename to examples/Planetary/EarthImpact/earth_impact.yml
index b98467616f659babb488352c87697a43fedcc0db..b1df0a6974e060074b66b2ab9c3021aa02d65190 100644
--- a/examples/GiantImpacts/GiantImpact/uranus_1e6.yml
+++ b/examples/Planetary/EarthImpact/earth_impact.yml
@@ -8,24 +8,22 @@ InternalUnitSystem:
# Parameters related to the initial conditions
InitialConditions:
- file_name: uranus_1e6.hdf5 # The initial conditions file to read
+ file_name: earth_impact.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).
+ 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_max: 100 # The maximal 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
Snapshots:
- basename: uranus_1e6 # Common part of the name of output files
+ basename: earth_impact # 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:
@@ -41,7 +39,7 @@ 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).
+ h_max: 0.2 # 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
@@ -53,11 +51,8 @@ Gravity:
# Parameters for the task scheduling
Scheduler:
- max_top_level_cells: 64 # Maximal number of top-level cells in any dimension. The nu
+ max_top_level_cells: 16 # 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
+ planetary_use_Til: 1 # Whether to initialise the Tillotson EOS
\ No newline at end of file
diff --git a/examples/Planetary/EarthImpact/get_init_cond.sh b/examples/Planetary/EarthImpact/get_init_cond.sh
new file mode 100755
index 0000000000000000000000000000000000000000..c8428c528afd61f7d25f0f0ccf4305eced63e429
--- /dev/null
+++ b/examples/Planetary/EarthImpact/get_init_cond.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/Planetary/earth_impact.hdf5
diff --git a/examples/Planetary/EarthImpact/make_anim.sh b/examples/Planetary/EarthImpact/make_anim.sh
new file mode 100755
index 0000000000000000000000000000000000000000..6331c0d1c9da6830e7dc649a21503c1fef2b7663
--- /dev/null
+++ b/examples/Planetary/EarthImpact/make_anim.sh
@@ -0,0 +1,7 @@
+#!/bin/bash
+
+# Make a simple animation of the snapshots
+out="earth_impact.mp4"
+ffmpeg -framerate 5 -i earth_impact_%?%?%?%?%?%?.png $out -y
+
+echo Saved $out
\ No newline at end of file
diff --git a/examples/Planetary/EarthImpact/plot_solution.py b/examples/Planetary/EarthImpact/plot_solution.py
new file mode 100644
index 0000000000000000000000000000000000000000..faf6612bada9f8d5f11c3c1f93843d3d8224ae90
--- /dev/null
+++ b/examples/Planetary/EarthImpact/plot_solution.py
@@ -0,0 +1,147 @@
+###############################################################################
+# 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 the proto-Earth, showing
+# the particles in a thin slice near z=0, coloured by their material.
+
+import matplotlib
+import matplotlib.pyplot as plt
+import numpy as np
+import swiftsimio as sw
+import unyt
+
+font_size = 20
+params = {
+ "axes.labelsize": font_size,
+ "font.size": font_size,
+ "xtick.labelsize": font_size,
+ "ytick.labelsize": font_size,
+ "font.family": "serif",
+}
+matplotlib.rcParams.update(params)
+
+# Material IDs ( = type_id * type_factor + unit_id )
+type_factor = 100
+type_Til = 1
+id_body = 10000
+# Name and ID
+Di_mat_id = {
+ "Til_iron": type_Til * type_factor,
+ "Til_iron_2": type_Til * type_factor + id_body,
+ "Til_granite": type_Til * type_factor + 1,
+ "Til_granite_2": type_Til * type_factor + 1 + id_body,
+}
+# Colour
+Di_mat_colour = {
+ "Til_iron": "darkgray",
+ "Til_granite": "orangered",
+ "Til_iron_2": "saddlebrown",
+ "Til_granite_2": "gold",
+}
+Di_id_colour = {Di_mat_id[mat]: colour for mat, colour in Di_mat_colour.items()}
+
+
+def load_snapshot(snapshot_time, ax_lim):
+ """ Select and load the particles to plot. """
+ # Snapshot to load
+ snapshot = "earth_impact_%06d.hdf5" % snapshot_time
+
+ # Only load data with the axis limits and below z=0
+ ax_lim = 0.1
+ mask = sw.mask(snapshot)
+ box_mid = 0.5 * mask.metadata.boxsize[0].to(unyt.Rearth)
+ x_min = box_mid - ax_lim * unyt.Rearth
+ x_max = box_mid + ax_lim * unyt.Rearth
+ load_region = [[x_min, x_max], [x_min, x_max], [x_min, box_mid]]
+ mask.constrain_spatial(load_region)
+
+ # Load
+ data = sw.load(snapshot, mask=mask)
+ pos = data.gas.coordinates.to(unyt.Rearth) - box_mid
+ id = data.gas.particle_ids
+ mat_id = data.gas.material_ids.value
+
+ # Restrict to z < 0
+ sel = np.where(pos[:, 2] < 0)[0]
+ pos = pos[sel]
+ id = id[sel]
+ mat_id = mat_id[sel]
+
+ # Sort in z order so higher particles are plotted on top
+ sort = np.argsort(pos[:, 2])
+ pos = pos[sort]
+ id = id[sort]
+ mat_id = mat_id[sort]
+
+ # Edit material IDs for particles in the impactor
+ num_in_target = 99740
+ mat_id[num_in_target <= id] += id_body
+
+ return pos, mat_id
+
+
+def plot_snapshot(pos, mat_id, ax_lim):
+ """ Plot the particles, coloured by their material. """
+ plt.figure(figsize=(7, 7))
+ ax = plt.gca()
+ ax.set_aspect("equal")
+
+ colour = np.empty(len(pos), dtype=object)
+ for id_c, c in Di_id_colour.items():
+ colour[mat_id == id_c] = c
+
+ ax.scatter(
+ pos[:, 0],
+ pos[:, 1],
+ c=colour,
+ edgecolors="none",
+ marker=".",
+ s=10,
+ alpha=0.5,
+ )
+
+ ax.set_xlim(-ax_lim, ax_lim)
+ ax.set_yticks(ax.get_xticks())
+ ax.set_ylim(-ax_lim, ax_lim)
+ ax.set_xlabel(r"x Position ($R_\oplus$)")
+ ax.set_ylabel(r"y Position ($R_\oplus$)")
+
+ plt.tight_layout()
+
+
+if __name__ == "__main__":
+ print()
+ # Axis limits (Earth radii)
+ ax_lim = 3.4
+
+ # Plot each snapshot
+ for snapshot_time in range(0, 36000 + 1, 1000):
+ # Load the data
+ pos, mat_id = load_snapshot(snapshot_time, ax_lim)
+
+ # Plot the data
+ plot_snapshot(pos, mat_id, ax_lim)
+
+ # Save the figure
+ save = "earth_impact_%06d.png" % snapshot_time
+ plt.savefig(save, dpi=100)
+
+ print("\rSaved %s" % save)
+
+ print()
diff --git a/examples/Planetary/EarthImpact/run.sh b/examples/Planetary/EarthImpact/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..40ddcbb63a74751fc9a623e17a0d9bbe94360cb2
--- /dev/null
+++ b/examples/Planetary/EarthImpact/run.sh
@@ -0,0 +1,17 @@
+#!/bin/bash
+
+# Get the initial conditions if they are not present.
+if [ ! -e earth_impact.hdf5 ]
+then
+ echo "Fetching initial conditions file for the Earth impact example..."
+ ./get_init_cond.sh
+fi
+
+# Run SWIFT
+../../swift -s -G -t 8 earth_impact.yml 2>&1 | tee output.log
+
+# Plot the snapshots
+python3 plot_solution.py
+
+# Make a simple animation
+./make_anim.sh
diff --git a/examples/Planetary/EarthImpact/system.yml b/examples/Planetary/EarthImpact/system.yml
new file mode 100644
index 0000000000000000000000000000000000000000..8e942aa1fcc99da28f15d212fc1792cfcd99c417
--- /dev/null
+++ b/examples/Planetary/EarthImpact/system.yml
@@ -0,0 +1,8 @@
+input:
+ - earth_impact.yml
+
+output:
+ - earth_impact.pdf
+
+swift_parameters: -s -G
+swift_threads: 8
diff --git a/examples/GiantImpacts/EoSTables/get_eos_tables.sh b/examples/Planetary/EoSTables/get_eos_tables.sh
similarity index 100%
rename from examples/GiantImpacts/EoSTables/get_eos_tables.sh
rename to examples/Planetary/EoSTables/get_eos_tables.sh
diff --git a/src/hydro/Planetary/hydro_io.h b/src/hydro/Planetary/hydro_io.h
index 7829fd6376f603d4786489504773de35bf0dd1b5..b0599ceba9101be21c843838c600e2e971cb3e69 100644
--- a/src/hydro/Planetary/hydro_io.h
+++ b/src/hydro/Planetary/hydro_io.h
@@ -61,15 +61,15 @@ INLINE static void hydro_read_particles(struct part* parts,
UNIT_CONV_SPEED, parts, v);
list[2] = io_make_input_field("Masses", FLOAT, 1, COMPULSORY, UNIT_CONV_MASS,
parts, mass);
- list[3] = io_make_input_field("SmoothingLength", FLOAT, 1, COMPULSORY,
+ list[3] = io_make_input_field("SmoothingLengths", FLOAT, 1, COMPULSORY,
UNIT_CONV_LENGTH, parts, h);
- list[4] = io_make_input_field("InternalEnergy", FLOAT, 1, COMPULSORY,
+ list[4] = io_make_input_field("InternalEnergies", FLOAT, 1, COMPULSORY,
UNIT_CONV_ENERGY_PER_UNIT_MASS, parts, u);
list[5] = io_make_input_field("ParticleIDs", ULONGLONG, 1, COMPULSORY,
UNIT_CONV_NO_UNITS, parts, id);
list[6] = io_make_input_field("Accelerations", FLOAT, 3, OPTIONAL,
UNIT_CONV_ACCELERATION, parts, a_hydro);
- list[7] = io_make_input_field("Density", FLOAT, 1, OPTIONAL,
+ list[7] = io_make_input_field("Densities", FLOAT, 1, OPTIONAL,
UNIT_CONV_DENSITY, parts, rho);
list[8] = io_make_input_field("MaterialIDs", INT, 1, COMPULSORY,
UNIT_CONV_NO_UNITS, parts, mat_id);