First really working version of the Minimal SPH implementation. Need more...

First really working version of the Minimal SPH implementation. Need more documentation and testing.

src/const.h

@@ -59,8 +59,8 @@
 //#define WENDLAND_C6_KERNEL
 
 /* SPH variant to use */
-//#define MINIMAL_SPH
-#define GADGET2_SPH
+#define MINIMAL_SPH
+//#define GADGET2_SPH
 //#define DEFAULT_SPH
 
 /* Self gravity stuff. */

src/hydro/Minimal/hydro.h

@@ -18,8 +18,10 @@
  ******************************************************************************/
 
 #include "adiabatic_index.h"
-#include "approx_math.h"
+#include "dimension.h"
 #include "equation_of_state.h"
+#include "hydro_properties.h"
+#include "kernel_hydro.h"
 
 /**
  * @brief Returns the internal energy of a particle
@@ -34,7 +36,7 @@
 __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
     const struct part *restrict p, float dt) {
 
-  return p->u;
+  return p->u + p->u_dt * dt;
 }
 
 /**
@@ -46,7 +48,9 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
 __attribute__((always_inline)) INLINE static float hydro_get_pressure(
     const struct part *restrict p, float dt) {
 
-  return gas_pressure_from_internal_energy(p->rho, p->u);
+  const float u = p->u + p->u_dt * dt;
+
+  return gas_pressure_from_internal_energy(p->rho, u);
 }
 
 /**
@@ -62,7 +66,9 @@ __attribute__((always_inline)) INLINE static float hydro_get_pressure(
 __attribute__((always_inline)) INLINE static float hydro_get_entropy(
     const struct part *restrict p, float dt) {
 
-  return gas_entropy_from_internal_energy(p->rho, p->u);
+  const float u = p->u + p->u_dt * dt;
+
+  return gas_entropy_from_internal_energy(p->rho, u);
 }
 
 /**
@@ -74,7 +80,9 @@ __attribute__((always_inline)) INLINE static float hydro_get_entropy(
 __attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
     const struct part *restrict p, float dt) {
 
-  return gas_soundspeed_from_internal_energy(p->rho, p->u);
+  const float u = p->u + p->u_dt * dt;
+
+  return gas_soundspeed_from_internal_energy(p->rho, u);
 }
 
 /**
@@ -150,7 +158,6 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
   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;
 }
 
 /**
@@ -164,6 +171,7 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
  */
 __attribute__((always_inline)) INLINE static void hydro_init_part(
     struct part *restrict p) {
+
   p->density.wcount = 0.f;
   p->density.wcount_dh = 0.f;
   p->rho = 0.f;
@@ -227,7 +235,10 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
     struct part *restrict p, struct xpart *restrict xp, int ti_current,
     double timeBase) {
 
-  p->force.pressure = p->rho * p->u * hydro_gamma_minus_one;
+  const float half_dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase;
+  const float pressure = hydro_get_pressure(p, half_dt);
+
+  p->force.pressure = pressure;
 }
 
 /**
@@ -268,10 +279,9 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
     struct part *restrict p, const struct xpart *restrict xp, int t0, int t1,
     double timeBase) {
 
-  p->u = xp->u_full;
-
-  /* Need to recompute the pressure as well */
-  p->force.pressure = p->rho * p->u * hydro_gamma_minus_one;
+  /* Drift the pressure */
+  const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase;
+  p->force.pressure = hydro_get_pressure(p, dt_entr);
 }
 
 /**
@@ -304,11 +314,15 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
     struct part *restrict p, struct xpart *restrict xp, float dt,
     float half_dt) {
 
-  /* Kick in momentum space */
-  xp->u_full += p->u_dt * dt;
+  /* Do not decrease the energy by more than a factor of 2*/
+  const float u_change = p->u_dt * dt;
+  if (u_change > -0.5f * p->u)
+    p->u += u_change;
+  else
+    p->u *= 0.5f;
 
-  /* Get the predicted internal energy */
-  p->u = xp->u_full - half_dt * p->u_dt;
+  /* Do not 'overcool' when timestep increases */
+  if (p->u + p->u_dt * half_dt < 0.5f * p->u) p->u_dt = -0.5f * p->u / half_dt;
 }
 
 /**

src/hydro/Minimal/hydro_debug.h

@@ -22,12 +22,12 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
   printf(
       "x=[%.3e,%.3e,%.3e], "
       "v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e], "
-      "u_full=%.3e, u=%.3e, du/dt=%.3e v_sig=%.3e, P=%.3e\n"
-      "h=%.3e, dh/dt=%.3e "
-      "wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, t_begin=%d, t_end=%d\n",
+      "u=%.3e, du/dt=%.3e v_sig=%.3e, P=%.3e\n"
+      "h=%.3e, dh/dt=%.3e wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, "
+      "t_begin=%d, t_end=%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],
-      xp->u_full, p->u, p->u_dt, p->force.v_sig, p->force.pressure, p->h,
-      p->force.h_dt, (int)p->density.wcount, p->mass, p->rho_dh, p->rho,
-      p->ti_begin, p->ti_end);
+      p->u, p->u_dt, p->force.v_sig, p->force.pressure, p->h, p->force.h_dt,
+      (int)p->density.wcount, p->mass, p->rho_dh, p->rho, p->ti_begin,
+      p->ti_end);
 }

src/hydro/Minimal/hydro_iact.h

@@ -115,6 +115,8 @@ runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj,
 __attribute__((always_inline)) INLINE static void runner_iact_force(
     float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
 
+  const float fac_mu = 1.f; /* Will change with cosmological integration */
+
   const float r = sqrtf(r2);
   const float r_inv = 1.0f / r;
 
@@ -143,34 +145,60 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
   const float wj_dr = hjd_inv * wj_dx;
 
   /* Compute gradient terms */
-  const float P_over_rho_i = pressurei / (rhoi * rhoi) * pi->rho_dh;
-  const float P_over_rho_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
+  const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->rho_dh;
+  const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
 
   /* Compute dv dot r. */
   const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
                      (pi->v[1] - pj->v[1]) * dx[1] +
                      (pi->v[2] - pj->v[2]) * dx[2];
 
+  /* Are the particles moving towards each others ? */
+  const float omega_ij = fminf(dvdr, 0.f);
+  const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
+
   /* Compute sound speeds */
   const float ci = sqrtf(hydro_gamma * pressurei / rhoi);
   const float cj = sqrtf(hydro_gamma * pressurej / rhoj);
-  const float v_sig = ci + cj;
+  const float v_sig = ci + cj - 3.f * mu_ij;
+
+  /* Construct the full viscosity term */
+  const float rho_ij = 0.5f * (rhoi + rhoj);
+  const float visc = -const_viscosity_alpha * v_sig * mu_ij / rho_ij;
+
+  /* Convolve with the kernel */
+  const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
 
   /* SPH acceleration term */
-  const float sph_term = (P_over_rho_i * wi_dr + P_over_rho_j * wj_dr) * r_inv;
+  const float sph_acc_term =
+      (P_over_rho2_i * wi_dr + P_over_rho2_j * wj_dr) * r_inv;
+
+  /* Assemble the acceleration */
+  const float acc = sph_acc_term + visc_acc_term;
 
   /* Use the force Luke ! */
-  pi->a_hydro[0] -= mj * sph_term * dx[0];
-  pi->a_hydro[1] -= mj * sph_term * dx[1];
-  pi->a_hydro[2] -= mj * sph_term * dx[2];
+  pi->a_hydro[0] -= mj * acc * dx[0];
+  pi->a_hydro[1] -= mj * acc * dx[1];
+  pi->a_hydro[2] -= mj * acc * dx[2];
 
-  pj->a_hydro[0] += mi * sph_term * dx[0];
-  pj->a_hydro[1] += mi * sph_term * dx[1];
-  pj->a_hydro[2] += mi * sph_term * dx[2];
+  pj->a_hydro[0] += mi * acc * dx[0];
+  pj->a_hydro[1] += mi * acc * dx[1];
+  pj->a_hydro[2] += mi * acc * dx[2];
 
   /* Get the time derivative for u. */
-  pi->u_dt += P_over_rho_i * mj * dvdr * r_inv * wi_dr;
-  pj->u_dt += P_over_rho_j * mi * dvdr * r_inv * wj_dr;
+  const float sph_du_term_i = P_over_rho2_i * dvdr * r_inv * wi_dr;
+  const float sph_du_term_j = P_over_rho2_j * dvdr * r_inv * wi_dr;
+
+  /* Viscosity term */
+  const float visc_du_term = 0.5f * visc_acc_term * dvdr;
+
+  /* Assemble the energy equation term */
+  const float du_dt_i = sph_du_term_i + visc_du_term;
+  const float du_dt_j = sph_du_term_j + visc_du_term;
+
+  /* Internal energy time derivatibe */
+  pi->u_dt += du_dt_i * mj;
+  pj->u_dt += du_dt_j * mi;
 
   /* Get the time derivative for h. */
   pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
@@ -198,6 +226,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_force(
 __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
     float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
 
+  const float fac_mu = 1.f; /* Will change with cosmological integration */
+
   const float r = sqrtf(r2);
   const float r_inv = 1.0f / r;
 
@@ -226,29 +256,53 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
   const float wj_dr = hjd_inv * wj_dx;
 
   /* Compute gradient terms */
-  const float P_over_rho_i = pressurei / (rhoi * rhoi) * pi->rho_dh;
-  const float P_over_rho_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
+  const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->rho_dh;
+  const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
 
   /* Compute dv dot r. */
   const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
                      (pi->v[1] - pj->v[1]) * dx[1] +
                      (pi->v[2] - pj->v[2]) * dx[2];
 
-  /* Compute sound speeds */
+  /* Are the particles moving towards each others ? */
+  const float omega_ij = fminf(dvdr, 0.f);
+  const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
+
+  /* Compute sound speeds and signal velocity */
   const float ci = sqrtf(hydro_gamma * pressurei / rhoi);
   const float cj = sqrtf(hydro_gamma * pressurej / rhoj);
-  const float v_sig = ci + cj;
+  const float v_sig = ci + cj - 3.f * mu_ij;
+
+  /* Construct the full viscosity term */
+  const float rho_ij = 0.5f * (rhoi + rhoj);
+  const float visc = -const_viscosity_alpha * v_sig * mu_ij / rho_ij;
+
+  /* Convolve with the kernel */
+  const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
 
   /* SPH acceleration term */
-  const float sph_term = (P_over_rho_i * wi_dr + P_over_rho_j * wj_dr) * r_inv;
+  const float sph_acc_term =
+      (P_over_rho2_i * wi_dr + P_over_rho2_j * wj_dr) * r_inv;
+
+  /* Assemble the acceleration */
+  const float acc = sph_acc_term + visc_acc_term;
 
   /* Use the force Luke ! */
-  pi->a_hydro[0] -= mj * sph_term * dx[0];
-  pi->a_hydro[1] -= mj * sph_term * dx[1];
-  pi->a_hydro[2] -= mj * sph_term * dx[2];
+  pi->a_hydro[0] -= mj * acc * dx[0];
+  pi->a_hydro[1] -= mj * acc * dx[1];
+  pi->a_hydro[2] -= mj * acc * dx[2];
 
   /* Get the time derivative for u. */
-  pi->u_dt += P_over_rho_i * mj * dvdr * r_inv * wi_dr;
+  const float sph_du_term_i = P_over_rho2_i * dvdr * r_inv * wi_dr;
+
+  /* Viscosity term */
+  const float visc_du_term = 0.5f * visc_acc_term * dvdr;
+
+  /* Assemble the energy equation term */
+  const float du_dt_i = sph_du_term_i + visc_du_term;
+
+  /* Internal energy time derivatibe */
+  pi->u_dt += du_dt_i * mj;
 
   /* Get the time derivative for h. */
   pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;

src/hydro/Minimal/hydro_io.h

@@ -17,6 +17,8 @@
  *
  ******************************************************************************/
 
+#include "adiabatic_index.h"
+#include "hydro.h"
 #include "io_properties.h"
 #include "kernel_hydro.h"
 
@@ -51,6 +53,16 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
                                 UNIT_CONV_DENSITY, parts, rho);
 }
 
+float convert_S(struct engine* e, struct part* p) {
+
+  return hydro_get_entropy(p, 0);
+}
+
+float convert_P(struct engine* e, struct part* p) {
+
+  return hydro_get_pressure(p, 0);
+}
+
 /**
  * @brief Specifies which particle fields to write to a dataset
  *
@@ -61,7 +73,7 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
 void hydro_write_particles(struct part* parts, struct io_props* list,
                            int* num_fields) {
 
-  *num_fields = 8;
+  *num_fields = 10;
 
   /* List what we want to write */
   list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
@@ -80,6 +92,11 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
                                  UNIT_CONV_ACCELERATION, parts, a_hydro);
   list[7] =
       io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts, rho);
+  list[8] = io_make_output_field_convert_part("Entropy", FLOAT, 1,
+                                              UNIT_CONV_ENTROPY_PER_UNIT_MASS,
+                                              parts, rho, convert_S);
+  list[9] = io_make_output_field_convert_part(
+      "Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, rho, convert_P);
 }
 
 /**

src/hydro/Minimal/hydro_part.h

@@ -31,8 +31,6 @@ struct xpart {
 
   float v_full[3]; /*!< Velocity at the last full step. */
 
-  float u_full; /*!< Thermal energy at the last full step. */
-
 } __attribute__((aligned(xpart_align)));
 
 /**