diff --git a/examples/GEAR/ZoomIn/zoom_in.yml b/examples/GEAR/ZoomIn/zoom_in.yml
index 1d01e6abd33454c4682bb2adfbfdb76a64e4948a..e0927d344f31b820743c83560893bce60982b1b9 100644
--- a/examples/GEAR/ZoomIn/zoom_in.yml
+++ b/examples/GEAR/ZoomIn/zoom_in.yml
@@ -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
diff --git a/examples/Logger/SimpleOrbits/makeIC.py b/examples/Logger/SimpleOrbits/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..a2418492584d0c79fa6eb900ca2b0f758c4c6084
--- /dev/null
+++ b/examples/Logger/SimpleOrbits/makeIC.py
@@ -0,0 +1,75 @@
+###############################################################################
+# 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)
diff --git a/examples/Logger/SimpleOrbits/mpl_style b/examples/Logger/SimpleOrbits/mpl_style
new file mode 100644
index 0000000000000000000000000000000000000000..9b85038b2490f0079ae09ecc45d24c0050863731
--- /dev/null
+++ b/examples/Logger/SimpleOrbits/mpl_style
@@ -0,0 +1,26 @@
+# 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
diff --git a/examples/Logger/SimpleOrbits/plotSolution.py b/examples/Logger/SimpleOrbits/plotSolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..b160e24964240b7a8a89fc03594aa4fb149b65bc
--- /dev/null
+++ b/examples/Logger/SimpleOrbits/plotSolution.py
@@ -0,0 +1,252 @@
+###############################################################################
+# 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()
diff --git a/examples/Logger/SimpleOrbits/run.sh b/examples/Logger/SimpleOrbits/run.sh
new file mode 100644
index 0000000000000000000000000000000000000000..e1742c357d445faf5fc92330c47300fe2dcff602
--- /dev/null
+++ b/examples/Logger/SimpleOrbits/run.sh
@@ -0,0 +1,11 @@
+#!/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
diff --git a/examples/Logger/SimpleOrbits/simple_orbits.yml b/examples/Logger/SimpleOrbits/simple_orbits.yml
new file mode 100644
index 0000000000000000000000000000000000000000..448f146f3cd38a97fbbf542e541a17438711e43b
--- /dev/null
+++ b/examples/Logger/SimpleOrbits/simple_orbits.yml
@@ -0,0 +1,47 @@
+# 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
diff --git a/examples/Planetary/EarthImpact/earth_impact.yml b/examples/Planetary/EarthImpact/earth_impact.yml
index 83aa4d027bbaf7320353d32ba44cde6c2b4ed252..336a0eb82d570a1544c1344ab2261e35dbf3d59e 100644
--- a/examples/Planetary/EarthImpact/earth_impact.yml
+++ b/examples/Planetary/EarthImpact/earth_impact.yml
@@ -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
diff --git a/examples/Planetary/EarthImpact/make_movie_logger.py b/examples/Planetary/EarthImpact/make_movie_logger.py
new file mode 100644
index 0000000000000000000000000000000000000000..5007543abe8e9984cf8cf86387c526e13b8c948a
--- /dev/null
+++ b/examples/Planetary/EarthImpact/make_movie_logger.py
@@ -0,0 +1,245 @@
+#!/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)
diff --git a/logger/examples/reader_example.py b/logger/examples/reader_example.py
index b15446a6cd3d5b159038b56910eede5f571a6d07..62cca38a94014253e70c6e2f3562e4db99b8cdae 100644
--- a/logger/examples/reader_example.py
+++ b/logger/examples/reader_example.py
@@ -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)
diff --git a/logger/logger_header.c b/logger/logger_header.c
index 2e2ffe07647dedfc204a0f2f31be1a08de1b4627..ea94f3907e94c9dade5449827daa4cead1ed8b0f 100644
--- a/logger/logger_header.c
+++ b/logger/logger_header.c
@@ -184,7 +184,6 @@ void header_read(struct header *h, struct logger_logfile *log) {
/* Check the logfile header's size. */
if (map != (char *)log->log.map + h->offset_first_record) {
- header_print(h);
size_t offset = (char *)map - (char *)log->log.map;
error("Wrong header size (in header %zi, current %zi).",
h->offset_first_record, offset);
diff --git a/logger/logger_index.c b/logger/logger_index.c
index 98368a89b93c707a9e1004bf400e0b898b53b261..80a0858d99dd394964180ad46abcf648076920bb 100644
--- a/logger/logger_index.c
+++ b/logger/logger_index.c
@@ -60,7 +60,9 @@ void logger_index_read_header(struct logger_index *index,
const char *filename) {
/* Open the file */
- message("Reading %s", filename);
+ if (index->reader->verbose > 1) {
+ message("Reading %s", filename);
+ }
logger_index_map_file(index, filename, 0);
/* Read times */
@@ -135,7 +137,7 @@ void logger_index_map_file(struct logger_index *index, const char *filename,
/* Check if need to sort the file */
if (sorted && !index->is_sorted) {
- if (index->reader->verbose > 0) {
+ if (index->reader->verbose > 1) {
message("Sorting the index file.");
}
/* Map the index file */
@@ -155,7 +157,7 @@ void logger_index_map_file(struct logger_index *index, const char *filename,
/* Free the index file before opening it again in read only */
logger_index_free(index);
- if (index->reader->verbose > 0) {
+ if (index->reader->verbose > 1) {
message("Sorting done.");
}
}
diff --git a/logger/logger_particle.c b/logger/logger_particle.c
index 480b0f7f50502a6c9e3f73dfd0ddfd70ee164b07..07cddd89b4e39fef863dcbfd387cfda440c1ae62 100644
--- a/logger/logger_particle.c
+++ b/logger/logger_particle.c
@@ -211,6 +211,104 @@ size_t logger_particle_read(struct logger_particle *part,
return offset;
}
+/**
+ * @brief Compute the quintic hermite spline interpolation.
+ *
+ * @param t0 The time at the left of the interval.
+ * @param x0 The function at the left of the interval.
+ * @param v0 The first derivative at the left of the interval.
+ * @param a0 The second derivative at the left of the interval.
+ * @param t1 The time at the right of the interval.
+ * @param x1 The function at the right of the interval.
+ * @param v1 The first derivative at the right of the interval.
+ * @param a1 The second derivative at the right of the interval.
+ * @param t The time of the interpolation.
+ *
+ * @return The function evaluated at t.
+ */
+double logger_particle_quintic_hermite_spline(double t0, double x0, float v0,
+ float a0, double t1, double x1,
+ float v1, float a1, double t) {
+
+ /* Generates recurring variables */
+ /* Time differences */
+ const double dt = t1 - t0;
+ const double dt2 = dt * dt;
+ const double dt3 = dt2 * dt;
+ const double dt4 = dt3 * dt;
+ const double dt5 = dt4 * dt;
+
+ const double t_t0 = t - t0;
+ const double t_t0_2 = t_t0 * t_t0;
+ const double t_t0_3 = t_t0_2 * t_t0;
+ const double t_t1 = t - t1;
+ const double t_t1_2 = t_t1 * t_t1;
+
+ /* Derivatives */
+ const double v0_dt = v0 * dt;
+ const double a0_dt2 = 0.5 * a0 * dt2;
+ const double v1_dt = v1 * dt;
+ const double a1_dt2 = 0.5 * a1 * dt2;
+
+ /* Do the first 3 terms of the hermite spline */
+ double x = x0 + v0 * t_t0 + 0.5 * a0 * t_t0_2;
+
+ /* Cubic term */
+ x += (x1 - x0 - v0_dt - a0_dt2) * t_t0_3 / dt3;
+
+ /* Quartic term */
+ x += (3. * x0 - 3. * x1 + v1_dt + 2. * v0_dt + a0_dt2) * t_t0_3 * t_t1 / dt4;
+
+ /* Quintic term */
+ x += (6. * x1 - 6. * x0 - 3. * v0_dt - 3. * v1_dt + a1_dt2 - a0_dt2) *
+ t_t0_3 * t_t1_2 / dt5;
+
+ return x;
+}
+
+/**
+ * @brief Compute the cubic hermite spline interpolation.
+ *
+ * @param t0 The time at the left of the interval.
+ * @param v0 The first derivative at the left of the interval.
+ * @param a0 The second derivative at the left of the interval.
+ * @param t1 The time at the right of the interval.
+ * @param v1 The first derivative at the right of the interval.
+ * @param a1 The second derivative at the right of the interval.
+ * @param t The time of the interpolation.
+ *
+ * @return The function evaluated at t.
+ */
+float logger_particle_cubic_hermite_spline(double t0, float v0, float a0,
+ double t1, float v1, float a1,
+ double t) {
+
+ /* Generates recurring variables */
+ /* Time differences */
+ const float dt = t1 - t0;
+ const float dt2 = dt * dt;
+ const float dt3 = dt2 * dt;
+
+ const float t_t0 = t - t0;
+ const float t_t0_2 = t_t0 * t_t0;
+ const float t_t1 = t - t1;
+
+ /* Derivatives */
+ const float a0_dt = a0 * dt;
+ const float a1_dt = a1 * dt;
+
+ /* Do the first 2 terms of the hermite spline */
+ float x = v0 + a0 * t_t0;
+
+ /* Square term */
+ x += (v1 - v0 - a0_dt) * t_t0_2 / dt2;
+
+ /* Cubic term */
+ x += (2. * v0 - 2. * v1 + a1_dt + a0_dt) * t_t0_2 * t_t1 / dt3;
+
+ return x;
+}
+
/**
* @brief interpolate two particles at a given time
*
@@ -245,17 +343,22 @@ void logger_particle_interpolate(struct logger_particle *part_curr,
scaling = (time - part_curr->time) / scaling;
- double tmp;
float ftmp;
/* interpolate vectors. */
- for (size_t i = 0; i < DIM; i++) {
- tmp = (part_next->pos[i] - part_curr->pos[i]);
- part_curr->pos[i] += tmp * scaling;
-
- ftmp = (part_next->vel[i] - part_curr->vel[i]);
- part_curr->vel[i] += ftmp * scaling;
-
+ for (int i = 0; i < 3; i++) {
+ /* Positions */
+ part_curr->pos[i] = logger_particle_quintic_hermite_spline(
+ part_curr->time, part_curr->pos[i], part_curr->vel[i],
+ part_curr->acc[i], part_next->time, part_next->pos[i],
+ part_next->vel[i], part_next->acc[i], time);
+
+ /* Velocities */
+ part_curr->vel[i] = logger_particle_cubic_hermite_spline(
+ part_curr->time, part_curr->vel[i], part_curr->acc[i], part_next->time,
+ part_next->vel[i], part_next->acc[i], time);
+
+ /* Accelerations */
ftmp = (part_next->acc[i] - part_curr->acc[i]);
part_curr->acc[i] += ftmp * scaling;
}
@@ -264,6 +367,12 @@ void logger_particle_interpolate(struct logger_particle *part_curr,
ftmp = (part_next->entropy - part_curr->entropy);
part_curr->entropy += ftmp * scaling;
+ ftmp = (part_next->h - part_curr->h);
+ part_curr->h += ftmp * scaling;
+
+ ftmp = (part_next->density - part_curr->density);
+ part_curr->density += ftmp * scaling;
+
/* set time. */
part_curr->time = time;
}
diff --git a/logger/logger_python_wrapper.c b/logger/logger_python_wrapper.c
index ee8042df8e6dc9fb4914c51fbf729149e4e3a6cf..48973daf969b605b1deed4c55739ad008dbe6df6 100644
--- a/logger/logger_python_wrapper.c
+++ b/logger/logger_python_wrapper.c
@@ -57,7 +57,7 @@ static PyObject *loadSnapshotAtTime(__attribute__((unused)) PyObject *self,
char *basename = NULL;
double time = 0;
- int verbose = 2;
+ int verbose = 1;
/* parse arguments. */
if (!PyArg_ParseTuple(args, "sd|i", &basename, &time, &verbose)) return NULL;
@@ -82,17 +82,14 @@ static PyObject *loadSnapshotAtTime(__attribute__((unused)) PyObject *self,
n_tot += n_parts[i];
}
-#ifdef SWIFT_DEBUG_CHECKS
- message("Found %lu particles", n_tot);
-#endif // SWIFT_DEBUG_CHECKS
+ if (verbose > 0) {
+ message("Found %lu particles", n_tot);
+ }
/* Allocate the output memory */
PyArrayObject *out = (PyArrayObject *)PyArray_SimpleNewFromDescr(
1, &n_tot, logger_particle_descr);
- /* Reference is stolen, therefore need to take it into account */
- Py_INCREF(logger_particle_descr);
-
void *data = PyArray_DATA(out);
/* Allows to use threads */
Py_BEGIN_ALLOW_THREADS;
@@ -125,7 +122,7 @@ static PyObject *getTimeLimits(__attribute__((unused)) PyObject *self,
/* declare variables. */
char *basename = NULL;
- int verbose = 2;
+ int verbose = 1;
/* parse arguments. */
if (!PyArg_ParseTuple(args, "s|i", &basename, &verbose)) return NULL;
@@ -177,6 +174,97 @@ static PyObject *pyReverseOffset(__attribute__((unused)) PyObject *self,
return Py_BuildValue("");
}
+/**
+ * @brief Move forward in time an array of particles.
+ *
+ * <b>filename</b> string filename of the log file.
+ * <b>parts</b> Numpy array containing the particles to evolve.
+ * <b>time</b> Time requested for the particles.
+ * <b>verbose</b> Verbose level
+ *
+ * <b>returns</b> The evolved array of particles.
+ */
+static PyObject *pyMoveForwardInTime(__attribute__((unused)) PyObject *self,
+ PyObject *args) {
+ /* input variables. */
+ char *filename = NULL;
+
+ int verbose = 0;
+ PyArrayObject *parts = NULL;
+ double time = 0;
+ int new_array = 1;
+
+ /* parse the arguments. */
+ if (!PyArg_ParseTuple(args, "sOd|ii", &filename, &parts, &time, &verbose,
+ &new_array))
+ return NULL;
+
+ /* Check parts */
+ if (!PyArray_Check(parts)) {
+ error("Expecting a numpy array of particles.");
+ }
+
+ if (PyArray_NDIM(parts) != 1) {
+ error("Expecting a 1D array of particles.");
+ }
+
+ if (PyArray_TYPE(parts) != logger_particle_descr->type_num) {
+ error("Expecting an array of particles.");
+ }
+
+ /* Create the interpolated array. */
+ PyArrayObject *interp = NULL;
+ if (new_array) {
+ interp =
+ (PyArrayObject *)PyArray_NewLikeArray(parts, NPY_ANYORDER, NULL, 0);
+
+ /* Check if the allocation was fine */
+ if (interp == NULL) {
+ return NULL;
+ }
+
+ /* Reference stolen in PyArray_NewLikeArray => incref */
+ Py_INCREF(PyArray_DESCR(parts));
+ } else {
+ interp = parts;
+ // We return it, therefore one more reference exists.
+ Py_INCREF(interp);
+ }
+
+ /* initialize the reader. */
+ struct logger_reader reader;
+ logger_reader_init(&reader, filename, verbose);
+
+ /* Get the offset of the requested time. */
+ size_t offset = logger_reader_get_next_offset_from_time(&reader, time);
+
+ /* Loop over all the particles */
+ size_t N = PyArray_DIM(parts, 0);
+ for (size_t i = 0; i < N; i++) {
+
+ /* Obtain the required particle records. */
+ struct logger_particle *p = PyArray_GETPTR1(parts, i);
+
+ /* Check that we are really going forward in time. */
+ if (time < p->time) {
+ error("Requesting to go backward in time (%g < %g)", time, p->time);
+ }
+ struct logger_particle new;
+ logger_reader_get_next_particle(&reader, p, &new, offset);
+
+ /* Interpolate the particle. */
+ struct logger_particle *p_ret = PyArray_GETPTR1(interp, i);
+ *p_ret = *p;
+
+ logger_particle_interpolate(p_ret, &new, time);
+ }
+
+ /* Free the reader. */
+ logger_reader_free(&reader);
+
+ return PyArray_Return(interp);
+}
+
/* definition of the method table. */
static PyMethodDef libloggerMethods[] = {
@@ -214,6 +302,24 @@ static PyMethodDef libloggerMethods[] = {
"-------\n\n"
"times: tuple\n"
" time min, time max\n"},
+ {"moveForwardInTime", pyMoveForwardInTime, METH_VARARGS,
+ "Move the particles forward in time.\n\n"
+ "Parameters\n"
+ "----------\n\n"
+ "basename: str\n"
+ " The basename of the index files.\n\n"
+ "parts: np.array\n"
+ " The array of particles.\n\n"
+ "time: double\n"
+ " The requested time for the particles.\n\n"
+ "verbose: int, optional\n"
+ " The verbose level of the loader.\n\n"
+ "new_array: bool, optional\n"
+ " Does the function return a new array (or use the provided one)?\n\n"
+ "Returns\n"
+ "-------\n\n"
+ "parts: np.array\n"
+ " The particles at the requested time.\n"},
{NULL, NULL, 0, NULL} /* Sentinel */
};
diff --git a/logger/logger_reader.c b/logger/logger_reader.c
index 854cf084e1d6c9c6fa5b066ce83b7cec5ad820ad..ce7799c656a9a3a3072b144b9b4ba242db6a4c87 100644
--- a/logger/logger_reader.c
+++ b/logger/logger_reader.c
@@ -378,6 +378,10 @@ double logger_reader_get_time_end(struct logger_reader *reader) {
size_t logger_reader_get_next_offset_from_time(struct logger_reader *reader,
double time) {
size_t ind = time_array_get_index_from_time(&reader->log.times, time);
+ /* We do not want to have the sentiel */
+ if (reader->log.times.size - 2 == ind) {
+ ind -= 1;
+ }
return reader->log.times.records[ind + 1].offset;
}
@@ -427,7 +431,11 @@ void logger_reader_get_next_particle(struct logger_reader *reader,
/* Are we at the end of the file? */
if (next_offset == 0) {
time_array_print(&reader->log.times);
- error("End of file for offset %zi", prev_offset);
+ error(
+ "End of file for particle %lli offset %zi when requesting time %g "
+ "with offset %zi",
+ prev->id, prev_offset,
+ time_array_get_time(&reader->log.times, time_offset), time_offset);
}
next_offset += prev_offset;
diff --git a/src/hydro/Planetary/hydro_part.h b/src/hydro/Planetary/hydro_part.h
index 4c0c4b786c7010a1d8080dbc934d661ca61c6c75..ca0a5b908b3ae140b154c41ce328c0f970a3a6a5 100644
--- a/src/hydro/Planetary/hydro_part.h
+++ b/src/hydro/Planetary/hydro_part.h
@@ -38,6 +38,7 @@
#include "cooling_struct.h"
#include "equation_of_state.h" // For enum material_id
#include "feedback_struct.h"
+#include "logger.h"
#include "star_formation_struct.h"
#include "timestep_limiter_struct.h"
#include "tracers_struct.h"
@@ -78,6 +79,10 @@ struct xpart {
/* Additional data used by the feedback */
struct feedback_part_data feedback_data;
+#ifdef WITH_LOGGER
+ /* Additional data for the particle logger */
+ struct logger_part_data logger_data;
+#endif
} SWIFT_STRUCT_ALIGN;
/**
diff --git a/src/logger.c b/src/logger.c
index 38056dd680a3b6feeb7291196897722aad760f34..c970385e27b561dc466c89bcf41a4845843f69bb 100644
--- a/src/logger.c
+++ b/src/logger.c
@@ -249,6 +249,7 @@ void logger_copy_part_fields(const struct part *p, unsigned int mask,
buff += logger_mask_data[logger_u].size;
mask &= ~logger_mask_data[logger_u].mask;
}
+#endif
/* Particle smoothing length as a single float. */
if (mask & logger_mask_data[logger_h].mask) {
@@ -275,8 +276,6 @@ void logger_copy_part_fields(const struct part *p, unsigned int mask,
mask &= ~logger_mask_data[logger_consts].mask;
}
-#endif
-
/* Special flags */
if (mask & logger_mask_data[logger_special_flags].mask) {
memcpy(buff, &special_flags, logger_mask_data[logger_special_flags].size);