Commit 23efe7cd authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Moved the constants related to self-gravity out of const.h

parent 72127be9
......@@ -41,9 +41,6 @@
/* Self gravity stuff. */
#define const_gravity_multipole_order 1
#define const_gravity_a_smooth FLT_MAX
#define const_gravity_r_cut 4.5f
#define const_gravity_eta 0.025f
/* Type of gradients to use (GIZMO_SPH only) */
/* If no option is chosen, no gradients are used (first order scheme) */
......
......@@ -21,15 +21,18 @@
#define SWIFT_DEFAULT_GRAVITY_H
#include <float.h>
#include "gravity_properties.h"
#include "minmax.h"
/**
* @brief Computes the gravity time-step of a given particle due to self-gravity
*
* @param gp Pointer to the g-particle data.
* @param grav_props Constants used in the gravity scheme.
*/
__attribute__((always_inline)) INLINE static float
gravity_compute_timestep_self(const struct gpart* const gp) {
gravity_compute_timestep_self(const struct gpart* const gp,
const struct gravity_props* restrict grav_props) {
const float ac2 = gp->a_grav[0] * gp->a_grav[0] +
gp->a_grav[1] * gp->a_grav[1] +
......@@ -37,7 +40,7 @@ gravity_compute_timestep_self(const struct gpart* const gp) {
const float ac = (ac2 > 0.f) ? sqrtf(ac2) : FLT_MIN;
const float dt = sqrtf(2.f * const_gravity_eta * gp->epsilon / ac);
const float dt = sqrtf(2.f * grav_props->eta * gp->epsilon / ac);
return dt;
}
......
......@@ -69,7 +69,8 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pm(
const int gcount = ci->gcount;
struct gpart *restrict gparts = ci->gparts;
const struct multipole *multi = cj->multipole;
const float rlr_inv = 1. / (const_gravity_a_smooth * ci->super->width[0]);
const float a_smooth = e->gravity_properties->a_smooth;
const float rlr_inv = 1. / (a_smooth * ci->super->width[0]);
TIMER_TIC;
......@@ -138,7 +139,8 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pp(
const int gcount_j = cj->gcount;
struct gpart *restrict gparts_i = ci->gparts;
struct gpart *restrict gparts_j = cj->gparts;
const float rlr_inv = 1. / (const_gravity_a_smooth * ci->super->width[0]);
const float a_smooth = e->gravity_properties->a_smooth;
const float rlr_inv = 1. / (a_smooth * ci->super->width[0]);
TIMER_TIC;
......@@ -244,7 +246,8 @@ __attribute__((always_inline)) INLINE static void runner_doself_grav_pp(
const struct engine *e = r->e;
const int gcount = c->gcount;
struct gpart *restrict gparts = c->gparts;
const float rlr_inv = 1. / (const_gravity_a_smooth * c->super->width[0]);
const float a_smooth = e->gravity_properties->a_smooth;
const float rlr_inv = 1. / (a_smooth * c->super->width[0]);
TIMER_TIC;
......
......@@ -78,7 +78,8 @@ __attribute__((always_inline)) INLINE static integertime_t get_gpart_timestep(
e->physical_constants, gp));
if (e->policy & engine_policy_self_gravity)
new_dt = min(new_dt, gravity_compute_timestep_self(gp));
new_dt =
min(new_dt, gravity_compute_timestep_self(gp, e->gravity_properties));
/* Limit timestep within the allowed range */
new_dt = min(new_dt, e->dt_max);
......@@ -121,7 +122,8 @@ __attribute__((always_inline)) INLINE static integertime_t get_part_timestep(
e->physical_constants, p->gpart));
if (e->policy & engine_policy_self_gravity)
new_dt_grav = min(new_dt_grav, gravity_compute_timestep_self(p->gpart));
new_dt_grav = min(new_dt_grav, gravity_compute_timestep_self(
p->gpart, e->gravity_properties));
}
/* Final time-step is minimum of hydro and gravity */
......@@ -163,7 +165,8 @@ __attribute__((always_inline)) INLINE static integertime_t get_spart_timestep(
e->physical_constants, sp->gpart));
if (e->policy & engine_policy_self_gravity)
new_dt = min(new_dt, gravity_compute_timestep_self(sp->gpart));
new_dt = min(new_dt, gravity_compute_timestep_self(sp->gpart,
e->gravity_properties));
/* Limit timestep within the allowed range */
new_dt = min(new_dt, e->dt_max);
......
......@@ -732,19 +732,24 @@ int compare_particles(struct part a, struct part b, double threshold) {
#endif
}
/** @brief Computes the forces between all g-particles using the N^2 algorithm
/**
* @brief Computes the forces between all g-particles using the N^2 algorithm
*
* Overwrites the accelerations of the gparts with the values.
* Do not use for actual runs.
*
* @brief gparts The array of particles.
* @brief gcount The number of particles.
* @brief constants Physical constants in internal units.
* @brief gravity_properties Constants governing the gravity scheme.
*/
void gravity_n2(struct gpart *gparts, const int gcount,
const struct phys_const *constants, float rlr) {
const struct phys_const *constants,
const struct gravity_props *gravity_properties, float rlr) {
const float rlr_inv = 1. / rlr;
const float max_d = const_gravity_r_cut * rlr;
const float r_cut = gravity_properties->r_cut;
const float max_d = r_cut * rlr;
const float max_d2 = max_d * max_d;
message("rlr_inv= %f", rlr_inv);
......
......@@ -23,6 +23,7 @@
#define SWIFT_TOOL_H
#include "cell.h"
#include "gravity_properties.h"
#include "part.h"
#include "physical_constants.h"
#include "runner.h"
......@@ -45,8 +46,8 @@ void pairs_n2(double *dim, struct part *restrict parts, int N, int periodic);
double random_uniform(double a, double b);
void shuffle_particles(struct part *parts, const int count);
void gravity_n2(struct gpart *gparts, const int gcount,
const struct phys_const *constants, float rlr);
const struct phys_const *constants,
const struct gravity_props *gravity_properties, float rlr);
int compare_values(double a, double b, double threshold, double *absDiff,
double *absSum, double *relDiff);
int compare_particles(struct part a, struct part b, double threshold);
......
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