diff --git a/src/const.h b/src/const.h
index 82f5b74604ba315efbbe01c922cb30b973e56965..01326e3a22b3f608dfa9c0cf25a1dffe41697edc 100644
--- a/src/const.h
+++ b/src/const.h
@@ -38,7 +38,7 @@
1.f /* Value taken from (Price,2008), not used in legacy gadget mode */
/* Time integration constants. */
-#define const_cfl 0.3f
+#define const_cfl 0.2f
#define const_ln_max_h_change \
0.231111721f /* Particle can't change volume by more than a factor of \
2=1.26^3 over one time step */
@@ -47,7 +47,7 @@
/* Neighbour search constants. */
#define const_eta_kernel \
1.2349f /* Corresponds to 48 ngbs with the cubic spline kernel */
-#define const_delta_nwneigh 1.f
+#define const_delta_nwneigh 0.1f
#define const_smoothing_max_iter 30
#define CUBIC_SPLINE_KERNEL
diff --git a/src/debug.c b/src/debug.c
index c6f37f3b2651c12cd1c016b6a27bdd097c92ebae..aa9edfe5b173ccf50e7bd82b4f2643fe30e7040a 100644
--- a/src/debug.c
+++ b/src/debug.c
@@ -50,7 +50,7 @@ void printParticle(struct part *parts, long long int id, int N) {
"v=[%.3e,%.3e,%.3e], a=[%.3e,%.3e,%.3e],\n h=%.3e, "
"wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, P=%.3e, S=%.3e, dS/dt=%.3e,\n"
"divV=%.3e, curlV=%.3e, rotV=[%.3e,%.3e,%.3e] \n "
- "t_begin=%.3e, t_end=%.3e\n",
+ "v_sig=%e t_begin=%.3e, t_end=%.3e\n",
i, parts[i].id, parts[i].x[0], parts[i].x[1], parts[i].x[2],
parts[i].v[0], parts[i].v[1], parts[i].v[2], parts[i].a[0],
parts[i].a[1], parts[i].a[2], 2.*parts[i].h,
@@ -65,6 +65,7 @@ void printParticle(struct part *parts, long long int id, int N) {
parts[i].rot_v[0],
parts[i].rot_v[1],
parts[i].rot_v[2],
+ parts[i].v_sig,
parts[i].t_begin, parts[i].t_end);
found = 1;
}
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index bc18400927b1dfa1d3e81f347cfca08c64e5a707..cbfcea09c4d59f153fa09c0e12adc039c2e42d3e 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -27,21 +27,16 @@
__attribute__((always_inline)) INLINE static float hydro_compute_timestep(
struct part* p, struct xpart* xp) {
- /* CFL condition */
- const float dt_cfl = const_cfl * p->h / p->v_sig;
-
- /* /\* Limit change in h *\/ */
- /* float dt_h_change = (p->h_dt != 0.0f) */
- /* ? fabsf(const_ln_max_h_change * p->h / p->h_dt) */
- /* : FLT_MAX; */
+ /* Acceleration */
+ float ac = sqrtf(p->a[0]*p->a[0] + p->a[1]*p->a[1] + p->a[2]*p->a[2]);
+ ac = fmaxf(ac, 1e-30);
- /* /\* Limit change in u *\/ */
- /* float dt_u_change = (p->force.u_dt != 0.0f) */
- /* ? fabsf(const_max_u_change * p->u / p->force.u_dt) */
- /* : FLT_MAX; */
+ const float dt_accel = sqrtf(2.f);
+
+ /* CFL condition */
+ const float dt_cfl = 2.f * const_cfl * p->h / p->v_sig;
- // return fminf(dt_cfl, fminf(dt_h_change, dt_u_change));
- return dt_cfl;
+ return fminf(dt_cfl, dt_accel);
}
@@ -176,8 +171,8 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(struc
/* Reset the time derivatives. */
p->entropy_dt = 0.0f;
- /* Reset minimal signal velocity */
- p->v_sig = FLT_MAX;
+ /* Reset maximal signal velocity */
+ p->v_sig = 0.0f;
}
@@ -228,6 +223,6 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(struct part* p
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(struct part* p) {
- p->entropy = (const_hydro_gamma - 1.f) * p->entropy * pow(p->rho, -(const_hydro_gamma - 1.f));
+ p->entropy = (const_hydro_gamma - 1.f) * p->entropy * powf(p->rho, -(const_hydro_gamma - 1.f));
}
diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index 8018807cb54a84cee8d57345e5ea37a9c521a304..d190501d9fd421162ca5051bcc06ecf868d78a2b 100644
--- a/src/hydro/Gadget2/hydro_iact.h
+++ b/src/hydro/Gadget2/hydro_iact.h
@@ -113,13 +113,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
pj->rot_v[0] += facj * curlvr[0];
pj->rot_v[1] += facj * curlvr[1];
pj->rot_v[2] += facj * curlvr[2];
-
- /* if(pi->id == 1000) */
- /* message("Interacting with %lld. r=%f\n", pj->id, r); */
-
- /* if(pj->id == 1000) */
- /* message("Interacting with %lld. r=%f\n", pi->id, r); */
-
}
@@ -139,8 +132,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
/* Get r and r inverse. */
const float r = sqrtf(r2);
const float ri = 1.0f / r;
-
- if(pi->id == 1000 && pj->id == 1103) hi = 0.2234976 / 2.f;
/* Compute the kernel function */
const float h_inv = 1.0f / hi;
@@ -167,7 +158,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2];
pi->div_v -= fac * dvdr;
- if(pi->id == 1000 && pj->id == 1103)
+ if(pi->id == 1000 && pj->id == 1203)
message("Interacting with %lld. r=%e hi=%e u=%e W=%e dW/dx=%e dh_drho1=%e dh_drho2=%e\n fac=%e dvdr=%e",
pj->id,
r,
@@ -260,6 +251,23 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Eventually got the acceleration */
const float acc = visc_term + sph_term;
+
+ if(pi->id == 1000 && pj->id == 1203)
+ message("Interacting with %lld. r=%e hi=%e hj=%e dWi/dx=%e dWj/dx=%3e dvdr=%e visc=%e sph=%e",
+ pj->id,
+ r,
+ 2*hi,
+ 2*hj,
+ wi_dr,
+ wj_dr,
+ dvdr,
+ visc_term,
+ sph_term
+ );
+ if(pi->id == 1203 && pj->id == 1000)
+ message("oO");
+
+
/* Use the force Luke ! */
pi->a[0] -= acc * dx[0];
pi->a[1] -= acc * dx[1];
@@ -274,8 +282,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
pj->v_sig = fmaxf(pj->v_sig, v_sig) ;
/* Change in entropy */
- pi->entropy_dt += 0.5f * visc_term * dvdr;
- pj->entropy_dt -= 0.5f * visc_term * dvdr;
+ pi->entropy_dt -= 0.5f * visc_term * dvdr;
+ pj->entropy_dt += 0.5f * visc_term * dvdr;
}
@@ -345,7 +353,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Eventually got the acceleration */
const float acc = visc_term + sph_term;
-
+
/* Use the force Luke ! */
pi->a[0] -= acc * dx[0];
pi->a[1] -= acc * dx[1];
@@ -355,7 +363,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
pi->v_sig = fmaxf(pi->v_sig, v_sig) ;
/* Change in entropy */
- pi->entropy_dt += 0.5f * visc_term * dvdr;
+ pi->entropy_dt -= 0.5f * visc_term * dvdr;
}
diff --git a/src/runner.c b/src/runner.c
index bbfe78032113f905f1fb737eafd743dcb3b673fb..848076515a18cccd614f2a9c1b290fe4c292dd5b 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -725,7 +725,7 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) {
p->x[0] += xp->v_full[0] * dt;
p->x[1] += xp->v_full[1] * dt;
p->x[2] += xp->v_full[2] * dt;
-
+
/* Predict velocities */
p->v[0] += p->a[0] * dt;
p->v[1] += p->a[1] * dt;
@@ -860,6 +860,13 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
new_dt = fminf(new_dt_hydro, new_dt_grav);
+ if(p->id == 1000)
+ message("1000 dt_hydro=%e", new_dt_hydro);
+
+ if(p->id == 515050)
+ message("515050 dt_hydro=%e", new_dt_hydro);
+
+
/* Recover the current timestep */
const float current_dt = p->t_end - p->t_begin;