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/examples/GiantImpactUranus/README.md b/examples/GiantImpactUranus/README.md new file mode 100644 index 0000000000000000000000000000000000000000..fb5c67bda631dbde57485190a31471e268c04e40 --- /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. \ No newline at end of file diff --git a/examples/GiantImpactUranus/configuration.yml b/examples/GiantImpactUranus/configuration.yml new file mode 100644 index 0000000000000000000000000000000000000000..2cc6b16f6bdb8ebf5c6c1d459a1c3e2d7465052a --- /dev/null +++ b/examples/GiantImpactUranus/configuration.yml @@ -0,0 +1,2 @@ +with-hydro: planetary +with-equation-of-state: planetary \ No newline at end of file diff --git a/examples/GiantImpactUranus/get_eos_tables.sh b/examples/GiantImpactUranus/get_eos_tables.sh new file mode 100755 index 0000000000000000000000000000000000000000..8299127340755d8047e0207ca96cf7904cc29823 --- /dev/null +++ b/examples/GiantImpactUranus/get_eos_tables.sh @@ -0,0 +1,4 @@ +#!/bin/bash +wget http://astro.dur.ac.uk/~cklv53/swift_planet_examples/HM80_HHe.txt +wget http://astro.dur.ac.uk/~cklv53/swift_planet_examples/HM80_ice.txt +wget http://astro.dur.ac.uk/~cklv53/swift_planet_examples/HM80_rock.txt diff --git a/examples/GiantImpactUranus/get_init_cond.sh b/examples/GiantImpactUranus/get_init_cond.sh new file mode 100755 index 0000000000000000000000000000000000000000..8db49f9a7ef89d735f76897cdc360034ac36a67f --- /dev/null +++ b/examples/GiantImpactUranus/get_init_cond.sh @@ -0,0 +1,2 @@ +#!/bin/bash +wget http://astro.dur.ac.uk/~cklv53/swift_planet_examples/uranus_1e6.hdf5 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..0b9c75f118fb801ae46896411efc45abc900c35f --- /dev/null +++ b/examples/GiantImpactUranus/run.sh @@ -0,0 +1,21 @@ +#!/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. +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 + +# 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..bba472bbc173f2610871c2a9946f7c813f5f77fd --- /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: HM80_HHe.txt # Table file paths + planetary_HM80_ice_table_file: HM80_ice.txt + planetary_HM80_rock_table_file: HM80_rock.txt