From f2cb03e33c01a11d43dbaf3f7b16d61d84ea9e52 Mon Sep 17 00:00:00 2001 From: Loic Hausammann <loic.hausammann@protonmail.ch> Date: Fri, 11 Jun 2021 08:17:33 +0000 Subject: [PATCH] Make the CSDS indepedent from SWIFT's module --- .clang-format | 19 + doc/source/NewScheme/index.rst | 28 +- examples/create_snapshot.py | 4 +- examples/reader_example.py | 2 +- format.sh | 80 ++++ src/Makefile.am | 54 +-- src/chemistry/GEAR/csds_chemistry.c | 39 -- src/chemistry/GEAR/csds_chemistry.h | 182 -------- src/chemistry/none/csds_chemistry.c | 35 -- src/chemistry/none/csds_chemistry.h | 124 ----- src/csds_chemistry.h | 40 -- src/csds_cosmology.c | 66 ++- src/csds_cosmology.h | 15 +- src/csds_fields.c | 39 ++ src/csds_fields.h | 427 ++++++++++++++++++ src/csds_gravity.h | 38 -- src/csds_header.c | 274 ++++++----- src/csds_header.h | 62 +-- src/csds_hydro.h | 61 --- src/csds_index.c | 2 +- src/csds_index.h | 2 +- src/csds_interpolation.h | 81 ++-- src/csds_loader_io.c | 2 +- src/csds_loader_io.h | 2 +- src/csds_logfile.c | 2 +- src/csds_logfile.h | 2 +- src/csds_parameters.c | 2 +- src/csds_parameters.h | 4 +- src/csds_particle.c | 241 ++++++---- src/csds_particle.h | 28 +- src/csds_python_tools.h | 2 +- src/csds_python_wrapper.c | 204 +++------ src/csds_reader.c | 283 +++++------- src/csds_reader.h | 18 +- src/csds_reader_generate_index.c | 131 ++---- src/csds_star_formation.h | 38 -- src/csds_stars.h | 39 -- src/csds_time.c | 23 +- src/csds_time.h | 2 +- src/csds_tools.c | 140 +----- src/csds_tools.h | 63 +-- src/gravity/MultiSoftening/csds_gravity.c | 34 -- src/gravity/MultiSoftening/csds_gravity.h | 141 ------ src/hydro/Gadget2/csds_hydro.c | 37 -- src/hydro/Gadget2/csds_hydro.h | 147 ------ src/hydro/SPHENIX/csds_hydro.c | 39 -- src/hydro/SPHENIX/csds_hydro.h | 209 --------- src/quick_sort.c | 2 +- src/quick_sort.h | 2 +- src/star_formation/GEAR/csds_star_formation.c | 32 -- src/star_formation/GEAR/csds_star_formation.h | 119 ----- src/star_formation/none/csds_star_formation.c | 28 -- src/star_formation/none/csds_star_formation.h | 77 ---- src/stars/Basic/csds_stars.c | 35 -- src/stars/Basic/csds_stars.h | 139 ------ src/stars/GEAR/csds_stars.c | 36 -- src/stars/GEAR/csds_stars.h | 142 ------ tests/generate_log.h | 4 +- tests/testLogfileHeader.c | 43 +- tests/testLogfileReader.c | 54 ++- tests/testQuickSort.c | 2 + tests/testVirtualReality.c | 24 +- 62 files changed, 1304 insertions(+), 2942 deletions(-) create mode 100644 .clang-format create mode 100755 format.sh delete mode 100644 src/chemistry/GEAR/csds_chemistry.c delete mode 100644 src/chemistry/GEAR/csds_chemistry.h delete mode 100644 src/chemistry/none/csds_chemistry.c delete mode 100644 src/chemistry/none/csds_chemistry.h delete mode 100644 src/csds_chemistry.h create mode 100644 src/csds_fields.c create mode 100644 src/csds_fields.h delete mode 100644 src/csds_gravity.h delete mode 100644 src/csds_hydro.h delete mode 100644 src/csds_star_formation.h delete mode 100644 src/csds_stars.h delete mode 100644 src/gravity/MultiSoftening/csds_gravity.c delete mode 100644 src/gravity/MultiSoftening/csds_gravity.h delete mode 100644 src/hydro/Gadget2/csds_hydro.c delete mode 100644 src/hydro/Gadget2/csds_hydro.h delete mode 100644 src/hydro/SPHENIX/csds_hydro.c delete mode 100644 src/hydro/SPHENIX/csds_hydro.h delete mode 100644 src/star_formation/GEAR/csds_star_formation.c delete mode 100644 src/star_formation/GEAR/csds_star_formation.h delete mode 100644 src/star_formation/none/csds_star_formation.c delete mode 100644 src/star_formation/none/csds_star_formation.h delete mode 100644 src/stars/Basic/csds_stars.c delete mode 100644 src/stars/Basic/csds_stars.h delete mode 100644 src/stars/GEAR/csds_stars.c delete mode 100644 src/stars/GEAR/csds_stars.h diff --git a/.clang-format b/.clang-format new file mode 100644 index 0000000..0ec6295 --- /dev/null +++ b/.clang-format @@ -0,0 +1,19 @@ +--- +Language: Cpp +BasedOnStyle: Google +KeepEmptyLinesAtTheStartOfBlocks: true +PenaltyBreakAssignment: 2 +IncludeCategories: + - Regex: '^"(llvm|llvm-c|clang|clang-c)/' + Priority: 3 + SortPriority: 3 + - Regex: '^(<|"(gtest|gmock|isl|json)/)' + Priority: 4 + - Regex: '<[[:alnum:].]+>' + Priority: 5 + - Regex: '.*' + Priority: 2 + SortPriority: 0 + - Regex: config.h + Priority: 1 +... diff --git a/doc/source/NewScheme/index.rst b/doc/source/NewScheme/index.rst index 87326ed..a7b8656 100644 --- a/doc/source/NewScheme/index.rst +++ b/doc/source/NewScheme/index.rst @@ -8,25 +8,17 @@ Implementing new Scheme To add the CSDS to a new scheme, a few files need to be created. Let's use the Gadget2 module for the illustration. -First, the files :file:`src/hydro/Gadget2/hydro_csds.c` and its header :file:`hydro_csds.h` describe +First, the file :file:`src/hydro/Gadget2/hydro_csds.h` describes how to write the particles inside the file. -They define the variable :command:`hydro_csds_field_names` that provides the name of the fields. -The names will be used within the python API. -The enum :command:`hydro_csds_fields` allows to loop over all the fields when writing / reading them. -Three functions are required within the header. -One initializes the mask, one computes the mask for a given step and the last one writes the particles -into the logfile. +The most important function is :command:`csds_hydro_define_fields` that defines all the fields. +The names should be the same than within the reader and correspond to the python API. +A conversion function is given for the acceleration :command:`csds_hydro_convert_acc`. For the compilation, the files need to be included in :file:`src/hydro_csds.h` and :file:`src/Makefile.am`. -Now in the reader, more or less the same is needed. -The files are :file:`csds/src/hydro/Gadget2/csds_hydro.c` and its header :file:`csds_hydro.h`. -In the files, we define two arrays. -The first one, :command:`hydro_csds_local_to_global`, is set by the reader and -allows a quick access between the module's fields and the global ones. -The second one, :command:`hydro_csds_field_size` gives the size of each field. -Three functions are implemented and consists in linking the derivatives -together (e.g. linking the velocity to the position), interpolating a field and -defining the data type for the numpy arrays. -For the compilation, it is exactly the same than for the writer and the files are -:file:`csds/src/csds_hydro.h` and :file:`csds/src/Makefile.am`. +The reader is independent from the modules in order to simplify its usage. +Therefore the new fields need to be added within :file:`csds/src/csds_fields.h`, +:file:`csds/src/csds_fields.c` and :file:`csds/src/csds_particles.c`. +The two first files define the properties of the fields (name, size, derivatives +and the python interface). +In the particle file, the interpolation is defined for each field. diff --git a/examples/create_snapshot.py b/examples/create_snapshot.py index d4876da..29d5b04 100644 --- a/examples/create_snapshot.py +++ b/examples/create_snapshot.py @@ -73,8 +73,8 @@ def write_particle_type(snap, part_type, args): if fields_name is None: fields_name = reader.get_list_fields(part_type) - # Abort if only requesting a type not implemented - if "SpecialFlags" in fields_name: + # Skip non implemented types + if len(fields_name) == 0: print("Part type %i not implemented, skipping it" % part_type) continue diff --git a/examples/reader_example.py b/examples/reader_example.py index 1a9d0a7..ef4a6b8 100644 --- a/examples/reader_example.py +++ b/examples/reader_example.py @@ -62,7 +62,7 @@ for f in args.files: filename = f[:-5] else: raise Exception("It seems that you are not providing a logfile (.dump)") - with csds.Reader(filename, verbose=0) as reader: + with csds.Reader(filename, verbose=0, number_threads=1, number_index=10) as reader: # Ensure that the fields are present fields = ["Coordinates", "InternalEnergies", "ParticleIDs"] diff --git a/format.sh b/format.sh new file mode 100755 index 0000000..b279ad1 --- /dev/null +++ b/format.sh @@ -0,0 +1,80 @@ +#!/bin/bash + +# Clang format command, can be overridden using CLANG_FORMAT_CMD. +# We currrently use version 10.0 so any overrides should provide that. +clang=${CLANG_FORMAT_CMD:="clang-format-10"} + +# Formatting command +cmd="$clang -style=file $(git ls-files | grep '\.[ch]$')" + +# Test if `clang-format-10` works +command -v $clang > /dev/null +if [[ $? -ne 0 ]] +then + echo "ERROR: cannot find $clang" + exit 1 +fi + +# Print the help +function show_help { + echo -e "This script formats CSDS according to Google style" + echo -e " -h, --help \t Show this help" + echo -e " -t, --test \t Test if CSDS is well formatted" +} + +# Parse arguments (based on https://stackoverflow.com/questions/192249) +TEST=0 +while [[ $# -gt 0 ]] +do + key="$1" + + case $key in + # print the help and exit + -h|--help) + show_help + exit + ;; + # check if the code is well formatted + -t|--test) + TEST=1 + shift + ;; + # unknown option + *) + echo "Argument '$1' not implemented" + show_help + exit + ;; + esac +done + +# Run the required commands +if [[ $TEST -eq 1 ]] +then + # Note trapping the exit status from both commands in the pipe. Also note + # do not use -q in grep as that closes the pipe on first match and we get + # a SIGPIPE error. + echo "Testing if CSDS is correctly formatted" + $cmd -output-replacements-xml | grep "<replacement " > /dev/null + status=("${PIPESTATUS[@]}") + + # Trap if first command failed. Note 141 is SIGPIPE, that happens when no + # output + if [[ ${status[0]} -ne 0 ]] + then + echo "ERROR: $clang command failed" + exit 1 + fi + + # Check formatting + if [[ ${status[1]} -eq 0 ]] + then + echo "ERROR: needs formatting" + exit 1 + else + echo "...is correctly formatted" + fi +else + echo "Formatting CSDS" + $cmd -i +fi diff --git a/src/Makefile.am b/src/Makefile.am index 07cac9c..e5bdf19 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -32,60 +32,19 @@ GIT_CMD = @GIT_CMD@ EXTRA_LIBS = $(PROFILER_LIBS) $(TCMALLOC_LIBS) $(JEMALLOC_LIBS) $(TBBMALLOC_LIBS) $(HDF5_LIBS) $(FFTW_LIBS) $(GRACKLE_LIBS) \ $(VELOCIRAPTOR_LIBS) $(GSL_LIBS) -L../../src/.libs -lswiftsim -# MPI libraries. -# MPI_LIBS = $(MPI_THREAD_LIBS) -# MPI_FLAGS = -DWITH_MPI - # Build the libcsds library lib_LTLIBRARIES = libcsds.la -# Build a MPI-enabled version too? -# if HAVEMPI -# lib_LTLIBRARIES += libcsds_mpi.la -# endif - GRAVITY_SRC = gravity/MultiSoftening/csds_gravity.c -# Include the stars modules -if HAVE_STARS_BASIC -STARS_SRC = stars/Basic/csds_stars.c -endif -if HAVE_STARS_GEAR -STARS_SRC = stars/GEAR/csds_stars.c -endif - -# Include the hydro modules -if HAVE_SPHENIX -HYDRO_SRC = hydro/SPHENIX/csds_hydro.c -endif -if HAVE_GADGET2 -HYDRO_SRC = hydro/Gadget2/csds_hydro.c -endif - -# Include the chemistry modules -if HAVE_CHEMISTRY_NONE -CHEMISTRY_SRC = chemistry/none/csds_chemistry.c -endif -if HAVE_CHEMISTRY_GEAR -CHEMISTRY_SRC = chemistry/GEAR/csds_chemistry.c -endif - -# Include the star formation modules -if HAVE_STAR_FORMATION_DEFAULT -STAR_FORMATION_SRC = star_formation/none/csds_star_formation.c -endif -if HAVE_STAR_FORMATION_GEAR -STAR_FORMATION_SRC = star_formation/GEAR/csds_star_formation.c -endif # List required headers include_HEADERS = csds_header.h csds_loader_io.h csds_particle.h csds_time.h csds_tools.h include_HEADERS += csds_reader.h csds_logfile.h csds_index.h quick_sort.h csds_python_tools.h -include_HEADERS += csds_interpolation.h csds_parameters.h csds_cosmology.h +include_HEADERS += csds_interpolation.h csds_parameters.h csds_cosmology.h csds_fields.h # Common source files AM_SOURCES = csds_header.c csds_loader_io.c csds_particle.c csds_time.c csds_tools.c csds_reader.c AM_SOURCES += csds_logfile.c csds_index.c quick_sort.c csds_parameters.c csds_reader_generate_index.c -AM_SOURCES += csds_cosmology.c -AM_SOURCES += $(GRAVITY_SRC) $(STARS_SRC) $(HYDRO_SRC) $(CHEMISTRY_SRC) $(STAR_FORMATION_SRC) +AM_SOURCES += csds_cosmology.c csds_fields.c if HAVEPYTHON AM_SOURCES += csds_python_wrapper.c @@ -98,15 +57,6 @@ if HAVEPYTHON PYTHON_EXTRA_COMPILER_FLAG = -Wno-missing-field-initializers -Wno-cast-function-type -Wno-unused-function -Wno-unused-function endif -# Include files for distribution, not installation. -nobase_noinst_HEADERS = csds_hydro.h hydro/Gadget2/csds_hydro.h hydro/SPHENIX/csds_hydro.h -nobase_noinst_HEADERS += csds_stars.h stars/Basic/csds_stars.h stars/GEAR/csds_stars.h -nobase_noinst_HEADERS += csds_gravity.h gravity/MultiSoftening/csds_gravity.h -nobase_noinst_HEADERS += csds_chemistry.h chemistry/none/csds_chemistry.h chemistry/GEAR/csds_chemistry.h -nobase_noinst_HEADERS += csds_star_formation.h star_formation/none/csds_star_formation.h -nobase_noinst_HEADERS += star_formation/GEAR/csds_star_formation.h - - # Sources and flags for regular library libcsds_la_SOURCES = $(AM_SOURCES) libcsds_la_CFLAGS = $(AM_CFLAGS) $(PYTHON_EXTRA_COMPILER_FLAG) diff --git a/src/chemistry/GEAR/csds_chemistry.c b/src/chemistry/GEAR/csds_chemistry.c deleted file mode 100644 index a06ba2e..0000000 --- a/src/chemistry/GEAR/csds_chemistry.c +++ /dev/null @@ -1,39 +0,0 @@ -/******************************************************************************* - * 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 this object's header */ -#include "csds_chemistry.h" - -/* Local headers */ -#include "csds_tools.h" - -/* Hydro part */ -const int chemistry_csds_field_size_part[chemistry_csds_field_part_count] = { - member_size(struct part, chemistry_data.smoothed_metal_mass_fraction) + - member_size(struct part, chemistry_data.metal_mass), -}; - -int chemistry_csds_local_to_global_part[chemistry_csds_field_part_count]; - -/* Stellar part */ -const int chemistry_csds_field_size_spart[chemistry_csds_field_spart_count] = { - member_size(struct spart, chemistry_data.metal_mass_fraction), -}; - -int chemistry_csds_local_to_global_spart[chemistry_csds_field_spart_count]; diff --git a/src/chemistry/GEAR/csds_chemistry.h b/src/chemistry/GEAR/csds_chemistry.h deleted file mode 100644 index bf5997f..0000000 --- a/src/chemistry/GEAR/csds_chemistry.h +++ /dev/null @@ -1,182 +0,0 @@ -/******************************************************************************* - * 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_GEAR_CSDS_CHEMISTRY_H -#define SWIFT_GEAR_CSDS_CHEMISTRY_H - -#include "../config.h" - -/* local includes */ -#include "chemistry_csds.h" -#include "csds_interpolation.h" -#include "csds_loader_io.h" -#include "csds_parameters.h" -#include "csds_python_tools.h" - -/* Index of the mask in the header mask array */ -extern int chemistry_csds_local_to_global_part[chemistry_csds_field_part_count]; -extern int - chemistry_csds_local_to_global_spart[chemistry_csds_field_spart_count]; - -/* Size for each mask */ -extern const int - chemistry_csds_field_size_part[chemistry_csds_field_part_count]; -extern const int - chemistry_csds_field_size_spart[chemistry_csds_field_spart_count]; - -/** - * @brief Create the link between the fields and their derivatives for #part. - * - * @param head The #header. - */ -__attribute__((always_inline)) INLINE static void -chemistry_csds_reader_link_derivatives_part(struct header *head) {} - -/** - * @brief Create the link between the fields and their derivatives for #spart. - * - * @param head The #header. - */ -__attribute__((always_inline)) INLINE static void -chemistry_csds_reader_link_derivatives_spart(struct header *head) {} - -/** - * @brief Interpolate a field of the #part 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 #csds_field at a time < t. - * @param after Pointer to the #csds_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 - * #chemistry_csds_fields). - * @param params The simulation's #csds_parameters. - */ -__attribute__((always_inline)) INLINE static void -chemistry_csds_interpolate_field_part(const double t_before, - const struct csds_field *restrict before, - const double t_after, - const struct csds_field *restrict after, - void *restrict output, const double t, - const int field, - const struct csds_parameters *params) { -#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 - - switch (field) { - case chemistry_csds_field_part_all: - interpolate_linear_double_ND(t_before, before, t_after, after, output, t, - 2 * GEAR_CHEMISTRY_ELEMENT_COUNT, - /* periodic */ 0, params); - break; - - default: - error("Not implemented"); - } -} - -/** - * @brief Interpolate a field of the #spart 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 #csds_field at a time < t. - * @param after Pointer to the #csds_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 - * #chemistry_csds_fields). - * @param params The simulation's #csds_parameters. - */ -__attribute__((always_inline)) INLINE static void -chemistry_csds_interpolate_field_spart(const double t_before, - const struct csds_field *restrict before, - const double t_after, - const struct csds_field *restrict after, - void *restrict output, const double t, - const int field, - const struct csds_parameters *params) { - -#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 - - switch (field) { - case chemistry_csds_field_spart_metal_mass_fractions: - interpolate_linear_double_ND(t_before, before, t_after, after, output, t, - GEAR_CHEMISTRY_ELEMENT_COUNT, - /* periodic= */ 0, params); - break; - - default: - error("Not implemented"); - } -} - -#ifdef HAVE_PYTHON -/** - * @brief Defines the different arrays for python. - */ -__attribute__((always_inline)) INLINE static void -chemistry_csds_generate_python_part(struct csds_python_field *fields) { - fields[chemistry_csds_field_part_all] = - csds_loader_python_field(/* Dimension */ 2, CUSTOM_NPY_TYPE); - - csds_loader_python_field_add_subfield( - &fields[chemistry_csds_field_part_all], "SmoothedMetalMassFractions", - /* offset */ 0, GEAR_CHEMISTRY_ELEMENT_COUNT, NPY_DOUBLE); - - csds_loader_python_field_add_subfield( - &fields[chemistry_csds_field_part_all], "MetalMassFractions", - GEAR_CHEMISTRY_ELEMENT_COUNT * sizeof(double), - GEAR_CHEMISTRY_ELEMENT_COUNT, NPY_DOUBLE); -} - -/** - * @brief Defines the different arrays for python. - */ -__attribute__((always_inline)) INLINE static void -chemistry_csds_generate_python_spart(struct csds_python_field *fields) { - fields[chemistry_csds_field_spart_metal_mass_fractions] = - csds_loader_python_field(GEAR_CHEMISTRY_ELEMENT_COUNT, NPY_DOUBLE); -} - -#endif // HAVE_PYTHON -#endif // SWIFT_GEAR_CSDS_CHEMISTRY_H diff --git a/src/chemistry/none/csds_chemistry.c b/src/chemistry/none/csds_chemistry.c deleted file mode 100644 index 98e364f..0000000 --- a/src/chemistry/none/csds_chemistry.c +++ /dev/null @@ -1,35 +0,0 @@ -/******************************************************************************* - * 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 this object's header */ -#include "csds_chemistry.h" - -/* Local headers */ -#include "csds_tools.h" - -/* Hydro part */ -const int chemistry_csds_field_size_part[chemistry_csds_field_part_count] = {}; - -int chemistry_csds_local_to_global_part[chemistry_csds_field_part_count]; - -/* Stellar part */ -const int chemistry_csds_field_size_spart[chemistry_csds_field_spart_count] = - {}; - -int chemistry_csds_local_to_global_spart[chemistry_csds_field_spart_count]; diff --git a/src/chemistry/none/csds_chemistry.h b/src/chemistry/none/csds_chemistry.h deleted file mode 100644 index 4e44f2d..0000000 --- a/src/chemistry/none/csds_chemistry.h +++ /dev/null @@ -1,124 +0,0 @@ -/******************************************************************************* - * 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_NONE_CSDS_CHEMISTRY_H -#define SWIFT_NONE_CSDS_CHEMISTRY_H - -#include "../config.h" - -/* local includes */ -#include "chemistry_csds.h" -#include "csds_interpolation.h" -#include "csds_loader_io.h" -#include "csds_parameters.h" -#include "csds_python_tools.h" - -/* Index of the mask in the header mask array */ -extern int chemistry_csds_local_to_global_part[chemistry_csds_field_part_count]; -extern int - chemistry_csds_local_to_global_spart[chemistry_csds_field_spart_count]; - -/* Size for each mask */ -extern const int - chemistry_csds_field_size_part[chemistry_csds_field_part_count]; -extern const int - chemistry_csds_field_size_spart[chemistry_csds_field_spart_count]; - -/** - * @brief Create the link between the fields and their derivatives for #part. - * - * @param head The #header. - */ -__attribute__((always_inline)) INLINE static void -chemistry_csds_reader_link_derivatives_part(struct header *head) {} - -/** - * @brief Create the link between the fields and their derivatives for #spart. - * - * @param head The #header. - */ -__attribute__((always_inline)) INLINE static void -chemistry_csds_reader_link_derivatives_spart(struct header *head) {} - -/** - * @brief Interpolate a field of the #part 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 #csds_field at a time < t. - * @param after Pointer to the #csds_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 - * #chemistry_csds_fields). - * @param params The simulation's #csds_parameters. - */ -__attribute__((always_inline)) INLINE static void -chemistry_csds_interpolate_field_part(const double t_before, - const struct csds_field *restrict before, - const double t_after, - const struct csds_field *restrict after, - void *restrict output, const double t, - const int field, - const struct csds_parameters *params) {} - -/** - * @brief Interpolate a field of the #spart 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 #csds_field at a time < t. - * @param after Pointer to the #csds_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 - * #chemistry_csds_fields). - * @param params The simulation's #csds_parameters. - */ -__attribute__((always_inline)) INLINE static void -chemistry_csds_interpolate_field_spart(const double t_before, - const struct csds_field *restrict before, - const double t_after, - const struct csds_field *restrict after, - void *restrict output, const double t, - const int field, - const struct csds_parameters *params) {} - -#ifdef HAVE_PYTHON -/** - * @brief Defines the different arrays for python. - */ -__attribute__((always_inline)) INLINE static void -chemistry_csds_generate_python_part(struct csds_python_field *fields) {} - -/** - * @brief Defines the different arrays for python. - */ -__attribute__((always_inline)) INLINE static void -chemistry_csds_generate_python_spart(struct csds_python_field *fields) {} - -#endif // HAVE_PYTHON -#endif // SWIFT_NONE_CSDS_CHEMISTRY_H diff --git a/src/csds_chemistry.h b/src/csds_chemistry.h deleted file mode 100644 index 18c7a86..0000000 --- a/src/csds_chemistry.h +++ /dev/null @@ -1,40 +0,0 @@ -/******************************************************************************* - * 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_CSDS_CHEMISTRY_H -#define SWIFT_CSDS_CHEMISTRY_H - -/* Config parameters. */ -#include "../config.h" - -/* Import the right functions */ -#if defined(CHEMISTRY_NONE) -#include "./chemistry/none/csds_chemistry.h" -#elif defined(CHEMISTRY_GEAR) -#include "./chemistry/GEAR/csds_chemistry.h" -#elif defined(CHEMISTRY_GEAR_DIFFUSION) -#error TODO -#elif defined(CHEMISTRY_QLA) -#error TODO -#elif defined(CHEMISTRY_EAGLE) -#error TODO -#else -#error "Invalid choice of chemistry function." -#endif - -#endif /* SWIFT_CHEMISTRY_H */ diff --git a/src/csds_cosmology.c b/src/csds_cosmology.c index 4f94c48..2f531d5 100644 --- a/src/csds_cosmology.c +++ b/src/csds_cosmology.c @@ -1,5 +1,5 @@ /******************************************************************************* - * This file is part of SWIFT. + * This file is part of CSDS. * Copyright (c) 2021 Loic Hausammann (loic.hausammann@epfl.ch) * * This program is free software: you can redistribute it and/or modify @@ -29,7 +29,6 @@ /* Some standard headers */ #include <math.h> - /** * @brief Add the cosmological factors for the position. * See Hausammann et al. 2021 for details. @@ -39,14 +38,15 @@ * @param vel The velocity to modify. * @param acc The acceleration to modify. */ -void csds_cosmology_add_factor_position( - const struct csds_cosmology *cosmo, - const double scale, float *vel, float *acc) { +void csds_cosmology_add_factor_position(const struct csds_cosmology *cosmo, + const double scale, float *vel, + float *acc) { const double a_dot = csds_cosmology_get_scale_factor_derivative(cosmo, scale); const double scale2 = scale * scale; if (acc) { - const double a_ddot = csds_cosmology_get_scale_factor_second_derivative(cosmo, scale); + const double a_ddot = + csds_cosmology_get_scale_factor_second_derivative(cosmo, scale); *acc = (*acc - 2 * a_dot * *vel) * scale; *acc -= scale2 * a_ddot * *vel / a_dot; *acc /= scale2 * scale2 * a_dot * a_dot; @@ -62,8 +62,10 @@ void csds_cosmology_add_factor_position( * * @return The equation of state. */ -double csds_cosmology_get_w_tilde(const struct csds_cosmology *cosmo, const double scale) { - const double w_tilde = (scale - 1.) * cosmo->w_a - (1. + cosmo->w_0 + cosmo->w_a) * log(scale); +double csds_cosmology_get_w_tilde(const struct csds_cosmology *cosmo, + const double scale) { + const double w_tilde = + (scale - 1.) * cosmo->w_a - (1. + cosmo->w_0 + cosmo->w_a) * log(scale); return exp(3. * w_tilde); } @@ -75,7 +77,8 @@ double csds_cosmology_get_w_tilde(const struct csds_cosmology *cosmo, const doub * * @return The first derivative of the scale factor. */ -double csds_cosmology_get_E(const struct csds_cosmology *cosmo, const double scale) { +double csds_cosmology_get_E(const struct csds_cosmology *cosmo, + const double scale) { const double a_inv = 1. / scale; const double a_inv2 = a_inv * a_inv; const double Omega_m = cosmo->Omega_b + cosmo->Omega_cdm; @@ -95,7 +98,8 @@ double csds_cosmology_get_E(const struct csds_cosmology *cosmo, const double sca * * @return The first derivative of the scale factor. */ -double csds_cosmology_get_scale_factor_derivative(const struct csds_cosmology *cosmo, const double scale) { +double csds_cosmology_get_scale_factor_derivative( + const struct csds_cosmology *cosmo, const double scale) { const double E = csds_cosmology_get_E(cosmo, scale); return cosmo->H0 * scale * E; } @@ -108,7 +112,8 @@ double csds_cosmology_get_scale_factor_derivative(const struct csds_cosmology *c * * @return The first derivative of the scale factor. */ -double csds_cosmology_get_scale_factor_second_derivative(const struct csds_cosmology *cosmo, const double scale) { +double csds_cosmology_get_scale_factor_second_derivative( + const struct csds_cosmology *cosmo, const double scale) { /* Compute some powers of a */ const double a_inv = 1. / scale; @@ -122,10 +127,9 @@ double csds_cosmology_get_scale_factor_second_derivative(const struct csds_cosmo double dark_eos_der = (cosmo->w_a) - (1 + cosmo->w_0 + cosmo->w_a) * a_inv; dark_eos_der *= csds_cosmology_get_w_tilde(cosmo, scale); const double Omega_m = cosmo->Omega_b + cosmo->Omega_cdm; - const double E_der = (-4 * cosmo->Omega_r * a_inv4 * a_inv - - 3 * Omega_m * a_inv4 - - 2 * cosmo->Omega_k * a_inv2 * a_inv + - dark_eos_der); + const double E_der = + (-4 * cosmo->Omega_r * a_inv4 * a_inv - 3 * Omega_m * a_inv4 - + 2 * cosmo->Omega_k * a_inv2 * a_inv + dark_eos_der); /* Compute the full expression */ const double a_dot = csds_cosmology_get_scale_factor_derivative(cosmo, scale); @@ -140,44 +144,38 @@ double csds_cosmology_get_scale_factor_second_derivative(const struct csds_cosmo * * @return The first derivative of the scale factor. */ -void csds_cosmology_init(struct csds_cosmology *cosmo, struct swift_params *params) { +void csds_cosmology_init(struct csds_cosmology *cosmo, + struct swift_params *params) { /* Read Omega_cdm */ - cosmo->Omega_cdm = parser_get_param_double( - params, "Cosmology:Omega_cdm"); + cosmo->Omega_cdm = parser_get_param_double(params, "Cosmology:Omega_cdm"); /* Read Omega_lambda */ - cosmo->Omega_lambda = parser_get_param_double( - params, "Cosmology:Omega_lambda"); + cosmo->Omega_lambda = + parser_get_param_double(params, "Cosmology:Omega_lambda"); /* Read Omega_b */ - cosmo->Omega_b = parser_get_param_double( - params, "Cosmology:Omega_b"); + cosmo->Omega_b = parser_get_param_double(params, "Cosmology:Omega_b"); /* Read Omega_r */ - cosmo->Omega_r = parser_get_param_double( - params, "Cosmology:Omega_r"); + cosmo->Omega_r = parser_get_param_double(params, "Cosmology:Omega_r"); /* Read Omega_k */ - cosmo->Omega_k = parser_get_param_double( - params, "Cosmology:Omega_k"); + cosmo->Omega_k = parser_get_param_double(params, "Cosmology:Omega_k"); /* Read Omega_nu_0 */ - double Omega_nu_0 = parser_get_opt_param_double( - params, "Cosmology:Omega_nu_0", 0); + double Omega_nu_0 = + parser_get_opt_param_double(params, "Cosmology:Omega_nu_0", 0); if (Omega_nu_0 != 0) error_python("Neutrinos are not implemented in the cosmology."); /* Read w_0 */ - cosmo->w_0 = parser_get_param_double( - params, "Cosmology:w_0"); + cosmo->w_0 = parser_get_param_double(params, "Cosmology:w_0"); /* Read w_a */ - cosmo->w_a = parser_get_param_double( - params, "Cosmology:w_a"); + cosmo->w_a = parser_get_param_double(params, "Cosmology:w_a"); /* Read the Hubble constant at z=0 */ - cosmo->H0 = parser_get_param_double( - params, "Cosmology:Hubble0"); + cosmo->H0 = parser_get_param_double(params, "Cosmology:Hubble0"); } diff --git a/src/csds_cosmology.h b/src/csds_cosmology.h index b5ba6fa..5c103cd 100644 --- a/src/csds_cosmology.h +++ b/src/csds_cosmology.h @@ -1,5 +1,5 @@ /******************************************************************************* - * This file is part of SWIFT. + * This file is part of CSDS. * Copyright (c) 2021 Loic Hausammann (loic.hausammann@epfl.ch) * * This program is free software: you can redistribute it and/or modify @@ -57,11 +57,14 @@ struct csds_cosmology { double H0; }; - -double csds_cosmology_get_scale_factor_derivative(const struct csds_cosmology *cosmo, const double scale); -double csds_cosmology_get_scale_factor_second_derivative(const struct csds_cosmology *cosmo, const double scale); -void csds_cosmology_init(struct csds_cosmology *cosmo, struct swift_params *params); +double csds_cosmology_get_scale_factor_derivative( + const struct csds_cosmology *cosmo, const double scale); +double csds_cosmology_get_scale_factor_second_derivative( + const struct csds_cosmology *cosmo, const double scale); +void csds_cosmology_init(struct csds_cosmology *cosmo, + struct swift_params *params); void csds_cosmology_add_factor_position(const struct csds_cosmology *cosmo, - const double scale, float *vel, float *acc); + const double scale, float *vel, + float *acc); #endif // CSDS_CSDS_PARAMETERS_H diff --git a/src/csds_fields.c b/src/csds_fields.c new file mode 100644 index 0000000..861b1a9 --- /dev/null +++ b/src/csds_fields.c @@ -0,0 +1,39 @@ +/******************************************************************************* + * This file is part of CSDS. + * Copyright (c) 2021 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 this file header */ +#include "csds_fields.h" + +/* Task type names. */ +const char *fields_names[task_type_count] = {SPECIAL_FLAGS_NAME, + "Coordinates", + "Velocities", + "Accelerations", + "Masses", + "SmoothingLengths", + "InternalEnergies", + "Entropies", + "ParticleIDs", + "Densities", + "SPHENIXSecondaryFields", + "BirthTimes", + "GEARStarFormation", + "GEARChemistryParts", + "GEARChemistrySparts", + "None"}; diff --git a/src/csds_fields.h b/src/csds_fields.h new file mode 100644 index 0000000..ba359ff --- /dev/null +++ b/src/csds_fields.h @@ -0,0 +1,427 @@ +/******************************************************************************* + * This file is part of CSDS. + * Copyright (c) 2021 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/>. + * + ******************************************************************************/ +/** + * @brief This file contains functions that deals with the fields + */ +#ifndef CSDS_FIELDS_H +#define CSDS_FIELDS_H + +/* Include CSDS headers */ +#include "csds_python_tools.h" +#include "csds_tools.h" + +#define SPECIAL_FLAGS_NAME "SpecialFlags" +// The fields arrays are constructed such as the special +// flag is always first (see csds_header.c) +#define SPECIAL_FLAGS_INDEX 0 +#define SPECIAL_FLAGS_SIZE sizeof(uint32_t) +#define SPECIAL_FLAGS_MASK (1 << 0) + +#define TIMESTAMP_MASK (1 << 1) +#define TIMESTAMP_NAME "Timestamp" +#define TIMESTAMP_SIZE sizeof(integertime_t) + sizeof(double) + +#ifndef GEAR_CHEMISTRY_ELEMENT_COUNT +#define GEAR_CHEMISTRY_ELEMENT_COUNT 10 +#endif + +/** + * @brief All the fields that are implemented */ +enum field_enum { + field_enum_special_flags = 0, + field_enum_coordinates, + field_enum_velocities, + field_enum_accelerations, + field_enum_masses, + field_enum_smoothing_lengths, + field_enum_internal_energies, + field_enum_entropies, + field_enum_particles_ids, + field_enum_densities, + field_enum_sphenix_secondary_fields, + field_enum_birth_time, + field_enum_gear_star_formation, + field_enum_gear_chemistry_part, + field_enum_gear_chemistry_spart, + field_enum_none, /* Should be just before count */ + field_enum_count +}; + +/** + * @brief Names of each field. + */ +extern const char *fields_names[]; + +/** + * @brief Structure defining a single field without its derivatives + */ +struct single_field { + /* The enum corresponding to the field */ + enum field_enum field; + + /* Its mask */ + unsigned int mask; + + /* Its size */ + int size; + + /* Its position within the array h->fields */ + int position; +}; + +/** + * @brief Structure containing all the information + * of a field including the derivatives. + */ +struct field_information { + /* The field */ + struct single_field field; + + /* Its first derivative */ + struct single_field first; + + /* Its second derivative */ + struct single_field second; +}; + +/** + * @brief Get the #field_enum from the field name. + */ +__attribute__((always_inline)) INLINE static enum field_enum field_enum_get( + const char *name) { + for (enum field_enum i = 0; i < field_enum_count; i++) { + if (strcmp(name, fields_names[i]) == 0) return i; + } + + error_python("Field %s not found", name); +} + +/** + * @brief Get the name of the field. + */ +__attribute__((always_inline)) INLINE static const char *field_enum_get_name( + enum field_enum field) { + return fields_names[field]; +} + +/** + * @brief Get the #field_enum for the first derivative. + */ +__attribute__((always_inline)) INLINE static enum field_enum +field_enum_get_first_derivative(enum field_enum field) { + switch (field) { + case field_enum_coordinates: + return field_enum_velocities; + case field_enum_velocities: + return field_enum_accelerations; + + default: + return field_enum_none; + } +} + +/** + * @brief Get the #field_enum for the second derivative. + */ +__attribute__((always_inline)) INLINE static enum field_enum +field_enum_get_second_derivative(enum field_enum field) { + switch (field) { + case field_enum_coordinates: + return field_enum_accelerations; + + default: + return field_enum_none; + } +} + +/** + * @brief Get the size of the field. + */ +__attribute__((always_inline)) INLINE static unsigned int field_enum_get_size( + enum field_enum field) { + switch (field) { + + case field_enum_special_flags: + return SPECIAL_FLAGS_SIZE; + + case field_enum_coordinates: + return 3 * sizeof(double); + + case field_enum_velocities: + case field_enum_accelerations: + return 3 * sizeof(float); + + case field_enum_masses: + case field_enum_smoothing_lengths: + case field_enum_internal_energies: + case field_enum_entropies: + case field_enum_birth_time: + case field_enum_densities: + return sizeof(float); + + case field_enum_particles_ids: + return sizeof(long long); + + case field_enum_sphenix_secondary_fields: + return 7 * sizeof(float); + + case field_enum_gear_star_formation: + return 2 * sizeof(float) + sizeof(long long); + + case field_enum_gear_chemistry_part: + return 2 * GEAR_CHEMISTRY_ELEMENT_COUNT * sizeof(double); + + case field_enum_gear_chemistry_spart: + return GEAR_CHEMISTRY_ELEMENT_COUNT * sizeof(double); + + default: + error_python("The field %s is not implemented", fields_names[field]); + } +} + +#ifdef HAVE_PYTHON +/** + * @brief Set a #csds_python_field for the field. + */ +__attribute__((always_inline)) INLINE static struct csds_python_field +field_enum_get_python_field(enum field_enum field) { + switch (field) { + case field_enum_coordinates: + return csds_loader_python_field(/* Dimension */ 3, NPY_DOUBLE); + + case field_enum_velocities: + case field_enum_accelerations: + return csds_loader_python_field(/* Dimension */ 3, NPY_FLOAT32); + + case field_enum_masses: + case field_enum_smoothing_lengths: + case field_enum_internal_energies: + case field_enum_entropies: + case field_enum_birth_time: + case field_enum_densities: + return csds_loader_python_field(/* Dimension */ 1, NPY_FLOAT32); + + case field_enum_particles_ids: + return csds_loader_python_field(/* Dimension */ 1, NPY_LONGLONG); + + case field_enum_sphenix_secondary_fields: { + struct csds_python_field python_field = + csds_loader_python_field(/* Dimension */ 7, CUSTOM_NPY_TYPE); + + csds_loader_python_field_add_subfield(&python_field, "Entropies", + /* offset */ 0, + /* Dimension */ 1, NPY_FLOAT32); + + csds_loader_python_field_add_subfield(&python_field, "Pressures", + /* offset */ sizeof(float), + /* Dimension */ 1, NPY_FLOAT32); + + csds_loader_python_field_add_subfield( + &python_field, "ViscosityParameters", /* offset */ 2 * sizeof(float), + /* Dimension */ 1, NPY_FLOAT32); + + csds_loader_python_field_add_subfield( + &python_field, "DiffusionParameters", /* offset */ 3 * sizeof(float), + /* Dimension */ 1, NPY_FLOAT32); + + csds_loader_python_field_add_subfield(&python_field, + "LaplacianInternalEnergies", + /* offset */ 4 * sizeof(float), + /* Dimension */ 1, NPY_FLOAT32); + + csds_loader_python_field_add_subfield( + &python_field, "VelocityDivergences", /* offset */ 5 * sizeof(float), + /* Dimension */ 1, NPY_FLOAT32); + + csds_loader_python_field_add_subfield( + &python_field, "VelocityDivergenceTimeDifferentials", + /* offset */ 6 * sizeof(float), + /* Dimension */ 1, NPY_FLOAT32); + return python_field; + } + + case field_enum_gear_star_formation: { + struct csds_python_field python_field = + csds_loader_python_field(/* Dimension */ 3, CUSTOM_NPY_TYPE); + + csds_loader_python_field_add_subfield(&python_field, "BirthDensities", + /* offset */ 0, + /* Dimension */ 1, NPY_FLOAT32); + + csds_loader_python_field_add_subfield(&python_field, "BirthMasses", + /* offset */ sizeof(float), + /* Dimension */ 1, NPY_FLOAT32); + + csds_loader_python_field_add_subfield(&python_field, "ProgenitorIDs", + /* offset */ 2 * sizeof(float), + /* Dimension */ 1, NPY_LONGLONG); + + return python_field; + } + + case field_enum_gear_chemistry_part: { + struct csds_python_field python_field = + csds_loader_python_field(/* Dimension */ 2, CUSTOM_NPY_TYPE); + + csds_loader_python_field_add_subfield( + &python_field, "SmoothedMetalMassFractions", /* offset */ 0, + GEAR_CHEMISTRY_ELEMENT_COUNT, NPY_DOUBLE); + + csds_loader_python_field_add_subfield( + &python_field, "MetalMassFractions", + GEAR_CHEMISTRY_ELEMENT_COUNT * sizeof(double), + GEAR_CHEMISTRY_ELEMENT_COUNT, NPY_DOUBLE); + return python_field; + } + + case field_enum_gear_chemistry_spart: + return csds_loader_python_field(GEAR_CHEMISTRY_ELEMENT_COUNT, NPY_DOUBLE); + + default: + error_python("The field %s is not implemented", fields_names[field]); + } +} +#endif + +struct mask_data { + /* Number of bytes for a mask. */ + int size; + + /* Mask value. */ + unsigned int mask; + + /* Name of the mask. */ + char name[STRING_SIZE]; +}; + +/** + * @brief Initialize a derivative of a #field_information. + * + * @param der The derivative to initialize. + * @param field The field corresponding to the derivative. + * @param all_data All the #mask_data from the header. + * @param size_all_data Size of all_data. + */ +__attribute__((always_inline)) INLINE static void +field_information_set_derivative(struct single_field *der, + enum field_enum field, + struct mask_data *all_data, + int size_all_data) { + + der->field = field; + der->mask = 0; + der->size = 0; + der->position = -1; + /* Set the first derivative */ + if (der->field != field_enum_none) { + struct mask_data *data = NULL; + const char *name_first = field_enum_get_name(der->field); + /* Get the first derivative in all_data */ + for (int i = 0; i < size_all_data; i++) { + if (strcmp(name_first, all_data[i].name) == 0) { + data = &all_data[i]; + break; + } + } + /* The first derivative does not exist in the header. */ + if (data == NULL) { + der->field = field_enum_none; + return; + } + + der->mask = data->mask; + der->size = field_enum_get_size(der->field); + if (der->size != data->size) { + error_python("Inconsistent size for %s (%i != %i)", fields_names[field], + der->size, data->size); + } + } +} + +/** + * @brief Initialize a #field_information. + * + * @param info The #field_information to initialize. + * @param data The #mask_data of the field. + * @param all_data All the #mask_data from the header. + * @param size_all_data Size of all_data. + */ +__attribute__((always_inline)) INLINE static void field_information_init( + struct field_information *info, struct mask_data *data, + struct mask_data *all_data, int size_all_data) { + + /* Set the field */ + enum field_enum field = field_enum_get(data->name); + info->field.field = field; + info->field.mask = data->mask; + info->field.size = field_enum_get_size(field); + if (info->field.size != data->size) { + error_python("Inconsistent size for %s (%i != %i)", fields_names[field], + info->field.size, data->size); + } + info->field.position = -1; + + /* Get the derivatives */ + enum field_enum first = field_enum_get_first_derivative(field); + enum field_enum second = field_enum_get_second_derivative(field); + + /* Set the two derivatives */ + field_information_set_derivative(&info->first, first, all_data, + size_all_data); + field_information_set_derivative(&info->second, second, all_data, + size_all_data); +} + +/** + * @brief Set the position of the field and its derivatives + * according to the array. + * + * @param info An array of #field_information. + * @param size The size of the array. + */ +__attribute__((always_inline)) INLINE static void +field_information_set_positions(struct field_information *info, + const unsigned int size) { + for (unsigned int i = 0; i < size; i++) { + /* Set the field */ + info[i].field.position = i; + + /* Set the first derivative */ + if (info[i].first.field != field_enum_none) { + for (unsigned int j = 0; j < size; j++) { + if (info[i].first.field == info[j].field.field) { + info[i].first.position = j; + break; + } + } + } + + /* Set the second derivative */ + if (info[i].second.field != field_enum_none) { + for (unsigned int j = 0; j < size; j++) { + if (info[i].second.field == info[j].field.field) { + info[i].second.position = j; + break; + } + } + } + } +} + +#endif // CSDS_FIELDS_H diff --git a/src/csds_gravity.h b/src/csds_gravity.h deleted file mode 100644 index bc4f477..0000000 --- a/src/csds_gravity.h +++ /dev/null @@ -1,38 +0,0 @@ -/******************************************************************************* - * 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_CSDS_GRAVITY_H -#define SWIFT_CSDS_GRAVITY_H - -/* Config parameters. */ -#include "../config.h" - -/* Import the right functions */ -#if defined(DEFAULT_GRAVITY) -#error TODO -#include "./gravity/Default/csds_gravity.h" -#elif defined(POTENTIAL_GRAVITY) -#error TODO -#include "./gravity/Potential/csds_gravity.h" -#elif defined(MULTI_SOFTENING_GRAVITY) -#include "./gravity/MultiSoftening/csds_gravity.h" -#else -#error "Invalid choice of gravity variant" -#endif - -#endif diff --git a/src/csds_header.c b/src/csds_header.c index e75c313..8dc6042 100644 --- a/src/csds_header.c +++ b/src/csds_header.c @@ -1,5 +1,5 @@ /******************************************************************************* - * This file is part of SWIFT. + * This file is part of CSDS. * Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch) * * This program is free software: you can redistribute it and/or modify @@ -45,13 +45,28 @@ void header_print(const struct header *h) { #endif message("First Offset: %lu.", h->offset_first_record); message("Offset direction: %s.", csds_offset_name[h->offset_direction]); - message("Number masks: %i.", h->masks_count); - for (int i = 0; i < h->masks_count; i++) { - message(" Mask: %s.", h->masks[i].name); - message(" Value: %u.", h->masks[i].mask); - message(" Size: %i.", h->masks[i].size); - message(""); + for (int type = 0; type < swift_type_count; type++) { + message("Masks for particle type %i:", type); + for (int k = 0; k < h->number_fields[type]; k++) { + /* Do the field */ + struct field_information *field = &h->fields[type][k]; + const char *name = field_enum_get_name(field->field.field); + message(" %20s:\t mask=%03u\t size=%i", name, field->field.mask, + field->field.size); + + /* Do the first derivative */ + if (field->first.field != field_enum_none) { + const char *name_first = field_enum_get_name(field->first.field); + message(" First derivative: %s", name_first); + } + + /* Do the second derivative */ + if (field->second.field != field_enum_none) { + const char *name_second = field_enum_get_name(field->second.field); + message(" Second derivative: %s", name_second); + } + } } }; @@ -60,23 +75,13 @@ void header_print(const struct header *h) { * * @param h The #header. */ -void header_free(struct header *h) { free(h->masks); }; - -/** - * @brief Check if a field is present in the header. - * - * @param h The #header. - * @param field name of the requested field. - * @return Index of the field (-1 if not found). - */ -int header_get_field_index(const struct header *h, const char *field) { - for (int i = 0; i < h->masks_count; i++) { - if (strcmp(h->masks[i].name, field) == 0) { - return i; +void header_free(struct header *h) { + for (int i = 0; i < swift_type_count; i++) { + if (h->fields[i]) { + free(h->fields[i]); + h->fields[i] = NULL; } } - - return -1; }; /** @@ -97,47 +102,6 @@ void header_change_offset_direction(struct header *h, &new_value); } -/** - * @brief Set the links between the local (module) and global (header) indexes. - */ -#define set_links_local_global_internal(MODULE, PART) \ - for (int j = 0; j < MODULE##_csds_field##PART##_count; j++) { \ - MODULE##_csds_local_to_global##PART[j] = -1; \ - for (int i = 0; i < h->masks_count; i++) { \ - if (strcmp(h->masks[i].name, MODULE##_csds_field_names##PART[j]) == 0) { \ - MODULE##_csds_local_to_global##PART[j] = i; \ - break; \ - } \ - } \ - \ - /* Check if everything is fine. */ \ - const int index = MODULE##_csds_local_to_global##PART[j]; \ - if (index == -1) { \ - error("Field %s in " #MODULE " is not set", \ - MODULE##_csds_field_names##PART[j]); \ - } \ - if (h->masks[index].size != MODULE##_csds_field_size##PART[j]) { \ - error("Field %s in " #MODULE " does not have the correct size", \ - MODULE##_csds_field_names##PART[j]); \ - } \ - } \ - \ - MODULE##_csds_reader_link_derivatives##PART(h); - -/** - * Same function as set_links_local_global_internal before but with a single - * argument. - */ -#define set_links_local_global_single_particle_type(MODULE) \ - set_links_local_global_internal(MODULE, ) - -/** - * Same function as set_links_local_global_internal before but with a cleaner - * argument. - */ -#define set_links_local_global(MODULE, PART) \ - set_links_local_global_internal(MODULE, _##PART) - /** * @brief read the csds header. * @@ -190,51 +154,111 @@ void header_read(struct header *h, struct csds_logfile *log) { } /* Read the number of masks. */ - map = csds_loader_io_read_data(map, sizeof(unsigned int), &h->masks_count); + unsigned int masks_count = 0; + map = csds_loader_io_read_data(map, sizeof(unsigned int), &masks_count); /* Allocate the masks memory. */ - h->masks = malloc(sizeof(struct mask_data) * h->masks_count); - if (h->masks == NULL) { + struct mask_data *masks = malloc(sizeof(struct mask_data) * masks_count); + if (masks == NULL) { error_python("Failed to allocate the memory for the masks."); } /* Loop over all masks. */ - h->timestamp_mask = 0; - for (int i = 0; i < h->masks_count; i++) { + h->mask_max = 0; + for (unsigned int i = 0; i < masks_count; i++) { /* Read the mask name. */ - map = csds_loader_io_read_data(map, h->string_length, h->masks[i].name); + map = csds_loader_io_read_data(map, h->string_length, masks[i].name); /* Set the mask value. */ - h->masks[i].mask = 1 << i; + masks[i].mask = 1 << i; + + /* Set the maximal mask */ + h->mask_max = max(2 * masks[i].mask, h->mask_max); /* Read the mask data size. */ - map = - csds_loader_io_read_data(map, sizeof(unsigned int), &h->masks[i].size); + map = csds_loader_io_read_data(map, sizeof(unsigned int), &masks[i].size); /* Print the information. */ if (reader->verbose > 1) { - message("Field %s found in the logfile", h->masks[i].name); + message("Field %s found in the logfile", masks[i].name); } + } - /* Keep the timestamp mask in memory */ - if (strcmp(h->masks[i].name, "Timestamp") == 0) { - h->timestamp_mask = h->masks[i].mask; - } + /* Check that the special flag is at the expected place. */ + const int special_flag = 0; + if (strcmp(masks[special_flag].name, SPECIAL_FLAGS_NAME) != 0) { + error_python("The special flag should be the first mask '%s'", + masks[special_flag].name); + } + + /* Check the special flag mask */ + if (masks[special_flag].mask != SPECIAL_FLAGS_MASK) { + error_python("The mask of the special flag should be %i", + SPECIAL_FLAGS_MASK); } - /* Check that the timestamp mask exists */ - if (h->timestamp_mask == 0) { - error_python("Unable to find the timestamp mask."); + /* Check the size of the special flag. */ + if (masks[special_flag].size != SPECIAL_FLAGS_SIZE) { + error_python("The special flag is expected to be %li bytes", + SPECIAL_FLAGS_SIZE); } - /* check the special flag. */ - if (strcmp(h->masks[csds_index_special_flags].name, "SpecialFlags") != 0) { - error("The special flag is expected to be at position %i", - csds_index_special_flags); + /* Check that the timestamp is at the expected place */ + const int timestamp = 1; + if (strcmp(masks[timestamp].name, TIMESTAMP_NAME) != 0) { + error_python("The time stamp should be the second mask"); } - if (h->masks[csds_index_special_flags].size != sizeof(uint32_t)) { - error("The special flag is expected to be 32 bits"); + /* Check the timestamp mask */ + if (masks[timestamp].mask != TIMESTAMP_MASK) { + error_python("The mask of the timestamp should be %i", TIMESTAMP_MASK); + } + + /* Check the size of the timestamp. */ + if (masks[timestamp].size != TIMESTAMP_SIZE) { + error_python( + "The timestamp is expected to be an integertime_t and a double."); + } + + /* Read the number of fields per particle */ + map = + csds_loader_io_read_data(map, sizeof(h->number_fields), h->number_fields); + + /* Read the order of the fields */ + for (int type = 0; type < swift_type_count; type++) { + if (h->number_fields[type] == 0) { + h->fields[type] = NULL; + continue; + } + + /* Allocate and read the order */ + size_t size = h->number_fields[type] * sizeof(int); + int *order = (int *)malloc(size); + map = csds_loader_io_read_data(map, size, order); + + /* Add the special flag */ + size = (h->number_fields[type] + 1) * sizeof(struct field_information); + + /* Allocate the memory in the header */ + h->fields[type] = (struct field_information *)malloc(size); + + /* Set the fields */ + for (int k = 0; k < h->number_fields[type]; k++) { + /* Special flag is at 0, thus skip it */ + field_information_init(&h->fields[type][k + 1], &masks[order[k]], masks, + masks_count); + } + + /* Add the contribution from the special flag */ + field_information_init(&h->fields[type][0], &masks[SPECIAL_FLAGS_INDEX], + masks, masks_count); + h->number_fields[type]++; + + /* Set the positions */ + field_information_set_positions(h->fields[type], h->number_fields[type]); + + /* Free the memory */ + free(order); } /* Check the logfile header's size. */ @@ -244,19 +268,8 @@ void header_read(struct header *h, struct csds_logfile *log) { h->offset_first_record, offset); } - /* Set the first and second derivatives as non existent. */ - for (int i = 0; i < h->masks_count; i++) { - h->masks[i].reader.first_deriv = -1; - h->masks[i].reader.second_deriv = -1; - } - - /* Set the link between local and global */ - set_links_local_global_single_particle_type(hydro); - set_links_local_global_single_particle_type(gravity); - set_links_local_global_single_particle_type(stars); - set_links_local_global(chemistry, part); - set_links_local_global(chemistry, spart); - set_links_local_global_single_particle_type(star_formation); + /* Cleanup the memory */ + free(masks); }; /** @@ -267,14 +280,67 @@ void header_read(struct header *h, struct csds_logfile *log) { * * @return number of bits in mask. */ -size_t header_get_record_size_from_mask(const struct header *h, - const size_t mask) { +size_t header_get_record_size_from_mask(const struct header *h, size_t mask) { size_t count = 0; + + /* Do the time stamp */ + if (mask & TIMESTAMP_MASK) { + count += TIMESTAMP_SIZE; + mask &= ~TIMESTAMP_MASK; + } + /* Loop over each masks. */ - for (int i = 0; i < h->masks_count; i++) { - if (mask & h->masks[i].mask) { - count += h->masks[i].size; + for (int part_type = 0; part_type < swift_type_count; part_type++) { + for (int i = 0; i < h->number_fields[part_type]; i++) { + if (mask & h->fields[part_type][i].field.mask) { + count += h->fields[part_type][i].field.size; + mask &= ~h->fields[part_type][i].field.mask; + + /* Early exit */ + if (mask == 0) { + part_type = swift_type_count; + break; + } + } } } + + /* Ensure that we found all the masks */ + if (mask != 0) { + error("A mask was not found"); + } return count; } + +/** + * @brief Find the #field_information from the field name. + * + * @param h The #header + * @param name The name of the field. + * @param part_type The particle type. + * + * @return The first #field_information found. + */ +const struct field_information *header_get_field_from_name( + const struct header *h, const char *name, enum part_type part_type) { + + const struct field_information *info = NULL; + for (int i = 0; i < h->number_fields[part_type]; i++) { + const char *current_name = + field_enum_get_name(h->fields[part_type][i].field.field); + if (strcmp(name, current_name) == 0) { + info = &h->fields[part_type][i]; + break; + } + } + + /* Check if we obtained a field. */ + if (info == NULL) { + error( + "The logfile does not contain the field %s" + "for particle type %i", + name, part_type); + } + + return info; +} diff --git a/src/csds_header.h b/src/csds_header.h index 40a2eae..e1e9ed7 100644 --- a/src/csds_header.h +++ b/src/csds_header.h @@ -1,5 +1,5 @@ /******************************************************************************* - * This file is part of SWIFT. + * This file is part of CSDS. * Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch) * * This program is free software: you can redistribute it and/or modify @@ -19,6 +19,7 @@ #ifndef CSDS_CSDS_HEADER_H #define CSDS_CSDS_HEADER_H +#include "csds_fields.h" #include "csds_tools.h" #include <stdio.h> @@ -71,30 +72,31 @@ struct header { /* Number of bytes for strings. */ unsigned int string_length; - /* Number of masks. */ - int masks_count; - - /* List of masks. */ - struct mask_data *masks; - /* Direction of the offset in the records. */ enum csds_offset_direction offset_direction; /* The corresponding log. */ struct csds_logfile *log; - /* Timestamp mask */ - size_t timestamp_mask; + /* Largest available mask (for debugging) */ + unsigned int mask_max; + + /* Number of mask per particle type */ + int number_fields[swift_type_count]; + + /* The fields for each particle type in the correct order. + By construction, the special flag is the first element in each array. */ + struct field_information *fields[swift_type_count]; }; void header_print(const struct header *h); void header_free(struct header *h); -int header_get_field_index(const struct header *h, const char *field); void header_read(struct header *h, struct csds_logfile *log); -size_t header_get_record_size_from_mask(const struct header *h, - const size_t mask); +size_t header_get_record_size_from_mask(const struct header *h, size_t mask); void header_change_offset_direction(struct header *h, enum csds_offset_direction new_value); +const struct field_information *header_get_field_from_name( + const struct header *h, const char *name, enum part_type part_type); /** * @brief Check if the offset are forward. @@ -123,40 +125,4 @@ __attribute__((always_inline)) INLINE static int header_is_corrupted( return h->offset_direction == csds_offset_corrupted; } -/** - * @brief Set the first derivative of a field. - * - * @param h The #header. - * @param field The index of the field. - * @param field_first The index of the first derivative of the field. - */ -__attribute__((always_inline)) INLINE static void header_set_first_derivative( - struct header *h, const int field, const int field_first) { - if (h->masks[field].reader.first_deriv != -1 && - h->masks[field].reader.first_deriv != field_first) { - error( - "The first derivative is not consistent with the one from a different " - "module."); - } - h->masks[field].reader.first_deriv = field_first; -} - -/** - * @brief Set the second derivative of a field. - * - * @param h The #header. - * @param field The index of the field. - * @param field_second The index of the second derivative of the field. - */ -__attribute__((always_inline)) INLINE static void header_set_second_derivative( - struct header *h, const int field, const int field_second) { - if (h->masks[field].reader.second_deriv != -1 && - h->masks[field].reader.second_deriv != field_second) { - error( - "The second derivative is not consistent with the one from a different " - "module."); - } - h->masks[field].reader.second_deriv = field_second; -} - #endif // CSDS_CSDS_HEADER_H diff --git a/src/csds_hydro.h b/src/csds_hydro.h deleted file mode 100644 index a4a1c15..0000000 --- a/src/csds_hydro.h +++ /dev/null @@ -1,61 +0,0 @@ -/******************************************************************************* - * 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_CSDS_HYDRO_H -#define SWIFT_CSDS_HYDRO_H - -/* Config parameters. */ -#include "../config.h" - -/* Import the right functions */ -#if defined(MINIMAL_SPH) -#error TODO -#include "./hydro/Minimal/csds_hydro.h" -#elif defined(GADGET2_SPH) -#include "./hydro/Gadget2/csds_hydro.h" -#elif defined(HOPKINS_PE_SPH) -#error TODO -#include "./hydro/PressureEntropy/csds_hydro.h" -#elif defined(HOPKINS_PU_SPH) -#error TODO -#include "./hydro/PressureEnergy/csds_hydro.h" -#elif defined(HOPKINS_PU_SPH_MONAGHAN) -#error TODO -#include "./hydro/PressureEnergyMorrisMonaghanAV/csds_hydro.h" -#elif defined(PHANTOM_SPH) -#error TODO -#include "./hydro/Phantom/csds_hydro.h" -#elif defined(GIZMO_MFV_SPH) || defined(GIZMO_MFM_SPH) -#error TODO -#include "./hydro/Gizmo/csds_hydro.h" -#elif defined(SHADOWFAX_SPH) -#error TODO -#include "./hydro/Shadowswift/csds_hydro.h" -#elif defined(PLANETARY_SPH) -#error TODO -#include "./hydro/Planetary/csds_hydro.h" -#elif defined(SPHENIX_SPH) -#include "./hydro/SPHENIX/csds_hydro.h" -#elif defined(ANARCHY_PU_SPH) -#error TODO -#include "./hydro/AnarchyPU/csds_hydro.h" -#else -#error "Invalid choice of SPH variant" -#endif - -#endif /* SWIFT_HYDRO_H */ diff --git a/src/csds_index.c b/src/csds_index.c index 628d92d..fa07793 100644 --- a/src/csds_index.c +++ b/src/csds_index.c @@ -1,5 +1,5 @@ /******************************************************************************* - * This file is part of SWIFT. + * This file is part of CSDS. * Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch) * * This program is free software: you can redistribute it and/or modify diff --git a/src/csds_index.h b/src/csds_index.h index 128cc2b..ba7c1df 100644 --- a/src/csds_index.h +++ b/src/csds_index.h @@ -1,5 +1,5 @@ /******************************************************************************* - * This file is part of SWIFT. + * This file is part of CSDS. * Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch) * * This program is free software: you can redistribute it and/or modify diff --git a/src/csds_interpolation.h b/src/csds_interpolation.h index 5499616..a51e83e 100644 --- a/src/csds_interpolation.h +++ b/src/csds_interpolation.h @@ -1,5 +1,5 @@ /******************************************************************************* - * This file is part of SWIFT. + * This file is part of CSDS. * Copyright (c) 2020 Loic Hausammann (loic.hausammann@epfl.ch) * * This program is free software: you can redistribute it and/or modify @@ -146,8 +146,8 @@ interpolate_cubic_hermite_spline(const double t0, const float v0, * * The field is in double and both derivatives in float. * - * @param before Pointer to the #csds_field at a time < t. - * @param after Pointer to the #csds_field at a time > t. + * @param before Pointer to the #csds_reader_field at a time < t. + * @param after Pointer to the #csds_reader_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). @@ -156,13 +156,11 @@ interpolate_cubic_hermite_spline(const double t0, const float v0, * @param params The simulation's #csds_parameters. */ __attribute__((always_inline)) INLINE static void -interpolate_quintic_double_float_ND(const double t_before, - const struct csds_field *restrict before, - const double t_after, - const struct csds_field *restrict after, - void *restrict output, const double t, - const int dimension, int periodic, - const struct csds_parameters *params) { +interpolate_quintic_double_float_ND( + const double t_before, const struct csds_reader_field *restrict before, + const double t_after, const struct csds_reader_field *restrict after, + void *restrict output, const double t, const int dimension, int periodic, + const struct csds_parameters *params) { /* Get the arrays */ const double *x_bef = (double *)before->field; @@ -187,10 +185,10 @@ interpolate_quintic_double_float_ND(const double t_before, /* Add cosmological factor for velocity and acceleration */ if (params->with_cosmology && v_bef && v_aft) { - csds_cosmology_add_factor_position(¶ms->cosmo, t_before, - &v0, a_bef ? &a0 : NULL); - csds_cosmology_add_factor_position(¶ms->cosmo, t_after, - &v1, a_aft ? &a1 : NULL); + csds_cosmology_add_factor_position(¶ms->cosmo, t_before, &v0, + a_bef ? &a0 : NULL); + csds_cosmology_add_factor_position(¶ms->cosmo, t_after, &v1, + a_aft ? &a1 : NULL); } /* Apply the periodic boundaries */ @@ -245,8 +243,8 @@ interpolate_quintic_double_float_ND(const double t_before, * * The field and the first derivatives are floats. * - * @param before Pointer to the #csds_field at a time < t. - * @param after Pointer to the #csds_field at a time > t. + * @param before Pointer to the #csds_reader_field at a time < t. + * @param after Pointer to the #csds_reader_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). @@ -255,8 +253,8 @@ interpolate_quintic_double_float_ND(const double t_before, * @param params The simulation's #csds_parameters. */ __attribute__((always_inline)) INLINE static void interpolate_cubic_float_ND( - const double t_before, const struct csds_field *restrict before, - const double t_after, const struct csds_field *restrict after, + const double t_before, const struct csds_reader_field *restrict before, + const double t_after, const struct csds_reader_field *restrict after, void *restrict output, const double t, const int dimension, int periodic, const struct csds_parameters *params) { @@ -296,8 +294,8 @@ __attribute__((always_inline)) INLINE static void interpolate_cubic_float_ND( * * The field is in float. * - * @param before Pointer to the #csds_field at a time < t. - * @param after Pointer to the #csds_field at a time > t. + * @param before Pointer to the #csds_reader_field at a time < t. + * @param after Pointer to the #csds_reader_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). @@ -306,8 +304,8 @@ __attribute__((always_inline)) INLINE static void interpolate_cubic_float_ND( * @param params The simulation's #csds_parameters. */ __attribute__((always_inline)) INLINE static void interpolate_linear_float_ND( - const double t_before, const struct csds_field *restrict before, - const double t_after, const struct csds_field *restrict after, + const double t_before, const struct csds_reader_field *restrict before, + const double t_after, const struct csds_reader_field *restrict after, void *restrict output, const double t, const int dimension, int periodic, const struct csds_parameters *params) { @@ -336,8 +334,8 @@ __attribute__((always_inline)) INLINE static void interpolate_linear_float_ND( * * The field is in double. * - * @param before Pointer to the #csds_field at a time < t. - * @param after Pointer to the #csds_field at a time > t. + * @param before Pointer to the #csds_reader_field at a time < t. + * @param after Pointer to the #csds_reader_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). @@ -346,8 +344,8 @@ __attribute__((always_inline)) INLINE static void interpolate_linear_float_ND( * @param params The simulation's #csds_parameters. */ __attribute__((always_inline)) INLINE static void interpolate_linear_double_ND( - const double t_before, const struct csds_field *restrict before, - const double t_after, const struct csds_field *restrict after, + const double t_before, const struct csds_reader_field *restrict before, + const double t_after, const struct csds_reader_field *restrict after, void *restrict output, const double t, const int dimension, int periodic, const struct csds_parameters *params) { @@ -375,8 +373,8 @@ __attribute__((always_inline)) INLINE static void interpolate_linear_double_ND( /** * @brief Interpolates a field stored as a float. * - * @param before Pointer to the #csds_field at a time < t. - * @param after Pointer to the #csds_field at a time > t. + * @param before Pointer to the #csds_reader_field at a time < t. + * @param after Pointer to the #csds_reader_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). @@ -385,8 +383,8 @@ __attribute__((always_inline)) INLINE static void interpolate_linear_double_ND( * @param params The simulation's #csds_parameters. */ __attribute__((always_inline)) INLINE static void interpolate_linear_float( - const double t_before, const struct csds_field *restrict before, - const double t_after, const struct csds_field *restrict after, + const double t_before, const struct csds_reader_field *restrict before, + const double t_after, const struct csds_reader_field *restrict after, void *restrict output, const double t, int periodic, const struct csds_parameters *params) { @@ -404,27 +402,4 @@ __attribute__((always_inline)) INLINE static void interpolate_linear_float( wa * ((float *)after->field)[0] + wb * ((float *)before->field)[0]; } -/** - * @brief Interpolation function for the ids. - * - * As the IDs should not change during the simulation, we ensure - * that the IDs are the same between two records. - * - * @param before Pointer to the #csds_field at a time < t. - * @param after Pointer to the #csds_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. - */ -__attribute__((always_inline)) INLINE static void interpolate_ids( - const double t_before, const struct csds_field *restrict before, - const double t_after, const struct csds_field *restrict after, - void *restrict output, const double t) { - if (*(long long *)after->field != *(long long *)before->field) { - error("Interpolating different particles"); - } - *(long long *)output = *(long long *)after->field; -} - #endif // CSDS_CSDS_INTERPOLATION_H diff --git a/src/csds_loader_io.c b/src/csds_loader_io.c index d6874bb..4281767 100644 --- a/src/csds_loader_io.c +++ b/src/csds_loader_io.c @@ -1,5 +1,5 @@ /******************************************************************************* - * This file is part of SWIFT. + * This file is part of CSDS. * Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch) * * This program is free software: you can redistribute it and/or modify diff --git a/src/csds_loader_io.h b/src/csds_loader_io.h index 56ba179..3022de9 100644 --- a/src/csds_loader_io.h +++ b/src/csds_loader_io.h @@ -1,5 +1,5 @@ /******************************************************************************* - * This file is part of SWIFT. + * This file is part of CSDS. * Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch) * * This program is free software: you can redistribute it and/or modify diff --git a/src/csds_logfile.c b/src/csds_logfile.c index 39541b5..f369391 100644 --- a/src/csds_logfile.c +++ b/src/csds_logfile.c @@ -1,5 +1,5 @@ /******************************************************************************* - * This file is part of SWIFT. + * This file is part of CSDS. * Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch) * * This program is free software: you can redistribute it and/or modify diff --git a/src/csds_logfile.h b/src/csds_logfile.h index 6e4c6ca..b48d49b 100644 --- a/src/csds_logfile.h +++ b/src/csds_logfile.h @@ -1,5 +1,5 @@ /******************************************************************************* - * This file is part of SWIFT. + * This file is part of CSDS. * Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch) * * This program is free software: you can redistribute it and/or modify diff --git a/src/csds_parameters.c b/src/csds_parameters.c index bc40493..e33a680 100644 --- a/src/csds_parameters.c +++ b/src/csds_parameters.c @@ -1,5 +1,5 @@ /******************************************************************************* - * This file is part of SWIFT. + * This file is part of CSDS. * Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch) * * This program is free software: you can redistribute it and/or modify diff --git a/src/csds_parameters.h b/src/csds_parameters.h index d65c137..431fcf5 100644 --- a/src/csds_parameters.h +++ b/src/csds_parameters.h @@ -1,5 +1,5 @@ /******************************************************************************* - * This file is part of SWIFT. + * This file is part of CSDS. * Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch) * * This program is free software: you can redistribute it and/or modify @@ -31,8 +31,6 @@ /* SWIFT include */ #include "part_type.h" - - struct csds_reader; /** diff --git a/src/csds_particle.c b/src/csds_particle.c index d63c9d1..2f4cedf 100644 --- a/src/csds_particle.c +++ b/src/csds_particle.c @@ -1,5 +1,5 @@ /******************************************************************************* - * This file is part of SWIFT. + * This file is part of CSDS. * Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch) * * This program is free software: you can redistribute it and/or modify @@ -19,6 +19,7 @@ #include "csds_particle.h" #include "csds_header.h" +#include "csds_interpolation.h" #include "csds_loader_io.h" #include "csds_reader.h" #include "csds_time.h" @@ -33,10 +34,10 @@ * * @param reader The #csds_reader. * @param offset offset of the record to read. - * @param output Buffer for the requested field.. - * @param all_fields The list of fields for the particle type. - * @param all_fields_count The number of element in all_fields. - * @param global_wanted The global index of the requested field. + * @param output Buffer for the requested field. + * @param part_type The particle type. + * @param field_read The field to read. + * @param derivative The derivative to read (0 for the field). * @param mask (out) The mask of the record. * @param h_offset (out) Difference of offset with the next record. * @@ -44,11 +45,12 @@ */ __attribute__((always_inline)) INLINE size_t csds_particle_read_field( const struct csds_reader *reader, size_t offset, void *output, - const struct field_information *all_fields, const int all_fields_count, - const int global_wanted, size_t *mask, size_t *h_offset) { + const enum part_type part_type, const struct field_information *field_read, + const int derivative, size_t *mask, size_t *h_offset) { /* Get a few pointers. */ const struct header *h = &reader->log.header; + const struct field_information *fields = h->fields[part_type]; char *map = reader->log.log.map; *mask = 0; @@ -58,35 +60,51 @@ __attribute__((always_inline)) INLINE size_t csds_particle_read_field( map = csds_loader_io_read_mask(h, map + offset, mask, h_offset); /* Check that the mask is within the limits. */ - if (*mask > (unsigned int)(1 << h->masks_count)) { + if (*mask > h->mask_max) { error_python("Found an unexpected mask %zi", *mask); } /* Check if it is not a time record. */ - if (*mask == h->timestamp_mask) { + if (*mask == TIMESTAMP_MASK) { error_python("Cannot read a particle from timestep record."); } - /* Skip the special field. */ - if (*mask & h->masks[csds_index_special_flags].mask) { - map += h->masks[csds_index_special_flags].size; +#ifdef SWIFT_DEBUG_CHECKS + if (derivative == 1 && field_read->first.field == field_enum_none) { + error("Trying to read the non existing first derivative."); } + if (derivative == 2 && field_read->second.field == field_enum_none) { + error("Trying to read the non existing second derivative."); + } +#endif + + /* Get the position of the field to read */ + const int position = derivative == 0 + ? field_read->field.position + : derivative == 1 ? field_read->first.position + : field_read->second.position; + +#ifdef SWIFT_DEBUG_CHECKS + if (position < 0) { + error("Position not set (derivative %i of %s)", derivative, + field_enum_get_name(field_read->field.field)); + } +#endif /* Read the record and copy it to a particle. */ - for (int i = 0; i < all_fields_count; i++) { - const int global = all_fields[i].global_index; + for (int i = 0; i < h->number_fields[part_type]; i++) { /* Is the mask present? */ - if (!(*mask & h->masks[global].mask)) { + if (!(*mask & fields[i].field.mask)) { continue; } - if (global_wanted == global) { + if (position == i) { /* Read the data. */ - map = csds_loader_io_read_data(map, h->masks[global].size, output); + map = csds_loader_io_read_data(map, fields[i].field.size, output); } else { /* Update the buffer's position. */ - map += h->masks[global].size; + map += fields[i].field.size; } } @@ -121,19 +139,18 @@ csds_particle_read_special_flag(const struct csds_reader *reader, size_t offset, map = csds_loader_io_read_mask(h, map + offset, mask, h_offset); /* Check that the mask is within the limits. */ - if (*mask > (unsigned int)(1 << h->masks_count)) { + if (*mask > h->mask_max) { error_python("Found an unexpected mask %zi", *mask); } /* Check if it is not a time record. */ - if (*mask == h->timestamp_mask) { + if (*mask == TIMESTAMP_MASK) { error_python("Cannot read a particle from timestep record."); } /* Read the special flag */ uint32_t packed_data = 0; - map = csds_loader_io_read_data(map, h->masks[csds_index_special_flags].size, - &packed_data); + map = csds_loader_io_read_data(map, SPECIAL_FLAGS_SIZE, &packed_data); return csds_unpack_flags_and_data(packed_data, data, part_type); } @@ -141,8 +158,8 @@ csds_particle_read_special_flag(const struct csds_reader *reader, size_t offset, /** * @brief Interpolate a field of the particle at the given time. * - * @param before Pointer to the #csds_field at a time < t. - * @param after Pointer to the #csds_field at a time > t. + * @param before Pointer to the #csds_reader_field at a time < t. + * @param after Pointer to the #csds_reader_field at a time > t. * @param otuput Pointer to the output value. * @param time_before Time of field_before (< t). * @param time_after Time of field_after (> t). @@ -150,76 +167,128 @@ csds_particle_read_special_flag(const struct csds_reader *reader, size_t offset, * @param field The field to reconstruct. * @param params The simulation's #csds_parameters. */ -void csds_particle_interpolate_field(const double time_before, - const struct csds_field *restrict before, - const double time_after, - const struct csds_field *restrict after, - void *restrict output, const double time, - const struct field_information *field, - enum part_type type, - const struct csds_parameters *params) { - - /* Select the correct interpolation */ - switch (type) { - - /* Hydro */ - case swift_type_gas: - switch (field->module) { - case field_module_default: - hydro_csds_interpolate_field(time_before, before, time_after, after, - output, time, field->local_index, - params); - break; - - case field_module_chemistry: - chemistry_csds_interpolate_field_part(time_before, before, time_after, - after, output, time, - field->local_index, params); - break; - - default: - error("Module not implemented"); +void csds_particle_interpolate_field( + const double time_before, const struct csds_reader_field *restrict before, + const double time_after, const struct csds_reader_field *restrict after, + void *restrict output, const double time, + const struct field_information *field, enum part_type type, + const struct csds_parameters *params) { + + switch (field->field.field) { + /* Coordinates */ + case field_enum_coordinates: + interpolate_quintic_double_float_ND( + time_before, before, time_after, after, output, time, + /* dimension */ 3, params->periodic, params); + break; + + /* Velocities */ + case field_enum_velocities: + interpolate_cubic_float_ND(time_before, before, time_after, after, output, + time, + /* dimension= */ 3, /* periodic= */ 0, params); + break; + + /* Accelerations */ + case field_enum_accelerations: + interpolate_linear_float_ND( + time_before, before, time_after, after, output, time, + /* dimension= */ 3, /* periodic= */ 0, params); + break; + + /* Interpolation of float */ + case field_enum_masses: + case field_enum_smoothing_lengths: + case field_enum_internal_energies: + case field_enum_entropies: + case field_enum_densities: + interpolate_linear_float(time_before, before, time_after, after, output, + time, /* periodic= */ 0, params); + break; + + /* Constant float term */ + case field_enum_birth_time: + if (*(float *)after->field != *(float *)before->field) { + error("Interpolating different particles."); + } + *(float *)output = *(float *)after->field; + break; + + /* Particle IDs */ + case field_enum_particles_ids: + /* Ensure to have the same particle */ + if (*(long long *)after->field != *(long long *)before->field) { + error("Interpolating different particles."); + } + *(long long *)output = *(long long *)after->field; + break; + + /* Gear star formation */ + case field_enum_gear_star_formation: { + /* Get some variables */ + const float wa = (time - time_before) / (time_after - time_before); + const float wb = 1. - wa; + float *x = (float *)output; + const float *bef = (float *)before->field; + const float *aft = (float *)after->field; + + /* Birth densities + Birth masses */ + x[0] = wa * aft[0] + wb * bef[0]; + x[1] = wa * aft[1] + wb * bef[1]; + + /* ProgenitorIDs */ + long long *id = (long long *)(x + 2); + *id = *(long long *)(aft + 2); + if (*id != *(long long *)(bef + 2)) { + error("Interpolating different particles."); } break; + } - /* Dark matter */ - case swift_type_dark_matter: - case swift_type_dark_matter_background: - case swift_type_neutrino: - if (field->module != field_module_default) - error("Module not implemented"); - gravity_csds_interpolate_field(time_before, before, time_after, after, - output, time, field->local_index, params); + /* GEAR chemistry (part) */ + case field_enum_gear_chemistry_part: + interpolate_linear_double_ND(time_before, before, time_after, after, + output, time, + 2 * GEAR_CHEMISTRY_ELEMENT_COUNT, + /* periodic */ 0, params); break; - /* Stars */ - case swift_type_stars: - switch (field->module) { - case field_module_default: - stars_csds_interpolate_field(time_before, before, time_after, after, - output, time, field->local_index, - params); - break; - - case field_module_chemistry: - chemistry_csds_interpolate_field_spart( - time_before, before, time_after, after, output, time, - field->local_index, params); - break; - - case field_module_star_formation: - star_formation_csds_interpolate_field(time_before, before, time_after, - after, output, time, - field->local_index, params); - break; - - default: - error("Module not implemented"); + /* GEAR chemistry (spart) */ + case field_enum_gear_chemistry_spart: + interpolate_linear_double_ND(time_before, before, time_after, after, + output, time, GEAR_CHEMISTRY_ELEMENT_COUNT, + /* periodic */ 0, params); + break; + + /* Sphenix grouped fields */ + case field_enum_sphenix_secondary_fields: { + /* Get some variables */ + const float wa = (time - time_before) / (time_after - time_before); + const float wb = 1. - wa; + float *x = (float *)output; + const float *bef = (float *)before->field; + const float *aft = (float *)after->field; + + /* Entropy + pressure + Viscosity + Diffusion + Laplacian*/ + const int n_linear = 5; + for (int i = 0; i < n_linear; i++) { + x[i] = wa * aft[i] + wb * bef[i]; } + + /* Div v */ + x[n_linear] = interpolate_cubic_hermite_spline( + time_before, bef[n_linear], bef[n_linear + 1], time_after, + aft[n_linear], aft[n_linear + 1], time); + + /* d Div v / dt */ + x[n_linear + 1] = wa * aft[n_linear + 1] + wb * bef[n_linear + 1]; + /* Use the linear interpolation */ + x[1] = wa * aft[n_linear + 2] + wb * bef[n_linear + 2]; break; + } - /* Default */ default: - error_python("Particle type not implemented"); + error("Field %s not implemented", + field_enum_get_name(field->field.field)); } } diff --git a/src/csds_particle.h b/src/csds_particle.h index 189b2da..eef28eb 100644 --- a/src/csds_particle.h +++ b/src/csds_particle.h @@ -1,5 +1,5 @@ /******************************************************************************* - * This file is part of SWIFT. + * This file is part of CSDS. * Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch) * * This program is free software: you can redistribute it and/or modify @@ -23,13 +23,8 @@ #include "csds_tools.h" /* Include the other local files. */ -#include "csds_chemistry.h" -#include "csds_gravity.h" #include "csds_header.h" -#include "csds_hydro.h" #include "csds_parameters.h" -#include "csds_star_formation.h" -#include "csds_stars.h" #include "csds_time.h" #include <stdio.h> @@ -44,20 +39,17 @@ enum csds_reader_type { }; size_t csds_particle_read_field(const struct csds_reader *reader, size_t offset, - void *output, - const struct field_information *all_fields, - const int all_fields_count, - const int global_index, size_t *mask, + void *output, const enum part_type part_type, + const struct field_information *field, + const int derivative, size_t *mask, size_t *h_offset); -void csds_particle_interpolate_field(const double t_before, - const struct csds_field *restrict before, - const double t_after, - const struct csds_field *restrict after, - void *restrict output, const double t, - const struct field_information *field, - enum part_type type, - const struct csds_parameters *params); +void csds_particle_interpolate_field( + const double t_before, const struct csds_reader_field *restrict before, + const double t_after, const struct csds_reader_field *restrict after, + void *restrict output, const double t, + const struct field_information *field, enum part_type type, + const struct csds_parameters *params); enum csds_special_flags csds_particle_read_special_flag( const struct csds_reader *reader, size_t offset, size_t *mask, diff --git a/src/csds_python_tools.h b/src/csds_python_tools.h index a7da56d..755c196 100644 --- a/src/csds_python_tools.h +++ b/src/csds_python_tools.h @@ -1,5 +1,5 @@ /******************************************************************************* - * This file is part of SWIFT. + * This file is part of CSDS. * Copyright (c) 2020 Loic Hausammann (loic.hausammann@epfl.ch) * * This program is free software: you can redistribute it and/or modify diff --git a/src/csds_python_wrapper.c b/src/csds_python_wrapper.c index 97dbb4b..c0695a0 100644 --- a/src/csds_python_wrapper.c +++ b/src/csds_python_wrapper.c @@ -1,5 +1,5 @@ /******************************************************************************* - * This file is part of SWIFT. + * This file is part of CSDS. * Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch) * * This program is free software: you can redistribute it and/or modify @@ -16,6 +16,7 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. * ******************************************************************************/ +#include "csds_fields.h" #include "csds_header.h" #include "csds_loader_io.h" #include "csds_particle.h" @@ -216,42 +217,6 @@ static PyObject *get_box_size(PyObject *self, PyObject *Py_UNUSED(ignored)) { return (PyObject *)out; } -#define find_field_in_module_internal(MODULE, PART) \ - for (int local = 0; local < MODULE##_csds_field##PART##_count; local++) { \ - const int global = MODULE##_csds_local_to_global##PART[local]; \ - const int local_shifted = local + total_number_fields; \ - if (field_indices[i] == global) { \ - /* Check if we have the same fields for the different modules */ \ - if (current_field != NULL) { \ - if (current_field->dimension != \ - python_fields[local_shifted].dimension || \ - current_field->typenum != python_fields[local_shifted].typenum) { \ - error_python( \ - "The python definition of the field %s does not correspond " \ - "between the modules.", \ - MODULE##_csds_field_names##PART[local]); \ - } \ - } \ - strcpy(field_name, MODULE##_csds_field_names##PART[local]); \ - current_field = &python_fields[local_shifted]; \ - break; \ - } \ - } \ - total_number_fields += MODULE##_csds_field##PART##_count; - -/** - * Same function as find_field_in_module_internal before but with a single - * argument. - */ -#define find_field_in_module_single_particle_type(MODULE) \ - find_field_in_module_internal(MODULE, ) - -/** - * Same function as find_field_in_module_internal but with a cleaner argument. - */ -#define find_field_in_module(MODULE, PART) \ - find_field_in_module_internal(MODULE, _##PART) - /** * @brief Create a list of numpy array containing the fields. * @@ -264,85 +229,50 @@ static PyObject *get_box_size(PyObject *self, PyObject *Py_UNUSED(ignored)) { * @return The python list of numpy array. */ __attribute__((always_inline)) INLINE static PyObject * -csds_loader_create_output(void **output, const int *field_indices, +csds_loader_create_output(void **output, const enum field_enum *fields, const int n_fields, const uint64_t *n_part, uint64_t n_tot) { - struct csds_python_field python_fields[100]; - /* Create the python dictionary */ PyObject *dict = PyDict_New(); - struct csds_python_field *current_field; - - /* Get the hydro fields */ - hydro_csds_generate_python(python_fields); - int total_number_fields = hydro_csds_field_count; - /* Get the gravity fields */ - gravity_csds_generate_python(python_fields + total_number_fields); - total_number_fields += gravity_csds_field_count; - /* Get the stars fields */ - stars_csds_generate_python(python_fields + total_number_fields); - total_number_fields += stars_csds_field_count; - /* Get the chemistry (part) fields */ - chemistry_csds_generate_python_part(python_fields + total_number_fields); - total_number_fields += chemistry_csds_field_part_count; - /* Get the chemistry (spart) fields */ - chemistry_csds_generate_python_spart(python_fields + total_number_fields); - total_number_fields += chemistry_csds_field_spart_count; - /* Get the star formation fields */ - star_formation_csds_generate_python(python_fields + total_number_fields); - total_number_fields += star_formation_csds_field_count; /* Get all the requested fields */ for (int i = 0; i < n_fields; i++) { - /* Reset the variables. */ - current_field = NULL; - total_number_fields = 0; - char field_name[STRING_SIZE]; - - /* Find the fields in the different modules. */ - find_field_in_module_single_particle_type(hydro); - find_field_in_module_single_particle_type(gravity); - find_field_in_module_single_particle_type(stars); - find_field_in_module(chemistry, part); - find_field_in_module(chemistry, spart); - find_field_in_module_single_particle_type(star_formation); - - /* Check if we got a field */ - if (current_field == NULL) { - error_python("Failed to find the required field"); - } + /* Get all the required informations */ + struct csds_python_field current_field = + field_enum_get_python_field(fields[i]); + PyObject *array = NULL; /* Simple case where we only have normal type */ - if (current_field->typenum != CUSTOM_NPY_TYPE) { - if (current_field->dimension > 1) { - npy_intp dims[2] = {n_tot, current_field->dimension}; - array = PyArray_SimpleNewFromData(2, dims, current_field->typenum, + if (current_field.typenum != CUSTOM_NPY_TYPE) { + if (current_field.dimension > 1) { + npy_intp dims[2] = {n_tot, current_field.dimension}; + array = PyArray_SimpleNewFromData(2, dims, current_field.typenum, output[i]); } else { npy_intp dims = n_tot; - array = PyArray_SimpleNewFromData(1, &dims, current_field->typenum, + array = PyArray_SimpleNewFromData(1, &dims, current_field.typenum, output[i]); } } /* Now the complex types */ else { /* Ensure that the field is properly defined */ - if (current_field->subfields_registred != current_field->dimension) { + if (current_field.subfields_registred != current_field.dimension) { error_python( "It seems that you forgot to register a field (found %i and " "expecting %i)", - current_field->subfields_registred, current_field->dimension); + current_field.subfields_registred, current_field.dimension); } /* Creates the dtype */ - PyObject *names = PyList_New(current_field->dimension); - PyObject *formats = PyList_New(current_field->dimension); - PyObject *offsets = PyList_New(current_field->dimension); + PyObject *names = PyList_New(current_field.dimension); + PyObject *formats = PyList_New(current_field.dimension); + PyObject *offsets = PyList_New(current_field.dimension); /* Gather the data into the lists */ - for (int k = 0; k < current_field->dimension; k++) { - struct csds_python_subfield *subfield = ¤t_field->subfields[k]; + for (int k = 0; k < current_field.dimension; k++) { + struct csds_python_subfield *subfield = ¤t_field.subfields[k]; /* Transform the information into python */ PyObject *name = PyUnicode_FromString(subfield->name); @@ -380,7 +310,8 @@ csds_loader_create_output(void **output, const int *field_indices, NPY_ARRAY_CARRAY, NULL); } - csds_loader_python_field_free(current_field); + csds_loader_python_field_free(¤t_field); + const char *field_name = field_enum_get_name(fields[i]); PyDict_SetItemString(dict, field_name, array); } @@ -487,69 +418,54 @@ static PyObject *pyGetListFields(__attribute__((unused)) PyObject *self, const struct header *h = &reader->log.header; /* Create the array to check if a field is present. */ - int *field_present = (int *)malloc(h->masks_count * sizeof(int)); + int *field_present = (int *)malloc(field_enum_count * sizeof(int)); if (field_present == NULL) { error("Failed to allocate the memory for the fields present."); } - /* Initialize the array */ - for (int i = 0; i < h->masks_count; i++) { + /* Find if the fields are present in all the types. */ + for (enum field_enum i = 0; i < field_enum_count; i++) { field_present[i] = 1; - } - - /* Check all the fields */ - for (int i = 0; i < swift_type_count; i++) { - /* Skip the types that are not required */ - if (read_types[i] == 0) continue; - /* Get the number of fields for the current particle type. */ - int number_fields = tools_get_number_fields(i); - if (number_fields == 0) { - continue; - } - - /* Get the list of fields for the current particle type. */ - struct field_information *fields = (struct field_information *)malloc( - number_fields * sizeof(struct field_information)); - if (fields == NULL) error("Failed to initialize the field information"); - tools_get_list_fields(fields, i, h); + /* Check all the types */ + for (int type = 0; type < swift_type_count; type++) { + /* Skip the types that are not required + and the field not found */ + if (read_types[type] == 0 || field_present[i] == 0) continue; - for (int j = 0; j < h->masks_count; j++) { - /* Skip the fields not found in previous type. */ - if (field_present[j] == 0) continue; - - /* Check if the field is present */ + /* Search among all the fields of the particle type */ int found = 0; - for (int k = 0; k < number_fields; k++) { - if (strcmp(h->masks[j].name, fields[k].name) == 0) { + for (int j = 0; j < h->number_fields[type]; j++) { + if (i == h->fields[type][j].field.field) { found = 1; break; } } /* Set the field as not found */ - if (!found) field_present[j] = 0; + if (!found) field_present[i] = 0; } - - /* Free the memory */ - free(fields); } + /* Remove the special flags */ + field_present[field_enum_special_flags] = 0; + /* Count the number of fields found */ int number_fields = 0; - for (int i = 0; i < h->masks_count; i++) { + for (int i = 0; i < field_enum_count; i++) { number_fields += field_present[i]; } /* Create the python list for the output*/ PyObject *list = PyList_New(number_fields); int current = 0; - for (int i = 0; i < h->masks_count; i++) { + for (int i = 0; i < field_enum_count; i++) { /* Keep only the field present. */ if (field_present[i] == 0) continue; - PyObject *name = PyUnicode_FromString(h->masks[i].name); - PyList_SetItem(list, current, name); + const char *name = field_enum_get_name(i); + PyObject *python_name = PyUnicode_FromString(name); + PyList_SetItem(list, current, python_name); current += 1; } @@ -720,14 +636,12 @@ static PyObject *pyGetData(__attribute__((unused)) PyObject *self, /* initialize the reader. */ struct csds_reader *reader = &self_reader->reader; - const struct header *h = &reader->log.header; - /* Get the fields indexes from the header. */ + /* Get the field enum from the header. */ const int n_fields = PyList_Size(fields); - int *field_indices = (int *)malloc(n_fields * sizeof(int)); + enum field_enum *fields_wanted = + (enum field_enum *)malloc(n_fields * sizeof(enum field_enum)); for (int i = 0; i < n_fields; i++) { - field_indices[i] = -1; - /* Get the an item in the list. */ PyObject *field = PyList_GetItem(fields, i); if (!PyUnicode_Check(field)) { @@ -738,18 +652,8 @@ static PyObject *pyGetData(__attribute__((unused)) PyObject *self, Py_ssize_t size = 0; const char *field_name = PyUnicode_AsUTF8AndSize(field, &size); - /* Find the equivalent field inside the header. */ - for (int j = 0; j < h->masks_count; j++) { - if (strcmp(field_name, h->masks[j].name) == 0) { - field_indices[i] = j; - break; - } - } - - /* Check if we found a field. */ - if (field_indices[i] == -1) { - error_python("Failed to find the field %s", field_name); - } + /* Set the enum */ + fields_wanted[i] = field_enum_get(field_name); } /* Set the time. */ @@ -775,10 +679,16 @@ static PyObject *pyGetData(__attribute__((unused)) PyObject *self, n_tot += n_part[i]; } + /* Print the number of particles */ + if (reader->verbose > 0) { + message("Found %lu particles", n_tot); + } + /* Allocate the output memory. */ void **output = malloc(n_fields * sizeof(void *)); for (int i = 0; i < n_fields; i++) { - output[i] = malloc(n_tot * h->masks[field_indices[i]].size); + unsigned int size = field_enum_get_size(fields_wanted[i]); + output[i] = malloc(n_tot * size); } /* Read the particles. */ @@ -786,7 +696,7 @@ static PyObject *pyGetData(__attribute__((unused)) PyObject *self, /* Enable multithreading */ Py_BEGIN_ALLOW_THREADS; - csds_reader_read_all_particles(reader, time, csds_reader_lin, field_indices, + csds_reader_read_all_particles(reader, time, csds_reader_lin, fields_wanted, n_fields, output, n_part); /* Disable multithreading */ @@ -816,7 +726,7 @@ static PyObject *pyGetData(__attribute__((unused)) PyObject *self, /* Read the particles. */ csds_reader_read_particles_from_ids(reader, time, csds_reader_lin, - field_indices, n_fields, output, n_part, + fields_wanted, n_fields, output, n_part, c_part_ids); /* Disable multithreading */ @@ -833,10 +743,10 @@ static PyObject *pyGetData(__attribute__((unused)) PyObject *self, /* Create the python output. */ PyObject *array = - csds_loader_create_output(output, field_indices, n_fields, n_part, n_tot); + csds_loader_create_output(output, fields_wanted, n_fields, n_part, n_tot); /* Free the reader. */ - free(field_indices); + free(fields_wanted); return array; } diff --git a/src/csds_reader.c b/src/csds_reader.c index 395bf11..86b4142 100644 --- a/src/csds_reader.c +++ b/src/csds_reader.c @@ -1,5 +1,5 @@ /******************************************************************************* - * This file is part of SWIFT. + * This file is part of CSDS. * Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch) * * This program is free software: you can redistribute it and/or modify @@ -20,6 +20,8 @@ /* Include corresponding header */ #include "csds_reader.h" +/* Include CSDS headers */ +#include "csds_fields.h" #include "csds_parameters.h" /* Include standard library */ @@ -373,8 +375,6 @@ void csds_reader_get_number_particles(struct csds_reader *reader, * @param offset_last_full_record The offset of this particle last record containing all the fields. * @param field_wanted The field to read. - * @param all_fields The list of all the fields in the particle type. - * @param all_fields_count The number of elements in all_fields. * @param output Pointers to the output array * @param type The #part_type. * @@ -385,34 +385,16 @@ int csds_reader_read_field(struct csds_reader *reader, double time, enum csds_reader_type interp_type, const size_t offset_last_full_record, const struct field_information *field_wanted, - const struct field_information *all_fields, - const int all_fields_count, void *output, - enum part_type type) { + void *output, enum part_type type) { const struct header *h = &reader->log.header; size_t offset = offset_last_full_record; - /* Get the indexes */ - const int global = field_wanted->global_index; - const int global_first = field_wanted->global_index_first; - const int global_second = field_wanted->global_index_second; - /* Check if the offsets are correct. */ if (offset > offset_time) { error_python("Last offset is larger than the requested time."); } - /* Get the masks. */ - struct mask_data *mask_field = &h->masks[global]; - struct mask_data *mask_first = NULL; - if (global_first >= 0) { - mask_first = &h->masks[global_first]; - } - struct mask_data *mask_second = NULL; - if (global_second >= 0) { - mask_second = &h->masks[global_second]; - } - /* Offset of the field read before the requested time. */ size_t offset_before = offset_last_full_record; @@ -428,7 +410,7 @@ int csds_reader_read_field(struct csds_reader *reader, double time, csds_loader_io_read_mask(h, reader->log.log.map + offset, &mask, &h_offset); /* Is the particle removed from the logfile? */ - if (mask & h->masks[csds_index_special_flags].mask) { + if (mask & SPECIAL_FLAGS_MASK) { int data = 0; int part_type = 0; enum csds_special_flags flag = csds_particle_read_special_flag( @@ -452,7 +434,7 @@ int csds_reader_read_field(struct csds_reader *reader, double time, } /* Check if the field is present */ - if (mask & mask_field->mask) { + if (mask & field_wanted->field.mask) { offset_before = offset; } @@ -461,29 +443,32 @@ int csds_reader_read_field(struct csds_reader *reader, double time, } /* Read the field */ - csds_particle_read_field(reader, offset_before, output, all_fields, - all_fields_count, global, &mask, &h_offset); + csds_particle_read_field(reader, offset_before, output, type, field_wanted, + /* derivative */ 0, &mask, &h_offset); /* Deal with the first derivative. */ - const size_t size_first = mask_first == NULL ? 0 : mask_first->size; - char first_deriv[size_first]; + int first_found = field_wanted->first.field != field_enum_none && + mask & field_wanted->first.mask; + char first_deriv[field_wanted->first.size]; - int first_found = mask_first != NULL && mask & mask_first->mask; if (first_found) { /* Read the first derivative */ - csds_particle_read_field(reader, offset_before, first_deriv, all_fields, - all_fields_count, global_first, &mask, &h_offset); + csds_particle_read_field(reader, offset_before, first_deriv, type, + field_wanted, /* derivative */ 1, &mask, + &h_offset); } /* Deal with the second derivative. */ - const size_t size_second = mask_second == NULL ? 0 : mask_second->size; - char second_deriv[size_second]; - int second_found = mask_second != NULL && mask & mask_second->mask; + int second_found = field_wanted->second.field != field_enum_none && + mask & field_wanted->second.mask; + char second_deriv[field_wanted->second.size]; + if (second_found) { /* Read the first derivative */ - csds_particle_read_field(reader, offset_before, second_deriv, all_fields, - all_fields_count, global_second, &mask, &h_offset); + csds_particle_read_field(reader, offset_before, second_deriv, type, + field_wanted, /* derivative */ 2, &mask, + &h_offset); } /* Get the time. */ @@ -503,7 +488,7 @@ int csds_reader_read_field(struct csds_reader *reader, double time, &h_offset); /* Do we have the field? */ - if (mask & mask_field->mask) { + if (mask & field_wanted->field.mask) { break; } @@ -517,34 +502,35 @@ int csds_reader_read_field(struct csds_reader *reader, double time, } /* Output after the requested time. */ - char output_after[mask_field->size]; + char output_after[field_wanted->field.size]; /* Read the field */ - csds_particle_read_field(reader, offset, output_after, all_fields, - all_fields_count, global, &mask, &h_offset); + csds_particle_read_field(reader, offset, output_after, type, field_wanted, + /* derivative */ 0, &mask, &h_offset); /* Deal with the first derivative. */ - char first_deriv_after[size_first]; + char first_deriv_after[field_wanted->first.size]; /* Did we find the derivative before and in this record? */ - first_found = mask_first != NULL && first_found && mask & mask_first->mask; + first_found = field_wanted->first.field != field_enum_none && first_found && + mask & field_wanted->first.mask; if (first_found) { /* Read the first derivative */ - csds_particle_read_field(reader, offset, first_deriv_after, all_fields, - all_fields_count, global_first, &mask, + csds_particle_read_field(reader, offset, first_deriv_after, type, + field_wanted, /*derivative*/ 1, &mask, &h_offset); } /* Deal with the second derivative. */ - char second_deriv_after[size_second]; + char second_deriv_after[field_wanted->second.size]; /* Did we find the derivative before and in this record? */ - second_found = - mask_second != NULL && second_found && mask & mask_second->mask; + second_found = field_wanted->second.field != field_enum_none && + second_found && mask & field_wanted->second.mask; if (second_found) { /* Read the second derivative */ - csds_particle_read_field(reader, offset, second_deriv_after, all_fields, - all_fields_count, global_second, &mask, + csds_particle_read_field(reader, offset, second_deriv_after, type, + field_wanted, /* derivative */ 2, &mask, &h_offset); } @@ -553,8 +539,8 @@ int csds_reader_read_field(struct csds_reader *reader, double time, double time_after = time_array_get_time(&reader->log.times, offset); /* Deal with the derivatives */ - struct csds_field before; - struct csds_field after; + struct csds_reader_field before; + struct csds_reader_field after; before.field = output; before.first_deriv = first_found ? first_deriv : NULL; before.second_deriv = second_found ? second_deriv : NULL; @@ -574,36 +560,32 @@ int csds_reader_read_field(struct csds_reader *reader, double time, * @brief Convert the fields from global indexes to local. * * @param reader The #csds_reader. - * @param global_fields_wanted The fields to convert. - * @param fields_wanted (out) Fields wanted (need to be allocated). - * @param n_fields_wanted Number of elements in global_fields_wanted, + * @param fields_wanted_wanted The fields to convert. + * @param fields_info_wanted (out) #field_information for the wanted fields + * (need to be allocated). + * @param n_fields_wanted Number of elements in fields_wanted, * local_fields_wanted and its derivatives. - * @param all_fields The list of all the fields. - * @param all_fields_count The number of elements in all_fields. * @param type The type of the particle. */ -void csds_reader_get_fields_wanted(const struct csds_reader *reader, - const int *global_fields_wanted, - struct field_information *fields_wanted, - const int n_fields_wanted, - struct field_information *all_fields, - int all_fields_count, enum part_type type) { +void csds_reader_get_fields_wanted( + const struct csds_reader *reader, const enum field_enum *fields_wanted, + const struct field_information **fields_found, const int n_fields_wanted, + enum part_type type) { const struct header *h = &reader->log.header; - /* Mark fields_wanted in order to check that the fields are found. */ + /* Mark fields_found in order to check that the fields are found. */ for (int i = 0; i < n_fields_wanted; i++) { - fields_wanted[i].local_index = -1; + fields_found[i] = NULL; } /* Find the corresponding fields */ for (int i = 0; i < n_fields_wanted; i++) { - const int global = global_fields_wanted[i]; - for (int j = 0; j < all_fields_count; j++) { + for (int j = 0; j < h->number_fields[type]; j++) { /* Copy the structure if the field is found. */ - if (all_fields[j].global_index == global) { - fields_wanted[i] = all_fields[j]; + if (h->fields[type][j].field.field == fields_wanted[i]) { + fields_found[i] = &h->fields[type][j]; break; } } @@ -611,9 +593,8 @@ void csds_reader_get_fields_wanted(const struct csds_reader *reader, /* Check that we found the fields */ for (int i = 0; i < n_fields_wanted; i++) { - if (fields_wanted[i].local_index < 0) { - const int global = global_fields_wanted[i]; - const char *name = h->masks[global].name; + if (fields_found[i] == NULL) { + const char *name = field_enum_get_name(fields_wanted[i]); error_python("Field %s not found in particle type %s", name, part_type_names[type]); } @@ -636,7 +617,7 @@ struct extra_data_read_all { enum csds_reader_type interp_type; /*! The #field_information for all the fields requested. */ - const struct field_information *fields_wanted; + const struct field_information **fields_wanted; /*! The number of fields requested. */ int n_fields_wanted; @@ -650,12 +631,6 @@ struct extra_data_read_all { /*! The type of particle to read. */ enum part_type type; - /*! The number of fields available for this type of particles. */ - int all_fields_count; - - /*! The #field_information of all the fields available. */ - const struct field_information *all_fields; - /*! The current position in the index file. */ size_t current_index; @@ -681,17 +656,13 @@ void csds_reader_read_all_particles_single_type_mapper(void *map_data, struct csds_reader *reader = data_tp->reader; double time = data_tp->time; enum csds_reader_type interp_type = data_tp->interp_type; - const struct field_information *fields_wanted = data_tp->fields_wanted; + const struct field_information **fields_wanted = data_tp->fields_wanted; int n_fields_wanted = data_tp->n_fields_wanted; void **output = data_tp->output; const uint64_t *n_part = data_tp->n_part; enum part_type type = data_tp->type; - int all_fields_count = data_tp->all_fields_count; - const struct field_information *all_fields = data_tp->all_fields; /* Create a few variables for later */ - const struct header *h = &reader->log.header; - const size_t size_index = reader->index.index_prev.nparts[type]; const size_t size_history = csds_reader_count_number_new_particles(reader, type); @@ -713,8 +684,7 @@ void csds_reader_read_all_particles_single_type_mapper(void *map_data, error_python("Failed to allocate the temporary output buffer"); } for (int i = 0; i < n_fields_wanted; i++) { - const int global = fields_wanted[i].global_index; - output_tmp[i] = malloc(h->masks[global].size); + output_tmp[i] = malloc(fields_wanted[i]->field.size); if (output_tmp[i] == NULL) { error_python("Failed to allocate the temporary output buffer"); } @@ -747,8 +717,7 @@ void csds_reader_read_all_particles_single_type_mapper(void *map_data, /* Read the field. */ particle_removed = csds_reader_read_field( reader, time, reader->time.time_offset, interp_type, offset, - &fields_wanted[field], all_fields, all_fields_count, - output_tmp[field], type); + fields_wanted[field], output_tmp[field], type); /* Should we continue to read the fields of this particle? */ if (particle_removed) { @@ -763,13 +732,12 @@ void csds_reader_read_all_particles_single_type_mapper(void *map_data, for (int field = 0; field < n_fields_wanted; field++) { /* Get the array */ - const int global = fields_wanted[field].global_index; + const size_t size = fields_wanted[field]->field.size; void *output_single = - (char *)output[field] + - (current_output + prev_npart) * h->masks[global].size; + (char *)output[field] + (current_output + prev_npart) * size; /* Copy the temporary buffer into the global one */ - memcpy(output_single, output_tmp[field], h->masks[global].size); + memcpy(output_single, output_tmp[field], size); } } } @@ -787,7 +755,7 @@ void csds_reader_read_all_particles_single_type_mapper(void *map_data, * @param reader The #csds_reader. * @param time The requested time for the particle. * @param interp_type The type of interpolation. - * @param global_fields_wanted The fields requested (global index). + * @param fields_wanted The fields requested. * @param n_fields_wanted Number of field requested. * @param output Pointer to the output array. Size: (n_fields_wanted, * sum(n_part)). @@ -796,39 +764,30 @@ void csds_reader_read_all_particles_single_type_mapper(void *map_data, */ void csds_reader_read_all_particles_single_type( struct csds_reader *reader, double time, enum csds_reader_type interp_type, - const int *global_fields_wanted, const int n_fields_wanted, void **output, - const uint64_t *n_part, enum part_type type) { - const struct header *h = &reader->log.header; + const enum field_enum *fields_wanted, const int n_fields_wanted, + void **output, const uint64_t *n_part, enum part_type type) { /* Allocate temporary memory. */ - struct field_information *fields_wanted = (struct field_information *)malloc( - sizeof(struct field_information) * n_fields_wanted); - if (fields_wanted == NULL) { + const struct field_information **fields_info_wanted = + (const struct field_information **)malloc( + sizeof(struct field_information *) * n_fields_wanted); + if (fields_info_wanted == NULL) { error_python("Failed to allocate the field information."); } - /* Get the list of fields. */ - const int all_fields_count = tools_get_number_fields(type); - struct field_information *all_fields = (struct field_information *)malloc( - all_fields_count * sizeof(struct field_information)); - tools_get_list_fields(all_fields, type, h); - /* Convert fields into the local array. */ - csds_reader_get_fields_wanted(reader, global_fields_wanted, fields_wanted, - n_fields_wanted, all_fields, all_fields_count, - type); + csds_reader_get_fields_wanted(reader, fields_wanted, fields_info_wanted, + n_fields_wanted, type); /* Compact the data */ struct extra_data_read_all extra = {reader, time, interp_type, - fields_wanted, + fields_info_wanted, n_fields_wanted, output, n_part, type, - all_fields_count, - all_fields, /* current_index */ 0, /* current_output */ 0}; @@ -847,8 +806,7 @@ void csds_reader_read_all_particles_single_type( /* Free the memory */ threadpool_clean(&tp); - free(all_fields); - free(fields_wanted); + free(fields_info_wanted); } /** @@ -857,7 +815,7 @@ void csds_reader_read_all_particles_single_type( * @param reader The #csds_reader. * @param time The requested time for the particle. * @param interp_type The type of interpolation. - * @param global_fields_wanted The fields requested (global index). + * @param fields_wanted The fields requested. * @param n_fields_wanted Number of field requested. * @param output Pointer to the output array. Size: (n_fields_wanted, * sum(n_part)). @@ -865,40 +823,40 @@ void csds_reader_read_all_particles_single_type( */ void csds_reader_read_all_particles(struct csds_reader *reader, double time, enum csds_reader_type interp_type, - const int *global_fields_wanted, + const enum field_enum *fields_wanted, const int n_fields_wanted, void **output, const uint64_t *n_part) { /* Read the gas */ if (n_part[swift_type_gas] != 0) { - csds_reader_read_all_particles_single_type( - reader, time, interp_type, global_fields_wanted, n_fields_wanted, - output, n_part, swift_type_gas); + csds_reader_read_all_particles_single_type(reader, time, interp_type, + fields_wanted, n_fields_wanted, + output, n_part, swift_type_gas); } /* Read the dark matter. */ if (n_part[swift_type_dark_matter] != 0) { csds_reader_read_all_particles_single_type( - reader, time, interp_type, global_fields_wanted, n_fields_wanted, - output, n_part, swift_type_dark_matter); + reader, time, interp_type, fields_wanted, n_fields_wanted, output, + n_part, swift_type_dark_matter); } /* Read the dark matter background. */ if (n_part[swift_type_dark_matter_background] != 0) { csds_reader_read_all_particles_single_type( - reader, time, interp_type, global_fields_wanted, n_fields_wanted, - output, n_part, swift_type_dark_matter_background); + reader, time, interp_type, fields_wanted, n_fields_wanted, output, + n_part, swift_type_dark_matter_background); } /* Read the neutrino dark matter. */ if (n_part[swift_type_neutrino] != 0) { csds_reader_read_all_particles_single_type( - reader, time, interp_type, global_fields_wanted, n_fields_wanted, - output, n_part, swift_type_neutrino); + reader, time, interp_type, fields_wanted, n_fields_wanted, output, + n_part, swift_type_neutrino); } /* Read the stars. */ if (n_part[swift_type_stars] != 0) { csds_reader_read_all_particles_single_type( - reader, time, interp_type, global_fields_wanted, n_fields_wanted, - output, n_part, swift_type_stars); + reader, time, interp_type, fields_wanted, n_fields_wanted, output, + n_part, swift_type_stars); } } @@ -908,7 +866,7 @@ void csds_reader_read_all_particles(struct csds_reader *reader, double time, * @param reader The #csds_reader. * @param time The requested time for the particle. * @param interp_type The type of interpolation. - * @param global_fields_wanted The fields requested (global index). + * @param fields_wanted The fields requested. * @param n_fields_wanted Number of field requested. * @param output Pointer to the output array. Size: (n_fields_wanted, * sum(n_part)). @@ -920,10 +878,8 @@ void csds_reader_read_all_particles(struct csds_reader *reader, double time, */ void csds_reader_read_particles_from_ids_single_type( struct csds_reader *reader, double time, enum csds_reader_type interp_type, - const int *global_fields_wanted, const int n_fields_wanted, void **output, - uint64_t *n_part, enum part_type type, long long *ids) { - - const struct header *h = &reader->log.header; + const enum field_enum *fields_wanted, const int n_fields_wanted, + void **output, uint64_t *n_part, enum part_type type, long long *ids) { /* Count the number of previous parts for the shift in output */ uint64_t prev_npart = 0; @@ -932,22 +888,16 @@ void csds_reader_read_particles_from_ids_single_type( } /* Allocate temporary memory. */ - struct field_information *fields_wanted = (struct field_information *)malloc( - sizeof(struct field_information) * n_fields_wanted); - if (fields_wanted == NULL) { + const struct field_information **fields_info_wanted = + (const struct field_information **)malloc( + sizeof(struct field_information *) * n_fields_wanted); + if (fields_info_wanted == NULL) { error_python("Failed to allocate the field information."); } - /* Get the list of fields. */ - const int all_fields_count = tools_get_number_fields(type); - struct field_information *all_fields = (struct field_information *)malloc( - all_fields_count * sizeof(struct field_information)); - tools_get_list_fields(all_fields, type, h); - /* Convert fields into the local array. */ - csds_reader_get_fields_wanted(reader, global_fields_wanted, fields_wanted, - n_fields_wanted, all_fields, all_fields_count, - type); + csds_reader_get_fields_wanted(reader, fields_wanted, fields_info_wanted, + n_fields_wanted, type); /* Read the particles */ for (size_t i = 0; i < n_part[type]; i++) { @@ -967,15 +917,14 @@ void csds_reader_read_particles_from_ids_single_type( /* Loop over each field. */ for (int field = 0; field < n_fields_wanted; field++) { - const int global = fields_wanted[field].global_index; void *output_single = - output[field] + (i + prev_npart) * h->masks[global].size; + output[field] + + (i + prev_npart) * fields_info_wanted[field]->field.size; /* Read the field. */ int particle_removed = csds_reader_read_field( reader, time, reader->time.time_offset, interp_type, offset, - &fields_wanted[field], all_fields, all_fields_count, output_single, - type); + fields_info_wanted[field], output_single, type); /* Is the particle still present? */ if (particle_removed) { @@ -988,8 +937,7 @@ void csds_reader_read_particles_from_ids_single_type( } /* Free the memory */ - free(all_fields); - free(fields_wanted); + free(fields_info_wanted); } /** @@ -998,7 +946,7 @@ void csds_reader_read_particles_from_ids_single_type( * @param reader The #csds_reader. * @param time The requested time for the particle. * @param interp_type The type of interpolation. - * @param global_fields_wanted The fields requested (global index). + * @param fields_wanted The fields requested. * @param n_fields_wanted Number of field requested. * @param output Pointer to the output array. Size: (n_fields_wanted, * sum(n_part)). @@ -1009,34 +957,34 @@ void csds_reader_read_particles_from_ids_single_type( */ void csds_reader_read_particles_from_ids( struct csds_reader *reader, double time, enum csds_reader_type interp_type, - const int *global_fields_wanted, const int n_fields_wanted, void **output, - uint64_t *n_part, long long **ids) { + const enum field_enum *fields_wanted, const int n_fields_wanted, + void **output, uint64_t *n_part, long long **ids) { /* Read the gas */ if (n_part[swift_type_gas] != 0) { csds_reader_read_particles_from_ids_single_type( - reader, time, interp_type, global_fields_wanted, n_fields_wanted, - output, n_part, swift_type_gas, ids[swift_type_gas]); + reader, time, interp_type, fields_wanted, n_fields_wanted, output, + n_part, swift_type_gas, ids[swift_type_gas]); } /* Read the dark matter. */ if (n_part[swift_type_dark_matter] != 0) { csds_reader_read_particles_from_ids_single_type( - reader, time, interp_type, global_fields_wanted, n_fields_wanted, - output, n_part, swift_type_dark_matter, ids[swift_type_dark_matter]); + reader, time, interp_type, fields_wanted, n_fields_wanted, output, + n_part, swift_type_dark_matter, ids[swift_type_dark_matter]); } /* Read the dark matter background. */ if (n_part[swift_type_dark_matter_background] != 0) { csds_reader_read_particles_from_ids_single_type( - reader, time, interp_type, global_fields_wanted, n_fields_wanted, - output, n_part, swift_type_dark_matter_background, + reader, time, interp_type, fields_wanted, n_fields_wanted, output, + n_part, swift_type_dark_matter_background, ids[swift_type_dark_matter_background]); } /* Read the stars. */ if (n_part[swift_type_stars] != 0) { csds_reader_read_particles_from_ids_single_type( - reader, time, interp_type, global_fields_wanted, n_fields_wanted, - output, n_part, swift_type_stars, ids[swift_type_stars]); + reader, time, interp_type, fields_wanted, n_fields_wanted, output, + n_part, swift_type_stars, ids[swift_type_stars]); } } @@ -1109,24 +1057,15 @@ size_t csds_reader_read_record(struct csds_reader *reader, void **output, /* Read the record's mask. */ map = csds_loader_io_read_mask(h, map + offset, &mask, &h_offset); - *is_particle = !(mask & h->timestamp_mask); + *is_particle = !(mask & TIMESTAMP_MASK); /* The record is a particle. */ if (*is_particle) { - /* Get the list of fields. */ - const int all_fields_count = tools_get_number_fields(swift_type_gas); - struct field_information *all_fields = (struct field_information *)malloc( - all_fields_count * sizeof(struct field_information)); - tools_get_list_fields(all_fields, swift_type_gas, h); - size_t offset_tmp = offset; - for (int i = 0; i < all_fields_count; i++) { + for (int i = 0; i < h->number_fields[swift_type_gas]; i++) { offset = csds_particle_read_field( - reader, offset_tmp, output[i], all_fields, all_fields_count, - all_fields[i].global_index, &mask, &h_offset); + reader, offset_tmp, output[i], swift_type_gas, + &h->fields[swift_type_gas][i], /* derivative */ 0, &mask, &h_offset); } - - /* Cleanup */ - free(all_fields); } /* The record is a timestamp. */ else { diff --git a/src/csds_reader.h b/src/csds_reader.h index 6f5fbe4..224a72f 100644 --- a/src/csds_reader.h +++ b/src/csds_reader.h @@ -1,5 +1,5 @@ /******************************************************************************* - * This file is part of SWIFT. + * This file is part of CSDS. * Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch) * * This program is free software: you can redistribute it and/or modify @@ -131,23 +131,21 @@ void csds_reader_get_number_particles(struct csds_reader *reader, void csds_reader_read_all_particles(struct csds_reader *reader, double time, enum csds_reader_type interp_type, - const int *id_masks_wanted, - const int n_mask_wanted, void **output, + const enum field_enum *fields_wanted, + const int n_fields_wanted, void **output, const uint64_t *n_part); void csds_reader_read_particles_from_ids(struct csds_reader *reader, double time, enum csds_reader_type interp_type, - const int *id_masks_wanted, + const enum field_enum *id_masks_wanted, const int n_mask_wanted, void **output, uint64_t *n_part, long long **ids); size_t csds_reader_read_record(struct csds_reader *reader, void **output, double *time, int *is_particle, size_t offset); void csds_reader_generate_index_files(const struct csds_reader *reader, int number_index); -void csds_reader_get_fields_wanted(const struct csds_reader *reader, - const int *global_fields_wanted, - struct field_information *fields_wanted, - const int n_fields_wanted, - struct field_information *all_fields, - int all_fields_count, enum part_type type); +void csds_reader_get_fields_wanted( + const struct csds_reader *reader, const enum field_enum *fields_wanted, + const struct field_information **fields_found, const int n_fields_wanted, + enum part_type type); #endif // CSDS_CSDS_READER_H diff --git a/src/csds_reader_generate_index.c b/src/csds_reader_generate_index.c index b17c9aa..b2bb4dd 100644 --- a/src/csds_reader_generate_index.c +++ b/src/csds_reader_generate_index.c @@ -1,5 +1,5 @@ /******************************************************************************* - * This file is part of SWIFT. + * This file is part of CSDS. * Copyright (c) 2021 Loic Hausammann (loic.hausammann@epfl.ch) * * This program is free software: you can redistribute it and/or modify @@ -376,10 +376,6 @@ void csds_reader_write_index(const struct csds_reader *reader, * @param reader The #csds_reader. * @param current_state The arrays that will be updated with (size * swift_type_count). - * @param all_fields The complete list of fields for the different particles - * (size swift_type_count). - * @param all_fields_count The number of fields for each particle type (size - * swift_type_count). * @param time_record (output) The first #time_record in the logfile (for * writing the index file). * @@ -388,21 +384,13 @@ void csds_reader_write_index(const struct csds_reader *reader, */ size_t csds_reader_get_initial_state(const struct csds_reader *reader, struct index_writer *current_state, - struct field_information **all_fields, - int *all_fields_count, struct time_record *time_record) { /* Get a few variables. */ const struct csds_logfile *log = &reader->log; - const struct header *header = &log->header; + const struct header *h = &log->header; const int size_record_header = CSDS_MASK_SIZE + CSDS_OFFSET_SIZE; - /* Get the mask for the IDs */ - const int index_mask_id = header_get_field_index(header, "ParticleIDs"); - if (index_mask_id == -1) { - error("The logfile does not contain the field ParticleIDs"); - } - /* Get the offset after the dump of all the particles and the time information of the first time record. */ if (log->times.size < 2) { @@ -412,13 +400,12 @@ size_t csds_reader_get_initial_state(const struct csds_reader *reader, *time_record = log->times.records[0]; /* Get the offset of the first particle record */ - size_t offset_first = header->offset_first_record; + size_t offset_first = h->offset_first_record; /* Skip the time record */ size_t time_mask = 0; - csds_loader_io_read_mask(header, log->log.map + offset_first, &time_mask, - NULL); - const int time_size = header_get_record_size_from_mask(header, time_mask); + csds_loader_io_read_mask(h, log->log.map + offset_first, &time_mask, NULL); + const int time_size = header_get_record_size_from_mask(h, time_mask); offset_first += time_size + size_record_header; /* Get the initial state */ @@ -441,22 +428,25 @@ size_t csds_reader_get_initial_state(const struct csds_reader *reader, error("Particle type not implemented (%s)", part_type_names[part_type]); } + /* Get the mask for the IDs */ + const struct field_information *field_id = + header_get_field_from_name(h, "ParticleIDs", part_type); + /* Get the particle ID */ - if (!(mask & header->masks[index_mask_id].mask)) { + if (!(mask & field_id->field.mask)) { error("The particle ID is missing in the first log"); } /* Read the particle ID */ int64_t id = 0; - csds_particle_read_field(reader, offset, &id, all_fields[part_type], - all_fields_count[part_type], index_mask_id, &mask, - &prev_offset); + csds_particle_read_field(reader, offset, &id, part_type, field_id, + /* derivative */ 0, &mask, &prev_offset); /* Log the particle */ index_writer_log(¤t_state[part_type], id, offset); /* Increment the offset */ - const int record_size = header_get_record_size_from_mask(header, mask); + const int record_size = header_get_record_size_from_mask(h, mask); offset += record_size + size_record_header; } @@ -470,10 +460,6 @@ size_t csds_reader_get_initial_state(const struct csds_reader *reader, * @param reader The #csds_reader. * @param init_offset The initial offset to read. * @param time_record The #time_record of the next index file. - * @param all_fields The complete list of fields for the different particles - * (size swift_type_count). - * @param all_fields_count The number of fields for each particle type (size - * swift_type_count). * @param current_state The arrays containing the state of the simulation (size * swift_type_count). * @param parts_created The arrays containing the particles created since last @@ -485,38 +471,30 @@ size_t csds_reader_get_initial_state(const struct csds_reader *reader, */ size_t csds_reader_update_state_to_next_index( const struct csds_reader *reader, size_t init_offset, - struct time_record time_record, struct field_information **all_fields, - int *all_fields_count, struct index_writer *current_state, + struct time_record time_record, struct index_writer *current_state, struct index_writer *parts_created, struct index_writer *parts_removed) { const struct csds_logfile *log = &reader->log; - const struct header *header = &log->header; + const struct header *h = &log->header; const int size_record_header = CSDS_MASK_SIZE + CSDS_OFFSET_SIZE; - /* Get the mask for the IDs */ - const int index_mask_id = header_get_field_index(header, "ParticleIDs"); - if (index_mask_id == -1) { - error("The logfile does not contain the field ParticleIDs"); - } - /* Look for all the created / removed particles */ size_t offset = init_offset; while (offset < time_record.offset) { size_t mask = 0; size_t h_offset = 0; - int part_type = 0; // only available if the record is flagged. + int part_type = -1; // only available if the record is flagged. int data = 0; /* Get the mask */ - csds_loader_io_read_mask(header, log->log.map + offset, &mask, &h_offset); + csds_loader_io_read_mask(h, log->log.map + offset, &mask, &h_offset); /* Go to the next record */ const size_t old_offset = offset; - offset += header_get_record_size_from_mask(header, mask); + offset += header_get_record_size_from_mask(h, mask); offset += size_record_header; /* Check if we have a particle with a flag */ - if (mask & header->masks[csds_index_timestamp].mask || - !(mask & header->masks[csds_index_special_flags].mask)) { + if (mask & TIMESTAMP_MASK || !(mask & SPECIAL_FLAGS_MASK)) { continue; } @@ -531,12 +509,19 @@ size_t csds_reader_update_state_to_next_index( "for a flag and the flag set to 0"); } #endif + /* Check if we have a meaningful particle type */ + if (part_type < 0) { + error("Found a special flag without a particle type"); + } + + /* Get the mask for the IDs */ + const struct field_information *field_id = + header_get_field_from_name(h, "ParticleIDs", part_type); /* Read the ID */ int64_t id = 0; - csds_particle_read_field(reader, old_offset, &id, all_fields[part_type], - all_fields_count[part_type], index_mask_id, &mask, - &h_offset); + csds_particle_read_field(reader, old_offset, &id, part_type, field_id, + /* derivative */ 0, &mask, &h_offset); /* Add the particle to the arrays */ if (flag == csds_flag_change_type || flag == csds_flag_mpi_exit || @@ -562,11 +547,11 @@ size_t csds_reader_update_state_to_next_index( /* Get the full mask */ size_t full_mask = 0; size_t h_offset = 0; - csds_loader_io_read_mask(header, log->log.map + current_offset, - &full_mask, &h_offset); + csds_loader_io_read_mask(h, log->log.map + current_offset, &full_mask, + &h_offset); /* Remove the special mask */ - full_mask = full_mask & !header->masks[csds_index_special_flags].mask; + full_mask = full_mask & !SPECIAL_FLAGS_MASK; /* Find the last offset before the current time */ size_t last_full_offset = current_offset; @@ -575,7 +560,7 @@ size_t csds_reader_update_state_to_next_index( /* Get the mask */ size_t mask = 0; h_offset = 0; - csds_loader_io_read_mask(header, log->log.map + current_offset, &mask, + csds_loader_io_read_mask(h, log->log.map + current_offset, &mask, &h_offset); /* update the offset */ @@ -586,7 +571,7 @@ size_t csds_reader_update_state_to_next_index( /* The particle should not have a special flag due to the previous loop */ - if (mask & header->masks[csds_index_special_flags].mask) { + if (mask & SPECIAL_FLAGS_MASK) { error("Found a special flag when updating the particles."); } @@ -619,7 +604,7 @@ void csds_reader_generate_index_files(const struct csds_reader *reader, int number_index) { /* Get a few pointers */ const struct csds_logfile *log = &reader->log; - const struct header *header = &log->header; + const struct header *h = &log->header; /* Write a quick message */ if (reader->verbose > 0) { @@ -632,39 +617,10 @@ void csds_reader_generate_index_files(const struct csds_reader *reader, } /* Ensure that the offset are in the assumed direction */ - if (!header_is_forward(header)) { + if (!header_is_forward(h)) { error("The offset are not in the expected direction"); } - /* Get the mask for the IDs */ - const int index_mask_id = header_get_field_index(header, "ParticleIDs"); - if (index_mask_id == -1) { - error("The logfile does not contain the field ParticleIDs"); - } - - /* Get the list of fields. */ - int global_fields_wanted = index_mask_id; - int all_fields_count[swift_type_count]; - struct field_information *all_fields[swift_type_count]; - struct field_information field_id[swift_type_count]; - for (int type = 0; type < swift_type_count; type++) { - /* TODO Implement missing particle types */ - if (type == swift_type_neutrino || type == swift_type_sink || - type == swift_type_black_hole) - continue; - - all_fields_count[type] = tools_get_number_fields(type); - all_fields[type] = (struct field_information *)malloc( - all_fields_count[type] * sizeof(struct field_information)); - tools_get_list_fields(all_fields[type], type, header); - - /* Convert fields into the local array. */ - csds_reader_get_fields_wanted(reader, &global_fields_wanted, - &field_id[type], - /* n_fields_wanted */ 1, all_fields[type], - all_fields_count[type], type); - } - /* Create the different arrays that will store the information */ struct index_writer current_state[swift_type_count]; struct index_writer parts_created[swift_type_count]; @@ -684,8 +640,8 @@ void csds_reader_generate_index_files(const struct csds_reader *reader, /* Get the initial state */ struct time_record time_record; - size_t offset = csds_reader_get_initial_state( - reader, current_state, all_fields, all_fields_count, &time_record); + size_t offset = + csds_reader_get_initial_state(reader, current_state, &time_record); /* Write the first index file */ csds_reader_write_index(reader, current_state, parts_created, parts_removed, &time_record, /* file_number */ 0); @@ -709,8 +665,8 @@ void csds_reader_generate_index_files(const struct csds_reader *reader, /* Update the state until the next index file. */ offset = csds_reader_update_state_to_next_index( - reader, offset, time_record, all_fields, all_fields_count, - current_state, parts_created, parts_removed); + reader, offset, time_record, current_state, parts_created, + parts_removed); /* Write the index file */ csds_reader_write_index(reader, current_state, parts_created, parts_removed, @@ -722,13 +678,6 @@ void csds_reader_generate_index_files(const struct csds_reader *reader, index_writer_free(¤t_state[type]); index_writer_free(&parts_created[type]); index_writer_free(&parts_removed[type]); - - /* TODO Implement missing particle types */ - if (type == swift_type_neutrino || type == swift_type_sink || - type == swift_type_black_hole) - continue; - - free(all_fields[type]); } if (reader->verbose > 0) { diff --git a/src/csds_star_formation.h b/src/csds_star_formation.h deleted file mode 100644 index c0bd8ad..0000000 --- a/src/csds_star_formation.h +++ /dev/null @@ -1,38 +0,0 @@ -/******************************************************************************* - * 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_CSDS_STARA_FORMATION_H -#define SWIFT_CSDS_STAR_FORMATION_H - -/* Config parameters. */ -#include "../config.h" - -/* Import the right functions */ -#if defined(STAR_FORMATION_NONE) -#include "./star_formation/none/csds_star_formation.h" -#elif defined(STAR_FORMATION_QLA) -#error TODO -#elif defined(STAR_FORMATION_EAGLE) -#error TODO -#elif defined(STAR_FORMATION_GEAR) -#include "./star_formation/GEAR/csds_star_formation.h" -#else -#error "Invalid choice of star formation law" -#endif - -#endif /* SWIFT_STAR_FORMATION_H */ diff --git a/src/csds_stars.h b/src/csds_stars.h deleted file mode 100644 index 43b3ecf..0000000 --- a/src/csds_stars.h +++ /dev/null @@ -1,39 +0,0 @@ -/******************************************************************************* - * 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_CSDS_STARS_H -#define SWIFT_CSDS_STARS_H - -/* Config parameters. */ -#include "../config.h" - -/* Select the correct star model */ -#if defined(STARS_NOME) -#error TODO -#elif defined(STARS_BASIC) -#include "./stars/Basic/csds_stars.h" -#elif defined(STARS_EAGLE) -#error TODO -#include "./stars/EAGLE/csds_stars.h" -#elif defined(STARS_GEAR) -#include "./stars/GEAR/csds_stars.h" -#else -#error "Invalid choice of star model" -#endif - -#endif /* SWIFT_STARS_H */ diff --git a/src/csds_time.c b/src/csds_time.c index df93f07..5dc46b1 100644 --- a/src/csds_time.c +++ b/src/csds_time.c @@ -1,5 +1,5 @@ /******************************************************************************* - * This file is part of SWIFT. + * This file is part of CSDS. * Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch) * * This program is free software: you can redistribute it and/or modify @@ -97,15 +97,8 @@ size_t time_read(struct time_record *time_record, /* read record header. */ map = csds_loader_io_read_mask(h, map + offset, &mask, &prev_offset); -#ifdef SWIFT_DEBUG_CHECKS - - /* check if time mask is present in log file header. */ - int ind = header_get_field_index(h, "Timestamp"); - if (ind == -1) error_python("File header does not contain a mask for time."); - /* check if reading a time record. */ - if (h->masks[ind].mask != mask) error_python("Not a time record."); -#endif + if (TIMESTAMP_MASK != mask) error_python("Not a time record."); /* read the record. */ map = csds_loader_io_read_data(map, sizeof(unsigned long long int), @@ -128,15 +121,10 @@ size_t time_offset_first_record(const struct header *h) { size_t offset = h->offset_first_record; char *map = h->log->log.map; - /* Check that the first record is really a time record. */ - int i = header_get_field_index(h, "Timestamp"); - - if (i == -1) error_python("Time mask not present in the log file header."); - size_t mask = 0; csds_loader_io_read_mask(h, map + offset, &mask, NULL); - if (mask != h->masks[i].mask) + if (mask != TIMESTAMP_MASK) error_python("Log file should begin by timestep."); return h->offset_first_record; @@ -400,9 +388,8 @@ size_t time_array_get_index_from_time(const struct time_array *t, if (!t) error_python("NULL pointer."); if (time < t->records[0].time || time > t->records[t->size - 1].time) - error_python("Time outside of range (%g < %g < %g).", - t->records[0].time, time, - t->records[t->size - 1].time); + error_python("Time outside of range (%g < %g < %g).", t->records[0].time, + time, t->records[t->size - 1].time); #endif /* right will contain the index at the end of the loop */ diff --git a/src/csds_time.h b/src/csds_time.h index 71ae20d..f104931 100644 --- a/src/csds_time.h +++ b/src/csds_time.h @@ -1,5 +1,5 @@ /******************************************************************************* - * This file is part of SWIFT. + * This file is part of CSDS. * Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch) * * This program is free software: you can redistribute it and/or modify diff --git a/src/csds_tools.c b/src/csds_tools.c index 8a12429..4e1078a 100644 --- a/src/csds_tools.c +++ b/src/csds_tools.c @@ -1,5 +1,5 @@ /******************************************************************************* - * This file is part of SWIFT. + * This file is part of CSDS. * Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch) * * This program is free software: you can redistribute it and/or modify @@ -25,131 +25,6 @@ #include <stdio.h> -/** - * @brief Compute the number of fields for a given particle type. - * - * @param type The type of particles. - * - * @return The number of fields. - */ -int tools_get_number_fields(enum part_type type) { - int number_fields = 0; - switch (type) { - case swift_type_gas: - number_fields = hydro_csds_field_count; - number_fields += chemistry_csds_field_part_count; - return number_fields; - - case swift_type_dark_matter: - case swift_type_dark_matter_background: - case swift_type_neutrino: - number_fields = gravity_csds_field_count; - return number_fields; - - case swift_type_stars: - number_fields = stars_csds_field_count; - number_fields += chemistry_csds_field_spart_count; - number_fields += star_formation_csds_field_count; - return number_fields; - - default: - message("Particle type %i not implemented, skipping it", type); - return 0; - } -} - -#define copy_field_to_struct_internal(MODULE, PART, TYPE) \ - for (int j = 0; j < MODULE##_csds_field##PART##_count; j++) { \ - \ - /* Save the main properties */ \ - fields[i].module = TYPE; \ - fields[i].name = MODULE##_csds_field_names##PART[j]; \ - \ - /* Get the indexes */ \ - const int global = MODULE##_csds_local_to_global##PART[j]; \ - int first = h->masks[global].reader.first_deriv; \ - int second = h->masks[global].reader.second_deriv; \ - \ - /* Save the global indexes */ \ - fields[i].global_index = global; \ - fields[i].global_index_first = first; \ - fields[i].global_index_second = second; \ - \ - /* Convert the first derivatives into local index */ \ - if (first != -1) { \ - for (int k = 0; k < MODULE##_csds_field##PART##_count; k++) { \ - if (MODULE##_csds_local_to_global##PART[k] == first) { \ - first = k; \ - break; \ - } \ - } \ - } \ - \ - /* Convert the second derivatives into local index */ \ - if (second != -1) { \ - for (int k = 0; k < MODULE##_csds_field##PART##_count; k++) { \ - if (MODULE##_csds_local_to_global##PART[k] == second) { \ - second = k; \ - break; \ - } \ - } \ - } \ - \ - /* Initialize the structure */ \ - fields[i].local_index = j; \ - fields[i].local_index_first = first; \ - fields[i].local_index_second = second; \ - i++; \ - } - -/** - * Same function as set_links_local_global_internal before but with only two - * arguments. - */ -#define copy_field_to_struct_single_particle_type(MODULE, TYPE) \ - copy_field_to_struct_internal(MODULE, , TYPE) - -/** - * Same function as set_links_local_global_internal before but with a cleaner - * argument. - */ -#define copy_field_to_struct(MODULE, PART, TYPE) \ - copy_field_to_struct_internal(MODULE, _##PART, TYPE) - -/** - * @brief Construct the list of fields for a given particle type. - * - * @param fields (output) The list of fields (need to be already allocated). - * @param type The type of particle. - * @param h The #header - */ -void tools_get_list_fields(struct field_information *fields, - enum part_type type, const struct header *h) { - int i = 0; - switch (type) { - case swift_type_gas: - copy_field_to_struct_single_particle_type(hydro, field_module_default); - copy_field_to_struct(chemistry, part, field_module_chemistry); - break; - - case swift_type_dark_matter: - case swift_type_dark_matter_background: - case swift_type_neutrino: - copy_field_to_struct_single_particle_type(gravity, field_module_default); - break; - - case swift_type_stars: - copy_field_to_struct_single_particle_type(stars, field_module_default); - copy_field_to_struct(chemistry, spart, field_module_chemistry); - copy_field_to_struct_single_particle_type(star_formation, - field_module_star_formation); - break; - - default: - error("Particle type %i not implemented", type); - } -} - /** * @brief get the offset of the next corresponding record. * @@ -275,8 +150,8 @@ size_t tools_reverse_offset(const struct header *h, char *file_map, csds_loader_io_read_mask(h, map, &prev_mask, NULL); /* Check if we are not mixing timestamp and particles */ - if ((prev_mask != h->timestamp_mask && mask == h->timestamp_mask) || - (prev_mask == h->timestamp_mask && mask != h->timestamp_mask)) + if ((prev_mask != TIMESTAMP_MASK && mask == TIMESTAMP_MASK) || + (prev_mask == TIMESTAMP_MASK && mask != TIMESTAMP_MASK)) error_python("Unexpected mask: %lu, got %lu.", mask, prev_mask); #endif // SWIFT_DEBUG_CHECKS @@ -308,9 +183,6 @@ size_t tools_check_record_consistency(const struct csds_reader *reader, size_t mask; size_t pointed_offset; - const size_t mask_special_flag = - h->masks[header_get_field_index(h, "SpecialFlags")].mask; - /* read mask + offset. */ map = csds_loader_io_read_mask(h, map, &mask, &pointed_offset); @@ -319,7 +191,7 @@ size_t tools_check_record_consistency(const struct csds_reader *reader, const size_t offset_ret = (size_t)(map - file_init); /* If something happened, skip the check. */ - if (mask & mask_special_flag) { + if (mask & SPECIAL_FLAGS_MASK) { return offset_ret; } @@ -342,8 +214,8 @@ size_t tools_check_record_consistency(const struct csds_reader *reader, csds_loader_io_read_mask(h, file_init + pointed_offset, &pointed_mask, NULL); /* check if not mixing timestamp and particles. */ - if ((pointed_mask != h->timestamp_mask && mask == h->timestamp_mask) || - (pointed_mask == h->timestamp_mask && mask != h->timestamp_mask)) + if ((pointed_mask != TIMESTAMP_MASK && mask == TIMESTAMP_MASK) || + (pointed_mask == TIMESTAMP_MASK && mask != TIMESTAMP_MASK)) error_python("Error in the offset (mask %lu at %lu != %lu at %lu).", mask, offset, pointed_mask, pointed_offset); diff --git a/src/csds_tools.h b/src/csds_tools.h index 74c2974..4b01f6a 100644 --- a/src/csds_tools.h +++ b/src/csds_tools.h @@ -1,5 +1,5 @@ /******************************************************************************* - * This file is part of SWIFT. + * This file is part of CSDS. * Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch) * * This program is free software: you can redistribute it and/or modify @@ -53,7 +53,7 @@ struct csds_reader; /** * @brief Structure dealing with reading / interpolating a field. */ -struct csds_field { +struct csds_reader_field { /* Value of the field. */ void *field; @@ -64,51 +64,6 @@ struct csds_field { void *second_deriv; }; -/** - * @brief Defines a type of modules (e.g. gravity, stars, - * chemistry, star_formation, ...) - */ -enum module_field { - field_module_default = 0, - field_module_chemistry, - field_module_star_formation, -}; - -/** - * @brief Structure containing all the information - * of a field. - */ -struct field_information { - - /* Module containing the field */ - enum module_field module; - - /* Name of the field */ - const char *name; - - /* Index of the field in the local array. */ - int local_index; - - /* Index of the field in the global array. */ - int global_index; - - /* Index of the first derivative in the local array. */ - int local_index_first; - - /* Index of the first derivative in the global array. */ - int global_index_first; - - /* Index of the second derivative in the local array. */ - int local_index_second; - - /* Index of the second derivative in the global array. */ - int global_index_second; -}; - -int tools_get_number_fields(enum part_type type); -void tools_get_list_fields(struct field_information *fields, - enum part_type type, const struct header *h); - int tools_get_next_record(const struct header *h, char *map, size_t *offset, size_t file_size); int _tools_get_next_record_backward(const struct header *h, char *map, @@ -125,6 +80,8 @@ void tools_print_progress(float percentage, int remaining_time, #ifndef HAVE_PYTHON #define error_python(s, ...) error(s, ##__VA_ARGS__); #else +#include "frameobject.h" + /** * @brief Print the python trace back */ @@ -132,11 +89,15 @@ __attribute__((always_inline)) INLINE static void csds_loader_print_traceback( void) { PyGILState_STATE gstate = PyGILState_Ensure(); - /* Import the traceback module */ - PyObject *pyth_module = PyImport_ImportModule("traceback"); - PyObject_CallMethod(pyth_module, "print_stack", ""); - Py_DECREF(pyth_module); + PyFrameObject *frame = PyEval_GetFrame(); + int lineno = PyFrame_GetLineNumber(frame); + PyObject *py_filename = frame->f_code->co_filename; + const char *filename = PyUnicode_AsUTF8(py_filename); + message("File %s, line: %i", filename, lineno); + + Py_DECREF(py_filename); + Py_DECREF(frame); PyGILState_Release(gstate); } diff --git a/src/gravity/MultiSoftening/csds_gravity.c b/src/gravity/MultiSoftening/csds_gravity.c deleted file mode 100644 index 972bd10..0000000 --- a/src/gravity/MultiSoftening/csds_gravity.c +++ /dev/null @@ -1,34 +0,0 @@ -/******************************************************************************* - * 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 this object's header */ -#include "csds_gravity.h" - -/* Local headers */ -#include "csds_tools.h" - -const int gravity_csds_field_size[gravity_csds_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_csds_local_to_global[gravity_csds_field_count]; diff --git a/src/gravity/MultiSoftening/csds_gravity.h b/src/gravity/MultiSoftening/csds_gravity.h deleted file mode 100644 index c0f4c79..0000000 --- a/src/gravity/MultiSoftening/csds_gravity.h +++ /dev/null @@ -1,141 +0,0 @@ -/******************************************************************************* - * 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_CSDS_GRAVITY_H -#define SWIFT_MULTISOFTENING_CSDS_GRAVITY_H - -#include "../config.h" - -/* local includes */ -#include "csds_interpolation.h" -#include "csds_loader_io.h" -#include "csds_parameters.h" -#include "csds_python_tools.h" -#include "gravity_csds.h" - -/* Index of the mask in the header mask array */ -extern int gravity_csds_local_to_global[gravity_csds_field_count]; - -/* Size for each mask */ -extern const int gravity_csds_field_size[gravity_csds_field_count]; - -/** - * @brief Create the link between the fields and their derivatives. - * - * @param head The #header. - */ -__attribute__((always_inline)) INLINE static void -gravity_csds_reader_link_derivatives(struct header *head) { - - /* Set the first and second derivatives */ - const int pos_id = - gravity_csds_local_to_global[gravity_csds_field_coordinates]; - const int vel_id = - gravity_csds_local_to_global[gravity_csds_field_velocities]; - const int acc_id = - gravity_csds_local_to_global[gravity_csds_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 #csds_field at a time < t. - * @param after Pointer to the #csds_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_csds_fields). - * @param params The simulation's #csds_parameters. - */ -__attribute__((always_inline)) INLINE static void -gravity_csds_interpolate_field(const double t_before, - const struct csds_field *restrict before, - const double t_after, - const struct csds_field *restrict after, - void *restrict output, const double t, - const int field, - const struct csds_parameters *params) { - -#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 - - switch (field) { - case gravity_csds_field_coordinates: - interpolate_quintic_double_float_ND(t_before, before, t_after, after, - output, t, /* dimension= */ 3, - params->periodic, params); - break; - case gravity_csds_field_velocities: - interpolate_cubic_float_ND(t_before, before, t_after, after, output, t, - /* dimension= */ 3, /* periodic= */ 0, params); - break; - case gravity_csds_field_accelerations: - interpolate_linear_float_ND(t_before, before, t_after, after, output, t, - /* dimension= */ 3, /* periodic= */ 0, - params); - break; - case gravity_csds_field_masses: - interpolate_linear_float(t_before, before, t_after, after, output, t, - /* periodic= */ 0, params); - break; - case gravity_csds_field_particle_ids: - interpolate_ids(t_before, before, t_after, after, output, t); - break; - default: - error("Not implemented"); - } -} - -#ifdef HAVE_PYTHON -__attribute__((always_inline)) INLINE static void gravity_csds_generate_python( - struct csds_python_field *fields) { - - fields[gravity_csds_field_coordinates] = - csds_loader_python_field(/* Dimension */ 3, NPY_DOUBLE); - fields[gravity_csds_field_velocities] = - csds_loader_python_field(/* Dimension */ 3, NPY_FLOAT32); - fields[gravity_csds_field_accelerations] = - csds_loader_python_field(/* Dimension */ 3, NPY_FLOAT32); - fields[gravity_csds_field_masses] = - csds_loader_python_field(/* Dimension */ 1, NPY_FLOAT32); - fields[gravity_csds_field_particle_ids] = - csds_loader_python_field(/* Dimension */ 1, NPY_LONGLONG); -} - -#endif // HAVE_PYTHON -#endif // SWIFT_MULTISOFTENING_CSDS_GRAVITY_H diff --git a/src/hydro/Gadget2/csds_hydro.c b/src/hydro/Gadget2/csds_hydro.c deleted file mode 100644 index 37b2c5c..0000000 --- a/src/hydro/Gadget2/csds_hydro.c +++ /dev/null @@ -1,37 +0,0 @@ -/******************************************************************************* - * 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 this object's header */ -#include "csds_hydro.h" - -/* Local headers */ -#include "csds_tools.h" - -int hydro_csds_local_to_global[hydro_csds_field_count]; - -const int hydro_csds_field_size[hydro_csds_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 -}; diff --git a/src/hydro/Gadget2/csds_hydro.h b/src/hydro/Gadget2/csds_hydro.h deleted file mode 100644 index 5aad453..0000000 --- a/src/hydro/Gadget2/csds_hydro.h +++ /dev/null @@ -1,147 +0,0 @@ -/******************************************************************************* - * 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_CSDS_HYDRO_H -#define SWIFT_GADGET2_CSDS_HYDRO_H - -#include "../config.h" - -/* local includes */ -#include "csds_interpolation.h" -#include "csds_loader_io.h" -#include "csds_parameters.h" -#include "csds_python_tools.h" -#include "hydro_csds.h" - -/* Index of the mask in the header mask array */ -extern int hydro_csds_local_to_global[hydro_csds_field_count]; - -/* Size for each mask */ -extern const int hydro_csds_field_size[hydro_csds_field_count]; - -/** - * @brief Create the link between the fields and their derivatives. - * - * @param head The #header. - */ -__attribute__((always_inline)) INLINE static void -hydro_csds_reader_link_derivatives(struct header *head) { - - /* Set the first and second derivatives */ - const int pos_id = hydro_csds_local_to_global[hydro_csds_field_coordinates]; - const int vel_id = hydro_csds_local_to_global[hydro_csds_field_velocities]; - const int acc_id = hydro_csds_local_to_global[hydro_csds_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 #csds_field at a time < t. - * @param after Pointer to the #csds_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_csds_fields). - * @param params The simulation's #csds_parameters. - */ -__attribute__((always_inline)) INLINE static void hydro_csds_interpolate_field( - const double t_before, const struct csds_field *restrict before, - const double t_after, const struct csds_field *restrict after, - void *restrict output, const double t, const int field, - const struct csds_parameters *params) { - -#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 - - switch (field) { - /* Do the position */ - case hydro_csds_field_coordinates: - interpolate_quintic_double_float_ND(t_before, before, t_after, after, - output, t, /* dimesion= */ 3, - params->periodic, params); - break; - /* Do the velocity */ - case hydro_csds_field_velocities: - interpolate_cubic_float_ND(t_before, before, t_after, after, output, t, - /* dimesion= */ 3, /* periodic= */ 0, params); - break; - case hydro_csds_field_accelerations: - interpolate_linear_float_ND(t_before, before, t_after, after, output, t, - /* dimesion= */ 3, /* periodic= */ 0, params); - break; - /* Do the linear interpolation of float. */ - case hydro_csds_field_smoothing_lengths: - case hydro_csds_field_entropies: - case hydro_csds_field_densities: - case hydro_csds_field_masses: - interpolate_linear_float(t_before, before, t_after, after, output, t, - /* periodic= */ 0, params); - break; - /* Check the ids */ - case hydro_csds_field_particle_ids: - interpolate_ids(t_before, before, t_after, after, output, t); - break; - default: - error("Not implemented"); - } -} - -#ifdef HAVE_PYTHON -__attribute__((always_inline)) INLINE static void hydro_csds_generate_python( - struct csds_python_field *fields) { - - fields[hydro_csds_field_coordinates] = - csds_loader_python_field(/* Dimension */ 3, NPY_DOUBLE); - fields[hydro_csds_field_velocities] = - csds_loader_python_field(/* Dimension */ 3, NPY_FLOAT32); - fields[hydro_csds_field_accelerations] = - csds_loader_python_field(/* Dimension */ 3, NPY_FLOAT32); - fields[hydro_csds_field_masses] = - csds_loader_python_field(/* Dimension */ 1, NPY_FLOAT32); - fields[hydro_csds_field_smoothing_lengths] = - csds_loader_python_field(/* Dimension */ 1, NPY_FLOAT32); - fields[hydro_csds_field_entropies] = - csds_loader_python_field(/* Dimension */ 1, NPY_FLOAT32); - fields[hydro_csds_field_particle_ids] = - csds_loader_python_field(/* Dimension */ 1, NPY_LONGLONG); - fields[hydro_csds_field_densities] = - csds_loader_python_field(/* Dimension */ 1, NPY_FLOAT32); -} - -#endif // HAVE_PYTHON -#endif /* SWIFT_GADGET2_CSDS_HYDRO_H */ diff --git a/src/hydro/SPHENIX/csds_hydro.c b/src/hydro/SPHENIX/csds_hydro.c deleted file mode 100644 index 926702f..0000000 --- a/src/hydro/SPHENIX/csds_hydro.c +++ /dev/null @@ -1,39 +0,0 @@ -/******************************************************************************* - * 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 this object's header */ -#include "csds_hydro.h" - -/* Local headers */ -#include "csds_tools.h" - -int hydro_csds_local_to_global[hydro_csds_field_count]; - -const int hydro_csds_field_size[hydro_csds_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, u), // InternalEnergies - member_size(struct part, id), // IDs - member_size(struct part, rho), // density - 7 * sizeof(float), // entropy + pressure + - // viscosity + diffusion + laplacian u + Velocity divergence + deriv -}; diff --git a/src/hydro/SPHENIX/csds_hydro.h b/src/hydro/SPHENIX/csds_hydro.h deleted file mode 100644 index 2560deb..0000000 --- a/src/hydro/SPHENIX/csds_hydro.h +++ /dev/null @@ -1,209 +0,0 @@ -/******************************************************************************* - * 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_SPHENIX_CSDS_HYDRO_H -#define SWIFT_SPHENIX_CSDS_HYDRO_H - -#include "../config.h" - -/* local includes */ -#include "csds_interpolation.h" -#include "csds_loader_io.h" -#include "csds_parameters.h" -#include "csds_python_tools.h" -#include "hydro.h" -#include "hydro_csds.h" - -/* Index of the mask in the header mask array */ -extern int hydro_csds_local_to_global[hydro_csds_field_count]; - -/* Size for each mask */ -extern const int hydro_csds_field_size[hydro_csds_field_count]; - -/** - * @brief Create the link between the fields and their derivatives. - * - * @param head The #header. - */ -__attribute__((always_inline)) INLINE static void -hydro_csds_reader_link_derivatives(struct header *head) { - - /* Set the first and second derivatives */ - const int pos_id = hydro_csds_local_to_global[hydro_csds_field_coordinates]; - const int vel_id = hydro_csds_local_to_global[hydro_csds_field_velocities]; - const int acc_id = hydro_csds_local_to_global[hydro_csds_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 #csds_field at a time < t. - * @param after Pointer to the #csds_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_csds_fields). - * @param params The simulation's #csds_parameters. - */ -__attribute__((always_inline)) INLINE static void hydro_csds_interpolate_field( - const double t_before, const struct csds_field *restrict before, - const double t_after, const struct csds_field *restrict after, - void *restrict output, const double t, const int field, - const struct csds_parameters *params) { - -#ifdef SWIFT_DEBUG_CHECKS - /* Check the times */ - if (t_before > t || t_after < t) { - error( - "The times for the interpolation are not " - "monotonically increasing %g < %g < %g.", - t_before, t, t_after); - } -#endif - - switch (field) { - /* Do the position */ - case hydro_csds_field_coordinates: - interpolate_quintic_double_float_ND(t_before, before, t_after, after, - output, t, /* dimension= */ 3, - params->periodic, params); - break; - /* Do the velocity */ - case hydro_csds_field_velocities: - interpolate_cubic_float_ND(t_before, before, t_after, after, output, t, - /* dimension= */ 3, /* periodic= */ 0, params); - break; - case hydro_csds_field_accelerations: - interpolate_linear_float_ND(t_before, before, t_after, after, output, t, - /* dimension= */ 3, /* periodic= */ 0, - params); - break; - /* Do the linear interpolation of float. */ - case hydro_csds_field_masses: - case hydro_csds_field_smoothing_lengths: - case hydro_csds_field_internal_energies: - case hydro_csds_field_densities: - interpolate_linear_float(t_before, before, t_after, after, output, t, - /* periodic= */ 0, params); - break; - /* Check the ids */ - case hydro_csds_field_particle_ids: - interpolate_ids(t_before, before, t_after, after, output, t); - break; - - case hydro_csds_field_secondary_fields: { - /* Get some variables */ - const float wa = (t - t_before) / (t_after - t_before); - const float wb = 1. - wa; - float *x = (float *)output; - const float *bef = (float *)before->field; - const float *aft = (float *)after->field; - - /* Entropy + pressure + Viscosity + Diffusion + Laplacian*/ - const int n_linear = 5; - for (int i = 0; i < n_linear; i++) { - x[i] = wa * aft[i] + wb * bef[i]; - } - - /* Div v */ - x[n_linear] = interpolate_cubic_hermite_spline( - t_before, bef[n_linear], bef[n_linear + 1], t_after, aft[n_linear], - aft[n_linear + 1], t); - - /* d Div v / dt */ - x[n_linear + 1] = wa * aft[n_linear + 1] + wb * bef[n_linear + 1]; - /* Use the linear interpolation */ - x[1] = wa * aft[n_linear + 2] + wb * bef[n_linear + 2]; - break; - } - - default: - error("Not implemented"); - } -} - -#ifdef HAVE_PYTHON -__attribute__((always_inline)) INLINE static void hydro_csds_generate_python( - struct csds_python_field *fields) { - - fields[hydro_csds_field_coordinates] = - csds_loader_python_field(/* Dimension */ 3, NPY_DOUBLE); - fields[hydro_csds_field_velocities] = - csds_loader_python_field(/* Dimension */ 3, NPY_FLOAT32); - fields[hydro_csds_field_accelerations] = - csds_loader_python_field(/* Dimension */ 3, NPY_FLOAT32); - fields[hydro_csds_field_masses] = - csds_loader_python_field(/* Dimension */ 1, NPY_FLOAT32); - fields[hydro_csds_field_smoothing_lengths] = - csds_loader_python_field(/* Dimension */ 1, NPY_FLOAT32); - fields[hydro_csds_field_internal_energies] = - csds_loader_python_field(/* Dimension */ 1, NPY_FLOAT32); - fields[hydro_csds_field_particle_ids] = - csds_loader_python_field(/* Dimension */ 1, NPY_LONGLONG); - fields[hydro_csds_field_densities] = - csds_loader_python_field(/* Dimension */ 1, NPY_FLOAT32); - - /* Now separate the secondary fields */ - fields[hydro_csds_field_secondary_fields] = - csds_loader_python_field(/* Dimension */ 7, CUSTOM_NPY_TYPE); - - csds_loader_python_field_add_subfield( - &fields[hydro_csds_field_secondary_fields], "Entropies", /* offset */ 0, - /* Dimension */ 1, NPY_FLOAT32); - - csds_loader_python_field_add_subfield( - &fields[hydro_csds_field_secondary_fields], "Pressures", - /* offset */ sizeof(float), /* Dimension */ 1, NPY_FLOAT32); - - csds_loader_python_field_add_subfield( - &fields[hydro_csds_field_secondary_fields], "ViscosityParameters", - /* offset */ 2 * sizeof(float), /* Dimension */ 1, NPY_FLOAT32); - - csds_loader_python_field_add_subfield( - &fields[hydro_csds_field_secondary_fields], "DiffusionParameters", - /* offset */ 3 * sizeof(float), /* Dimension */ 1, NPY_FLOAT32); - - csds_loader_python_field_add_subfield( - &fields[hydro_csds_field_secondary_fields], "LaplacianInternalEnergies", - /* offset */ 4 * sizeof(float), /* Dimension */ 1, NPY_FLOAT32); - - csds_loader_python_field_add_subfield( - &fields[hydro_csds_field_secondary_fields], "VelocityDivergences", - /* offset */ 5 * sizeof(float), /* Dimension */ 1, NPY_FLOAT32); - - csds_loader_python_field_add_subfield( - &fields[hydro_csds_field_secondary_fields], - "VelocityDivergenceTimeDifferentials", /* offset */ 6 * sizeof(float), - /* Dimension */ 1, NPY_FLOAT32); -} - -#endif // HAVE_PYTHON -#endif /* SWIFT_SPHENIX_CSDS_HYDRO_H */ diff --git a/src/quick_sort.c b/src/quick_sort.c index 885db88..ddda5ef 100644 --- a/src/quick_sort.c +++ b/src/quick_sort.c @@ -1,5 +1,5 @@ /******************************************************************************* - * This file is part of SWIFT. + * This file is part of CSDS. * Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch) * * This program is free software: you can redistribute it and/or modify diff --git a/src/quick_sort.h b/src/quick_sort.h index 27a9100..a5101a6 100644 --- a/src/quick_sort.h +++ b/src/quick_sort.h @@ -1,5 +1,5 @@ /******************************************************************************* - * This file is part of SWIFT. + * This file is part of CSDS. * Copyright (c) 2019 Loic Hausammann (loic.hausammann@epfl.ch) * * This program is free software: you can redistribute it and/or modify diff --git a/src/star_formation/GEAR/csds_star_formation.c b/src/star_formation/GEAR/csds_star_formation.c deleted file mode 100644 index e8e41f3..0000000 --- a/src/star_formation/GEAR/csds_star_formation.c +++ /dev/null @@ -1,32 +0,0 @@ -/******************************************************************************* - * 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 this object's header */ -#include "csds_star_formation.h" - -/* Local headers */ -#include "csds_tools.h" - -const int star_formation_csds_field_size[star_formation_csds_field_count] = { - member_size(struct spart, sf_data.birth_density) + - member_size(struct spart, sf_data.birth_mass) + - member_size(struct spart, sf_data.progenitor_id), -}; - -int star_formation_csds_local_to_global[star_formation_csds_field_count]; diff --git a/src/star_formation/GEAR/csds_star_formation.h b/src/star_formation/GEAR/csds_star_formation.h deleted file mode 100644 index 0c79210..0000000 --- a/src/star_formation/GEAR/csds_star_formation.h +++ /dev/null @@ -1,119 +0,0 @@ -/******************************************************************************* - * 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_GEAR_CSDS_STAR_FORMATION_H -#define SWIFT_GEAR_CSDS_STAR_FORMATION_H - -#include "../config.h" - -/* local includes */ -#include "csds_interpolation.h" -#include "csds_loader_io.h" -#include "csds_parameters.h" -#include "csds_python_tools.h" -#include "star_formation_csds.h" - -/* Index of the mask in the header mask array */ -extern int star_formation_csds_local_to_global[star_formation_csds_field_count]; - -/* Size for each mask */ -extern const int - star_formation_csds_field_size[star_formation_csds_field_count]; - -/** - * @brief Create the link between the fields and their derivatives. - * - * @param head The #header. - */ -__attribute__((always_inline)) INLINE static void -star_formation_csds_reader_link_derivatives(struct header *head) {} - -/** - * @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 #csds_field at a time < t. - * @param after Pointer to the #csds_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 - * #star_formation_csds_fields). - * @param params The simulation's #csds_parameters. - */ -__attribute__((always_inline)) INLINE static void -star_formation_csds_interpolate_field(const double t_before, - const struct csds_field *restrict before, - const double t_after, - const struct csds_field *restrict after, - void *restrict output, const double t, - const int field, - const struct csds_parameters *params) { - -#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 - - switch (field) { - case star_formation_csds_field_all: - interpolate_linear_float_ND(t_before, before, t_after, after, output, t, - /* Dimension */ 2, /* periodic= */ 0, params); - - /* Deal with the progenitor id */ - long long *proj_id = - (long long *)((void *)before->field + 2 * sizeof(float)); - output += 2 * sizeof(float); - *(long long *)output = *proj_id; - break; - - default: - error("Not implemented"); - } -} - -#ifdef HAVE_PYTHON -__attribute__((always_inline)) INLINE static void -star_formation_csds_generate_python(struct csds_python_field *fields) { - fields[star_formation_csds_field_all] = - csds_loader_python_field(/* Dimension */ 2, NPY_FLOAT); - - csds_loader_python_field_add_subfield( - &fields[star_formation_csds_field_all], "BirthDensities", - /* offset */ 0, /* Dimension */ 1, NPY_FLOAT32); - - csds_loader_python_field_add_subfield( - &fields[star_formation_csds_field_all], "BirthMasses", - /* offset */ sizeof(float), /* Dimension */ 1, NPY_FLOAT32); - - csds_loader_python_field_add_subfield( - &fields[star_formation_csds_field_all], "ProgenitorIDs", - /* offset */ 2 * sizeof(float), /* Dimension */ 1, NPY_LONGLONG); -} - -#endif // HAVE_PYTHON -#endif // SWIFT_GEAR_CSDS_STAR_FORMATION_H diff --git a/src/star_formation/none/csds_star_formation.c b/src/star_formation/none/csds_star_formation.c deleted file mode 100644 index 53f5b7c..0000000 --- a/src/star_formation/none/csds_star_formation.c +++ /dev/null @@ -1,28 +0,0 @@ -/******************************************************************************* - * 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 this object's header */ -#include "csds_star_formation.h" - -/* Local headers */ -#include "csds_tools.h" - -const int star_formation_csds_field_size[star_formation_csds_field_count] = {}; - -int star_formation_csds_local_to_global[star_formation_csds_field_count]; diff --git a/src/star_formation/none/csds_star_formation.h b/src/star_formation/none/csds_star_formation.h deleted file mode 100644 index 2e129a7..0000000 --- a/src/star_formation/none/csds_star_formation.h +++ /dev/null @@ -1,77 +0,0 @@ -/******************************************************************************* - * 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_NONE_CSDS_STAR_FORMATION_H -#define SWIFT_NONE_CSDS_STAR_FORMATION_H - -#include "../config.h" - -/* local includes */ -#include "csds_interpolation.h" -#include "csds_loader_io.h" -#include "csds_parameters.h" -#include "csds_python_tools.h" -#include "star_formation_csds.h" - -/* Index of the mask in the header mask array */ -extern int star_formation_csds_local_to_global[star_formation_csds_field_count]; - -/* Size for each mask */ -extern const int - star_formation_csds_field_size[star_formation_csds_field_count]; - -/** - * @brief Create the link between the fields and their derivatives. - * - * @param head The #header. - */ -__attribute__((always_inline)) INLINE static void -star_formation_csds_reader_link_derivatives(struct header *head) {} - -/** - * @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 #csds_field at a time < t. - * @param after Pointer to the #csds_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 - * #star_formation_csds_fields). - * @param params The simulation's #csds_parameters. - */ -__attribute__((always_inline)) INLINE static void -star_formation_csds_interpolate_field(const double t_before, - const struct csds_field *restrict before, - const double t_after, - const struct csds_field *restrict after, - void *restrict output, const double t, - const int field, - const struct csds_parameters *params) {} - -#ifdef HAVE_PYTHON -__attribute__((always_inline)) INLINE static void -star_formation_csds_generate_python(struct csds_python_field *fields) {} - -#endif // HAVE_PYTHON -#endif // SWIFT_NONE_CSDS_STAR_FORMATION_H diff --git a/src/stars/Basic/csds_stars.c b/src/stars/Basic/csds_stars.c deleted file mode 100644 index 6af98c0..0000000 --- a/src/stars/Basic/csds_stars.c +++ /dev/null @@ -1,35 +0,0 @@ -/******************************************************************************* - * 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 this object's header */ -#include "csds_stars.h" - -/* Local headers */ -#include "csds_tools.h" - -const int stars_csds_field_size[stars_csds_field_count] = { - member_size(struct spart, x), // coordinates - member_size(struct spart, v), // velocities - member_size(struct gpart, a_grav), // accelerations -> stored inside gparts - member_size(struct spart, mass), // massses - member_size(struct spart, h), // Smoothing Length - member_size(struct spart, id), // IDs -}; - -int stars_csds_local_to_global[stars_csds_field_count]; diff --git a/src/stars/Basic/csds_stars.h b/src/stars/Basic/csds_stars.h deleted file mode 100644 index d3dbc42..0000000 --- a/src/stars/Basic/csds_stars.h +++ /dev/null @@ -1,139 +0,0 @@ -/******************************************************************************* - * 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_BASIC_CSDS_STARS_H -#define SWIFT_BASIC_CSDS_STARS_H - -#include "../config.h" - -/* local includes */ -#include "csds_interpolation.h" -#include "csds_loader_io.h" -#include "csds_parameters.h" -#include "csds_python_tools.h" -#include "stars_csds.h" - -/* Index of the mask in the header mask array */ -extern int stars_csds_local_to_global[stars_csds_field_count]; - -/* Size for each mask */ -extern const int stars_csds_field_size[stars_csds_field_count]; - -/** - * @brief Create the link between the fields and their derivatives. - * - * @param head The #header. - */ -__attribute__((always_inline)) INLINE static void -stars_csds_reader_link_derivatives(struct header *head) { - - /* Set the first and second derivatives */ - const int pos_id = stars_csds_local_to_global[stars_csds_field_coordinates]; - const int vel_id = stars_csds_local_to_global[stars_csds_field_velocities]; - const int acc_id = stars_csds_local_to_global[stars_csds_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 #csds_field at a time < t. - * @param after Pointer to the #csds_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 - * #stars_csds_fields). - * @param params The simulation's #csds_parameters. - */ -__attribute__((always_inline)) INLINE static void stars_csds_interpolate_field( - const double t_before, const struct csds_field *restrict before, - const double t_after, const struct csds_field *restrict after, - void *restrict output, const double t, const int field, - const struct csds_parameters *params) { - -#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 - - switch (field) { - case stars_csds_field_coordinates: - interpolate_quintic_double_float_ND(t_before, before, t_after, after, - output, t, /* dimension= */ 3, - params->periodic, params); - break; - case stars_csds_field_velocities: - interpolate_cubic_float_ND(t_before, before, t_after, after, output, t, - /* dimension= */ 3, /* periodic= */ 0, params); - break; - case stars_csds_field_accelerations: - interpolate_linear_float_ND(t_before, before, t_after, after, output, t, - /* dimension= */ 3, /* periodic= */ 0, - params); - break; - case stars_csds_field_smoothing_lengths: - case stars_csds_field_masses: - interpolate_linear_float(t_before, before, t_after, after, output, t, - /* periodic= */ 0, params); - break; - case stars_csds_field_particle_ids: - interpolate_ids(t_before, before, t_after, after, output, t); - break; - default: - error("Not implemented"); - } -} - -#ifdef HAVE_PYTHON -__attribute__((always_inline)) INLINE static void stars_csds_generate_python( - struct csds_python_field *fields) { - - fields[stars_csds_field_coordinates] = - csds_loader_python_field(/* Dimension */ 3, NPY_DOUBLE); - fields[stars_csds_field_velocities] = - csds_loader_python_field(/* Dimension */ 3, NPY_FLOAT32); - fields[stars_csds_field_accelerations] = - csds_loader_python_field(/* Dimension */ 3, NPY_FLOAT32); - fields[stars_csds_field_masses] = - csds_loader_python_field(/* Dimension */ 1, NPY_FLOAT32); - fields[stars_csds_field_smoothing_lengths] = - csds_loader_python_field(/* Dimension */ 1, NPY_FLOAT32); - fields[stars_csds_field_particle_ids] = - csds_loader_python_field(/* Dimension */ 1, NPY_LONGLONG); -} - -#endif // HAVE_PYTHON - -#endif /* SWIFT_BASIC_CSDS_STARS_H */ diff --git a/src/stars/GEAR/csds_stars.c b/src/stars/GEAR/csds_stars.c deleted file mode 100644 index a80f909..0000000 --- a/src/stars/GEAR/csds_stars.c +++ /dev/null @@ -1,36 +0,0 @@ -/******************************************************************************* - * 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 this object's header */ -#include "csds_stars.h" - -/* Local headers */ -#include "csds_tools.h" - -const int stars_csds_field_size[stars_csds_field_count] = { - member_size(struct spart, x), // coordinates - member_size(struct spart, v), // velocities - member_size(struct gpart, a_grav), // accelerations -> stored inside gparts - member_size(struct spart, mass), // massses - member_size(struct spart, h), // Smoothing Length - member_size(struct spart, id), // IDs - member_size(struct spart, birth_scale_factor), // birth scale factor -}; - -int stars_csds_local_to_global[stars_csds_field_count]; diff --git a/src/stars/GEAR/csds_stars.h b/src/stars/GEAR/csds_stars.h deleted file mode 100644 index 65f7831..0000000 --- a/src/stars/GEAR/csds_stars.h +++ /dev/null @@ -1,142 +0,0 @@ -/******************************************************************************* - * 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_GEAR_CSDS_STARS_H -#define SWIFT_GEAR_CSDS_STARS_H - -#include "../config.h" - -/* local includes */ -#include "csds_interpolation.h" -#include "csds_loader_io.h" -#include "csds_parameters.h" -#include "csds_python_tools.h" -#include "stars_csds.h" - -/* Index of the mask in the header mask array */ -extern int stars_csds_local_to_global[stars_csds_field_count]; - -/* Size for each mask */ -extern const int stars_csds_field_size[stars_csds_field_count]; - -/** - * @brief Create the link between the fields and their derivatives. - * - * @param head The #header. - */ -__attribute__((always_inline)) INLINE static void -stars_csds_reader_link_derivatives(struct header *head) { - - /* Set the first and second derivatives */ - const int pos_id = stars_csds_local_to_global[stars_csds_field_coordinates]; - const int vel_id = stars_csds_local_to_global[stars_csds_field_velocities]; - const int acc_id = stars_csds_local_to_global[stars_csds_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 #csds_field at a time < t. - * @param after Pointer to the #csds_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 - * #stars_csds_fields). - * @param params The simulation's #csds_parameters. - */ -__attribute__((always_inline)) INLINE static void stars_csds_interpolate_field( - const double t_before, const struct csds_field *restrict before, - const double t_after, const struct csds_field *restrict after, - void *restrict output, const double t, const int field, - const struct csds_parameters *params) { - -#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 - - switch (field) { - case stars_csds_field_coordinates: - interpolate_quintic_double_float_ND(t_before, before, t_after, after, - output, t, /* dimension= */ 3, - params->periodic, params); - break; - case stars_csds_field_velocities: - interpolate_cubic_float_ND(t_before, before, t_after, after, output, t, - /* dimension= */ 3, /* periodic= */ 0, params); - break; - case stars_csds_field_accelerations: - interpolate_linear_float_ND(t_before, before, t_after, after, output, t, - /* dimension= */ 3, /* periodic= */ 0, - params); - break; - case stars_csds_field_smoothing_lengths: - case stars_csds_field_masses: - case stars_csds_field_birth_scale_factors: - interpolate_linear_float(t_before, before, t_after, after, output, t, - /* periodic= */ 0, params); - break; - case stars_csds_field_particle_ids: - interpolate_ids(t_before, before, t_after, after, output, t); - break; - default: - error("Not implemented"); - } -} - -#ifdef HAVE_PYTHON -__attribute__((always_inline)) INLINE static void stars_csds_generate_python( - struct csds_python_field *fields) { - - fields[stars_csds_field_coordinates] = - csds_loader_python_field(/* Dimension */ 3, NPY_DOUBLE); - fields[stars_csds_field_velocities] = - csds_loader_python_field(/* Dimension */ 3, NPY_FLOAT32); - fields[stars_csds_field_accelerations] = - csds_loader_python_field(/* Dimension */ 3, NPY_FLOAT32); - fields[stars_csds_field_masses] = - csds_loader_python_field(/* Dimension */ 1, NPY_FLOAT32); - fields[stars_csds_field_smoothing_lengths] = - csds_loader_python_field(/* Dimension */ 1, NPY_FLOAT32); - fields[stars_csds_field_particle_ids] = - csds_loader_python_field(/* Dimension */ 1, NPY_LONGLONG); - fields[stars_csds_field_birth_scale_factors] = - csds_loader_python_field(/* Dimension */ 1, NPY_FLOAT32); -} - -#endif // HAVE_PYTHON - -#endif // SWIFT_GEAR_CSDS_STARS_H diff --git a/tests/generate_log.h b/tests/generate_log.h index 270aba9..8a0ca0d 100644 --- a/tests/generate_log.h +++ b/tests/generate_log.h @@ -104,7 +104,7 @@ void write_particles(struct csds_writer *log, struct engine *e) { csds_log_timestamp(log, ti_int, ti_double, &log->timestamp_offset); /* Make sure that we have enough space in the particle csds file * to store the particles in current time step. */ - csds_ensure_size(log, nparts, /* number gpart */ 0, 0); + csds_ensure_size(log, e); /* Loop over all the particles. */ for (size_t j = 0; j < nparts; j++) { @@ -181,7 +181,7 @@ void generate_log(struct swift_params *params, struct part *parts, csds_log_timestamp(&log, e.ti_current, e.time, &log.timestamp_offset); /* Make sure that we have enough space in the particle csds file * to store the particles in current time step. */ - csds_ensure_size(&log, nparts, /* number gpart */ 0, 0); + csds_ensure_size(&log, &e); /* Log all the particles before starting */ csds_log_all_particles(&log, &e, /* first_log */ 1); diff --git a/tests/testLogfileHeader.c b/tests/testLogfileHeader.c index 3890a76..6d1e0ad 100644 --- a/tests/testLogfileHeader.c +++ b/tests/testLogfileHeader.c @@ -54,11 +54,11 @@ int main(int argc, char *argv[]) { csds_write_file_header(&log); /* Copy the masks */ - const int csds_count_mask_tmp = log.csds_count_mask; - struct mask_data *mask_data = (struct mask_data *)malloc( - log.csds_count_mask * sizeof(struct mask_data)); - memcpy(mask_data, log.csds_mask_data, - log.csds_count_mask * sizeof(struct mask_data)); + const int total_number_fields = log.total_number_fields; + struct csds_field *list_fields = (struct csds_field *)malloc( + total_number_fields * sizeof(struct csds_field)); + memcpy(list_fields, log.list_fields, + total_number_fields * sizeof(struct csds_field)); /* clean memory. */ csds_free(&log); @@ -91,36 +91,19 @@ int main(int argc, char *argv[]) { assert(h->offset_first_record == logfile->log.mmap_size); message("Checking number of masks"); - /* Get info about the masks */ - int count_mask = 0; - size_t mask = 0; - /* Skip the mask that appears multiple times */ - for (int i = 0; i < csds_count_mask_tmp; i++) { - if (mask_data[i].mask & ~mask) { - mask |= mask_data[i].mask; - count_mask += 1; - } - } - assert(h->masks_count == count_mask); - - message("Checking masks"); - mask = 0; - int j = 0; - for (int i = 0; i < count_mask; i++) { - /* skip the duplicate masks */ - if (mask_data[i].mask & ~mask) { - mask |= mask_data[i].mask; - assert(mask_data[j].size == h->masks[i].size); - assert(mask_data[j].mask == h->masks[i].mask); - assert(strcmp(mask_data[j].name, h->masks[i].name) == 0); - j += 1; - } + /* Compute the number of masks */ + int count_fields = 2; // Timestamp + special flags + for (int type = 0; type < swift_type_count; type++) { + count_fields += h->number_fields[type]; + /* Remove the copies of the special flags */ + if (h->number_fields[type] != 0) count_fields -= 1; } + assert(total_number_fields == count_fields); message("Checking offset direction"); assert(h->offset_direction == csds_offset_backward); - free(mask_data); + free(list_fields); csds_logfile_free(logfile); return 0; } diff --git a/tests/testLogfileReader.c b/tests/testLogfileReader.c index 5d33ce9..2d3f09b 100644 --- a/tests/testLogfileReader.c +++ b/tests/testLogfileReader.c @@ -54,15 +54,11 @@ void check_data(struct csds_reader *reader, struct part *parts, struct header *h = &reader->log.header; /* Create a particle */ - const int all_fields_count = tools_get_number_fields(swift_type_gas); - struct field_information *all_fields = (struct field_information *)malloc( - all_fields_count * sizeof(struct field_information)); - tools_get_list_fields(all_fields, swift_type_gas, h); - - void **output = malloc(all_fields_count * sizeof(void *)); - for (int i = 0; i < all_fields_count; i++) { - const int global = all_fields[i].global_index; - output[i] = malloc(h->masks[global].size); + void **output = malloc(h->number_fields[swift_type_gas] * sizeof(void *)); + for (int i = 0; i < h->number_fields[swift_type_gas]; i++) { + output[i] = malloc(h->fields[swift_type_gas][i].field.size); + message("%i %i %i", i, h->fields[swift_type_gas][i].field.field, + h->fields[swift_type_gas][i].field.position); } /* Define a few variables */ @@ -77,6 +73,24 @@ void check_data(struct csds_reader *reader, struct part *parts, const uint64_t id_flag = 5 * number_parts; uint64_t previous_id = id_flag; + /* Get all the fields */ + const struct field_information *field_id = + header_get_field_from_name(h, "ParticleIDs", swift_type_gas); + const struct field_information *field_pos = + header_get_field_from_name(h, "Coordinates", swift_type_gas); + const struct field_information *field_vel = + header_get_field_from_name(h, "Velocities", swift_type_gas); + const struct field_information *field_acc = + header_get_field_from_name(h, "Accelerations", swift_type_gas); + const struct field_information *field_mass = + header_get_field_from_name(h, "Masses", swift_type_gas); + const struct field_information *field_u = + header_get_field_from_name(h, "InternalEnergies", swift_type_gas); + const struct field_information *field_h = + header_get_field_from_name(h, "SmoothingLengths", swift_type_gas); + const struct field_information *field_rho = + header_get_field_from_name(h, "Densities", swift_type_gas); + /* Loop over each record. */ for (size_t offset = csds_reader_read_record(reader, output, &time, &is_particle, @@ -93,8 +107,8 @@ void check_data(struct csds_reader *reader, struct part *parts, Check that we are really increasing the id in the logfile. See the writing part to see that we are always increasing the id. */ - const uint64_t current_id = - *(long long *)output[hydro_csds_field_particle_ids]; + + const uint64_t current_id = *(uint64_t *)output[field_id->field.position]; if (previous_id != id_flag && previous_id >= current_id) { error("Wrong particle found"); previous_id = current_id; @@ -104,9 +118,9 @@ void check_data(struct csds_reader *reader, struct part *parts, if (current_id >= number_parts) error("Wrong id %li", current_id); struct part *p = &parts[current_id]; - const double *pos = (double *)output[hydro_csds_field_coordinates]; - const float *vel = (float *)output[hydro_csds_field_velocities]; - const float *acc = (float *)output[hydro_csds_field_accelerations]; + const double *pos = (double *)output[field_pos->field.position]; + const float *vel = (float *)output[field_vel->field.position]; + const float *acc = (float *)output[field_acc->field.position]; /* Check the record's data. */ for (int i = 0; i < 3; i++) { @@ -124,23 +138,22 @@ void check_data(struct csds_reader *reader, struct part *parts, assert(p->a_hydro[i] == acc[i]); } - const float energy = *(float *)output[hydro_csds_field_internal_energies]; + const float energy = *(float *)output[field_u->field.position]; assert(p->u == energy); - const float mass = *(float *)output[hydro_csds_field_masses]; + const float mass = *(float *)output[field_mass->field.position]; assert(p->mass == mass); /* Check optional fields. */ // int number_steps = step / p->time_bin; // TODO check only every few steps - const float current_h = - *(float *)output[hydro_csds_field_smoothing_lengths]; + const float current_h = *(float *)output[field_h->field.position]; assert(p->h == current_h); /* if (number_steps % period_h == 0 || step > max_step) { */ /* assert(p->h == lp.h); */ /* } else { */ /* assert(-1 == lp.h); */ /* } */ - const float rho = *(float *)output[hydro_csds_field_densities]; + const float rho = *(float *)output[field_rho->field.position]; assert(p->rho == rho); /* if (number_steps % period_rho == 0 || step > max_step) { */ /* assert(p->rho == lp.rho); */ @@ -176,11 +189,10 @@ void check_data(struct csds_reader *reader, struct part *parts, } /* Cleanup */ - for (int i = 0; i < all_fields_count; i++) { + for (int i = 0; i < h->number_fields[swift_type_gas]; i++) { free(output[i]); } free(output); - free(all_fields); } int main(int argc, char *argv[]) { diff --git a/tests/testQuickSort.c b/tests/testQuickSort.c index 85e2a13..f8ad6d4 100644 --- a/tests/testQuickSort.c +++ b/tests/testQuickSort.c @@ -93,5 +93,7 @@ int main(int argc, char *argv[]) { /* Check the results */ check_sort(data); + free(data); + return 0; } diff --git a/tests/testVirtualReality.c b/tests/testVirtualReality.c index 183e62e..2899336 100644 --- a/tests/testVirtualReality.c +++ b/tests/testVirtualReality.c @@ -66,7 +66,7 @@ int main(int argc, char *argv[]) { parser_get_param_string(¶ms, "CSDS:basename", basename); strcat(basename, "_0000"); csds_reader_init(&reader, basename, - /* Verbose */ 0, /* number_threads */ 1, + /* Verbose */ 2, /* number_threads */ 1, /* number_index */ 5); /* Read the time limits */ @@ -104,24 +104,10 @@ int main(int argc, char *argv[]) { /* Create the list of fields. */ const int n_fields = 2; - int *required_fields = (int *)malloc(n_fields * sizeof(int)); - const struct header *h = &reader.log.header; - for (int i = 0; i < n_fields; i++) { - required_fields[i] = -1; - } - for (int j = 0; j < h->masks_count; j++) { - if (strcmp(h->masks[j].name, "Coordinates") == 0) { - required_fields[0] = j; - } else if (strcmp(h->masks[j].name, "ParticleIDs") == 0) { - required_fields[1] = j; - } - } - if (required_fields[0] == -1) { - error("Coordinates not found"); - } - if (required_fields[1] == -1) { - error("ParticleIDs not found."); - } + enum field_enum *required_fields = + (enum field_enum *)malloc(n_fields * sizeof(enum field_enum)); + required_fields[0] = field_enum_coordinates; + required_fields[1] = field_enum_particles_ids; /* Create the output */ void **output = malloc(n_fields * sizeof(void *)); -- GitLab