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

More consistent notation for the mesh-smoothing scale r_s

parent 307a7b40
Hide whitespace changes
Inline Side-by-side
src/engine.c
View file @ f3e310b5
......@@ -6268,11 +6268,11 @@ void engine_recompute_displacement_constraint(struct engine *e) {
const float vel_norm_b = vel_norm[0] + vel_norm[4];
/* Mesh forces smoothing scale */
float a_smooth;
float r_s;
if ((e->policy & engine_policy_self_gravity) && e->s->periodic == 1)
a_smooth = e->mesh->a_smooth;
r_s = e->mesh->r_s;
else
a_smooth = FLT_MAX;
r_s = FLT_MAX;
float dt_dm = FLT_MAX, dt_b = FLT_MAX;
......@@ -6289,7 +6289,7 @@ void engine_recompute_displacement_constraint(struct engine *e) {
const float rms_vel_dm = vel_norm_dm / N_dm;
/* Time-step based on maximum displacement */
dt_dm = a * a * min(a_smooth, d_dm) / sqrtf(rms_vel_dm);
dt_dm = a * a * min(r_s, d_dm) / sqrtf(rms_vel_dm);
}
/* Baryon case */
......@@ -6305,7 +6305,7 @@ void engine_recompute_displacement_constraint(struct engine *e) {
const float rms_vel_b = vel_norm_b / N_b;
/* Time-step based on maximum displacement */
dt_b = a * a * min(a_smooth, d_b) / sqrtf(rms_vel_b);
dt_b = a * a * min(r_s, d_b) / sqrtf(rms_vel_b);
}
/* Use the minimum */
......
src/gravity_derivatives.h
View file @ f3e310b5
......@@ -132,7 +132,7 @@ struct potential_derivatives_M2P {
* @param eps Softening length.
* @param eps_inv Inverse of softening length.
* @param periodic Is the calculation periodic ?
* @param rs_inv
* @param r_s_inv Inverse of the long-range gravity mesh smoothing length.
* @param pot (return) The structure containing all the derivatives.
*/
__attribute__((always_inline)) INLINE static void
......@@ -140,7 +140,7 @@ compute_potential_derivatives_M2L(const float r_x, const float r_y,
const float r_z, const float r2,
const float r_inv, const float eps,
const float eps_inv, const int periodic,
const float rs_inv,
const float r_s_inv,
struct potential_derivatives_M2L *pot) {
float Dt_1;
......@@ -189,8 +189,8 @@ compute_potential_derivatives_M2L(const float r_x, const float r_y,
/* Get the derivatives of the truncated potential */
const float r = r2 * r_inv;
struct truncated_derivatives derivs;
kernel_long_grav_derivatives(r, rs_inv, &derivs);
struct chi_derivatives derivs;
kernel_long_grav_derivatives(r, r_s_inv, &derivs);
Dt_1 = derivs.chi_0 * r_inv;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
......@@ -388,7 +388,7 @@ compute_potential_derivatives_M2L(const float r_x, const float r_y,
* @param eps Softening length.
* @param eps_inv Inverse of softening length.
* @param periodic Is the calculation using periodic BCs?
* @param rs_inv The inverse of the gravity mesh-smoothing scale.
* @param r_s_inv The inverse of the gravity mesh-smoothing scale.
* @param pot (return) The structure containing all the derivatives.
*/
__attribute__((always_inline)) INLINE static void
......@@ -396,7 +396,7 @@ compute_potential_derivatives_M2P(const float r_x, const float r_y,
const float r_z, const float r2,
const float r_inv, const float eps,
const float eps_inv, const int periodic,
const float rs_inv,
const float r_s_inv,
struct potential_derivatives_M2P *pot) {
float Dt_1;
......@@ -421,8 +421,8 @@ compute_potential_derivatives_M2P(const float r_x, const float r_y,
/* Get the derivatives of the truncated potential */
const float r = r2 * r_inv;
struct truncated_derivatives d;
kernel_long_grav_derivatives(r, rs_inv, &d);
struct chi_derivatives d;
kernel_long_grav_derivatives(r, r_s_inv, &d);
const float r_inv2 = r_inv * r_inv;
const float r_inv3 = r_inv2 * r_inv;
......
src/kernel_long_gravity.h
View file @ f3e310b5
......@@ -33,41 +33,63 @@
#define GADGET2_LONG_RANGE_CORRECTION
struct truncated_derivatives {
/**
* @brief Derivatives of the long-range truncation function \f$\chi(r,r_s)\f$ up
* to 5th order.
*/
struct chi_derivatives {
/*! 0th order derivative \f$\chi(r,r_s)\f$ */
float chi_0;
/*! 1st order derivative \f$\partial_{r}\chi(r,r_s)\f$ */
float chi_1;
/*! 2nd order derivative \f$\partial_{rr}\chi(r,r_s)\f$ */
float chi_2;
/*! 3rd order derivative \f$\partial_{rrr}\chi(r,r_s)\f$ */
float chi_3;
/*! 4th order derivative \f$\partial_{rrrr}\chi(r,r_s)\f$ */
float chi_4;
/*! 5th order derivative \f$\partial_{rrrrr}\chi(r,r_s)\f$ */
float chi_5;
};
/**
* @brief Compute the derivatives of the long-range truncation function
* \f$\chi(r,r_s)\f$ up to 5th order.
*
* @param r The distance.
* @param r_s_inv The inverse of the long-range gravity mesh scale.
* @param derivs (return) The computed #chi_derivatives.
*/
__attribute__((always_inline)) INLINE static void kernel_long_grav_derivatives(
const float r, const float rs_inv,
struct truncated_derivatives *const derivs) {
const float r, const float r_s_inv, struct chi_derivatives *const derivs) {
#ifdef GADGET2_LONG_RANGE_CORRECTION
/* Powers of u=r/2r_s */
const float u = 0.5f * r * rs_inv;
const float u = 0.5f * r * r_s_inv;
const float u2 = u * u;
const float u3 = u2 * u;
const float u4 = u3 * u;
/* Powers of (1/r_s) */
const float rs_inv2 = rs_inv * rs_inv;
const float rs_inv3 = rs_inv2 * rs_inv;
const float rs_inv4 = rs_inv3 * rs_inv;
const float rs_inv5 = rs_inv4 * rs_inv;
const float r_s_inv2 = r_s_inv * r_s_inv;
const float r_s_inv3 = r_s_inv2 * r_s_inv;
const float r_s_inv4 = r_s_inv3 * r_s_inv;
const float r_s_inv5 = r_s_inv4 * r_s_inv;
/* Derivatives of \chi */
derivs->chi_0 = erfcf(u);
derivs->chi_1 = -rs_inv;
derivs->chi_2 = rs_inv2 * u;
derivs->chi_3 = -rs_inv3 * (u2 - 0.5f);
derivs->chi_4 = rs_inv4 * (u3 - 1.5f * u);
derivs->chi_5 = -rs_inv5 * (u4 - 3.f * u2 + 0.75f);
derivs->chi_1 = -r_s_inv;
derivs->chi_2 = r_s_inv2 * u;
derivs->chi_3 = -r_s_inv3 * (u2 - 0.5f);
derivs->chi_4 = r_s_inv4 * (u3 - 1.5f * u);
derivs->chi_5 = -r_s_inv5 * (u4 - 3.f * u2 + 0.75f);
const float one_over_sqrt_pi = ((float)(M_2_SQRTPI * 0.5));
const float common_factor = one_over_sqrt_pi * expf(-u2);
......@@ -83,7 +105,7 @@ __attribute__((always_inline)) INLINE static void kernel_long_grav_derivatives(
/* Powers of 2/r_s */
const float c0 = 1.f;
const float c1 = 2.f * rs_inv;
const float c1 = 2.f * r_s_inv;
const float c2 = c1 * c1;
const float c3 = c2 * c1;
const float c4 = c3 * c1;
......@@ -93,7 +115,7 @@ __attribute__((always_inline)) INLINE static void kernel_long_grav_derivatives(
const float x = c1 * r;
/* e^(2r / r_s) */
const float exp_x = expf(x); // good_approx_expf(x);
const float exp_x = good_approx_expf(x);
/* 1 / alpha(w) */
const float a_inv = 1.f + exp_x;
......@@ -136,7 +158,7 @@ __attribute__((always_inline)) INLINE static void kernel_long_grav_pot_eval(
#else
const float x = 2.f * u;
const float exp_x = expf(x); // good_approx_expf(x);
const float exp_x = good_approx_expf(x);
const float alpha = 1.f / (1.f + exp_x);
/* We want 2 - 2 exp(x) * alpha */
......@@ -169,7 +191,7 @@ __attribute__((always_inline)) INLINE static void kernel_long_grav_force_eval(
#else
const float x = 2.f * u;
const float exp_x = expf(x); // good_approx_expf(x);
const float exp_x = good_approx_expf(x);
const float alpha = 1.f / (1.f + exp_x);
/* We want 2*(x*alpha - x*alpha^2 - exp(x)*alpha + 1) */
......@@ -177,13 +199,6 @@ __attribute__((always_inline)) INLINE static void kernel_long_grav_force_eval(
*W = *W * x - exp_x;
*W = *W * alpha + 1.f;
*W *= 2.f;
/* const float arg = 2.f * u; */
/* const float exp_arg = good_approx_expf(arg); */
/* const float term = 1.f / (1.f + exp_arg); */
/* *W = arg * exp_arg * term * term - exp_arg * term + 1.f; */
/* *W *= 2.f; */
#endif
}
......
src/mesh_gravity.c
View file @ f3e310b5
......@@ -327,13 +327,13 @@ void pm_mesh_compute_potential(struct pm_mesh* mesh, const struct engine* e) {
#ifdef HAVE_FFTW
const struct space* s = e->s;
const double a_smooth = mesh->a_smooth;
const double r_s = mesh->r_s;
const double box_size = s->dim[0];
const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]};
const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
if (cdim[0] != cdim[1] || cdim[0] != cdim[2]) error("Non-square mesh");
// if (a_smooth <= 0.) error("Invalid value of a_smooth");
if (r_s <= 0.) error("Invalid value of a_smooth");
/* Some useful constants */
const int N = mesh->N;
......@@ -389,8 +389,7 @@ void pm_mesh_compute_potential(struct pm_mesh* mesh, const struct engine* e) {
/* Some common factors */
const double green_fac = -1. / (M_PI * box_size);
const double a_smooth2 =
4. * M_PI * M_PI * a_smooth * a_smooth / (box_size * box_size);
const double a_smooth2 = 4. * M_PI * M_PI * r_s * r_s / (box_size * box_size);
const double k_fac = M_PI / (double)N;
/* Now de-convolve the CIC kernel and apply the Green function */
......@@ -516,12 +515,12 @@ void pm_mesh_init(struct pm_mesh* mesh, const struct gravity_props* props,
#ifdef HAVE_FFTW
/* Initi the mesh size */
const int N = props->mesh_size;
mesh->N = N;
mesh->box_size = box_size;
mesh->cell_fac = N / box_size;
mesh->a_smooth = props->a_smooth * box_size / N;
mesh->r_s = props->a_smooth * box_size / N;
mesh->r_s_inv = 1. / mesh->r_s;
/* Allocate the memory for the combined density and potential array */
mesh->potential = (double*)fftw_malloc(sizeof(double) * N * N * N);
......
src/mesh_gravity.h
View file @ f3e310b5
......@@ -44,7 +44,10 @@ struct pm_mesh {
double box_size;
/*! Scale over which we smooth the forces */
double a_smooth;
double r_s;
/*! Inverse of the scale over which we smooth the forces */
double r_s_inv;
/*! Potential field */
double *potential;
......
src/runner_doiact_grav.h
View file @ f3e310b5
......@@ -175,8 +175,7 @@ static INLINE void runner_dopair_grav_mm(struct runner *r,
const int periodic = s->periodic;
const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
const struct gravity_props *props = e->gravity_properties;
const float rlr = e->mesh->a_smooth;
const float rlr_inv = 1.f / rlr;
const float r_s_inv = e->mesh->r_s_inv;
TIMER_TIC;
......@@ -209,7 +208,7 @@ static INLINE void runner_dopair_grav_mm(struct runner *r,
/* Let's interact at this level */
if (0)
gravity_M2L(&ci->multipole->pot, multi_j, ci->multipole->CoM,
cj->multipole->CoM, props, periodic, dim, rlr_inv);
cj->multipole->CoM, props, periodic, dim, r_s_inv);
runner_dopair_grav_pp(r, ci, cj, 0);
......@@ -323,7 +322,7 @@ static INLINE void runner_dopair_grav_pp_full(const struct engine *e,
}
static INLINE void runner_dopair_grav_pp_truncated(
const struct engine *e, const float rlr_inv, struct gravity_cache *ci_cache,
const struct engine *e, const float r_s_inv, struct gravity_cache *ci_cache,
struct gravity_cache *cj_cache, int gcount_i, int gcount_j,
int gcount_padded_j, struct gpart *restrict gparts_i,
struct gpart *restrict gparts_j) {
......@@ -403,7 +402,7 @@ static INLINE void runner_dopair_grav_pp_truncated(
/* Interact! */
float f_ij, pot_ij;
runner_iact_grav_pp_truncated(r2, h2_i, h_inv_i, h_inv3_i, mass_j,
rlr_inv, &f_ij, &pot_ij);
r_s_inv, &f_ij, &pot_ij);
/* if (id_i == 1000 && id_j == 900 && pjd < gcount_j) { */
/* message("--- Interacting part"); */
......@@ -540,7 +539,7 @@ static INLINE void runner_dopair_grav_pm_truncated(
struct cell *restrict cj) {
const float dim[3] = {e->s->dim[0], e->s->dim[1], e->s->dim[2]};
const float rlr_inv = 1. / e->mesh->a_smooth;
const float r_s_inv = e->mesh->r_s_inv;
TIMER_TIC;
......@@ -608,7 +607,7 @@ static INLINE void runner_dopair_grav_pm_truncated(
/* Interact! */
float f_x, f_y, f_z, pot_ij;
runner_iact_grav_pm_truncated(dx, dy, dz, r2, h_i, h_inv_i, rlr_inv,
runner_iact_grav_pm_truncated(dx, dy, dz, r2, h_i, h_inv_i, r_s_inv,
multi_j, &f_x, &f_y, &f_z, &pot_ij);
/* Store it back */
......@@ -656,9 +655,9 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
struct space *s = e->s;
const int periodic = s->periodic;
const double r_cut_min = e->gravity_properties->r_cut_min;
const double rlr = e->mesh->a_smooth;
const double min_trunc = rlr * r_cut_min;
const float rlr_inv = 1. / rlr;
const double r_s = e->mesh->r_s;
const float r_s_inv = e->mesh->r_s_inv;
const double min_trunc = r_s * r_cut_min;
/* Caches to play with */
struct gravity_cache *const ci_cache = &r->ci_gravity_cache;
......@@ -773,7 +772,7 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
if (ci_active) {
/* First the (truncated) P2P */
runner_dopair_grav_pp_truncated(e, rlr_inv, ci_cache, cj_cache,
runner_dopair_grav_pp_truncated(e, r_s_inv, ci_cache, cj_cache,
gcount_i, gcount_j, gcount_padded_j,
ci->gparts, cj->gparts);
......@@ -784,7 +783,7 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
if (cj_active && symmetric) {
/* First the (truncated) P2P */
runner_dopair_grav_pp_truncated(e, rlr_inv, cj_cache, ci_cache,
runner_dopair_grav_pp_truncated(e, r_s_inv, cj_cache, ci_cache,
gcount_j, gcount_i, gcount_padded_i,
cj->gparts, ci->gparts);
......@@ -913,6 +912,7 @@ static INLINE void runner_doself_grav_pp_full(struct runner *r,
const float mass_j = ci_cache->m[pjd];
/* Compute the pairwise (square) distance. */
/* Note: no need for periodic wrapping inside a cell */
const float dx = x_j - x_i;
const float dy = y_j - y_i;
const float dz = z_j - z_i;
......@@ -971,8 +971,7 @@ static INLINE void runner_doself_grav_pp_truncated(struct runner *r,
/* Some constants */
const struct engine *const e = r->e;
const double rlr = e->mesh->a_smooth;
const float rlr_inv = 1. / rlr;
const float r_s_inv = e->mesh->r_s_inv;
/* Caches to play with */
struct gravity_cache *const ci_cache = &r->ci_gravity_cache;
......@@ -1042,13 +1041,10 @@ static INLINE void runner_doself_grav_pp_truncated(struct runner *r,
const float mass_j = ci_cache->m[pjd];
/* Compute the pairwise (square) distance. */
float dx = x_j - x_i;
float dy = y_j - y_i;
float dz = z_j - z_i;
/* dx = nearestf(dx, e->s->dim[0]); */
/* dy = nearestf(dy, e->s->dim[1]); */
/* dz = nearestf(dz, e->s->dim[2]); */
/* Note: no need for periodic wrapping inside a cell */
const float dx = x_j - x_i;
const float dy = y_j - y_i;
const float dz = z_j - z_i;
const float r2 = dx * dx + dy * dy + dz * dz;
......@@ -1063,38 +1059,10 @@ static INLINE void runner_doself_grav_pp_truncated(struct runner *r,
error("gpj not drifted to current time");
#endif
/* long long id_i = e->s->parts[-gparts[pid].id_or_neg_offset].id; */
/* long long id_j = e->s->parts[-gparts[pjd].id_or_neg_offset].id; */
/* Interact! */
float f_ij, pot_ij;
runner_iact_grav_pp_truncated(r2, h2_i, h_inv_i, h_inv3_i, mass_j,
rlr_inv, &f_ij, &pot_ij);
/* if (id_i == 1000 && id_j == 901) { */
/* message("--- Interacting part"); */
/* message("id=%lld mass=%e r=%e h=%e", id_j, mass_j, sqrtf(r2), */
/* sqrtf(h2_i)); */
/* /\* message("e->s->dim[0]=%e", e->s->dim[0]); *\/ */
/* message("dx=[%e %e %e]", dx, dy, dz); */
/* message("pos_i=[%e %e %e]", x_i, y_i, z_i); */
/* message("pos_j=[%e %e %e]", x_j, y_j, z_j); */
/* message("fac=%e corr=%e fac2=%e", f1_ij, corr, pot_ij); */
/* message("a=[%e %e %e]", f_ij * dx, f_ij * dy, f_ij * dz); */
/* message("pot=%e", pot_ij); */
/* } */
/* if (0 && id_i == 1000) { */
/* message("--- Interacting part"); */
/* message("id=%lld mass=%e r=%e h=%e", id_j, mass_j, sqrtf(r2), */
/* sqrtf(h2_i)); */
/* /\* message("e->s->dim[0]=%e", e->s->dim[0]); *\/ */
/* message("dx=[%e %e %e]", dx, dy, dz); */
/* message("pos_i=[%e %e %e]", x_i, y_i, z_i); */
/* message("pos_j=[%e %e %e]", x_j, y_j, z_j); */
/* message("fac=%e corr=%e fac2=%e", f1_ij, corr, f_ij); */
/* message("a=[%e %e %e]", f_ij * dx, f_ij * dy, f_ij * dz); */
/* } */
r_s_inv, &f_ij, &pot_ij);
/* Store it back */
a_x += f_ij * dx;
......@@ -1132,9 +1100,8 @@ static INLINE void runner_doself_grav_pp(struct runner *r, struct cell *c) {
const struct engine *e = r->e;
const struct space *s = e->s;
const int periodic = s->periodic;
const double a_smooth = e->mesh->a_smooth;
const double r_cut_min = e->gravity_properties->r_cut_min;
const double min_trunc = r_cut_min * a_smooth;
const double r_s = e->mesh->r_s;
const double min_trunc = e->gravity_properties->r_cut_min * r_s;
TIMER_TIC;
......@@ -1189,9 +1156,8 @@ static INLINE void runner_dopair_grav(struct runner *r, struct cell *ci,
const int nodeID = e->nodeID;
const int periodic = s->periodic;
const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
const struct gravity_props *props = e->gravity_properties;
const double theta_crit2 = props->theta_crit2;
const double max_distance = e->mesh->a_smooth * props->r_cut_max;
const double theta_crit2 = e->gravity_properties->theta_crit2;
const double max_distance = e->mesh->r_s * e->gravity_properties->r_cut_max;
/* Anything to do here? */
if (!((cell_is_active_gravity(ci, e) && ci->nodeID == nodeID) ||
......@@ -1389,13 +1355,11 @@ static INLINE void runner_do_grav_long_range(struct runner *r, struct cell *ci,
/* Some constants */
const struct engine *e = r->e;
const struct space *s = e->s;
const struct gravity_props *props = e->gravity_properties;
const int periodic = s->periodic;
const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]};
const double theta_crit2 = props->theta_crit2;
const double max_distance = e->mesh->a_smooth * props->r_cut_max;
// const double max_distance2 = max_distance * max_distance;
const double theta_crit2 = e->gravity_properties->theta_crit2;
const double max_distance = e->mesh->r_s * e->gravity_properties->r_cut_max;
/* const float u = 1.; */
/* float W; */
......
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