Commit 07f1e5b1 authored by Peter W. Draper's avatar Peter W. Draper
Browse files

Merge branch 'cosmo_io' into 'master'

More cosmology work

Closes #407 and #404

See merge request !516
parents a727dda0 4f946aa2
......@@ -278,7 +278,7 @@ struct engine {
const struct hydro_props *hydro_properties;
/* Properties of the self-gravity scheme */
const struct gravity_props *gravity_properties;
struct gravity_props *gravity_properties;
/* Properties of external gravitational potential */
const struct external_potential *external_potential;
......@@ -335,7 +335,7 @@ void engine_init(
long long Ngas, long long Ndm, int policy, int verbose,
struct repartition *reparttype, const struct unit_system *internal_units,
const struct phys_const *physical_constants, struct cosmology *cosmo,
const struct hydro_props *hydro, const struct gravity_props *gravity,
const struct hydro_props *hydro, struct gravity_props *gravity,
const struct external_potential *potential,
const struct cooling_function_data *cooling_func,
const struct chemistry_data *chemistry, struct sourceterms *sourceterms);
......
......@@ -24,16 +24,17 @@
/* Local headers. */
#include "const.h"
#include "engine.h"
#include "inline.h"
#include "part.h"
#include "space.h"
/* So far only one model here */
/* Straight-forward import */
#include "./gravity/Default/gravity.h"
#include "./gravity/Default/gravity_iact.h"
struct engine;
struct space;
void gravity_exact_force_ewald_init(double boxSize);
void gravity_exact_force_ewald_free();
void gravity_exact_force_compute(struct space *s, const struct engine *e);
......
......@@ -21,8 +21,10 @@
#define SWIFT_DEFAULT_GRAVITY_H
#include <float.h>
#include "cosmology.h"
#include "gravity_properties.h"
#include "kernel_gravity.h"
#include "minmax.h"
/**
......@@ -40,11 +42,12 @@ __attribute__((always_inline)) INLINE static float gravity_get_mass(
* @brief Returns the softening of a particle
*
* @param gp The particle of interest
* @param grav_props The global gravity properties.
*/
__attribute__((always_inline)) INLINE static float gravity_get_softening(
const struct gpart* restrict gp) {
const struct gpart* gp, const struct gravity_props* restrict grav_props) {
return gp->epsilon;
return grav_props->epsilon_cur;
}
/**
......@@ -102,12 +105,10 @@ gravity_compute_timestep_self(const struct gpart* const gp,
const float ac_inv = (ac2 > 0.f) ? 1.f / sqrtf(ac2) : FLT_MAX;
// MATTHIEU cosmological evolution of the softening?
const float epsilon = gravity_get_softening(gp);
const float epsilon = gravity_get_softening(gp, grav_props);
/* Note that 0.66666667 = 2. (from Gadget) / 3. (Plummer softening) */
const float dt =
sqrtf(0.66666667f * cosmo->a * grav_props->eta * epsilon * ac_inv);
const float dt = sqrtf(2. * kernel_gravity_softening_plummer_equivalent_inv *
cosmo->a * grav_props->eta * epsilon * ac_inv);
return dt;
}
......@@ -183,7 +184,6 @@ __attribute__((always_inline)) INLINE static void gravity_first_init_gpart(
struct gpart* gp, const struct gravity_props* grav_props) {
gp->time_bin = 0;
gp->epsilon = grav_props->epsilon;
gravity_init_gpart(gp);
}
......
......@@ -22,9 +22,9 @@
__attribute__((always_inline)) INLINE static void gravity_debug_particle(
const struct gpart* p) {
printf(
"mass=%.3e epsilon=%.5e time_bin=%d\n"
"mass=%.3e time_bin=%d\n"
"x=[%.5e,%.5e,%.5e], v_full=[%.5e,%.5e,%.5e], a=[%.5e,%.5e,%.5e]\n",
p->mass, p->epsilon, p->time_bin, p->x[0], p->x[1], p->x[2], p->v_full[0],
p->mass, p->time_bin, p->x[0], p->x[1], p->x[2], p->v_full[0],
p->v_full[1], p->v_full[2], p->a_grav[0], p->a_grav[1], p->a_grav[2]);
#ifdef SWIFT_DEBUG_CHECKS
printf("num_interacted=%lld ti_drift=%lld ti_kick=%lld\n", p->num_interacted,
......
......@@ -35,6 +35,38 @@ void convert_gpart_pos(const struct engine* e, const struct gpart* gp,
}
}
void convert_gpart_vel(const struct engine* e, const struct gpart* gp,
float* ret) {
const int with_cosmology = (e->policy & engine_policy_cosmology);
const struct cosmology* cosmo = e->cosmology;
const integertime_t ti_current = e->ti_current;
const double time_base = e->time_base;
const integertime_t ti_beg = get_integer_time_begin(ti_current, gp->time_bin);
const integertime_t ti_end = get_integer_time_end(ti_current, gp->time_bin);
/* Get time-step since the last kick */
float dt_kick_grav;
if (with_cosmology) {
dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, ti_beg, ti_current);
dt_kick_grav -=
cosmology_get_grav_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
} else {
dt_kick_grav = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
}
/* Extrapolate the velocites to the current time */
ret[0] = gp->v_full[0] + gp->a_grav[0] * dt_kick_grav;
ret[1] = gp->v_full[1] + gp->a_grav[1] * dt_kick_grav;
ret[2] = gp->v_full[2] + gp->a_grav[2] * dt_kick_grav;
/* Conversion from internal units to peculiar velocities */
ret[0] *= cosmo->a2_inv;
ret[1] *= cosmo->a2_inv;
ret[2] *= cosmo->a2_inv;
}
/**
* @brief Specifies which g-particle fields to read from a dataset
*
......@@ -69,20 +101,20 @@ void darkmatter_read_particles(struct gpart* gparts, struct io_props* list,
void darkmatter_write_particles(const struct gpart* gparts,
struct io_props* list, int* num_fields) {
/* Say how much we want to read */
/* Say how much we want to write */
*num_fields = 5;
/* List what we want to read */
/* List what we want to write */
list[0] = io_make_output_field_convert_gpart(
"Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, gparts, convert_gpart_pos);
list[1] = io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED,
gparts, v_full);
list[1] = io_make_output_field_convert_gpart(
"Velocities", FLOAT, 3, UNIT_CONV_SPEED, gparts, convert_gpart_vel);
list[2] =
io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, gparts, mass);
list[3] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
UNIT_CONV_NO_UNITS, gparts, id_or_neg_offset);
list[4] = io_make_output_field("Acceleration", FLOAT, 3,
UNIT_CONV_ACCELERATION, gparts, a_grav);
list[4] = io_make_output_field("Potential", FLOAT, 1, UNIT_CONV_POTENTIAL,
gparts, potential);
}
#endif /* SWIFT_DEFAULT_GRAVITY_IO_H */
......@@ -22,35 +22,32 @@
/* Gravity particle. */
struct gpart {
/* Particle ID. If negative, it is the negative offset of the #part with
/*! Particle ID. If negative, it is the negative offset of the #part with
which this gpart is linked. */
long long id_or_neg_offset;
/* Particle position. */
/*! Particle position. */
double x[3];
/* Offset between current position and position at last tree rebuild. */
/*! Offset between current position and position at last tree rebuild. */
float x_diff[3];
/* Particle velocity. */
/*! Particle velocity. */
float v_full[3];
/* Particle acceleration. */
/*! Particle acceleration. */
float a_grav[3];
/* Particle mass. */
/*! Particle mass. */
float mass;
/* Gravitational potential */
/*! Gravitational potential */
float potential;
/* Softening length */
float epsilon;
/* Time-step length */
/*! Time-step length */
timebin_t time_bin;
/* Type of the #gpart (DM, gas, star, ...) */
/*! Type of the #gpart (DM, gas, star, ...) */
enum part_type type;
#ifdef SWIFT_DEBUG_CHECKS
......
......@@ -149,14 +149,16 @@ static INLINE void gravity_cache_init(struct gravity_cache *c, int count) {
* @param shift A shift to apply to all the particles.
* @param CoM The position of the multipole.
* @param r_max2 The square of the multipole radius.
* @param theta_crit2 The square of the opening angle.
* @param cell The cell we play with (to get reasonable padding positions).
* @param grav_props The global gravity properties.
*/
__attribute__((always_inline)) INLINE static void gravity_cache_populate(
timebin_t max_active_bin, struct gravity_cache *c,
const struct gpart *restrict gparts, int gcount, int gcount_padded,
const double shift[3], const float CoM[3], float r_max2, float theta_crit2,
const struct cell *cell) {
const double shift[3], const float CoM[3], float r_max2,
const struct cell *cell, const struct gravity_props *grav_props) {
const float theta_crit2 = grav_props->theta_crit2;
/* Make the compiler understand we are in happy vectorization land */
swift_declare_aligned_ptr(float, x, c->x, SWIFT_CACHE_ALIGNMENT);
......@@ -174,7 +176,7 @@ __attribute__((always_inline)) INLINE static void gravity_cache_populate(
x[i] = (float)(gparts[i].x[0] - shift[0]);
y[i] = (float)(gparts[i].x[1] - shift[1]);
z[i] = (float)(gparts[i].x[2] - shift[2]);
epsilon[i] = gparts[i].epsilon;
epsilon[i] = gravity_get_softening(&gparts[i], grav_props);
m[i] = gparts[i].mass;
active[i] = (int)(gparts[i].time_bin <= max_active_bin);
......@@ -220,13 +222,15 @@ __attribute__((always_inline)) INLINE static void gravity_cache_populate(
* multiple of the vector length.
* @param shift A shift to apply to all the particles.
* @param cell The cell we play with (to get reasonable padding positions).
* @param grav_props The global gravity properties.
*/
__attribute__((always_inline)) INLINE static void
gravity_cache_populate_no_mpole(timebin_t max_active_bin,
struct gravity_cache *c,
const struct gpart *restrict gparts, int gcount,
int gcount_padded, const double shift[3],
const struct cell *cell) {
const struct cell *cell,
const struct gravity_props *grav_props) {
/* Make the compiler understand we are in happy vectorization land */
swift_declare_aligned_ptr(float, x, c->x, SWIFT_CACHE_ALIGNMENT);
......@@ -242,7 +246,7 @@ gravity_cache_populate_no_mpole(timebin_t max_active_bin,
x[i] = (float)(gparts[i].x[0] - shift[0]);
y[i] = (float)(gparts[i].x[1] - shift[1]);
z[i] = (float)(gparts[i].x[2] - shift[2]);
epsilon[i] = gparts[i].epsilon;
epsilon[i] = gravity_get_softening(&gparts[i], grav_props);
m[i] = gparts[i].mass;
active[i] = (int)(gparts[i].time_bin <= max_active_bin);
}
......
......@@ -37,7 +37,8 @@
#define gravity_props_default_r_cut_min 0.1f
void gravity_props_init(struct gravity_props *p,
const struct swift_params *params) {
const struct swift_params *params,
const struct cosmology *cosmo) {
/* Tree-PM parameters */
p->a_smooth = parser_get_opt_param_float(params, "Gravity:a_smooth",
......@@ -56,11 +57,35 @@ void gravity_props_init(struct gravity_props *p,
p->theta_crit2 = p->theta_crit * p->theta_crit;
p->theta_crit_inv = 1. / p->theta_crit;
/* Softening lengths */
p->epsilon = 3. * parser_get_param_double(params, "Gravity:epsilon");
p->epsilon2 = p->epsilon * p->epsilon;
p->epsilon_inv = 1.f / p->epsilon;
p->epsilon_inv3 = p->epsilon_inv * p->epsilon_inv * p->epsilon_inv;
/* Softening parameters */
p->epsilon_comoving =
parser_get_param_double(params, "Gravity:comoving_softening");
p->epsilon_max_physical =
parser_get_param_double(params, "Gravity:max_physical_softening");
/* Set the softening to the current time */
gravity_update(p, cosmo);
}
void gravity_update(struct gravity_props *p, const struct cosmology *cosmo) {
/* Current softening lengths */
double softening;
if (p->epsilon_comoving * cosmo->a > p->epsilon_max_physical)
softening = p->epsilon_max_physical / cosmo->a;
else
softening = p->epsilon_comoving;
/* Plummer equivalent -> internal */
softening *= kernel_gravity_softening_plummer_equivalent;
/* Store things */
p->epsilon_cur = softening;
/* Other factors */
p->epsilon_cur2 = softening * softening;
p->epsilon_cur_inv = 1. / softening;
p->epsilon_cur_inv3 = 1. / (softening * softening * softening);
}
void gravity_props_print(const struct gravity_props *p) {
......@@ -72,8 +97,17 @@ void gravity_props_print(const struct gravity_props *p) {
message("Self-gravity opening angle: theta=%.4f", p->theta_crit);
message("Self-gravity softening: epsilon=%.4f (Plummer equivalent: %.4f)",
p->epsilon, p->epsilon / 3.);
message(
"Self-gravity comoving softening: epsilon=%.4f (Plummer equivalent: "
"%.4f)",
p->epsilon_comoving * kernel_gravity_softening_plummer_equivalent,
p->epsilon_comoving);
message(
"Self-gravity maximal physical softening: epsilon=%.4f (Plummer "
"equivalent: %.4f)",
p->epsilon_max_physical * kernel_gravity_softening_plummer_equivalent,
p->epsilon_max_physical);
message("Self-gravity mesh smoothing-scale: a_smooth=%f", p->a_smooth);
......@@ -86,9 +120,18 @@ void gravity_props_print_snapshot(hid_t h_grpgrav,
const struct gravity_props *p) {
io_write_attribute_f(h_grpgrav, "Time integration eta", p->eta);
io_write_attribute_f(h_grpgrav, "Softening length", p->epsilon);
io_write_attribute_f(h_grpgrav, "Softening length (Plummer equivalent)",
p->epsilon / 3.);
io_write_attribute_f(
h_grpgrav, "Comoving softening length",
p->epsilon_comoving * kernel_gravity_softening_plummer_equivalent);
io_write_attribute_f(h_grpgrav,
"Comoving Softening length (Plummer equivalent)",
p->epsilon_comoving);
io_write_attribute_f(
h_grpgrav, "Maximal physical softening length",
p->epsilon_max_physical * kernel_gravity_softening_plummer_equivalent);
io_write_attribute_f(h_grpgrav,
"Maximal physical softening length (Plummer equivalent)",
p->epsilon_max_physical);
io_write_attribute_f(h_grpgrav, "Opening angle", p->theta_crit);
io_write_attribute_d(h_grpgrav, "MM order", SELF_GRAVITY_MULTIPOLE_ORDER);
io_write_attribute_f(h_grpgrav, "Mesh a_smooth", p->a_smooth);
......
......@@ -27,6 +27,7 @@
#endif
/* Local includes. */
#include "cosmology.h"
#include "parser.h"
#include "restart.h"
......@@ -58,22 +59,30 @@ struct gravity_props {
/*! Inverse of opening angle */
double theta_crit_inv;
/*! Softening length */
float epsilon;
/*! Comoving softening */
double epsilon_comoving;
/*! Square of softening length */
float epsilon2;
/*! Maxium physical softening */
double epsilon_max_physical;
/*! Inverse of softening length */
float epsilon_inv;
/*! Current sftening length */
float epsilon_cur;
/*! Cube of the inverse of softening length */
float epsilon_inv3;
/*! Square of current softening length */
float epsilon_cur2;
/*! Inverse of current softening length */
float epsilon_cur_inv;
/*! Cube of the inverse of current oftening length */
float epsilon_cur_inv3;
};
void gravity_props_print(const struct gravity_props *p);
void gravity_props_init(struct gravity_props *p,
const struct swift_params *params);
const struct swift_params *params,
const struct cosmology *cosmo);
void gravity_update(struct gravity_props *p, const struct cosmology *cosmo);
#if defined(HAVE_HDF5)
void gravity_props_print_snapshot(hid_t h_grpsph,
......
......@@ -160,6 +160,18 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
return p->mass;
}
/**
* @brief Sets the mass of a particle
*
* @param p The particle of interest
* @param m The mass to set.
*/
__attribute__((always_inline)) INLINE static void hydro_set_mass(
struct part *restrict p, float m) {
p->mass = m;
}
/**
* @brief Returns the velocities drifted to the current time of a particle.
*
......@@ -492,7 +504,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
* @param dt_therm The time-step for this kick (for thermodynamic quantities)
*/
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
struct part *restrict p, struct xpart *restrict xp, float dt_therm) {}
struct part *restrict p, struct xpart *restrict xp, float dt_therm,
const struct cosmology *cosmo, const struct hydro_props *hydro_props) {}
/**
* @brief Converts hydro quantity of a particle at the start of a run
......@@ -500,9 +513,12 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
* Requires the density to be known
*
* @param p The particle to act upon
* @param xp The extended data.
* @param cosmo The cosmological model.
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part *restrict p, struct xpart *restrict xp) {}
struct part *restrict p, struct xpart *restrict xp,
const struct cosmology *cosmo) {}
/**
* @brief Initialises the particles for the first time
......@@ -529,4 +545,21 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
hydro_init_part(p, NULL);
}
/**
* @brief Overwrite the initial internal energy of a particle.
*
* Note that in the cases where the thermodynamic variable is not
* internal energy but gets converted later, we must overwrite that
* field. The conversion to the actual variable happens later after
* the initial fake time-step.
*
* @param p The #part to write to.
* @param u_init The new initial internal energy.
*/
__attribute__((always_inline)) INLINE static void
hydro_set_init_internal_energy(struct part *p, float u_init) {
p->u = u_init;
}
#endif /* SWIFT_DEFAULT_HYDRO_H */
......@@ -54,7 +54,7 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
}
void convert_part_pos(const struct engine* e, const struct part* p,
double* ret) {
const struct xpart* xp, double* ret) {
if (e->s->periodic) {
ret[0] = box_wrap(p->x[0], 0.0, e->s->dim[0]);
......@@ -67,6 +67,49 @@ void convert_part_pos(const struct engine* e, const struct part* p,
}
}
void convert_part_vel(const struct engine* e, const struct part* p,
const struct xpart* xp, float* ret) {
const int with_cosmology = (e->policy & engine_policy_cosmology);
const struct cosmology* cosmo = e->cosmology;
const integertime_t ti_current = e->ti_current;
const double time_base = e->time_base;
const integertime_t ti_beg = get_integer_time_begin(ti_current, p->time_bin);
const integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin);
/* Get time-step since the last kick */
float dt_kick_grav, dt_kick_hydro;
if (with_cosmology) {
dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, ti_beg, ti_current);
dt_kick_grav -=
cosmology_get_grav_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
dt_kick_hydro = cosmology_get_hydro_kick_factor(cosmo, ti_beg, ti_current);
dt_kick_hydro -=
cosmology_get_hydro_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
} else {
dt_kick_grav = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
dt_kick_hydro = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
}
/* Extrapolate the velocites to the current time */
hydro_get_drifted_velocities(p, xp, dt_kick_hydro, dt_kick_grav, ret);
/* Conversion from internal units to peculiar velocities */
ret[0] *= cosmo->a2_inv;
ret[1] *= cosmo->a2_inv;
ret[2] *= cosmo->a2_inv;
}
void convert_part_potential(const struct engine* e, const struct part* p,
const struct xpart* xp, float* ret) {
if (p->gpart != NULL)
ret[0] = gravity_get_comoving_potential(p->gpart);
else
ret[0] = 0.f;
}
/**
* @brief Specifies which particle fields to write to a dataset
*
......@@ -74,16 +117,17 @@ void convert_part_pos(const struct engine* e, const struct part* p,
* @param list The list of i/o properties to write.
* @param num_fields The number of i/o fields to write.
*/
void hydro_write_particles(struct part* parts, struct io_props* list,
int* num_fields) {
void hydro_write_particles(const struct part* parts, const struct xpart* xparts,
struct io_props* list, int* num_fields) {
*num_fields = 7;
*num_fields = 8;
/* List what we want to write */
list[0] = io_make_output_field_convert_part(
"Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, parts, convert_part_pos);
list[1] =
io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, v);
list[0] = io_make_output_field_convert_part("Coordinates", DOUBLE, 3,
UNIT_CONV_LENGTH, parts, xparts,
convert_part_pos);
list[1] = io_make_output_field_convert_part(
"Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, xparts, convert_part_vel);
list[2] =
io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts, mass);
list[3] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
......@@ -94,6 +138,9 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
UNIT_CONV_NO_UNITS, parts, id);
list[6] =
io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts, rho);
list[7] = io_make_output_field_convert_part("Potential", FLOAT, 1,
UNIT_CONV_POTENTIAL, parts,
xparts, convert_part_potential);
}
/**
......
......@@ -170,6 +170,18 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
return p->mass;
}
/**
* @brief Sets the mass of a particle
*
* @param p The particle of interest
* @param m The mass to set.
*/
__attribute__((always_inline)) INLINE static void hydro_set_mass(
struct part *restrict p, float m) {
p->mass = m;
}
/**
* @brief Returns the velocities drifted to the current time of a particle.
*
......@@ -513,19 +525,29 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
* @param p The particle to act upon
* @param xp The particle extended data to act upon
* @param dt_therm The time-step for this kick (for thermodynamic quantities)
* @param cosmo The cosmological model.
* @param hydro_props The constants used in the scheme
*/
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
struct part *restrict p, struct xpart *restrict xp, float dt_therm) {
struct part *restrict p, struct xpart *restrict xp, float dt_therm,
const struct cosmology *cosmo, const struct hydro_props *hydro_props) {
/* Do not decrease the entropy by more than a factor of 2 */
if (dt_therm > 0. && p->entropy_dt * dt_therm < -0.5f * xp->entropy_full) {
/* message("Warning! Limiting entropy_dt. Possible cooling error.\n
* entropy_full = %g \n entropy_dt * dt =%g \n", */
/* xp->entropy_full,p->entropy_dt * dt); */
p->entropy_dt = -0.5f * xp->entropy_full / dt_therm;
}
xp->entropy_full += p->entropy_dt * dt_therm;
/* Apply the minimal energy limit */
const float density = p->rho * cosmo->a3_inv;
const float min_energy = hydro_props->minimal_internal_energy;
const float min_entropy =
gas_entropy_from_internal_energy(density, min_energy);
if (xp->entropy_full < min_entropy) {