Skip to content
Snippets Groups Projects
Commit 087254a4 authored by Loic Hausammann's avatar Loic Hausammann
Browse files

Initial commit for logger loader

parent 2cf116fc
No related branches found
No related tags found
1 merge request!685Logger loader
......@@ -23,6 +23,9 @@ SUBDIRS = src argparse examples doc tests tools
if HAVEEAGLECOOLING
SUBDIRS += examples/Cooling/CoolingRates
endif
if HAVELOGGER
SUBDIRS += logger
endif
# Non-standard files that should be part of the distribution.
EXTRA_DIST = INSTALL.swift .clang-format format.sh
......@@ -86,6 +86,7 @@ AC_ARG_ENABLE([logger],
if test "$with_logger" = "yes"; then
AC_DEFINE([WITH_LOGGER], 1, [logger enabled])
fi
AM_CONDITIONAL([HAVELOGGER],[test $with_logger = "yes"])
# Interprocedural optimization support. Needs special handling for linking and
# archiving as well as compilation with Intels, needs to be done before
......@@ -996,6 +997,50 @@ fi
AC_SUBST([TBBMALLOC_LIBS])
AM_CONDITIONAL([HAVETBBMALLOC],[test -n "$TBBMALLOC_LIBS"])
# Check for python.
have_python="no"
AC_ARG_WITH([python],
[AS_HELP_STRING([--with-python=PATH],
[root directory where python is installed @<:@yes/no@:>@]
)],
[with_python="$withval"],
[with_python="no"]
)
if test "x$with_python" != "xno"; then
AC_CACHE_CHECK(
"python version",
ac_cv_ver_python,
[ac_cv_ver_python=`python -c 'import sys;print(sys.version[[:3]])' 2> /dev/null`])
if test "x$with_python" != "xyes" -a "x$with_python" != "x"; then
PYTHON_LIBS=""
PYTHON_INCS="-I$with_python/include/python$ac_cv_ver_python"
else
PYTHON_LIBS=""
PYTHON_INCS=""
fi
have_python="yes"
AC_CHECK_PROGS(
[PYTHON_BIN],
[python$ac_cv_ver_python],
[AC_MSG_ERROR(Cannot find python binary!)],
[$with_python/bin]
)
AC_CHECK_LIB(
[python${ac_cv_ver_python}m],
[PyArg_ParseTuple],
[AC_DEFINE([HAVE_PYTHON],1,[The python library appears to be present.])
],
[AC_MSG_ERROR(Cannot find python library!)]
)
fi
AC_SUBST([PYTHON_LIBS])
AC_SUBST([PYTHON_INCS])
AM_CONDITIONAL([HAVEPYTHON],[test -n "$PYTHON_LIBS"])
# Check for HDF5. This is required.
AX_LIB_HDF5
if test "$with_hdf5" != "yes"; then
......@@ -1991,7 +2036,7 @@ AM_CONDITIONAL([HAVEEAGLEFEEDBACK], [test $with_feedback = "EAGLE"])
# Handle .in files.
AC_CONFIG_FILES([Makefile src/Makefile examples/Makefile examples/Cooling/CoolingRates/Makefile doc/Makefile doc/Doxyfile tests/Makefile])
AC_CONFIG_FILES([argparse/Makefile tools/Makefile])
AC_CONFIG_FILES([argparse/Makefile tools/Makefile logger/Makefile])
AC_CONFIG_FILES([tests/testReading.sh], [chmod +x tests/testReading.sh])
AC_CONFIG_FILES([tests/testActivePair.sh], [chmod +x tests/testActivePair.sh])
AC_CONFIG_FILES([tests/test27cells.sh], [chmod +x tests/test27cells.sh])
......@@ -2046,6 +2091,7 @@ AC_MSG_RESULT([
VELOCIraptor enabled : $have_velociraptor
Particle Logger : $with_logger
FoF activated: : $enable_fof
Python enabled : $have_python
Hydro scheme : $with_hydro
Dimensionality : $with_dimension
......
# This file is part of SWIFT.
# Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
# Matthieu Schaller (matthieu.schaller@durham.ac.uk).
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU 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 General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
# Add the non-standard paths to the included library headers
AM_CFLAGS = $(PYTHON_INCS)
# Assign a "safe" version number
AM_LDFLAGS = -version-info 0:0:0
# The git command, if available.
GIT_CMD = @GIT_CMD@
# Additional dependencies for shared libraries.
EXTRA_LIBS = $(PROFILER_LIBS) $(TCMALLOC_LIBS) $(JEMALLOC_LIBS) $(TBBMALLOC_LIBS) $(PYTHON_LIBS)
# MPI libraries.
# MPI_LIBS = $(MPI_THREAD_LIBS)
# MPI_FLAGS = -DWITH_MPI
# Build the libswiftsim library
lib_LTLIBRARIES = libswiftlogger.la
# Build a MPI-enabled version too?
# if HAVEMPI
# lib_LTLIBRARIES += libswiftlogger_mpi.la
# endif
# List required headers
include_HEADERS = header.h io.h particle.h timeline.h tools.h
# Common source files
AM_SOURCES = header.c io.c logger_loader.c particle.c \
timeline.c tools.c
# Include files for distribution, not installation.
nobase_noinst_HEADERS =
# Sources and flags for regular library
libswiftlogger_la_SOURCES = $(AM_SOURCES)
libswiftlogger_la_CFLAGS = $(AM_CFLAGS)
libswiftlogger_la_LDFLAGS = $(AM_LDFLAGS) $(EXTRA_LIBS)
# Sources and flags for MPI library
# libswiftlogger_mpi_la_SOURCES = $(AM_SOURCES)
# libswiftlogger_mpi_la_CFLAGS = $(AM_CFLAGS) $(MPI_FLAGS)
# libswiftlogger_mpi_la_LDFLAGS = $(AM_LDFLAGS) $(MPI_LIBS) $(EXTRA_LIBS)
# libswiftlogger_mpi_la_SHORTNAME = mpi
# libswiftlogger_mpi_la_LIBADD =
# Versioning. If any sources change then update the version_string.h file with
# the current git revision and package version.
# May have a checkout without a version_string.h file and no git command (tar/zip
# download), allow that, but make sure we know it.
version_string.h: ../src/version_string.h.in $(AM_SOURCES) $(include_HEADERS) $(noinst_HEADERS)
if test "X$(GIT_CMD)" != "X"; then \
GIT_REVISION=`$(GIT_CMD) describe --abbrev=8 --always --tags --dirty`; \
GIT_BRANCH=`$(GIT_CMD) branch | sed -n 's/^\* \(.*\)/\1/p'`; \
GIT_DATE=`$(GIT_CMD) log -1 --format=%ci`; \
sed -e "s,@PACKAGE_VERSION\@,$(PACKAGE_VERSION)," \
-e "s,@GIT_REVISION\@,$${GIT_REVISION}," \
-e "s|@GIT_BRANCH\@|$${GIT_BRANCH}|" \
-e "s|@GIT_DATE\@|$${GIT_DATE}|" \
-e "s|@SWIFT_CFLAGS\@|$(CFLAGS)|" $< > version_string.h; \
else \
if test ! -f version_string.h; then \
sed -e "s,@PACKAGE_VERSION\@,$(PACKAGE_VERSION)," \
-e "s,@GIT_REVISION\@,unknown," \
-e "s,@GIT_BRANCH\@,unknown," \
-e "s,@GIT_DATE\@,unknown," \
-e "s|@SWIFT_CFLAGS\@|$(CFLAGS)|" $< > version_string.h; \
fi; \
fi
# Make sure version_string.h is built first.
BUILT_SOURCES = version_string.h
# And distribute the built files.
EXTRA_DIST = version_string.h ../src/version_string.h.in
#include "header.h"
#include "io.h"
#include "tools.h"
#include <Python.h>
#include <stdio.h>
#include <stdlib.h>
void header_print(const struct header *h) {
#ifdef SWIFT_DEBUG_CHECKS
printf("Debug checks enabled\n");
#endif
printf("Version: %s\n", h->version);
printf("First Offset: %lu\n", h->offset_first);
char direction[20];
if (h->forward_offset)
strcpy(direction, "Forward");
else
strcpy(direction, "Backward");
printf("Offset direction: %s\n", direction);
printf("Number masks: %lu\n", h->nber_mask);
for (size_t i = 0; i < h->nber_mask; i++) {
printf("\tMask: %s\n", h->masks_name[i]);
printf("\tValue: %lu\n", h->masks[i]);
printf("\tSize: %lu\n", h->masks_size[i]);
printf("\n");
}
/* mask contains... TODO */
};
void header_free(struct header *h) {
for (size_t i = 0; i < h->nber_mask; i++) {
free(h->masks_name[i]);
}
free(h->masks_name);
free(h->masks);
free(h->masks_size);
};
int header_is_present_and_get_index(const struct header *h, const char *field,
size_t *ind) {
for (size_t i = 0; i < h->nber_mask; i++) {
if (strcmp(h->masks_name[i], field) == 0) {
if (ind != NULL) {
*ind = i;
}
return 1;
}
}
return 0;
};
int header_is_present(const struct header *h, const char *field) {
return header_is_present_and_get_index(h, field, NULL);
};
int header_change_offset_direction(struct header *h, void *map) {
h->forward_offset = !h->forward_offset;
size_t offset = LOGGER_VERSION_SIZE;
return io_write_data(map, LOGGER_NBER_SIZE, &h->forward_offset,
&offset);
}
int header_read(struct header *h, void *map) {
size_t offset = 0;
/* read version */
io_read_data(map, LOGGER_VERSION_SIZE, &h->version, &offset);
/* read offset direction */
h->forward_offset = 0;
io_read_data(map, LOGGER_NBER_SIZE, &h->forward_offset, &offset);
if (h->forward_offset != 0 && h->forward_offset != 1)
error(EIO, "Non boolean value for the offset direction (%i)", h->forward_offset);
/* read offset to first data */
h->offset_first = 0;
io_read_data(map, LOGGER_OFFSET_SIZE, &h->offset_first, &offset);
/* read name size */
h->name = 0;
io_read_data(map, LOGGER_NBER_SIZE, &h->name, &offset);
/* check if value defined in this file is large enough */
if (STRING_SIZE < h->name) {
error(EOVERFLOW, "Name too large in dump file");
}
/* read number of masks */
h->nber_mask = 0;
io_read_data(map, LOGGER_NBER_SIZE, &h->nber_mask, &offset);
/* allocate memory */
h->masks = malloc(sizeof(size_t) * h->nber_mask);
h->masks_name = malloc(sizeof(void *) * h->nber_mask);
h->masks_size = malloc(sizeof(size_t) * h->nber_mask);
/* loop over all masks */
for (size_t i = 0; i < h->nber_mask; i++) {
/* read mask name */
h->masks_name[i] = malloc(h->name);
io_read_data(map, h->name, h->masks_name[i], &offset);
/* get mask value */
h->masks[i] = 1 << i;
/* read mask data size */
h->masks_size[i] = 0;
io_read_data(map, LOGGER_NBER_SIZE, &h->masks_size[i], &offset);
}
if (offset != h->offset_first) {
#ifdef SWIFT_DEBUG_CHECKS
header_print(h);
#endif
error(EIO, "Wrong header size (in header %li, current %li)",
h->offset_first, offset);
}
return 0;
};
size_t header_get_mask_size(const struct header *h, const size_t mask) {
size_t count = 0;
for (size_t i = 0; i < h->nber_mask; i++) {
if (mask & h->masks[i]) {
count += h->masks_size[i];
}
}
return count;
}
#ifndef __LOGGER_HEADER_H__
#define __LOGGER_HEADER_H__
#include "tools.h"
#include <Python.h>
#include <stdio.h>
#include <stdlib.h>
#define LOGGER_VERSION_SIZE 20
#define LOGGER_NBER_SIZE 4
#define LOGGER_OFFSET_SIZE 7
#define LOGGER_MASK_SIZE 1
struct header {
/* Logger version */
char version[STRING_SIZE];
/* offset of the first header */
size_t offset_first;
/* Number of bytes for names */
size_t name;
/* number of masks */
size_t nber_mask;
/* list of masks */
size_t *masks;
/* list of mask name */
char **masks_name;
/* size of data in a mask */
size_t *masks_size;
/* offset direction */
int forward_offset;
};
/**
* @brief print a header struct
*/
void header_print(const struct header *h);
/**
* @brief free allocated memory
*/
void header_free(struct header *h);
/**
* @brief check if field is present in header
*
* @param field name of the requested field
* @param ind (return value) indice of the requested field
*/
int header_is_present_and_get_index(const struct header *h, const char *field,
size_t *ind);
/**
* @brief check if field is present in header
*
* @param field name of the requested field
*/
int header_is_present(const struct header *h, const char *field);
/**
* @brief read the logger header
*
* @param h out: header
* @param map file mapping
*/
int header_read(struct header *h, void *map);
/**
* @brief count number of bits in a given mask
*
* @param h #header file structure
* @param mask mask to compute
*
* @return number of bits in mask
*/
size_t header_get_mask_size(const struct header *h, const size_t mask);
/**
* @brief Inverse the offset direction
*
* @param h #header file structure
* @param map file mapping
*
* @return error code
*/
int header_change_offset_direction(struct header *h, void *map);
#endif // __LOGGER_HEADER_H__
#include "io.h"
#include "header.h"
#include "tools.h"
/* file size */
#include <sys/stat.h>
/* mapping */
#include <sys/mman.h>
/* open */
#include <fcntl.h>
int io_get_file_size(int fd, size_t *size) {
struct stat s;
int status = fstat(fd, &s);
if (status != 0) error(errno, "Unable to get file size");
*size = s.st_size;
return 0;
}
int io_open_file(char *filename, int *fd, void **map) {
/* open file */
*fd = open(filename, O_RDWR);
if (*fd == -1) error(errno, "Unable to open file %s", filename);
/* get file size */
size_t size;
int status = io_get_file_size(*fd, &size);
if (status != 0) return status;
/* map memory */
*map = mmap(NULL, size, PROT_WRITE | PROT_READ, MAP_SHARED, *fd, 0);
if (map == MAP_FAILED)
error(errno, "Failed to allocate map of size %zi bytes.", size);
return 0;
}
int io_close_file(int *fd, void **map) {
/* get file size */
size_t size;
int status = io_get_file_size(*fd, &size);
if (status != 0) return status;
/* unmap */
if (munmap(*map, size) != 0) error(errno, "Unable to unmap the file");
close(*fd);
return 0;
}
int io_read_mask(const struct header *h, void *map, size_t *offset,
size_t *mask, size_t *diff_offset) {
/* read mask */
if (mask) {
*mask = 0;
memcpy(mask, map + *offset, LOGGER_MASK_SIZE);
}
*offset += LOGGER_MASK_SIZE;
/* read offset */
if (diff_offset) {
*diff_offset = 0;
memcpy(diff_offset, map + *offset, LOGGER_OFFSET_SIZE);
}
*offset += LOGGER_OFFSET_SIZE;
return 0;
}
int io_read_data(void *map, const size_t size, void *p, size_t *offset) {
memcpy(p, map + *offset, size);
*offset += size;
return 0;
};
int io_write_data(void *map, const size_t size, const void *p, size_t *offset) {
memcpy(map + *offset, p, size);
*offset += size;
return 0;
};
#ifndef __SWIFT_LOGGER_IO_H__
#define __SWIFT_LOGGER_IO_H__
#include "header.h"
#include "tools.h"
#include <stdio.h>
#include <stdlib.h>
/**
* @brief get the size of a file
*
* @param fd file id
* @param size out: file size
*
* @return error code
*/
int io_get_file_size(int fd, size_t *size);
/**
* @brief Open a file and map it
*
* @param filename file to read
* @param fd out: file id
* @param map out: file mapping
*
* @return error code
*/
int io_open_file(char *filename, int *fd, void **map);
/**
* @brief Close a file and unmap it
*
* @param fd file id
* @param map file mapping
*
* @return error code
*/
int io_close_file(int *fd, void **map);
/**
* @brief read a single value in a file
*
* @param map file mapping
* @param size size of the chunk to read
* @param p pointer where to store the data
* @param offset In: position to read, Out: shifted by size
*
* @return error code
*/
int io_read_data(void *map, const size_t size, void *p, size_t *offset);
/**
* @brief write a single value in a file
*
* @param map file mapping
* @param size size of the chunk to write
* @param p pointer to the data
* @param offset In: position to write, Out: shifted by size
*
* @return error code
*/
int io_write_data(void *map, const size_t size, const void *p, size_t *offset);
/**
* @brief read a maks with its offset
*
* @param h #header file structure
* @param map file mapping
* @param offset In: position in the file, Out: shifted by the mask + offset
* size
* @param mask mask read
* @param diff_offset offset difference to previous/next corresponding chunk
*
* @return error code
*/
int io_read_mask(const struct header *h, void *map, size_t *offset,
size_t *mask, size_t *diff_offset);
#endif // __SWIFT_LOGGER_IO_H__
#include "header.h"
#include "io.h"
#include "particle.h"
#include "timeline.h"
#include <Python.h>
#include <errno.h>
#include <numpy/arrayobject.h>
#include <stdio.h>
#include <stdlib.h>
/**
* @brief Reverse offset in dump file
*
* @param filename string filename of the dump file
*/
int reverseOffset(char *filename) {
struct header h;
/* open file */
int fd;
void *map;
if (io_open_file(filename, &fd, &map) != 0) return 1;
/* read header */
if (header_read(&h, map) != 0) return 1;
header_print(&h);
/* check offset direction */
if (h.forward_offset) {
error_no_return(EIO, "Offset are already reversed");
return 1;
}
/* compute file size */
size_t sz;
int status = io_get_file_size(fd, &sz);
if (status != 0) return 1;
size_t offset;
int error_code;
#ifdef SWIFT_DEBUG_CHECKS
/* check offset */
printf("Check offsets...\n");
offset = h.offset_first;
while (offset < sz) {
error_code = tools_check_offset(&h, map, &offset);
if (error_code != 0) return 1;
}
printf("Check done\n");
#endif
/* reverse header offset */
header_change_offset_direction(&h, map);
offset = h.offset_first;
/* reverse chunks */
printf("Reversing offsets...\n");
while (offset < sz) {
error_code = tools_reverse_offset(&h, map, &offset);
if (error_code != 0) return 1;
}
printf("Reversing done\n");
#ifdef SWIFT_DEBUG_CHECKS
/* check offset */
printf("Check offsets...\n");
offset = h.offset_first;
while (offset < sz) {
error_code = tools_check_offset(&h, map, &offset);
if (error_code != 0) return 1;
}
printf("Check done\n");
#endif
/* free internal variables */
header_free(&h);
if (io_close_file(&fd, &map) != 0) return 1;
return 0;
}
/**
* @brief load data from the offset without any interpolation
*
* @param offset PyArrayObject list of offset for each particle
* @param filename string filename of the dump file
* @return dictionnary containing the data read
*/
static PyObject *loadFromIndex(__attribute__((unused)) PyObject *self,
PyObject *args) {
struct header h;
/* input */
PyArrayObject *offset = NULL;
char *filename = NULL;
/* output */
PyArrayObject *pos = NULL;
PyArrayObject *vel = NULL;
PyArrayObject *acc = NULL;
PyArrayObject *entropy = NULL;
PyArrayObject *h_sph = NULL;
PyArrayObject *rho = NULL;
PyArrayObject *mass = NULL;
PyArrayObject *id = NULL;
size_t time_offset;
/* parse arguments */
if (!PyArg_ParseTuple(args, "OsL", &offset, &filename, &time_offset))
return NULL;
if (!PyArray_Check(offset)) {
error_no_return(ENOTSUP, "Offset is not a numpy array");
return NULL;
}
if (PyArray_NDIM(offset) != 1) {
error_no_return(ENOTSUP, "Offset is not a 1 dimensional array");
return NULL;
}
if (PyArray_TYPE(offset) != NPY_UINT64) {
error_no_return(ENOTSUP, "Offset does not contain unsigned int");
return NULL;
}
/* open file */
int fd;
void *map;
if (io_open_file(filename, &fd, &map) != 0) return NULL;
/* read header */
if (header_read(&h, map) != 0) return NULL;
/* reverse offset if needed */
if (!h.forward_offset) {
if (io_close_file(&fd, &map) != 0) return NULL;
if (reverseOffset(filename) != 0) return NULL;
if (io_open_file(filename, &fd, &map) != 0) return NULL;
/* Reset header */
header_free(&h);
if (header_read(&h, map) != 0) return NULL;
}
/* read timestamps */
struct time_array times;
if (time_array_init(&times, &h, map, fd) != 0) return NULL;
time_array_print(&times);
/* get required time */
double time = time_array_get_time(&times, time_offset);
/* init array */
npy_intp dim[2];
dim[0] = PyArray_DIMS(offset)[0];
dim[1] = DIM;
/* init output */
if (header_is_present(&h, "position")) {
pos = (PyArrayObject *)PyArray_SimpleNew(2, dim, NPY_DOUBLE);
}
if (header_is_present(&h, "velocity")) {
vel = (PyArrayObject *)PyArray_SimpleNew(2, dim, NPY_FLOAT);
}
if (header_is_present(&h, "acceleration")) {
acc = (PyArrayObject *)PyArray_SimpleNew(2, dim, NPY_FLOAT);
}
if (header_is_present(&h, "entropy")) {
entropy =
(PyArrayObject *)PyArray_SimpleNew(1, PyArray_DIMS(offset), NPY_FLOAT);
}
if (header_is_present(&h, "cutoff radius")) {
h_sph =
(PyArrayObject *)PyArray_SimpleNew(1, PyArray_DIMS(offset), NPY_FLOAT);
}
if (header_is_present(&h, "density")) {
rho =
(PyArrayObject *)PyArray_SimpleNew(1, PyArray_DIMS(offset), NPY_FLOAT);
}
if (header_is_present(&h, "consts")) {
mass =
(PyArrayObject *)PyArray_SimpleNew(1, PyArray_DIMS(offset), NPY_FLOAT);
id = (PyArrayObject *)PyArray_SimpleNew(1, PyArray_DIMS(offset), NPY_ULONG);
}
int error_code = 0;
/* check data type in particles */
if (!particle_check_data_type(&h)) {
error_no_return(ENOTSUP, "Particle data type are not compatible");
return NULL;
}
/* loop over all particles */
for (npy_intp i = 0; i < PyArray_DIMS(offset)[0]; i++) {
struct particle part;
size_t *offset_particle = (size_t *)PyArray_GETPTR1(offset, i);
error_code = particle_read(&part, &h, map, offset_particle, time,
reader_lin, &times);
if (error_code != 0) return NULL;
double *dtmp;
float *ftmp;
size_t *stmp;
/* copy data */
for (size_t k = 0; k < DIM; k++) {
if (pos) {
dtmp = PyArray_GETPTR2(pos, i, k);
*dtmp = part.pos[k];
}
if (vel) {
ftmp = PyArray_GETPTR2(vel, i, k);
*ftmp = part.vel[k];
}
if (acc) {
ftmp = PyArray_GETPTR2(acc, i, k);
*ftmp = part.acc[k];
}
}
if (entropy) {
ftmp = PyArray_GETPTR1(entropy, i);
*ftmp = part.entropy;
}
if (rho) {
ftmp = PyArray_GETPTR1(rho, i);
*ftmp = part.density;
}
if (h_sph) {
ftmp = PyArray_GETPTR1(h_sph, i);
*ftmp = part.h;
}
if (mass) {
ftmp = PyArray_GETPTR1(mass, i);
*ftmp = part.mass;
}
if (id) {
stmp = PyArray_GETPTR1(id, i);
*stmp = part.id;
}
}
header_free(&h);
/* construct return */
PyObject *dict = PyDict_New();
PyObject *key = PyUnicode_FromString("position");
PyDict_SetItem(dict, key, PyArray_Return(pos));
if (vel) {
key = PyUnicode_FromString("velocity");
PyDict_SetItem(dict, key, PyArray_Return(vel));
}
if (acc) {
key = PyUnicode_FromString("acceleration");
PyDict_SetItem(dict, key, PyArray_Return(acc));
}
if (entropy) {
key = PyUnicode_FromString("entropy");
PyDict_SetItem(dict, key, PyArray_Return(entropy));
}
if (rho) {
key = PyUnicode_FromString("rho");
PyDict_SetItem(dict, key, PyArray_Return(rho));
}
if (h_sph) {
key = PyUnicode_FromString("h_sph");
PyDict_SetItem(dict, key, PyArray_Return(h_sph));
}
if (mass) {
key = PyUnicode_FromString("mass");
PyDict_SetItem(dict, key, PyArray_Return(mass));
}
if (id) {
key = PyUnicode_FromString("id");
PyDict_SetItem(dict, key, PyArray_Return(id));
}
io_close_file(&fd, &map);
return dict;
}
/**
* @brief Reverse offset in dump file
*
* @param filename string filename of the dump file
*/
static PyObject *pyReverseOffset(__attribute__((unused)) PyObject *self,
PyObject *args) {
/* input */
char *filename = NULL;
/* parse arguments */
if (!PyArg_ParseTuple(args, "s", &filename)) return NULL;
if (reverseOffset(filename) != 0) return NULL;
return Py_BuildValue("");
}
/* definition of the method table */
static PyMethodDef logger_loaderMethods[] = {
{"loadFromIndex", loadFromIndex, METH_VARARGS,
"Load snapshot directly from the offset in an index file."},
{"reverseOffset", pyReverseOffset, METH_VARARGS,
"Reverse the offset (from pointing backward to forward)."},
{NULL, NULL, 0, NULL} /* Sentinel */
};
static struct PyModuleDef logger_loadermodule = {
PyModuleDef_HEAD_INIT,
"logger_loader",
"Module reading a SWIFTsim logger snapshot",
-1,
logger_loaderMethods,
NULL, /* m_slots */
NULL, /* m_traverse */
NULL, /* m_clear */
NULL /* m_free */
};
PyMODINIT_FUNC PyInit_logger_loader(void) {
PyObject *m;
m = PyModule_Create(&logger_loadermodule);
if (m == NULL) return NULL;
import_array();
return m;
}
#include "particle.h"
#include "header.h"
#include "io.h"
#include "timeline.h"
#include "tools.h"
#include <stdio.h>
#include <stdlib.h>
void particle_print(const struct particle *p) {
printf("ID: %lu\n", p->id);
printf("Mass: %g\n", p->mass);
printf("Time: %g\n", p->time);
printf("Cutoff Radius: %g\n", p->h);
printf("Position: (%g, %g, %g)\n", p->pos[0], p->pos[1], p->pos[2]);
printf("Velocity: (%g, %g, %g)\n", p->vel[0], p->vel[1], p->vel[2]);
printf("Acceleration: (%g, %g, %g)\n", p->acc[0], p->acc[1], p->acc[2]);
printf("Entropy: %g\n", p->entropy);
printf("Density: %g\n", p->density);
}
int particle_check_data_type(__attribute__((unused)) const struct header *h) {
printf("TODO check_data_type\n");
return 1;
}
void particle_init(struct particle *part) {
for (size_t k = 0; k < DIM; k++) {
part->pos[k] = 0;
part->vel[k] = 0;
part->acc[k] = 0;
}
part->entropy = -1;
part->density = -1;
part->h = -1;
part->mass = -1;
part->id = SIZE_MAX;
}
int particle_read_field(struct particle *part, void *map, size_t *offset,
const char *field, const size_t size) {
void *p = NULL;
if (strcmp("position", field) == 0) {
p = &part->pos;
} else if (strcmp("velocity", field) == 0) {
p = &part->vel;
} else if (strcmp("acceleration", field) == 0) {
p = &part->acc;
} else if (strcmp("entropy", field) == 0) {
p = &part->entropy;
} else if (strcmp("cutoff radius", field) == 0) {
p = &part->h;
} else if (strcmp("density", field) == 0) {
p = &part->density;
} else if (strcmp("consts", field) == 0) {
p = malloc(size);
} else {
error(ENOTSUP, "Type %s not defined", field);
}
int error_code = io_read_data(map, size, p, offset);
if (error_code != 0) return error_code;
if (strcmp("consts", field) == 0) {
part->mass = 0;
part->id = 0;
memcpy(&part->mass, p, sizeof(float));
p += sizeof(float);
memcpy(&part->id, p, sizeof(size_t));
p -= sizeof(float);
free(p);
}
return error_code;
}
int particle_read(struct particle *part, const struct header *h, void *map,
size_t *offset, const double time, const int reader,
struct time_array *times) {
int error_code = 0;
size_t mask = 0;
size_t h_offset = 0;
particle_init(part);
error_code = io_read_mask(h, map, offset, &mask, &h_offset);
if (error_code != 0) return error_code;
if (mask != 127) error(EIO, "Unexpected mask: %lu", mask);
for (size_t i = 0; i < h->nber_mask; i++) {
if (mask & h->masks[i]) {
error_code = particle_read_field(part, map, offset, h->masks_name[i],
h->masks_size[i]);
if (error_code != 0) return error_code;
}
}
if (times) /* move offset by 1 in order to be in the required chunk */
part->time = time_array_get_time(times, *offset - 1);
else
part->time = -1;
/* end of const case */
if (reader == reader_const) return 0;
/* read next particle */
struct particle part_next;
if (!h->forward_offset) error(ENOSYS, "TODO");
if (h_offset == 0) return 0;
/* get absolute offset of next particle */
h_offset += *offset - header_get_mask_size(h, mask) - LOGGER_MASK_SIZE - LOGGER_OFFSET_SIZE;
part_next.time = time_array_get_time(times, h_offset);
/* previous part exists */
error_code = particle_read(&part_next, h, map, &h_offset, part_next.time,
reader_const, times);
if (error_code != 0) return error_code;
error_code = particle_interpolate(part, &part_next, time);
if (error_code != 0) return error_code;
return 0;
}
int particle_interpolate(struct particle *part_curr,
const struct particle *part_next,
const double time) {
if (!part_curr) error(EFAULT, "part_curr is NULL");
if (!part_next) error(EFAULT, "part_next is NULL");
#ifdef SWIFT_DEBUG_CHECKS
if (part_next->time <= part_curr->time)
error(EIO, "Wrong particle order (next before current)");
if ((time < part_curr->time) || (part_next->time < time))
error(EIO,
"Interpolating, not extrapolating (particle time: %f, "
"interpolating time: %f, next particle time: %f)",
part_curr->time, time, part_next->time);
#endif
double scaling = part_next->time - part_curr->time;
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;
ftmp = (part_next->acc[i] - part_curr->acc[i]);
part_curr->acc[i] += ftmp * scaling;
}
/* interpolate scalars */
ftmp = (part_next->entropy - part_curr->entropy);
part_curr->entropy += ftmp * scaling;
/* set time */
part_curr->time = time;
return 0;
}
#ifndef __PARTICLE_H__
#define __PARTICLE_H__
#include "header.h"
#include "timeline.h"
#include "tools.h"
#include <stdio.h>
#include <stdlib.h>
#define DIM 3
struct particle {
/* position */
double pos[DIM];
/* velocity */
float vel[DIM];
/* acceleration */
float acc[DIM];
/* entropy */
float entropy;
/* cutoff radius */
float h;
/* density */
float density;
/* mass */
float mass;
/* id */
size_t id;
/* time */
double time;
};
enum reader_type {
reader_const,
reader_lin,
};
/**
* @brief print a particle
*
* @param p #particle particle to print
*/
void particle_print(const struct particle *p);
/**
* @brief read a particle in the dump file
*
* @param part #particle particle to update
* @param h #header structure of the file
* @param map file mapping
* @param offset offset of the chunk to read (update it to the end of the chunk)
* @param time time to interpolate (if #reader_type is an interpolating one)
* @param reader #reader_type
* @param times #time_array times in the dump
*
* @return error code
*/
int particle_read(struct particle *part, const struct header *h, void *map,
size_t *offset, const double time, const int reader,
struct time_array *times);
/**
* @brief Check if dump data type are compatible with the particle type
*
* @param h #header structure of the file
*
* @return error code
*/
int particle_check_data_type(const struct header *h);
/**
* @brief initialize a particle
*
* @param part #particle particle to initialize
*/
void particle_init(struct particle *part);
/**
* @brief read a single field for a particle
*
* @param part @particle particle to update
* @param map file mapping
* @param offset In: read position, Out: input shifted by the required amount of
* data
* @param field field to read
* @param size number of bits to read
*
* @return error code
*/
int particle_read_field(struct particle *part, void *map, size_t *offset,
const char *field, const size_t size);
/**
* @brief interpolate two particles at a given time
*
* @param part_curr #particle In: current particle (before time), Out:
* interpolated particle
* @param part_next #particle next particle (after time)
* @param time interpolation time
*
* @return error code
*/
int particle_interpolate(struct particle *part_curr,
const struct particle *part_next,
const double time);
#endif //__PARTICLE_H__
#include "timeline.h"
#include "io.h"
double time_convert_to_double(const integertime_t ti, const double timeBase) {
return ti * timeBase; // should add timebegin
}
integertime_t time_convert_to_integer(const double_t ti,
const double timeBase) {
return ti / timeBase; // should add timebegin
}
int time_read(integertime_t *timestamp, double *time, const struct header *h, void *map,
size_t *offset) {
int error_code = 0;
size_t mask = 0;
size_t ind = 0;
size_t prev_offset = 0;
*timestamp = 0;
*time = 0;
/* read chunck header */
error_code = io_read_mask(h, map, offset, &mask, &prev_offset);
if (error_code != 0) return error_code;
#ifdef SWIFT_DEBUG_CHECKS
/* check if timestamp is present */
if (!header_is_present_and_get_index(h, "timestamp", &ind))
error(EIO, "Header does not contain a timestamp");
/* check if timestamp */
if (h->masks[ind] != mask) error(EIO, "Not a timestamp");
#endif
/* read data */
// TODO
error_code = io_read_data(map, sizeof(unsigned long long int), timestamp, offset);
error_code = io_read_data(map, sizeof(double), time, offset);
return error_code;
}
int time_first_timestamp(const struct header *h, void *map, size_t *offset) {
*offset = h->offset_first;
int test;
size_t i;
if (!header_is_present_and_get_index(h, "timestamp", &i))
error(EIO, "Time stamp not present in header");
size_t tmp = *offset;
size_t mask = 0;
test = io_read_mask(h, map, offset, &mask, NULL);
if (test != 0) return test;
if (mask != h->masks[i]) error(EIO, "Dump should begin by timestep");
*offset = tmp;
return 0;
}
int time_array_init(struct time_array *t, const struct header *h, void *map,
int fd) {
t->next = NULL;
t->prev = NULL;
/* get first time stamp */
size_t offset = 0;
time_first_timestamp(h, map, &offset);
integertime_t timestamp = 0;
double time = 0;
/* get file size */
size_t file_size;
int test = io_get_file_size(fd, &file_size);
while (offset < file_size) {
/* read time */
t->offset = offset;
size_t tmp_offset = offset;
time_read(&timestamp, &time, h, map, &tmp_offset);
t->timestamp = timestamp;
t->time = time;
/* get next chunk */
test = tools_get_next_chunk(h, map, &offset, fd);
if (test == -1)
break;
else if (test != 0)
return test;
/* allocate next time_array */
struct time_array *tmp = malloc(sizeof(struct time_array));
tmp->prev = t;
tmp->next = NULL;
t->next = tmp;
t = tmp;
}
/* free unused time_array */
struct time_array *tmp = t->prev;
tmp->next = NULL;
free(t);
return 0;
}
integertime_t time_array_get_integertime(struct time_array *t, const size_t offset) {
const struct time_array *tmp = time_array_get_time_array(t, offset);
return tmp->timestamp;
}
double time_array_get_time(struct time_array *t, const size_t offset) {
const struct time_array *tmp = time_array_get_time_array(t, offset);
return tmp->time;
}
struct time_array *time_array_get_time_array(struct time_array *t,
const size_t offset) {
#ifdef SWIFT_DEBUG_CHECKS
if (!t) return 0;
#endif
while (t->next) {
if (t->offset > offset) break;
t = t->next;
}
if (t->prev == NULL) return t;
struct time_array *tmp = t->prev;
return tmp;
}
void time_array_free(struct time_array *t) {
struct time_array *tmp;
while (t->next) {
tmp = t->next;
free(t);
t = tmp;
}
}
void time_array_print(const struct time_array *t) {
size_t threshold = 4;
size_t n = time_array_count(t);
size_t i = 0;
size_t up_threshold = n - threshold;
printf("Times (size %lu): [%lli (%g)", n, t->timestamp, t->time);
while (t->next) {
i += 1;
t = t->next;
if (i < threshold || i > up_threshold) printf(", %lli (%g)", t->timestamp, t->time);
if (i == threshold) printf(", ...");
}
printf("]\n");
}
void time_array_print_offset(const struct time_array *t) {
size_t threshold = 4;
size_t n = time_array_count(t);
size_t i = 0;
size_t up_threshold = n - threshold;
printf("Times (size %lu): [%lu", n, t->offset);
while (t->next) {
i += 1;
t = t->next;
if (i < threshold || i > up_threshold) printf(", %lu", t->offset);
if (i == threshold) printf(", ...");
}
printf("]\n");
}
size_t time_array_count(const struct time_array *t) {
size_t count = 1;
while (t->next) {
t = t->next;
count += 1;
}
return count;
}
#ifndef __TIMELINE_H__
#define __TIMELINE_H__
#include "header.h"
#include "tools.h"
typedef char timebin_t;
typedef long long integertime_t;
struct time_array {
void *next;
void *prev;
integertime_t timestamp;
double time;
size_t offset;
};
/**
* @brief convert an integer time to a real time
*
* @param ti integer time to convert
* @param timeBase time base of the integer time
*
* @return converted time
*/
double time_convert_to_double(const integertime_t ti, const double timeBase);
/**
* @brief convert an integer time to a real time
*
* @param ti integer time to convert
* @param timeBase time base of the integer time
*
* @return converted time
*/
integertime_t time_convert_to_integer(const double ti, const double timeBase);
/**
* @brief read a time stamp
*
* @param timestamp timestamp read
* @param time time read
* @param h #header file structure
* @param map file mapping
* @param offset In: position in the file, Out: shifted by the timestamp
*/
int time_read(integertime_t *timestamp, double *time, const struct header *h, void *map,
size_t *offset);
/**
* @brief Initialize a time array
*
* @param t #time_array to initialize
* @param h #header file structure
* @param map file mapping
* @param fd file id
*/
int time_array_init(struct time_array *t, const struct header *h, void *map,
int fd);
/**
* @brief access the time of a given chunk (by its offset)
*
* @param t #time_array to access
* @param offset offset of the chunk
*
* @return integer time of the chunk
*/
integertime_t time_array_get_integertime(struct time_array *t, const size_t offset);
/**
* @brief access the time of a given chunk (by its offset)
*
* @param t #time_array to access
* @param offset offset of the chunk
*
* @return time of the chunk
*/
double time_array_get_time(struct time_array *t, const size_t offset);
/**
* @brief access the #time_array of a given chunk (by its offset)
*
* @param t #time_array to access
* @param offset offset of the chunk
*
* @return pointer to the requested #time_array
*/
struct time_array *time_array_get_time_array(struct time_array *t,
const size_t offset);
/**
* @brief free memory of a #time_array
*
* @param t #time_array to free
*/
void time_array_free(struct time_array *t);
/**
* @brief print a #time_array
*
* @param t #time_array to print
*/
void time_array_print(const struct time_array *t);
/**
* @brief print a #time_array (offset)
*
* @param t #time_array to print
*/
void time_array_print_offset(const struct time_array *t);
/**
* @brief count number of element in #time_array
*
* @param t #time_array to count
*
* @return number of element
*/
size_t time_array_count(const struct time_array *t);
/**
* @brief get offset of first timestamp
*
* @param h file #header
* @param map file mapping
* @param offset out: offset of first timestamp
*
* @return error code
*/
int time_first_timestamp(const struct header *h, void *map, size_t *offset);
#endif // __TIMELINE_H__
#include "tools.h"
#include "header.h"
#include "io.h"
#include "particle.h"
#include <stdio.h>
int tools_get_next_chunk(const struct header *h, void *map, size_t *offset,
int fd) {
if (h->forward_offset)
return _tools_get_next_chunk_forward(h, map, offset);
else
return _tools_get_next_chunk_backward(h, map, offset, fd);
}
int _tools_get_next_chunk_forward(const struct header *h, void *map,
size_t *offset) {
size_t diff_offset = 0;
/* read offset */
int test = io_read_mask(h, map, offset, NULL, &diff_offset);
if (test != 0) return test;
if (diff_offset == 0) return -1;
/* set absolute offset */
*offset += diff_offset - LOGGER_MASK_SIZE - LOGGER_OFFSET_SIZE;
return 0;
}
int _tools_get_next_chunk_backward(const struct header *h, void *map,
size_t *offset, int fd) {
#ifndef SWIFT_DEBUG_CHECKS
error(ENOTSUP, "Should not be used, method too slow");
#endif
size_t current_offset = *offset;
size_t chunk_header = LOGGER_MASK_SIZE + LOGGER_OFFSET_SIZE;
size_t file_size = 0;
int test = io_get_file_size(fd, &file_size);
if (test != 0) return test;
while (current_offset < file_size) {
size_t mask = 0;
size_t prev_offset;
int test = io_read_mask(h, map, &current_offset, &mask, &prev_offset);
if (test != 0) return test;
prev_offset = current_offset - prev_offset - chunk_header;
if (*offset == prev_offset) {
*offset = current_offset - chunk_header;
return 0;
}
current_offset += header_get_mask_size(h, mask);
}
return -1;
}
int tools_reverse_offset(const struct header *h, void *map, size_t *offset) {
size_t mask = 0;
size_t prev_offset = 0;
const size_t cur_offset = *offset;
/* read mask + offset */
int error_code = io_read_mask(h, map, offset, &mask, &prev_offset);
if (error_code != 0) return error_code;
/* write offset of zero (in case it is the last chunk) */
const size_t zero = 0;
*offset -= LOGGER_OFFSET_SIZE;
error_code = io_write_data(map, LOGGER_OFFSET_SIZE, &zero, offset);
if (error_code != 0) error(errno, "Unable to write data while reversing");
/* set offset after current chunk */
*offset += header_get_mask_size(h, mask);
/* first chunks do not have a previous partner */
if (prev_offset == cur_offset) return 0;
if (prev_offset > cur_offset)
error(EIO, "Unexpected offset, header %lu, current %lu", prev_offset,
cur_offset);
/* modify previous offset */
size_t abs_prev_offset = cur_offset - prev_offset + LOGGER_MASK_SIZE;
error_code = io_write_data(map, LOGGER_OFFSET_SIZE, &prev_offset, &abs_prev_offset);
if (error_code != 0) error(errno, "Unable to write offset");
#ifdef SWIFT_DEBUG_CHECKS
size_t prev_mask = 0;
abs_prev_offset -= LOGGER_MASK_SIZE + LOGGER_OFFSET_SIZE;
error_code = io_read_mask(h, map, &abs_prev_offset, &prev_mask, NULL);
if (prev_mask != mask)
error(EIO, "Unexpected mask: %lu (at %lu), got %lu (at %lu)", mask, *offset,
prev_mask, prev_offset);
#endif // SWIFT_DEBUG_CHECKS
return 0;
}
int tools_check_offset(const struct header *h, void *map, size_t *offset) {
#ifndef SWIFT_DEBUG_CHECKS
error(ENOTSUP, "Should not check in non debug mode");
#endif
size_t tmp = *offset;
size_t mask;
size_t pointed_offset;
/* read mask + offset */
int error_code = io_read_mask(h, map, offset, &mask, &pointed_offset);
if (error_code != 0) return error_code;
/* get absolute offset */
if (h->forward_offset)
pointed_offset += tmp;
else {
if (tmp < pointed_offset)
error(EIO, "Offset too large (%lu) at %lu with mask %lu", pointed_offset,
tmp, mask);
pointed_offset = tmp - pointed_offset;
}
/* set offset after current chunk */
*offset += header_get_mask_size(h, mask);
if (pointed_offset == tmp || pointed_offset == 0) return 0;
/* read mask of the pointed chunk */
size_t pointed_mask = 0;
error_code = io_read_mask(h, map, &pointed_offset, &pointed_mask, NULL);
/* check masks */
if (pointed_mask != mask)
error(EIO, "Error in the offset (mask %lu != %lu) at %lu and %lu", mask,
pointed_mask,
*offset - header_get_mask_size(h, mask) - LOGGER_MASK_SIZE - LOGGER_OFFSET_SIZE,
pointed_offset - LOGGER_MASK_SIZE - LOGGER_OFFSET_SIZE);
if (pointed_mask == 128) return 0;
struct particle part;
error_code = particle_read(&part, h, map, &tmp, 0, reader_const, NULL);
size_t id = part.id;
tmp = pointed_offset - LOGGER_MASK_SIZE - LOGGER_OFFSET_SIZE;
error_code = particle_read(&part, h, map, &tmp, 0, reader_const, NULL);
if (id != part.id)
error(EIO, "Offset wrong, id incorrect (%lu != %lu) at %lu", id, part.id,
tmp);
return 0;
}
#ifndef __TOOLS_H__
#define __TOOLS_H__
#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
#include <Python.h>
#include <stdio.h>
#include <stdlib.h>
#define STRING_SIZE 200
struct header;
/* Set the error message for python and return the error code
* WARNING for python, you need to return NULL and not the error code
*/
#define error(err, s, ...) \
({ \
error_no_return(err, s, ##__VA_ARGS__); \
return err; \
})
/* Same as error but does not return the error code.
* Should be used for python functions */
#define error_no_return(err, s, ...) \
({ \
char error_msg[STRING_SIZE]; \
sprintf(error_msg, "%s:%s():%i: " s ": %s\n", __FILE__, __FUNCTION__, \
__LINE__, ##__VA_ARGS__, strerror(err)); \
PyErr_SetString(PyExc_RuntimeError, error_msg); \
})
#define message(s, ...) \
({ \
printf("%s:%s():%i: " s "\n", __FILE__, __FUNCTION__, __LINE__, \
##__VA_ARGS__); \
})
/**
* @brief get the offset of the next corresponding chunk
*
* @param h #header structure of the file
* @param map file mapping
* @param offset In: initial offset, Out: offset of the next chunk
* @param fd file id
*
* @return error code, -1 if no next chunk
*/
int tools_get_next_chunk(const struct header *h, void *map, size_t *offset,
int fd);
/**
* @brief internal function of #tools_get_next_chunk. Should not be used (very
* slow)
*
* @param h #header structure of the file
* @param map file mapping
* @param offset In: initial offset, Out: offset of the next chunk
* @param fd file id
*
* @return error code, -1 if no next chunk
*/
int _tools_get_next_chunk_backward(const struct header *h, void *map,
size_t *offset, int fd);
/**
* @brief internal function of #tools_get_next_chunk. Should not be used outside
*
* @param h #header structure of the file
* @param map file mapping
* @param offset In: initial offset, Out: offset of the next chunk
*
* @return error code, -1 if no next chunk
*/
int _tools_get_next_chunk_forward(const struct header *h, void *map,
size_t *offset);
/**
* @brief switch side offset
*
* From current chunk, switch side of the offset of the previous one.
* @param h #header structure of the file
* @param map file mapping
* @param offset In: initial offset, Out: offset of next chunk
*
* @return error code
*/
int tools_reverse_offset(const struct header *h, void *map, size_t *offset);
/**
* @brief debugging function checking the offset of a chunk
*
* Compare the mask with the one pointed by the header.
* if the chunk is a particle, check the id too.
* @param h #header structure of the file
* @param map file mapping
* @param offset In: chunk offset, Out: offset after the chunk
*
* @return error code
*/
int tools_check_offset(const struct header *h, void *map, size_t *offset);
#endif //__TOOLS_H__
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment