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

Merge branch 'master' into PressureEntropy_SPH

parents ba607295 f4cd515f
......@@ -131,8 +131,10 @@ def parse_files():
file_list = [ file_list[j] for j in sorted_indices]
parse_header(file_list[0])
branch[i] = branch[i].replace("_", "\\_")
version.append(branch[i] + " " + revision[i] + "\n" + hydro_scheme[i] +
version.append("$\\textrm{%s}$"%str(branch[i]) + " " + revision[i] + "\n" + hydro_scheme[i] +
"\n" + hydro_kernel[i] + r", $N_{ngb}=%d$"%float(hydro_neighbours[i]) +
r", $\eta=%.3f$"%float(hydro_eta[i]))
times.append([])
......@@ -215,14 +217,21 @@ def plot_results(times,totalTime,speedUp,parallelEff):
pts = [1, 10**np.floor(np.log10(threadList[i][-1])+1)]
totalTimePlot.loglog(pts,totalTime[i][0]/pts, 'k--', lw=1., color='0.2')
totalTimePlot.loglog(threadList[i],totalTime[i],linestyle[i],label=version[i])
y_min = 10**np.floor(np.log10(np.min(totalTime[:][-1])*0.6))
y_max = 1.2*10**np.floor(np.log10(np.max(totalTime[:][0]) * 1.5)+1)
totalTimePlot.set_xscale('log')
totalTimePlot.set_xlabel("${\\rm Threads}$", labelpad=0.)
totalTimePlot.set_ylabel("${\\rm Time~to~solution}~[{\\rm ms}]$", labelpad=0.)
totalTimePlot.set_xlim([0.9, 10**(np.floor(np.log10(threadList[i][-1]))+0.5)])
totalTimePlot.set_ylim([10**np.floor(np.log10(np.min(totalTime)*0.6)), 1.2*10**np.floor(np.log10(np.max(totalTime) * 1.5)+1)])
totalTimePlot.set_ylim(y_min, y_max)
ax2 = totalTimePlot.twinx()
ax2.set_yscale('log')
ax2.set_ylabel("${\\rm Time~to~solution}~[{\\rm hr}]$", labelpad=0.)
ax2.set_ylim(y_min / 3.6e6, y_max / 3.6e6)
totalTimePlot.legend(bbox_to_anchor=(1.14, 0.97), loc=2, borderaxespad=0.,prop={'size':12}, frameon=False)
totalTimePlot.legend(bbox_to_anchor=(1.21, 0.97), loc=2, borderaxespad=0.,prop={'size':12}, frameon=False)
emptyPlot.axis('off')
for i, txt in enumerate(threadList[0]):
......
......@@ -2818,6 +2818,8 @@ void engine_step(struct engine *e) {
submask |= 1 << task_subtype_tend;
}
if (e->verbose) engine_print_task_counts(e);
/* Send off the runners. */
TIMER_TIC;
engine_launch(e, e->nr_threads, mask, submask);
......
......@@ -231,18 +231,18 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
p->density.wcount *= kernel_norm;
p->density.wcount_dh *= h_inv * kernel_gamma * kernel_norm;
const float irho = 1.f / p->rho;
const float rho_inv = 1.f / p->rho;
/* Compute the derivative term */
p->rho_dh = 1.f / (1.f + hydro_dimension_inv * p->h * p->rho_dh * irho);
p->rho_dh = 1.f / (1.f + hydro_dimension_inv * p->h * p->rho_dh * rho_inv);
/* Finish calculation of the velocity curl components */
p->density.rot_v[0] *= h_inv_dim_plus_one * irho;
p->density.rot_v[1] *= h_inv_dim_plus_one * irho;
p->density.rot_v[2] *= h_inv_dim_plus_one * irho;
p->density.rot_v[0] *= h_inv_dim_plus_one * rho_inv;
p->density.rot_v[1] *= h_inv_dim_plus_one * rho_inv;
p->density.rot_v[2] *= h_inv_dim_plus_one * rho_inv;
/* Finish calculation of the velocity divergence */
p->density.div_v *= h_inv_dim_plus_one * irho;
p->density.div_v *= h_inv_dim_plus_one * rho_inv;
}
/**
......@@ -273,13 +273,13 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
const float half_dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase;
const float pressure = hydro_get_pressure(p, half_dt);
const float irho = 1.f / p->rho;
const float rho_inv = 1.f / p->rho;
/* Divide the pressure by the density and density gradient */
const float P_over_rho2 = pressure * irho * irho * p->rho_dh;
const float P_over_rho2 = pressure * rho_inv * rho_inv * p->rho_dh;
/* Compute the sound speed */
const float soundspeed = sqrtf(hydro_gamma * pressure * irho);
const float soundspeed = sqrtf(hydro_gamma * pressure * rho_inv);
/* Compute the Balsara switch */
const float balsara =
......@@ -349,13 +349,13 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase;
const float pressure = hydro_get_pressure(p, dt_entr);
const float irho = 1.f / p->rho;
const float rho_inv = 1.f / p->rho;
/* Divide the pressure by the density and density gradient */
const float P_over_rho2 = pressure * irho * irho * p->rho_dh;
const float P_over_rho2 = pressure * rho_inv * rho_inv * p->rho_dh;
/* Compute the new sound speed */
const float soundspeed = sqrtf(hydro_gamma * pressure * irho);
const float soundspeed = sqrtf(hydro_gamma * pressure * rho_inv);
/* Update variables */
p->force.P_over_rho2 = P_over_rho2;
......
......@@ -248,17 +248,17 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
const float ri = 1.0f / r;
/* Compute the kernel function */
const float h_inv = 1.0f / hi;
const float u = r * h_inv;
kernel_deval(u, &wi, &wi_dx);
const float hi_inv = 1.0f / hi;
const float ui = r * hi_inv;
kernel_deval(ui, &wi, &wi_dx);
/* Compute contribution to the density */
pi->rho += mj * wi;
pi->rho_dh -= mj * (hydro_dimension * wi + u * wi_dx);
pi->rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
/* Compute contribution to the number of neighbours */
pi->density.wcount += wi;
pi->density.wcount_dh -= u * wi_dx;
pi->density.wcount_dh -= ui * wi_dx;
const float fac = mj * wi_dx * ri;
......
......@@ -1269,6 +1269,12 @@ void space_split_mapper(void *map_data, int num_elements, void *extra_data) {
atomic_cas(&s->maxdepth, maxdepth, c->depth);
}
/* If the depth is too large, we have a problem and should stop. */
if (maxdepth > space_cell_maxdepth) {
error("Exceeded maximum depth (%d) when splitting cells, aborting",
space_cell_maxdepth);
}
/* Split or let it be? */
if (count > space_splitsize || gcount > space_splitsize) {
......
......@@ -44,6 +44,9 @@
#define space_stretch 1.10f
#define space_maxreldx 0.25f
/* Maximum allowed depth of cell splits. */
#define space_cell_maxdepth 52
/* Split size. */
extern int space_splitsize;
extern int space_maxsize;
......
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