diff --git a/doc/RTD/source/ExternalPotentials/index.rst b/doc/RTD/source/ExternalPotentials/index.rst index ca33eb8189eea216863feb02579344aa22916696..329d33987209e986a5c3ca0675a285b73fc2e692 100644 --- a/doc/RTD/source/ExternalPotentials/index.rst +++ b/doc/RTD/source/ExternalPotentials/index.rst @@ -1,5 +1,6 @@ .. External potentials in SWIFT Folkert Nobels, 25th October 2018 + Alejandro Benitez-Llambay, October 2019 External Potentials =================== @@ -38,9 +39,17 @@ give a short overview of the potentials that are implemented in the code: This potential has as free parameters the concentration of the DM halo, the virial mass (:math:`M_{200}`) and the critical density. -7. Sine wave (sine-wave) -8. Point mass ring (point-mass-ring) -9. Disc Patch (disc-patch) +7. NFW poential + Miyamoto-Nagai potential (nfw_mn): This includes and NFW potential (identical to nfw) + plus an axisymmetric Miyamoto-Nagai potential. The Miyamoto-Nagai potential is given by: + + :math:`\Phi(R,z) = - \frac{G M_{d}}{\sqrt{R^2 + \left ( R_d + \sqrt{z^2 + Z_d^2} \right )^2}}`, + + where :math:`R^2 = x^2 + y^2` is the projected radius and :math:`M_d`, :math:`R_d`, :math:`Z_d` are the + mass, scalelength and scaleheight of the disk (in internal units), respectively. + +8. Sine wave (sine-wave) +9. Point mass ring (point-mass-ring) +10. Disc Patch (disc-patch) How to implement your own potential diff --git a/src/black_holes.h b/src/black_holes.h index d0c572061101eb8992b6420f4bc4b8acfb207dba..c45a8f112e8c449643f6502da259e232dda7ea6d 100644 --- a/src/black_holes.h +++ b/src/black_holes.h @@ -22,7 +22,7 @@ /* Config parameters. */ #include "../config.h" -/* Select the correct star model */ +/* Select the correct BH model */ #if defined(BLACK_HOLES_NONE) #include "./black_holes/Default/black_holes.h" #include "./black_holes/Default/black_holes_iact.h" diff --git a/src/black_holes_io.h b/src/black_holes_io.h index babe30656083d1e7e7b493850c421c562a2ce1ee..e82bcc813984a13b76f036edfea202dbaa5a07ce 100644 --- a/src/black_holes_io.h +++ b/src/black_holes_io.h @@ -24,13 +24,13 @@ /* Local includes */ #include "engine.h" -/* Load the correct star type */ +/* Load the correct BH model */ #if defined(BLACK_HOLES_NONE) #include "./black_holes/Default/black_holes_io.h" #elif defined(BLACK_HOLES_EAGLE) #include "./black_holes/EAGLE/black_holes_io.h" #else -#error "Invalid choice of star model" +#error "Invalid choice of BH model" #endif #endif /* SWIFT_BLACK_HOLES_IO_H */ diff --git a/src/black_holes_struct.h b/src/black_holes_struct.h index fdc70ebb7e73dab16bcf86b099a8ea3b61f51e6a..5e656d7876b168d45c4bd4bfe017d4b7f4371f4f 100644 --- a/src/black_holes_struct.h +++ b/src/black_holes_struct.h @@ -35,7 +35,7 @@ #elif defined(BLACK_HOLES_EAGLE) #include "./black_holes/EAGLE/black_holes_struct.h" #else -#error "Invalid choice of black holes function." +#error "Invalid choice of black hole model." #endif #endif /* SWIFT_BLACK_HOLES_STRUCT_H */ diff --git a/src/drift.h b/src/drift.h index b0b3e972995eefb57f3807034c26e2e1a6dbc200..74ac46346f8b2038f65793cb172f5138bbfa1545 100644 --- a/src/drift.h +++ b/src/drift.h @@ -173,7 +173,7 @@ __attribute__((always_inline)) INLINE static void drift_bpart( #ifdef SWIFT_DEBUG_CHECKS if (bp->ti_drift != ti_old) error( - "s-particle has not been drifted to the current time " + "b-particle has not been drifted to the current time " "bp->ti_drift=%lld, " "c->ti_old=%lld, ti_current=%lld", bp->ti_drift, ti_old, ti_current); diff --git a/src/feedback/EAGLE/feedback.c b/src/feedback/EAGLE/feedback.c index b25e71432c8ae9f825fbacf0b80131007fa35eb5..8f5b3e48652c9711fb4ec15ab13328c573bc9535 100644 --- a/src/feedback/EAGLE/feedback.c +++ b/src/feedback/EAGLE/feedback.c @@ -1171,7 +1171,7 @@ void feedback_props_init(struct feedback_props* fp, if (!fp->SNII_sampled_delay && fp->stellar_evolution_age_cut < fp->SNII_wind_delay) error( - "Time at which the enrichment downsampling stars is lower than the " + "Time at which the enrichment downsampling starts is lower than the " "SNII wind delay!"); /* Gather common conversion factors --------------------------------------- */ diff --git a/src/feedback/EAGLE/feedback.h b/src/feedback/EAGLE/feedback.h index 9a3b33aa5cfd102b3d429d9c928df70559edbb05..f306473efa3ecd46c799ef7007c9364a10601329 100644 --- a/src/feedback/EAGLE/feedback.h +++ b/src/feedback/EAGLE/feedback.h @@ -297,7 +297,7 @@ INLINE static void feedback_write_flavour(struct feedback_props* feedback, hid_t h_grp) { io_write_attribute_s(h_grp, "Feedback Model", "EAGLE"); -}; +} #endif // HAVE_HDF5 #endif /* SWIFT_FEEDBACK_EAGLE_H */ diff --git a/src/potential.h b/src/potential.h index 9d1418cadafcabe1d66658c43f2c1ef15778fc2b..f5ab1b9bc9c160f94ceff9138f5a8a8f46d2521f 100644 --- a/src/potential.h +++ b/src/potential.h @@ -38,6 +38,8 @@ #include "./potential/hernquist/potential.h" #elif defined(EXTERNAL_POTENTIAL_NFW) #include "./potential/nfw/potential.h" +#elif defined(EXTERNAL_POTENTIAL_NFW_MN) +#include "./potential/nfw_mn/potential.h" #elif defined(EXTERNAL_POTENTIAL_DISC_PATCH) #include "./potential/disc_patch/potential.h" #elif defined(EXTERNAL_POTENTIAL_SINE_WAVE) diff --git a/src/potential/nfw_mn/potential.h b/src/potential/nfw_mn/potential.h new file mode 100644 index 0000000000000000000000000000000000000000..355f1cbcd358d5fc7017886a433f547ab6399e85 --- /dev/null +++ b/src/potential/nfw_mn/potential.h @@ -0,0 +1,301 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2019 Alejandro Benitez-Llambay + * + * 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 <http://www.gnu.org/licenses/>. + * + ******************************************************************************/ + +#ifndef SWIFT_POTENTIAL_NFW_MN_H +#define SWIFT_POTENTIAL_NFW_MN_H + +/* Config parameters. */ +#include "../config.h" + +/* Some standard headers. */ +#include <float.h> +#include <math.h> + +/* Local includes. */ +#include "error.h" +#include "parser.h" +#include "part.h" +#include "physical_constants.h" +#include "space.h" +#include "units.h" + +/** + * @brief External Potential Properties - NFW Potential + Miyamoto-Nagai + * + * halo --> rho(r) = rho_0 / ( (r/R_s)*(1+r/R_s)^2 ) + * disk --> phi(R,z) = -G * Mdisk / (R^2 + (Rdisk + (z^2+Zdisk^2)^1/2)^2)^(1/2) + * + * We however parameterise this in terms of c and virial_mass, Mdisk, Rdisk + * and Zdisk + */ +struct external_potential { + + /*! Position of the centre of potential */ + double x[3]; + + /*! The scale radius of the NFW potential */ + double r_s; + + /*! The pre-factor \f$ 4 \pi G \rho_0 \r_s^3 \f$ */ + double pre_factor; + + /*! The critical density of the universe */ + double rho_c; + + /*! The concentration parameter */ + double c_200; + + /*! The virial mass */ + double M_200; + + /*! Disk Size */ + double Rdisk; + + /*! Disk height */ + double Zdisk; + + /*! Disk Mass */ + double Mdisk; + + /*! 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; + + /*! Softening length */ + double eps; +}; + +/** + * @brief Computes the time-step due to the acceleration from the NFW + MN + * 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 R2 = dx * dx + dy * dy; + const float r = sqrtf(R2 + 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 Vcirc_NFW = sqrtf((phys_const->const_newton_G * mr) / r); + const float f1 = + potential->Rdisk + sqrtf(potential->Zdisk * potential->Zdisk + dz * dz); + const float Vcirc_MN = sqrtf(phys_const->const_newton_G * potential->Mdisk * + R2 / pow(R2 + f1 * f1, 3.0 / 2.0)); + const float Vcirc = sqrtf(Vcirc_NFW * Vcirc_NFW + Vcirc_MN * Vcirc_MN); + + const float period = 2 * M_PI * r / Vcirc; + + /* 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 + + * MN disk. + * + * Note that the accelerations are multiplied by Newton's G constant + * later on. + * + * a_x = 4 pi \rho_0 r_s^3 ( 1/((r+rs)*r^2) - log(1+r/rs)/r^3) * x + * a_y = 4 pi \rho_0 r_s^3 ( 1/((r+rs)*r^2) - log(1+r/rs)/r^3) * y + * a_z = 4 pi \rho_0 r_s^3 ( 1/((r+rs)*r^2) - log(1+r/rs)/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) { + + 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]; + + /* First for the NFW part */ + const float R2 = dx * dx + dy * dy; + const float r = sqrtf(R2 + dz * dz + potential->eps * potential->eps); + const float term1 = potential->pre_factor; + const float term2 = (1.0f / ((r + potential->r_s) * r * r) - + logf(1.0f + r / potential->r_s) / (r * r * r)); + + g->a_grav[0] += term1 * term2 * dx; + g->a_grav[1] += term1 * term2 * dy; + g->a_grav[2] += term1 * term2 * dz; + + /* Now the the MN disk */ + const float f1 = sqrtf(potential->Zdisk * potential->Zdisk + dz * dz); + const float f2 = potential->Rdisk + f1; + const float f3 = pow(R2 + f2 * f2, -3.0 / 2.0); + + g->a_grav[0] -= potential->Mdisk * f3 * dx; + g->a_grav[1] -= potential->Mdisk * f3 * dy; + g->a_grav[2] -= potential->Mdisk * f3 * (f2 / f1) * dz; +} + +/** + * @brief Computes the gravitational potential energy of a particle in an + * NFW potential + MN potential. + * + * phi = -4 * pi * G * rho_0 * r_s^3 * ln(1+r/r_s) - G * Mdisk / sqrt(R^2 + + * (Rdisk + sqrt(z^2 + Zdisk^2))^2) + * + * @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]; + + /* First for the NFW profile */ + const float R2 = dx * dx + dy * dy; + const float r = sqrtf(R2 + dz * dz + potential->eps * potential->eps); + const float term1 = -potential->pre_factor / r; + const float term2 = logf(1.0f + r / potential->r_s); + + /* Now for the MN disk */ + const float MN_term = potential->Rdisk + sqrtf(potential->Zdisk + dz * dz); + const float MN_pot = -potential->Mdisk / sqrtf(R2 + MN_term * MN_term); + + return term1 * term2 + MN_pot; +} + +/** + * @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, "NFW_MNPotential: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, "NFW_MNPotential: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, "NFW_MNPotential:timestep_mult"); + potential->c_200 = + parser_get_param_double(parameter_file, "NFW_MNPotential:concentration"); + potential->M_200 = + parser_get_param_double(parameter_file, "NFW_MNPotential:M_200"); + potential->rho_c = parser_get_param_double( + parameter_file, "NFW_MNPotential:critical_density"); + potential->Mdisk = + parser_get_param_double(parameter_file, "NFW_MNPotential:Mdisk"); + potential->Rdisk = + parser_get_param_double(parameter_file, "NFW_MNPotential:Rdisk"); + potential->Zdisk = + parser_get_param_double(parameter_file, "NFW_MNPotential:Zdisk"); + + potential->eps = 0.05; + + /* Compute R_200 */ + const double R_200 = + cbrtf(3.0 * potential->M_200 / (4. * M_PI * 200.0 * potential->rho_c)); + + /* NFW scale-radius */ + potential->r_s = R_200 / potential->c_200; + const double r_s3 = potential->r_s * potential->r_s * potential->r_s; + + /* Log(c_200) term appearing in many expressions */ + potential->log_c200_term = + log(1. + potential->c_200) - potential->c_200 / (1. + potential->c_200); + + const double rho_0 = + potential->M_200 / (4.f * M_PI * r_s3 * potential->log_c200_term); + + /* Pre-factor for the accelerations (note G is multiplied in later on) */ + potential->pre_factor = 4.0f * M_PI * rho_0 * r_s3; + + /* 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 + MN disk' 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); +} + +#endif /* SWIFT_POTENTIAL_NFW_MN_H */ diff --git a/src/units.c b/src/units.c index e6bba3ac6651d9dc6f8ef0eda5766a90bd3710e2..75b26e5561d458515b4f422949e801e70220df7f 100644 --- a/src/units.c +++ b/src/units.c @@ -308,6 +308,7 @@ void units_get_base_unit_exponents_array(float baseUnitsExp[5], break; case UNIT_CONV_ENERGY_PER_UNIT_MASS: + case UNIT_CONV_VELOCITY_SQUARED: baseUnitsExp[UNIT_LENGTH] = 2.f; baseUnitsExp[UNIT_TIME] = -2.f; break; @@ -408,6 +409,16 @@ void units_get_base_unit_exponents_array(float baseUnitsExp[5], baseUnitsExp[UNIT_TIME] = -1.f; break; + case UNIT_CONV_DIFF_RATE: + baseUnitsExp[UNIT_TIME] = -1.f; + break; + + case UNIT_CONV_DIFF_COEFF: + baseUnitsExp[UNIT_MASS] = 1.f; + baseUnitsExp[UNIT_LENGTH] = -1.f; + baseUnitsExp[UNIT_TIME] = -1.f; + break; + default: error("Invalid choice of pre-defined units"); break; diff --git a/src/units.h b/src/units.h index 898221d95b1a74bca1318ec6a2dbf71f8f3a4760..59ff31dff80553d47da8a49ec24b2792f8f9b4d1 100644 --- a/src/units.h +++ b/src/units.h @@ -100,7 +100,10 @@ enum unit_conversion_factor { UNIT_CONV_INV_VOLUME, UNIT_CONV_SFR, UNIT_CONV_SSFR, - UNIT_CONV_MASS_PER_UNIT_TIME + UNIT_CONV_DIFF_RATE, + UNIT_CONV_DIFF_COEFF, + UNIT_CONV_MASS_PER_UNIT_TIME, + UNIT_CONV_VELOCITY_SQUARED }; void units_init_cgs(struct unit_system*);