Skip to content
GitLab
Commit c8be4935 authored by Pedro Gonnet's avatar Pedro Gonnet
Browse files

Merge branch 'new_timeline_mpi' of gitlab.cosma.dur.ac.uk:swift/swiftsim into new_timeline_mpi 

Conflicts:
	src/task.c
parents 73da4be8 8e3f584c
Hide whitespace changes
Inline Side-by-side
examples/CoolingHaloWithSpin/cooling_halo.yml
View file @ c8be4935
......@@ -10,7 +10,7 @@ InternalUnitSystem:
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
time_end: 10. # The end time of the simulation (in internal units).
dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
dt_min: 1e-10 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-1 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the conserved quantities statistics
......@@ -35,12 +35,12 @@ InitialConditions:
# External potential parameters
IsothermalPotential:
position_x: 0. # location of centre of isothermal potential in internal units
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
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
# Cooling parameters
LambdaCooling:
......@@ -48,4 +48,4 @@ LambdaCooling:
minimum_temperature: 1.0e4 # Minimal temperature (Kelvin)
mean_molecular_weight: 0.59 # Mean molecular weight
hydrogen_mass_abundance: 0.75 # Hydrogen mass abundance (dimensionless)
cooling_tstep_mult: 0.1 # Dimensionless pre-factor for the time-step condition
cooling_tstep_mult: 1.0 # Dimensionless pre-factor for the time-step condition
examples/CoolingHaloWithSpin/run.sh
View file @ c8be4935
......@@ -4,7 +4,8 @@
echo "Generating initial conditions for the isothermal potential box example..."
python makeIC.py 10000
../swift -g -s -C -t 16 cooling_halo.yml 2>&1 | tee output.log
# Run SWIFT with external potential, SPH and cooling
../swift -g -s -C -t 1 cooling_halo.yml 2>&1 | tee output.log
# python radial_profile.py 10
......
src/debug.c
View file @ c8be4935
......@@ -194,7 +194,8 @@ int checkSpacehmax(struct space *s) {
/**
* @brief Check if the h_max and dx_max values of a cell's hierarchy are
* consistent with the particles. Report verbosely if not.
* consistent with the particles. Also checks if particles are correctly
* in a cell. Report verbosely if not.
*
* @param c the top cell of the hierarchy.
* @param depth the recursion depth for use in messages. Set to 0 initially.
......@@ -206,43 +207,62 @@ int checkCellhdxmax(const struct cell *c, int *depth) {
float h_max = 0.0f;
float dx_max = 0.0f;
if (!c->split) {
const size_t nr_parts = c->count;
struct part *parts = c->parts;
for (size_t k = 0; k < nr_parts; k++) {
h_max = (h_max > parts[k].h) ? h_max : parts[k].h;
int result = 1;
const double loc_min[3] = {c->loc[0], c->loc[1], c->loc[2]};
const double loc_max[3] = {c->loc[0] + c->width[0], c->loc[1] + c->width[1],
c->loc[2] + c->width[2]};
const size_t nr_parts = c->count;
struct part *parts = c->parts;
struct xpart *xparts = c->xparts;
for (size_t k = 0; k < nr_parts; k++) {
struct part *const p = &parts[k];
struct xpart *const xp = &xparts[k];
if (p->x[0] < loc_min[0] || p->x[0] > loc_max[0] || p->x[1] < loc_min[1] ||
p->x[1] > loc_max[1] || p->x[2] < loc_min[2] || p->x[2] > loc_max[2]) {
message(
"Inconsistent part position p->x=[%e %e %e], c->loc=[%e %e %e] "
"c->width=[%e %e %e]",
p->x[0], p->x[1], p->x[2], c->loc[0], c->loc[1], c->loc[2],
c->width[0], c->width[1], c->width[2]);
result = 0;
}
} else {
for (int k = 0; k < 8; k++)
const float dx2 = xp->x_diff[0] * xp->x_diff[0] +
xp->x_diff[1] * xp->x_diff[1] +
xp->x_diff[2] * xp->x_diff[2];
h_max = max(h_max, p->h);
dx_max = max(dx_max, sqrt(dx2));
}
if (c->split) {
for (int k = 0; k < 8; k++) {
if (c->progeny[k] != NULL) {
struct cell *cp = c->progeny[k];
checkCellhdxmax(cp, depth);
dx_max = max(dx_max, cp->dx_max);
h_max = max(h_max, cp->h_max);
}
}
}
/* Check. */
int result = 1;
if (c->h_max != h_max) {
message("%d Inconsistent h_max: cell %f != parts %f", *depth, c->h_max,
h_max);
error("location: %f %f %f", c->loc[0], c->loc[1], c->loc[2]);
message("location: %f %f %f", c->loc[0], c->loc[1], c->loc[2]);
result = 0;
}
if (c->dx_max != dx_max) {
message("%d Inconsistent dx_max: %f != %f", *depth, c->dx_max, dx_max);
error("location: %f %f %f", c->loc[0], c->loc[1], c->loc[2]);
message("location: %f %f %f", c->loc[0], c->loc[1], c->loc[2]);
result = 0;
}
/* Check rebuild criterion. */
/* if (h_max > c->dmin) {
message("%d Inconsistent c->dmin: %f > %f", *depth, h_max, c->dmin);
error("location: %f %f %f", c->loc[0], c->loc[1], c->loc[2]);
result = 0;
} */
return result;
}
......
src/engine.c
View file @ c8be4935
......@@ -2291,15 +2291,16 @@ void engine_rebuild(struct engine *e) {
* @param drift_all Whether to drift particles before rebuilding or not. Will
* not be necessary if all particles have already been
* drifted (before repartitioning for instance).
* @param deferskip Whether to defer the skip until after the rebuild.
* Needed after a repartition.
* @param postrepart If we have just repartitioned, if so we need to defer the
* skip until after the rebuild and not check the if all
* cells have been drifted.
*/
void engine_prepare(struct engine *e, int drift_all, int deferskip) {
void engine_prepare(struct engine *e, int drift_all, int postrepart) {
TIMER_TIC;
/* Unskip active tasks and check for rebuild */
if (!deferskip) engine_unskip(e);
if (!postrepart) engine_unskip(e);
/* Run through the tasks and mark as skip or not. */
int rebuild = e->forcerebuild;
......@@ -2320,13 +2321,15 @@ void engine_prepare(struct engine *e, int drift_all, int deferskip) {
if (drift_all) engine_drift_all(e);
#ifdef SWIFT_DEBUG_CHECKS
/* Check that all cells have been drifted to the current time */
space_check_drift_point(e->s, e->ti_current);
/* Check that all cells have been drifted to the current time, unless
* we have just repartitioned, that can include cells that have not
* previously been active on this rank. */
if (!postrepart) space_check_drift_point(e->s, e->ti_current);
#endif
engine_rebuild(e);
}
if (deferskip) engine_unskip(e);
if (postrepart) engine_unskip(e);
/* Re-rank the tasks every now and then. */
if (e->tasks_age % engine_tasksreweight == 1) {
......@@ -2807,6 +2810,11 @@ void engine_drift_all(struct engine *e) {
threadpool_map(&e->threadpool, runner_do_drift_mapper, e->s->cells_top,
e->s->nr_cells, sizeof(struct cell), 1, e);
#ifdef SWIFT_DEBUG_CHECKS
/* Check that all cells have been drifted to the current time. */
space_check_drift_point(e->s, e->ti_current);
#endif
if (e->verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
clocks_getunit());
......
src/engine.h
View file @ c8be4935
......@@ -225,7 +225,7 @@ void engine_init(struct engine *e, struct space *s,
const struct cooling_function_data *cooling,
struct sourceterms *sourceterms);
void engine_launch(struct engine *e, int nr_runners);
void engine_prepare(struct engine *e, int drift_all, int deferskip);
void engine_prepare(struct engine *e, int drift_all, int postrepart);
void engine_print(struct engine *e);
void engine_init_particles(struct engine *e, int flag_entropy_ICs);
void engine_step(struct engine *e);
......
src/hydro/Default/hydro.h
View file @ c8be4935
......@@ -33,7 +33,7 @@
* @param dt Time since the last kick
*/
__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
const struct part *restrict p, float dt) {
const struct part *restrict p) {
return p->u;
}
......@@ -45,7 +45,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
* @param dt Time since the last kick
*/
__attribute__((always_inline)) INLINE static float hydro_get_pressure(
const struct part *restrict p, float dt) {
const struct part *restrict p) {
return gas_pressure_from_internal_energy(p->rho, p->u);
}
......@@ -57,7 +57,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_pressure(
* @param dt Time since the last kick
*/
__attribute__((always_inline)) INLINE static float hydro_get_entropy(
const struct part *restrict p, float dt) {
const struct part *restrict p) {
return gas_entropy_from_internal_energy(p->rho, p->u);
}
......@@ -69,7 +69,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_entropy(
* @param dt Time since the last kick
*/
__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
const struct part *restrict p, float dt) {
const struct part *restrict p) {
return p->force.soundspeed;
}
......@@ -97,34 +97,30 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
}
/**
* @brief Modifies the thermal state of a particle to the imposed internal
* energy
* @brief Returns the time derivative of internal energy of a particle
*
* This overrides the current state of the particle but does *not* change its
* time-derivatives
* We assume a constant density.
*
* @param p The particle
* @param u The new internal energy
* @param p The particle of interest
*/
__attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
struct part *restrict p, float u) {
__attribute__((always_inline)) INLINE static float hydro_get_internal_energy_dt(
const struct part *restrict p) {
p->u = u;
return p->force.u_dt;
}
/**
* @brief Modifies the thermal state of a particle to the imposed entropy
* @brief Returns the time derivative of internal energy of a particle
*
* This overrides the current state of the particle but does *not* change its
* time-derivatives
* We assume a constant density.
*
* @param p The particle
* @param S The new entropy
* @param p The particle of interest.
* @param du_dt The new time derivative of the internal energy.
*/
__attribute__((always_inline)) INLINE static void hydro_set_entropy(
struct part *restrict p, float S) {
__attribute__((always_inline)) INLINE static void hydro_set_internal_energy_dt(
struct part *restrict p, float du_dt) {
p->u = gas_internal_energy_from_entropy(p->rho, S);
p->force.u_dt = du_dt;
}
/**
......@@ -152,26 +148,6 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
return min(dt_cfl, dt_u_change);
}
/**
* @brief Initialises the particles for the first time
*
* This function is called only once just after the ICs have been
* read in to do some conversions.
*
* @param p The particle to act upon
* @param xp The extended particle data to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_first_init_part(
struct part *restrict p, struct xpart *restrict xp) {
p->ti_begin = 0;
p->ti_end = 0;
xp->v_full[0] = p->v[0];
xp->v_full[1] = p->v[1];
xp->v_full[2] = p->v[2];
xp->u_full = p->u;
}
/**
* @brief Prepares a particle for the density calculation.
*
......@@ -244,8 +220,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
* @param time The current time
*/
__attribute__((always_inline)) INLINE static void hydro_prepare_force(
struct part *restrict p, struct xpart *restrict xp, int ti_current,
double timeBase) {
struct part *restrict p, struct xpart *restrict xp) {
/* Some smoothing length multiples. */
const float h = p->h;
......@@ -270,17 +245,18 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
p->force.balsara = normDiv_v / (normDiv_v + normRot_v + 0.0001f * fc * h_inv);
/* Viscosity parameter decay time */
const float tau = h / (2.f * const_viscosity_length * p->force.soundspeed);
/* const float tau = h / (2.f * const_viscosity_length * p->force.soundspeed);
*/
/* Viscosity source term */
const float S = max(-normDiv_v, 0.f);
/* const float S = max(-normDiv_v, 0.f); */
/* Compute the particle's viscosity parameter time derivative */
const float alpha_dot = (const_viscosity_alpha_min - p->alpha) / tau +
(const_viscosity_alpha_max - p->alpha) * S;
/* const float alpha_dot = (const_viscosity_alpha_min - p->alpha) / tau + */
/* (const_viscosity_alpha_max - p->alpha) * S; */
/* Update particle's viscosity paramter */
p->alpha += alpha_dot * (p->ti_end - p->ti_begin) * timeBase;
/* p->alpha += alpha_dot * (p->ti_end - p->ti_begin) * timeBase; */ // MATTHIEU
}
/**
......@@ -305,6 +281,22 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
p->force.v_sig = 0.0f;
}
/**
* @brief Sets the values to be predicted in the drifts to their values at a
* kick time
*
* @param p The particle.
* @param xp The extended data of this particle.
*/
__attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
struct part *restrict p, const struct xpart *restrict xp) {
/* Re-set the predicted velocities */
p->v[0] = xp->v_full[0];
p->v[1] = xp->v_full[1];
p->v[2] = xp->v_full[2];
}
/**
* @brief Predict additional particle fields forward in time when drifting
*
......@@ -316,8 +308,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
* @param timeBase The minimal time-step size
*/
__attribute__((always_inline)) INLINE static void hydro_predict_extra(
struct part *restrict p, struct xpart *restrict xp, float dt, int t0,
int t1, double timeBase) {
struct part *restrict p, struct xpart *restrict xp, float dt) {
float u, w;
const float h_inv = 1.f / p->h;
......@@ -368,8 +359,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
* @param half_dt The half time-step for this kick
*/
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
struct part *restrict p, struct xpart *restrict xp, float dt,
float half_dt) {}
struct part *restrict p, struct xpart *restrict xp, float dt) {}
/**
* @brief Converts hydro quantity of a particle at the start of a run
......@@ -379,6 +369,28 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
* @param p The particle to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part *restrict p) {}
struct part *restrict p, struct xpart *restrict xp) {}
/**
* @brief Initialises the particles for the first time
*
* This function is called only once just after the ICs have been
* read in to do some conversions.
*
* @param p The particle to act upon
* @param xp The extended particle data to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_first_init_part(
struct part *restrict p, struct xpart *restrict xp) {
p->time_bin = 0;
xp->v_full[0] = p->v[0];
xp->v_full[1] = p->v[1];
xp->v_full[2] = p->v[2];
xp->u_full = p->u;
hydro_reset_acceleration(p);
hydro_init_part(p);
}
#endif /* SWIFT_DEFAULT_HYDRO_H */
src/hydro/Default/hydro_debug.h
View file @ c8be4935
......@@ -25,11 +25,10 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
"x=[%.3e,%.3e,%.3e], "
"v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
"h=%.3e, "
"wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, t_begin=%d, t_end=%d\n",
"wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, time_bin=%d\n",
p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
p->h, (int)p->density.wcount, p->mass, p->rho_dh, p->rho, p->ti_begin,
p->ti_end);
p->h, (int)p->density.wcount, p->mass, p->rho_dh, p->rho, p->time_bin);
}
#endif /* SWIFT_DEFAULT_HYDRO_DEBUG_H */
src/hydro/Default/hydro_part.h
View file @ c8be4935
......@@ -55,12 +55,6 @@ struct part {
/* Particle cutoff radius. */
float h;
/* Particle time of beginning of time-step. */
int ti_begin;
/* Particle time of end of time-step. */
int ti_end;
/* Particle internal energy. */
float u;
......@@ -125,6 +119,9 @@ struct part {
/* Pointer to corresponding gravity part. */
struct gpart* gpart;
/* Particle time-bin */
timebin_t time_bin;
} SWIFT_STRUCT_ALIGN;
#endif /* SWIFT_DEFAULT_HYDRO_PART_H */
src/runner.c
View file @ c8be4935
......@@ -734,8 +734,16 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
}
}
#ifdef SWIFT_DEBUG_CHECKS
if (count) {
message("Smoothing length failed to converge on %i particles.", count);
error("Aborting....");
}
#else
if (count)
message("Smoothing length failed to converge on %i particles.", count);
#endif
/* Be clean */
free(pid);
......@@ -1319,7 +1327,13 @@ void *runner_main(void *data) {
#ifdef SWIFT_DEBUG_CHECKS
t->ti_run = e->ti_current;
#ifndef WITH_MPI
if (cj == NULL) { /* self */
if (ci == NULL && cj == NULL) {
if (t->type != task_type_grav_gather_m && t->type != task_type_grav_fft)
error("Task not associated with cells!");
} else if (cj == NULL) { /* self */
if (!cell_is_active(ci, e) && t->type != task_type_sort &&
t->type != task_type_send && t->type != task_type_recv)
error(
......
src/scheduler.c
View file @ c8be4935
......@@ -1056,8 +1056,6 @@ void scheduler_enqueue_mapper(void *map_data, int num_elements,
*/
void scheduler_start(struct scheduler *s) {
message("launching %i active tasks.", s->active_count);
/* Re-wait the tasks. */
if (s->active_count > 1000) {
threadpool_map(s->threadpool, scheduler_rewait_mapper, s->tid_active,
......@@ -1082,7 +1080,12 @@ void scheduler_start(struct scheduler *s) {
/* Don't check MPI stuff */
if (t->type == task_type_send || t->type == task_type_recv) continue;
if (cj == NULL) { /* self */
if (ci == NULL && cj == NULL) {
if (t->type != task_type_grav_gather_m && t->type != task_type_grav_fft)
error("Task not associated with cells!");
} else if (cj == NULL) { /* self */
if (ci->ti_end_min == ti_current && t->skip &&
t->type != task_type_sort && t->type)
......
src/timestep.h
View file @ c8be4935
......@@ -70,14 +70,15 @@ make_integer_timestep(float new_dt, timebin_t old_bin, integertime_t ti_current,
__attribute__((always_inline)) INLINE static integertime_t get_gpart_timestep(
const struct gpart *restrict gp, const struct engine *restrict e) {
const float new_dt_external = external_gravity_timestep(
e->time, e->external_potential, e->physical_constants, gp);
float new_dt = FLT_MAX;
/* const float new_dt_self = */
/* gravity_compute_timestep_self(e->physical_constants, gp); */
const float new_dt_self = FLT_MAX; // MATTHIEU
if (e->policy & engine_policy_external_gravity)
new_dt =
min(new_dt, external_gravity_timestep(e->time, e->external_potential,
e->physical_constants, gp));
float new_dt = min(new_dt_external, new_dt_self);
if (e->policy & engine_policy_self_gravity)
new_dt = min(new_dt, gravity_compute_timestep_self(gp));
/* Limit timestep within the allowed range */
new_dt = min(new_dt, e->dt_max);
......@@ -114,14 +115,13 @@ __attribute__((always_inline)) INLINE static integertime_t get_part_timestep(
float new_dt_grav = FLT_MAX;
if (p->gpart != NULL) {
const float new_dt_external = external_gravity_timestep(
e->time, e->external_potential, e->physical_constants, p->gpart);
if (e->policy & engine_policy_external_gravity)
new_dt_grav = min(new_dt_grav, external_gravity_timestep(
e->time, e->external_potential,
e->physical_constants, p->gpart));
/* const float new_dt_self = */
/* gravity_compute_timestep_self(e->physical_constants, p->gpart); */
const float new_dt_self = FLT_MAX; // MATTHIEU
new_dt_grav = min(new_dt_external, new_dt_self);
if (e->policy & engine_policy_self_gravity)
new_dt_grav = min(new_dt_grav, gravity_compute_timestep_self(p->gpart));
}
/* Final time-step is minimum of hydro and gravity */
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment