Commit 0442c30d authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'yaml_arrays' into 'master'

Yaml arrays

See merge request !564
parents b2308cc0 aa3105fe
......@@ -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
......@@ -95,9 +95,7 @@ InitialConditions:
cleanup_h_factors: 0 # (Optional) Clean up the h-factors used 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
......@@ -114,9 +112,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",
......@@ -147,17 +143,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)
......
This diff is collapsed.
......@@ -80,6 +80,32 @@ double parser_get_opt_param_double(struct swift_params *params,
const char *name, double def);
void parser_get_opt_param_string(struct swift_params *params, const char *name,
char *retParam, const char *def);
void parser_get_param_char_array(struct swift_params *params, const char *name,
int nval, char *values);
void parser_get_param_int_array(struct swift_params *params, const char *name,
int nval, int *values);
void parser_get_param_float_array(struct swift_params *params, const char *name,
int nval, float *values);
void parser_get_param_double_array(struct swift_params *params,
const char *name, int nval, double *values);
int parser_get_opt_param_char_array(struct swift_params *params,
const char *name, int nval, char *values);
int parser_get_opt_param_int_array(struct swift_params *params,
const char *name, int nval, int *values);
int parser_get_opt_param_float_array(struct swift_params *params,
const char *name, int nval, float *values);
int parser_get_opt_param_double_array(struct swift_params *params,
const char *name, int nval,
double *values);
void parser_get_param_string_array(struct swift_params *params,
const char *name, int *nval, char ***values);
int parser_get_opt_param_string_array(struct swift_params *params,
const char *name, int *nval,
char ***values, int ndef,
const char *def[]);
void parser_free_param_string_array(int nval, char **values);
#if defined(HAVE_HDF5)
void parser_write_params_to_hdf5(const struct swift_params *params, hid_t grp,
......
......@@ -1095,12 +1095,8 @@ void partition_init(struct partition *partition,
/* In case of grid, read more parameters */
if (part_type[0] == 'g') {
partition->grid[0] = parser_get_opt_param_int(
params, "DomainDecomposition:initial_grid_x", partition->grid[0]);
partition->grid[1] = parser_get_opt_param_int(
params, "DomainDecomposition:initial_grid_y", partition->grid[1]);
partition->grid[2] = parser_get_opt_param_int(
params, "DomainDecomposition:initial_grid_z", partition->grid[2]);
parser_get_opt_param_int_array(params, "DomainDecomposition:initial_grid",
3, partition->grid);
}
/* Now let's check what the user wants as a repartition strategy */
......
......@@ -42,7 +42,7 @@
struct external_potential {
/*! Position of the centre of potential */
double x, y, z;
double x[3];
/*! Rotation velocity */
double vrot;
......@@ -73,9 +73,9 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep(
const struct phys_const* restrict phys_const,
const struct gpart* restrict g) {
const float dx = g->x[0] - potential->x;
const float dy = g->x[1] - potential->y;
const float dz = g->x[2] - potential->z;
const float dx = g->x[0] - potential->x[0];
const float dy = g->x[1] - potential->x[1];
const float dz = g->x[2] - potential->x[2];
const float r2_plus_epsilon2_inv =
1.f / (dx * dx + dy * dy + dz * dz + potential->epsilon2);
......@@ -115,9 +115,9 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
double time, const struct external_potential* potential,
const struct phys_const* const phys_const, struct gpart* g) {
const float dx = g->x[0] - potential->x;
const float dy = g->x[1] - potential->y;
const float dz = g->x[2] - potential->z;
const float dx = g->x[0] - potential->x[0];
const float dy = g->x[1] - potential->x[1];
const float dz = g->x[2] - potential->x[2];
const float r2_plus_epsilon2_inv =
1.f / (dx * dx + dy * dy + dz * dz + potential->epsilon2);
......@@ -144,9 +144,9 @@ external_gravity_get_potential_energy(
double time, const struct external_potential* potential,
const struct phys_const* const phys_const, const struct gpart* g) {
const float dx = g->x[0] - potential->x;
const float dy = g->x[1] - potential->y;
const float dz = g->x[2] - potential->z;
const float dx = g->x[0] - potential->x[0];
const float dy = g->x[1] - potential->x[1];
const float dz = g->x[2] - potential->x[2];
return -0.5f * potential->vrot * potential->vrot *
logf(dx * dx + dy * dy + dz * dz + potential->epsilon2);
......@@ -166,15 +166,12 @@ static INLINE void potential_init_backend(
const struct unit_system* us, const struct space* s,
struct external_potential* potential) {
potential->x =
s->dim[0] / 2. +
parser_get_param_double(parameter_file, "IsothermalPotential:position_x");
potential->y =
s->dim[1] / 2. +
parser_get_param_double(parameter_file, "IsothermalPotential:position_y");
potential->z =
s->dim[2] / 2. +
parser_get_param_double(parameter_file, "IsothermalPotential:position_z");
parser_get_param_double_array(parameter_file, "IsothermalPotential:position",
3, potential->x);
potential->x[0] += s->dim[0] / 2.;
potential->x[1] += s->dim[1] / 2.;
potential->x[2] += s->dim[2] / 2.;
potential->vrot =
parser_get_param_double(parameter_file, "IsothermalPotential:vrot");
potential->timestep_mult = parser_get_param_float(
......@@ -198,7 +195,7 @@ static INLINE void potential_print_backend(
"External potential is 'Isothermal' with properties are (x,y,z) = (%e, "
"%e, %e), vrot = %e "
"timestep multiplier = %e, epsilon = %e",
potential->x, potential->y, potential->z, potential->vrot,
potential->x[0], potential->x[1], potential->x[2], potential->vrot,
potential->timestep_mult, sqrtf(potential->epsilon2));
}
......
......@@ -41,7 +41,7 @@
struct external_potential {
/*! Position of the point mass */
double x, y, z;
double x[3];
/*! Mass */
double mass;
......@@ -66,15 +66,15 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep(
const struct gpart* restrict g) {
const float G_newton = phys_const->const_newton_G;
const float dx = g->x[0] - potential->x;
const float dy = g->x[1] - potential->y;
const float dz = g->x[2] - potential->z;
const float dx = g->x[0] - potential->x[0];
const float dy = g->x[1] - potential->x[1];
const float dz = g->x[2] - potential->x[2];
const float rinv = 1.f / sqrtf(dx * dx + dy * dy + dz * dz);
const float rinv2 = rinv * rinv;
const float rinv3 = rinv2 * rinv;
const float drdv = (g->x[0] - potential->x) * (g->v_full[0]) +
(g->x[1] - potential->y) * (g->v_full[1]) +
(g->x[2] - potential->z) * (g->v_full[2]);
const float drdv = (g->x[0] - potential->x[0]) * (g->v_full[0]) +
(g->x[1] - potential->x[1]) * (g->v_full[1]) +
(g->x[2] - potential->x[2]) * (g->v_full[2]);
const float dota_x = G_newton * potential->mass * rinv3 *
(-g->v_full[0] + 3.f * rinv2 * drdv * dx);
const float dota_y = G_newton * potential->mass * rinv3 *
......@@ -109,9 +109,9 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
double time, const struct external_potential* restrict potential,
const struct phys_const* restrict phys_const, struct gpart* restrict g) {
const float dx = g->x[0] - potential->x;
const float dy = g->x[1] - potential->y;
const float dz = g->x[2] - potential->z;
const float dx = g->x[0] - potential->x[0];
const float dy = g->x[1] - potential->x[1];
const float dz = g->x[2] - potential->x[2];
const float rinv = 1.f / sqrtf(dx * dx + dy * dy + dz * dz);
const float rinv3 = rinv * rinv * rinv;
......@@ -134,9 +134,9 @@ external_gravity_get_potential_energy(
double time, const struct external_potential* potential,
const struct phys_const* const phys_const, const struct gpart* g) {
const float dx = g->x[0] - potential->x;
const float dy = g->x[1] - potential->y;
const float dz = g->x[2] - potential->z;
const float dx = g->x[0] - potential->x[0];
const float dy = g->x[1] - potential->x[1];
const float dz = g->x[2] - potential->x[2];
const float rinv = 1. / sqrtf(dx * dx + dy * dy + dz * dz);
return -phys_const->const_newton_G * potential->mass * rinv;
}
......@@ -156,12 +156,8 @@ static INLINE void potential_init_backend(
const struct unit_system* us, const struct space* s,
struct external_potential* potential) {
potential->x =
parser_get_param_double(parameter_file, "PointMassPotential:position_x");
potential->y =
parser_get_param_double(parameter_file, "PointMassPotential:position_y");
potential->z =
parser_get_param_double(parameter_file, "PointMassPotential:position_z");
parser_get_param_double_array(parameter_file, "PointMassPotential:position",
3, potential->x);
potential->mass =
parser_get_param_double(parameter_file, "PointMassPotential:mass");
potential->timestep_mult = parser_get_param_float(
......@@ -179,7 +175,7 @@ static INLINE void potential_print_backend(
message(
"External potential is 'Point mass' with properties (x,y,z) = (%e, %e, "
"%e), M = %e timestep multiplier = %e.",
potential->x, potential->y, potential->z, potential->mass,
potential->x[0], potential->x[1], potential->x[2], potential->mass,
potential->timestep_mult);
}
......
......@@ -41,7 +41,7 @@
struct external_potential {
/*! Position of the point mass */
double x, y, z;
double x[3];
/*! Mass */
double mass;
......@@ -68,16 +68,16 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep(
const struct gpart* restrict g) {
const float G_newton = phys_const->const_newton_G;
const float dx = g->x[0] - potential->x;
const float dy = g->x[1] - potential->y;
const float dz = g->x[2] - potential->z;
const float dx = g->x[0] - potential->x[0];
const float dy = g->x[1] - potential->x[1];
const float dz = g->x[2] - potential->x[2];
const float r = sqrtf(dx * dx + dy * dy + dz * dz);
const float rinv = 1.f / r;
const float rinv2 = rinv * rinv;
const float rinv3 = rinv2 * rinv;
const float drdv = (g->x[0] - potential->x) * (g->v_full[0]) +
(g->x[1] - potential->y) * (g->v_full[1]) +
(g->x[2] - potential->z) * (g->v_full[2]);
const float drdv = (g->x[0] - potential->x[0]) * (g->v_full[0]) +
(g->x[1] - potential->x[1]) * (g->v_full[1]) +
(g->x[2] - potential->x[2]) * (g->v_full[2]);
float factor;
if (r < 0.175) {
......@@ -129,9 +129,9 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
double time, const struct external_potential* restrict potential,
const struct phys_const* restrict phys_const, struct gpart* restrict g) {
const float dx = g->x[0] - potential->x;
const float dy = g->x[1] - potential->y;
const float dz = g->x[2] - potential->z;
const float dx = g->x[0] - potential->x[0];
const float dy = g->x[1] - potential->x[1];
const float dz = g->x[2] - potential->x[2];
const float r = sqrtf(dx * dx + dy * dy + dz * dz);
const float rinv = 1.f / r;
const float rinv2 = rinv * rinv;
......@@ -174,9 +174,9 @@ external_gravity_get_potential_energy(
double time, const struct external_potential* potential,
const struct phys_const* const phys_const, const struct gpart* g) {
const float dx = g->x[0] - potential->x;
const float dy = g->x[1] - potential->y;
const float dz = g->x[2] - potential->z;
const float dx = g->x[0] - potential->x[0];
const float dy = g->x[1] - potential->x[1];
const float dz = g->x[2] - potential->x[2];
const float rinv = 1. / sqrtf(dx * dx + dy * dy + dz * dz);
return -phys_const->const_newton_G * potential->mass * rinv;
}
......@@ -196,12 +196,8 @@ static INLINE void potential_init_backend(
const struct unit_system* us, const struct space* s,
struct external_potential* potential) {
potential->x =
parser_get_param_double(parameter_file, "PointMassPotential:position_x");
potential->y =
parser_get_param_double(parameter_file, "PointMassPotential:position_y");
potential->z =
parser_get_param_double(parameter_file, "PointMassPotential:position_z");
parser_get_param_double_array(parameter_file, "PointMassPotential:position",
3, potential->x);
potential->mass =
parser_get_param_double(parameter_file, "PointMassPotential:mass");
potential->timestep_mult = parser_get_param_float(
......@@ -219,7 +215,7 @@ static INLINE void potential_print_backend(
message(
"External potential is 'Point mass' with properties (x,y,z) = (%e, %e, "
"%e), M = %e timestep multiplier = %e.",
potential->x, potential->y, potential->z, potential->mass,
potential->x[0], potential->x[1], potential->x[2], potential->mass,
potential->timestep_mult);
}
......
......@@ -50,7 +50,7 @@
struct external_potential {
/*! Position of the point mass */
double x, y, z;
double x[3];
/*! Mass */
double mass;
......@@ -77,9 +77,9 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep(
const struct phys_const* restrict phys_const,
const struct gpart* restrict g) {
const float dx = g->x[0] - potential->x;
const float dy = g->x[1] - potential->y;
const float dz = g->x[2] - potential->z;
const float dx = g->x[0] - potential->x[0];
const float dy = g->x[1] - potential->x[1];
const float dz = g->x[2] - potential->x[2];
const float softening2 = potential->softening * potential->softening;
const float r2 = dx * dx + dy * dy + dz * dz;
......@@ -133,9 +133,9 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
double time, const struct external_potential* restrict potential,
const struct phys_const* restrict phys_const, struct gpart* restrict g) {
const float dx = g->x[0] - potential->x;
const float dy = g->x[1] - potential->y;
const float dz = g->x[2] - potential->z;
const float dx = g->x[0] - potential->x[0];
const float dy = g->x[1] - potential->x[1];
const float dz = g->x[2] - potential->x[2];
const float rinv = 1.f / sqrtf(dx * dx + dy * dy + dz * dz +
potential->softening * potential->softening);
const float rinv3 = rinv * rinv * rinv;
......@@ -159,9 +159,9 @@ external_gravity_get_potential_energy(
double time, const struct external_potential* potential,
const struct phys_const* const phys_const, const struct gpart* g) {
const float dx = g->x[0] - potential->x;
const float dy = g->x[1] - potential->y;
const float dz = g->x[2] - potential->z;
const float dx = g->x[0] - potential->x[0];
const float dy = g->x[1] - potential->x[1];
const float dz = g->x[2] - potential->x[2];
const float rinv = 1. / sqrtf(dx * dx + dy * dy + dz * dz +
potential->softening * potential->softening);
......@@ -183,12 +183,8 @@ static INLINE void potential_init_backend(
const struct unit_system* us, const struct space* s,
struct external_potential* potential) {
potential->x =
parser_get_param_double(parameter_file, "PointMassPotential:position_x");
potential->y =
parser_get_param_double(parameter_file, "PointMassPotential:position_y");
potential->z =
parser_get_param_double(parameter_file, "PointMassPotential:position_z");
parser_get_param_double_array(parameter_file, "PointMassPotential:position",
3, potential->x);
potential->mass =
parser_get_param_double(parameter_file, "PointMassPotential:mass");
potential->timestep_mult = parser_get_param_float(
......@@ -208,7 +204,7 @@ static INLINE void potential_print_backend(
message(
"External potential is 'Point mass' with properties (x,y,z) = (%e, %e, "
"%e), M = %e timestep multiplier = %e, softening = %e.",
potential->x, potential->y, potential->z, potential->mass,
potential->x[0], potential->x[1], potential->x[2], potential->mass,
potential->timestep_mult, potential->softening);
}
......
......@@ -2758,12 +2758,8 @@ void space_init(struct space *s, struct swift_params *params,
/* Apply shift */
double shift[3] = {0.0, 0.0, 0.0};
shift[0] =
parser_get_opt_param_double(params, "InitialConditions:shift_x", 0.0);
shift[1] =
parser_get_opt_param_double(params, "InitialConditions:shift_y", 0.0);
shift[2] =
parser_get_opt_param_double(params, "InitialConditions:shift_z", 0.0);
parser_get_opt_param_double_array(params, "InitialConditions:shift", 3,
shift);
if ((shift[0] != 0. || shift[1] != 0. || shift[2] != 0.) && !dry_run) {
message("Shifting particles by [%e %e %e]", shift[0], shift[1], shift[2]);
for (size_t k = 0; k < Npart; k++) {
......
......@@ -23,6 +23,7 @@
#include "../config.h"
/* Some standard headers. */
#include <ctype.h>
#include <math.h>
#include <stddef.h>
#include <stdio.h>
......@@ -728,3 +729,46 @@ long get_maxrss() {
getrusage(RUSAGE_SELF, &usage);
return usage.ru_maxrss;
}
/**
* @brief trim leading white space from a string.
*
* Returns pointer to first character.
*
* @param s the string.
* @result the result.
*/
char *trim_leading(char *s) {
if (s == NULL || strlen(s) < 2) return s;
while (isspace(*s)) s++;
return s;
}
/**
* @brief trim trailing white space from a string.
*
* Modifies the string by adding a NULL to the end.