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

Added a 'sound speed' variable to the Minimal SPH particle to save redundant calculations.

parent 6f8b4ff5
...@@ -63,8 +63,8 @@ ...@@ -63,8 +63,8 @@
//#define WENDLAND_C6_KERNEL //#define WENDLAND_C6_KERNEL
/* SPH variant to use */ /* SPH variant to use */
//#define MINIMAL_SPH #define MINIMAL_SPH
#define GADGET2_SPH //#define GADGET2_SPH
//#define HOPKINS_PE_SPH //#define HOPKINS_PE_SPH
//#define DEFAULT_SPH //#define DEFAULT_SPH
//#define GIZMO_SPH //#define GIZMO_SPH
......
...@@ -141,6 +141,12 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy( ...@@ -141,6 +141,12 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
/* Compute the pressure */ /* Compute the pressure */
const float pressure = gas_pressure_from_internal_energy(p->rho, p->u); const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
/* Compute the sound speed from the pressure*/
const float rho_inv = 1.f / p->rho;
const float soundspeed = sqrtf(hydro_gamma * pressure * rho_inv);
p->force.soundspeed = soundspeed;
p->force.pressure = pressure; p->force.pressure = pressure;
} }
...@@ -160,6 +166,12 @@ __attribute__((always_inline)) INLINE static void hydro_set_entropy( ...@@ -160,6 +166,12 @@ __attribute__((always_inline)) INLINE static void hydro_set_entropy(
/* Compute the pressure */ /* Compute the pressure */
const float pressure = gas_pressure_from_internal_energy(p->rho, p->u); const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
/* Compute the sound speed from the pressure*/
const float rho_inv = 1.f / p->rho;
const float soundspeed = sqrtf(hydro_gamma * pressure * rho_inv);
p->force.soundspeed = soundspeed;
p->force.pressure = pressure; p->force.pressure = pressure;
} }
...@@ -277,10 +289,15 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( ...@@ -277,10 +289,15 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
struct part *restrict p, struct xpart *restrict xp, int ti_current, struct part *restrict p, struct xpart *restrict xp, int ti_current,
double timeBase) { double timeBase) {
/* Compute the pressure */
const float half_dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase; 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 pressure = hydro_get_pressure(p, half_dt);
const float rho_inv = 1.f / p->rho; const float rho_inv = 1.f / p->rho;
/* Compute the sound speed */
const float soundspeed = sqrtf(hydro_gamma * pressure * rho_inv);
/* Compute the "grad h" term */ /* Compute the "grad h" term */
const float grad_h_term = const float grad_h_term =
1.f / (1.f + hydro_dimension_inv * p->h * p->density.rho_dh * rho_inv); 1.f / (1.f + hydro_dimension_inv * p->h * p->density.rho_dh * rho_inv);
...@@ -288,6 +305,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( ...@@ -288,6 +305,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Update variables. */ /* Update variables. */
p->force.f = grad_h_term; p->force.f = grad_h_term;
p->force.pressure = pressure; p->force.pressure = pressure;
p->force.soundspeed = soundspeed;
} }
/** /**
...@@ -347,7 +365,15 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( ...@@ -347,7 +365,15 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* Drift the pressure */ /* Drift the pressure */
const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase; const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase;
p->force.pressure = hydro_get_pressure(p, dt_entr); const float pressure = hydro_get_pressure(p, dt_entr);
const float rho_inv = 1.f / p->rho;
/* Compute the new sound speed */
const float soundspeed = sqrtf(hydro_gamma * pressure * rho_inv);
p->force.pressure = pressure;
p->force.soundspeed = soundspeed;
} }
/** /**
...@@ -392,6 +418,12 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra( ...@@ -392,6 +418,12 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
/* Compute the pressure */ /* Compute the pressure */
const float pressure = gas_pressure_from_internal_energy(p->rho, p->u); const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
/* Compute the sound speed from the pressure*/
const float rho_inv = 1.f / p->rho;
const float soundspeed = sqrtf(hydro_gamma * pressure * rho_inv);
p->force.soundspeed = soundspeed;
p->force.pressure = pressure; p->force.pressure = pressure;
} }
......
...@@ -165,8 +165,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( ...@@ -165,8 +165,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */ const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
/* Compute sound speeds and signal velocity */ /* Compute sound speeds and signal velocity */
const float ci = sqrtf(hydro_gamma * pressurei / rhoi); const float ci = pi->force.soundspeed;
const float cj = sqrtf(hydro_gamma * pressurej / rhoj); const float cj = pj->force.soundspeed;
const float v_sig = ci + cj - 3.f * mu_ij; const float v_sig = ci + cj - 3.f * mu_ij;
/* Construct the full viscosity term */ /* Construct the full viscosity term */
...@@ -276,8 +276,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( ...@@ -276,8 +276,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */ const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
/* Compute sound speeds and signal velocity */ /* Compute sound speeds and signal velocity */
const float ci = sqrtf(hydro_gamma * pressurei / rhoi); const float ci = pi->force.soundspeed;
const float cj = sqrtf(hydro_gamma * pressurej / rhoj); const float cj = pj->force.soundspeed;
const float v_sig = ci + cj - 3.f * mu_ij; const float v_sig = ci + cj - 3.f * mu_ij;
/* Construct the full viscosity term */ /* Construct the full viscosity term */
......
...@@ -131,6 +131,9 @@ struct part { ...@@ -131,6 +131,9 @@ struct part {
/*! Particle pressure. */ /*! Particle pressure. */
float pressure; float pressure;
/*! Particle soundspeed. */
float soundspeed;
/*! Particle signal velocity */ /*! Particle signal velocity */
float v_sig; float v_sig;
......
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