Commit 48d8f0fe authored by Matthieu Schaller's avatar Matthieu Schaller

Merge branch 'logger_example' into 'master'

Logger example See merge request !1019
parents 453b1ef9 8d4d8fc4
......@@ -55,7 +55,7 @@ SPH:
# Parameters related to the initial conditions
InitialConditions:
file_name: ./zoom_in.hdf5 # The file to read
file_name: ./h177.hdf5 # The file to read
periodic: 1
cleanup_h_factors: 1 # Remove the h-factors inherited from Gadget
cleanup_velocity_factors: 1 # Remove the sqrt(a) factor in the velocities inherited from Gadget
......
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2020 Loic Hausammann (loic.hausammann@epfl.ch)
#
# 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 numpy as np
from astropy import units
from swiftsimio import Writer
import unyt
np.random.seed(50)
# parameters
filename = "simple_orbits.hdf5"
num_part = 1
masses = 1.
# If changed, need to update simple_orbits.yml
M = units.solMass.to("earthMass")
M = float(M)
min_r = 0.2
max_r = 5
boxsize = 2 * max_r
u_l = 1.49e13 # AU
u_m = 5.97e27 # Earth mass
u_v = 474047. # AU / yr
u_t = u_l / u_v
G = 1.189972e-04 # grav. const.
# generate the coordinates
dist = np.random.rand(num_part) * (max_r - min_r) + min_r
angle = np.random.rand(num_part) * 2 * np.pi
# more easy to do it in 2D, therefore coords[:, 2] == 0
coords = np.zeros((num_part, 3))
coords[:, 0] = dist * np.sin(angle)
coords[:, 1] = dist * np.cos(angle)
coords += boxsize * 0.5
# generate the masses
m = np.ones(num_part) * masses
# generate the velocities
sign = np.random.rand(num_part)
sign[sign < 0.5] = -1
sign[sign >= 0.5] = 1
v = np.zeros((num_part, 3))
v[:, 0] = sign * np.sqrt(G * M / (dist * (1 + np.tan(angle)**2)))
v[:, 1] = - np.tan(angle) * v[:, 0]
# Write the snapshot
units = unyt.UnitSystem("Planets", unyt.AU, unyt.mearth, unyt.yr)
snapshot = Writer(units, boxsize * unyt.AU)
snapshot.dark_matter.coordinates = coords * unyt.AU
snapshot.dark_matter.velocities = v * unyt.AU / unyt.yr
snapshot.dark_matter.masses = m * unyt.mearth
snapshot.write(filename)
# Line Properties
lines.linewidth: 2.
# Font options
font.size: 10
font.family: STIXGeneral
mathtext.fontset: stix
# LaTeX options
text.usetex: True
# Legend options
legend.frameon: True
legend.fontsize: 10
# Figure options for publications
figure.dpi: 300
# Histogram options
hist.bins: auto
# Ticks inside plots; more space devoted to plot.
xtick.direction: in
ytick.direction: in
# Always plot ticks on top of data
axes.axisbelow: False
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2020 Loic Hausammann (loic.hausammann@epfl.ch)
#
# 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/>.
#
##############################################################################
from h5py import File as HDF5File
import numpy as np
from glob import glob
import matplotlib.pyplot as plt
import makeIC
import sys
sys.path.append("../../../logger/.libs/")
import liblogger as logger
# Plot parameters
plt.style.use("mpl_style")
center = np.array([0.5 * makeIC.boxsize]*3)
id_focus = 0
# Defines the constants
G = 1.189972e-04 # gravitational constant
M = 3.329460e+05 # mass of the sun
# generate the figures
fig_1 = plt.figure()
fig_2 = plt.figure()
fig_3 = plt.figure()
def gravity(t, y):
"""
Compute the equation of motion.
Parameters
----------
t: float
Time of the step
y: np.array
Variable to evolve [x, y, vx, vy].
"""
dy = np.zeros(4)
dy[:2] = y[2:]
r = np.sum((y[:2] - center[:2])**2)**0.5
a = - G * M / r**3
dy[2] = y[2] * a
dy[3] = y[3] * a
return dy
def plotRelative(t, p, *args, **kwargs):
"""
Wrapper around the function plot from matplotlib, but
plot the relative evolution of the variable.
"""
p = (p - p[0]) / p[0]
plt.plot(t, p, *args, **kwargs)
def doSnapshots():
"""
Read the snapshots and plot the corresponding variables.
"""
# get all the filenames
filenames = glob("simple_orbits_*.hdf5")
N = len(filenames)
filenames.sort()
# generate the output arrays
E = np.zeros((N, makeIC.num_part))
t = np.zeros(N)
p = np.zeros((N, 3))
v = np.zeros((N, 3))
for i, f in enumerate(filenames):
# get the data from the file
f = HDF5File(f, "r")
ids = f["PartType1/ParticleIDs"][:]
sort = np.argsort(ids)
ids = ids[sort]
pos = f["PartType1/Coordinates"][sort, :]
pos -= center
vel = f["PartType1/Velocities"][sort, :]
t[i] = f["Header"].attrs["Time"]
r = np.sum(pos**2, axis=1)**0.5
v2 = np.sum(vel**2, axis=1)
E[i, :] = 0.5 * v2 - G * M / r
# Get the pos / vel of the required particle
ind = ids == id_focus
p[i, :] = pos[ind, :]
v[i, :] = vel[ind, :]
# Compute the solution
y0 = np.zeros(4)
y0[:2] = p[0, :2]
y0[2:] = v[0, :2]
# compute the plotting variables
plt.figure(fig_1.number)
plotRelative(t, E, ".", label="Snapshot")
plt.figure(fig_2.number)
plt.plot(p[:, 0], p[:, 1], "-", label="Snapshot", lw=1.)
plt.figure(fig_3.number)
plt.plot(v[:, 0], v[:, 1], "-", label="Snapshot", lw=1.)
def doStatistics():
"""
Do the plots with the energy output file.
"""
data = np.genfromtxt("energy.txt", names=True)
times = data["Time"]
E = data["E_tot"]
plt.figure(fig_1.number)
plotRelative(times, E, "-", label="Statistics")
def doLogger():
"""
Read the logfile and plot the corresponding variables.
"""
basename = "index"
N = 1000
verbose = 0
# Get time limits
t_min, t_max = logger.getTimeLimits(basename, verbose)
times = np.linspace(t_min, t_max, N)
# Create output arrays
E = np.zeros((N, makeIC.num_part))
E_parts = np.zeros((N, makeIC.num_part))
p = np.zeros((N, 3))
v = np.zeros((N, 3))
t_parts = np.zeros((N, makeIC.num_part))
p_parts = np.zeros((N, 3))
v_parts = np.zeros((N, 3))
# Read the particles
parts = logger.loadSnapshotAtTime(
basename, times[0], verbose)
for i, t in enumerate(times):
# Get the next particles
interp = logger.moveForwardInTime(
basename, parts, t, verbose)
ids = interp["ids"]
sort = np.argsort(ids)
ids = ids[sort]
rel_pos = interp["positions"][sort, :] - center
vel = interp["velocities"][sort, :]
rel_pos_parts = parts["positions"][sort, :] - center
vel_parts = parts["velocities"][sort, :]
# Compute the interpolated variables
r = np.sum(rel_pos**2, axis=1)**0.5
v2 = np.sum(vel**2, axis=1)
E[i, :] = 0.5 * v2 - G * M / r
ind = ids == id_focus
p[i, :] = rel_pos[ind, :]
v[i, :] = vel[ind, :]
# Compute the variables of the last record
r = np.sum(rel_pos_parts**2, axis=1)**0.5
v2 = np.sum(vel_parts**2, axis=1)
E_parts[i, :] = 0.5 * v2 - G * M / r
t_parts[i, :] = parts["times"][sort]
ind = ids == id_focus
p_parts[i, :] = rel_pos_parts[ind, :]
v_parts[i, :] = vel_parts[ind, :]
# compute the plotting variables
plt.figure(fig_1.number)
plotRelative(t_parts, E_parts, "x", label="Logger")
plotRelative(times, E, "--", label="Logger (Interpolation)")
# Compute the solution
y0 = np.zeros(4)
y0[:2] = p[0, :2]
y0[2:] = v[0, :2]
t_parts, ind = np.unique(t_parts[:, 0], return_index=True)
# plot the solution
plt.figure(fig_2.number)
plt.plot(p_parts[:, 0], p_parts[:, 1], "x", label="Logger")
plt.figure(fig_3.number)
plt.plot(v_parts[:, 0], v_parts[:, 1], "x", label="Logger")
# Compute the solution
y0 = np.zeros(4)
y0[:2] = p[0, :2]
y0[2:] = v[0, :2]
plt.figure(fig_2.number)
plt.plot(p[:, 0], p[:, 1], ":r", label="Logger (Interpolation)")
plt.figure(fig_3.number)
plt.plot(v[:, 0], v[:, 1], ":r", label="Logger (Interpolation)")
# do all the plots
doStatistics()
doSnapshots()
doLogger()
# add text
plt.figure(fig_1.number)
plt.xlabel("Time [yr]")
plt.ylabel(r"$\frac{E - E(t=0)}{E(t=0)}$")
plt.legend(ncol=2)
plt.savefig("Energy.png")
plt.figure(fig_2.number)
plt.xlabel("Position [AU]")
plt.ylabel("Position [AU]")
plt.axis("equal")
plt.legend()
plt.savefig("Positions.pdf")
plt.figure(fig_3.number)
plt.xlabel("Velocity [AU / yr]")
plt.ylabel("Velocity [AU / yr]")
plt.axis("equal")
plt.legend()
plt.savefig("Velocities.png")
plt.show()
#!/bin/bash
# Generate the initial conditions.
echo "Generating initial conditions for the Simple Orbits example..."
python makeIC.py
# Run SWIFT
../../swift --external-gravity --threads=1 simple_orbits.yml 2>&1 | tee output.log
# Plot the solution
python3 plotSolution.py
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 5.97e27 # M_earth
UnitLength_in_cgs: 1.49e13 # AU
UnitVelocity_in_cgs: 474047. # AU / yr
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: 4.2 # The end time of the simulation (in internal units).
dt_min: 1e-12 # The minimal time-step size of the simulation (in internal units).
dt_max: 1 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the snapshots
Snapshots:
basename: simple_orbits # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 5e-3 # Time difference between consecutive outputs (in internal units)
# Parameters governing the conserved quantities statistics
Statistics:
delta_time: 1 # Time between statistics output
# Parameters related to the initial conditions
InitialConditions:
file_name: simple_orbits.hdf5 # The file to read
periodic: 0
# Point mass external potentials
PointMassPotential:
useabspos: 0 # 0 -> positions based on centre, 1 -> absolute positions
position: [0.,0,0.] # location of external point mass (internal units)
mass: 332946 # mass of external point mass (solar mass in internal units)
timestep_mult: 1e-3 # Dimensionless pre-factor for the time-step condition
# Parameters governing the logger snapshot system
Logger:
delta_step: 4000 # Update the particle log every this many updates
basename: index # Common part of the filenames
initial_buffer_size: 0.01 # (Optional) Buffer size in GB
buffer_scale: 10 # (Optional) When buffer size is too small, update it with required memory times buffer_scale
index_mem_frac: 1e-6
......@@ -54,4 +54,11 @@ Scheduler:
# Parameters related to the equation of state
EoS:
planetary_use_Til: 1 # Whether to initialise the Tillotson EOS
\ No newline at end of file
planetary_use_Til: 1 # Whether to initialise the Tillotson EOS
# Parameters governing the logger snapshot system
Logger:
delta_step: 10 # Update the particle log every this many updates
basename: index # Common part of the filenames
initial_buffer_size: 0.01 # (Optional) Buffer size in GB
buffer_scale: 10 # (Optional) When buffer size is too small, update it with required memory times buffer_scale
#!/usr/bin/env python3
import numpy as np
import matplotlib.pyplot as plt
import sys
import swiftsimio.visualisation.projection as vis
from copy import deepcopy
from shutil import copyfile
import os
from subprocess import call
sys.path.append("../../../logger/.libs/")
import liblogger as logger
boxsize = 80.
large_width = 15
small_width = 0.4
alpha_particle = 0
width = 0
basename = "index"
resolution = 1080
resolution_mini = resolution * 100 // 512
id_foc = 99839
v_max = 8.
traj = None
def selectParticles(pos, width):
ind1 = np.logical_and(pos[:, 0] > 0, pos[:, 0] < width)
ind2 = np.logical_and(pos[:, 1] > 0, pos[:, 1] < width)
ind3 = np.logical_and(ind1, ind2)
# avoid the zoom
size = resolution_mini * width / resolution
box = np.logical_and(pos[:, 0] > (width - size), pos[:, 1] < size)
ind3 = np.logical_and(~box, ind3)
return ind3
def getImage(parts, width, center, res):
pos = parts["positions"]
pos = pos - center + 0.5 * width
pos /= width
h = parts["smoothing_lengths"] / width
m = np.ones(pos.shape[0])
# Do the projection
img = vis.scatter(pos[:, 0], pos[:, 1], m, h, res).T
img /= width**3
ind = img == 0
img[ind] = img[~ind].min()
img = np.log10(img)
return img
def doProjection(parts, t, skip):
plt.figure(figsize=(8, 8))
global traj, id_foc
# evolve in time the particles
interp = logger.moveForwardInTime(basename, parts, t)
# Check if some particles where removed
ind = parts["smoothing_lengths"] == 0
ind = np.arange(parts.shape[0])[ind]
if len(ind) != 0:
parts = np.delete(parts, ind)
interp = np.delete(interp, ind)
id_foc -= np.sum(ind < id_foc)
# Get the position and the image center
pos = interp["positions"]
position = pos[id_foc, :]
# save the trajectory
if traj is None:
traj = position[np.newaxis, :]
else:
traj = np.append(traj, position[np.newaxis, :], axis=0)
# Do we want to generate the image?
if skip:
return parts
# Compute the images
img = getImage(interp, width, position, resolution)
img[:resolution_mini, -resolution_mini:] = getImage(
interp, 0.2 * boxsize, position, resolution_mini)
box = width * np.array([0, 1, 0, 1])
plt.imshow(img, origin="lower", extent=box, cmap="plasma",
vmin=0, vmax=v_max)
# plot the particles
pos = interp["positions"] - position + 0.5 * width
ind = selectParticles(pos, width)
ms = 2
plt.plot(pos[ind, 0], pos[ind, 1], "o", markersize=ms,
alpha=alpha_particle, color="silver")
plt.plot(pos[id_foc, 0], pos[id_foc, 1], "or", markersize=2*ms,
alpha=alpha_particle)
# plot time
plt.text(0.5 * width, 0.95 * width,
"Time = %0.2f Hours" % (t / (60 * 60)), color="w",
horizontalalignment="center")
# plot trajectory
tr = deepcopy(traj) - position + 0.5 * width
plt.plot(tr[:, 0], tr[:, 1], "-", color="w", alpha=alpha_particle)
# plot scale
plt.plot([0.05 * width, 0.15 * width], [0.05 * width, 0.05 * width],
"-w")
plt.text(0.1 * width, 0.06 * width, "%.2f R$_\oplus$" % (0.1 * width),
horizontalalignment='center', color="w")
# style
plt.axis("off")
plt.xlim(box[:2])
plt.ylim(box[2:])
plt.tight_layout()
plt.style.use('dark_background')
return parts
def skipImage(i):
if os.path.isfile("output/image_%04i.png" % i):
print("Skipping image %i" % i)
return True
else:
return False
def doMovie(parts, t0, t1, N, init):
times = np.linspace(t0, t1, N)
for i, t in enumerate(times):
print("Image %i / %i" % (i+1, N))
skip = skipImage(i + init)
parts = doProjection(parts, t, skip)
if not skip:
plt.savefig("output/image_%04i.png" % (i + init))
plt.close()
return init + N
def doZoom(parts, w_init, w_end, N, init, t0, t1, increase_alpha):
global width, alpha_particle
widths = np.linspace(w_init, w_end, N)
alpha = np.linspace(-8, 0.2, N)
if not increase_alpha:
alpha = alpha[::-1]
k = 5 # parameter for the steepness
alpha = 1. / (1 + np.exp(- 2 * k * alpha))
times = np.linspace(t0, t1, N)
for i, w in enumerate(widths):
print("Image %i / %i" % (i+1, N))
skip = skipImage(i + init)
width = w
alpha_particle = alpha[i]
parts = doProjection(parts, times[i], skip)
if not skip:
plt.savefig("output/image_%04i.png" % (i + init))
plt.close()
return init + N
def doStatic(init, number, ref):
# copy the first picture
for i in range(number):
copyfile("output/image_%04i.png" % ref,
"output/image_%04i.png" % (i + init))
return init + number
def doTitle(frames):
plt.figure()
plt.axis("off")
box = np.array([0, 1, 0, 1])
plt.xlim(box[:2])
plt.ylim(box[2:])
plt.tight_layout()
plt.style.use('dark_background')
style = {
"verticalalignment": "center",
"horizontalalignment": "center",
"fontweight": "bold"
}
plt.text(0.5, 0.6, "Planetary Impact with the Particle Logger",
**style)
plt.text(0.5, 0.5, "L. Hausammann, J. Kegerreis and P. Gonnet 2020",
**style)
for i in range(frames):
plt.savefig("output/image_%04i.png" % i)
plt.close()
return frames
def main():
t0, t1 = logger.getTimeLimits(basename)
parts = logger.loadSnapshotAtTime(basename, t0)
# Do a few title frames
init = doTitle(40)
after_title = init
frames_after_title = 10
init += frames_after_title
# do a zoom while moving forward in time
init = doZoom(parts, large_width, small_width, 50,
init=init, t0=t0, t1=0.07 * t1,
increase_alpha=True)
# Copy a few frames
doStatic(after_title, frames_after_title,
after_title + frames_after_title)
init = doStatic(init, 10, init-1)
# Do the main part of the movie
init = doMovie(parts, 0.07 * t1, 0.15 * t1, N=250,
init=init)
# copy a few frames
init = doStatic(init, 10, init-1)
# zoom out and finish the movie
init = doZoom(parts, small_width, large_width, 607,
init=init, t0=0.15 * t1, t1=t1,
increase_alpha=False)
# copy a few frames
init = doStatic(init, 10, init-1)
if __name__ == "__main__":
main()
convert = "ffmpeg -i output/image_%04d.png -y -vcodec libx264 "
convert += "-profile:v high444 -refs 16 -crf 0 "
convert += "-preset ultrafast movie.mp4"
call(convert, shell=True)
......@@ -50,7 +50,7 @@ print("time: %g" % time)
# read dump
t = logger.getTimeLimits(basename)
data = logger.loadSnapshotAtTime(basename, time)
data = logger.loadSnapshotAtTime(basename, time, 2)
print("The data contains the following elements:")
print(data.dtype.names)
......
......@@ -184,7 +184,6 @@ void header_read(struct header *h, struct logger_logfile *log) {