Logger example

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

examples/Logger/SimpleOrbits/makeIC.py

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

examples/Logger/SimpleOrbits/mpl_style

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

examples/Logger/SimpleOrbits/plotSolution.py

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

examples/Logger/SimpleOrbits/run.sh

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

examples/Logger/SimpleOrbits/simple_orbits.yml

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

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
+    planetary_use_Til:  1               # Whether to initialise the Tillotson EOS
\ No newline at end of file
+
+# 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

examples/Planetary/EarthImpact/make_movie_logger.py

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

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)