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

Made the potential file compile by removing duplicated lines.

parent a80b68ca
......@@ -54,6 +54,7 @@ struct external_potential {
float timestep_mult;
} isothermal_potential;
#endif
#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
struct {
double surface_density;
......@@ -65,8 +66,11 @@ struct external_potential {
#endif
};
/* ------------------------------------------------------------------------- */
/* external potential: disk_patch */
#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
/**
* @brief Computes the time-step due to the acceleration from a hydrostatic disk
* (Creasy, Theuns & Bower, 2013)
......@@ -122,12 +126,13 @@ external_gravity_disk_patch_timestep(const struct external_potential* potential,
return potential->disk_patch_potential.timestep_mult * dt;
}
/**
* @brief Computes the gravitational acceleration from a hydrostatic disk
* (Creasy, Theuns & Bower, 2013)
*
* @param phys_cont The physical constants in internal units.
* @param gp Pointer to the g-particle data.
* @param g Pointer to the g-particle data.
*/
__attribute__((always_inline)) INLINE static void
external_gravity_disk_patch_potential(
......@@ -169,6 +174,8 @@ external_gravity_disk_patch_potential(
}
#endif /* EXTERNAL_POTENTIAL_DISK_PATCH */
/* ------------------------------------------------------------------------- */
#ifdef EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
/**
......@@ -182,147 +189,136 @@ __attribute__((always_inline)) INLINE static float
external_gravity_isothermalpotential_timestep(
const struct external_potential* potential,
const struct phys_const* const phys_const, const struct gpart* const g) {
const struct gpart* const g) {
const float dx = g->x[0] - potential->isothermal_potential.x;
const float dy = g->x[1] - potential->isothermal_potential.y;
const float dz = g->x[2] - potential->isothermal_potential.z;
const float rinv2 = 1.f / (dx * dx + dy * dy + dz * dz);
const float drdv =
dx * (g->v_full[0]) + dy * (g->v_full[1]) + dz * (g->v_full[2]);
const double vrot = potential->isothermal_potential.vrot;
const float dota_x =
vrot * vrot * rinv2 * (g->v_full[0] - 2 * drdv * dx * rinv2);
const float dota_y =
vrot * vrot * rinv2 * (g->v_full[1] - 2 * drdv * dy * rinv2);
const float dota_z =
vrot * vrot * rinv2 * (g->v_full[2] - 2 * drdv * dz * rinv2);
const float dota_2 =
dota_x * dota_x + dota_y * dota_y + dota_z * dota_z;
const float a_2 = g->a_grav[0] * g->a_grav[0] +
g->a_grav[1] * g->a_grav[1] +
g->a_grav[2] * g->a_grav[2];
return potential->isothermal_potential.timestep_mult *
sqrtf(a_2 / dota_2);
}
/**
* @brief Computes the gravitational acceleration of a particle due to a
* point
* mass
*
* Note that the accelerations are multiplied by Newton's G constant
* later on.
*
* @param potential The #external_potential used in the run.
* @param phys_const The physical constants in internal units.
* @param g Pointer to the g-particle data.
*/
__attribute__((always_inline)) INLINE static void
external_gravity_isothermalpotential(
const struct external_potential* potential,
const struct phys_const* const phys_const, struct gpart* g) {
__attribute__((always_inline)) INLINE static void
external_gravity_isothermalpotential(
const struct external_potential* potential,
const struct phys_const* const phys_const, struct gpart* g) {
const float G_newton = phys_const->const_newton_G;
const float dx = g->x[0] - potential->isothermal_potential.x;
const float dy = g->x[1] - potential->isothermal_potential.y;
const float dz = g->x[2] - potential->isothermal_potential.z;
const float rinv2 = 1.f / (dx * dx + dy * dy + dz * dz);
const double vrot = potential->isothermal_potential.vrot;
const double term = -vrot * vrot * rinv2 / G_newton;
g->a_grav[0] += term * dx;
g->a_grav[1] += term * dy;
g->a_grav[2] += term * dz;
}
const float dx = g->x[0] - potential->isothermal_potential.x;
const float dy = g->x[1] - potential->isothermal_potential.y;
const float dz = g->x[2] - potential->isothermal_potential.z;
const float rinv2 = 1.f / (dx * dx + dy * dy + dz * dz);
const float drdv =
dx * (g->v_full[0]) + dy * (g->v_full[1]) + dz * (g->v_full[2]);
const double vrot = potential->isothermal_potential.vrot;
const float dota_x =
vrot * vrot * rinv2 * (g->v_full[0] - 2 * drdv * dx * rinv2);
const float dota_y =
vrot * vrot * rinv2 * (g->v_full[1] - 2 * drdv * dy * rinv2);
const float dota_z =
vrot * vrot * rinv2 * (g->v_full[2] - 2 * drdv * dz * rinv2);
const float dota_2 = dota_x * dota_x + dota_y * dota_y + dota_z * dota_z;
const float a_2 = g->a_grav[0] * g->a_grav[0] + g->a_grav[1] * g->a_grav[1] +
g->a_grav[2] * g->a_grav[2];
return potential->isothermal_potential.timestep_mult * sqrtf(a_2 / dota_2);
}
/**
* @brief Computes the gravitational acceleration of a particle due to a
* point
* mass
*
* Note that the accelerations are multiplied by Newton's G constant
* later on.
*
* @param potential The #external_potential used in the run.
* @param phys_const The physical constants in internal units.
* @param g Pointer to the g-particle data.
*/
__attribute__((always_inline)) INLINE static void
external_gravity_isothermalpotential(const struct external_potential* potential,
const struct phys_const* const phys_const,
struct gpart* g) {
const float G_newton = phys_const->const_newton_G;
const float dx = g->x[0] - potential->isothermal_potential.x;
const float dy = g->x[1] - potential->isothermal_potential.y;
const float dz = g->x[2] - potential->isothermal_potential.z;
const float rinv2 = 1.f / (dx * dx + dy * dy + dz * dz);
const double vrot = potential->isothermal_potential.vrot;
const double term = -vrot * vrot * rinv2 / G_newton;
g->a_grav[0] += term * dx;
g->a_grav[1] += term * dy;
g->a_grav[2] += term * dz;
}
#endif /* EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL */
/* ------------------------------------------------------------------------- */
/* Include exteral pointmass potential */
#ifdef EXTERNAL_POTENTIAL_POINTMASS
/**
* @brief Computes the time-step due to the acceleration from a point mass
*
* @param potential The properties of the externa potential.
* @param phys_const The physical constants in internal units.
* @param g Pointer to the g-particle data.
*/
__attribute__((always_inline)) INLINE static float
external_gravity_pointmass_timestep(
const struct external_potential* potential,
const struct phys_const* const phys_const,
const struct gpart* const g) {
const float G_newton = phys_const->const_newton_G;
const float dx = g->x[0] - potential->point_mass.x;
const float dy = g->x[1] - potential->point_mass.y;
const float dz = g->x[2] - potential->point_mass.z;
const float rinv = 1.f / sqrtf(dx * dx + dy * dy + dz * dz);
const float drdv = (g->x[0] - potential->point_mass.x) * (g->v_full[0]) +
(g->x[1] - potential->point_mass.y) * (g->v_full[1]) +
(g->x[2] - potential->point_mass.z) * (g->v_full[2]);
const float dota_x = G_newton * potential->point_mass.mass * rinv * rinv *
rinv *
(-g->v_full[0] + 3.f * rinv * rinv * drdv * dx);
const float dota_y = G_newton * potential->point_mass.mass * rinv * rinv *
rinv *
(-g->v_full[1] + 3.f * rinv * rinv * drdv * dy);
const float dota_z = G_newton * potential->point_mass.mass * rinv * rinv *
rinv *
(-g->v_full[2] + 3.f * rinv * rinv * drdv * dz);
const float dota_2 = dota_x * dota_x + dota_y * dota_y + dota_z * dota_z;
const float a_2 = g->a_grav[0] * g->a_grav[0] +
g->a_grav[1] * g->a_grav[1] +
g->a_grav[2] * g->a_grav[2];
return potential->point_mass.timestep_mult * sqrtf(a_2 / dota_2);
}
/**
* @brief Computes the gravitational acceleration of a particle due to a
* point
* mass
*
* Note that the accelerations are multiplied by Newton's G constant later
* on.
*
* @param potential The proerties of the external potential.
* @param phys_const The physical constants in internal units.
* @param g Pointer to the g-particle data.
*/
__attribute__((always_inline)) INLINE static void
external_gravity_pointmass(const struct external_potential* potential,
const struct phys_const* const phys_const,
struct gpart* g) {
const float dx = g->x[0] - potential->point_mass.x;
const float dy = g->x[1] - potential->point_mass.y;
const float dz = g->x[2] - potential->point_mass.z;
const float rinv = 1.f / sqrtf(dx * dx + dy * dy + dz * dz);
const float rinv3 = rinv * rinv * rinv;
g->a_grav[0] += -potential->point_mass.mass * dx * rinv3;
g->a_grav[1] += -potential->point_mass.mass * dy * rinv3;
g->a_grav[2] += -potential->point_mass.mass * dz * rinv3;
}
/**
* @brief Computes the time-step due to the acceleration from a point mass
*
* @param potential The properties of the externa potential.
* @param phys_const The physical constants in internal units.
* @param g Pointer to the g-particle data.
*/
__attribute__((always_inline)) INLINE static float
external_gravity_pointmass_timestep(const struct external_potential* potential,
const struct phys_const* const phys_const,
const struct gpart* const g) {
const float G_newton = phys_const->const_newton_G;
const float dx = g->x[0] - potential->point_mass.x;
const float dy = g->x[1] - potential->point_mass.y;
const float dz = g->x[2] - potential->point_mass.z;
const float rinv = 1.f / sqrtf(dx * dx + dy * dy + dz * dz);
const float drdv = (g->x[0] - potential->point_mass.x) * (g->v_full[0]) +
(g->x[1] - potential->point_mass.y) * (g->v_full[1]) +
(g->x[2] - potential->point_mass.z) * (g->v_full[2]);
const float dota_x = G_newton * potential->point_mass.mass * rinv * rinv *
rinv * (-g->v_full[0] + 3.f * rinv * rinv * drdv * dx);
const float dota_y = G_newton * potential->point_mass.mass * rinv * rinv *
rinv * (-g->v_full[1] + 3.f * rinv * rinv * drdv * dy);
const float dota_z = G_newton * potential->point_mass.mass * rinv * rinv *
rinv * (-g->v_full[2] + 3.f * rinv * rinv * drdv * dz);
const float dota_2 = dota_x * dota_x + dota_y * dota_y + dota_z * dota_z;
const float a_2 = g->a_grav[0] * g->a_grav[0] + g->a_grav[1] * g->a_grav[1] +
g->a_grav[2] * g->a_grav[2];
return potential->point_mass.timestep_mult * sqrtf(a_2 / dota_2);
}
/**
* @brief Computes the gravitational acceleration of a particle due to a
* point mass
*
* Note that the accelerations are multiplied by Newton's G constant later
* on.
*
* @param potential The proerties of the external potential.
* @param phys_const The physical constants in internal units.
* @param g Pointer to the g-particle data.
*/
__attribute__((always_inline)) INLINE static void external_gravity_pointmass(
const struct external_potential* potential,
const struct phys_const* const phys_const, struct gpart* g) {
const float dx = g->x[0] - potential->point_mass.x;
const float dy = g->x[1] - potential->point_mass.y;
const float dz = g->x[2] - potential->point_mass.z;
const float rinv = 1.f / sqrtf(dx * dx + dy * dy + dz * dz);
const float rinv3 = rinv * rinv * rinv;
g->a_grav[0] += -potential->point_mass.mass * dx * rinv3;
g->a_grav[1] += -potential->point_mass.mass * dy * rinv3;
g->a_grav[2] += -potential->point_mass.mass * dz * rinv3;
}
#endif /* EXTERNAL_POTENTIAL_POINTMASS */
/* Now, some generic functions, defined in the source file */
void potential_init(const struct swift_params* parameter_file,
const struct phys_const* const phys_const,
struct UnitSystem* us,
struct external_potential* potential);
/* ------------------------------------------------------------------------- */
/* Now, some generic functions, defined in the source file */
void potential_init(const struct swift_params* parameter_file,
const struct phys_const* const phys_const,
struct UnitSystem* us,
struct external_potential* potential);
void potential_print(const struct external_potential* potential);
void potential_print(const struct external_potential* potential);
#endif /* SWIFT_POTENTIALS_H */
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