Skip to content
GitLab
Commit 2bd3be76 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

First working version of PE-SPH

parent cb5ebdeb
Hide whitespace changes
Inline Side-by-side
src/hydro/PressureEntropy/hydro.h
View file @ 2bd3be76
......@@ -48,7 +48,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
const float entropy = p->entropy + p->entropy_dt * dt;
return gas_internal_energy_from_entropy(p->rho, entropy);
return gas_internal_energy_from_entropy(p->rho_bar, entropy);
}
/**
......@@ -62,7 +62,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_pressure(
const float entropy = p->entropy + p->entropy_dt * dt;
return gas_pressure_from_entropy(p->rho, entropy);
return gas_pressure_from_entropy(p->rho_bar, entropy);
}
/**
......@@ -124,7 +124,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
__attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
struct part *restrict p, float u) {
p->entropy = gas_entropy_from_internal_energy(p->rho, u);
p->entropy = gas_entropy_from_internal_energy(p->rho_bar, u);
}
/**
......@@ -194,9 +194,11 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
p->density.wcount = 0.f;
p->density.wcount_dh = 0.f;
p->rho = 0.f;
p->weightedPressure = 0.f;
p->density.weightedPressure_dh = 0.f;
// p->rho_dh = 0.f;
p->rho_dh = 0.f;
p->rho_bar = 0.f;
p->pressure_dh = 0.f;
// p->density.weightedPressure_dh = 0.f;
/* p->density.div_v = 0.f; */
/* p->density.rot_v[0] = 0.f; */
/* p->density.rot_v[1] = 0.f; */
......@@ -217,32 +219,36 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
/* Some smoothing length multiples. */
const float h = p->h;
const float h_inv = 1.0f / h; /* 1/h */
const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
// const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
const float h_inv = 1.0f / h; /* 1/h */
const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
/* Final operation on the density (add self-contribution). */
p->rho += p->mass * kernel_root;
p->rho_dh -= hydro_dimension * p->mass * kernel_root;
p->rho_bar += p->mass * p->entropy_one_over_gamma * kernel_root;
p->pressure_dh -=
hydro_dimension * p->mass * p->entropy_one_over_gamma * kernel_root;
p->density.wcount += kernel_root;
// p->density.wcount_dh -= hydro_dimension * kernel_root;
p->weightedPressure += p->mass * pow_one_over_gamma(p->entropy) * kernel_root;
p->density.weightedPressure_dh -=
hydro_dimension * p->mass * pow_one_over_gamma(p->entropy) * kernel_root;
/* Finish the calculation by inserting the missing h-factors */
p->rho *= h_inv_dim;
p->rho_dh *= h_inv_dim;
p->rho_bar *= h_inv_dim;
p->pressure_dh *= h_inv_dim_plus_one;
p->density.wcount *= kernel_norm;
p->density.wcount_dh *= h_inv * kernel_gamma * kernel_norm;
p->weightedPressure *= h_inv_dim;
p->density.weightedPressure_dh *= h_inv * kernel_gamma * kernel_norm;
/* Final operation on the weighted pressure */
p->weightedPressure = pow_gamma(p->weightedPressure);
const float rho_inv = 1.f / p->rho;
const float entropy_minus_one_over_gamma = 1.f / p->entropy_one_over_gamma;
/* const float irho = 1.f / p->rho; */
/* Final operation on the weighted density */
p->rho_bar *= entropy_minus_one_over_gamma;
/* 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);
p->pressure_dh *=
p->h * rho_inv * hydro_dimension_inv * entropy_minus_one_over_gamma;
/* Finish calculation of the velocity curl components */
// p->density.rot_v[0] *= h_inv_dim_plus_one * irho;
......@@ -269,25 +275,15 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Compute the pressure */
const float half_dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase;
/* Compute the sound speed from the actual pressure*/
const float pressure = hydro_get_pressure(p, half_dt);
const float irho = 1.f / p->rho;
const float soundspeed = sqrtf(hydro_gamma * pressure * irho);
/* Compute the "i" part of the f_ij term */
const float number_density = p->density.wcount / kernel_norm;
const float factor = p->h / (hydro_dimension * number_density);
const float f_ij =
factor * p->density.weightedPressure_dh / (1.f + factor * number_density);
/* Comput the pressure term */
const float pressure_term = pow_one_minus_two_over_gamma(p->weightedPressure);
/* Compute the sound speed from the pressure*/
const float rho_bar_inv = 1.f / p->rho_bar;
const float soundspeed = sqrtf(hydro_gamma * pressure * rho_bar_inv);
/* Update variables. */
p->force.soundspeed = soundspeed;
p->force.f_ij = f_ij;
p->force.pressure_term = pressure_term;
p->force.P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
}
/**
......@@ -328,26 +324,40 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
struct part *restrict p, const struct xpart *restrict xp, float dt, int t0,
int t1, double timeBase) {
/* Drift the pressure */
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 h_inv = 1.f / p->h;
/* Compute the new sound speed */
const float soundspeed = sqrtf(hydro_gamma * pressure * irho);
const float w1 = -hydro_dimension * p->force.h_dt * dt_entr / p->h;
/* Predict smoothing length */
const float w1 = p->force.h_dt * h_inv * dt;
if (fabsf(w1) < 0.2f)
p->weightedPressure *= approx_expf(w1); /* 4th order expansion of exp(w) */
p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
else
p->weightedPressure *= expf(w1);
p->h *= expf(w1);
/* Predict density */
const float w2 = -hydro_dimension * w1;
if (fabsf(w2) < 0.2f) {
p->rho *= approx_expf(w2); /* 4th order expansion of exp(w) */
p->rho_bar *= approx_expf(w2);
} else {
p->rho *= expf(w2);
p->rho_bar *= expf(w2);
}
/* Drift the entropy */
const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase;
const float entropy = hydro_get_entropy(p, dt_entr);
/* Compute the pressure */
const float pressure = gas_pressure_from_entropy(p->rho_bar, entropy);
const float pressure_term = pow_one_minus_two_over_gamma(p->weightedPressure);
/* Compute the sound speed from the pressure*/
const float rho_bar_inv = 1.f / p->rho_bar;
const float soundspeed = sqrtf(hydro_gamma * pressure * rho_bar_inv);
/* Update variables */
/* Update the variables */
p->entropy_one_over_gamma = pow_one_over_gamma(entropy);
p->force.soundspeed = soundspeed;
p->force.pressure_term = pressure_term;
p->force.P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
}
/**
......@@ -363,7 +373,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
p->force.h_dt *= p->h * hydro_dimension_inv;
p->entropy_dt *=
0.5f * hydro_gamma_minus_one * pow_minus_gamma_minus_one(p->rho);
0.5f * hydro_gamma_minus_one * pow_minus_gamma_minus_one(p->rho_bar);
}
/**
......@@ -388,6 +398,17 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
/* Do not 'overcool' when timestep increases */
if (p->entropy + p->entropy_dt * half_dt < 0.5f * p->entropy)
p->entropy_dt = -0.5f * p->entropy / half_dt;
/* Compute the pressure */
const float pressure = gas_pressure_from_entropy(p->rho_bar, p->entropy);
/* Compute the sound speed from the pressure*/
const float rho_bar_inv = 1.f / p->rho_bar;
const float soundspeed = sqrtf(hydro_gamma * pressure * rho_bar_inv);
p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);
p->force.soundspeed = soundspeed;
p->force.P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
}
/**
......@@ -402,6 +423,17 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
/* We read u in the entropy field. We now get S from u */
p->entropy = gas_entropy_from_internal_energy(p->rho, p->entropy);
p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);
/* Compute the pressure */
const float entropy = p->entropy;
const float pressure = gas_pressure_from_entropy(p->rho_bar, entropy);
/* Compute the sound speed from the pressure*/
const float rho_bar_inv = 1.f / p->rho_bar;
const float soundspeed = sqrtf(hydro_gamma * pressure * rho_bar_inv);
p->force.soundspeed = soundspeed;
p->force.P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
}
#endif /* SWIFT_PRESSURE_ENTROPY_HYDRO_H */
src/hydro/PressureEntropy/hydro_iact.h
View file @ 2bd3be76
......@@ -43,9 +43,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
const float mi = pi->mass;
const float mj = pj->mass;
/* Get the entropies. */
const float entropy_i = pi->entropy;
const float entropy_j = pj->entropy;
/* /\* Get the entropies. *\/ */
/* const float entropy_i = pi->entropy; */
/* const float entropy_j = pj->entropy; */
/* Get r and r inverse. */
const float r = sqrtf(r2);
......@@ -58,15 +58,16 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
/* Compute contribution to the density */
pi->rho += mj * wi;
/* Weighted pressure */
pi->weightedPressure += mj * pow_one_over_gamma(pj->entropy) * wi;
pi->density.weightedPressure_dh -=
mj * pow_one_over_gamma(entropy_j) * (hydro_dimension * wi + ui * 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 -= ui * wi_dx; //(hydro_dimension * wi + ui * wi_dx);
pi->density.wcount_dh -= ui * wi_dx;
/* Compute contribution to the weighted density */
pi->rho_bar += mj * pj->entropy_one_over_gamma * wi;
pi->pressure_dh -=
mj * pj->entropy_one_over_gamma * (hydro_dimension * wi + ui * wi_dx);
/* Compute the kernel function for pj */
const float hj_inv = 1.f / hj;
......@@ -75,16 +76,16 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
/* Compute contribution to the density */
pj->rho += mi * wj;
// pj->rho_dh -= mi * (hydro_dimension * wj + uj * wj_dx);
/* Weighted pressure */
pj->weightedPressure += mi * pow_one_over_gamma(pi->entropy) * wj;
pj->density.weightedPressure_dh -=
mi * pow_one_over_gamma(entropy_i) * (hydro_dimension * wj + uj * wj_dx);
pj->rho_dh -= mi * (hydro_dimension * wj + uj * wj_dx);
/* Compute contribution to the number of neighbours */
pj->density.wcount += wj;
pj->density.wcount_dh -= uj * wj_dx; //(hydro_dimension * wj + uj * wj_dx);
pj->density.wcount_dh -= uj * wj_dx;
/* Compute contribution to the weighted density */
pj->rho_bar += mi * pi->entropy_one_over_gamma * wj;
pj->pressure_dh -=
mi * pi->entropy_one_over_gamma * (hydro_dimension * wj + uj * wj_dx);
/* const float faci = mj * wi_dx * r_inv; */
/* const float facj = mi * wj_dx * r_inv; */
......@@ -135,7 +136,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
const float mj = pj->mass;
/* Get the entropies. */
const float entropy_j = pj->entropy;
// const float entropy_j = pj->entropy;
/* Get r and r inverse. */
const float r = sqrtf(r2);
......@@ -143,20 +144,21 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
/* 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 ui = r * h_inv;
kernel_deval(ui, &wi, &wi_dx);
/* Compute contribution to the density */
pi->rho += mj * wi;
/* Weighted pressure */
pi->weightedPressure += mj * pow_one_over_gamma(pj->entropy) * wi;
pi->density.weightedPressure_dh -=
mj * pow_one_over_gamma(entropy_j) * (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; //(hydro_dimension * wi + u * wi_dx);
pi->density.wcount_dh -= ui * wi_dx; //(hydro_dimension * wi + u * wi_dx);
/* Compute contribution to the weighted density */
pi->rho_bar += mj * pj->entropy_one_over_gamma * wi;
pi->pressure_dh -=
mj * pj->entropy_one_over_gamma * (hydro_dimension * wi + ui * wi_dx);
/* const float fac = mj * wi_dx * ri; */
......@@ -205,9 +207,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float mj = pj->mass;
const float rhoi = pi->rho;
const float rhoj = pj->rho;
const float entropy_i = pi->entropy;
const float entropy_j = pj->entropy;
const float entropy_product = pow_one_over_gamma(entropy_i * entropy_j);
/* Get the kernel for hi. */
const float hi_inv = 1.0f / hi;
......@@ -224,10 +223,16 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float wj_dr = hjd_inv * wj_dx;
/* Compute gradient terms */
const float f_ij =
1.f; // - pi->force.f_ij / (mj * pow_one_over_gamma(entropy_j));
const float f_ji =
1.f; // - pj->force.f_ij / (mi * pow_one_over_gamma(entropy_i));
const float f_i = pi->rho_dh * pi->pressure_dh;
const float f_j = pj->rho_dh * pj->pressure_dh;
/* Compute Pressure terms */
const float P_over_rho2_i = pi->force.P_over_rho2;
const float P_over_rho2_j = pj->force.P_over_rho2;
/* Compute entropy terms */
const float S_gamma_i = pi->entropy_one_over_gamma;
const float S_gamma_j = pj->entropy_one_over_gamma;
/* Compute sound speeds */
const float ci = pi->force.soundspeed;
......@@ -254,22 +259,12 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float visc = -0.5f * const_viscosity_alpha * v_sig * mu_ij / rho_ij;
/* Now, convolve with the kernel */
const float visc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
const float sph_term = entropy_product *
(f_ij * pi->force.pressure_term * wi_dr +
f_ji * pj->force.pressure_term * wj_dr) *
r_inv;
/* if(pi->entropy != pj->entropy) { */
/* message("Si = %f Pi = %f, Pi_term = %f", pi->entropy, pi->weightedPressure,
* pi->force.pressure_term); */
/* message("Sj = %f Pj = %f, Pj_term = %f", pj->entropy, pj->weightedPressure,
* pj->force.pressure_term); */
/* error("done"); */
/* } */
const float visc_term = 0.5f * visc * (wi_dr + wj_dr);
const float sph_term = (S_gamma_j / S_gamma_i - f_i) * P_over_rho2_i * wi_dr +
(S_gamma_i / S_gamma_j - f_j) * P_over_rho2_j * wj_dr;
/* Eventually got the acceleration */
const float acc = visc_term + sph_term;
const float acc = (visc_term + sph_term) * r_inv;
/* Use the force Luke ! */
pi->a_hydro[0] -= mj * acc * dx[0];
......@@ -289,8 +284,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
pj->force.v_sig = fmaxf(pj->force.v_sig, v_sig);
/* Change in entropy */
pi->entropy_dt += mj * visc_term * dvdr;
pj->entropy_dt += mi * visc_term * dvdr;
pi->entropy_dt += mj * visc_term * r_inv * dvdr;
pj->entropy_dt += mi * visc_term * r_inv * dvdr;
}
/**
......@@ -321,9 +316,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float mj = pj->mass;
const float rhoi = pi->rho;
const float rhoj = pj->rho;
const float entropy_i = pi->entropy;
const float entropy_j = pj->entropy;
const float entropy_product = pow_one_over_gamma(entropy_i * entropy_j);
/* Get the kernel for hi. */
const float hi_inv = 1.0f / hi;
......@@ -340,10 +332,16 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float wj_dr = hjd_inv * wj_dx;
/* Compute gradient terms */
const float f_ij =
1.f; // - pi->force.f_ij / (mj * pow_one_over_gamma(entropy_j));
const float f_ji =
1.f; // - pj->force.f_ij / (mi * pow_one_over_gamma(entropy_i));
const float f_i = pi->rho_dh * pi->pressure_dh;
const float f_j = pj->rho_dh * pj->pressure_dh;
/* Compute Pressure terms */
const float P_over_rho2_i = pi->force.P_over_rho2;
const float P_over_rho2_j = pj->force.P_over_rho2;
/* Compute entropy terms */
const float S_gamma_i = pi->entropy_one_over_gamma;
const float S_gamma_j = pj->entropy_one_over_gamma;
/* Compute sound speeds */
const float ci = pi->force.soundspeed;
......@@ -370,14 +368,12 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float visc = -0.5f * const_viscosity_alpha * v_sig * mu_ij / rho_ij;
/* Now, convolve with the kernel */
const float visc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
const float sph_term = entropy_product *
(f_ij * pi->force.pressure_term * wi_dr +
f_ji * pj->force.pressure_term * wj_dr) *
r_inv;
const float visc_term = 0.5f * visc * (wi_dr + wj_dr);
const float sph_term = (S_gamma_j / S_gamma_i - f_i) * P_over_rho2_i * wi_dr +
(S_gamma_i / S_gamma_j - f_j) * P_over_rho2_j * wj_dr;
/* Eventually got the acceleration */
const float acc = visc_term + sph_term;
const float acc = (visc_term + sph_term) * r_inv;
/* Use the force Luke ! */
pi->a_hydro[0] -= mj * acc * dx[0];
......@@ -391,7 +387,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
pi->force.v_sig = fmaxf(pi->force.v_sig, v_sig);
/* Change in entropy */
pi->entropy_dt += mj * visc_term * dvdr;
pi->entropy_dt += mj * visc_term * r_inv * dvdr;
}
/**
......
src/hydro/PressureEntropy/hydro_io.h
View file @ 2bd3be76
......@@ -113,8 +113,8 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
"Entropy", FLOAT, 1, UNIT_CONV_ENTROPY_PER_UNIT_MASS, parts, entropy);
list[9] = io_make_output_field_convert_part(
"Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, rho, convert_P);
list[10] = io_make_output_field("WeightedPressure", FLOAT, 1,
UNIT_CONV_PRESSURE, parts, weightedPressure);
list[10] = io_make_output_field("WeightedDensity", FLOAT, 1,
UNIT_CONV_DENSITY, parts, rho_bar);
}
/**
......
src/hydro/PressureEntropy/hydro_part.h
View file @ 2bd3be76
......@@ -72,8 +72,17 @@ struct part {
/*! Particle density. */
float rho;
/*! Particle weighted pressure. */
float weightedPressure;
/*! Derivative of density with respect to h */
float rho_dh;
/*! Particle pressure. */
float pressure;
/*! Derivative of pressure with respect to h */
float pressure_dh;
/*! Particle weighted density */
float rho_bar;
/*! Particle entropy. */
float entropy;
......@@ -81,6 +90,9 @@ struct part {
/*! Entropy time derivative */
float entropy_dt;
/*! Particle entropy to the power 1/gamma. */
float entropy_one_over_gamma;
union {
struct {
......@@ -92,23 +104,23 @@ struct part {
float wcount_dh;
/*! Derivative of particle weighted pressure with h. */
float weightedPressure_dh;
// float weightedPressure_dh;
/*! Particle velocity curl. */
// float rot_v[3];
/*! Particle velocity divergence. */
// float div_v;
float div_v;
} density;
struct {
/*! "Grad h" term */
float f_ij;
// float f_ij;
/*! Pressure term */
float pressure_term;
float P_over_rho2;
/*! Particle sound speed. */
float soundspeed;
......
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