Skip to content
Snippets Groups Projects
Commit cbcd454b authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'hernquist-sdmh05-fix' into 'master' 

Fix factor of 2 in equation for the scale radius and tiny change in plotting script.

See merge request !1646
parents 66109803 af35dbd9
No related branches found
No related tags found
5 merge requests!1715Update planetary strength after planetary plus's master rebase,!1693More threapool plotting tweaks,!1668before Mag.Egy in all the flavors,!1662Initial sync from previous months,!1646Fix factor of 2 in equation for the scale radius and tiny change in plotting script.
Hide whitespace changes
Inline Side-by-side
doc/RTD/source/ExternalPotentials/index.rst
+ 63
− 29
View file @ cbcd454b
...@@ -150,22 +150,22 @@ The parameters of the model are: ...@@ -150,22 +150,22 @@ The parameters of the model are:
6. Hernquist potential (``hernquist``) 6. Hernquist potential (``hernquist``)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
A potential that is given by the Hernquist (1990) We can set up a potential as given by Hernquist (1990):
potential:
:math:`\Phi(r) = - \frac{G_{\rm N}M}{r+a}.` :math:`\Phi(r) = - \frac{G_{\rm N}M}{r+a},`
The free parameters of the Hernquist potential are mass, scale length, where :math:`M` is the total Hernquist mass and :math: `a` is the Hernquist-
and softening. The potential can be set at any position in the box. equivalent scale radius of the potential. The potential can be set at any
The potential add an additional time step constrain that limits the position in the box. It adds an additional time step constraint that limits
time step to a maximum of a specified fraction of the circular orbital the time step to a maximum of a specified fraction of the circular orbital
time at the current radius of the particle. The other criteria time at the current radius of the particle. The other criteria
(CFL, self-gravity, ...) are applied on top of this criterion. (CFL, self-gravity, ...) are applied on top of this criterion. For example, a
fraction of 0.01 means that an orbital period would be resolved by 100 time steps.
The Hernquist potential can be run in the most basic version, then only the In the most basic version, the Hernquist potential can be run by providing
Hernquist mass, scale length, softening length and fraction of the only the Hernquist mass, scale length, softening length and fraction of the
orbital time for the time step limit are used, the parameters of the orbital time for the time stepping. In this case the model parameters may
model in this case are: look something like:
.. code:: YAML .. code:: YAML
...@@ -173,20 +173,23 @@ model in this case are: ...@@ -173,20 +173,23 @@ model in this case are:
useabspos: 0 # 0 -> positions based on centre, 1 -> absolute positions useabspos: 0 # 0 -> positions based on centre, 1 -> absolute positions
position: [0.,0.,0.] # Location of centre of isothermal potential with respect to centre of the box (if 0) otherwise absolute (if 1) (internal units) position: [0.,0.,0.] # Location of centre of isothermal potential with respect to centre of the box (if 0) otherwise absolute (if 1) (internal units)
mass: 1e12 # Mass of the Hernquist potential mass: 1e12 # Mass of the Hernquist potential
scalelength: 2.0 # scale length a scalelength: 2.0 # scale length a (internal units)
timestep_mult: 0.01 # Dimensionless pre-factor for the time-step condition, basically determines the fraction of the orbital time we use to do the time integration timestep_mult: 0.01 # Dimensionless pre-factor for the time-step condition, determines the fraction of the orbital time we use to do the time integration; fraction of 0.01 means we resolve an orbit with 100 timesteps
epsilon: 0.2 # Softening size (internal units) epsilon: 0.2 # Softening size (internal units)
Besides the basic version, it is also possible to run the Hernquist Besides the basic version, it is also possible to run the Hernquist
potential for idealised disk galaxies that follow the approach of potential for idealised disk galaxies that follow the approach of
Springel+ 2005. The default Hernquist potential uses a corrected value `Springel, Di Matteo & Hernquist (2005)
for the formulation that improves the match with the NFW (below) with <https://ui.adsabs.harvard.edu/abs/2005MNRAS.361..776S/abstract>`_. This
the same M200 (Nobels+ in prep). Contrary to above, the idealised disk potential, however, uses a corrected value of the formulation that improves
setup runs with a specified M200, concentration and reduced Hubble the match with the NFW profile (below) with the same M200 (Nobels+ in prep).
constant that set both the mass and scale length parameter. The reduced Contrary to the above (idealizeddisk: 0 setup), the idealised disk setup runs
Hubble constant is only used to determine R200. by specifying one out of :math:`M_{200}`, :math:`V_{200}`, or :math:`R_{200}`,
plus concentration and reduced Hubble constant.
The parameters of the model are:
In this case, we don't provide the 'mass' and 'scalelength' parameters, but
'M200' (or 'V200', or 'R200') and 'concentration' :math:`c`, as well as reduced Hubble
constant :math:`h` to define the potential. The parameters of the model may look something like:
.. code:: YAML .. code:: YAML
...@@ -199,21 +202,52 @@ The parameters of the model are: ...@@ -199,21 +202,52 @@ The parameters of the model are:
concentration: 9.0 # concentration of the Halo concentration: 9.0 # concentration of the Halo
diskfraction: 0.040 # Disk mass fraction diskfraction: 0.040 # Disk mass fraction
bulgefraction: 0.0 # Bulge mass fraction bulgefraction: 0.0 # Bulge mass fraction
timestep_mult: 0.01 # Dimensionless pre-factor for the time-step condition, basically determines the fraction of the orbital time we use to do the time integration timestep_mult: 0.01 # Dimensionless pre-factor for the time-step condition, determines the fraction of the orbital time we use to do the time integration; fraction of 0.01 means we resolve an orbit with 100 timesteps
epsilon: 0.2 # Softening size (internal units) epsilon: 0.2 # Softening size (internal units)
The user should specify one out of 'M200', 'V200', or 'R200' to define
the potential. The reduced Hubble constant is then used to determine the
other two. Then, the scale radius is calculated as :math:`R_s = R_{200}/c`,
where :math:`c` is the concentration, and the Hernquist-equivalent scale-length
is calculated as:
:math:`a = \frac{b+\sqrt{b}}{1-b} R_{200}`,
where :math:`b = \frac{2}{c^2}(\ln(1+c) - \frac{c}{1+c})`.
Two examples using the Hernquist potential can be found in ``swiftsim/examples/GravityTests/``.
The ``Hernquist_radialinfall`` example puts 5 particles with zero velocity in a Hernquist
potential, resulting in radial orbits. The ``Hernquist_circularorbit``example puts three
particles on a circular orbit in a Hernquist potential, one at the inner region, one at the
maximal circular velocity, and one in the outer region. To run these examples, SWIFT must
be configured with the flag ``--with-ext-potential=hernquist``, or ``hernquist-sdmh05``
(see below).
7. Hernquist potential (``hernquist-sdmh05``) 7. Hernquist SDMH05 potential (``hernquist-sdmh05``)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
This is the same potential as Hernquist with the difference being the This is the same potential as Hernquist with the difference being the
way that the idealised disk part is calculated. This potential uses way that the idealised disk part is calculated. This potential uses
exactly the same approach as `Springel, Di Matteo & Hernquist (2005) exactly the same approach as `Springel, Di Matteo & Hernquist (2005)
<https://ui.adsabs.harvard.edu/abs/2005MNRAS.361..776S/abstract>`_, <https://ui.adsabs.harvard.edu/abs/2005MNRAS.361..776S/abstract>`_,
this means that ICs generated with the original `MakeNewDisk` code can which means that ICs generated with the original `MakeNewDisk` code can
be used when using this potential. Contrary to the updated Hernquist be used with this potential. Contrary to the updated Hernquist
potential (above) it is not possible to have an identically matched potential (above), it is not possible to have an identically matched
NFW potential. NFW potential in this case.
The parameters needed for this potential are the same set of variables as
above, i.e. 'mass' and 'scalelength' when we don't use the idealised
disk, and 'concentration' plus one out of 'M200', 'V200', or 'R200' if
we do. As one of the three is provided, the reduced Hubble constant is
used to calculate the other two. Then, the scale radius is calculated
using the NFW definition, :math:`R_s = R_{200}/c`, and the Hernquist-
equivalent scale length is given by
:math:`a = R_s \sqrt{2(\ln(1+c) - \frac{c}{1+c})}.`
Runs that provide e.g. M200 and c (using idealised disk) are thus equivalent
to providing mass and scale length if calculated by the above equation (without
idealized disk).
8. Navarro-Frenk-White potential (``nfw``): 8. Navarro-Frenk-White potential (``nfw``):
......
examples/GravityTests/Hernquist_circularorbit/plotprog.py
+ 1
− 1
View file @ cbcd454b
...@@ -22,7 +22,7 @@ import h5py ...@@ -22,7 +22,7 @@ import h5py
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from scipy.integrate import odeint from scipy.integrate import odeint
t = np.linspace(0, 40, 1e5) t = np.linspace(0, 40, 100000)
y0 = [0, 10] y0 = [0, 10]
a = 30.0 a = 30.0
G = 4.300927e-06 G = 4.300927e-06
......
src/potential/hernquist/potential.h
+ 19
− 17
View file @ cbcd454b
...@@ -73,6 +73,9 @@ struct external_potential { ...@@ -73,6 +73,9 @@ struct external_potential {
* time to get the time steps */ * time to get the time steps */
double timestep_mult; double timestep_mult;
/*! Inverse of the sqrt of G*M, a common factor */
double sqrtgm_inv;
/*! Mode to use 0 for simplest form of potential purely 1 for idealized /*! Mode to use 0 for simplest form of potential purely 1 for idealized
* galaxies */ * galaxies */
int usedisk; int usedisk;
...@@ -92,21 +95,18 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep( ...@@ -92,21 +95,18 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep(
const struct phys_const* restrict phys_const, const struct phys_const* restrict phys_const,
const struct gpart* restrict g) { const struct gpart* restrict g) {
const float G_newton = phys_const->const_newton_G;
/* Calculate the relative potential with respect to the centre of the /* Calculate the relative potential with respect to the centre of the
* potential */ * potential */
const float dx = g->x[0] - potential->x[0]; const float dx = g->x[0] - potential->x[0];
const float dy = g->x[1] - potential->x[1]; const float dy = g->x[1] - potential->x[1];
const float dz = g->x[2] - potential->x[2]; const float dz = g->x[2] - potential->x[2];
/* calculate the radius */ /* Calculate the radius */
const float r = sqrtf(dx * dx + dy * dy + dz * dz + potential->epsilon2); const float r = sqrtf(dx * dx + dy * dy + dz * dz + potential->epsilon2);
const float sqrtgm_inv = 1.f / sqrtf(G_newton * potential->mass);
/* Calculate the circular orbital period */ /* Calculate the circular orbital period */
const float period = 2.f * M_PI * sqrtf(r) * potential->al * const float period =
(1 + r / potential->al) * sqrtgm_inv; 2.f * M_PI * sqrtf(r) * (potential->al + r) * potential->sqrtgm_inv;
/* Time-step as a fraction of the cirecular orbital time */ /* Time-step as a fraction of the cirecular orbital time */
const float time_step = potential->timestep_mult * period; const float time_step = potential->timestep_mult * period;
...@@ -211,7 +211,7 @@ static INLINE void potential_init_backend( ...@@ -211,7 +211,7 @@ static INLINE void potential_init_backend(
potential->x[2] += s->dim[2] / 2.; potential->x[2] += s->dim[2] / 2.;
} }
/* check whether we use the more advanced idealized disk setting */ /* Check whether we use the more advanced idealized disk setting */
potential->usedisk = parser_get_opt_param_int( potential->usedisk = parser_get_opt_param_int(
parameter_file, "HernquistPotential:idealizeddisk", parameter_file, "HernquistPotential:idealizeddisk",
idealized_disk_default); idealized_disk_default);
...@@ -269,7 +269,7 @@ static INLINE void potential_init_backend( ...@@ -269,7 +269,7 @@ static INLINE void potential_init_backend(
potential->M200 = M200; potential->M200 = M200;
potential->R200 = R200; potential->R200 = R200;
/* get the concentration from the parameter file */ /* Get the concentration from the parameter file */
potential->c = parser_get_param_double(parameter_file, potential->c = parser_get_param_double(parameter_file,
"HernquistPotential:concentration"); "HernquistPotential:concentration");
...@@ -286,7 +286,7 @@ static INLINE void potential_init_backend( ...@@ -286,7 +286,7 @@ static INLINE void potential_init_backend(
const double b = 2. * cc_inv * cc_inv * (log(1. + cc) - cc / (1. + cc)); const double b = 2. * cc_inv * cc_inv * (log(1. + cc) - cc / (1. + cc));
/* Calculate the Hernquist equivalent scale length */ /* Calculate the Hernquist equivalent scale length */
potential->al = (b + sqrt(b)) / (1 - b) * R200; potential->al = R200 * (b + sqrt(b)) / (1 - b);
/* Define R200 inv*/ /* Define R200 inv*/
const double R200_inv = 1. / R200; const double R200_inv = 1. / R200;
...@@ -296,9 +296,9 @@ static INLINE void potential_init_backend( ...@@ -296,9 +296,9 @@ static INLINE void potential_init_backend(
(potential->R200 + potential->al) * R200_inv * (potential->R200 + potential->al) * R200_inv *
R200_inv * potential->M200; R200_inv * potential->M200;
/* Depending on the disk mass and and the bulge mass the halo /* Depending on the disk mass and and the bulge mass, the halo
* gets a different mass, because of this we read the fractions * gets a different mass. Because of this, we read the fractions
* from the parameter file and calculate the absolute mass*/ * from the parameter file and calculate the absolute mass */
const double diskfraction = parser_get_param_double( const double diskfraction = parser_get_param_double(
parameter_file, "HernquistPotential:diskfraction"); parameter_file, "HernquistPotential:diskfraction");
const double bulgefraction = parser_get_param_double( const double bulgefraction = parser_get_param_double(
...@@ -318,12 +318,14 @@ static INLINE void potential_init_backend( ...@@ -318,12 +318,14 @@ static INLINE void potential_init_backend(
parser_get_param_double(parameter_file, "HernquistPotential:epsilon"); parser_get_param_double(parameter_file, "HernquistPotential:epsilon");
potential->epsilon2 = epsilon * epsilon; potential->epsilon2 = epsilon * epsilon;
/* Compute the minimal time-step. */ /* Calculate a common factor in the calculation, i.e. 1/sqrt(GM)*/
/* This is the circular orbital time at the softened radius */
const float sqrtgm = sqrtf(phys_const->const_newton_G * potential->mass); const float sqrtgm = sqrtf(phys_const->const_newton_G * potential->mass);
potential->mintime = 2.f * sqrtf(epsilon) * potential->al * M_PI * potential->sqrtgm_inv = 1. / sqrtgm;
(1. + epsilon / potential->al) / sqrtgm *
potential->timestep_mult; /* Compute the minimal time-step. */
/* This is a fraction of the circular orbital time at the softened radius */
potential->mintime = potential->timestep_mult * 2.f * sqrtf(epsilon) * M_PI *
(potential->al + epsilon) * potential->sqrtgm_inv;
} }
/** /**
......
src/potential/hernquist_sdmh05/potential.h
+ 18
− 16
View file @ cbcd454b
...@@ -63,6 +63,9 @@ struct external_potential { ...@@ -63,6 +63,9 @@ struct external_potential {
/*! Time-step condition pre-factor, is multiplied times the circular orbital /*! Time-step condition pre-factor, is multiplied times the circular orbital
* time to get the time steps */ * time to get the time steps */
double timestep_mult; double timestep_mult;
/* Additional common parameter inverse of sqrt(GM)*/
double sqrtgm_inv;
}; };
/** /**
...@@ -79,8 +82,6 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep( ...@@ -79,8 +82,6 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep(
const struct phys_const* restrict phys_const, const struct phys_const* restrict phys_const,
const struct gpart* restrict g) { const struct gpart* restrict g) {
const float G_newton = phys_const->const_newton_G;
/* Calculate the relative potential with respect to the centre of the /* Calculate the relative potential with respect to the centre of the
* potential */ * potential */
const float dx = g->x[0] - potential->x[0]; const float dx = g->x[0] - potential->x[0];
...@@ -89,11 +90,10 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep( ...@@ -89,11 +90,10 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep(
/* calculate the radius */ /* calculate the radius */
const float r = sqrtf(dx * dx + dy * dy + dz * dz + potential->epsilon2); const float r = sqrtf(dx * dx + dy * dy + dz * dz + potential->epsilon2);
const float sqrtgm_inv = 1.f / sqrtf(G_newton * potential->mass);
/* Calculate the circular orbital period */ /* Calculate the circular orbital period */
const float period = 2.f * M_PI * sqrtf(r) * potential->al * const float period =
(1 + r / potential->al) * sqrtgm_inv; 2.f * M_PI * sqrtf(r) * (potential->al + r) * potential->sqrtgm_inv;
/* Time-step as a fraction of the cirecular orbital time */ /* Time-step as a fraction of the cirecular orbital time */
const float time_step = potential->timestep_mult * period; const float time_step = potential->timestep_mult * period;
...@@ -254,7 +254,7 @@ static INLINE void potential_init_backend( ...@@ -254,7 +254,7 @@ static INLINE void potential_init_backend(
/* message("M200 = %g, R200 = %g, V200 = %g", M200, R200, V200); */ /* message("M200 = %g, R200 = %g, V200 = %g", M200, R200, V200); */
/* message("H0 = %g", H0); */ /* message("H0 = %g", H0); */
/* get the concentration from the parameter file */ /* Get the concentration from the parameter file */
const double concentration = parser_get_param_double( const double concentration = parser_get_param_double(
parameter_file, "HernquistPotential:concentration"); parameter_file, "HernquistPotential:concentration");
...@@ -262,12 +262,12 @@ static INLINE void potential_init_backend( ...@@ -262,12 +262,12 @@ static INLINE void potential_init_backend(
const double RS = R200 / concentration; const double RS = R200 / concentration;
/* Calculate the Hernquist equivalent scale length */ /* Calculate the Hernquist equivalent scale length */
potential->al = RS * sqrt(1. * (log(1. + concentration) - potential->al = RS * sqrt(2. * (log(1. + concentration) -
concentration / (1. + concentration))); concentration / (1. + concentration)));
/* Depending on the disk mass and and the bulge mass the halo /* Depending on the disk mass and the bulge mass, the halo
* gets a different mass, because of this we read the fractions * gets a different mass. Because of this, we read the fractions
* from the parameter file and calculate the absolute mass*/ * from the parameter file and calculate the absolute mass */
const double diskfraction = parser_get_param_double( const double diskfraction = parser_get_param_double(
parameter_file, "HernquistPotential:diskfraction"); parameter_file, "HernquistPotential:diskfraction");
const double bulgefraction = parser_get_param_double( const double bulgefraction = parser_get_param_double(
...@@ -287,12 +287,14 @@ static INLINE void potential_init_backend( ...@@ -287,12 +287,14 @@ static INLINE void potential_init_backend(
parser_get_param_double(parameter_file, "HernquistPotential:epsilon"); parser_get_param_double(parameter_file, "HernquistPotential:epsilon");
potential->epsilon2 = epsilon * epsilon; potential->epsilon2 = epsilon * epsilon;
/* Compute the minimal time-step. */ /* Calculate a common factor in the calculation, i.e. 1/sqrt(GM)*/
/* This is the circular orbital time at the softened radius */
const float sqrtgm = sqrtf(phys_const->const_newton_G * potential->mass); const float sqrtgm = sqrtf(phys_const->const_newton_G * potential->mass);
potential->mintime = 2.f * sqrtf(epsilon) * potential->al * M_PI * potential->sqrtgm_inv = 1. / sqrtgm;
(1. + epsilon / potential->al) / sqrtgm *
potential->timestep_mult; /* Compute the minimal time-step. */
/* This is a fraction of the circular orbital time at the softened radius */
potential->mintime = potential->timestep_mult * 2.f * sqrtf(epsilon) * M_PI *
(potential->al + epsilon) * potential->sqrtgm_inv;
} }
/** /**
...@@ -304,7 +306,7 @@ static inline void potential_print_backend( ...@@ -304,7 +306,7 @@ static inline void potential_print_backend(
const struct external_potential* potential) { const struct external_potential* potential) {
message( message(
"external potential is 'hernquist Springel, Di Matteo & Herquist 2005' " "external potential is 'hernquist Springel, Di Matteo & Hernquist 2005' "
"with properties are (x,y,z) = (%e, %e, %e), mass = %e scale length = %e " "with properties are (x,y,z) = (%e, %e, %e), mass = %e scale length = %e "
", minimum time = %e timestep multiplier = %e", ", minimum time = %e timestep multiplier = %e",
potential->x[0], potential->x[1], potential->x[2], potential->mass, potential->x[0], potential->x[1], potential->x[2], potential->mass,
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment