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

Finish the force calculation before triggereing the cooling.

parent 226ff24a
......@@ -9,7 +9,7 @@ InternalUnitSystem:
# Parameters governing the time integration
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
time_end: 0.4 # The end time of the simulation (in internal units).
time_end: 0.25 # The end time of the simulation (in internal units).
dt_min: 1e-5 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
......@@ -17,7 +17,7 @@ TimeIntegration:
Snapshots:
basename: coolingBox # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 1.0e-1 # Time difference between consecutive outputs (in internal units)
delta_time: 1e-2 # Time difference between consecutive outputs (in internal units)
# Parameters governing the conserved quantities statistics
Statistics:
......
......@@ -83,7 +83,7 @@ rho_cgs = rho * unit_mass / (unit_length)**3
time_cgs = time * unit_time
total_energy_cgs = total_energy / total_mass[0] * unit_length**2 / (unit_time)**2
kinetic_energy_cgs = kinetic_energy / total_mass[0] * unit_length**2 / (unit_time)**2
internal_energy_cgs = kinetic_energy / total_mass[0] * unit_length**2 / (unit_time)**2
internal_energy_cgs = internal_energy / total_mass[0] * unit_length**2 / (unit_time)**2
radiated_energy_cgs = radiated_energy / total_mass[0] * unit_length**2 / (unit_time)**2
# Find the energy floor
......@@ -100,19 +100,27 @@ u_analytic_cgs[u_analytic_cgs < u_floor_cgs] = u_floor_cgs
print "Cooling time:", cooling_time_cgs, "[s]"
# Read snapshots
u_snapshots_cgs = zeros(25)
t_snapshots_cgs = zeros(25)
for i in range(25):
snap = h5.File("coolingBox_%0.3d.hdf5"%i,'r')
u_snapshots_cgs[i] = sum(snap["/PartType0/InternalEnergy"][:] * snap["/PartType0/Masses"][:]) / total_mass[0] * unit_length**2 / (unit_time)**2
t_snapshots_cgs[i] = snap["/Header"].attrs["Time"] * unit_time
figure()
#plot(time_cgs, internal_energy_cgs, 'b-', lw=1.5, label="Gas internal energy")
#plot(time_cgs, kinetic_energy_cgs, 'y-', lw=1.5, label="Gas kinetic energy")
plot(time_cgs, total_energy_cgs, 'r-', lw=1.6, label="Gas total energy")
plot(t_snapshots_cgs, u_snapshots_cgs, 'rD', ms=3)
plot(time_cgs, radiated_energy_cgs, 'g-', lw=1.6, label="Radiated energy")
plot(time_cgs, total_energy_cgs + radiated_energy_cgs, 'b-', lw=0.6, label="Gas total + radiated")
plot(analytic_time_cgs, u_analytic_cgs, '--', color='k', alpha=0.8, lw=1.0, label="Analytic solution")
legend(loc="upper right", fontsize=8)
legend(loc="upper right", fontsize=8, frameon=False, handlelength=3, ncol=1)
xlabel("${\\rm{Time~[s]}}$", labelpad=0)
ylabel("${\\rm{Gas~internal~energy~[erg]}}$")
xlim(0, 1.8*cooling_time_cgs)
ylabel("${\\rm{Energy~[erg]}}$")
xlim(0, 1.5*cooling_time_cgs)
ylim(0, 1.5*u_analytic_cgs[0])
savefig("energy.png", dpi=200)
......
......@@ -89,7 +89,7 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
float cooling_du_dt = cooling_rate(phys_const, us, cooling, p);
/* Integrate cooling equation to enforce energy floor */
if (u_old + hydro_du_dt < u_floor) {
if (u_old + hydro_du_dt * dt < u_floor) {
cooling_du_dt = 0.f;
} else if (u_old + (hydro_du_dt + cooling_du_dt) * dt < u_floor) {
cooling_du_dt = (u_old + dt * hydro_du_dt - u_floor) / dt;
......
......@@ -176,13 +176,18 @@ void engine_make_hierarchical_tasks(struct engine *e, struct cell *c) {
#endif
/* Cooling task */
if (is_with_cooling)
if (is_with_cooling) {
c->cooling = scheduler_addtask(s, task_type_cooling, task_subtype_none,
0, 0, c, NULL, 0);
scheduler_addunlock(s, c->cooling, c->kick2);
}
/* add source terms */
if (is_with_sourceterms)
if (is_with_sourceterms) {
c->sourceterms = scheduler_addtask(s, task_type_sourceterms,
task_subtype_none, 0, 0, c, NULL, 0);
}
}
} else { /* We are above the super-cell so need to go deeper */
......@@ -1634,11 +1639,18 @@ static inline void engine_make_hydro_loops_dependencies(struct scheduler *sched,
struct task *density,
struct task *force,
struct cell *c) {
/* init --> density loop --> ghost --> force loop --> kick2 */
/* init --> density loop --> ghost --> force loop */
scheduler_addunlock(sched, c->super->init, density);
scheduler_addunlock(sched, density, c->super->ghost);
scheduler_addunlock(sched, c->super->ghost, force);
scheduler_addunlock(sched, force, c->super->kick2);
if (c->super->cooling != NULL) {
/* force loop --> cooling (--> kick2) */
scheduler_addunlock(sched, force, c->super->cooling);
} else {
/* force loop --> kick2 */
scheduler_addunlock(sched, force, c->super->kick2);
}
}
#endif
......@@ -1832,12 +1844,6 @@ void engine_make_extra_hydroloop_tasks(struct engine *e) {
}
#endif
}
/* Cooling tasks take place after the second kick */
else if (t->type == task_type_cooling) {
scheduler_addunlock(sched, t->ci->kick2, t);
scheduler_addunlock(sched, t, t->ci->timestep);
}
}
}
......@@ -2543,9 +2549,7 @@ void engine_skip_drift_and_kick(struct engine *e) {
struct task *t = &tasks[i];
/* Skip everything that updates the particles */
if (t->type == task_type_drift || t->type == task_type_kick1 ||
t->type == task_type_cooling || t->type == task_type_sourceterms)
t->skip = 1;
if (t->type == task_type_drift || t->type == task_type_kick1) t->skip = 1;
}
}
......
......@@ -924,7 +924,6 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
struct part *restrict parts = c->parts;
struct xpart *restrict xparts = c->xparts;
struct gpart *restrict gparts = c->gparts;
const double const_G = e->physical_constants->const_newton_G;
TIMER_TIC;
......@@ -947,10 +946,6 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
/* If particle needs to be kicked */
if (part_is_active(p, e)) {
/* First, finish the force loop */
hydro_end_force(p);
if (p->gpart != NULL) gravity_end_force(p->gpart, const_G);
const integertime_t ti_step = get_integer_timestep(p->time_bin);
const integertime_t ti_begin =
get_integer_time_begin(ti_current, p->time_bin);
......@@ -977,9 +972,6 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
/* If the g-particle has no counterpart and needs to be kicked */
if (gp->id_or_neg_offset > 0 && gpart_is_active(gp, e)) {
/* First, finish the force loop */
gravity_end_force(gp, const_G);
const integertime_t ti_step = get_integer_timestep(gp->time_bin);
const integertime_t ti_begin =
get_integer_time_begin(ti_current, gp->time_bin);
......@@ -1140,6 +1132,63 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
if (timer) TIMER_TOC(timer_timestep);
}
/**
* @brief End the force calculation of all active particles in a cell
* by multiplying the acccelerations by the relevant constants
*
* @param r The #runner thread.
* @param c The #cell.
* @param timer Are we timing this ?
*/
void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
const struct engine *e = r->e;
const int count = c->count;
const int gcount = c->gcount;
struct part *restrict parts = c->parts;
struct gpart *restrict gparts = c->gparts;
const double const_G = e->physical_constants->const_newton_G;
TIMER_TIC;
/* Anything to do here? */
if (!cell_is_active(c, e)) return;
/* Recurse? */
if (c->split) {
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) runner_do_kick2(r, c->progeny[k], 0);
} else {
/* Loop over the particles in this cell. */
for (int k = 0; k < count; k++) {
/* Get a handle on the part. */
struct part *restrict p = &parts[k];
if (part_is_active(p, e)) {
/* First, finish the force loop */
hydro_end_force(p);
if (p->gpart != NULL) gravity_end_force(p->gpart, const_G);
}
}
/* Loop over the g-particles in this cell. */
for (int k = 0; k < gcount; k++) {
/* Get a handle on the gpart. */
struct gpart *restrict gp = &gparts[k];
if (gp->id_or_neg_offset > 0 && gpart_is_active(gp, e)) {
gravity_end_force(gp, const_G);
}
}
}
if (timer) TIMER_TOC(timer_endforce);
}
/**
* @brief Construct the cell properties from the received particles
*
......@@ -1379,6 +1428,8 @@ void *runner_main(void *data) {
runner_do_kick1(r, ci, 1);
break;
case task_type_kick2:
if (!(e->policy & engine_policy_cooling))
runner_do_end_force(r, ci, 1);
runner_do_kick2(r, ci, 1);
break;
case task_type_timestep:
......@@ -1411,6 +1462,7 @@ void *runner_main(void *data) {
runner_do_grav_fft(r);
break;
case task_type_cooling:
if (e->policy & engine_policy_cooling) runner_do_end_force(r, ci, 1);
runner_do_cooling(r, t->ci, 1);
break;
case task_type_sourceterms:
......
......@@ -36,6 +36,7 @@ enum {
timer_kick1,
timer_kick2,
timer_timestep,
timer_endforce,
timer_dosort,
timer_doself_density,
timer_doself_gradient,
......
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