/******************************************************************************* * This file is part of SWIFT. * Copyright (c) 2018 Ashley Kelly () * Folkert Nobels (nobels@strw.leidenuniv.nl) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published * by the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU Lesser General Public License * along with this program. If not, see . * ******************************************************************************/ #ifndef SWIFT_POTENTIAL_NFW_H #define SWIFT_POTENTIAL_NFW_H /* Config parameters. */ #include /* Some standard headers. */ #include #include /* Local includes. */ #include "error.h" #include "gravity.h" #include "parser.h" #include "part.h" #include "physical_constants.h" #include "space.h" #include "units.h" /** * @brief External Potential Properties - NFW Potential rho(r) = rho_0 / ( (r/R_s)*(1+r/R_s)^2 ) We however parameterise this in terms of c and virial_mass */ struct external_potential { /*! Position of the centre of potential */ double x[3]; /*! The scale radius of the NFW potential */ double r_s; /*! The critical density of the universe */ double rho_c; /*! The concentration parameter */ double c_200; /*! The mass at R200 */ double M_200; /*! R200 */ double R_200; /*! Time-step condition pre_factor, this factor is used to multiply times the * orbital time, so in the case of 0.01 we take 1% of the orbital time as * the time integration steps */ double timestep_mult; /*! Minimum time step based on the orbital time at the softening times * the timestep_mult */ double mintime; /*! Common log term \f$ \ln(1+c_{200}) - \frac{c_{200}}{1 + c_{200}} \f$ */ double log_c200_term; /*! M_200 times inverse of common log term \f$ \ln(1+c_{200}) - * \frac{c_{200}}{1 + c_{200}} \f$ */ double M_200_times_log_c200_term_inv; /*! Softening length */ double eps; /*! Bulge fraction */ double bulgefraction; /*! disk fraction */ double diskfraction; }; /** * @brief Computes the enclosed mass due to the NFW potential * * @param potential The #external_potential used in the run. * @param radius The radius of the particle */ __attribute__((always_inline)) INLINE static float enclosed_mass_NFW( const struct external_potential* restrict potential, const double r) { const double r_over_Rs = r / potential->r_s; return potential->M_200_times_log_c200_term_inv * (log(1 + r_over_Rs) - r / (r + potential->r_s)); } /** * @brief Computes the time-step due to the acceleration from the NFW potential * as a fraction (timestep_mult) of the circular orbital time of that * particle. * * @param time The current time. * @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 float external_gravity_timestep( double time, const struct external_potential* restrict potential, const struct phys_const* restrict phys_const, const struct gpart* restrict g) { const float dx = g->x[0] - potential->x[0]; const float dy = g->x[1] - potential->x[1]; const float dz = g->x[2] - potential->x[2]; const float r = sqrtf(dx * dx + dy * dy + dz * dz + potential->eps * potential->eps); const float mr = potential->M_200 * (logf(1.f + r / potential->r_s) - r / (r + potential->r_s)) / potential->log_c200_term; const float period = 2 * M_PI * r * sqrtf(r / (phys_const->const_newton_G * mr)); /* Time-step as a fraction of the circular period */ const float time_step = potential->timestep_mult * period; return max(time_step, potential->mintime); } /** * @brief Computes the gravitational acceleration from an NFW Halo potential. * * Note that the accelerations are multiplied by Newton's G constant * later on. * * a_x = M_encl(r) /r^3 * x * a_y = M_encl(r) /r^3 * y * a_z = M_encl(r) /r^3 * z * * @param time The current time. * @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_acceleration( double time, const struct external_potential* restrict potential, const struct phys_const* restrict phys_const, struct gpart* restrict g) { /* Determine the position relative to the centre of the potential */ const float dx = g->x[0] - potential->x[0]; const float dy = g->x[1] - potential->x[1]; const float dz = g->x[2] - potential->x[2]; /* Calculate the acceleration */ const float r2 = dx * dx + dy * dy + dz * dz + potential->eps * potential->eps; const float r = sqrtf(r2); const float r_inv = 1.f / r; const float M_encl = enclosed_mass_NFW(potential, r); const float acc = -M_encl * r_inv * r_inv * r_inv; const float pot = -potential->M_200_times_log_c200_term_inv * r_inv * logf(1.f + r / potential->r_s); g->a_grav[0] += acc * dx; g->a_grav[1] += acc * dy; g->a_grav[2] += acc * dz; gravity_add_comoving_potential(g, pot); } /** * @brief Computes the gravitational potential energy of a particle in an * NFW potential. * * phi = -4 * pi * G * rho_0 * r_s^3 * ln(1+r/r_s) * * @param time The current time (unused here). * @param potential The #external_potential used in the run. * @param phys_const Physical constants in internal units. * @param g Pointer to the particle data. */ __attribute__((always_inline)) INLINE static float external_gravity_get_potential_energy( double time, const struct external_potential* potential, const struct phys_const* const phys_const, const struct gpart* g) { const float dx = g->x[0] - potential->x[0]; const float dy = g->x[1] - potential->x[1]; const float dz = g->x[2] - potential->x[2]; const float r = sqrtf(dx * dx + dy * dy + dz * dz + potential->eps * potential->eps); const float term1 = -potential->M_200_times_log_c200_term_inv / r; const float term2 = logf(1.0f + r / potential->r_s); return phys_const->const_newton_G * term1 * term2; } /** * @brief Initialises the external potential properties in the internal system * of units. * * @param parameter_file The parsed parameter file * @param phys_const Physical constants in internal units * @param us The current internal system of units * @param potential The external potential properties to initialize */ static INLINE void potential_init_backend( struct swift_params* parameter_file, const struct phys_const* phys_const, const struct unit_system* us, const struct space* s, struct external_potential* potential) { /* Read in the position of the centre of potential */ parser_get_param_double_array(parameter_file, "NFWPotential:position", 3, potential->x); /* Is the position absolute or relative to the centre of the box? */ const int useabspos = parser_get_param_int(parameter_file, "NFWPotential:useabspos"); if (!useabspos) { potential->x[0] += s->dim[0] / 2.; potential->x[1] += s->dim[1] / 2.; potential->x[2] += s->dim[2] / 2.; } /* Read the other parameters of the model */ potential->timestep_mult = parser_get_param_double(parameter_file, "NFWPotential:timestep_mult"); potential->c_200 = parser_get_param_double(parameter_file, "NFWPotential:concentration"); potential->M_200 = parser_get_param_double(parameter_file, "NFWPotential:M_200"); potential->eps = parser_get_param_double(parameter_file, "NFWPotential:epsilon"); const double h = parser_get_param_double(parameter_file, "NFWPotential:h"); potential->bulgefraction = parser_get_opt_param_double( parameter_file, "NFWPotential:bulgefraction", 0.0); potential->diskfraction = parser_get_opt_param_double( parameter_file, "NFWPotential:diskfraction", 0.0); /* Some constants we need to calculate the critical density */ const double G_newton = phys_const->const_newton_G; const double kmoversoverMpc = phys_const->const_reduced_hubble; /* Hubble constant assumed for halo masses conversion */ const double H0 = h * kmoversoverMpc; /* Compute R_200 for this use the parameter critical density*/ potential->R_200 = cbrt(10. * potential->M_200 * G_newton * H0) / (10 * H0); /* NFW scale-radius */ potential->r_s = potential->R_200 / potential->c_200; /* Log(c_200) term appearing in many expressions */ potential->log_c200_term = log(1. + potential->c_200) - potential->c_200 / (1. + potential->c_200); potential->M_200_times_log_c200_term_inv = potential->M_200 * (1 - potential->bulgefraction - potential->diskfraction) / potential->log_c200_term; /* Compute the orbital time at the softening radius */ const double sqrtgm = sqrt(phys_const->const_newton_G * potential->M_200); const double epslnthing = log(1.f + potential->eps / potential->r_s) - potential->eps / (potential->eps + potential->r_s); potential->mintime = 2. * M_PI * potential->eps * sqrtf(potential->eps) * sqrtf(potential->log_c200_term / epslnthing) / sqrtgm * potential->timestep_mult; } /** * @brief Prints the properties of the external potential to stdout. * * @param potential The external potential properties. */ static INLINE void potential_print_backend( const struct external_potential* potential) { message( "External potential is 'NFW' with properties are (x,y,z) = (%e, " "%e, %e), scale radius = %e " "timestep multiplier = %e, mintime = %e", potential->x[0], potential->x[1], potential->x[2], potential->r_s, potential->timestep_mult, potential->mintime); message("Properties of the halo M200 = %e, R200 = %e, c = %e", potential->M_200, potential->R_200, potential->c_200); if ((potential->bulgefraction > 0.) || (potential->diskfraction > 0.)) { message("bulge fraction = %e, disk fraction = %e", potential->bulgefraction, potential->diskfraction); } } #endif /* SWIFT_POTENTIAL_NFW_H */