From 3c3d6586d5138c155e35576346c150383c3ac995 Mon Sep 17 00:00:00 2001 From: loikki Date: Fri, 7 Feb 2020 22:50:50 +0100 Subject: [PATCH 01/12] Logger: implement simple orbits example --- examples/Logger/SimpleOrbits/makeIC.py | 99 +++++++++++++++++++ examples/Logger/SimpleOrbits/mpl_style | 26 +++++ examples/Logger/SimpleOrbits/plotSolution.py | 99 +++++++++++++++++++ examples/Logger/SimpleOrbits/run.sh | 11 +++ .../Logger/SimpleOrbits/simple_orbits.yml | 46 +++++++++ 5 files changed, 281 insertions(+) create mode 100644 examples/Logger/SimpleOrbits/makeIC.py create mode 100644 examples/Logger/SimpleOrbits/mpl_style create mode 100644 examples/Logger/SimpleOrbits/plotSolution.py create mode 100644 examples/Logger/SimpleOrbits/run.sh create mode 100644 examples/Logger/SimpleOrbits/simple_orbits.yml diff --git a/examples/Logger/SimpleOrbits/makeIC.py b/examples/Logger/SimpleOrbits/makeIC.py new file mode 100644 index 000000000..9642d4145 --- /dev/null +++ b/examples/Logger/SimpleOrbits/makeIC.py @@ -0,0 +1,99 @@ +############################################################################### +# 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 . +# +############################################################################## + +from h5py import File +import numpy as np +from astropy import units, constants + +np.random.seed(50) + +# parameters + +filename = "simple_orbits.hdf5" +num_part = 3 +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 ids +ids = np.arange(num_part) + +# 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] + +# File +file = File(filename, 'w') + +# Header +grp = file.create_group("/Header") +grp.attrs["BoxSize"] = [boxsize]*3 +grp.attrs["NumPart_Total"] = [0, num_part, 0, 0, 0, 0] +grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0] +grp.attrs["NumPart_ThisFile"] = [0, num_part, 0, 0, 0, 0] +grp.attrs["Time"] = 0.0 +grp.attrs["NumFilesPerSnapshot"] = 1 +grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0] +grp.attrs["Flag_Entropy_ICs"] = 0 +grp.attrs["Dimension"] = 3 + +# Units +grp = file.create_group("/Units") +grp.attrs["Unit length in cgs (U_L)"] = u_l +grp.attrs["Unit mass in cgs (U_M)"] = u_m +grp.attrs["Unit time in cgs (U_t)"] = u_t +grp.attrs["Unit current in cgs (U_I)"] = 1. +grp.attrs["Unit temperature in cgs (U_T)"] = 1. + +# Particle group +grp = file.create_group("/PartType1") +grp.create_dataset('Coordinates', data=coords, dtype='d') +grp.create_dataset('Velocities', data=v, dtype='f') +grp.create_dataset('Masses', data=m, dtype='f') +grp.create_dataset('ParticleIDs', data=ids, dtype='L') + +file.close() diff --git a/examples/Logger/SimpleOrbits/mpl_style b/examples/Logger/SimpleOrbits/mpl_style new file mode 100644 index 000000000..8e4296931 --- /dev/null +++ b/examples/Logger/SimpleOrbits/mpl_style @@ -0,0 +1,26 @@ +# Line Properties +lines.linewidth: 1. + +# Font options +font.size: 8 +font.family: STIXGeneral +mathtext.fontset: stix + +# LaTeX options +text.usetex: False + +# Legend options +legend.frameon: False +legend.fontsize: 8 + +# 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 000000000..a1c7e399e --- /dev/null +++ b/examples/Logger/SimpleOrbits/plotSolution.py @@ -0,0 +1,99 @@ +############################################################################### +# 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 . +# +############################################################################## + +from h5py import File +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) + + +def getErrorSnapshots(): + filenames = glob("simple_orbits_*.hdf5") + N = len(filenames) + filenames.sort() + + r2 = np.zeros((N, makeIC.num_part)) + t = np.zeros(N) + + for i, f in enumerate(filenames): + if i % 10 == 0: + print("File %i / %i" % (i, N)) + f = File(f, "r") + pos = f["PartType1/Coordinates"][:] + pos -= center + ids = f["PartType1/ParticleIDs"][:] + print(ids) + + t[i] = f["Header"].attrs["Time"] + for j in range(makeIC.num_part): + ind = ids == j + print(ind) + r2[i, j] = np.sum(pos[ind, :]**2) + + # compute the plotting variables + r = np.sqrt(r2) + av = r.mean(axis=0) + eps = (r - av) / r[0, :] + return t, eps + + +def getErrorLogger(): + basename = "index" + N = 1000 + t_min, t_max = logger.getTimeLimits(basename) + times = np.linspace(t_min, t_max, N) + r2 = np.zeros((N, makeIC.num_part)) + + for i, t in enumerate(times): + parts = logger.loadSnapshotAtTime(basename, t) + pos = parts["positions"] + pos -= center + + for j in range(makeIC.num_part): + ind = parts["ids"] == j + r2[i, j] = np.sum(pos[ind, :]**2) + + # compute the plotting variables + r = np.sqrt(r2) + av = r.mean(axis=0) + eps = (r - av) / r[0, :] + return times, eps + + +t_log, eps_log = getErrorLogger() +t_snap, eps_snap = getErrorSnapshots() + +plt.figure() +colors = ["b", "r", "m", "g", "k"] +for i in range(makeIC.num_part): + plt.plot(t_snap, eps_snap[:, i], ":", color=colors[i]) + plt.plot(t_log, eps_log[:, i], "-", color=colors[i]) +plt.xlabel("Time [yr]") +plt.ylabel("Relative error on the radius") +plt.legend(["Snapshot", "Logger"]) +plt.show() diff --git a/examples/Logger/SimpleOrbits/run.sh b/examples/Logger/SimpleOrbits/run.sh new file mode 100644 index 000000000..768c40df1 --- /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=4 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 000000000..93ce1fbf2 --- /dev/null +++ b/examples/Logger/SimpleOrbits/simple_orbits.yml @@ -0,0 +1,46 @@ +# 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: 6. # 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: 1e-2 # 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: 0.05 # Time difference between consecutive outputs (in internal units) + +# Parameters governing the conserved quantities statistics +Statistics: + delta_time: 0.01 # 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: 0.03 # Dimensionless pre-factor for the time-step condition + + +# 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 -- GitLab From 647558c56e8b315e29f347619777c9ed78f13e1f Mon Sep 17 00:00:00 2001 From: loikki Date: Sun, 9 Feb 2020 18:46:26 +0100 Subject: [PATCH 02/12] Logger: messages are shown only if verbose level > 0 --- examples/Logger/SimpleOrbits/plotSolution.py | 34 ++++++++++++++------ logger/logger_header.c | 1 - logger/logger_index.c | 8 +++-- logger/logger_python_wrapper.c | 10 +++--- 4 files changed, 34 insertions(+), 19 deletions(-) diff --git a/examples/Logger/SimpleOrbits/plotSolution.py b/examples/Logger/SimpleOrbits/plotSolution.py index a1c7e399e..3f5352775 100644 --- a/examples/Logger/SimpleOrbits/plotSolution.py +++ b/examples/Logger/SimpleOrbits/plotSolution.py @@ -30,15 +30,17 @@ import liblogger as logger plt.style.use("mpl_style") center = np.array([0.5 * makeIC.boxsize]*3) +id_focus = 0 -def getErrorSnapshots(): +def getSnapshots(): filenames = glob("simple_orbits_*.hdf5") N = len(filenames) filenames.sort() r2 = np.zeros((N, makeIC.num_part)) t = np.zeros(N) + p = np.zeros((N, 3)) for i, f in enumerate(filenames): if i % 10 == 0: @@ -47,46 +49,50 @@ def getErrorSnapshots(): pos = f["PartType1/Coordinates"][:] pos -= center ids = f["PartType1/ParticleIDs"][:] - print(ids) t[i] = f["Header"].attrs["Time"] for j in range(makeIC.num_part): ind = ids == j - print(ind) r2[i, j] = np.sum(pos[ind, :]**2) + if j == id_focus: + p[i, :] = pos[ind, :] # compute the plotting variables r = np.sqrt(r2) av = r.mean(axis=0) eps = (r - av) / r[0, :] - return t, eps + return t, eps, p -def getErrorLogger(): +def getLogger(): basename = "index" N = 1000 - t_min, t_max = logger.getTimeLimits(basename) + verbose = 0 + t_min, t_max = logger.getTimeLimits(basename, verbose) times = np.linspace(t_min, t_max, N) r2 = np.zeros((N, makeIC.num_part)) + p = np.zeros((N, 3)) for i, t in enumerate(times): - parts = logger.loadSnapshotAtTime(basename, t) + parts = logger.loadSnapshotAtTime(basename, t, verbose) pos = parts["positions"] pos -= center for j in range(makeIC.num_part): ind = parts["ids"] == j r2[i, j] = np.sum(pos[ind, :]**2) + if j == id_focus: + p[i, :] = pos[ind, :] # compute the plotting variables r = np.sqrt(r2) av = r.mean(axis=0) eps = (r - av) / r[0, :] - return times, eps + return times, eps, p -t_log, eps_log = getErrorLogger() -t_snap, eps_snap = getErrorSnapshots() +t_log, eps_log, p_log = getLogger() +t_snap, eps_snap, p_snap = getSnapshots() plt.figure() colors = ["b", "r", "m", "g", "k"] @@ -96,4 +102,12 @@ for i in range(makeIC.num_part): plt.xlabel("Time [yr]") plt.ylabel("Relative error on the radius") plt.legend(["Snapshot", "Logger"]) + + +plt.figure() +plt.plot(p_log[:, 0], p_log[:, 1], label="Logger") +plt.plot(p_snap[:, 0], p_snap[:, 1], label="Snapshot") +plt.xlabel("Position [AU]") +plt.ylabel("Position [AU]") +plt.legend() plt.show() diff --git a/logger/logger_header.c b/logger/logger_header.c index 2e2ffe076..ea94f3907 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 98368a89b..80a0858d9 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_python_wrapper.c b/logger/logger_python_wrapper.c index ee8042df8..fe3081b59 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,9 +82,9 @@ 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( @@ -125,7 +125,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; -- GitLab From 7bb24c5f5a48cae7cb9095961c615e1c318dd7fe Mon Sep 17 00:00:00 2001 From: loikki Date: Wed, 12 Feb 2020 11:47:07 +0100 Subject: [PATCH 03/12] Logger: implement hermite interpolation and add first logger example --- examples/Logger/SimpleOrbits/makeIC.py | 2 +- examples/Logger/SimpleOrbits/mpl_style | 2 +- examples/Logger/SimpleOrbits/plotSolution.py | 87 +++++++++---- .../Logger/SimpleOrbits/simple_orbits.yml | 11 +- logger/logger_particle.c | 117 ++++++++++++++++-- logger/logger_python_wrapper.c | 104 ++++++++++++++++ logger/logger_reader.c | 7 +- 7 files changed, 289 insertions(+), 41 deletions(-) diff --git a/examples/Logger/SimpleOrbits/makeIC.py b/examples/Logger/SimpleOrbits/makeIC.py index 9642d4145..6d12ea46e 100644 --- a/examples/Logger/SimpleOrbits/makeIC.py +++ b/examples/Logger/SimpleOrbits/makeIC.py @@ -26,7 +26,7 @@ np.random.seed(50) # parameters filename = "simple_orbits.hdf5" -num_part = 3 +num_part = 1 masses = 1. # If changed, need to update simple_orbits.yml M = units.solMass.to("earthMass") diff --git a/examples/Logger/SimpleOrbits/mpl_style b/examples/Logger/SimpleOrbits/mpl_style index 8e4296931..80e125f52 100644 --- a/examples/Logger/SimpleOrbits/mpl_style +++ b/examples/Logger/SimpleOrbits/mpl_style @@ -10,7 +10,7 @@ mathtext.fontset: stix text.usetex: False # Legend options -legend.frameon: False +legend.frameon: True legend.fontsize: 8 # Figure options for publications diff --git a/examples/Logger/SimpleOrbits/plotSolution.py b/examples/Logger/SimpleOrbits/plotSolution.py index 3f5352775..03702eca1 100644 --- a/examples/Logger/SimpleOrbits/plotSolution.py +++ b/examples/Logger/SimpleOrbits/plotSolution.py @@ -32,13 +32,16 @@ plt.style.use("mpl_style") center = np.array([0.5 * makeIC.boxsize]*3) id_focus = 0 +G = 1.189972e-04 +M = 3.329460e+05 + def getSnapshots(): filenames = glob("simple_orbits_*.hdf5") N = len(filenames) filenames.sort() - r2 = np.zeros((N, makeIC.num_part)) + E = np.zeros((N, makeIC.num_part)) t = np.zeros(N) p = np.zeros((N, 3)) @@ -48,65 +51,99 @@ def getSnapshots(): f = File(f, "r") pos = f["PartType1/Coordinates"][:] pos -= center + vel = f["PartType1/Velocities"][:] ids = f["PartType1/ParticleIDs"][:] t[i] = f["Header"].attrs["Time"] for j in range(makeIC.num_part): ind = ids == j - r2[i, j] = np.sum(pos[ind, :]**2) + r = np.sum(pos[ind, :]**2)**0.5 + v2 = np.sum(vel[ind, :]**2) + E[i, ind] = 0.5 * v2 - G * M / r if j == id_focus: p[i, :] = pos[ind, :] # compute the plotting variables - r = np.sqrt(r2) - av = r.mean(axis=0) - eps = (r - av) / r[0, :] + eps = (E - E[0, :]) / E[0, :] return t, eps, p +def getStatistics(): + data = np.genfromtxt("energy.txt", names=True) + + times = data["Time"] + E = data["E_tot"] + eps = (E - E[0]) / E[0] + return times, eps + + def getLogger(): 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) - r2 = np.zeros((N, makeIC.num_part)) + + # Create output arrays + E = np.zeros((N, makeIC.num_part)) + E_parts = np.zeros((N, makeIC.num_part)) p = np.zeros((N, 3)) + t_parts = np.zeros((N, makeIC.num_part)) + + # Read the particles + parts = logger.loadSnapshotAtTime( + basename, times[0], verbose) for i, t in enumerate(times): - parts = logger.loadSnapshotAtTime(basename, t, verbose) - pos = parts["positions"] - pos -= center + interp = logger.moveForwardInTime( + basename, parts, t, verbose) + rel_pos = interp["positions"] - center + v = interp["velocities"] + + rel_pos_parts = parts["positions"] - center + v_parts = parts["velocities"] for j in range(makeIC.num_part): - ind = parts["ids"] == j - r2[i, j] = np.sum(pos[ind, :]**2) + # do interp + ind = interp["ids"] == j + r = np.sum(rel_pos[ind, :]**2)**0.5 + v2 = np.sum(v[ind, :]**2) + E[i, ind] = 0.5 * v2 - G * M / r if j == id_focus: - p[i, :] = pos[ind, :] + p[i, :] = rel_pos[ind, :] + + # do parts + ind = parts["ids"] == j + r = np.sum(rel_pos_parts[ind, :]**2)**0.5 + v2 = np.sum(v_parts[ind, :]**2) + E_parts[i, ind] = 0.5 * v2 - G * M / r + t_parts[i, ind] = parts["times"][ind] # compute the plotting variables - r = np.sqrt(r2) - av = r.mean(axis=0) - eps = (r - av) / r[0, :] - return times, eps, p + eps = (E - E[0, :]) / E[0, :] + eps_parts = (E_parts - E_parts[0, :]) / E_parts[0, :] + return times, eps, p, eps_parts, t_parts -t_log, eps_log, p_log = getLogger() +t_log, eps_log, p_log, eps_parts, t_parts = getLogger() t_snap, eps_snap, p_snap = getSnapshots() +t_stat, eps_stat = getStatistics() plt.figure() -colors = ["b", "r", "m", "g", "k"] -for i in range(makeIC.num_part): - plt.plot(t_snap, eps_snap[:, i], ":", color=colors[i]) - plt.plot(t_log, eps_log[:, i], "-", color=colors[i]) +plt.plot(t_snap, eps_snap, ".", label="Snapshot") +plt.plot(t_parts, eps_parts, "x", label="Logger") +plt.plot(t_log, eps_log, "--", label="Logger (Interpolation)") +plt.plot(t_stat, eps_stat, "-", label="Statistics") plt.xlabel("Time [yr]") -plt.ylabel("Relative error on the radius") -plt.legend(["Snapshot", "Logger"]) +plt.ylabel("Relative error on the specific energy") +plt.legend(ncol=2) plt.figure() -plt.plot(p_log[:, 0], p_log[:, 1], label="Logger") -plt.plot(p_snap[:, 0], p_snap[:, 1], label="Snapshot") +plt.plot(p_snap[:, 0], p_snap[:, 1], ".", label="Snapshot") +plt.plot(p_log[:, 0], p_log[:, 1], "--", label="Logger") plt.xlabel("Position [AU]") plt.ylabel("Position [AU]") plt.legend() diff --git a/examples/Logger/SimpleOrbits/simple_orbits.yml b/examples/Logger/SimpleOrbits/simple_orbits.yml index 93ce1fbf2..d6023e8b9 100644 --- a/examples/Logger/SimpleOrbits/simple_orbits.yml +++ b/examples/Logger/SimpleOrbits/simple_orbits.yml @@ -9,7 +9,7 @@ InternalUnitSystem: # Parameters governing the time integration TimeIntegration: time_begin: 0. # The starting time of the simulation (in internal units). - time_end: 6. # The end time of the simulation (in internal units). + time_end: 0.1 # 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: 1e-2 # The maximal time-step size of the simulation (in internal units). @@ -17,11 +17,11 @@ TimeIntegration: 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: 0.05 # Time difference between consecutive outputs (in internal units) + delta_time: 3e-3 # Time difference between consecutive outputs (in internal units) # Parameters governing the conserved quantities statistics Statistics: - delta_time: 0.01 # Time between statistics output + delta_time: 1e-5 # Time between statistics output # Parameters related to the initial conditions @@ -35,12 +35,13 @@ 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: 0.03 # Dimensionless pre-factor for the time-step condition + timestep_mult: 1e-5 # Dimensionless pre-factor for the time-step condition # Parameters governing the logger snapshot system Logger: - delta_step: 10 # Update the particle log every this many updates + delta_step: 500 # 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/logger/logger_particle.c b/logger/logger_particle.c index 480b0f7f5..60ba32556 100644 --- a/logger/logger_particle.c +++ b/logger/logger_particle.c @@ -211,6 +211,103 @@ 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 = a0 * dt2; + const double v1_dt = v1 * dt; + const double a1_dt2 = a1 * dt2; + + /* Do the first 3 terms of the hermite spline */ + double x = x0 + v0 * t_t0 + 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 +342,21 @@ 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; } diff --git a/logger/logger_python_wrapper.c b/logger/logger_python_wrapper.c index fe3081b59..1774d01b1 100644 --- a/logger/logger_python_wrapper.c +++ b/logger/logger_python_wrapper.c @@ -177,6 +177,91 @@ static PyObject *pyReverseOffset(__attribute__((unused)) PyObject *self, return Py_BuildValue(""); } + +/** + * @brief Move forward in time an array of particles. + * + * filename string filename of the log file. + * parts Numpy array containing the particles to evolve. + * time Time requested for the particles. + * verbose Verbose level + * + * returns 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); + } + 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 +299,25 @@ 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 854cf084e..c9783b1e3 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,8 @@ 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 offset %zi when requesting time with offset %zi", + prev_offset, time_offset); } next_offset += prev_offset; -- GitLab From 9e1b3fab27370cb06d5ca3624d7839a34575279f Mon Sep 17 00:00:00 2001 From: loikki Date: Thu, 13 Feb 2020 10:00:52 +0100 Subject: [PATCH 04/12] Logger: fix interpolation + update plots --- examples/Logger/SimpleOrbits/mpl_style | 2 +- examples/Logger/SimpleOrbits/plotSolution.py | 204 +++++++++++++++--- .../Logger/SimpleOrbits/simple_orbits.yml | 4 +- logger/logger_particle.c | 6 +- 4 files changed, 177 insertions(+), 39 deletions(-) diff --git a/examples/Logger/SimpleOrbits/mpl_style b/examples/Logger/SimpleOrbits/mpl_style index 80e125f52..b8fc86777 100644 --- a/examples/Logger/SimpleOrbits/mpl_style +++ b/examples/Logger/SimpleOrbits/mpl_style @@ -7,7 +7,7 @@ font.family: STIXGeneral mathtext.fontset: stix # LaTeX options -text.usetex: False +text.usetex: True # Legend options legend.frameon: True diff --git a/examples/Logger/SimpleOrbits/plotSolution.py b/examples/Logger/SimpleOrbits/plotSolution.py index 03702eca1..ad32bd8f1 100644 --- a/examples/Logger/SimpleOrbits/plotSolution.py +++ b/examples/Logger/SimpleOrbits/plotSolution.py @@ -21,6 +21,7 @@ from h5py import File import numpy as np from glob import glob import matplotlib.pyplot as plt +from scipy.integrate import solve_ivp import makeIC import sys sys.path.append("../../../logger/.libs/") @@ -31,23 +32,87 @@ plt.style.use("mpl_style") center = np.array([0.5 * makeIC.boxsize]*3) id_focus = 0 +# Plot the absolute positions (+vel) or the difference +# to a solution computed with a RK solver. +diff_to_sol = False -G = 1.189972e-04 -M = 3.329460e+05 +# 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 getSnapshots(): + +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 Evolve(y0, t): + """ + Compute the solution to the gravitational problem. + + Parameters + ---------- + + y0: np.array + Initial conditions [x, y, vx, vy] + + t: np.array + Time of the outputs. + """ + y = solve_ivp(gravity, (t[0], t[-1]), y0, t_eval=t) + return y.y + + +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): - if i % 10 == 0: - print("File %i / %i" % (i, N)) + # get the data from the file f = File(f, "r") pos = f["PartType1/Coordinates"][:] pos -= center @@ -55,29 +120,60 @@ def getSnapshots(): ids = f["PartType1/ParticleIDs"][:] t[i] = f["Header"].attrs["Time"] + + # loop over all the particles + # (avoids comparing different particles) for j in range(makeIC.num_part): + # Get the current id ind = ids == j + + # Compute the energy r = np.sum(pos[ind, :]**2)**0.5 v2 = np.sum(vel[ind, :]**2) E[i, ind] = 0.5 * v2 - G * M / r + + # Get the pos / vel of the required particle if j == 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] + y = Evolve(y0, t) # compute the plotting variables - eps = (E - E[0, :]) / E[0, :] - return t, eps, p + plt.figure(fig_1.number) + plotRelative(t, E, ".", label="Snapshot") + + plt.figure(fig_2.number) + if diff_to_sol: + p[:, :2] = p[:, :2] - y[:2, :].transpose() + plt.plot(p[:, 0], p[:, 1], ".", label="Snapshot") + plt.figure(fig_3.number) + if diff_to_sol: + v[:, :2] = v[:, :2] - y[2:, :].transpose() + plt.plot(v[:, 0], v[:, 1], ".", label="Snapshot") -def getStatistics(): + +def doStatistics(): + """ + Do the plots with the energy output file. + """ data = np.genfromtxt("energy.txt", names=True) times = data["Time"] E = data["E_tot"] - eps = (E - E[0]) / E[0] - return times, eps + plt.figure(fig_1.number) + plotRelative(times, E, "-", label="Statistics") -def getLogger(): +def doLogger(): + """ + Read the logfile and plot the corresponding variables. + """ basename = "index" N = 1000 verbose = 0 @@ -90,61 +186,103 @@ def getLogger(): 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) rel_pos = interp["positions"] - center - v = interp["velocities"] + vel = interp["velocities"] rel_pos_parts = parts["positions"] - center - v_parts = parts["velocities"] + vel_parts = parts["velocities"] for j in range(makeIC.num_part): - # do interp + # Compute the interpolated variables ind = interp["ids"] == j r = np.sum(rel_pos[ind, :]**2)**0.5 - v2 = np.sum(v[ind, :]**2) + v2 = np.sum(vel[ind, :]**2) E[i, ind] = 0.5 * v2 - G * M / r if j == id_focus: p[i, :] = rel_pos[ind, :] + v[i, :] = vel[ind, :] - # do parts + # Compute the variables of the last record ind = parts["ids"] == j r = np.sum(rel_pos_parts[ind, :]**2)**0.5 - v2 = np.sum(v_parts[ind, :]**2) + v2 = np.sum(vel_parts[ind, :]**2) E_parts[i, ind] = 0.5 * v2 - G * M / r t_parts[i, ind] = parts["times"][ind] + if j == id_focus: + p_parts[i, :] = rel_pos_parts[ind, :] + v_parts[i, :] = vel_parts[ind, :] # compute the plotting variables - eps = (E - E[0, :]) / E[0, :] - eps_parts = (E_parts - E_parts[0, :]) / E_parts[0, :] - return times, eps, p, eps_parts, t_parts + 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) + y = Evolve(y0, t_parts) + if diff_to_sol: + p_parts = p_parts[ind, :2] - y[:2, :].transpose() + # plot the solution + plt.figure(fig_2.number) + plt.plot(p_parts[:, 0], p_parts[:, 1], "x", label="Logger") -t_log, eps_log, p_log, eps_parts, t_parts = getLogger() -t_snap, eps_snap, p_snap = getSnapshots() -t_stat, eps_stat = getStatistics() + plt.figure(fig_3.number) + if diff_to_sol: + v_parts = v_parts[ind, :2] - y[2:, :].transpose() + plt.plot(v_parts[:, 0], v_parts[:, 1], "x", label="Logger") -plt.figure() -plt.plot(t_snap, eps_snap, ".", label="Snapshot") -plt.plot(t_parts, eps_parts, "x", label="Logger") -plt.plot(t_log, eps_log, "--", label="Logger (Interpolation)") -plt.plot(t_stat, eps_stat, "-", label="Statistics") + # Compute the solution + y0 = np.zeros(4) + y0[:2] = p[0, :2] + y0[2:] = v[0, :2] + y = Evolve(y0, times) + if diff_to_sol: + p[:, :2] = p[:, :2] - y[:2, :].transpose() + plt.figure(fig_2.number) + plt.plot(p[:, 0], p[:, 1], "--", label="Logger (Interpolation)") + + plt.figure(fig_3.number) + if diff_to_sol: + v[:, :2] = v[:, :2] - y[2:, :].transpose() + plt.plot(v[:, 0], v[:, 1], "--", label="Logger (Interpolation)") + + +# do all the plots +doStatistics() +doSnapshots() +doLogger() + +# add text +plt.figure(fig_1.number) plt.xlabel("Time [yr]") -plt.ylabel("Relative error on the specific energy") +plt.ylabel(r"$\frac{E - E(t=0)}{E(t=0)}$") plt.legend(ncol=2) - -plt.figure() -plt.plot(p_snap[:, 0], p_snap[:, 1], ".", label="Snapshot") -plt.plot(p_log[:, 0], p_log[:, 1], "--", label="Logger") +plt.figure(fig_2.number) plt.xlabel("Position [AU]") plt.ylabel("Position [AU]") plt.legend() + +plt.figure(fig_3.number) +plt.xlabel("Velocity [AU / yr]") +plt.ylabel("Velocity [AU / yr]") +plt.legend() + plt.show() diff --git a/examples/Logger/SimpleOrbits/simple_orbits.yml b/examples/Logger/SimpleOrbits/simple_orbits.yml index d6023e8b9..dca898c4a 100644 --- a/examples/Logger/SimpleOrbits/simple_orbits.yml +++ b/examples/Logger/SimpleOrbits/simple_orbits.yml @@ -17,7 +17,7 @@ TimeIntegration: 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: 3e-3 # Time difference between consecutive outputs (in internal units) + delta_time: 4e-3 # Time difference between consecutive outputs (in internal units) # Parameters governing the conserved quantities statistics Statistics: @@ -40,7 +40,7 @@ PointMassPotential: # Parameters governing the logger snapshot system Logger: - delta_step: 500 # Update the particle log every this many updates + delta_step: 2000 # 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/logger/logger_particle.c b/logger/logger_particle.c index 60ba32556..b07b9a570 100644 --- a/logger/logger_particle.c +++ b/logger/logger_particle.c @@ -246,12 +246,12 @@ double logger_particle_quintic_hermite_spline( /* Derivatives */ const double v0_dt = v0 * dt; - const double a0_dt2 = a0 * dt2; + const double a0_dt2 = 0.5 * a0 * dt2; const double v1_dt = v1 * dt; - const double a1_dt2 = a1 * dt2; + const double a1_dt2 = 0.5 * a1 * dt2; /* Do the first 3 terms of the hermite spline */ - double x = x0 + v0 * t_t0 + a0 * t_t0_2; + 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; -- GitLab From 61ed61505e0031a93a4d6e886dd674bd2d29b80f Mon Sep 17 00:00:00 2001 From: loikki Date: Thu, 13 Feb 2020 12:55:50 +0100 Subject: [PATCH 05/12] Logger: save figure --- examples/Logger/SimpleOrbits/plotSolution.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/examples/Logger/SimpleOrbits/plotSolution.py b/examples/Logger/SimpleOrbits/plotSolution.py index ad32bd8f1..0531d8d4a 100644 --- a/examples/Logger/SimpleOrbits/plotSolution.py +++ b/examples/Logger/SimpleOrbits/plotSolution.py @@ -274,15 +274,16 @@ 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.legend() +plt.savefig("Positions.png") plt.figure(fig_3.number) plt.xlabel("Velocity [AU / yr]") plt.ylabel("Velocity [AU / yr]") plt.legend() - -plt.show() +plt.savefig("Velocities.png") -- GitLab From 971234a00dd54bf03369894d9568776edc8a3c81 Mon Sep 17 00:00:00 2001 From: loikki Date: Thu, 20 Feb 2020 16:25:59 +0100 Subject: [PATCH 06/12] Logger: add movie for the planetary impact example --- .../EarthImpact/make_movie_logger.py | 255 ++++++++++++++++++ 1 file changed, 255 insertions(+) create mode 100644 examples/Planetary/EarthImpact/make_movie_logger.py diff --git a/examples/Planetary/EarthImpact/make_movie_logger.py b/examples/Planetary/EarthImpact/make_movie_logger.py new file mode 100644 index 000000000..24c03db77 --- /dev/null +++ b/examples/Planetary/EarthImpact/make_movie_logger.py @@ -0,0 +1,255 @@ +#!/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) + +""" +TODO: + - change speed (back in time?) + - zoom in / out + - trace particle + - zoom on different area (different time scale and show it?) + - Add length scale + - show large scale +""" -- GitLab From 3f2298d501d0c56664aa81c0dbab6138afc55a24 Mon Sep 17 00:00:00 2001 From: loikki Date: Thu, 20 Feb 2020 16:27:16 +0100 Subject: [PATCH 07/12] Logger: remove comments --- examples/Planetary/EarthImpact/make_movie_logger.py | 12 +----------- 1 file changed, 1 insertion(+), 11 deletions(-) diff --git a/examples/Planetary/EarthImpact/make_movie_logger.py b/examples/Planetary/EarthImpact/make_movie_logger.py index 24c03db77..5007543ab 100644 --- a/examples/Planetary/EarthImpact/make_movie_logger.py +++ b/examples/Planetary/EarthImpact/make_movie_logger.py @@ -237,19 +237,9 @@ def main(): if __name__ == "__main__": - # 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) - -""" -TODO: - - change speed (back in time?) - - zoom in / out - - trace particle - - zoom on different area (different time scale and show it?) - - Add length scale - - show large scale -""" -- GitLab From b2715eb2f1c00444bfcc00f480be90f08a66d9e5 Mon Sep 17 00:00:00 2001 From: loikki Date: Wed, 26 Feb 2020 16:22:16 +0100 Subject: [PATCH 08/12] Logger: cleanup simpleorbits example --- examples/GEAR/ZoomIn/zoom_in.yml | 2 +- examples/Logger/SimpleOrbits/makeIC.py | 44 +++++-------------- examples/Logger/SimpleOrbits/mpl_style | 6 +-- examples/Logger/SimpleOrbits/plotSolution.py | 18 +++++--- examples/Logger/SimpleOrbits/run.sh | 2 +- .../Logger/SimpleOrbits/simple_orbits.yml | 12 ++--- .../Planetary/EarthImpact/earth_impact.yml | 9 +++- logger/examples/reader_example.py | 2 +- logger/logger_particle.c | 13 ++++-- logger/logger_particle.h | 2 +- logger/logger_python_wrapper.c | 12 +++-- logger/logger_reader.c | 10 +++-- src/hydro/Planetary/hydro_part.h | 5 +++ src/logger.c | 5 +-- 14 files changed, 74 insertions(+), 68 deletions(-) diff --git a/examples/GEAR/ZoomIn/zoom_in.yml b/examples/GEAR/ZoomIn/zoom_in.yml index 1d01e6abd..e0927d344 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 index 6d12ea46e..a24184925 100644 --- a/examples/Logger/SimpleOrbits/makeIC.py +++ b/examples/Logger/SimpleOrbits/makeIC.py @@ -17,9 +17,10 @@ # ############################################################################## -from h5py import File import numpy as np -from astropy import units, constants +from astropy import units +from swiftsimio import Writer +import unyt np.random.seed(50) @@ -54,9 +55,6 @@ coords += boxsize * 0.5 # generate the masses m = np.ones(num_part) * masses -# generate the ids -ids = np.arange(num_part) - # generate the velocities sign = np.random.rand(num_part) sign[sign < 0.5] = -1 @@ -66,34 +64,12 @@ 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] -# File -file = File(filename, 'w') - -# Header -grp = file.create_group("/Header") -grp.attrs["BoxSize"] = [boxsize]*3 -grp.attrs["NumPart_Total"] = [0, num_part, 0, 0, 0, 0] -grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0] -grp.attrs["NumPart_ThisFile"] = [0, num_part, 0, 0, 0, 0] -grp.attrs["Time"] = 0.0 -grp.attrs["NumFilesPerSnapshot"] = 1 -grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0] -grp.attrs["Flag_Entropy_ICs"] = 0 -grp.attrs["Dimension"] = 3 - -# Units -grp = file.create_group("/Units") -grp.attrs["Unit length in cgs (U_L)"] = u_l -grp.attrs["Unit mass in cgs (U_M)"] = u_m -grp.attrs["Unit time in cgs (U_t)"] = u_t -grp.attrs["Unit current in cgs (U_I)"] = 1. -grp.attrs["Unit temperature in cgs (U_T)"] = 1. +# Write the snapshot +units = unyt.UnitSystem("Planets", unyt.AU, unyt.mearth, unyt.yr) -# Particle group -grp = file.create_group("/PartType1") -grp.create_dataset('Coordinates', data=coords, dtype='d') -grp.create_dataset('Velocities', data=v, dtype='f') -grp.create_dataset('Masses', data=m, dtype='f') -grp.create_dataset('ParticleIDs', data=ids, dtype='L') +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 -file.close() +snapshot.write(filename) diff --git a/examples/Logger/SimpleOrbits/mpl_style b/examples/Logger/SimpleOrbits/mpl_style index b8fc86777..9b85038b2 100644 --- a/examples/Logger/SimpleOrbits/mpl_style +++ b/examples/Logger/SimpleOrbits/mpl_style @@ -1,8 +1,8 @@ # Line Properties -lines.linewidth: 1. +lines.linewidth: 2. # Font options -font.size: 8 +font.size: 10 font.family: STIXGeneral mathtext.fontset: stix @@ -11,7 +11,7 @@ text.usetex: True # Legend options legend.frameon: True -legend.fontsize: 8 +legend.fontsize: 10 # Figure options for publications figure.dpi: 300 diff --git a/examples/Logger/SimpleOrbits/plotSolution.py b/examples/Logger/SimpleOrbits/plotSolution.py index 0531d8d4a..ee06b6a5b 100644 --- a/examples/Logger/SimpleOrbits/plotSolution.py +++ b/examples/Logger/SimpleOrbits/plotSolution.py @@ -17,7 +17,7 @@ # ############################################################################## -from h5py import File +from h5py import HDF5File import numpy as np from glob import glob import matplotlib.pyplot as plt @@ -113,7 +113,7 @@ def doSnapshots(): for i, f in enumerate(filenames): # get the data from the file - f = File(f, "r") + f = HDF5File(f, "r") pos = f["PartType1/Coordinates"][:] pos -= center vel = f["PartType1/Velocities"][:] @@ -150,12 +150,12 @@ def doSnapshots(): plt.figure(fig_2.number) if diff_to_sol: p[:, :2] = p[:, :2] - y[:2, :].transpose() - plt.plot(p[:, 0], p[:, 1], ".", label="Snapshot") + plt.plot(p[:, 0], p[:, 1], "-", label="Snapshot", lw=1.) plt.figure(fig_3.number) if diff_to_sol: v[:, :2] = v[:, :2] - y[2:, :].transpose() - plt.plot(v[:, 0], v[:, 1], ".", label="Snapshot") + plt.plot(v[:, 0], v[:, 1], "-", label="Snapshot", lw=1.) def doStatistics(): @@ -256,12 +256,12 @@ def doLogger(): if diff_to_sol: p[:, :2] = p[:, :2] - y[:2, :].transpose() plt.figure(fig_2.number) - plt.plot(p[:, 0], p[:, 1], "--", label="Logger (Interpolation)") + plt.plot(p[:, 0], p[:, 1], ":r", label="Logger (Interpolation)") plt.figure(fig_3.number) if diff_to_sol: v[:, :2] = v[:, :2] - y[2:, :].transpose() - plt.plot(v[:, 0], v[:, 1], "--", label="Logger (Interpolation)") + plt.plot(v[:, 0], v[:, 1], ":r", label="Logger (Interpolation)") # do all the plots @@ -279,11 +279,15 @@ 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.png") +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 index 768c40df1..e1742c357 100644 --- a/examples/Logger/SimpleOrbits/run.sh +++ b/examples/Logger/SimpleOrbits/run.sh @@ -5,7 +5,7 @@ echo "Generating initial conditions for the Simple Orbits example..." python makeIC.py # Run SWIFT -../../swift --external-gravity --threads=4 simple_orbits.yml 2>&1 | tee output.log +../../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 index dca898c4a..448f146f3 100644 --- a/examples/Logger/SimpleOrbits/simple_orbits.yml +++ b/examples/Logger/SimpleOrbits/simple_orbits.yml @@ -9,19 +9,19 @@ InternalUnitSystem: # Parameters governing the time integration TimeIntegration: time_begin: 0. # The starting time of the simulation (in internal units). - time_end: 0.1 # The end 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: 1e-2 # The maximal 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: 4e-3 # Time difference between consecutive outputs (in internal units) + delta_time: 5e-3 # Time difference between consecutive outputs (in internal units) # Parameters governing the conserved quantities statistics Statistics: - delta_time: 1e-5 # Time between statistics output + delta_time: 1 # Time between statistics output # Parameters related to the initial conditions @@ -35,12 +35,12 @@ 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-5 # Dimensionless pre-factor for the time-step condition + timestep_mult: 1e-3 # Dimensionless pre-factor for the time-step condition # Parameters governing the logger snapshot system Logger: - delta_step: 2000 # Update the particle log every this many updates + 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 diff --git a/examples/Planetary/EarthImpact/earth_impact.yml b/examples/Planetary/EarthImpact/earth_impact.yml index b1df0a697..e98f7a072 100644 --- a/examples/Planetary/EarthImpact/earth_impact.yml +++ b/examples/Planetary/EarthImpact/earth_impact.yml @@ -55,4 +55,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/logger/examples/reader_example.py b/logger/examples/reader_example.py index b15446a6c..62cca38a9 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_particle.c b/logger/logger_particle.c index b07b9a570..c912c40dd 100644 --- a/logger/logger_particle.c +++ b/logger/logger_particle.c @@ -40,7 +40,7 @@ void logger_particle_print(const struct logger_particle *p) { message("Positions: (%g, %g, %g).", p->pos[0], p->pos[1], p->pos[2]); message("Velocities: (%g, %g, %g).", p->vel[0], p->vel[1], p->vel[2]); message("Accelerations: (%g, %g, %g).", p->acc[0], p->acc[1], p->acc[2]); - message("Entropy: %g.", p->entropy); + message("Entropy: %i.", p->entropy); message("Density: %g.", p->density); message("Type: %i.", p->type); } @@ -361,9 +361,14 @@ void logger_particle_interpolate(struct logger_particle *part_curr, part_curr->acc[i] += ftmp * scaling; } - /* interpolate scalars. */ - ftmp = (part_next->entropy - part_curr->entropy); - part_curr->entropy += ftmp * scaling; + /* /\* interpolate scalars. *\/ */ + /* 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_particle.h b/logger/logger_particle.h index 57cab1f6f..ede181a1a 100644 --- a/logger/logger_particle.h +++ b/logger/logger_particle.h @@ -62,7 +62,7 @@ struct logger_particle { float acc[3]; /* entropy. */ - float entropy; + int entropy; /* smoothing length. */ float h; diff --git a/logger/logger_python_wrapper.c b/logger/logger_python_wrapper.c index 1774d01b1..985234b9c 100644 --- a/logger/logger_python_wrapper.c +++ b/logger/logger_python_wrapper.c @@ -90,9 +90,6 @@ static PyObject *loadSnapshotAtTime(__attribute__((unused)) PyObject *self, 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; @@ -221,6 +218,14 @@ static PyObject *pyMoveForwardInTime(__attribute__((unused)) PyObject *self, 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; @@ -241,6 +246,7 @@ static PyObject *pyMoveForwardInTime(__attribute__((unused)) PyObject *self, /* 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)", diff --git a/logger/logger_reader.c b/logger/logger_reader.c index c9783b1e3..173aaff12 100644 --- a/logger/logger_reader.c +++ b/logger/logger_reader.c @@ -430,9 +430,13 @@ 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 when requesting time with offset %zi", - prev_offset, time_offset); + bzero(prev, sizeof(struct logger_particle)); + bzero(next, sizeof(struct logger_particle)); + return; + /* time_array_print(&reader->log.times); */ + /* 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 4c0c4b786..ca0a5b908 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 38056dd68..06299fc93 100644 --- a/src/logger.c +++ b/src/logger.c @@ -70,7 +70,7 @@ const struct mask_data logger_mask_data[logger_count_mask] = { /* Particle's acceleration. */ {3 * sizeof(float), 1 << logger_a, "accelerations"}, /* Particle's entropy. */ - {sizeof(float), 1 << logger_u, "entropy"}, + {sizeof(int), 1 << logger_u, "entropy"}, /* Particle's smoothing length. */ {sizeof(float), 1 << logger_h, "smoothing length"}, /* Particle's density. */ @@ -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); -- GitLab From 91150012d568b18241cc24be916488f2a58f0e56 Mon Sep 17 00:00:00 2001 From: loikki Date: Thu, 27 Feb 2020 07:38:19 +0100 Subject: [PATCH 09/12] Logger: sort the particles by ids when plotting the solution --- examples/Logger/SimpleOrbits/plotSolution.py | 81 +++++++++----------- 1 file changed, 38 insertions(+), 43 deletions(-) diff --git a/examples/Logger/SimpleOrbits/plotSolution.py b/examples/Logger/SimpleOrbits/plotSolution.py index ee06b6a5b..c49965176 100644 --- a/examples/Logger/SimpleOrbits/plotSolution.py +++ b/examples/Logger/SimpleOrbits/plotSolution.py @@ -17,7 +17,7 @@ # ############################################################################## -from h5py import HDF5File +from h5py import File as HDF5File import numpy as np from glob import glob import matplotlib.pyplot as plt @@ -114,28 +114,23 @@ def doSnapshots(): for i, f in enumerate(filenames): # get the data from the file f = HDF5File(f, "r") - pos = f["PartType1/Coordinates"][:] - pos -= center - vel = f["PartType1/Velocities"][:] 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"] - # loop over all the particles - # (avoids comparing different particles) - for j in range(makeIC.num_part): - # Get the current id - ind = ids == j - - # Compute the energy - r = np.sum(pos[ind, :]**2)**0.5 - v2 = np.sum(vel[ind, :]**2) - E[i, ind] = 0.5 * v2 - G * M / r + 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 - if j == id_focus: - p[i, :] = pos[ind, :] - v[i, :] = vel[ind, :] + # 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) @@ -199,31 +194,31 @@ def doLogger(): # Get the next particles interp = logger.moveForwardInTime( basename, parts, t, verbose) - rel_pos = interp["positions"] - center - vel = interp["velocities"] - - rel_pos_parts = parts["positions"] - center - vel_parts = parts["velocities"] - - for j in range(makeIC.num_part): - # Compute the interpolated variables - ind = interp["ids"] == j - r = np.sum(rel_pos[ind, :]**2)**0.5 - v2 = np.sum(vel[ind, :]**2) - E[i, ind] = 0.5 * v2 - G * M / r - if j == id_focus: - p[i, :] = rel_pos[ind, :] - v[i, :] = vel[ind, :] - - # Compute the variables of the last record - ind = parts["ids"] == j - r = np.sum(rel_pos_parts[ind, :]**2)**0.5 - v2 = np.sum(vel_parts[ind, :]**2) - E_parts[i, ind] = 0.5 * v2 - G * M / r - t_parts[i, ind] = parts["times"][ind] - if j == id_focus: - p_parts[i, :] = rel_pos_parts[ind, :] - v_parts[i, :] = vel_parts[ind, :] + 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) -- GitLab From 9ef9696a8e31ae974468f225350e52e23283b472 Mon Sep 17 00:00:00 2001 From: loikki Date: Wed, 25 Mar 2020 14:50:19 +0100 Subject: [PATCH 10/12] Logger example: Remove IVP --- examples/Logger/SimpleOrbits/plotSolution.py | 36 -------------------- 1 file changed, 36 deletions(-) diff --git a/examples/Logger/SimpleOrbits/plotSolution.py b/examples/Logger/SimpleOrbits/plotSolution.py index c49965176..b160e2496 100644 --- a/examples/Logger/SimpleOrbits/plotSolution.py +++ b/examples/Logger/SimpleOrbits/plotSolution.py @@ -21,7 +21,6 @@ from h5py import File as HDF5File import numpy as np from glob import glob import matplotlib.pyplot as plt -from scipy.integrate import solve_ivp import makeIC import sys sys.path.append("../../../logger/.libs/") @@ -32,9 +31,6 @@ plt.style.use("mpl_style") center = np.array([0.5 * makeIC.boxsize]*3) id_focus = 0 -# Plot the absolute positions (+vel) or the difference -# to a solution computed with a RK solver. -diff_to_sol = False # Defines the constants G = 1.189972e-04 # gravitational constant @@ -69,23 +65,6 @@ def gravity(t, y): return dy -def Evolve(y0, t): - """ - Compute the solution to the gravitational problem. - - Parameters - ---------- - - y0: np.array - Initial conditions [x, y, vx, vy] - - t: np.array - Time of the outputs. - """ - y = solve_ivp(gravity, (t[0], t[-1]), y0, t_eval=t) - return y.y - - def plotRelative(t, p, *args, **kwargs): """ Wrapper around the function plot from matplotlib, but @@ -136,20 +115,15 @@ def doSnapshots(): y0 = np.zeros(4) y0[:2] = p[0, :2] y0[2:] = v[0, :2] - y = Evolve(y0, t) # compute the plotting variables plt.figure(fig_1.number) plotRelative(t, E, ".", label="Snapshot") plt.figure(fig_2.number) - if diff_to_sol: - p[:, :2] = p[:, :2] - y[:2, :].transpose() plt.plot(p[:, 0], p[:, 1], "-", label="Snapshot", lw=1.) plt.figure(fig_3.number) - if diff_to_sol: - v[:, :2] = v[:, :2] - y[2:, :].transpose() plt.plot(v[:, 0], v[:, 1], "-", label="Snapshot", lw=1.) @@ -230,32 +204,22 @@ def doLogger(): y0[:2] = p[0, :2] y0[2:] = v[0, :2] t_parts, ind = np.unique(t_parts[:, 0], return_index=True) - y = Evolve(y0, t_parts) - if diff_to_sol: - p_parts = p_parts[ind, :2] - y[:2, :].transpose() # plot the solution plt.figure(fig_2.number) plt.plot(p_parts[:, 0], p_parts[:, 1], "x", label="Logger") plt.figure(fig_3.number) - if diff_to_sol: - v_parts = v_parts[ind, :2] - y[2:, :].transpose() 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] - y = Evolve(y0, times) - if diff_to_sol: - p[:, :2] = p[:, :2] - y[:2, :].transpose() plt.figure(fig_2.number) plt.plot(p[:, 0], p[:, 1], ":r", label="Logger (Interpolation)") plt.figure(fig_3.number) - if diff_to_sol: - v[:, :2] = v[:, :2] - y[2:, :].transpose() plt.plot(v[:, 0], v[:, 1], ":r", label="Logger (Interpolation)") -- GitLab From 09ac4847e6c9fc4efc270efdf37d5be2e481b4c4 Mon Sep 17 00:00:00 2001 From: loikki Date: Thu, 26 Mar 2020 08:08:47 +0100 Subject: [PATCH 11/12] Logger: fix entropy --- logger/logger_particle.c | 9 +++++---- logger/logger_particle.h | 2 +- logger/logger_reader.c | 11 ++++------- src/logger.c | 2 +- 4 files changed, 11 insertions(+), 13 deletions(-) diff --git a/logger/logger_particle.c b/logger/logger_particle.c index c912c40dd..9a9ab624c 100644 --- a/logger/logger_particle.c +++ b/logger/logger_particle.c @@ -40,7 +40,7 @@ void logger_particle_print(const struct logger_particle *p) { message("Positions: (%g, %g, %g).", p->pos[0], p->pos[1], p->pos[2]); message("Velocities: (%g, %g, %g).", p->vel[0], p->vel[1], p->vel[2]); message("Accelerations: (%g, %g, %g).", p->acc[0], p->acc[1], p->acc[2]); - message("Entropy: %i.", p->entropy); + message("Entropy: %g.", p->entropy); message("Density: %g.", p->density); message("Type: %i.", p->type); } @@ -361,9 +361,10 @@ void logger_particle_interpolate(struct logger_particle *part_curr, part_curr->acc[i] += ftmp * scaling; } - /* /\* interpolate scalars. *\/ */ - /* ftmp = (part_next->entropy - part_curr->entropy); */ - /* part_curr->entropy += ftmp * scaling; */ + /* interpolate scalars. */ + ftmp = (part_next->entropy - part_curr->entropy); + part_curr->entropy += ftmp * scaling; + ftmp = (part_next->h - part_curr->h); part_curr->h += ftmp * scaling; diff --git a/logger/logger_particle.h b/logger/logger_particle.h index ede181a1a..57cab1f6f 100644 --- a/logger/logger_particle.h +++ b/logger/logger_particle.h @@ -62,7 +62,7 @@ struct logger_particle { float acc[3]; /* entropy. */ - int entropy; + float entropy; /* smoothing length. */ float h; diff --git a/logger/logger_reader.c b/logger/logger_reader.c index 173aaff12..fc93ca7bc 100644 --- a/logger/logger_reader.c +++ b/logger/logger_reader.c @@ -430,13 +430,10 @@ void logger_reader_get_next_particle(struct logger_reader *reader, /* Are we at the end of the file? */ if (next_offset == 0) { - bzero(prev, sizeof(struct logger_particle)); - bzero(next, sizeof(struct logger_particle)); - return; - /* time_array_print(&reader->log.times); */ - /* 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); */ + time_array_print(&reader->log.times); + 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/logger.c b/src/logger.c index 06299fc93..c970385e2 100644 --- a/src/logger.c +++ b/src/logger.c @@ -70,7 +70,7 @@ const struct mask_data logger_mask_data[logger_count_mask] = { /* Particle's acceleration. */ {3 * sizeof(float), 1 << logger_a, "accelerations"}, /* Particle's entropy. */ - {sizeof(int), 1 << logger_u, "entropy"}, + {sizeof(float), 1 << logger_u, "entropy"}, /* Particle's smoothing length. */ {sizeof(float), 1 << logger_h, "smoothing length"}, /* Particle's density. */ -- GitLab From e3bba1805181cd2bbbe3a86b9a76eb280e252a15 Mon Sep 17 00:00:00 2001 From: loikki Date: Thu, 26 Mar 2020 08:09:16 +0100 Subject: [PATCH 12/12] Format --- logger/logger_particle.c | 24 +++++++++++++----------- logger/logger_python_wrapper.c | 24 ++++++++++-------------- logger/logger_reader.c | 8 +++++--- 3 files changed, 28 insertions(+), 28 deletions(-) diff --git a/logger/logger_particle.c b/logger/logger_particle.c index 9a9ab624c..07cddd89b 100644 --- a/logger/logger_particle.c +++ b/logger/logger_particle.c @@ -226,9 +226,9 @@ size_t logger_particle_read(struct logger_particle *part, * * @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) { +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 */ @@ -260,7 +260,8 @@ double logger_particle_quintic_hermite_spline( 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; + x += (6. * x1 - 6. * x0 - 3. * v0_dt - 3. * v1_dt + a1_dt2 - a0_dt2) * + t_t0_3 * t_t1_2 / dt5; return x; } @@ -278,9 +279,9 @@ double logger_particle_quintic_hermite_spline( * * @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) { +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 */ @@ -348,13 +349,14 @@ void logger_particle_interpolate(struct logger_particle *part_curr, 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); + 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); + 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]); diff --git a/logger/logger_python_wrapper.c b/logger/logger_python_wrapper.c index 985234b9c..48973daf9 100644 --- a/logger/logger_python_wrapper.c +++ b/logger/logger_python_wrapper.c @@ -174,7 +174,6 @@ static PyObject *pyReverseOffset(__attribute__((unused)) PyObject *self, return Py_BuildValue(""); } - /** * @brief Move forward in time an array of particles. * @@ -196,9 +195,9 @@ static PyObject *pyMoveForwardInTime(__attribute__((unused)) PyObject *self, int new_array = 1; /* parse the arguments. */ - if (!PyArg_ParseTuple(args, "sOd|ii", &filename, - &parts, &time, &verbose, &new_array)) return NULL; - + if (!PyArg_ParseTuple(args, "sOd|ii", &filename, &parts, &time, &verbose, + &new_array)) + return NULL; /* Check parts */ if (!PyArray_Check(parts)) { @@ -216,8 +215,8 @@ static PyObject *pyMoveForwardInTime(__attribute__((unused)) PyObject *self, /* Create the interpolated array. */ PyArrayObject *interp = NULL; if (new_array) { - interp = (PyArrayObject *) PyArray_NewLikeArray( - parts, NPY_ANYORDER, NULL, 0); + interp = + (PyArrayObject *)PyArray_NewLikeArray(parts, NPY_ANYORDER, NULL, 0); /* Check if the allocation was fine */ if (interp == NULL) { @@ -226,8 +225,7 @@ static PyObject *pyMoveForwardInTime(__attribute__((unused)) PyObject *self, /* Reference stolen in PyArray_NewLikeArray => incref */ Py_INCREF(PyArray_DESCR(parts)); - } - else { + } else { interp = parts; // We return it, therefore one more reference exists. Py_INCREF(interp); @@ -242,15 +240,14 @@ static PyObject *pyMoveForwardInTime(__attribute__((unused)) PyObject *self, /* Loop over all the particles */ size_t N = PyArray_DIM(parts, 0); - for(size_t i = 0; i < N; i++) { + 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); + 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); @@ -306,7 +303,7 @@ static PyMethodDef libloggerMethods[] = { "times: tuple\n" " time min, time max\n"}, {"moveForwardInTime", pyMoveForwardInTime, METH_VARARGS, - "Move the particles forward in time.\n\n" + "Move the particles forward in time.\n\n" "Parameters\n" "----------\n\n" "basename: str\n" @@ -322,8 +319,7 @@ static PyMethodDef libloggerMethods[] = { "Returns\n" "-------\n\n" "parts: np.array\n" - " The particles at the requested time.\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 fc93ca7bc..ce7799c65 100644 --- a/logger/logger_reader.c +++ b/logger/logger_reader.c @@ -431,9 +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 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); + 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; -- GitLab