Merge branch 'planetary_remix_merge_examples' into 'master'

Add new examples from REMIX paper and remove duplicate Planetary hydro tests

See merge request !2117

examples/Planetary/DemoImpactInitCond/plot_profiles.py

@@ -22,7 +22,7 @@ Note that, for standard SPH hydro schemes, especially at low resolution, the
 standard issues that arise at discontinuities in material and density lead the
 SPH particles to settle at slightly different densities near discontinuties.
 For more info and to explore ways to resolve these issues, check out e.g.
-Sandnes et al. (2024), Ruiz-Bonilla el al. (2022), and Kegerreis et al. (2019).
+Sandnes et al. (2025), Ruiz-Bonilla el al. (2022), and Kegerreis et al. (2019).
 The overall profile of the settled SPH planet should still align with the input.
 """
 

examples/Planetary/EvrardCollapse_3D/README.md

@@ -2,13 +2,20 @@ Evrard Collapse 3D (Planetary)
 ===============
 
 This is a copy of `/examples/HydroTests/EvrardCollapse_3D` for testing the 
-Planetary hydro scheme with the planetary ideal gas equation of state. 
+Planetary and REMIX hydro schemes with the planetary ideal gas equation of state.
 
-The results should be highly similar to the Minimal hydro scheme, though 
-slightly different because of the higher viscosity beta used here. To recover 
-the Minimal scheme behaviour, edit `const_viscosity_beta` from 4 to 3 in 
-`src/hydro/Planetary/hydro_parameters.h`.
+The Planetary hydro scheme results should be highly similar to the Minimal
+hydro scheme, though  slightly different because of the higher viscosity beta
+used here. To recover  the Minimal scheme behaviour, edit `const_viscosity_beta`
+from 4 to 3 in `src/hydro/Planetary/hydro_parameters.h`.
 
-This requires the code to be configured to use the planetary hydrodynamics 
-scheme and equations of state: 
-`--with-hydro=planetary --with-equation-of-state=planetary`
+Note that the default resolution_eta parameter is consistent with the use of a
+Wendland C2 kernel with ~100 neighbours.
+
+Some examples of configuration options with different hydro schemes:
+
+REMIX:
+`--with-hydro=remix --with-equation-of-state=planetary --with-kernel=wendland-C2`
+
+"Traditional" SPH (tSPH):
+`--with-hydro=planetary --with-equation-of-state=planetary --with-kernel=wendland-C2`

examples/Planetary/EvrardCollapse_3D/evrard.yml

@@ -26,9 +26,8 @@ Statistics:
 
 # 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).
+  resolution_eta:        1.487    # Target smoothing length in units of the mean inter-particle separation (1.487 == 100Ngbs with the Wendland C2 kernel).
   CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
-  viscosity_alpha:       0.8      # Override for the initial value of the artificial viscosity. In schemes that have a fixed AV, this remains as alpha throughout the run.
 
 # Parameters for the self-gravity scheme
 Gravity:

examples/Planetary/EvrardCollapse_3D/makeIC.py

@@ -67,6 +67,7 @@ h = ones(numPart) * 2.0 * R / numPart ** (1.0 / 3.0)
 v = zeros((numPart, 3))
 ids = linspace(1, numPart, numPart)
 m = ones(numPart) * M / numPart
+rho = M / (2 * pi * R ** 2 * r)
 u = ones(numPart) * u0
 mat = zeros(numPart)
 
@@ -100,6 +101,7 @@ grp = file.create_group("/PartType0")
 grp.create_dataset("Coordinates", data=pos, dtype="d")
 grp.create_dataset("Velocities", data=v, dtype="f")
 grp.create_dataset("Masses", data=m, dtype="f")
+grp.create_dataset("Density", data=rho, dtype="f")
 grp.create_dataset("SmoothingLength", data=h, dtype="f")
 grp.create_dataset("InternalEnergy", data=u, dtype="f")
 grp.create_dataset("ParticleIDs", data=ids, dtype="L")

examples/Planetary/GreshoVortex_3D/README.md

-Gresho Vortex 3D (Planetary)
-=============
-
-This is a copy of `/examples/HydroTests/GreshoVortex_3D` for testing the 
-Planetary hydro scheme with the planetary ideal gas equation of state. 
-
-The results should be highly similar to the Minimal hydro scheme, though 
-slightly different because of the higher viscosity beta used here. To recover 
-the Minimal scheme behaviour, edit `const_viscosity_beta` from 4 to 3 in 
-`src/hydro/Planetary/hydro_parameters.h`.
-
-This requires the code to be configured to use the planetary hydrodynamics 
-scheme and equations of state: 
-`--with-hydro=planetary --with-equation-of-state=planetary`

examples/Planetary/GreshoVortex_3D/makeIC.py

-################################################################################
-# This file is part of SWIFT.
-# Copyright (c) 2016 Matthieu Schaller (schaller@strw.leidenuniv.nl)
-#               2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
-#
-# 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/>.
-#
-################################################################################
-
-import h5py
-from numpy import *
-
-# Generates a swift IC file for the Gresho-Chan vortex in a periodic box
-
-# Parameters
-gamma = 5.0 / 3.0  # Gas adiabatic index
-rho0 = 1  # Gas density
-P0 = 0.0  # Constant additional pressure (should have no impact on the dynamics)
-fileOutputName = "greshoVortex.hdf5"
-fileGlass = "glassCube_64.hdf5"
-# ---------------------------------------------------
-
-# Get position and smoothing lengths from the glass
-fileInput = h5py.File(fileGlass, "r")
-coords = fileInput["/PartType0/Coordinates"][:, :]
-h = fileInput["/PartType0/SmoothingLength"][:]
-ids = fileInput["/PartType0/ParticleIDs"][:]
-boxSize = fileInput["/Header"].attrs["BoxSize"][0]
-numPart = size(h)
-fileInput.close()
-
-# Now generate the rest
-m = ones(numPart) * rho0 * boxSize ** 3 / numPart
-u = zeros(numPart)
-v = zeros((numPart, 3))
-mat = zeros(numPart)
-
-for i in range(numPart):
-
-    x = coords[i, 0]
-    y = coords[i, 1]
-
-    r2 = (x - boxSize / 2) ** 2 + (y - boxSize / 2) ** 2
-    r = sqrt(r2)
-
-    v_phi = 0.0
-    if r < 0.2:
-        v_phi = 5.0 * r
-    elif r < 0.4:
-        v_phi = 2.0 - 5.0 * r
-    else:
-        v_phi = 0.0
-    v[i, 0] = -v_phi * (y - boxSize / 2) / r
-    v[i, 1] = v_phi * (x - boxSize / 2) / r
-    v[i, 2] = 0.0
-
-    P = P0
-    if r < 0.2:
-        P = P + 5.0 + 12.5 * r2
-    elif r < 0.4:
-        P = P + 9.0 + 12.5 * r2 - 20.0 * r + 4.0 * log(r / 0.2)
-    else:
-        P = P + 3.0 + 4.0 * log(2.0)
-    u[i] = P / ((gamma - 1.0) * rho0)
-
-
-# File
-fileOutput = h5py.File(fileOutputName, "w")
-
-# Header
-grp = fileOutput.create_group("/Header")
-grp.attrs["BoxSize"] = [boxSize, boxSize, boxSize]
-grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
-grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
-grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
-grp.attrs["Time"] = 0.0
-grp.attrs["NumFileOutputsPerSnapshot"] = 1
-grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
-grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
-grp.attrs["Dimension"] = 3
-
-# Units
-grp = fileOutput.create_group("/Units")
-grp.attrs["Unit length in cgs (U_L)"] = 1.0
-grp.attrs["Unit mass in cgs (U_M)"] = 1.0
-grp.attrs["Unit time in cgs (U_t)"] = 1.0
-grp.attrs["Unit current in cgs (U_I)"] = 1.0
-grp.attrs["Unit temperature in cgs (U_T)"] = 1.0
-
-# Particle group
-grp = fileOutput.create_group("/PartType0")
-ds = grp.create_dataset("Coordinates", (numPart, 3), "d")
-ds[()] = coords
-ds = grp.create_dataset("Velocities", (numPart, 3), "f")
-ds[()] = v
-ds = grp.create_dataset("Masses", (numPart, 1), "f")
-ds[()] = m.reshape((numPart, 1))
-ds = grp.create_dataset("SmoothingLength", (numPart, 1), "f")
-ds[()] = h.reshape((numPart, 1))
-ds = grp.create_dataset("InternalEnergy", (numPart, 1), "f")
-ds[()] = u.reshape((numPart, 1))
-ds = grp.create_dataset("ParticleIDs", (numPart, 1), "L")
-ds[()] = ids.reshape((numPart, 1))
-ds = grp.create_dataset("MaterialIDs", (numPart, 1), "i")
-ds[()] = mat.reshape((numPart, 1))
-
-fileOutput.close()

examples/Planetary/GreshoVortex_3D/run.sh

-#!/bin/bash
-
- # Generate the initial conditions if they are not present.
-if [ ! -e glassCube_64.hdf5 ]
-then
-    echo "Fetching initial glass file for the Gresho-Chan vortex example..."
-    ../../HydroTests/GreshoVortex_3D/getGlass.sh
-fi
-if [ ! -e greshoVortex.hdf5 ]
-then
-    echo "Generating initial conditions for the Gresho-Chan vortex example..."
-    python3 makeIC.py
-fi
-
-# Run SWIFT
-../../../swift --hydro --threads=4 gresho.yml 2>&1 | tee output.log
-
-# Plot the solution
-python3 ../../HydroTests/GreshoVortex_3D/plotSolution.py 11

examples/Planetary/JupiterLikePlanet/demo_target_n70.yml

+# Define the system of units to use internally.
+InternalUnitSystem:
+    UnitMass_in_cgs:        1e27        # Sets Earth mass = 5.972
+    UnitLength_in_cgs:      1e8         # Sets Earth radius = 6.371
+    UnitVelocity_in_cgs:    1e8         # Sets time in seconds
+    UnitCurrent_in_cgs:     1           # Amperes
+    UnitTemp_in_cgs:        1           # Kelvin
+
+# Parameters related to the initial conditions
+InitialConditions:
+    file_name:  demo_target_n70.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:       20000               # The end time 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
+Snapshots:
+    subdir:             snapshots       # Sub-directory in which to write the snapshots. Defaults to "" (i.e. the directory where SWIFT is run).
+    basename:           demo_target_n70 # Common part of the name of output files
+    time_first:         0               # Time of the first output (in internal units)
+    delta_time:         2000            # Time difference between consecutive outputs (in internal units)
+
+# 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:             1               # Whether to enable dumping restarts at fixed intervals.
+    save:               1               # Whether to save copies of the previous set of restart files (named .prev)
+    subdir:             restart         # Name of subdirectory for restart files.
+    basename:           demo_target_n70 # Prefix used in naming restart files.
+    delta_hours:        10.0            # Decimal hours between dumps of restart files.
+
+# Parameters for the hydrodynamics scheme
+SPH:
+    resolution_eta:     1.487           # Target smoothing length in units of the mean inter-particle separation (1.487 == 100Ngbs with the Wendland C2 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:              10.             # 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.
+    MAC:                            adaptive    # Choice of mulitpole acceptance criterion: 'adaptive' OR 'geometric'.
+    epsilon_fmm:                    0.001       # Tolerance parameter for the adaptive multipole acceptance criterion.
+    theta_cr:                       0.5         # Opening angle for the purely gemoetric criterion.
+    max_physical_baryon_softening:  0.04        # 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.
+    
+# Parameters related to the equation of state
+EoS:
+    # Select which planetary EoS material(s) to enable for use.
+    planetary_use_CD21_HHe:  1     # Hydrogen--helium (Chabrier and Debras 2021), material ID 307
+    planetary_use_AQUA:      1     # AQUA (Haldemann et al. 2020), material ID 304
+    # Tablulated EoS file paths.
+    planetary_CD21_HHe_table_file: ../EoSTables/CD21_HHe.txt
+    planetary_AQUA_table_file:     ../EoSTables/AQUA_H20.txt

examples/Planetary/JupiterLikePlanet/make_init_cond.py

+################################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2025 Thomas Sandnes (thomas.d.sandnes@durham.ac.uk)
+# 		        2024 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/>.
+################################################################################
+
+"""Create initial conditions for settling, using WoMa. See README.md for more info."""
+
+import numpy as np
+import h5py
+import woma
+
+# Number of particles
+N = 10 ** 7
+N_label = "n%d" % (10 * np.log10(N))
+
+# Earth units
+M_E = 5.9724e24  # kg
+R_E = 6.3710e6  # m
+
+# Set profile inputs
+M_t = 308 * M_E
+target_prof = woma.Planet(
+    name="target",
+    A1_mat_layer=["AQUA", "CD21_HHe"],
+    A1_T_rho_type=["adiabatic", "adiabatic"],
+    M=M_t,
+    A1_M_layer=[10 * M_E, 298 * M_E],
+    P_s=1e5,
+    T_s=165,
+    num_prof=10000,
+)
+
+# Load material tables
+woma.load_eos_tables(np.unique(target_prof.A1_mat_layer))
+
+# Compute profiles
+target_prof.gen_prof_L2_find_R_R1_given_M1_M2(R_min=10.7 * R_E, R_max=11.2 * R_E)
+
+# Save profile data
+target_prof.save("demo_target_profile.hdf5")
+
+# Place particles
+target = woma.ParticlePlanet(target_prof, N, seed=12345)
+
+print()
+print("N_target     = %d" % target.N_particles)
+
+# Save the settling initial conditions
+file_to_SI = woma.Conversions(m=1e24, l=1e6, t=1)
+target.save(
+    "demo_target_%s.hdf5" % N_label,
+    boxsize=30 * R_E,
+    file_to_SI=file_to_SI,
+    do_entropies=True,
+)

examples/Planetary/JupiterLikePlanet/plot_profiles.py

+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2025 Thomas Sandnes (thomas.d.sandnes@durham.ac.uk)
+# 	            2023 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 post-settling planetary profiles of the DemoImpactInitCond simulations.
+
+Note that, for standard SPH hydro schemes, especially at low resolution, the
+standard issues that arise at discontinuities in material and density lead the
+SPH particles to settle at slightly different densities near discontinuties.
+For more info and to explore ways to resolve these issues, check out e.g.
+Sandnes et al. (2025), Ruiz-Bonilla el al. (2022), and Kegerreis et al. (2019).
+The overall profile of the settled SPH planet should still align with the input.
+"""
+
+import os
+import matplotlib
+import matplotlib.pyplot as plt
+import numpy as np
+import h5py
+import woma
+
+# Number of particles
+N = 10 ** 7
+N_label = "n%d" % (10 * np.log10(N))
+
+# Plotting options
+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)
+
+
+def plot_profile_and_particles(profile, A1_r, A1_rho):
+    """Plot the particles."""
+    plt.figure(figsize=(7, 7))
+    ax = plt.gca()
+
+    # Earth units
+    R_E = 6.3710e6  # m
+
+    # Profile
+    ax.plot(profile.A1_r / R_E, profile.A1_rho)
+
+    # Particles
+    ax.scatter(A1_r / R_E, A1_rho, c="k", marker=".", s=1 ** 2)
+
+    ax.set_xlim(0, None)
+    ax.set_xlabel(r"Radial distance ($R_\oplus$)")
+    ax.set_ylabel(r"Density (kg m$^{-3}$)")
+
+    plt.tight_layout()
+
+
+if __name__ == "__main__":
+    # Plot each snapshot
+    for body in ["target"]:
+        # Load profiles
+        profile = woma.Planet(load_file="demo_%s_profile.hdf5" % body)
+
+        # Load the data
+        snapshot_id = 5
+        filename = "snapshots/demo_%s_%s_%04d.hdf5" % (body, N_label, snapshot_id)
+        with h5py.File(filename, "r") as f:
+            # Units from file metadata
+            file_to_SI = woma.Conversions(
+                m=float(f["Units"].attrs["Unit mass in cgs (U_M)"]) * 1e-3,
+                l=float(f["Units"].attrs["Unit length in cgs (U_L)"]) * 1e-2,
+                t=float(f["Units"].attrs["Unit time in cgs (U_t)"]),
+            )
+
+            # Particle data
+            A2_pos = (
+                np.array(f["PartType0/Coordinates"][()])
+                - 0.5 * f["Header"].attrs["BoxSize"]
+            ) * file_to_SI.l
+            A1_r = np.sqrt(np.sum(A2_pos ** 2, axis=1))
+            A1_rho = np.array(f["PartType0/Densities"][()]) * file_to_SI.rho
+
+        # Plot the data
+        plot_profile_and_particles(profile, A1_r, A1_rho)
+
+        # Save the figure
+        save = "demo_%s_%s_%04d_prof.png" % (body, N_label, snapshot_id)
+        plt.savefig(save, dpi=200)
+        plt.close()
+
+        print("\rSaved %s" % save)

examples/Planetary/JupiterLikePlanet/plot_snapshots.py

+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2025 Thomas Sandnes (thomas.d.sandnes@durham.ac.uk)
+# 	            2023 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 particle positions from the DemoImpactInitCond settling simulations."""
+
+import os
+import matplotlib
+import matplotlib.pyplot as plt
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+import numpy as np
+import h5py
+import woma
+
+# Number of particles
+N = 10 ** 7
+N_label = "n%d" % (10 * np.log10(N))
+
+# Plotting options
+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 colours
+Di_mat_colour = {"AQUA": "orangered", "CD21_HHe": "gold"}
+Di_id_colour = {woma.Di_mat_id[mat]: colour for mat, colour in Di_mat_colour.items()}
+
+# Scale point size with resolution
+size = (1 * np.cbrt(10 ** 6 / N)) ** 2
+
+
+def load_snapshot(filename):
+    """Load and convert the particle data to plot."""
+    with h5py.File(filename, "r") as f:
+        # Units from file metadata
+        file_to_SI = woma.Conversions(
+            m=float(f["Units"].attrs["Unit mass in cgs (U_M)"]) * 1e-3,
+            l=float(f["Units"].attrs["Unit length in cgs (U_L)"]) * 1e-2,
+            t=float(f["Units"].attrs["Unit time in cgs (U_t)"]),
+        )
+
+        # Particle data
+        A2_pos = (
+            np.array(f["PartType0/Coordinates"][()])
+            - 0.5 * f["Header"].attrs["BoxSize"]
+        ) * file_to_SI.l
+        A1_u = np.array(f["PartType0/InternalEnergies"][()]) * file_to_SI.u
+
+    # Restrict to z < 0 for plotting
+    A1_sel = np.where(A2_pos[:, 2] < 0)[0]
+    A2_pos = A2_pos[A1_sel]
+    A1_u = A1_u[A1_sel]
+
+    return A2_pos, A1_u
+
+
+def plot_snapshot(A2_pos, A1_u):
+    """Plot the particles, coloured by their internal energy."""
+    plt.figure(figsize=(7, 7))
+    ax = plt.gca()
+    ax.set_aspect("equal")
+    cax = make_axes_locatable(ax).append_axes("right", size="5%", pad=0.05)
+
+    # Earth units
+    R_E = 6.3710e6  # m
+
+    # Plot
+    scat = ax.scatter(
+        A2_pos[:, 0] / R_E,
+        A2_pos[:, 1] / R_E,
+        c=A1_u,
+        edgecolors="none",
+        marker=".",
+        s=size,
+    )
+    cbar = plt.colorbar(scat, cax=cax)
+    cbar.set_label(r"Sp. Int. Energy (J kg$^{-1}$)")
+
+    ax_lim = 15
+    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$ ($R_\oplus$)")
+    ax.set_ylabel(r"$y$ ($R_\oplus$)")
+
+    plt.tight_layout()
+
+
+if __name__ == "__main__":
+    # Plot each snapshot
+    for body in ["target"]:
+        # Load the data
+        snapshot_id = 5
+        A2_pos, A1_u = load_snapshot(
+            "snapshots/demo_%s_%s_%04d.hdf5" % (body, N_label, snapshot_id)
+        )
+
+        # Plot the data
+        plot_snapshot(A2_pos, A1_u)
+
+        # Save the figure
+        save = "demo_%s_%s_%04d.png" % (body, N_label, snapshot_id)
+        plt.savefig(save, dpi=200)
+        plt.close()
+
+        print("\rSaved %s" % save)

examples/Planetary/JupiterLikePlanet/run.sh

+#!/bin/bash
+set -o xtrace
+
+# Resolution
+N_label=n70
+
+# Create the initial particle planets
+python3 make_init_cond.py
+
+# Download the equation of state tables if not already present
+if [ ! -e ../EoSTables/CD21_HHe.txt ]
+then
+    cd ../EoSTables
+    ./get_eos_tables.sh
+    cd -
+fi
+
+# Run SWIFT settling simulations
+../../../swift --hydro --self-gravity --threads=28 demo_target_"$N_label".yml \
+    2>&1 | tee output_"$N_label"_t.txt
+
+# Plot the settled particles
+python3 plot_snapshots.py
+python3 plot_profiles.py

examples/Planetary/KelvinHelmholtz_2D/README.md

-Kelvin Helmholtz 2D (Planetary)
-===================
-
-This is a copy of `/examples/HydroTests/KelvinHelmholtz_2D` for testing the 
-Planetary hydro scheme with the planetary ideal gas equation of state. 
-
-The results should be highly similar to the Minimal hydro scheme, though 
-slightly different because of the higher viscosity beta used here. To recover 
-the Minimal scheme behaviour, edit `const_viscosity_beta` from 4 to 3 in 
-`src/hydro/Planetary/hydro_parameters.h`.
-
-This requires the code to be configured to use the planetary hydrodynamics 
-scheme and equations of state: 
-`--with-hydro=planetary --with-equation-of-state=planetary --with-hydro-dimension=2`

examples/Planetary/KelvinHelmholtz_2D/makeIC.py

-###############################################################################
-# This file is part of SWIFT.
-# Copyright (c) 2016 Matthieu Schaller (schaller@strw.leidenuniv.nl)
-#
-# 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/>.
-#
-##############################################################################
-
-import h5py
-from numpy import *
-import sys
-
-# Generates a swift IC file for the Kelvin-Helmholtz vortex in a periodic box
-
-# Parameters
-L2 = 256  # Particles along one edge in the low-density region
-gamma = 5.0 / 3.0  # Gas adiabatic index
-P1 = 2.5  # Central region pressure
-P2 = 2.5  # Outskirts pressure
-v1 = 0.5  # Central region velocity
-v2 = -0.5  # Outskirts vlocity
-rho1 = 2  # Central density
-rho2 = 1  # Outskirts density
-omega0 = 0.1
-sigma = 0.05 / sqrt(2)
-fileOutputName = "kelvinHelmholtz.hdf5"
-# ---------------------------------------------------
-
-# Start by generating grids of particles at the two densities
-numPart2 = L2 * L2
-L1 = int(sqrt(numPart2 / rho2 * rho1))
-numPart1 = L1 * L1
-
-print("N2 =", numPart2, "N1 =", numPart1)
-print("L2 =", L2, "L1 = ", L1)
-print("rho2 =", rho2, "rho1 =", (float(L1 * L1)) / (float(L2 * L2)))
-
-coords1 = zeros((numPart1, 3))
-coords2 = zeros((numPart2, 3))
-h1 = ones(numPart1) * 1.2348 / L1
-h2 = ones(numPart2) * 1.2348 / L2
-m1 = zeros(numPart1)
-m2 = zeros(numPart2)
-u1 = zeros(numPart1)
-u2 = zeros(numPart2)
-vel1 = zeros((numPart1, 3))
-vel2 = zeros((numPart2, 3))
-
-# Particles in the central region
-for i in range(L1):
-    for j in range(L1):
-
-        index = i * L1 + j
-
-        x = i / float(L1) + 1.0 / (2.0 * L1)
-        y = j / float(L1) + 1.0 / (2.0 * L1)
-
-        coords1[index, 0] = x
-        coords1[index, 1] = y
-        u1[index] = P1 / (rho1 * (gamma - 1.0))
-        vel1[index, 0] = v1
-
-# Particles in the outskirts
-for i in range(L2):
-    for j in range(L2):
-
-        index = i * L2 + j
-
-        x = i / float(L2) + 1.0 / (2.0 * L2)
-        y = j / float(L2) + 1.0 / (2.0 * L2)
-
-        coords2[index, 0] = x
-        coords2[index, 1] = y
-        u2[index] = P2 / (rho2 * (gamma - 1.0))
-        vel2[index, 0] = v2
-
-
-# Now concatenate arrays
-where1 = abs(coords1[:, 1] - 0.5) < 0.25
-where2 = abs(coords2[:, 1] - 0.5) > 0.25
-
-coords = append(coords1[where1, :], coords2[where2, :], axis=0)
-
-# print L2*(L2/2), L1*(L1/2)
-# print shape(coords), shape(coords1[where1,:]), shape(coords2[where2,:])
-# print shape(coords), shape(logical_not(coords1[where1,:])), shape(logical_not(coords2[where2,:]))
-
-vel = append(vel1[where1, :], vel2[where2, :], axis=0)
-h = append(h1[where1], h2[where2], axis=0)
-m = append(m1[where1], m2[where2], axis=0)
-u = append(u1[where1], u2[where2], axis=0)
-numPart = size(h)
-ids = linspace(1, numPart, numPart)
-mat = zeros(numPart)
-m[:] = (0.5 * rho1 + 0.5 * rho2) / float(numPart)
-
-# Velocity perturbation
-vel[:, 1] = (
-    omega0
-    * sin(4 * pi * coords[:, 0])
-    * (
-        exp(-((coords[:, 1] - 0.25) ** 2) / (2 * sigma ** 2))
-        + exp(-((coords[:, 1] - 0.75) ** 2) / (2 * sigma ** 2))
-    )
-)
-
-# File
-fileOutput = h5py.File(fileOutputName, "w")
-
-# Header
-grp = fileOutput.create_group("/Header")
-grp.attrs["BoxSize"] = [1.0, 1.0, 0.1]
-grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
-grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
-grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
-grp.attrs["Time"] = 0.0
-grp.attrs["NumFileOutputsPerSnapshot"] = 1
-grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
-grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
-grp.attrs["Dimension"] = 2
-
-# Units
-grp = fileOutput.create_group("/Units")
-grp.attrs["Unit length in cgs (U_L)"] = 1.0
-grp.attrs["Unit mass in cgs (U_M)"] = 1.0
-grp.attrs["Unit time in cgs (U_t)"] = 1.0
-grp.attrs["Unit current in cgs (U_I)"] = 1.0
-grp.attrs["Unit temperature in cgs (U_T)"] = 1.0
-
-# Particle group
-grp = fileOutput.create_group("/PartType0")
-ds = grp.create_dataset("Coordinates", (numPart, 3), "d")
-ds[()] = coords
-ds = grp.create_dataset("Velocities", (numPart, 3), "f")
-ds[()] = vel
-ds = grp.create_dataset("Masses", (numPart, 1), "f")
-ds[()] = m.reshape((numPart, 1))
-ds = grp.create_dataset("SmoothingLength", (numPart, 1), "f")
-ds[()] = h.reshape((numPart, 1))
-ds = grp.create_dataset("InternalEnergy", (numPart, 1), "f")
-ds[()] = u.reshape((numPart, 1))
-ds = grp.create_dataset("ParticleIDs", (numPart, 1), "L")
-ds[()] = ids.reshape((numPart, 1))
-ds = grp.create_dataset("MaterialIDs", (numPart, 1), "i")
-ds[()] = mat.reshape((numPart, 1))
-
-fileOutput.close()

examples/Planetary/KelvinHelmholtz_2D/run.sh

-#!/bin/bash
-
- # Generate the initial conditions if they are not present.
-if [ ! -e kelvinHelmholtz.hdf5 ]
-then
-    echo "Generating initial conditions for the Kelvin-Helmholtz example..."
-    python3 makeIC.py
-fi
-
-# Run SWIFT
-../../../swift --hydro --threads=4 kelvinHelmholtz.yml 2>&1 | tee output.log
-
-
-# Plot the solution
-python3 ../../HydroTests/KelvinHelmholtz_2D/makeMovieSwiftsimIO.py

examples/Planetary/KelvinHelmholtz_EarthLike_3D/README.md

+Kelvin--Helmholtz Instabilty (Earth-like, equal mass, 3D)
+--------------
+
+This is a 3D version of the Kelvin--Helmholtz instability. These initial
+conditions are those used to produce the simulations presented by
+Sandnes et al. (2025), section 4.4.
+
+This test uses particles of equal mass and has sharp density and velocity
+discontinuities. Equations of state and conditions are representative of those
+within Earth's interior.
+
+Note that the default resolution_eta parameter is consistent with the use of a
+Wendland C2 kernel with ~100 neighbours.
+
+Some examples of configuration options with different hydro schemes:
+
+REMIX:
+`--with-hydro=remix --with-equation-of-state=planetary --with-kernel=wendland-C2`
+
+"Traditional" SPH (tSPH):
+`--with-hydro=planetary --with-equation-of-state=planetary --with-kernel=wendland-C2`

examples/Planetary/KelvinHelmholtz_EarthLike_3D/kelvin_helmholtz.yml

+# 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 governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   10520 # The end time of the simulation (in internal units). Corresponding to 2 \tau_{KH}.
+  dt_min:     1e-9  # The minimal time-step size of the simulation (in internal units).
+  dt_max:     1e2   # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:            kelvin_helmholtz # Common part of the name of output files
+  time_first:          0.               # Time of the first output (in internal units)
+  delta_time:          526              # Time difference between consecutive outputs (in internal units). Corresponding to 0.1 \tau_{KH}.
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          526 # Time between statistics output
+
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.487    # Target smoothing length in units of the mean inter-particle separation (1.487 == 100Ngbs with the Wendland C2 kernel).
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+  
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./kelvin_helmholtz.hdf5      # The file to read
+  periodic:   1
+
+Scheduler:
+    max_top_level_cells:     40
+
+# Parameters related to the equation of state
+EoS:
+    # Select which planetary EoS material(s) to enable for use.
+    planetary_use_ANEOS_forsterite:   1     # ANEOS forsterite (Stewart et al. 2019), material ID 400
+    planetary_use_ANEOS_Fe85Si15:     1     # ANEOS Fe85Si15 (Stewart 2020), material ID 402
+    # Tablulated EoS file paths.
+    planetary_ANEOS_forsterite_table_file:  ../EoSTables/ANEOS_forsterite_S19.txt
+    planetary_ANEOS_Fe85Si15_table_file:    ../EoSTables/ANEOS_Fe85Si15_S20.txt

examples/Planetary/KelvinHelmholtz_EarthLike_3D/makeIC.py

+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2025 Thomas Sandnes (thomas.d.sandnes@durham.ac.uk)
+#               2016 Matthieu Schaller (matthieu.schaller@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/>.
+#
+##############################################################################
+
+import h5py
+import numpy as np
+
+# Generates a swift IC file for the Kelvin-Helmholtz test in a periodic box
+
+# Constants
+R_earth = 6371000  # Earth radius
+
+# Parameters
+N2_l = 128  # Particles along one edge in the low-density region
+N2_depth = 18  # Particles in z direction in low-density region
+matID1 = 402  # Central region material ID: ANEOS Fe85Si15
+matID2 = 400  # Outskirts material ID: ANEOS forsterite
+P1 = 1.2e11  # Central region pressure
+P2 = 1.2e11  # Outskirts pressure
+u1 = 4069874  # Central region specific internal energy
+u2 = 9899952  # Outskirts specific internal energy
+rho1_approx = 10000  # Central region density. Readjusted later
+rho2 = 5000  # Outskirts density
+boxsize_l = R_earth  # size of simulation box in x and y dimension
+v1 = boxsize_l / 10000  # Central region velocity
+v2 = -boxsize_l / 10000  # Outskirts velocity
+boxsize_depth = boxsize_l * N2_depth / N2_l  # size of simulation box in z dimension
+mass = rho2 * (boxsize_l * boxsize_l * boxsize_depth) / (N2_l * N2_l * N2_depth)
+fileOutputName = "kelvin_helmholtz.hdf5"
+# ---------------------------------------------------
+
+# Start by calculating N1_l and rho1
+numPart2 = N2_l * N2_l * N2_depth
+numPart1_approx = int(numPart2 / rho2 * rho1_approx)
+
+# Consider numPart1 = N1_l * N1_l * N1_depth
+# Substituting boxsize_depth / boxsize_l = N1_depth / N1_l gives,
+# numPart1 = N1_l * N1_l * (boxsize_depth / boxsize_l) * N1_l, which ranges to:
+N1_l = int(np.cbrt(numPart1_approx * boxsize_l / boxsize_depth))
+# Make sure this is a multiple of 4 since this is the number of KH vortices
+N1_l -= N1_l % 4
+
+N1_depth = int(boxsize_depth * N1_l / boxsize_l)
+numPart1 = int(N1_l * N1_l * N1_depth)
+
+# The density of the central region can then be calculated
+rho1 = mass * (N1_l * N1_l * N1_depth) / (boxsize_l * boxsize_l * boxsize_depth)
+
+# Now construct two lattices of particles in the two regions
+A2_coords1 = np.empty((numPart1, 3))
+A2_coords2 = np.empty((numPart2, 3))
+A2_vel1 = np.zeros((numPart1, 3))
+A2_vel2 = np.zeros((numPart2, 3))
+A2_vel1[:, 0] = v1
+A2_vel2[:, 0] = v2
+
+A1_mat1 = np.full(numPart1, matID1)
+A1_mat2 = np.full(numPart2, matID2)
+A1_m1 = np.full(numPart1, mass)
+A1_m2 = np.full(numPart2, mass)
+A1_rho1 = np.full(numPart1, rho1)
+A1_rho2 = np.full(numPart2, rho2)
+A1_u1 = np.full(numPart1, u1)
+A1_u2 = np.full(numPart2, u2)
+A1_h1 = np.full(numPart1, boxsize_l / N1_l)
+A1_h2 = np.full(numPart2, boxsize_l / N2_l)
+
+# Particles in the central region
+for i in range(N1_depth):
+    for j in range(N1_l):
+        for k in range(N1_l):
+            index = i * N1_l ** 2 + j * N1_l + k
+            A2_coords1[index, 0] = (j / float(N1_l) + 1.0 / (2.0 * N1_l)) * boxsize_l
+            A2_coords1[index, 1] = (k / float(N1_l) + 1.0 / (2.0 * N1_l)) * boxsize_l
+            A2_coords1[index, 2] = (
+                i / float(N1_depth) + 1.0 / (2.0 * N1_depth)
+            ) * boxsize_depth
+
+# Particles in the outskirts
+for i in range(N2_depth):
+    for j in range(N2_l):
+        for k in range(N2_l):
+            index = i * N2_l ** 2 + j * N2_l + k
+            A2_coords2[index, 0] = (j / float(N2_l) + 1.0 / (2.0 * N2_l)) * boxsize_l
+            A2_coords2[index, 1] = (k / float(N2_l) + 1.0 / (2.0 * N2_l)) * boxsize_l
+            A2_coords2[index, 2] = (
+                i / float(N2_depth) + 1.0 / (2.0 * N2_depth)
+            ) * boxsize_depth
+
+
+# Masks for the particles to be selected for the outer and inner regions
+mask1 = abs(A2_coords1[:, 1] - 0.5 * boxsize_l) < 0.25 * boxsize_l
+mask2 = abs(A2_coords2[:, 1] - 0.5 * boxsize_l) > 0.25 * boxsize_l
+
+# The positions of the particles are now selected
+# and the placement of the lattices are adjusted to give appropriate interfaces
+A2_coords_inside = A2_coords1[mask1, :]
+A2_coords_outside = A2_coords2[mask2, :]
+
+# Calculate the separation of particles across the density discontinuity
+pcl_separation_1 = np.cbrt(mass / rho1)
+pcl_separation_2 = np.cbrt(mass / rho2)
+boundary_separation = 0.5 * (pcl_separation_1 + pcl_separation_2)
+
+# Shift all the "inside" particles to get boundary_separation across the bottom interface
+min_y_inside = np.min(A2_coords_inside[:, 1])
+max_y_outside_bot = np.max(
+    A2_coords_outside[A2_coords_outside[:, 1] - 0.5 * boxsize_l < -0.25 * boxsize_l, 1]
+)
+shift_distance_bot = boundary_separation - (min_y_inside - max_y_outside_bot)
+A2_coords_inside[:, 1] += shift_distance_bot
+
+# Shift the top section of the "outside" particles to get boundary_separation across the top interface
+max_y_inside = np.max(A2_coords_inside[:, 1])
+min_y_outside_top = np.min(
+    A2_coords_outside[A2_coords_outside[:, 1] - 0.5 * boxsize_l > 0.25 * boxsize_l, 1]
+)
+shift_distance_top = boundary_separation - (min_y_outside_top - max_y_inside)
+A2_coords_outside[
+    A2_coords_outside[:, 1] - 0.5 * boxsize_l > 0.25 * boxsize_l, 1
+] += shift_distance_top
+
+# Adjust box size in y direction based on the shifting of the lattices.
+new_box_y = boxsize_l + shift_distance_top
+
+# Now the two lattices can be combined
+A2_coords = np.append(A2_coords_inside, A2_coords_outside, axis=0)
+A2_vel = np.append(A2_vel1[mask1], A2_vel2[mask2], axis=0)
+A1_mat = np.append(A1_mat1[mask1], A1_mat2[mask2], axis=0)
+A1_m = np.append(A1_m1[mask1], A1_m2[mask2], axis=0)
+A1_rho = np.append(A1_rho1[mask1], A1_rho2[mask2], axis=0)
+A1_u = np.append(A1_u1[mask1], A1_u2[mask2], axis=0)
+A1_h = np.append(A1_h1[mask1], A1_h2[mask2], axis=0)
+numPart = np.size(A1_m)
+A1_ids = np.linspace(1, numPart, numPart)
+
+# Finally add the velocity perturbation
+vel_perturb_factor = 0.01 * (v1 - v2)
+A2_vel[:, 1] = vel_perturb_factor * np.sin(
+    2 * np.pi * A2_coords[:, 0] / (0.5 * boxsize_l)
+)
+
+# Write ICs to file
+with h5py.File(fileOutputName, "w") as f:
+    # Header
+    grp = f.create_group("/Header")
+    grp.attrs["BoxSize"] = [boxsize_l, new_box_y, boxsize_depth]
+    grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
+    grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+    grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
+    grp.attrs["Time"] = 0.0
+    grp.attrs["NumFileOutputsPerSnapshot"] = 1
+    grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+    grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
+    grp.attrs["Dimension"] = 3
+
+    # Units
+    grp = f.create_group("/Units")
+    grp.attrs["Unit length in cgs (U_L)"] = 100.0
+    grp.attrs["Unit mass in cgs (U_M)"] = 1000.0
+    grp.attrs["Unit time in cgs (U_t)"] = 1.0
+    grp.attrs["Unit current in cgs (U_I)"] = 1.0
+    grp.attrs["Unit temperature in cgs (U_T)"] = 1.0
+
+    # Particle group
+    grp = f.create_group("/PartType0")
+    ds = grp.create_dataset("Coordinates", (numPart, 3), "d")
+    ds[()] = A2_coords
+    ds = grp.create_dataset("Velocities", (numPart, 3), "f")
+    ds[()] = A2_vel
+    ds = grp.create_dataset("Masses", (numPart, 1), "f")
+    ds[()] = A1_m.reshape((numPart, 1))
+    ds = grp.create_dataset("Density", (numPart, 1), "f")
+    ds[()] = A1_rho.reshape((numPart, 1))
+    ds = grp.create_dataset("SmoothingLength", (numPart, 1), "f")
+    ds[()] = A1_h.reshape((numPart, 1))
+    ds = grp.create_dataset("InternalEnergy", (numPart, 1), "f")
+    ds[()] = A1_u.reshape((numPart, 1))
+    ds = grp.create_dataset("ParticleIDs", (numPart, 1), "L")
+    ds[()] = A1_ids.reshape((numPart, 1))
+    ds = grp.create_dataset("MaterialIDs", (numPart, 1), "i")
+    ds[()] = A1_mat.reshape((numPart, 1))

examples/Planetary/KelvinHelmholtz_EarthLike_3D/plotModeGrowth.py

+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2025 Thomas Sandnes (thomas.d.sandnes@durham.ac.uk)
+#               2025 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/>.
+#
+##############################################################################
+"""
+Generate plot Mode growth of the 3D Kelvin--Helmholtz instability.
+This is based on the quantity calculated in Eqns. 10--13 of McNally et al. 2012
+"""
+
+import sys
+import h5py
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+import numpy as np
+
+
+def calculate_mode_growth(snap, vy_init_amp, wavelength):
+
+    # Load data
+    snap_file = "kelvin_helmholtz_%04d.hdf5" % snap
+
+    with h5py.File(snap_file, "r") as f:
+        boxsize_l = f["Header"].attrs["BoxSize"][0]
+        A2_pos = f["/PartType0/Coordinates"][:, :] / boxsize_l
+        A2_vel = f["/PartType0/Velocities"][:, :] / (vy_init_amp * boxsize_l)
+        A1_h = f["/PartType0/SmoothingLengths"][:] / boxsize_l
+
+    # Adjust positions
+    A2_pos[:, 0] += 0.5
+    A2_pos[:, 1] += 0.5
+
+    # Masks to select the upper and lower halfs of the simulations
+    mask_up = A2_pos[:, 1] >= 0.5
+    mask_down = A2_pos[:, 1] < 0.5
+
+    # McNally et al. 2012 Eqn. 10
+    s = np.empty(len(A1_h))
+    s[mask_down] = (
+        A2_vel[mask_down, 1]
+        * A1_h[mask_down] ** 3
+        * np.sin(2 * np.pi * A2_pos[mask_down, 0] / wavelength)
+        * np.exp(-2 * np.pi * np.abs(A2_pos[mask_down, 1] - 0.25) / wavelength)
+    )
+    s[mask_up] = (
+        A2_vel[mask_up, 1]
+        * A1_h[mask_up] ** 3
+        * np.sin(2 * np.pi * A2_pos[mask_up, 0] / wavelength)
+        * np.exp(-2 * np.pi * np.abs((1 - A2_pos[mask_up, 1]) - 0.25) / wavelength)
+    )
+
+    # McNally et al. 2012 Eqn. 11
+    c = np.empty(len(A1_h))
+    c[mask_down] = (
+        A2_vel[mask_down, 1]
+        * A1_h[mask_down] ** 3
+        * np.cos(2 * np.pi * A2_pos[mask_down, 0] / wavelength)
+        * np.exp(-2 * np.pi * np.abs(A2_pos[mask_down, 1] - 0.25) / wavelength)
+    )
+    c[mask_up] = (
+        A2_vel[mask_up, 1]
+        * A1_h[mask_up] ** 3
+        * np.cos(2 * np.pi * A2_pos[mask_up, 0] / wavelength)
+        * np.exp(-2 * np.pi * np.abs((1 - A2_pos[mask_up, 1]) - 0.25) / wavelength)
+    )
+
+    # McNally et al. 2012 Eqn. 12
+    d = np.empty(len(A1_h))
+    d[mask_down] = A1_h[mask_down] ** 3 * np.exp(
+        -2 * np.pi * np.abs(A2_pos[mask_down, 1] - 0.25) / wavelength
+    )
+    d[mask_up] = A1_h[mask_up] ** 3 * np.exp(
+        -2 * np.pi * np.abs((1 - A2_pos[mask_up, 1]) - 0.25) / wavelength
+    )
+
+    # McNally et al. 2012 Eqn. 13
+    M = 2 * np.sqrt((np.sum(s) / np.sum(d)) ** 2 + (np.sum(c) / np.sum(d)) ** 2)
+
+    return M
+
+
+if __name__ == "__main__":
+
+    # Simulation paramerters for nomralisation of mode growth
+    vy_init_amp = (
+        2e-6
+    )  # Initial amplitude of y velocity perturbation in units of the boxsize per second
+    wavelength = 0.5  # wavelength of initial perturbation in units of the boxsize
+
+    # Use Gridspec to set up figure
+    n_gs_ax = 40
+    ax_len = 5
+    fig = plt.figure(figsize=(9, 6))
+    gs = mpl.gridspec.GridSpec(nrows=40, ncols=60)
+    ax = plt.subplot(gs[:, :])
+
+    # Snapshots and corresponding times
+    snaps = np.arange(21)
+    times = 0.1 * snaps
+
+    # Calculate mode mode growth
+    mode_growth = np.empty(len(snaps))
+    for i, snap in enumerate(snaps):
+        M = calculate_mode_growth(snap, vy_init_amp, wavelength)
+        mode_growth[i] = M
+
+    # Plot
+    ax.plot(times, np.log10(mode_growth), linewidth=1.5)
+
+    ax.set_ylabel(r"$\log( \, M \; / \; M^{}_0 \, )$", fontsize=18)
+    ax.set_xlabel(r"$t \; / \; \tau^{}_{\rm KH}$", fontsize=18)
+    ax.minorticks_on()
+    ax.tick_params(which="major", direction="in")
+    ax.tick_params(which="minor", direction="in")
+    ax.tick_params(labelsize=14)
+
+    plt.savefig("mode_growth.png", dpi=300, bbox_inches="tight")

examples/Planetary/KelvinHelmholtz_EarthLike_3D/plotSnapshots.py

+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2025 Thomas Sandnes (thomas.d.sandnes@durham.ac.uk)
+#               2025 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/>.
+#
+##############################################################################
+"""
+Generate plot of the 3D Kelvin--Helmholtz instability.
+"""
+
+import sys
+import h5py
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+import numpy as np
+
+
+def make_axes():
+    # Use Gridspec to set up figure
+    n_gs_ax = 40
+    n_gs_ax_gap = 1
+    n_gs_cbar_gap = 1
+    n_gs_cbar = 2
+
+    ax_len = 5
+    ax_gap_len = n_gs_ax_gap * ax_len / n_gs_ax
+    cbar_gap_len = n_gs_cbar_gap * ax_len / n_gs_ax
+    cbar_len = n_gs_cbar * ax_len / n_gs_ax
+
+    fig = plt.figure(
+        figsize=(3 * ax_len + 2 * ax_gap_len + cbar_gap_len + cbar_len, ax_len)
+    )
+    gs = mpl.gridspec.GridSpec(
+        nrows=n_gs_ax, ncols=3 * n_gs_ax + 2 * n_gs_ax_gap + n_gs_cbar_gap + n_gs_cbar
+    )
+
+    ax_0 = plt.subplot(gs[:n_gs_ax, :n_gs_ax])
+    ax_1 = plt.subplot(gs[:n_gs_ax, n_gs_ax + n_gs_ax_gap : 2 * n_gs_ax + n_gs_ax_gap])
+    ax_2 = plt.subplot(
+        gs[:n_gs_ax, 2 * n_gs_ax + 2 * n_gs_ax_gap : 3 * n_gs_ax + 2 * n_gs_ax_gap]
+    )
+
+    cax_0 = plt.subplot(
+        gs[
+            : int(0.5 * n_gs_ax) - 1,
+            3 * n_gs_ax
+            + 2 * n_gs_ax_gap
+            + n_gs_cbar_gap : 3 * n_gs_ax
+            + 2 * n_gs_ax_gap
+            + n_gs_cbar_gap
+            + n_gs_cbar,
+        ]
+    )
+    cax_1 = plt.subplot(
+        gs[
+            int(0.5 * n_gs_ax) + 1 :,
+            3 * n_gs_ax
+            + 2 * n_gs_ax_gap
+            + n_gs_cbar_gap : 3 * n_gs_ax
+            + 2 * n_gs_ax_gap
+            + n_gs_cbar_gap
+            + n_gs_cbar,
+        ]
+    )
+
+    axs = [ax_0, ax_1, ax_2]
+    caxs = [cax_0, cax_1]
+
+    return axs, caxs
+
+
+def plot_kh(ax, snap, mat_id1, mat_id2, cmap1, cmap2, norm1, norm2):
+
+    # Load data
+    snap_file = "kelvin_helmholtz_%04d.hdf5" % snap
+
+    with h5py.File(snap_file, "r") as f:
+        # Units from file metadata to SI
+        m = float(f["Units"].attrs["Unit mass in cgs (U_M)"][0]) * 1e-3
+        l = float(f["Units"].attrs["Unit length in cgs (U_L)"][0]) * 1e-2
+
+        boxsize_l = f["Header"].attrs["BoxSize"][0] * l
+        A1_x = f["/PartType0/Coordinates"][:, 0] * l
+        A1_y = f["/PartType0/Coordinates"][:, 1] * l
+        A1_z = f["/PartType0/Coordinates"][:, 2] * l
+        A1_rho = f["/PartType0/Densities"][:] * (m / l ** 3)
+        A1_m = f["/PartType0/Masses"][:] * m
+        A1_mat_id = f["/PartType0/MaterialIDs"][:]
+
+    # Sort arrays based on z position
+    sort_indices = np.argsort(A1_z)
+    A1_x = A1_x[sort_indices]
+    A1_y = A1_y[sort_indices]
+    A1_z = A1_z[sort_indices]
+    A1_rho = A1_rho[sort_indices]
+    A1_m = A1_m[sort_indices]
+    A1_mat_id = A1_mat_id[sort_indices]
+
+    # Mask to select slice
+    slice_thickness = 0.2 * (np.max(A1_z) - np.min(A1_z))
+    slice_pos_z = 0.5 * (np.max(A1_z) + np.min(A1_z))
+    mask_slice = np.logical_and(
+        A1_z > slice_pos_z - 0.5 * slice_thickness,
+        A1_z < slice_pos_z + 0.5 * slice_thickness,
+    )
+
+    # Select particles to plot
+    A1_x_slice = A1_x[mask_slice]
+    A1_y_slice = A1_y[mask_slice]
+    A1_rho_slice = A1_rho[mask_slice]
+    A1_m_slice = A1_m[mask_slice]
+    A1_mat_id_slice = A1_mat_id[mask_slice]
+
+    # Size of plotted particles
+    size_factor = 5e4
+    A1_size = size_factor * (A1_m_slice / A1_rho_slice) ** (2 / 3) / boxsize_l ** 2
+
+    mask_mat1 = A1_mat_id_slice == mat_id1
+    mask_mat2 = A1_mat_id_slice == mat_id2
+
+    # Plot
+    scatter = ax.scatter(
+        A1_x_slice[mask_mat1],
+        A1_y_slice[mask_mat1],
+        c=A1_rho_slice[mask_mat1],
+        norm=norm1,
+        cmap=cmap1,
+        s=A1_size[mask_mat1],
+        edgecolors="none",
+    )
+
+    scatter = ax.scatter(
+        A1_x_slice[mask_mat2],
+        A1_y_slice[mask_mat2],
+        c=A1_rho_slice[mask_mat2],
+        norm=norm2,
+        cmap=cmap2,
+        s=A1_size[mask_mat2],
+        edgecolors="none",
+    )
+
+    ax.set_xticks([])
+    ax.set_yticks([])
+    ax.set_facecolor((0.9, 0.9, 0.9))
+    ax.set_xlim((0, boxsize_l))
+    ax.set_ylim((0, boxsize_l))
+
+
+if __name__ == "__main__":
+
+    # Set colormap
+    cmap1 = plt.get_cmap("YlOrRd")
+    mat_id1 = 402
+    rho_min1 = 7950
+    rho_max1 = 10050
+    norm1 = mpl.colors.Normalize(vmin=rho_min1, vmax=rho_max1)
+
+    cmap2 = plt.get_cmap("Blues_r")
+    mat_id2 = 400
+    rho_min2 = 4950
+    rho_max2 = 5550
+    norm2 = mpl.colors.Normalize(vmin=rho_min2, vmax=rho_max2)
+
+    # Generate axes
+    axs, caxs = make_axes()
+
+    # The three snapshots to be plotted
+    snaps = [10, 15, 20]
+    times = ["1.0", "1.5", "2.0"]
+
+    # Plot
+    for i, snap in enumerate(snaps):
+        ax = axs[i]
+        time = times[i]
+
+        plot_kh(ax, snap, mat_id1, mat_id2, cmap1, cmap2, norm1, norm2)
+        ax.text(
+            0.5,
+            -0.1,
+            r"$t =\;$" + time + r"$\, \tau^{}_{\rm KH}$",
+            horizontalalignment="center",
+            size=18,
+            transform=ax.transAxes,
+        )
+
+    # Colour bar
+    sm1 = plt.cm.ScalarMappable(cmap=cmap1, norm=norm1)
+    cbar1 = plt.colorbar(sm1, caxs[0])
+    cbar1.ax.tick_params(labelsize=14)
+    cbar1.set_label(r"Iron density (kg/m$^3$)", rotation=90, labelpad=8, fontsize=12)
+
+    sm2 = plt.cm.ScalarMappable(cmap=cmap2, norm=norm2)
+    cbar2 = plt.colorbar(sm2, caxs[1])
+    cbar2.ax.tick_params(labelsize=14)
+    cbar2.set_label(r"Rock density (kg/m$^3$)", rotation=90, labelpad=16, fontsize=12)
+
+    plt.savefig("kelvin_helmholtz.png", dpi=300, bbox_inches="tight")