Skip to content
Snippets Groups Projects
Commit 733fed35 authored by Josh Borrow's avatar Josh Borrow
Browse files

Formatting

parent da1ca8a4
No related branches found
No related tags found
1 merge request!540Add Pressure-Energy (P-U) SPH
Hide whitespace changes
Inline Side-by-side
src/hydro/PressureEnergy/hydro.h
+ 3
− 2
View file @ 733fed35
......@@ -375,7 +375,8 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
/* Re-set problematic values */
p->rho = p->mass * kernel_root * h_inv_dim;
p->pressure_bar = p->mass * p->u * hydro_gamma_minus_one * kernel_root * h_inv_dim;
p->pressure_bar =
p->mass * p->u * hydro_gamma_minus_one * kernel_root * h_inv_dim;
p->density.wcount = kernel_root * kernel_norm * h_inv_dim;
p->density.rho_dh = 0.f;
p->density.wcount_dh = 0.f;
......@@ -409,7 +410,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Compute the "grad h" term */
const float common_factor = p->h / (hydro_dimension * p->density.wcount);
const float grad_h_term =
(p->density.pressure_bar_dh * common_factor / hydro_gamma_minus_one) /
(p->density.pressure_bar_dh * common_factor / hydro_gamma_minus_one) /
(1 + common_factor * p->density.wcount_dh);
/* Update variables. */
......
src/hydro/PressureEnergy/hydro_debug.h
+ 1
− 2
View file @ 733fed35
......@@ -45,8 +45,7 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
p->u, p->u_dt, p->force.v_sig, hydro_get_comoving_pressure(p), p->h,
p->force.h_dt, (int)p->density.wcount, p->mass, p->density.rho_dh, p->rho,
p->density.pressure_bar_dh, p->pressure_bar,
p->time_bin);
p->density.pressure_bar_dh, p->pressure_bar, p->time_bin);
}
#endif /* SWIFT_MINIMAL_HYDRO_DEBUG_H */
src/hydro/PressureEnergy/hydro_iact.h
+ 28
− 22
View file @ 733fed35
......@@ -48,7 +48,7 @@
* @param H Current Hubble parameter.
*/
__attribute__((always_inline)) INLINE static void runner_iact_density(
float r2, const float *dx, float hi, float hj, struct part* pi,
float r2, const float* dx, float hi, float hj, struct part* pi,
struct part* pj, float a, float H) {
float wi, wj, wi_dx, wj_dx;
......@@ -67,7 +67,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
pi->rho += mj * wi;
pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
pi->pressure_bar += mj * wi * pj->u;
pi->density.pressure_bar_dh -= mj * pj->u * (hydro_dimension * wi + ui * wi_dx);
pi->density.pressure_bar_dh -=
mj * pj->u * (hydro_dimension * wi + ui * wi_dx);
pi->density.wcount += wi;
pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
......@@ -79,7 +80,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
pj->rho += mi * wj;
pj->density.rho_dh -= mi * (hydro_dimension * wj + uj * wj_dx);
pj->pressure_bar += mi * wj * pi->u;
pj->density.pressure_bar_dh -= mi * pi->u * (hydro_dimension * wj + uj * wj_dx);
pj->density.pressure_bar_dh -=
mi * pi->u * (hydro_dimension * wj + uj * wj_dx);
pj->density.wcount += wj;
pj->density.wcount_dh -= (hydro_dimension * wj + uj * wj_dx);
}
......@@ -97,7 +99,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
* @param H Current Hubble parameter.
*/
__attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
float r2, const float *dx, float hi, float hj, struct part* pi,
float r2, const float* dx, float hi, float hj, struct part* pi,
const struct part* pj, float a, float H) {
float wi, wi_dx;
......@@ -115,7 +117,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
pi->rho += mj * wi;
pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
pi->pressure_bar += mj * wi * pj->u;
pi->density.pressure_bar_dh -= mj * pj->u * (hydro_dimension * wi + ui * wi_dx);
pi->density.pressure_bar_dh -=
mj * pj->u * (hydro_dimension * wi + ui * wi_dx);
pi->density.wcount += wi;
pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
}
......@@ -133,7 +136,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
* @param H Current Hubble parameter.
*/
__attribute__((always_inline)) INLINE static void runner_iact_force(
float r2, const float *dx, float hi, float hj, struct part* pi,
float r2, const float* dx, float hi, float hj, struct part* pi,
struct part* pj, float a, float H) {
/* Cosmological factors entering the EoMs */
......@@ -146,17 +149,16 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Recover some data */
const float mj = pj->mass;
const float mi = pi->mass;
const float miui = mi * pi->u;
const float mjuj = mj * pj->u;
const float rhoi = pi->rho;
const float rhoj = pj->rho;
/* Compute gradient terms */
const float f_ij = 1 - (pi->force.f / mjuj);
const float f_ji = 1 - (pj->force.f / miui);
/* Get the kernel for hi. */
const float hi_inv = 1.0f / hi;
const float hid_inv = pow_dimension_plus_one(hi_inv); /* 1/h^(d+1) */
......@@ -195,9 +197,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
/* SPH acceleration term */
const float sph_acc_term = pj->u * pi->u *
hydro_gamma_minus_one * hydro_gamma_minus_one *
((f_ij/pi->pressure_bar) * wi_dr + (f_ji/pj->pressure_bar) * wj_dr) * r_inv;
const float sph_acc_term =
pj->u * pi->u * hydro_gamma_minus_one * hydro_gamma_minus_one *
((f_ij / pi->pressure_bar) * wi_dr + (f_ji / pj->pressure_bar) * wj_dr) *
r_inv;
/* Assemble the acceleration */
const float acc = sph_acc_term + visc_acc_term;
......@@ -213,9 +216,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Get the time derivative for u. */
const float sph_du_term_i = hydro_gamma_minus_one * hydro_gamma_minus_one *
pj->u * pi->u * (f_ij / pi->pressure_bar) * wi_dr * dvdr * r_inv;
pj->u * pi->u * (f_ij / pi->pressure_bar) *
wi_dr * dvdr * r_inv;
const float sph_du_term_j = hydro_gamma_minus_one * hydro_gamma_minus_one *
pi->u * pj->u * (f_ji / pj->pressure_bar) * wj_dr * dvdr * r_inv;
pi->u * pj->u * (f_ji / pj->pressure_bar) *
wj_dr * dvdr * r_inv;
/* Viscosity term */
const float visc_du_term = 0.5f * visc_acc_term * dvdr;
......@@ -250,7 +255,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
* @param H Current Hubble parameter.
*/
__attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
float r2, const float *dx, float hi, float hj, struct part* pi,
float r2, const float* dx, float hi, float hj, struct part* pi,
const struct part* pj, float a, float H) {
/* Cosmological factors entering the EoMs */
......@@ -264,17 +269,16 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
// const float mi = pi->mass;
const float mj = pj->mass;
const float mi = pi->mass;
const float miui = mi * pi->u;
const float mjuj = mj * pj->u;
const float rhoi = pi->rho;
const float rhoj = pj->rho;
/* Compute gradient terms */
const float f_ij = 1 - (pi->force.f / mjuj);
const float f_ji = 1 - (pj->force.f / miui);
/* Get the kernel for hi. */
const float hi_inv = 1.0f / hi;
const float hid_inv = pow_dimension_plus_one(hi_inv); /* 1/h^(d+1) */
......@@ -313,9 +317,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
/* SPH acceleration term */
const float sph_acc_term = pj->u * pi->u *
hydro_gamma_minus_one * hydro_gamma_minus_one *
((f_ij/pi->pressure_bar) * wi_dr + (f_ji/pj->pressure_bar) * wj_dr) * r_inv;
const float sph_acc_term =
pj->u * pi->u * hydro_gamma_minus_one * hydro_gamma_minus_one *
((f_ij / pi->pressure_bar) * wi_dr + (f_ji / pj->pressure_bar) * wj_dr) *
r_inv;
/* Assemble the acceleration */
const float acc = sph_acc_term + visc_acc_term;
......@@ -327,7 +332,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Get the time derivative for u. */
const float sph_du_term_i = hydro_gamma_minus_one * hydro_gamma_minus_one *
pj->u * pi->u * (f_ij / pi->pressure_bar) * wi_dr * dvdr * r_inv;
pj->u * pi->u * (f_ij / pi->pressure_bar) *
wi_dr * dvdr * r_inv;
/* Viscosity term */
const float visc_du_term = 0.5f * visc_acc_term * dvdr;
......
src/hydro/PressureEnergy/hydro_io.h
+ 12
− 13
View file @ 733fed35
......@@ -69,13 +69,13 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
}
void convert_u(const struct engine* e, const struct part* p,
const struct xpart* xp, float* ret) {
const struct xpart* xp, float* ret) {
ret[0] = hydro_get_comoving_internal_energy(p);
}
void convert_P(const struct engine* e, const struct part* p,
const struct xpart* xp, float* ret) {
void convert_P(const struct engine* e, const struct part* p,
const struct xpart* xp, float* ret) {
ret[0] = hydro_get_comoving_pressure(p);
}
......@@ -135,31 +135,30 @@ void convert_part_vel(const struct engine* e, const struct part* p,
* @param list The list of i/o properties to write.
* @param num_fields The number of i/o fields to write.
*/
void hydro_write_particles(const struct part* parts,
const struct xpart* xparts,
struct io_props* list, int* num_fields) {
void hydro_write_particles(const struct part* parts, const struct xpart* xparts,
struct io_props* list, int* num_fields) {
*num_fields = 8;
/* List what we want to write */
list[0] = io_make_output_field_convert_part("Coordinates", DOUBLE, 3,
UNIT_CONV_LENGTH, parts, xparts,
convert_part_pos);
UNIT_CONV_LENGTH, parts, xparts,
convert_part_pos);
list[1] = io_make_output_field_convert_part(
"Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, xparts, convert_part_vel);
list[2] =
io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts, mass);
list[3] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
parts, h);
list[4] = io_make_output_field_convert_part(
"InternalEnergy", FLOAT, 1, UNIT_CONV_ENERGY_PER_UNIT_MASS, parts,
xparts, convert_u);
list[4] = io_make_output_field_convert_part("InternalEnergy", FLOAT, 1,
UNIT_CONV_ENERGY_PER_UNIT_MASS,
parts, xparts, convert_u);
list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
UNIT_CONV_NO_UNITS, parts, id);
list[6] =
io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts, rho);
list[7] = io_make_output_field(
"Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, pressure_bar);
list[7] = io_make_output_field("Pressure", FLOAT, 1, UNIT_CONV_PRESSURE,
parts, pressure_bar);
}
/**
......
src/hydro/PressureEnergy/hydro_part.h
+ 3
− 3
View file @ 733fed35
......@@ -101,7 +101,7 @@ struct part {
/*! Particle density. */
float rho;
/*! Particle pressure (weighted) */
float pressure_bar;
......@@ -126,9 +126,9 @@ struct part {
/*! Derivative of density with respect to h */
float rho_dh;
/*! Derivative of the weighted pressure with respect to h */
/*! Derivative of the weighted pressure with respect to h */
float pressure_bar_dh;
/*! Particle velocity curl. */
float rot_v[3];
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment