Commit 8872a1f5 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Added a truncation distance for the disc patch potential bezond which the...

Added a truncation distance for the disc patch potential bezond which the accelerations tend towards 0. Simplified the code in that potential.
parent df6eed1e
......@@ -40,5 +40,7 @@ DiscPatchPotential:
surface_density: 10.
scale_height: 100.
z_disc: 400.
z_trunc: 300.
z_max: 380.
timestep_mult: 0.03
growth_time: 5.
......@@ -40,4 +40,6 @@ DiscPatchPotential:
surface_density: 10.
scale_height: 100.
z_disc: 400.
z_trunc: 300.
z_max: 380.
timestep_mult: 0.03
......@@ -103,7 +103,9 @@ 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. # Position of the disc along the z-axis (internal units)
z_disc: 400. # Position of the disc along the z-axis (internal units)
z_trunc: 300. # (Optional) Distance from the disc along z-axis above which the potential gets truncated.
z_max: 380. # (Optional) Distance from the disc along z-axis above which the potential is set to 0.
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)
......
......@@ -39,27 +39,54 @@
/**
* @brief External Potential Properties - Disc patch case
*
* See Creasey, Theuns & Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948
* See Creasey, Theuns & Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948.
*
* We truncate the accelerations beyond z_trunc using a 1-cos(z) function
* that smoothly brings the accelerations to 0 at z_max.
*/
struct external_potential {
/*! Surface density of the disc */
double surface_density;
/*! Surface density of the disc (sigma) */
float surface_density;
/*! Disc scale-height (b) */
float scale_height;
/*! Disc scale-height */
double scale_height;
/*! Inverse of disc scale-height (1/b) */
float scale_height_inv;
/*! Position of the disc along the z-axis */
double z_disc;
float z_disc;
/*! Position above which the accelerations get truncated */
float z_trunc;
/*! Position above which the accelerations are zero */
float z_max;
/*! The truncated transition regime */
float z_trans;
/*! Inverse of the truncated transition regime */
float z_trans_inv;
/*! Dynamical time of the system */
double dynamical_time;
float dynamical_time;
/*! Time over which to grow the disk in units of the dynamical time */
double growth_time;
/*! Time over which to grow the disk */
float growth_time;
/*! Inverse of the growth time */
float growth_time_inv;
/*! Time-step condition pre-factor */
double timestep_mult;
float timestep_mult;
/*! Constant pre-factor (2 pi G sigma) */
float norm;
/*! Constant pre-factor (2 pi sigma)*/
float norm_over_G;
};
/**
......@@ -68,6 +95,7 @@ struct external_potential {
*
* See Creasey, Theuns & Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948,
* equations 17 and 20.
* We do not use the truncated potential here.
*
* @param time The current time.
* @param potential The properties of the potential.
......@@ -81,39 +109,41 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep(
/* initilize time step to disc dynamical time */
const float dt_dyn = potential->dynamical_time;
float dt = dt_dyn;
const float b = potential->scale_height;
const float b_inv = potential->scale_height_inv;
const float norm = potential->norm;
/* absolute value of height above disc */
const float dz = fabsf(g->x[2] - potential->z_disc);
/* vertical acceleration */
const float z_accel = 2.f * M_PI * phys_const->const_newton_G *
potential->surface_density *
tanhf(dz / potential->scale_height);
const float z_accel = norm * tanhf(dz * b_inv);
float dt = dt_dyn;
/* demand that dt * velocity < fraction of scale height of disc */
float dt1 = FLT_MAX;
if (g->v_full[2] != 0.f) {
dt1 = potential->scale_height / fabsf(g->v_full[2]);
if (dt1 < dt) dt = dt1;
const float dt1 = b / fabsf(g->v_full[2]);
dt = min(dt1, dt);
}
/* demand that dt^2 * acceleration < fraction of scale height of disc */
float dt2 = FLT_MAX;
if (z_accel != 0.f) {
dt2 = potential->scale_height / fabsf(z_accel);
const float dt2 = b / fabsf(z_accel);
if (dt2 < dt * dt) dt = sqrtf(dt2);
}
/* demand that dt^3 * jerk < fraction of scale height of disc */
float dt3 = FLT_MAX;
if (g->v_full[2] != 0.f) {
const float cosh_dz_inv = 1.f / coshf(dz * b_inv);
const float cosh_dz_inv2 = cosh_dz_inv * cosh_dz_inv;
const float dz_accel_over_dt =
2.f * M_PI * phys_const->const_newton_G * potential->surface_density /
potential->scale_height / coshf(dz / potential->scale_height) /
coshf(dz / potential->scale_height) * fabsf(g->v_full[2]);
norm * cosh_dz_inv2 * b_inv * fabsf(g->v_full[2]);
dt3 = potential->scale_height / fabsf(dz_accel_over_dt);
const float dt3 = b / fabsf(dz_accel_over_dt);
if (dt3 < dt * dt * dt) dt = cbrtf(dt3);
}
......@@ -126,6 +156,8 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep(
*
* See Creasey, Theuns & Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948,
* equation 17.
* We truncate the accelerations beyond z_trunc using a 1-cos(z) function
* that smoothly brings the accelerations to 0 at z_max.
*
* @param time The current time in internal units.
* @param potential The properties of the potential.
......@@ -137,19 +169,38 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
const struct phys_const* restrict phys_const, struct gpart* restrict g) {
const float dz = g->x[2] - potential->z_disc;
const float t_dyn = potential->dynamical_time;
float reduction_factor = 1.f;
if (time < potential->growth_time * t_dyn)
reduction_factor = time / (potential->growth_time * t_dyn);
/* Accelerations. Note that they are multiplied by G later on */
const float z_accel = reduction_factor * 2.f * M_PI *
potential->surface_density *
tanhf(fabsf(dz) / potential->scale_height);
const float abs_dz = fabsf(dz);
const float t_growth = potential->growth_time;
const float t_growth_inv = potential->growth_time_inv;
const float b_inv = potential->scale_height_inv;
const float z_trunc = potential->z_trunc;
const float z_max = potential->z_max;
const float z_trans_inv = potential->z_trans_inv;
const float norm_over_G = potential->norm_over_G;
/* Are we still growing the disc ? */
const float reduction_factor = time < t_growth ? time * t_growth_inv : 1.f;
/* Truncated or not ? */
float a_z;
if (abs_dz < z_trunc) {
/* Acc. 2 pi sigma tanh(z/b) */
a_z = reduction_factor * norm_over_G * tanhf(abs_dz * b_inv);
} else if (abs_dz < z_max) {
/* Acc. 2 pi sigma tanh(z/b) [1/2 + 1/2cos((z-zmax)/(pi z_trans))] */
a_z = reduction_factor * norm_over_G * tanhf(abs_dz * b_inv) *
(0.5f + 0.5f * cosf((abs_dz - z_trunc) * z_trans_inv));
} else {
/* Acc. 0 */
a_z = 0.f;
}
if (dz > 0) g->a_grav[2] -= z_accel;
if (dz < 0) g->a_grav[2] += z_accel;
/* Get the correct sign. Recall G is multipiled in later on */
if (dz > 0) g->a_grav[2] -= a_z;
if (dz < 0) g->a_grav[2] += a_z;
}
/**
......@@ -158,6 +209,8 @@ __attribute__((always_inline)) INLINE static void external_gravity_acceleration(
*
* See Creasey, Theuns & Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948,
* equation 22.
* We truncate the accelerations beyond z_trunc using a 1-cos(z) function
* that smoothly brings the accelerations to 0 at z_max.
*
* @param time The current time.
* @param potential The #external_potential used in the run.
......@@ -170,16 +223,35 @@ external_gravity_get_potential_energy(
const struct phys_const* const phys_const, const struct gpart* gp) {
const float dz = gp->x[2] - potential->z_disc;
const float t_dyn = potential->dynamical_time;
const float abs_dz = fabsf(dz);
const float t_growth = potential->growth_time;
const float t_growth_inv = potential->growth_time_inv;
const float b = potential->scale_height;
const float b_inv = potential->scale_height_inv;
const float norm = potential->norm;
const float z_trunc = potential->z_trunc;
const float z_max = potential->z_max;
/* Are we still growing the disc ? */
const float reduction_factor = time < t_growth ? time * t_growth_inv : 1.f;
/* Truncated or not ? */
float pot;
if (abs_dz < z_trunc) {
/* Potential (2 pi G sigma b ln(cosh(z/b)) */
pot = b * logf(coshf(dz * b_inv));
} else if (abs_dz < z_max) {
/* Potential. At z>>b, phi(z) = norm * z / b */
pot = 0.f;
float reduction_factor = 1.f;
if (time < potential->growth_time * t_dyn)
reduction_factor = time / (potential->growth_time * t_dyn);
} else {
/* Accelerations. Note that they are multiplied by G later on */
return reduction_factor * 2.f * M_PI * phys_const->const_newton_G *
potential->surface_density * potential->scale_height *
logf(coshf(dz / potential->scale_height));
pot = 0.f;
}
return pot * reduction_factor * norm;
}
/**
......@@ -205,13 +277,47 @@ static INLINE void potential_init_backend(
parameter_file, "DiscPatchPotential:scale_height");
potential->z_disc =
parser_get_param_double(parameter_file, "DiscPatchPotential:z_disc");
potential->z_trunc = parser_get_opt_param_double(
parameter_file, "DiscPatchPotential:z_trunc", FLT_MAX);
potential->z_max = parser_get_opt_param_double(
parameter_file, "DiscPatchPotential:z_max", FLT_MAX);
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.);
/* Compute the dynamical time */
potential->dynamical_time =
sqrt(potential->scale_height /
(phys_const->const_newton_G * potential->surface_density));
/* Convert the growth time multiplier to physical time */
potential->growth_time *= potential->dynamical_time;
/* Some cross-checks */
if (potential->z_trunc > potential->z_max)
error("Potential truncation z larger than maximal z");
if (potential->z_trunc < potential->scale_height)
error("Potential truncation z smaller than scale height");
/* Compute derived quantities */
potential->scale_height_inv = 1. / potential->scale_height;
potential->norm =
2. * M_PI * phys_const->const_newton_G * potential->surface_density;
potential->norm_over_G = 2 * M_PI * potential->surface_density;
potential->z_trans = potential->z_max - potential->z_trunc;
if (potential->z_trans != 0.f)
potential->z_trans_inv = 1. / (M_PI * potential->z_trans);
else
potential->z_trans_inv = FLT_MAX;
if (potential->growth_time != 0.)
potential->growth_time_inv = 1. / potential->growth_time;
else
potential->growth_time_inv = FLT_MAX;
}
/**
......@@ -223,13 +329,19 @@ static INLINE void potential_print_backend(
const struct external_potential* potential) {
message(
"External potential is 'Disk-patch' with properties surface_density = %e "
"disc height= %e scale height = %e timestep multiplier = %e.",
"External potential is 'Disk-patch' with Sigma=%f, z_disc=%f, b=%f and "
"dt_mult=%f.",
potential->surface_density, potential->z_disc, potential->scale_height,
potential->timestep_mult);
if (potential->z_max < FLT_MAX)
message("Potential will be truncated at z_trunc=%f and zeroed at z_max=%f",
potential->z_trunc, potential->z_max);
if (potential->growth_time > 0.)
message("Disc will grow for %f dynamical times.", potential->growth_time);
message("Disc will grow for %f [time_units]. (%f dynamical time)",
potential->growth_time,
potential->growth_time / potential->dynamical_time);
}
#endif /* SWIFT_DISC_PATCH_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