Skip to content
Snippets Groups Projects
Commit 65e79c9a authored by Jacob Kegerreis's avatar Jacob Kegerreis Committed by Matthieu Schaller
Browse files

Planetary examples

parent 02704bc5
Branches
Tags
No related merge requests found
Showing
with 203 additions and 197 deletions
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
#!/bin/bash
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/GiantImpacts/uranus_1e6.hdf5
# Time
4000
9000
14000
20000
30000
40000
\ No newline at end of file
###############################################################################
# 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
#!/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
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
......@@ -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
#!/bin/bash
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/Planetary/earth_impact.hdf5
#!/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
###############################################################################
# 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()
#!/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
input:
- uranus_1e6.yml
- output_list.txt
- earth_impact.yml
output:
- uranus_1e6.pdf
- earth_impact.pdf
swift_parameters: -s -G
swift_threads: 8
......@@ -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);
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment