Merge branch 'planetary_examples' into 'master'

Planetary examples

Closes #633

See merge request !976

examples/GiantImpacts/GiantImpact/README.md

-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

examples/GiantImpacts/GiantImpact/get_init_cond.sh

-#!/bin/bash
-wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/GiantImpacts/uranus_1e6.hdf5

examples/GiantImpacts/GiantImpact/output_list.txt

-# Time
-4000
-9000
-14000
-20000
-30000
-40000
\ No newline at end of file

examples/GiantImpacts/GiantImpact/plot_solution.py

-###############################################################################
- # 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

examples/GiantImpacts/GiantImpact/run.sh

-#!/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

examples/Planetary/EarthImpact/README.md

+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

examples/GiantImpacts/GiantImpact/uranus_1e6.yml → 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

examples/Planetary/EarthImpact/get_init_cond.sh

+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/Planetary/earth_impact.hdf5

examples/Planetary/EarthImpact/make_anim.sh

+#!/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

examples/Planetary/EarthImpact/plot_solution.py

+###############################################################################
+# 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()

examples/Planetary/EarthImpact/run.sh

+#!/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

examples/GiantImpacts/GiantImpact/system.yml → examples/Planetary/EarthImpact/system.yml

 input:
-  - uranus_1e6.yml
-  - output_list.txt
+  - earth_impact.yml
 
 output:
-  - uranus_1e6.pdf
+  - earth_impact.pdf
 
 swift_parameters:   -s -G
 swift_threads:      8

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);