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

Merge branch 'master' into 'mesh_force_task'

# Conflicts:
#   src/Makefile.am
parents 5dcb789c 0442c30d
......@@ -273,6 +273,18 @@ if test "$enable_debugging_checks" = "yes"; then
AC_DEFINE([SWIFT_DEBUG_CHECKS],1,[Enable expensive debugging])
fi
# Check if using our custom icbrtf is enalbled.
AC_ARG_ENABLE([custom-icbrtf],
[AS_HELP_STRING([--enable-custom-icbrtf],
[Use SWIFT's custom icbrtf function instead of the system cbrtf @<:@yes/no@:>@]
)],
[enable_custom_icbrtf="$enableval"],
[enable_custom_icbrtf="no"]
)
if test "$enable_custom_icbrtf" = "yes"; then
AC_DEFINE([WITH_ICBRTF],1,[Enable custom icbrtf])
fi
# Check whether we want to default to naive cell interactions
AC_ARG_ENABLE([naive-interactions],
[AS_HELP_STRING([--enable-naive-interactions],
......@@ -1354,5 +1366,6 @@ AC_MSG_RESULT([
Interaction debugging : $enable_debug_interactions
Naive interactions : $enable_naive_interactions
Gravity checks : $gravity_force_checks
Custom icbrtf : $enable_custom_icbrtf
------------------------])
......@@ -1957,7 +1957,7 @@ MACRO_EXPANSION = YES
# The default value is: NO.
# This tag requires that the tag ENABLE_PREPROCESSING is set to YES.
EXPAND_ONLY_PREDEF = YES
EXPAND_ONLY_PREDEF = NO
# If the SEARCH_INCLUDES tag is set to YES, the include files in the
# INCLUDE_PATH will be searched if a #include is found.
......
......@@ -31,15 +31,9 @@ SPH:
# Parameters related to the initial conditions
InitialConditions:
file_name: CoolingHalo.hdf5 # The file to read
shift_x: 0. # A shift to apply to all particles read from the ICs (in internal units).
shift_y: 0.
shift_z: 0.
# External potential parameters
IsothermalPotential:
position_x: 0. # location of centre of isothermal potential in internal units
position_y: 0.
position_z: 0.
vrot: 200. # rotation speed of isothermal potential in internal units
timestep_mult: 0.03 # controls time step
epsilon: 0.1 #softening for the isothermal potential
......
......@@ -34,9 +34,6 @@ InitialConditions:
# External potential parameters
IsothermalPotential:
position_x: 0. # Location of centre of isothermal potential in internal units
position_y: 0.
position_z: 0.
vrot: 200. # Rotation speed of isothermal potential in internal units
timestep_mult: 0.03 # Controls time step
epsilon: 1.0 # Softening for the isothermal potential
......
......@@ -31,15 +31,11 @@ SPH:
# Parameters related to the initial conditions
InitialConditions:
file_name: PointMass.hdf5 # The file to read
shift_x: 50. # A shift to apply to all particles read from the ICs (in internal units).
shift_y: 50.
shift_z: 50.
shift: [50.,50.,50.] # A shift to apply to all particles read from the ICs (in internal units).
# External potential parameters
PointMassPotential:
position_x: 50. # location of external point mass in internal units
position_y: 50.
position_z: 50.
position: [50.,50.,50.] # location of external point mass in internal units
mass: 1e10 # mass of external point mass in internal units
timestep_mult: 0.03 # controls time step
......@@ -34,9 +34,6 @@ InitialConditions:
# External potential parameters
IsothermalPotential:
position_x: 0. # location of centre of isothermal potential in internal units
position_y: 0.
position_z: 0.
vrot: 200. # rotation speed of isothermal potential in internal units
epsilon: 1.0
timestep_mult: 0.03 # controls time step
......
......@@ -26,15 +26,10 @@ Snapshots:
# Parameters related to the initial conditions
InitialConditions:
file_name: Isothermal.hdf5 # The file to read
shift_x: 200. # Shift all particles to be in the potential
shift_y: 200.
shift_z: 200.
shift: [200.,200.,200.] # Shift all particles to be in the potential
# External potential parameters
IsothermalPotential:
position_x: 0. # location of centre of isothermal potential in internal units
position_y: 0.
position_z: 0.
vrot: 200. # rotation speed of isothermal potential in internal units
timestep_mult: 0.01 # controls time step
epsilon: 0. # No softening at the centre of the halo
vrot: 200. # rotation speed of isothermal potential in internal units
timestep_mult: 0.01 # controls time step
epsilon: 0. # No softening at the centre of the halo
......@@ -32,15 +32,10 @@ SPH:
# Parameters related to the initial conditions
InitialConditions:
file_name: initial_conditions.hdf5 # The file to read
shift_x: 0. # A shift to apply to all particles read from the ICs (in internal units).
shift_y: 0.
shift_z: 0.
# External potential parameters
PointMassPotential:
position_x: 5. # location of external point mass in internal units
position_y: 5. # here just take this to be the centre of the ring
position_z: 5.
mass: 1.498351e7 # mass of external point mass in internal units
timestep_mult: 0.03 # controls time step
position: [5.,5.,5.] # location of external point mass in internal units
mass: 1.498351e7 # mass of external point mass in internal units
timestep_mult: 0.03 # controls time step
softening: 0.05
......@@ -300,7 +300,7 @@ def write_block(f, part_type, pos, vel, ids, mass, int_energy, smoothing, other=
else:
data = default_data
particles = f.create_group(f"PartType{part_type}")
particles = f.create_group("PartType" + str(part_type))
for name, value in data.items():
particles.create_dataset(name, data=value)
......
......@@ -35,8 +35,6 @@ InitialConditions:
# External potential parameters
PointMassPotential:
position_x: 50. # location of external point mass in internal units
position_y: 50.
position_z: 50.
position: [50.,50.,50.] # location of external point mass in internal units
mass: 1e10 # mass of external point mass in internal units
timestep_mult: 1e-2
......@@ -97,9 +97,7 @@ InitialConditions:
cleanup_velocity_factors: 0 # (Optional) Clean up the scale-factors used in the definition of the velocity variable in the ICs (e.g. in Gadget files).
cleanup_smoothing_lengths: 0 # (Optional) Clean the values of the smoothing lengths that are read in to remove stupid values. Set to 1 to activate.
smoothing_length_scaling: 1. # (Optional) A scaling factor to apply to all smoothing lengths in the ICs.
shift_x: 0. # (Optional) A shift to apply to all particles read from the ICs (in internal units).
shift_y: 0.
shift_z: 0.
shift: [0.0,0.0,0.0] # (Optional) A shift to apply to all particles read from the ICs (in internal units).
replicate: 2 # (Optional) Replicate all particles along each axis a given integer number of times. Default 1.
# Parameters controlling restarts
......@@ -116,9 +114,7 @@ Restarts:
DomainDecomposition:
initial_type: simple_metis # (Optional) The initial decomposition strategy: "grid",
# "simple_metis", "weighted_metis", or "vectorized".
initial_grid_x: 10 # (Optional) Grid size if the "grid" strategy is chosen.
initial_grid_y: 10 # ""
initial_grid_z: 10 # ""
initial_grid: [10,10,10] # (Optional) Grid sizes if the "grid" strategy is chosen.
repartition_type: costs/costs # (Optional) The re-decomposition strategy, one of:
# "none/none", "costs/costs", "counts/none", "none/costs", "counts/costs",
......@@ -149,17 +145,14 @@ EoS:
# Point mass external potentials
PointMassPotential:
position_x: 50. # location of external point mass (internal units)
position_y: 50.
position_z: 50.
mass: 1e10 # mass of external point mass (internal units)
timestep_mult: 0.03 # Dimensionless pre-factor for the time-step condition
position: [50.,50.0,50.] # location of external point mass (internal units)
mass: 1e10 # mass of external point mass (internal units)
timestep_mult: 0.03 # Dimensionless pre-factor for the time-step condition
softening: 0.05 # For point-mass-softened option
# Isothermal potential parameters
IsothermalPotential:
position_x: 100. # Location of centre of isothermal potential with respect to centre of the box (internal units)
position_y: 100.
position_z: 100.
position: [100.,100.,100.] # Location of centre of isothermal potential with respect to centre of the box (internal units)
vrot: 200. # Rotation speed of isothermal potential (internal units)
timestep_mult: 0.03 # Dimensionless pre-factor for the time-step condition
epsilon: 0.1 # Softening size (internal units)
......
......@@ -48,7 +48,7 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
dump.h logger.h active.h timeline.h xmf.h gravity_properties.h gravity_derivatives.h \
gravity_softened_derivatives.h vector_power.h collectgroup.h hydro_space.h sort_part.h \
chemistry.h chemistry_io.h chemistry_struct.h cosmology.h restart.h space_getsid.h utilities.h \
mesh_gravity.h
mesh_gravity.h cbrt.h
# Common source files
AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
......
......@@ -33,6 +33,7 @@
#include <math.h>
/* Local headers. */
#include "cbrt.h"
#include "error.h"
#include "inline.h"
......@@ -112,8 +113,13 @@ __attribute__((always_inline, const)) INLINE static float pow_gamma(float x) {
#if defined(HYDRO_GAMMA_5_3)
#ifdef WITH_ICBRTF
const float icbrt = icbrtf(x); /* x^(-1/3) */
return icbrt * x * x; /* x^(5/3) */
#else
const float cbrt = cbrtf(x); /* x^(1/3) */
return cbrt * cbrt * x; /* x^(5/3) */
#endif // WITH_ICBRTF
#elif defined(HYDRO_GAMMA_7_5)
......@@ -121,7 +127,12 @@ __attribute__((always_inline, const)) INLINE static float pow_gamma(float x) {
#elif defined(HYDRO_GAMMA_4_3)
#ifdef WITH_ICBRTF
const float icbrt = icbrtf(x); /* x^(-1/3) */
return icbrt * icbrt * x * x; /* x^(4/3) */
#else
return cbrtf(x) * x; /* x^(4/3) */
#endif // WITH_ICBRTF
#elif defined(HYDRO_GAMMA_2_1)
......@@ -146,8 +157,13 @@ __attribute__((always_inline, const)) INLINE static float pow_gamma_minus_one(
#if defined(HYDRO_GAMMA_5_3)
#ifdef WITH_ICBRTF
const float icbrt = icbrtf(x); /* x^(-1/3) */
return x * icbrt; /* x^(2/3) */
#else
const float cbrt = cbrtf(x); /* x^(1/3) */
return cbrt * cbrt; /* x^(2/3) */
#endif // WITH_ICBRTF
#elif defined(HYDRO_GAMMA_7_5)
......@@ -155,7 +171,12 @@ __attribute__((always_inline, const)) INLINE static float pow_gamma_minus_one(
#elif defined(HYDRO_GAMMA_4_3)
#ifdef WITH_ICBRTF
const float icbrt = icbrtf(x); /* x^(-1/3) */
return x * icbrt * icbrt; /* x^(1/3) */
#else
return cbrtf(x); /* x^(1/3) */
#endif // WITH_ICBRTF
#elif defined(HYDRO_GAMMA_2_1)
......@@ -180,8 +201,13 @@ pow_minus_gamma_minus_one(float x) {
#if defined(HYDRO_GAMMA_5_3)
#ifdef WITH_ICBRTF
const float icbrt = icbrtf(x); /* x^(-1/3) */
return icbrt * icbrt; /* x^(-2/3) */
#else
const float cbrt_inv = 1.f / cbrtf(x); /* x^(-1/3) */
return cbrt_inv * cbrt_inv; /* x^(-2/3) */
#endif // WITH_ICBRTF
#elif defined(HYDRO_GAMMA_7_5)
......@@ -189,7 +215,11 @@ pow_minus_gamma_minus_one(float x) {
#elif defined(HYDRO_GAMMA_4_3)
#ifdef WITH_ICBRTF
return icbrtf(x); /* x^(-1/3) */
#else
return 1.f / cbrtf(x); /* x^(-1/3) */
#endif // WITH_ICBRTF
#elif defined(HYDRO_GAMMA_2_1)
......@@ -217,9 +247,15 @@ __attribute__((always_inline, const)) INLINE static float pow_minus_gamma(
#if defined(HYDRO_GAMMA_5_3)
#ifdef WITH_ICBRTF
const float icbrt = icbrtf(x); /* x^(-1/3) */
const float icbrt2 = icbrt * icbrt; /* x^(-2/3) */
return icbrt * icbrt2 * icbrt2; /* x^(-5/3) */
#else
const float cbrt_inv = 1.f / cbrtf(x); /* x^(-1/3) */
const float cbrt_inv2 = cbrt_inv * cbrt_inv; /* x^(-2/3) */
return cbrt_inv * cbrt_inv2 * cbrt_inv2; /* x^(-5/3) */
#endif // WITH_ICBRTF
#elif defined(HYDRO_GAMMA_7_5)
......@@ -227,7 +263,11 @@ __attribute__((always_inline, const)) INLINE static float pow_minus_gamma(
#elif defined(HYDRO_GAMMA_4_3)
#ifdef WITH_ICBRTF
const float cbrt_inv = icbrtf(x); /* x^(-1/3) */
#else
const float cbrt_inv = 1.f / cbrtf(x); /* x^(-1/3) */
#endif // WITH_ICBRTF
const float cbrt_inv2 = cbrt_inv * cbrt_inv; /* x^(-2/3) */
return cbrt_inv2 * cbrt_inv2; /* x^(-4/3) */
......
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2016 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
* Matthieu Schaller (matthieu.schaller@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU 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_CBRT_H
#define SWIFT_CBRT_H
/* Config parameters. */
#include "../config.h"
/* Some standard headers. */
#include <math.h>
/* Local headers. */
#include "inline.h"
/**
* @brief Compute the inverse cube root of a single-precision floating-point
* number.
*
* This function does not care about non-finite inputs.
*
* @warning This function is faster than both gcc and Intel's `cbrtf()`
* functions on x86 systems. However, Other compilers or other architectures
* may have faster implementations of the standard function `cbrtf()` that
* will potentionally outperform this function.
*
* @param x_in The input value.
*
* @return The inverse cubic root of @c x_in (i.e. \f$x_{in}^{-1/3} \f$) .
*/
__attribute__((always_inline)) INLINE static float icbrtf(float x_in) {
union {
float as_float;
unsigned int as_uint;
int as_int;
} cast;
/* Extract the exponent. */
cast.as_float = x_in;
const int exponent = ((cast.as_int & 0x7f800000) >> 23) - 127;
/* Clear the exponent and sign to get the mantissa. */
cast.as_uint = (cast.as_uint & ~0xff800000) | 0x3f800000;
const float x_norm = cast.as_float;
/* Multiply by sqrt(1/2) and subtract one, should then be in the
range [sqrt(1/2) - 1, sqrt(2) - 1). */
const float x = x_norm * (float)M_SQRT1_2 - 1.0f;
/* Compute the polynomial interpolant. */
float res =
9.99976591940035e-01f +
x * (-3.32901212909283e-01f +
x * (2.24361110929912e-01f +
x * (-1.88913279594895e-01f + x * 1.28384036492344e-01f)));
/* Compute the new exponent and the correction factor. */
int exponent_new = exponent;
if (exponent_new < 0) exponent_new -= 2;
exponent_new = -exponent_new / 3;
const int exponent_rem = exponent + 3 * exponent_new;
cast.as_uint = (exponent_new + 127) << 23;
static const float scale[3] = {8.90898718140339e-01f, 7.07106781186548e-01f,
5.61231024154687e-01f};
const float exponent_scale = cast.as_float * scale[exponent_rem];
/* Scale the result and set the correct sign. */
res = copysignf(res * exponent_scale, x_in);
/* One step of Newton iteration to refine the result. */
res *= (1.0f / 3.0f) * (4.0f - x_in * res * res * res);
/* We're done. */
return res;
}
#endif /* SWIFT_CBRT_H */
......@@ -250,13 +250,18 @@ void io_write_attribute_s(hid_t grp, const char* name, const char* str) {
/**
* @brief Reads the Unit System from an IC file.
*
* If the 'Units' group does not exist in the ICs, we will use the internal
* system of units.
*
* @param h_file The (opened) HDF5 file from which to read.
* @param us The unit_system to fill.
* @param ic_units The unit_system to fill.
* @param internal_units The internal system of units to copy if needed.
* @param mpi_rank The MPI rank we are on.
*
* If the 'Units' group does not exist in the ICs, cgs units will be assumed
*/
void io_read_unit_system(hid_t h_file, struct unit_system* us, int mpi_rank) {
void io_read_unit_system(hid_t h_file, struct unit_system* ic_units,
const struct unit_system* internal_units,
int mpi_rank) {
/* First check if it exists as this is *not* required. */
const htri_t exists = H5Lexists(h_file, "/Units", H5P_DEFAULT);
......@@ -264,16 +269,12 @@ void io_read_unit_system(hid_t h_file, struct unit_system* us, int mpi_rank) {
if (exists == 0) {
if (mpi_rank == 0)
message("'Units' group not found in ICs. Assuming CGS unit system.");
message("'Units' group not found in ICs. Assuming internal unit system.");
/* Default to CGS */
us->UnitMass_in_cgs = 1.;
us->UnitLength_in_cgs = 1.;
us->UnitTime_in_cgs = 1.;
us->UnitCurrent_in_cgs = 1.;
us->UnitTemperature_in_cgs = 1.;
units_copy(ic_units, internal_units);
return;
} else if (exists < 0) {
error("Serious problem with 'Units' group in ICs. H5Lexists gives %d",
exists);
......@@ -284,15 +285,15 @@ void io_read_unit_system(hid_t h_file, struct unit_system* us, int mpi_rank) {
/* Ok, Read the damn thing */
io_read_attribute(h_grp, "Unit length in cgs (U_L)", DOUBLE,
&us->UnitLength_in_cgs);
&ic_units->UnitLength_in_cgs);
io_read_attribute(h_grp, "Unit mass in cgs (U_M)", DOUBLE,
&us->UnitMass_in_cgs);
&ic_units->UnitMass_in_cgs);
io_read_attribute(h_grp, "Unit time in cgs (U_t)", DOUBLE,
&us->UnitTime_in_cgs);
&ic_units->UnitTime_in_cgs);
io_read_attribute(h_grp, "Unit current in cgs (U_I)", DOUBLE,
&us->UnitCurrent_in_cgs);
&ic_units->UnitCurrent_in_cgs);
io_read_attribute(h_grp, "Unit temperature in cgs (U_T)", DOUBLE,
&us->UnitTemperature_in_cgs);
&ic_units->UnitTemperature_in_cgs);
/* Clean up */
H5Gclose(h_grp);
......
......@@ -75,7 +75,9 @@ void io_write_attribute_s(hid_t grp, const char* name, const char* str);
void io_write_code_description(hid_t h_file);
void io_write_engine_policy(hid_t h_file, const struct engine* e);
void io_read_unit_system(hid_t h_file, struct unit_system* us, int mpi_rank);
void io_read_unit_system(hid_t h_file, struct unit_system* ic_units,
const struct unit_system* internal_units,
int mpi_rank);
void io_write_unit_system(hid_t h_grp, const struct unit_system* us,
const char* groupName);
......
......@@ -54,7 +54,7 @@ enum eos_planetary_type_id {
eos_planetary_type_HM80 = 2,
eos_planetary_type_ANEOS = 3,
eos_planetary_type_SESAME = 4,
} __attribute__((packed));
};
/**
* @brief Minor type for the planetary equation of state.
......@@ -104,7 +104,7 @@ enum eos_planetary_material_id {
/*! SESAME iron */
eos_planetary_id_SESAME_iron =
eos_planetary_type_SESAME * eos_planetary_type_factor,
} __attribute__((packed));
};
/* Individual EOS function headers. */
#include "aneos.h"
......@@ -132,7 +132,8 @@ __attribute__((always_inline)) INLINE static float
gas_internal_energy_from_entropy(float density, float entropy,
enum eos_planetary_material_id mat_id) {
const enum eos_planetary_type_id type = mat_id / eos_planetary_type_factor;
const enum eos_planetary_type_id type =
(enum eos_planetary_type_id)(mat_id / eos_planetary_type_factor);
/* Select the material base type */
switch (type) {
......@@ -241,7 +242,8 @@ gas_internal_energy_from_entropy(float density, float entropy,
__attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
float density, float entropy, enum eos_planetary_material_id mat_id) {
const enum eos_planetary_type_id type = mat_id / eos_planetary_type_factor;
const enum eos_planetary_type_id type =
(enum eos_planetary_type_id)(mat_id / eos_planetary_type_factor);
/* Select the material base type */
switch (type) {
......@@ -344,7 +346,8 @@ __attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
__attribute__((always_inline)) INLINE static float gas_entropy_from_pressure(
float density, float P, enum eos_planetary_material_id mat_id) {
const enum eos_planetary_type_id type = mat_id / eos_planetary_type_factor;
const enum eos_planetary_type_id type =
(enum eos_planetary_type_id)(mat_id / eos_planetary_type_factor);
/* Select the material base type */
switch (type) {
......@@ -445,7 +448,8 @@ __attribute__((always_inline)) INLINE static float gas_entropy_from_pressure(
__attribute__((always_inline)) INLINE static float gas_soundspeed_from_entropy(
float density, float entropy, enum eos_planetary_material_id mat_id) {
const enum eos_planetary_type_id type = mat_id / eos_planetary_type_factor;
const enum eos_planetary_type_id type =
(enum eos_planetary_type_id)(mat_id / eos_planetary_type_factor);
/* Select the material base type */
switch (type) {
......@@ -549,8 +553,8 @@ __attribute__((always_inline)) INLINE static float gas_soundspeed_from_entropy(
__attribute__((always_inline)) INLINE static float
gas_entropy_from_internal_energy(float density, float u,
enum eos_planetary_material_id mat_id) {
const enum eos_planetary_type_id type = mat_id / eos_planetary_type_factor;
const enum eos_planetary_type_id type =
(enum eos_planetary_type_id)(mat_id / eos_planetary_type_factor);
/* Select the material base type */
switch (type) {
......@@ -654,7 +658,8 @@ __attribute__((always_inline)) INLINE static float
gas_pressure_from_internal_energy(float density, float u,
enum eos_planetary_material_id mat_id) {
const enum eos_planetary_type_id type = mat_id / eos_planetary_type_factor;
const enum eos_planetary_type_id type =
(enum eos_planetary_type_id)(mat_id / eos_planetary_type_factor);
/* Select the material base type */
switch (type) {
......@@ -762,7 +767,8 @@ __attribute__((always_inline)) INLINE static float
gas_internal_energy_from_pressure(float density, float P,
enum eos_planetary_material_id mat_id) {
const enum eos_planetary_type_id type = mat_id / eos_planetary_type_factor;
const enum eos_planetary_type_id type =
(enum eos_planetary_type_id)(mat_id / eos_planetary_type_factor);
/* Select the material base type */
switch (type) {
......@@ -867,7 +873,8 @@ __attribute__((always_inline)) INLINE static float
gas_soundspeed_from_internal_energy(float density, float u,
enum eos_planetary_material_id mat_id) {
const enum eos_planetary_type_id type = mat_id / eos_planetary_type_factor;
const enum eos_planetary_type_id type =
(enum eos_planetary_type_id)(mat_id / eos_planetary_type_factor);
/* Select the material base type */
switch (type) {
......@@ -975,7 +982,8 @@ gas_soundspeed_from_internal_energy(float density, float u,
__attribute__((always_inline)) INLINE static float gas_soundspeed_from_pressure(
float density, float P, enum eos_planetary_material_id mat_id) {
const enum eos_planetary_type_id type = mat_id / eos_planetary_type_factor;
const enum eos_planetary_type_id type =
(enum eos_planetary_type_id)(mat_id / eos_planetary_type_factor);
/* Select the material base type */
switch (type) {
......
......@@ -69,8 +69,8 @@ INLINE static void hydro_read_particles(struct part* parts,
UNIT_CONV_ACCELERATION, parts, a_hydro);
list[7] = io_make_input_field("Density", FLOAT, 1, OPTIONAL,
UNIT_CONV_DENSITY, parts, rho);
list[8] =
io_make_input_field("MaterialID", INT, 1, OPTIONAL, 1, parts, mat_id);
list[8] = io_make_input_field("MaterialID", INT, 1, COMPULSORY,
UNIT_CONV_NO_UNITS, parts, mat_id);
}
INLINE static void convert_S(const struct engine* e, const struct part* p,
......
......@@ -730,7 +730,7 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
struct unit_system* ic_units =
(struct unit_system*)malloc(sizeof(struct unit_system));
if (ic_units == NULL) error("Unable to allocate memory for IC unit system");
io_read_unit_system(h_file, ic_units, mpi_rank);
io_read_unit_system(h_file, ic_units, internal_units, mpi_rank);
/* Tell the user if a conversion will be needed */
if (mpi_rank == 0) {
......