diff --git a/examples/ExternalPointMass/externalPointMass.yml b/examples/ExternalPointMass/externalPointMass.yml index 2712cde82d43fb7d2babaf330bca37ce4ebb02a3..df68307fa36cd01488e0276630003fd822b294c3 100644 --- a/examples/ExternalPointMass/externalPointMass.yml +++ b/examples/ExternalPointMass/externalPointMass.yml @@ -6,9 +6,6 @@ InternalUnitSystem: UnitCurrent_in_cgs: 1 # Amperes UnitTemp_in_cgs: 1 # Kelvin -#Scheduler: -# max_top_level_cells: 20 - # Parameters governing the time integration TimeIntegration: time_begin: 0. # The starting time of the simulation (in internal units). diff --git a/examples/SedovBlast_3D/sedov.yml b/examples/SedovBlast_3D/sedov.yml index 1cc4aced9af3314cadde44f768016225426addf6..c7ea6bfb608d61b631af93c48804202fc0d8f50b 100644 --- a/examples/SedovBlast_3D/sedov.yml +++ b/examples/SedovBlast_3D/sedov.yml @@ -6,6 +6,9 @@ InternalUnitSystem: UnitCurrent_in_cgs: 1 # Amperes UnitTemp_in_cgs: 1 # Kelvin +#Scheduler: +# max_top_level_cells: 3 + # Parameters governing the time integration TimeIntegration: time_begin: 0. # The starting time of the simulation (in internal units). diff --git a/src/active.h b/src/active.h index 3c679fca7e85a0d8562b04da42d941a637ee77e8..7769c3adbdeff2026585a99c15a737f1a9212603 100644 --- a/src/active.h +++ b/src/active.h @@ -143,7 +143,6 @@ __attribute__((always_inline)) INLINE static int gpart_is_active( return (ti_end == ti_current); } - /* Are cells / particles active for kick1 tasks ? */ /** @@ -179,7 +178,8 @@ __attribute__((always_inline)) INLINE static int part_is_starting( const struct part *p, const struct engine *e) { const integertime_t ti_current = e->ti_current; - const integertime_t ti_beg = get_integer_time_begin(ti_current+1, p->time_bin); + const integertime_t ti_beg = + get_integer_time_begin(ti_current + 1, p->time_bin); #ifdef SWIFT_DEBUG_CHECKS if (ti_beg > ti_current) @@ -203,7 +203,8 @@ __attribute__((always_inline)) INLINE static int gpart_is_starting( const struct gpart *gp, const struct engine *e) { const integertime_t ti_current = e->ti_current; - const integertime_t ti_beg = get_integer_time_begin(ti_current+1, gp->time_bin); + const integertime_t ti_beg = + get_integer_time_begin(ti_current + 1, gp->time_bin); #ifdef SWIFT_DEBUG_CHECKS if (ti_beg > ti_current) @@ -216,5 +217,4 @@ __attribute__((always_inline)) INLINE static int gpart_is_starting( return (ti_beg == ti_current); } - #endif /* SWIFT_ACTIVE_H */ diff --git a/src/cell.c b/src/cell.c index 833961f574c58a5dd1bfd7cd3b8425d330485750..a8e23dee838f2bf00493ea876a2e5271f8ed3b0b 100644 --- a/src/cell.c +++ b/src/cell.c @@ -912,6 +912,11 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) { struct engine *e = s->space->e; #endif + //message("unskip! c=%p c->kick1=%p c->drift=%p", c, c->kick1, c->drift); + //task_print(c->drift); + + int ret = 0; + /* Un-skip the density tasks involved with this cell. */ for (struct link *l = c->density; l != NULL; l = l->next) { struct task *t = l->t; @@ -937,7 +942,7 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) { (max(ci->h_max, cj->h_max) + ci->dx_max + cj->dx_max > cj->dmin || ci->dx_max > space_maxreldx * ci->h_max || cj->dx_max > space_maxreldx * cj->h_max)) - return 1; + ret = 1; #ifdef WITH_MPI /* Activate the send/recv flags. */ @@ -1034,7 +1039,9 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) { if (c->cooling != NULL) scheduler_activate(s, c->cooling); if (c->sourceterms != NULL) scheduler_activate(s, c->sourceterms); - return 0; + //task_print(c->drift); + + return ret; } /** @@ -1077,7 +1084,8 @@ void cell_drift(struct cell *c, const struct engine *e) { float dx_max = 0.f, dx2_max = 0.f, h_max = 0.f; /* if (c->gcount > 0) */ - /* message("dt=%e, ti_old=%lld ti_current=%lld", dt, c->ti_old, e->ti_current); */ + /* message("dt=%e, ti_old=%lld ti_current=%lld", dt, c->ti_old, + * e->ti_current); */ /* Check that we are actually going to move forward. */ if (ti_current < ti_old) error("Attempt to drift to the past"); diff --git a/src/cell.h b/src/cell.h index b72f94570afd7c16f355d3b3220cdfa74561364a..63bd4da352e99ca5ceeae962e238977fc18b9ca1 100644 --- a/src/cell.h +++ b/src/cell.h @@ -226,7 +226,7 @@ struct cell { /*! Maximum beginning of (integer) time step in this cell. */ integertime_t ti_beg_max; - + /*! Last (integer) time the cell's content was drifted forward in time. */ integertime_t ti_old; diff --git a/src/drift.h b/src/drift.h index 98983cb98b59cccdd72597697813d1fd6ba1f234..9eb86b86973719eb78baf8209bbeed658ae72e7d 100644 --- a/src/drift.h +++ b/src/drift.h @@ -53,8 +53,8 @@ __attribute__((always_inline)) INLINE static void drift_gpart( gp->ti_drift = ti_current; #endif - //message("dt= %e", dt); - //fprintf(files_timestep[gp->id_or_neg_offset], "drift: dt=%e\n", dt); + // message("dt= %e", dt); + // fprintf(files_timestep[gp->id_or_neg_offset], "drift: dt=%e\n", dt); /* Drift... */ gp->x[0] += gp->v_full[0] * dt; diff --git a/src/engine.c b/src/engine.c index 5f74aea327785a05641d71e923719d3c2efc60f4..da9ea5059dfb28dccd142f7812602e5a33a36018 100644 --- a/src/engine.c +++ b/src/engine.c @@ -71,7 +71,7 @@ #include "units.h" #include "version.h" -FILE* files_timestep[num_files]; +FILE *files_timestep[num_files]; /* Particle cache size. */ #define CACHE_SIZE 512 @@ -2162,10 +2162,11 @@ void engine_marktasks_mapper(void *map_data, int num_elements, } /* Kick/Drift/Init? */ - else if (t->type == task_type_kick1 || t->type == task_type_kick2 || t->type == task_type_drift || t->type == task_type_init) { + else if (t->type == task_type_kick1 || t->type == task_type_kick2 || + t->type == task_type_drift || t->type == task_type_init) { if (cell_is_active(t->ci, e)) scheduler_activate(s, t); } - + /* Time-step? */ else if (t->type == task_type_timestep) { t->ci->updated = 0; @@ -2192,6 +2193,8 @@ int engine_marktasks(struct engine *e) { const ticks tic = getticks(); int rebuild_space = 0; + message("marktasks"); + /* Run through the tasks and mark as skip or not. */ size_t extra_data[3] = {(size_t)e, rebuild_space, (size_t)&e->sched}; threadpool_map(&e->threadpool, engine_marktasks_mapper, s->tasks, s->nr_tasks, @@ -2280,7 +2283,7 @@ void engine_rebuild(struct engine *e) { error("engine_marktasks failed after space_rebuild."); /* Print the status of the system */ - //if (e->verbose) engine_print_task_counts(e); + // if (e->verbose) engine_print_task_counts(e); if (e->verbose) message("took %.3f %s.", clocks_from_ticks(getticks() - tic), @@ -2318,60 +2321,6 @@ void engine_prepare(struct engine *e) { clocks_from_ticks(getticks() - tic), clocks_getunit()); } -/* void engine_prepare(struct engine *e, int drift_all, int postrepart) { */ - -/* TIMER_TIC; */ - -/* // message("e->forcerebuild: %d", e->forcerebuild); */ - -/* /\* Unskip active tasks and check for rebuild *\/ */ -/* if (!postrepart) engine_unskip(e); */ - -/* /\* Run through the tasks and mark as skip or not. *\/ */ -/* int rebuild = e->forcerebuild; */ - -/* // message("e->forcerebuild: %d", e->forcerebuild); */ - -/* /\* Collect the values of rebuild from all nodes. *\/ */ -/* #ifdef WITH_MPI */ -/* int buff = 0; */ -/* if (MPI_Allreduce(&rebuild, &buff, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD) != - */ -/* MPI_SUCCESS) */ -/* error("Failed to aggregate the rebuild flag across nodes."); */ -/* rebuild = buff; */ -/* #endif */ - -/* /\* And rebuild if necessary. *\/ */ -/* if (rebuild) { */ - -/* /\* Drift all particles to the current time if needed. *\/ */ -/* if (drift_all) engine_drift_all(e); */ - -/* #ifdef SWIFT_DEBUG_CHECKS */ -/* /\* 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 (postrepart) engine_unskip(e); */ - -/* /\* Re-rank the tasks every now and then. *\/ */ -/* if (e->tasks_age % engine_tasksreweight == 1) { */ -/* scheduler_reweight(&e->sched, e->verbose); */ -/* } */ -/* e->tasks_age += 1; */ - -/* TIMER_TOC(timer_prepare); */ - -/* if (e->verbose) */ -/* message("took %.3f %s (including drift all, rebuild and reweight).", */ -/* clocks_from_ticks(getticks() - tic), clocks_getunit()); */ -/* } */ - /** * @brief Implements a barrier for the #runner threads. * @@ -2468,7 +2417,7 @@ void engine_collect_timestep(struct engine *e) { const ticks tic = getticks(); int updates = 0, g_updates = 0; - integertime_t ti_end_min = max_nr_timesteps; + integertime_t ti_end_min = max_nr_timesteps, ti_end_max = 0, ti_beg_max = 0; const struct space *s = e->s; /* Collect the cell data. */ @@ -2481,6 +2430,8 @@ void engine_collect_timestep(struct engine *e) { /* And aggregate */ ti_end_min = min(ti_end_min, c->ti_end_min); + ti_end_max = max(ti_end_max, c->ti_end_max); + ti_beg_max = max(ti_beg_max, c->ti_beg_max); updates += c->updated; g_updates += c->g_updated; @@ -2514,6 +2465,8 @@ void engine_collect_timestep(struct engine *e) { #endif e->ti_end_min = ti_end_min; + e->ti_end_max = ti_end_max; + e->ti_beg_max = ti_beg_max; e->updates = updates; e->g_updates = g_updates; @@ -2668,7 +2621,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) { /* Print the number of active tasks ? */ if (e->verbose) engine_print_task_counts(e); - + /* Now, launch the calculation */ TIMER_TIC; engine_launch(e, e->nr_threads); @@ -2701,7 +2654,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) { /* Print the number of active tasks ? */ if (e->verbose) engine_print_task_counts(e); - + engine_launch(e, e->nr_threads); /* Recover the (integer) end of the next time-step */ @@ -2715,7 +2668,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) { /* Ready to go */ e->step = 0; - e->forcerebuild = 0; + e->forcerebuild = 1; e->wallclock_time = (float)clocks_diff(&time1, &time2); if (e->verbose) message("took %.3f %s.", e->wallclock_time, clocks_getunit()); @@ -2746,8 +2699,8 @@ void engine_step(struct engine *e) { if (e->nodeID == 0) { /* Print some information to the screen */ - printf(" %6d %14e %14e %10zu %10zu %21.3f\n", e->step, e->time, - e->timeStep, e->updates, e->g_updates, e->wallclock_time); + printf(" %6d %lld %14e %14e %10zu %10zu %21.3f\n", e->step, e->ti_current, + e->time, e->timeStep, e->updates, e->g_updates, e->wallclock_time); fflush(stdout); fprintf(e->file_timesteps, " %6d %14e %14e %10zu %10zu %21.3f\n", e->step, @@ -2755,24 +2708,26 @@ void engine_step(struct engine *e) { fflush(e->file_timesteps); /* for(int i=0; i<num_files; ++i) { */ - /* fprintf(files_timestep[i], " %6d %14e %14e %10zu %10zu %21.3f\n", e->step, */ - /* e->time, e->timeStep, e->updates, e->g_updates, e->wallclock_time); */ + /* fprintf(files_timestep[i], " %6d %14e %14e %10zu %10zu %21.3f\n", + * e->step, */ + /* e->time, e->timeStep, e->updates, e->g_updates, + * e->wallclock_time); */ /* fflush(files_timestep[i]); */ /* } */ } /* Do we need repartitioning ? */ - if(e->forcerepart) engine_repartition(e); + if (e->forcerepart != REPART_NONE) engine_repartition(e); /* Do we need rebuilding ? */ - if(e->forcerebuild) engine_rebuild(e); - + if (e->forcerebuild) engine_rebuild(e); + /* Prepare the tasks to be launched. */ engine_prepare(e); - if (e->step % 99 == 0) e->forcerebuild = 1; - + if (e->step == 1) e->forcerebuild = 1; + /* Print the number of active tasks ? */ if (e->verbose) engine_print_task_counts(e); @@ -2792,8 +2747,7 @@ void engine_step(struct engine *e) { /* Drift everybody (i.e. what has not yet been drifted) */ /* to the current time */ - if(e->dump_snapshot || e->forcerebuild || e->forcerepart != REPART_NONE) { - message("snap=%d, rebuild=%d repart=%d", e->dump_snapshot, e->forcerebuild, e->forcerepart); + if (e->dump_snapshot || e->forcerebuild || e->forcerepart != REPART_NONE) { engine_drift_all(e); @@ -2803,11 +2757,14 @@ void engine_step(struct engine *e) { * previously been active on this rank. */ space_check_drift_point(e->s, e->ti_current); #endif - } + message("snap=%d, rebuild=%d repart=%d", e->dump_snapshot, e->forcerebuild, + e->forcerepart); + + /* Write a snapshot ? */ - if(e->dump_snapshot) { + if (e->dump_snapshot) { /* Dump... */ engine_dump_snapshot(e); @@ -2815,49 +2772,11 @@ void engine_step(struct engine *e) { /* ... and find the next output time */ engine_compute_next_snapshot_time(e); } - + /* Recover the (integer) end of the next time-step */ engine_collect_timestep(e); - -/* /\* Check for snapshot dump *\/ */ -/* while { */ - -/* /\* Drift everybody (i.e. what has not yet been drifted) *\/ */ -/* /\* to the current time *\/ */ -/* engine_drift_all(e); */ - -/* /\* Dump... *\/ */ -/* engine_dump_snapshot(e); */ - -/* /\* ... and find the next output time *\/ */ -/* engine_compute_next_snapshot_time(e); */ -/* } */ - -/* /\* Repartition if necessary *\/ */ -/* if (e->forcerepart != REPART_NONE) { */ -/* /\* Drift all particles to the current time if needed. *\/ */ -/* engine_drift_all(e); */ - -/* engine_repartition(e); */ -/* } */ - -/* /\* Rebuild if necessary *\/ */ -/* if (e->forcerebuild) { */ - -/* /\* Drift all particles to the current time if needed. *\/ */ -/* engine_drift_all(e); */ - -/* #ifdef SWIFT_DEBUG_CHECKS */ -/* /\* Check that all cells have been drifted to the current time. */ -/* * That can include cells that have not */ -/* * previously been active on this rank. *\/ */ -/* space_check_drift_point(e->s, e->ti_current); */ -/* #endif */ - -/* engine_rebuild(e); */ -/* } */ TIMER_TOC2(timer_step); @@ -2867,91 +2786,6 @@ void engine_step(struct engine *e) { e->toc_step = getticks(); } -/* void engine_step(struct engine *e) { */ - -/* TIMER_TIC2; */ - -/* struct clocks_time time1, time2; */ -/* clocks_gettime(&time1); */ - -/* e->tic_step = getticks(); */ - -/* /\* Recover the (integer) end of the next time-step *\/ */ -/* engine_collect_timestep(e); */ - -/* /\* Move forward in time *\/ */ -/* e->ti_old = e->ti_current; */ -/* e->ti_current = e->ti_end_min; */ -/* e->step += 1; */ -/* e->time = e->ti_current * e->timeBase + e->timeBegin; */ -/* e->timeOld = e->ti_old * e->timeBase + e->timeBegin; */ -/* e->timeStep = (e->ti_current - e->ti_old) * e->timeBase; */ -/* // if(e->step % 100 == 0) e->forcerebuild = 1; */ - -/* if (e->nodeID == 0) { */ - -/* /\* Print some information to the screen *\/ */ -/* printf(" %6d %14e %14e %10zu %10zu %21.3f\n", e->step, e->time, */ -/* e->timeStep, e->updates, e->g_updates, e->wallclock_time); */ -/* fflush(stdout); */ - -/* fprintf(e->file_timesteps, " %6d %14e %14e %10zu %10zu %21.3f\n", - * e->step, */ -/* e->time, e->timeStep, e->updates, e->g_updates, - * e->wallclock_time); */ -/* fflush(e->file_timesteps); */ -/* } */ - -/* /\* Drift only the necessary particles, that means all particles */ -/* * if we are about to repartition. *\/ */ -/* const int repart = (e->forcerepart != REPART_NONE); */ -/* // message("repart: %d", repart); */ -/* const int drift_all = (e->policy & engine_policy_drift_all); */ -/* // message("drift_all: %d", drift_all); */ -/* if (repart || drift_all) engine_drift_all(e); */ - -/* /\* Re-distribute the particles amongst the nodes? *\/ */ -/* if (repart) engine_repartition(e); */ - -/* /\* Prepare the space. *\/ */ -/* engine_prepare(e, !(drift_all || repart), repart); */ - -/* if (e->verbose) engine_print_task_counts(e); */ - -/* /\* Save some statistics *\/ */ -/* if (e->time - e->timeLastStatistics >= e->deltaTimeStatistics) { */ -/* engine_print_stats(e); */ -/* e->timeLastStatistics += e->deltaTimeStatistics; */ -/* } */ - -/* /\* Send off the runners. *\/ */ -/* TIMER_TIC; */ -/* engine_launch(e, e->nr_threads); */ -/* TIMER_TOC(timer_runners); */ - -/* /\* Check for output *\/ */ -/* while (e->ti_end_min >= e->ti_nextSnapshot && e->ti_nextSnapshot > 0) { */ - -/* /\* Drift everybody (i.e. what has not yet been drifted) to the snapshot - */ -/* * position *\/ */ -/* engine_drift_all(e); */ - -/* /\* Dump... *\/ */ -/* engine_dump_snapshot(e); */ - -/* /\* ... and find the next output time *\/ */ -/* engine_compute_next_snapshot_time(e); */ -/* } */ - -/* TIMER_TOC2(timer_step); */ - -/* clocks_gettime(&time2); */ - -/* e->wallclock_time = (float)clocks_diff(&time1, &time2); */ -/* e->toc_step = getticks(); */ -/* } */ - /** * @brief Returns 1 if the simulation has reached its end point, 0 otherwise */ @@ -2966,6 +2800,8 @@ int engine_is_done(struct engine *e) { */ void engine_unskip(struct engine *e) { + message("unskip"); + const ticks tic = getticks(); threadpool_map(&e->threadpool, runner_do_unskip_mapper, e->s->cells_top, e->s->nr_cells, sizeof(struct cell), 1, e); @@ -2988,7 +2824,6 @@ void engine_drift_all(struct engine *e) { /* fprintf(files_timestep[i], "drift all\n"); */ /* } */ - const ticks tic = getticks(); threadpool_map(&e->threadpool, runner_do_drift_mapper, e->s->cells_top, e->s->nr_cells, sizeof(struct cell), 1, e); @@ -3210,7 +3045,6 @@ void engine_dump_snapshot(struct engine *e) { /* fprintf(files_timestep[i], "dump\n"); */ /* } */ - message("dump"); /* Dump... */ @@ -3230,7 +3064,7 @@ void engine_dump_snapshot(struct engine *e) { #endif e->dump_snapshot = 0; - + clocks_gettime(&time2); if (e->verbose) message("writing particle properties took %.3f %s.", @@ -3714,12 +3548,12 @@ void engine_init(struct engine *e, struct space *s, free(buf); #endif - for(int i=0; i<num_files; ++i) { + for (int i = 0; i < num_files; ++i) { char name[10]; sprintf(name, "dt_%d.txt", i); files_timestep[i] = fopen(name, "w"); } - + /* Wait for the runner threads to be in place. */ while (e->barrier_running || e->barrier_launch) if (pthread_cond_wait(&e->barrier_cond, &e->barrier_mutex) != 0) @@ -3787,9 +3621,8 @@ void engine_compute_next_snapshot_time(struct engine *e) { */ void engine_clean(struct engine *e) { - for(int i=0; i<num_files; ++i) - fclose(files_timestep[i]); - + for (int i = 0; i < num_files; ++i) fclose(files_timestep[i]); + free(e->snapshotUnits); free(e->links); scheduler_clean(&e->sched); diff --git a/src/engine.h b/src/engine.h index 414cc2966b1081280548124f3fda426b8aa31c1e..2c13f8040ea6e461b88a0c245ea89ab148024864 100644 --- a/src/engine.h +++ b/src/engine.h @@ -130,6 +130,12 @@ struct engine { /* Minimal ti_end for the next time-step */ integertime_t ti_end_min; + /* Maximal ti_end for the next time-step */ + integertime_t ti_end_max; + + /* Maximal ti_beg for the next time-step */ + integertime_t ti_beg_max; + /* Number of particles updated */ size_t updates, g_updates; @@ -182,7 +188,7 @@ struct engine { /* Need to dump a snapshot ? */ int dump_snapshot; - + /* How many steps have we done with the same set of tasks? */ int tasks_age; diff --git a/src/kick.h b/src/kick.h index fddac41efd17bf585eecb78630fd9ba801431a63..92f25461ad9743ad17f3a107478a930b24e8bc57 100644 --- a/src/kick.h +++ b/src/kick.h @@ -48,13 +48,12 @@ __attribute__((always_inline)) INLINE static void kick_gpart( "Particle has not been kicked to the current time gp->ti_kick=%lld, " "ti_start=%lld, ti_end=%lld", gp->ti_kick, ti_start, ti_end); - + gp->ti_kick = ti_end; #endif - - //message("dt= %e gp->ti_kick=%lld", dt, gp->ti_kick); - //fprintf(files_timestep[gp->id_or_neg_offset], "kick: dt=%e\n", dt); + // message("dt= %e gp->ti_kick=%lld", dt, gp->ti_kick); + // fprintf(files_timestep[gp->id_or_neg_offset], "kick: dt=%e\n", dt); /* Kick particles in momentum space */ gp->v_full[0] += gp->a_grav[0] * dt; diff --git a/src/potential/point_mass/potential.h b/src/potential/point_mass/potential.h index c05838b9f05e757b4ed9996441e28e42cea35a3e..6ccee4d8a7bfcd3212144a4a659e1dccc3e26ab3 100644 --- a/src/potential/point_mass/potential.h +++ b/src/potential/point_mass/potential.h @@ -80,14 +80,17 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep( /* (-g->v_full[1] + 3.f * rinv2 * drdv * dy); */ /* const float dota_z = G_newton * potential->mass * rinv3 * */ /* (-g->v_full[2] + 3.f * rinv2 * drdv * dz); */ - /* const float dota_2 = dota_x * dota_x + dota_y * dota_y + dota_z * dota_z; */ - /* const float a_2 = g->a_grav[0] * g->a_grav[0] + g->a_grav[1] * g->a_grav[1] */ + /* const float dota_2 = dota_x * dota_x + dota_y * dota_y + dota_z * dota_z; + */ + /* const float a_2 = g->a_grav[0] * g->a_grav[0] + g->a_grav[1] * g->a_grav[1] + */ /* + g->a_grav[2] * g->a_grav[2]; */ /* return potential->timestep_mult * sqrtf(a_2 / dota_2); */ const double r = sqrt(dx * dx + dy * dy + dz * dz); - return potential->timestep_mult * sqrt(r * r * r / (G_newton * potential->mass)); + return potential->timestep_mult * + sqrt(r * r * r / (G_newton * potential->mass)); } /** diff --git a/src/runner.c b/src/runner.c index 5b37b6ddba8b233a94b3b352f192225d5fd4de92..41b951a8faaa85af80370ce6ac4261c64d8bfc31 100644 --- a/src/runner.c +++ b/src/runner.c @@ -783,6 +783,8 @@ static void runner_do_unskip(struct cell *c, struct engine *e) { if (forcerebuild) atomic_inc(&e->forcerebuild); } + // message("c->depth=%d c->split=%d c->count=%.5d c->ti_end_min=%lld c->ti_end_max=%lld c->ti_beg_max=%lld ti_current=%lld", c->depth, c->split, c->count, c->ti_end_min, c->ti_end_max, c->ti_beg_max ,e->ti_current); + /* Recurse */ if (c->split) { for (int k = 0; k < 8; k++) { @@ -869,7 +871,7 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) { /* Anything to do here? */ if (!cell_is_starting(c, e)) return; - + /* Recurse? */ if (c->split) { for (int k = 0; k < 8; k++) @@ -890,16 +892,16 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) { const integertime_t ti_begin = get_integer_time_begin(ti_current + 1, p->time_bin); -/* #ifdef SWIFT_DEBUG_CHECKS */ -/* const integertime_t ti_end = */ -/* get_integer_time_end(ti_current, p->time_bin); */ +#ifdef SWIFT_DEBUG_CHECKS + const integertime_t ti_end = + get_integer_time_end(ti_current + 1, p->time_bin); -/* if (ti_end - ti_begin != ti_step) */ -/* error( */ -/* "Particle in wrong time-bin, ti_end=%lld, ti_begin=%lld, " */ -/* "ti_step=%lld time_bin=%d ti_current=%lld", */ -/* ti_end, ti_begin, ti_step, p->time_bin, ti_current); */ -/* #endif */ + if (ti_begin != ti_current) + error( + "Particle in wrong time-bin, ti_end=%lld, ti_begin=%lld, " + "ti_step=%lld time_bin=%d ti_current=%lld", + ti_end, ti_begin, ti_step, p->time_bin, ti_current); +#endif /* do the kick */ kick_part(p, xp, ti_begin, ti_begin + ti_step / 2, timeBase); @@ -919,24 +921,26 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) { const integertime_t ti_begin = get_integer_time_begin(ti_current + 1, gp->time_bin); -/* #ifdef SWIFT_DEBUG_CHECKS */ -/* const integertime_t ti_end = */ -/* get_integer_time_end(ti_current, gp->time_bin); */ +#ifdef SWIFT_DEBUG_CHECKS + const integertime_t ti_end = + get_integer_time_end(ti_current + 1, gp->time_bin); -/* if (ti_end - ti_begin != ti_step) error("Particle in wrong time-bin"); */ -/* #endif */ + if (ti_begin != ti_current) + error( + "Particle in wrong time-bin, ti_end=%lld, ti_begin=%lld, " + "ti_step=%lld time_bin=%d ti_current=%lld", + ti_end, ti_begin, ti_step, gp->time_bin, ti_current); +#endif /* do the kick */ kick_gpart(gp, ti_begin, ti_begin + ti_step / 2, timeBase); - + #ifdef ICHECK if (gp->id_or_neg_offset == ICHECK) { message("--- ti_current=%lld time=%e ---", e->ti_current, e->time); printgParticle_single(gp); } -#endif - } else { - //message("aaa"); +#endif } } } @@ -1097,8 +1101,8 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) { ti_end_min = min(ti_current + ti_new_step, ti_end_min); ti_end_max = max(ti_current + ti_new_step, ti_end_max); - /* What is the next starting point for this cell ? */ - ti_beg_max = max(ti_current, ti_beg_max); + /* What is the next starting point for this cell ? */ + ti_beg_max = max(ti_current, ti_beg_max); } else { /* part is inactive */ @@ -1110,8 +1114,11 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) { ti_end_min = min(ti_end, ti_end_min); ti_end_max = max(ti_end, ti_end_max); - /* What is the next starting point for this cell ? */ - ti_beg_max = max(ti_current, ti_beg_max); + const integertime_t ti_beg = + get_integer_time_begin(ti_current + 1, p->time_bin); + + /* What is the next starting point for this cell ? */ + ti_beg_max = max(ti_beg, ti_beg_max); } } @@ -1156,10 +1163,10 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) { ti_end_min = min(ti_current + ti_new_step, ti_end_min); ti_end_max = max(ti_current + ti_new_step, ti_end_max); - /* What is the next starting point for this cell ? */ - ti_beg_max = max(ti_current, ti_beg_max); + /* What is the next starting point for this cell ? */ + ti_beg_max = max(ti_current, ti_beg_max); - } else { /* gpart is inactive */ + } else { /* gpart is inactive */ const integertime_t ti_end = get_integer_time_end(ti_current, gp->time_bin); @@ -1168,8 +1175,11 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) { ti_end_min = min(ti_end, ti_end_min); ti_end_max = max(ti_end, ti_end_max); - /* What is the next starting point for this cell ? */ - ti_beg_max = max(ti_current, ti_beg_max); + const integertime_t ti_beg = + get_integer_time_begin(ti_current + 1, gp->time_bin); + + /* What is the next starting point for this cell ? */ + ti_beg_max = max(ti_beg, ti_beg_max); } } } @@ -1188,7 +1198,7 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) { g_updated += cp->g_updated; ti_end_min = min(cp->ti_end_min, ti_end_min); ti_end_max = max(cp->ti_end_max, ti_end_max); - ti_beg_max = max(cp->ti_beg_max, ti_beg_max); + ti_beg_max = max(cp->ti_beg_max, ti_beg_max); } } @@ -1456,8 +1466,8 @@ void *runner_main(void *data) { } 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 - && t->type != task_type_kick1) + t->type != task_type_send && t->type != task_type_recv && + t->type != task_type_kick1) error( "Task (type='%s/%s') should have been skipped ti_current=%lld " "c->ti_end_min=%lld", @@ -1481,7 +1491,7 @@ void *runner_main(void *data) { "c->ti_end_min=%lld t->flags=%d", taskID_names[t->type], subtaskID_names[t->subtype], e->ti_current, ci->ti_end_min, t->flags); - + } else { /* pair */ if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) @@ -1580,24 +1590,16 @@ void *runner_main(void *data) { runner_do_extra_ghost(r, ci, 1); break; #endif - case task_type_drift: { - /* integertime_t ti_current = e->ti_current; */ - /* integertime_t ti_old = ci->ti_old; */ - - /* e->ti_current = ci->ti_old + (ti_current - ti_old) / 2; */ - - /* runner_do_drift(r, ci, 1); */ - - /* e->ti_current = ti_current; */ + case task_type_drift: runner_do_drift(r, ci, 1); - } break; + break; case task_type_kick1: runner_do_kick1(r, ci, 1); break; case task_type_kick2: if (!(e->policy & engine_policy_cooling)) runner_do_end_force(r, ci, 1); - if (e->ti_current > 0) runner_do_kick2(r, ci, 1); + runner_do_kick2(r, ci, 1); break; case task_type_timestep: runner_do_timestep(r, ci, 1); diff --git a/src/space.c b/src/space.c index 06514b01e1a55af89889807ed3097b6fa929973d..b289f95b47b830e9192b7cd60dceeb69a24c0053 100644 --- a/src/space.c +++ b/src/space.c @@ -248,7 +248,7 @@ void space_regrid(struct space *s, int verbose) { const size_t nr_parts = s->nr_parts; const ticks tic = getticks(); - //const integertime_t ti_current = (s->e != NULL) ? s->e->ti_current : 0; + // const integertime_t ti_current = (s->e != NULL) ? s->e->ti_current : 0; const integertime_t ti_old = (s->e != NULL) ? s->e->ti_old : 0; /* Run through the cells and get the current h_max. */ @@ -480,14 +480,14 @@ void space_rebuild(struct space *s, int verbose) { /* for(int i=0; i<num_files; ++i) { */ /* fprintf(files_timestep[i], "REBUILD\n"); */ /* } */ - + /* Re-grid if necessary, or just re-set the cell data. */ space_regrid(s, verbose); size_t nr_parts = s->nr_parts; size_t nr_gparts = s->nr_gparts; struct cell *restrict cells_top = s->cells_top; - //const integertime_t ti_current = (s->e != NULL) ? s->e->ti_current : 0; + // const integertime_t ti_current = (s->e != NULL) ? s->e->ti_current : 0; const integertime_t ti_old = (s->e != NULL) ? s->e->ti_old : 0; /* Run through the particles and get their cell index. Allocates @@ -1587,7 +1587,7 @@ void space_split_recursive(struct space *s, struct cell *c, h_max = max(h_max, c->progeny[k]->h_max); ti_end_min = min(ti_end_min, c->progeny[k]->ti_end_min); ti_end_max = max(ti_end_max, c->progeny[k]->ti_end_max); - ti_end_max = max(ti_beg_max, c->progeny[k]->ti_beg_max); + ti_beg_max = max(ti_beg_max, c->progeny[k]->ti_beg_max); if (c->progeny[k]->maxdepth > maxdepth) maxdepth = c->progeny[k]->maxdepth; } @@ -1607,10 +1607,12 @@ void space_split_recursive(struct space *s, struct cell *c, struct part *p = &parts[k]; struct xpart *xp = &xparts[k]; const float h = p->h; + const integertime_t ti_end = + get_integer_time_end(e->ti_current, p->time_bin); /* const integertime_t ti_end = */ - /* e->ti_current + get_integer_timestep(p->time_bin); */ - const integertime_t ti_end = get_integer_time_end(e->ti_current, p->time_bin); - const integertime_t ti_beg = get_integer_time_begin(e->ti_current, p->time_bin); + /* get_integer_time_end(e->ti_current, p->time_bin); */ + const integertime_t ti_beg = + get_integer_time_begin(e->ti_current + 1, p->time_bin); xp->x_diff[0] = 0.f; xp->x_diff[1] = 0.f; xp->x_diff[2] = 0.f; @@ -1621,10 +1623,12 @@ void space_split_recursive(struct space *s, struct cell *c, } for (int k = 0; k < gcount; k++) { struct gpart *gp = &gparts[k]; + const integertime_t ti_end = + get_integer_time_end(e->ti_current, gp->time_bin); /* const integertime_t ti_end = */ - /* e->ti_current + get_integer_timestep(gp->time_bin); */ - const integertime_t ti_end = get_integer_time_end(e->ti_current, gp->time_bin); - const integertime_t ti_beg = get_integer_time_begin(e->ti_current, gp->time_bin); + /* get_integer_time_end(e->ti_current, gp->time_bin); */ + const integertime_t ti_beg = + get_integer_time_begin(e->ti_current + 1, gp->time_bin); gp->x_diff[0] = 0.f; gp->x_diff[1] = 0.f; gp->x_diff[2] = 0.f;