Commit 680340fa authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Improved documentation of the external potentials.

parent 21a4773f
......@@ -85,7 +85,7 @@ IsothermalPotential:
DiscPatchPotential:
surface_density: 10. # Surface density of the disc (internal units)
scale_height: 100. # Scale height of the disc (internal units)
z_disc: 200. # Disc height (internal units)
z_disc: 200. # Position of the disc along the z-axis (internal units)
timestep_mult: 0.03 # Dimensionless pre-factor for the time-step condition
growth_time: 5. # (Optional) Time for the disc to grow to its final size (multiple of the dynamical time)
......
......@@ -22,17 +22,14 @@
#include <float.h>
#include "minmax.h"
#include "physical_constants.h"
/**
* @brief Computes the gravity time-step of a given particle due to self-gravity
*
* @param phys_const The physical constants in internal units.
* @param gp Pointer to the g-particle data.
*/
__attribute__((always_inline)) INLINE static float
gravity_compute_timestep_self(const struct phys_const* const phys_const,
const struct gpart* const gp) {
gravity_compute_timestep_self(const struct gpart* const gp) {
const float ac2 = gp->a_grav[0] * gp->a_grav[0] +
gp->a_grav[1] * gp->a_grav[1] +
......
......@@ -35,14 +35,29 @@
#include "physical_constants.h"
#include "units.h"
/* External Potential Properties */
/**
* @brief External Potential Properties - Disc patch case
*
* See Creasey, Theuns & Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948
*/
struct external_potential {
/*! Surface density of the disc */
double surface_density;
/*! Disc scale-height */
double scale_height;
/*! Position of the disc along the z-axis */
double z_disc;
/*! Dynamical time of the system */
double dynamical_time;
/*! Time over which to grow the disk in units of the dynamical time */
double growth_time;
/*! Time-step condition pre-factor */
double timestep_mult;
};
......@@ -70,7 +85,7 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep(
const float dz = fabs(g->x[2] - potential->z_disc);
/* vertical cceleration */
const float z_accel = 2 * M_PI * phys_const->const_newton_G *
const float z_accel = 2.f * M_PI * phys_const->const_newton_G *
potential->surface_density *
tanh(dz / potential->scale_height);
......@@ -92,7 +107,7 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep(
float dt3 = FLT_MAX;
if (abs(g->v_full[2]) > 0) {
const float dz_accel_over_dt =
2 * M_PI * phys_const->const_newton_G * potential->surface_density /
2.f * M_PI * phys_const->const_newton_G * potential->surface_density /
potential->scale_height / cosh(dz / potential->scale_height) /
cosh(dz / potential->scale_height) * fabs(g->v_full[2]);
......@@ -118,21 +133,20 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
double time, const struct external_potential* restrict potential,
const struct phys_const* restrict phys_const, struct gpart* restrict g) {
const float G_newton = phys_const->const_newton_G;
const float dz = g->x[2] - potential->z_disc;
const float t_dyn = potential->dynamical_time;
float reduction_factor = 1.;
float reduction_factor = 1.f;
if (time < potential->growth_time * t_dyn)
reduction_factor = time / (potential->growth_time * t_dyn);
const float z_accel = reduction_factor * 2 * G_newton * M_PI *
/* Accelerations. Note that they are multiplied by G later on */
const float z_accel = reduction_factor * 2.f * M_PI *
potential->surface_density *
tanh(fabs(dz) / potential->scale_height);
/* Accelerations. Note that they are multiplied by G later on */
if (dz > 0) g->a_grav[2] -= z_accel / G_newton;
if (dz < 0) g->a_grav[2] += z_accel / G_newton;
if (dz > 0) g->a_grav[2] -= z_accel;
if (dz < 0) g->a_grav[2] += z_accel;
}
/**
......@@ -149,16 +163,16 @@ static INLINE void potential_init_backend(
const struct phys_const* phys_const, const struct UnitSystem* us,
struct external_potential* potential) {
potential->surface_density =
parser_get_param_double(parameter_file, "DiscPatchPotential:surface_density");
potential->scale_height =
parser_get_param_double(parameter_file, "DiscPatchPotential:scale_height");
potential->surface_density = parser_get_param_double(
parameter_file, "DiscPatchPotential:surface_density");
potential->scale_height = parser_get_param_double(
parameter_file, "DiscPatchPotential:scale_height");
potential->z_disc =
parser_get_param_double(parameter_file, "DiscPatchPotential:z_disc");
potential->timestep_mult =
parser_get_param_double(parameter_file, "DiscPatchPotential:timestep_mult");
potential->growth_time =
parser_get_opt_param_double(parameter_file, "DiscPatchPotential:growth_time", 0.);
potential->timestep_mult = parser_get_param_double(
parameter_file, "DiscPatchPotential:timestep_mult");
potential->growth_time = parser_get_opt_param_double(
parameter_file, "DiscPatchPotential:growth_time", 0.);
potential->dynamical_time =
sqrt(potential->scale_height /
(phys_const->const_newton_G * potential->surface_density));
......@@ -179,7 +193,7 @@ static INLINE void potential_print_backend(
potential->timestep_mult);
if (potential->growth_time > 0.)
message("Disc will grow for %f dynamiccal times.", potential->growth_time);
message("Disc will grow for %f dynamical times.", potential->growth_time);
}
#endif /* SWIFT_DISC_PATCH_H */
......@@ -33,12 +33,22 @@
#include "physical_constants.h"
#include "units.h"
/* External Potential Properties */
/**
* @brief External Potential Properties - Isothermal sphere case
*/
struct external_potential {
/*! Position of the centre of potential */
double x, y, z;
/*! Rotation velocity */
double vrot;
float timestep_mult;
/*! Square of vrot divided by G \f$ \frac{v_{rot}^2}{G} \f$ */
double vrot2_over_G;
/*! Time-step condition pre-factor */
double timestep_mult;
};
/**
......@@ -65,11 +75,11 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep(
const double vrot = potential->vrot;
const float dota_x =
vrot * vrot * rinv2 * (g->v_full[0] - 2 * drdv * dx * rinv2);
vrot * vrot * rinv2 * (g->v_full[0] - 2.f * drdv * dx * rinv2);
const float dota_y =
vrot * vrot * rinv2 * (g->v_full[1] - 2 * drdv * dy * rinv2);
vrot * vrot * rinv2 * (g->v_full[1] - 2.f * drdv * dy * rinv2);
const float dota_z =
vrot * vrot * rinv2 * (g->v_full[2] - 2 * drdv * dz * rinv2);
vrot * vrot * rinv2 * (g->v_full[2] - 2.f * 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];
......@@ -92,14 +102,12 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
double time, const struct external_potential* potential,
const struct phys_const* const phys_const, struct gpart* g) {
const double G_newton = phys_const->const_newton_G;
const float dx = g->x[0] - potential->x;
const float dy = g->x[1] - potential->y;
const float dz = g->x[2] - potential->z;
const float rinv2 = 1.f / (dx * dx + dy * dy + dz * dz);
const double vrot = potential->vrot;
const double term = -vrot * vrot * rinv2 / G_newton;
const double term = -potential->vrot2_over_G * rinv2;
g->a_grav[0] += term * dx;
g->a_grav[1] += term * dy;
......@@ -130,6 +138,9 @@ static INLINE void potential_init_backend(
parser_get_param_double(parameter_file, "IsothermalPotential:vrot");
potential->timestep_mult = parser_get_param_float(
parameter_file, "IsothermalPotential:timestep_mult");
potential->vrot2_over_G =
potential->vrot * potential->vrot / phys_const->const_newton_G;
}
/**
......
......@@ -34,12 +34,17 @@
#include "units.h"
/**
* @brief External Potential Properties
* @brief External Potential Properties - Point mass case
*/
struct external_potential {
/*! Position of the point mass */
double x, y, z;
/*! Mass */
double mass;
/*! Time-step condition pre-factor */
float timestep_mult;
};
......
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