gitlab is upgraded to version 13, please report any issues and enjoy

Commit c78c1aa2 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'logger_reader' into 'master'

Logger reader

See merge request !1127
parents dd2bff29 bcceec66
......@@ -43,22 +43,28 @@ lib_LTLIBRARIES = liblogger.la
# lib_LTLIBRARIES += liblogger_mpi.la
# endif
GRAVITY_SRC = gravity/MultiSoftening/logger_gravity.c
STARS_SRC = stars/Default/logger_stars.c
HYDRO_SRC = hydro/Gadget2/logger_hydro.c
# subdirectories
SUBDIRS = tests
# List required headers
include_HEADERS = logger_header.h logger_loader_io.h logger_particle.h logger_time.h logger_tools.h logger_reader.h \
logger_logfile.h logger_index.h quick_sort.h
logger_logfile.h logger_index.h quick_sort.h logger_python_tools.h
# Common source files
AM_SOURCES = logger_header.c logger_loader_io.c logger_particle.c logger_time.c logger_tools.c logger_reader.c \
logger_logfile.c logger_index.c quick_sort.c
logger_logfile.c logger_index.c quick_sort.c $(GRAVITY_SRC) $(STARS_SRC) $(HYDRO_SRC)
if HAVEPYTHON
AM_SOURCES += logger_python_wrapper.c
endif
# Include files for distribution, not installation.
nobase_noinst_HEADERS =
nobase_noinst_HEADERS = logger_hydro.h hydro/Gadget2/logger_hydro.h \
logger_stars.h stars/Default/logger_stars.h \
logger_gravity.h gravity/MultiSoftening/logger_gravity.h
# Sources and flags for regular library
liblogger_la_SOURCES = $(AM_SOURCES)
......
......@@ -6,34 +6,33 @@ Example: ./reader_example.py ../../examples/SedovBlast_3D/index 0.1
import sys
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
sys.path.append("../.libs/")
import liblogger as logger
def plot3D(pos):
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
ax.plot(pos[:, 0], pos[:, 1], pos[:, 2], ".", markersize=0.2)
ax.plot(pos[:, 0], pos[:, 1], pos[:, 2], ".", markersize=0.1)
def plot2D():
def plot2D(pos, entropies):
plt.figure()
center = np.array([0.5]*3)
r2 = np.sum((pos - center)**2, axis=1)
# plot entropy vs distance
plt.plot(np.sqrt(r2), data["entropies"], '.',
plt.plot(np.sqrt(r2), entropies, '.',
markersize=0.2)
plt.xlim(0., 0.5)
plt.ylim(-1, 50)
plt.xlabel("Radius")
plt.ylabel("Entropy")
basename = "../../examples/HydroTests/SedovBlast_3D/index"
time = 0.05
basename = "../../examples/HydroTests/SedovBlast_3D/index_0000"
time = 0.03
if len(sys.argv) >= 2:
basename = sys.argv[1]
else:
......@@ -47,14 +46,13 @@ if len(sys.argv) > 3:
print("basename: %s" % basename)
print("time: %g" % time)
# read dump
# read the logger
t = logger.getTimeLimits(basename)
data = logger.loadSnapshotAtTime(basename, time, 2)
print("The data contains the following elements:")
print(data.dtype.names)
pos = data["positions"]
pos, ent = logger.get_particle_data(basename, ["Coordinates", "Entropies"],
time)
print("Min/Max of the position:", pos.min(), pos.max())
print("Min/Max of the entropy:", ent.min(), ent.max())
plot3D(pos)
plot2D(pos, ent)
plt.show()
/*******************************************************************************
* 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/>.
*
******************************************************************************/
#include "logger_gravity.h"
/* Define the size of all the fields. */
#define member_size(type, member) sizeof(((type *)0)->member)
const int gravity_logger_field_size[gravity_logger_field_count] = {
member_size(struct gpart, x), // coordinates
member_size(struct gpart, v_full), // velocities
member_size(struct gpart, a_grav), // accelerations
member_size(struct gpart, mass), // massses
member_size(struct gpart, id_or_neg_offset), // IDs
};
int gravity_logger_local_to_global[gravity_logger_field_count];
/*******************************************************************************
* 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/>.
*
******************************************************************************/
#ifndef SWIFT_MULTISOFTENING_LOGGER_GRAVITY_H
#define SWIFT_MULTISOFTENING_LOGGER_GRAVITY_H
#include "../config.h"
/* local includes */
#include "gravity_logger.h"
#include "logger_loader_io.h"
#include "logger_python_tools.h"
/* Index of the mask in the header mask array */
extern int gravity_logger_local_to_global[gravity_logger_field_count];
/* Size for each mask */
extern const int gravity_logger_field_size[gravity_logger_field_count];
/**
* @brief Create the link between the fields and their derivatives.
*
* @param head The #header.
*/
__attribute__((always_inline)) INLINE static void
gravity_logger_reader_link_derivatives(struct header *head) {
/* Set the first and second derivatives */
const int pos_id =
gravity_logger_local_to_global[gravity_logger_field_coordinates];
const int vel_id =
gravity_logger_local_to_global[gravity_logger_field_velocities];
const int acc_id =
gravity_logger_local_to_global[gravity_logger_field_accelerations];
/* Coordinates */
header_set_first_derivative(head, pos_id, vel_id);
header_set_second_derivative(head, pos_id, acc_id);
/* Velocities */
header_set_first_derivative(head, vel_id, acc_id);
}
/**
* @brief Interpolate a field of the particle at the given time.
* Here we use a linear interpolation for most of the fields.
* For the position (velocity), we use a quintic (cubic) hermite interpolation
* based on the positions, velocities and accelerations at the time of the two
* particles.
*
* @param before Pointer to the #logger_field at a time < t.
* @param after Pointer to the #logger_field at a time > t.
* @param otuput Pointer to the output value.
* @param t_before Time of field_before (< t).
* @param t_after Time of field_after (> t).
* @param t Requested time.
* @param field The field to reconstruct (follows the order of
* #gravity_logger_fields).
*/
__attribute__((always_inline)) INLINE static void
gravity_logger_interpolate_field(const double t_before,
const struct logger_field *restrict before,
const double t_after,
const struct logger_field *restrict after,
void *restrict output, const double t,
const int field) {
#ifdef SWIFT_DEBUG_CHECKS
/* Check the times */
if (t_before > t || t_after < t) {
error(
"The times for the interpolation are not correct"
" %g < %g < %g.",
t_before, t, t_after);
}
#endif
/* Compute the interpolation scaling. */
const double wa = (t - t_before) / (t_after - t_before);
const double wb = 1. - wa;
switch (field) {
case gravity_logger_field_coordinates:
for (int i = 0; i < 3; i++) {
double *x = (double *)output;
const double *x_bef = (double *)before->field;
const float *v_bef = (float *)before->first_deriv;
const float *a_bef = (float *)before->second_deriv;
const double *x_aft = (double *)after->field;
const float *v_aft = (float *)after->first_deriv;
const float *a_aft = (float *)after->second_deriv;
/* Use quintic hermite spline. */
if (v_bef && v_aft && a_bef && a_aft) {
x[i] = logger_tools_quintic_hermite_spline(
t_before, x_bef[i], v_bef[i], a_bef[i], t_after, x_aft[i],
v_aft[i], a_aft[i], t);
}
/* Use cubic hermite spline. */
else if (v_bef && v_aft) {
x[i] = logger_tools_cubic_hermite_spline(
t_before, x_bef[i], v_bef[i], t_after, x_aft[i], v_aft[i], t);
}
/* Use linear interpolation. */
else {
x[i] = wa * x_aft[i] + wb * x_bef[i];
}
}
break;
case gravity_logger_field_velocities:
for (int i = 0; i < 3; i++) {
float *v = (float *)output;
const float *v_bef = (float *)before->field;
const float *a_bef = (float *)before->first_deriv;
const float *v_aft = (float *)after->field;
const float *a_aft = (float *)after->first_deriv;
/* Use a cubic hermite spline. */
if (a_bef && a_aft) {
v[i] = logger_tools_cubic_hermite_spline(
t_before, v_bef[i], a_bef[i], t_after, v_aft[i], a_aft[i], t);
}
/* Use linear interpolation. */
else {
v[i] = wa * v_aft[i] + wb * v_bef[i];
}
}
break;
case gravity_logger_field_accelerations:
/* interpolate vectors. */
for (int i = 0; i < 3; i++) {
float *a = (float *)output;
const float *a_bef = (float *)before->field;
const float *a_aft = (float *)after->field;
a[i] = wa * a_aft[i] + wb * a_bef[i];
}
break;
case gravity_logger_field_masses:
((float *)output)[0] =
wa * ((float *)after->field)[0] + wb * ((float *)before->field)[0];
break;
case gravity_logger_field_particle_ids:
if (*(long long *)after->field != *(long long *)before->field) {
error("Interpolating different particles");
}
*(long long *)output = *(long long *)after->field;
break;
default:
error("Not implemented");
}
}
#ifdef HAVE_PYTHON
__attribute__((always_inline)) INLINE static void
gravity_logger_generate_python(struct logger_python_field *fields) {
fields[gravity_logger_field_coordinates] =
logger_loader_python_field(/* Dimension */ 3, NPY_DOUBLE);
fields[gravity_logger_field_velocities] =
logger_loader_python_field(/* Dimension */ 3, NPY_FLOAT32);
fields[gravity_logger_field_accelerations] =
logger_loader_python_field(/* Dimension */ 3, NPY_FLOAT32);
fields[gravity_logger_field_masses] =
logger_loader_python_field(/* Dimension */ 1, NPY_FLOAT32);
fields[gravity_logger_field_particle_ids] =
logger_loader_python_field(/* Dimension */ 1, NPY_LONGLONG);
}
#endif // HAVE_PYTHON
#endif // SWIFT_MULTISOFTENING_LOGGER_GRAVITY_H
/*******************************************************************************
* 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/>.
*
******************************************************************************/
#include "logger_hydro.h"
int hydro_logger_local_to_global[hydro_logger_field_count];
/* Define the size of all the fields. */
#define member_size(type, member) sizeof(((type *)0)->member)
const int hydro_logger_field_size[hydro_logger_field_count] = {
member_size(struct part, x), // coordinates
member_size(struct part, v), // velocities
member_size(struct part, a_hydro), // accelerations
member_size(struct part, mass), // massses
member_size(struct part, h), // Smoothing Length
member_size(struct part, entropy), // Entropy
member_size(struct part, id), // IDs
member_size(struct part, rho), // density
};
/*******************************************************************************
* 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/>.
*
******************************************************************************/
#ifndef SWIFT_GADGET2_LOGGER_HYDRO_H
#define SWIFT_GADGET2_LOGGER_HYDRO_H
#include "../config.h"
/* local includes */
#include "hydro_logger.h"
#include "logger_loader_io.h"
#include "logger_python_tools.h"
/* Index of the mask in the header mask array */
extern int hydro_logger_local_to_global[hydro_logger_field_count];
/* Size for each mask */
extern const int hydro_logger_field_size[hydro_logger_field_count];
/**
* @brief Create the link between the fields and their derivatives.
*
* @param head The #header.
*/
__attribute__((always_inline)) INLINE static void
hydro_logger_reader_link_derivatives(struct header *head) {
/* Set the first and second derivatives */
const int pos_id =
hydro_logger_local_to_global[hydro_logger_field_coordinates];
const int vel_id =
hydro_logger_local_to_global[hydro_logger_field_velocities];
const int acc_id =
hydro_logger_local_to_global[hydro_logger_field_accelerations];
/* Coordinates */
header_set_first_derivative(head, pos_id, vel_id);
header_set_second_derivative(head, pos_id, acc_id);
/* Velocities */
header_set_first_derivative(head, vel_id, acc_id);
}
/**
* @brief Interpolate a field of the particle at the given time.
* Here we use a linear interpolation for most of the fields.
* For the position (velocity), we use a quintic (cubic) hermite interpolation
* based on the positions, velocities and accelerations at the time of the two
* particles.
*
* @param before Pointer to the #logger_field at a time < t.
* @param after Pointer to the #logger_field at a time > t.
* @param otuput Pointer to the output value.
* @param t_before Time of field_before (< t).
* @param t_after Time of field_after (> t).
* @param t Requested time.
* @param field The field to reconstruct (follows the order of
* #hydro_logger_fields).
*/
__attribute__((always_inline)) INLINE static void
hydro_logger_interpolate_field(const double t_before,
const struct logger_field *restrict before,
const double t_after,
const struct logger_field *restrict after,
void *restrict output, const double t,
const int field) {
#ifdef SWIFT_DEBUG_CHECKS
/* Check the times */
if (t_before > t || t_after < t) {
error(
"The times for the interpolation are not correct"
" %g < %g < %g.",
t_before, t, t_after);
}
#endif
/* Compute the interpolation scaling. */
const double wa = (t - t_before) / (t_after - t_before);
const double wb = 1. - wa;
switch (field) {
/* Do the position */
case hydro_logger_field_coordinates:
for (int i = 0; i < 3; i++) {
double *x = (double *)output;
const double *x_bef = (double *)before->field;
const float *v_bef = (float *)before->first_deriv;
const float *a_bef = (float *)before->second_deriv;
const double *x_aft = (double *)after->field;
const float *v_aft = (float *)after->first_deriv;
const float *a_aft = (float *)after->second_deriv;
/* Use quintic hermite spline. */
if (v_bef && v_aft && a_bef && a_aft) {
x[i] = logger_tools_quintic_hermite_spline(
t_before, x_bef[i], v_bef[i], a_bef[i], t_after, x_aft[i],
v_aft[i], a_aft[i], t);
}
/* Use cubic hermite spline. */
else if (v_bef && v_aft) {
x[i] = logger_tools_cubic_hermite_spline(
t_before, x_bef[i], v_bef[i], t_after, x_aft[i], v_aft[i], t);
}
/* Use linear interpolation. */
else {
x[i] = wa * x_aft[i] + wb * x_bef[i];
}
}
break;
/* Do the velocity */
case hydro_logger_field_velocities:
for (int i = 0; i < 3; i++) {
float *v = (float *)output;
const float *v_bef = (float *)before->field;
const float *a_bef = (float *)before->first_deriv;
const float *v_aft = (float *)after->field;
const float *a_aft = (float *)after->first_deriv;
/* Use a cubic hermite spline. */
if (a_bef && a_aft) {
v[i] = logger_tools_cubic_hermite_spline(
t_before, v_bef[i], a_bef[i], t_after, v_aft[i], a_aft[i], t);
}
/* Use linear interpolation. */
else {
v[i] = wa * v_aft[i] + wb * v_bef[i];
}
}
break;
case hydro_logger_field_accelerations:
/* interpolate vectors. */
for (int i = 0; i < 3; i++) {
float *a = (float *)output;
const float *a_bef = (float *)before->field;
const float *a_aft = (float *)after->field;
a[i] = wa * a_aft[i] + wb * a_bef[i];
}
break;
/* Do the linear interpolation of float. */
case hydro_logger_field_smoothing_lengths:
case hydro_logger_field_entropies:
case hydro_logger_field_densities:
case hydro_logger_field_masses:
((float *)output)[0] =
wa * ((float *)after->field)[0] + wb * ((float *)before->field)[0];
break;
/* Check the ids */
case hydro_logger_field_particle_ids:
if (*(long long *)after->field != *(long long *)before->field) {
error("Interpolating different particles");
}
*(long long *)output = *(long long *)after->field;
break;
default:
error("Not implemented");
}
}
#ifdef HAVE_PYTHON
__attribute__((always_inline)) INLINE static void hydro_logger_generate_python(
struct logger_python_field *fields) {
fields[hydro_logger_field_coordinates] =
logger_loader_python_field(/* Dimension */ 3, NPY_DOUBLE);
fields[hydro_logger_field_velocities] =
logger_loader_python_field(/* Dimension */ 3, NPY_FLOAT32);
fields[hydro_logger_field_accelerations] =
logger_loader_python_field(/* Dimension */ 3, NPY_FLOAT32);
fields[hydro_logger_field_masses] =
logger_loader_python_field(/* Dimension */ 1, NPY_FLOAT32);
fields[hydro_logger_field_smoothing_lengths] =
logger_loader_python_field(/* Dimension */ 1, NPY_FLOAT32);
fields[hydro_logger_field_entropies] =
logger_loader_python_field(/* Dimension */ 1, NPY_FLOAT32);
fields[hydro_logger_field_particle_ids] =
logger_loader_python_field(/* Dimension */ 1, NPY_LONGLONG);
fields[hydro_logger_field_densities] =
logger_loader_python_field(/* Dimension */ 1, NPY_FLOAT32);
}
#endif // HAVE_PYTHON
#endif /* SWIFT_GADGET2_LOGGER_HYDRO_H */
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (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/>.
*
******************************************************************************/
#ifndef SWIFT_LOGGER_GRAVITY_H
#define SWIFT_LOGGER_GRAVITY_H
/* Config parameters. */
#include "../config.h"
/* Import the right functions */
#if defined(DEFAULT_GRAVITY)
#error TODO
#include "./gravity/Default/logger_gravity.h"
#elif defined(POTENTIAL_GRAVITY)
#error TODO
#include "./gravity/Potential/logger_gravity.h"
#elif defined(MULTI_SOFTENING_GRAVITY)
#include "./gravity/MultiSoftening/logger_gravity.h"
#else
#error "Invalid choice of gravity variant"
#endif
#endif