Skip to content
Snippets Groups Projects
Commit af35dbd9 authored by Zorry Belcheva's avatar Zorry Belcheva Committed by Matthieu Schaller
Browse files

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

parent 66109803
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.
......@@ -150,22 +150,22 @@ The parameters of the model are:
6. Hernquist potential (``hernquist``)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
A potential that is given by the Hernquist (1990)
potential:
We can set up a potential as given by Hernquist (1990):
: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,
and softening. The potential can be set at any position in the box.
The potential add an additional time step constrain that limits the
time step to a maximum of a specified fraction of the circular orbital
where :math:`M` is the total Hernquist mass and :math: `a` is the Hernquist-
equivalent scale radius of the potential. The potential can be set at any
position in the box. It adds an additional time step constraint that limits
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
(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
Hernquist mass, scale length, softening length and fraction of the
orbital time for the time step limit are used, the parameters of the
model in this case are:
In the most basic version, the Hernquist potential can be run by providing
only the Hernquist mass, scale length, softening length and fraction of the
orbital time for the time stepping. In this case the model parameters may
look something like:
.. code:: YAML
......@@ -173,20 +173,23 @@ model in this case are:
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)
mass: 1e12 # Mass of the Hernquist potential
scalelength: 2.0 # scale length a
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
scalelength: 2.0 # scale length a (internal units)
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)
Besides the basic version, it is also possible to run the Hernquist
potential for idealised disk galaxies that follow the approach of
Springel+ 2005. The default Hernquist potential uses a corrected value
for the formulation that improves the match with the NFW (below) with
the same M200 (Nobels+ in prep). Contrary to above, the idealised disk
setup runs with a specified M200, concentration and reduced Hubble
constant that set both the mass and scale length parameter. The reduced
Hubble constant is only used to determine R200.
The parameters of the model are:
`Springel, Di Matteo & Hernquist (2005)
<https://ui.adsabs.harvard.edu/abs/2005MNRAS.361..776S/abstract>`_. This
potential, however, uses a corrected value of the formulation that improves
the match with the NFW profile (below) with the same M200 (Nobels+ in prep).
Contrary to the above (idealizeddisk: 0 setup), the idealised disk setup runs
by specifying one out of :math:`M_{200}`, :math:`V_{200}`, or :math:`R_{200}`,
plus concentration and reduced Hubble constant.
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
......@@ -199,21 +202,52 @@ The parameters of the model are:
concentration: 9.0 # concentration of the Halo
diskfraction: 0.040 # Disk 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)
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
way that the idealised disk part is calculated. This potential uses
exactly the same approach as `Springel, Di Matteo & Hernquist (2005)
<https://ui.adsabs.harvard.edu/abs/2005MNRAS.361..776S/abstract>`_,
this means that ICs generated with the original `MakeNewDisk` code can
be used when using this potential. Contrary to the updated Hernquist
potential (above) it is not possible to have an identically matched
NFW potential.
which means that ICs generated with the original `MakeNewDisk` code can
be used with this potential. Contrary to the updated Hernquist
potential (above), it is not possible to have an identically matched
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``):
......
......@@ -22,7 +22,7 @@ import h5py
import matplotlib.pyplot as plt
from scipy.integrate import odeint
t = np.linspace(0, 40, 1e5)
t = np.linspace(0, 40, 100000)
y0 = [0, 10]
a = 30.0
G = 4.300927e-06
......
......@@ -73,6 +73,9 @@ struct external_potential {
* time to get the time steps */
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
* galaxies */
int usedisk;
......@@ -92,21 +95,18 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep(
const struct phys_const* restrict phys_const,
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
* 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 radius */
/* Calculate the radius */
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 */
const float period = 2.f * M_PI * sqrtf(r) * potential->al *
(1 + r / potential->al) * sqrtgm_inv;
const float period =
2.f * M_PI * sqrtf(r) * (potential->al + r) * potential->sqrtgm_inv;
/* Time-step as a fraction of the cirecular orbital time */
const float time_step = potential->timestep_mult * period;
......@@ -211,7 +211,7 @@ static INLINE void potential_init_backend(
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(
parameter_file, "HernquistPotential:idealizeddisk",
idealized_disk_default);
......@@ -269,7 +269,7 @@ static INLINE void potential_init_backend(
potential->M200 = M200;
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,
"HernquistPotential:concentration");
......@@ -286,7 +286,7 @@ static INLINE void potential_init_backend(
const double b = 2. * cc_inv * cc_inv * (log(1. + cc) - cc / (1. + cc));
/* 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*/
const double R200_inv = 1. / R200;
......@@ -296,9 +296,9 @@ static INLINE void potential_init_backend(
(potential->R200 + potential->al) * R200_inv *
R200_inv * potential->M200;
/* Depending on the disk mass and and the bulge mass the halo
* gets a different mass, because of this we read the fractions
* from the parameter file and calculate the absolute mass*/
/* Depending on the disk mass and and the bulge mass, the halo
* gets a different mass. Because of this, we read the fractions
* from the parameter file and calculate the absolute mass */
const double diskfraction = parser_get_param_double(
parameter_file, "HernquistPotential:diskfraction");
const double bulgefraction = parser_get_param_double(
......@@ -318,12 +318,14 @@ static INLINE void potential_init_backend(
parser_get_param_double(parameter_file, "HernquistPotential:epsilon");
potential->epsilon2 = epsilon * epsilon;
/* Compute the minimal time-step. */
/* This is the circular orbital time at the softened radius */
/* Calculate a common factor in the calculation, i.e. 1/sqrt(GM)*/
const float sqrtgm = sqrtf(phys_const->const_newton_G * potential->mass);
potential->mintime = 2.f * sqrtf(epsilon) * potential->al * M_PI *
(1. + epsilon / potential->al) / sqrtgm *
potential->timestep_mult;
potential->sqrtgm_inv = 1. / sqrtgm;
/* 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;
}
/**
......
......@@ -63,6 +63,9 @@ struct external_potential {
/*! Time-step condition pre-factor, is multiplied times the circular orbital
* time to get the time steps */
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(
const struct phys_const* restrict phys_const,
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
* potential */
const float dx = g->x[0] - potential->x[0];
......@@ -89,11 +90,10 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep(
/* calculate the radius */
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 */
const float period = 2.f * M_PI * sqrtf(r) * potential->al *
(1 + r / potential->al) * sqrtgm_inv;
const float period =
2.f * M_PI * sqrtf(r) * (potential->al + r) * potential->sqrtgm_inv;
/* Time-step as a fraction of the cirecular orbital time */
const float time_step = potential->timestep_mult * period;
......@@ -254,7 +254,7 @@ static INLINE void potential_init_backend(
/* message("M200 = %g, R200 = %g, V200 = %g", M200, R200, V200); */
/* 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(
parameter_file, "HernquistPotential:concentration");
......@@ -262,12 +262,12 @@ static INLINE void potential_init_backend(
const double RS = R200 / concentration;
/* 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)));
/* Depending on the disk mass and and the bulge mass the halo
* gets a different mass, because of this we read the fractions
* from the parameter file and calculate the absolute mass*/
/* Depending on the disk mass and the bulge mass, the halo
* gets a different mass. Because of this, we read the fractions
* from the parameter file and calculate the absolute mass */
const double diskfraction = parser_get_param_double(
parameter_file, "HernquistPotential:diskfraction");
const double bulgefraction = parser_get_param_double(
......@@ -287,12 +287,14 @@ static INLINE void potential_init_backend(
parser_get_param_double(parameter_file, "HernquistPotential:epsilon");
potential->epsilon2 = epsilon * epsilon;
/* Compute the minimal time-step. */
/* This is the circular orbital time at the softened radius */
/* Calculate a common factor in the calculation, i.e. 1/sqrt(GM)*/
const float sqrtgm = sqrtf(phys_const->const_newton_G * potential->mass);
potential->mintime = 2.f * sqrtf(epsilon) * potential->al * M_PI *
(1. + epsilon / potential->al) / sqrtgm *
potential->timestep_mult;
potential->sqrtgm_inv = 1. / sqrtgm;
/* 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(
const struct external_potential* potential) {
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 "
", minimum time = %e timestep multiplier = %e",
potential->x[0], potential->x[1], potential->x[2], potential->mass,
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment