diff --git a/.gitignore b/.gitignore
index a6b0a250fb2818f954403ad6899dc39d0db16919..7ad337771989f94131f7091d4ed55e285582cf4e 100644
--- a/.gitignore
+++ b/.gitignore
@@ -172,6 +172,8 @@ tests/testCooling
tests/testComovingCooling
tests/testHashmap
tests/testNeutrinoCosmology
+tests/testNeutrinoFermiDirac
+tests/testLog
tests/*.png
tests/*.txt
diff --git a/doc/RTD/source/Neutrinos/index.rst b/doc/RTD/source/Neutrinos/index.rst
new file mode 100644
index 0000000000000000000000000000000000000000..df9af9e3ea19853eb4fa9f08ef1156ad9c6c90f5
--- /dev/null
+++ b/doc/RTD/source/Neutrinos/index.rst
@@ -0,0 +1,79 @@
+.. Neutrinos
+ Willem Elbers, 7 April 2021
+
+.. _neutrinos:
+
+Neutrino implementation
+=======================
+
+SWIFT can also accurately model the effects of massive neutrinos in
+cosmological simulations. At the background level, massive neutrinos
+and other relativistic species can be included by specifying their
+number and masses in the cosmology section of the parameter file
+(see :ref:`Parameters_cosmology`).
+
+At the perturbation level, neutrinos can be included as a separate particle
+species (``PartType6``). To facilitate this, SWIFT implements the
+:math:`\delta f` method for shot noise suppression (`Elbers et al. 2020
+<https://ui.adsabs.harvard.edu/abs/2020arXiv201007321E/>`_). The method
+works by statistically weighting the particles during the simulation,
+with weights computed from the Liouville equation using current and
+initial momenta. The method can be activated by specifying
+``Neutrino:use_delta_f`` in the parameter file.
+
+The implementation of the :math:`\delta f` method in SWIFT assumes a
+specific method for generating the initial neutrino momenta (see below).
+This makes it possible to reproduce the initial momentum when it is
+needed without increasing the memory footprint of the neutrino particles.
+If perturbed initial conditions are not needed, the initial momenta can
+be generated internally by specifying ``Neutrino:generate_ics`` in the
+parameter file. This will assign ``PartType6`` particles to each
+neutrino mass specified in the cosmology and generate new velocities
+based on the homogeneous (unperturbed) Fermi-Dirac distribution.
+
+Relativistic Drift
+------------------
+
+At high redshift, neutrino particles move faster than the speed of light
+if the usual Newtonian expressions are used. To rectify this, SWIFT
+implements a relativistic drift correction. In this convention, the
+internal velocity variable (see theory/Cosmology) is
+:math:`v^i=a^2u^i=a^2\dot{x}^i\gamma^{-1}`, where :math:`u^i` is the
+spatial part of the 4-velocity, :math:`a` the scale factor, and
+:math:`x^i` a comoving position vector. The conversion factor to the
+coordinate 3-velocity is :math:`\gamma=ac/\sqrt{a^2c^2+v^2}`. This
+factor is applied to the neutrino particles throughout the simulation.
+
+Generating Fermi-Dirac momenta
+------------------------------
+
+The implementation of the :math:`\delta f` method in SWIFT assumes that
+neutrinos were initially assigned a Fermi-Dirac momentum using the following
+method. Each particle has a fixed 64-bit unsigned integer :math:`\ell` given
+by the particle ID [#f1]_ (plus an optional seed: ``Neutrino:neutrino_seed``).
+This number is transformed into a floating point number :math:`u\in(0,1)`,
+using the following pseudo-code based on splitmix64:
+
+.. code-block:: none
+
+ m = l + 0x9E3779B97f4A7C15
+ m = (m ^ (m >> 30)) * 0xBF58476D1CE4E5B9;
+ m = (m ^ (m >> 27)) * 0x94D049BB133111EB;
+ m = m ^ (m >> 31);
+ u = (m + 0.5) / (UINT64_MAX + 1)
+
+This is subsequently transformed into a Fermi-Dirac momentum
+:math:`q = F^{-1}(u)` by evaluating the quantile function. To generate
+neutrino particle initial conditions with perturbations, one first generates
+momenta from the unperturbed Fermi-Dirac distribution using the above method
+and then applies perturbations in any suitable manner.
+
+When using the :math:`\delta f` method, SWIFT also assumes that ``PartType6``
+particles are assigned to all :math:`N_\nu` massive species present in the
+cosmology, such that the particle with fixed integer :math:`\ell` corresponds
+to species :math:`i = \ell\; \% \;N_\nu\in[0,N_\nu-1]`.
+
+The sampled Fermi-Dirac speeds and neutrino masses are written into the
+snapshot files as ``SampledSpeeds`` and ``MicroscopicMasses``.
+
+.. [#f1] Currently, it is not guaranteed that a particle ID is unique.
diff --git a/doc/RTD/source/ParameterFiles/parameter_description.rst b/doc/RTD/source/ParameterFiles/parameter_description.rst
index 5cdb52cdd74f1870c985d9feab012ec1620fb0c0..34c8f1b68deae0bb030b96202c6662c915156028 100644
--- a/doc/RTD/source/ParameterFiles/parameter_description.rst
+++ b/doc/RTD/source/ParameterFiles/parameter_description.rst
@@ -1481,3 +1481,24 @@ gparts are active after that snapshot output timestep.
.. [#f2] which would translate into a constant :math:`G_N=1.5517771\times10^{-9}~cm^{3}\,g^{-1}\,s^{-2}` if expressed in the CGS system.
+
+Neutrinos
+---------
+
+The ``Neutrino`` section of the parameter file controls the behaviour of
+neutrino particles (``PartType6``). This assumes that massive neutrinos have
+been specified in the ``Cosmology`` section described above. Random
+Fermi-Dirac momenta will be generated if ``generate_ics`` is used. The
+:math:`\delta f` method for shot noise reduction can be activated with
+``use_delta_f``. Finally, a random seed for the Fermi-Dirac momenta can
+be set with ``neutrino_seed``.
+
+For mode details on the neutrino implementation, refer to :ref:`Neutrinos`.
+A complete specification of the model looks like
+
+.. code:: YAML
+
+ Neutrino:
+ generate_ics: 1 # Replace neutrino particle velocities with random Fermi-Dirac momenta at the start
+ use_delta_f: 1 # Use the delta-f method for shot noise reduction
+ neutrino_seed: 1234 # A random seed used for the Fermi-Dirac momenta
diff --git a/doc/RTD/source/conf.py b/doc/RTD/source/conf.py
index 35d5b592800cc8556d074c6a5bbf68d9ec21e54c..f0e354eb06f49d437584f9f5fa2cec42bea5562d 100644
--- a/doc/RTD/source/conf.py
+++ b/doc/RTD/source/conf.py
@@ -19,7 +19,7 @@
# -- Project information -----------------------------------------------------
project = 'SWIFT: SPH With Inter-dependent Fine-grained Tasking'
-copyright = '2014-2020, SWIFT Collaboration'
+copyright = '2014-2021, SWIFT Collaboration'
author = 'SWIFT Team'
# The short X.Y version
diff --git a/doc/RTD/source/index.rst b/doc/RTD/source/index.rst
index 381f857728399024afe6a9b97c6e50935ea6d80d..d09c27fe5707f0ba576a5bfb7c7c9d2e11df1a4b 100644
--- a/doc/RTD/source/index.rst
+++ b/doc/RTD/source/index.rst
@@ -29,6 +29,7 @@ difference is the parameter file that will need to be adapted for SWIFT.
LineOfSights/index
EquationOfState/index
ExternalPotentials/index
+ Neutrinos/index
NewOption/index
Task/index
AnalysisTools/index
diff --git a/examples/Cosmology/NeutrinoCosmo/neutrino_cosmo.yml b/examples/Cosmology/NeutrinoCosmo/neutrino_cosmo.yml
index 9acf793e6f96f5b2f47f4814cd5f341c1c0483d7..f3d682d055427e8fdd3f2909371918a4e309f53c 100644
--- a/examples/Cosmology/NeutrinoCosmo/neutrino_cosmo.yml
+++ b/examples/Cosmology/NeutrinoCosmo/neutrino_cosmo.yml
@@ -43,6 +43,10 @@ InitialConditions:
file_name: particles.hdf5
periodic: 1
+Neutrino:
+ use_delta_f: 1 # Use the delta-f method for shot noise reduction
+ generate_ics: 1 # Generate new Fermi-Dirac velocities for the particles
+
Gravity:
mesh_side_length: 32
MAC: geometric
diff --git a/examples/main.c b/examples/main.c
index 1f96cd49239f0e44a89a17cc5ee6062b6ad2e755..37b11d4d0379405dfe765a3e1638e365a94f88a5 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -97,6 +97,7 @@ int main(int argc, char *argv[]) {
struct hydro_props hydro_properties;
struct stars_props stars_properties;
struct sink_props sink_properties;
+ struct neutrino_props neutrino_properties;
struct feedback_props feedback_properties;
struct rt_props rt_properties;
struct entropy_floor_properties entropy_floor;
@@ -1268,6 +1269,12 @@ int main(int argc, char *argv[]) {
/* Do we have neutrino DM particles? */
const int with_neutrinos = N_total[swift_type_neutrino] > 0;
+ /* Initialise the neutrino properties if we have neutrino particles */
+ bzero(&neutrino_properties, sizeof(struct neutrino_props));
+ if (with_neutrinos)
+ neutrino_props_init(&neutrino_properties, &prog_const, &us, params,
+ &cosmo);
+
/* Initialize the space with these data. */
if (myrank == 0) clocks_gettime(&tic);
space_init(&s, params, &cosmo, dim, &hydro_properties, parts, gparts, sinks,
@@ -1422,9 +1429,10 @@ int main(int argc, char *argv[]) {
N_total[swift_type_neutrino], engine_policies, talking,
&reparttype, &us, &prog_const, &cosmo, &hydro_properties,
&entropy_floor, &gravity_properties, &stars_properties,
- &black_holes_properties, &sink_properties, &feedback_properties,
- &rt_properties, &mesh, &potential, &cooling_func, &starform,
- &chemistry, &fof_properties, &los_properties);
+ &black_holes_properties, &sink_properties, &neutrino_properties,
+ &feedback_properties, &rt_properties, &mesh, &potential,
+ &cooling_func, &starform, &chemistry, &fof_properties,
+ &los_properties);
engine_config(/*restart=*/0, /*fof=*/0, &e, params, nr_nodes, myrank,
nr_threads, nr_pool_threads, with_aff, talking, restart_file);
diff --git a/examples/main_fof.c b/examples/main_fof.c
index 9002dde186f24be76100535c189c1f019e36c8ee..c7f40980b6b7d3311b9c24b8f46eb2000645fa2f 100644
--- a/examples/main_fof.c
+++ b/examples/main_fof.c
@@ -642,7 +642,7 @@ int main(int argc, char *argv[]) {
engine_policies, talking, &reparttype, &us, &prog_const, &cosmo,
/*hydro_properties=*/NULL, /*entropy_floor=*/NULL, &gravity_properties,
/*stars_properties=*/NULL, /*black_holes_properties=*/NULL,
- /*sink_properties=*/NULL,
+ /*sink_properties=*/NULL, /*neutrino_properties=*/NULL,
/*feedback_properties=*/NULL, /*rt_properties=*/NULL, &mesh,
/*potential=*/NULL,
/*cooling_func=*/NULL, /*starform=*/NULL, /*chemistry=*/NULL,
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 2f50e067f9501b07f63de4559837dad267f021ff..fda3c046a92eab8982ad618d3b675f63db9acf9e 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -638,6 +638,12 @@ EAGLEAGN:
merger_max_distance_ratio: 3.0 # Maximal distance over which two BHs can merge, in units of the softening length.
minimum_timestep_Myr: 1.0 # Minimum time-step size for BHs in Mega-years. The time-step size is computed based on the accretion rate and this serves as a minimum. Defaults to FLT_MAX.
+# Parameters related to the neutrino particles ----------------------------------
+Neutrino:
+ use_delta_f: 1 # Use the delta-f method for shot noise reduction
+ generate_ics: 0 # Automatically generate neutrino velocities at the start of the run
+ neutrino_seed: 1234 # Seed used in the delta-f weighting step for Fermi-Dirac momentum generation
+
# Parameters related to the sink particles ---------------------------------------
# Default sink particles
diff --git a/src/Makefile.am b/src/Makefile.am
index 7f127064436d052924c452447ecde44dcd80821f..f6013014caeab4bd584ec7f1f083b4f7e5071475 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -177,6 +177,8 @@ AM_SOURCES += chemistry.c cosmology.c mesh_gravity.c velociraptor_interface.c
AM_SOURCES += output_list.c velociraptor_dummy.c logger_io.c memuse.c mpiuse.c memuse_rnodes.c
AM_SOURCES += fof.c fof_catalogue_io.c
AM_SOURCES += hashmap.c pressure_floor.c logger_history.c
+AM_SOURCES += runner_neutrino.c
+AM_SOURCES += neutrino/Default/fermi_dirac.c
AM_SOURCES += $(QLA_COOLING_SOURCES)
AM_SOURCES += $(EAGLE_COOLING_SOURCES) $(EAGLE_FEEDBACK_SOURCES)
AM_SOURCES += $(GRACKLE_COOLING_SOURCES) $(GEAR_FEEDBACK_SOURCES)
@@ -186,7 +188,7 @@ AM_SOURCES += $(STAR_FORMATION_LOGGER)
# Include files for distribution, not installation.
nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h
-nobase_noinst_HEADERS += gravity_iact.h kernel_long_gravity.h vector.h accumulate.h cache.h exp.h
+nobase_noinst_HEADERS += gravity_iact.h kernel_long_gravity.h vector.h accumulate.h cache.h exp.h log.h
nobase_noinst_HEADERS += runner_doiact_nosort.h runner_doiact_hydro.h runner_doiact_stars.h runner_doiact_black_holes.h runner_doiact_grav.h
nobase_noinst_HEADERS += runner_doiact_functions_hydro.h runner_doiact_functions_stars.h runner_doiact_functions_black_holes.h
nobase_noinst_HEADERS += runner_doiact_functions_limiter.h runner_doiact_limiter.h units.h intrinsics.h minmax.h
@@ -400,7 +402,9 @@ nobase_noinst_HEADERS += pressure_floor/GEAR/pressure_floor_struct.h pressure_fl
nobase_noinst_HEADERS += sink/Default/sink.h sink/Default/sink_io.h sink/Default/sink_part.h sink/Default/sink_properties.h
nobase_noinst_HEADERS += sink/Default/sink_iact.h
nobase_noinst_HEADERS += sink.h sink_io.h sink_properties.h
-nobase_noinst_HEADERS += neutrino/neutrino.h neutrino/relativity.h
+nobase_noinst_HEADERS += neutrino.h neutrino_properties.h neutrino_io.h
+nobase_noinst_HEADERS += neutrino/Default/neutrino.h neutrino/Default/relativity.h neutrino/Default/fermi_dirac.h
+nobase_noinst_HEADERS += neutrino/Default/neutrino_properties.h neutrino/Default/neutrino_io.h
# Sources and special flags for the gravity library
libgrav_la_SOURCES = runner_doiact_grav.c
diff --git a/src/cell_drift.c b/src/cell_drift.c
index 1a128fe4caa1441a6f992b8458f453bd14fe5dfd..f541b4b89f5d0226a98d93949d958eaaf82ab55c 100644
--- a/src/cell_drift.c
+++ b/src/cell_drift.c
@@ -31,7 +31,7 @@
#include "feedback.h"
#include "gravity.h"
#include "multipole.h"
-#include "neutrino/relativity.h"
+#include "neutrino.h"
#include "pressure_floor.h"
#include "rt.h"
#include "star_formation.h"
@@ -281,7 +281,6 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) {
const double a = e->cosmology->a;
const double c_vel = e->physical_constants->const_speed_light_c;
const int with_neutrinos = e->s->with_neutrinos;
- const int with_relat_drift = with_neutrinos && (a < e->cosmology->a_nu_nr);
/* Drift irrespective of cell flags? */
force = (force || cell_get_flag(c, cell_flag_do_grav_drift));
@@ -344,7 +343,7 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) {
/* Relativistic drift correction for neutrinos */
double dt_drift_k = dt_drift;
- if (with_relat_drift && gp->type == swift_type_neutrino) {
+ if (with_neutrinos && gp->type == swift_type_neutrino) {
dt_drift_k *= relativistic_drift_factor(gp->v_full, a, c_vel);
}
diff --git a/src/cell_grav.h b/src/cell_grav.h
index ec690966efa7d3b79b0f085d52a190499eddbcd5..b0780c7d570694be0670960b2b3f26b6d62c4bfd 100644
--- a/src/cell_grav.h
+++ b/src/cell_grav.h
@@ -77,6 +77,9 @@ struct cell_grav {
/*! The task to end the force calculation */
struct task *end_force;
+ /*! Task for weighting neutrino particles */
+ struct task *neutrino_weight;
+
/*! Minimum end of (integer) time step in this cell for gravity tasks. */
integertime_t ti_end_min;
diff --git a/src/cell_unskip.c b/src/cell_unskip.c
index d4eb0525f97a74a1d184a7ca81ebb3ea32ed9361..d229fbffe550353cc3f3c0275044450299d01dcb 100644
--- a/src/cell_unskip.c
+++ b/src/cell_unskip.c
@@ -1961,6 +1961,8 @@ int cell_unskip_gravity_tasks(struct cell *c, struct scheduler *s) {
if (c->grav.down_in != NULL) scheduler_activate(s, c->grav.down_in);
if (c->grav.long_range != NULL) scheduler_activate(s, c->grav.long_range);
if (c->grav.end_force != NULL) scheduler_activate(s, c->grav.end_force);
+ if (c->grav.neutrino_weight != NULL)
+ scheduler_activate(s, c->grav.neutrino_weight);
#ifdef WITH_LOGGER
if (c->logger != NULL) scheduler_activate(s, c->logger);
#endif
diff --git a/src/common_io.c b/src/common_io.c
index e1ffd211a37d57c9a291f60ba70b953ae28175ef..53694875ff3c17b4e084b3a7b89647e33c8512df 100644
--- a/src/common_io.c
+++ b/src/common_io.c
@@ -42,6 +42,7 @@
#include "fof_io.h"
#include "gravity_io.h"
#include "hydro_io.h"
+#include "neutrino_io.h"
#include "particle_splitting.h"
#include "rt_io.h"
#include "sink_io.h"
@@ -1641,6 +1642,27 @@ void io_select_dm_fields(const struct gpart* const gparts,
}
}
+/**
+ * @brief Select the fields to write to snapshots for the neutrino particles.
+ *
+ * @param gparts The #gpart's
+ * @param with_fof Are we running FoF?
+ * @param with_stf Are we running with structure finding?
+ * @param e The #engine (to access scheme properties).
+ * @param num_fields (return) The number of fields to write.
+ * @param list (return) The list of fields to write.
+ */
+void io_select_neutrino_fields(
+ const struct gpart* const gparts,
+ const struct velociraptor_gpart_data* gpart_group_data, const int with_fof,
+ const int with_stf, const struct engine* const e, int* const num_fields,
+ struct io_props* const list) {
+
+ darkmatter_write_particles(gparts, list, num_fields);
+
+ *num_fields += neutrino_write_particles(gparts, list + *num_fields);
+}
+
/**
* @brief Select the fields to write to snapshots for the sink particles.
*
diff --git a/src/common_io.h b/src/common_io.h
index 1311e77e439ee75f57ce50e1f8ac5dc454c54a91..879a927fc6f04d857a7b1b9e351bfa2c4a60f5ec 100644
--- a/src/common_io.h
+++ b/src/common_io.h
@@ -227,6 +227,12 @@ void io_select_dm_fields(const struct gpart* const gparts,
const struct engine* const e, int* const num_fields,
struct io_props* const list);
+void io_select_neutrino_fields(
+ const struct gpart* const gparts,
+ const struct velociraptor_gpart_data* gpart_group_data, const int with_fof,
+ const int with_stf, const struct engine* const e, int* const num_fields,
+ struct io_props* const list);
+
void io_select_sink_fields(const struct sink* const sinks,
const int with_cosmology, const int with_fof,
const int with_stf, const struct engine* const e,
diff --git a/src/common_io_fields.c b/src/common_io_fields.c
index c7d8d4bf300c714350fcc60fe3c2520ce67056c2..00fa41aaf6c0f6fc824cbd7feea2c1a9778133e7 100644
--- a/src/common_io_fields.c
+++ b/src/common_io_fields.c
@@ -85,15 +85,16 @@ int io_get_ptype_fields(const int ptype, struct io_props* list,
case swift_type_dark_matter:
case swift_type_dark_matter_background:
- case swift_type_neutrino:
io_select_dm_fields(NULL, NULL, with_fof, with_stf, /*e=*/NULL,
&num_fields, list);
break;
-
+ case swift_type_neutrino:
+ io_select_neutrino_fields(NULL, NULL, with_fof, with_stf, /*e=*/NULL,
+ &num_fields, list);
+ break;
case swift_type_stars:
io_select_star_fields(NULL, with_cosmology, with_fof, with_stf,
- /*with_rt=*/1,
- /*e=*/NULL, &num_fields, list);
+ /*with_rt=*/1, /*e=*/NULL, &num_fields, list);
break;
case swift_type_sink:
diff --git a/src/cosmology.c b/src/cosmology.c
index 6603d5c1e77fc8bac7553087ed83490cfbad455b..23aa2877cdae08d9fc374b7ee0de3f29f11a58f9 100644
--- a/src/cosmology.c
+++ b/src/cosmology.c
@@ -801,7 +801,6 @@ void cosmology_init(struct swift_params *params, const struct unit_system *us,
c->Omega_nu_0 = 0.;
c->Omega_nu = 0.;
c->N_eff = 0.;
- c->a_nu_nr = 0.;
c->neutrino_density_early_table = NULL;
c->neutrino_density_late_table = NULL;
@@ -862,15 +861,6 @@ void cosmology_init(struct swift_params *params, const struct unit_system *us,
for (int i = 0; i < c->N_nu; i++) {
M_eV_min = fmin(M_eV_min, c->M_nu_eV[i]);
}
-
- /* Time when the relativistic drift correction is below 1% for all nu's */
- if (c->N_nu > 0 && M_eV_min > 0.) {
- const double five_sigma_velocity = 20.0 * c->T_nu_0_eV / M_eV_min;
- const double approx_treshold = 0.14249; // 1.0 / sqrt(1 + x^2) = 0.99
- c->a_nu_nr = five_sigma_velocity / approx_treshold;
- } else {
- c->a_nu_nr = 0.;
- }
}
/* Cold dark matter density */
@@ -927,7 +917,6 @@ void cosmology_init_no_cosmo(struct cosmology *c) {
c->N_nu = 0;
c->N_ur = 0.;
c->N_eff = 0.;
- c->a_nu_nr = 0.;
c->a_begin = 1.;
c->a_end = 1.;
@@ -1359,7 +1348,6 @@ void cosmology_write_model(hid_t h_grp, const struct cosmology *c) {
io_write_attribute_d(h_grp, "N_eff", c->N_eff);
io_write_attribute_d(h_grp, "N_ur", c->N_ur);
io_write_attribute_d(h_grp, "N_nu", c->N_nu);
- io_write_attribute_d(h_grp, "a_nu_nr", c->a_nu_nr);
if (c->N_nu > 0) {
io_write_attribute(h_grp, "M_nu_eV", DOUBLE, c->M_nu_eV, c->N_nu);
io_write_attribute(h_grp, "deg_nu", DOUBLE, c->deg_nu, c->N_nu);
diff --git a/src/cosmology.h b/src/cosmology.h
index 2f15a696c6781c55b022060e18e709d5c2b0c4b4..47a17e26b50e527ba8b0bea4070618096021057c 100644
--- a/src/cosmology.h
+++ b/src/cosmology.h
@@ -198,9 +198,6 @@ struct cosmology {
/*! Degeneracy of each massive neutrino species */
double *deg_nu;
- /*! Scale-factor after which all neutrinos move non-relativistically */
- double a_nu_nr;
-
/*! Log of starting expansion factor for neutrino interpolation tables */
double log_a_long_begin;
diff --git a/src/distributed_io.c b/src/distributed_io.c
index 8ceec30b20fb6fddcf0825d40d73dbb74266d861..15451fd85fa856675c6205c135a31772563c79a9 100644
--- a/src/distributed_io.c
+++ b/src/distributed_io.c
@@ -703,8 +703,8 @@ void write_output_distributed(struct engine* e,
gpart_group_data_written, Ntot, Ndm_neutrino, with_stf);
/* Select the fields to write */
- io_select_dm_fields(gparts_written, gpart_group_data_written, with_fof,
- with_stf, e, &num_fields, list);
+ io_select_neutrino_fields(gparts_written, gpart_group_data_written,
+ with_fof, with_stf, e, &num_fields, list);
} break;
case swift_type_sink: {
diff --git a/src/engine.c b/src/engine.c
index f59f328e33f4d0b63156be3139ea006f8c44a653..36b8f7048650b8197d623b146c692c29609b52d5 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -77,7 +77,8 @@
#include "minmax.h"
#include "mpiuse.h"
#include "multipole_struct.h"
-#include "neutrino/neutrino.h"
+#include "neutrino.h"
+#include "neutrino_properties.h"
#include "output_list.h"
#include "output_options.h"
#include "partition.h"
@@ -1572,6 +1573,7 @@ void engine_skip_force_and_kick(struct engine *e) {
t->type == task_type_bh_swallow_ghost3 || t->type == task_type_bh_in ||
t->type == task_type_bh_out || t->type == task_type_rt_ghost1 ||
t->type == task_type_rt_ghost2 || t->type == task_type_rt_tchem ||
+ t->type == task_type_neutrino_weight ||
t->subtype == task_subtype_force ||
t->subtype == task_subtype_limiter ||
t->subtype == task_subtype_gradient ||
@@ -2762,6 +2764,7 @@ void engine_unpin(void) {
* @param stars The #stars_props used for this run.
* @param black_holes The #black_holes_props used for this run.
* @param sinks The #sink_props used for this run.
+ * @param neutrinos The #neutrino_props used for this run.
* @param feedback The #feedback_props used for this run.
* @param mesh The #pm_mesh used for the long-range periodic forces.
* @param potential The properties of the external potential.
@@ -2782,7 +2785,8 @@ void engine_init(
const struct entropy_floor_properties *entropy_floor,
struct gravity_props *gravity, const struct stars_props *stars,
const struct black_holes_props *black_holes, const struct sink_props *sinks,
- struct feedback_props *feedback, struct rt_props *rt, struct pm_mesh *mesh,
+ const struct neutrino_props *neutrinos, struct feedback_props *feedback,
+ struct rt_props *rt, struct pm_mesh *mesh,
const struct external_potential *potential,
struct cooling_function_data *cooling_func,
const struct star_formation *starform,
@@ -2873,6 +2877,7 @@ void engine_init(
e->stars_properties = stars;
e->black_holes_properties = black_holes;
e->sink_properties = sinks;
+ e->neutrino_properties = neutrinos;
e->mesh = mesh;
e->external_potential = potential;
e->cooling_func = cooling_func;
@@ -3358,6 +3363,7 @@ void engine_struct_dump(struct engine *e, FILE *stream) {
rt_struct_dump(e->rt_props, stream);
black_holes_struct_dump(e->black_holes_properties, stream);
sink_struct_dump(e->sink_properties, stream);
+ neutrino_struct_dump(e->neutrino_properties, stream);
chemistry_struct_dump(e->chemistry, stream);
#ifdef WITH_FOF
fof_struct_dump(e->fof_properties, stream);
@@ -3487,6 +3493,11 @@ void engine_struct_restore(struct engine *e, FILE *stream) {
sink_struct_restore(sink_properties, stream);
e->sink_properties = sink_properties;
+ struct neutrino_props *neutrino_properties =
+ (struct neutrino_props *)malloc(sizeof(struct neutrino_props));
+ neutrino_struct_restore(neutrino_properties, stream);
+ e->neutrino_properties = neutrino_properties;
+
struct chemistry_global_data *chemistry =
(struct chemistry_global_data *)malloc(
sizeof(struct chemistry_global_data));
diff --git a/src/engine.h b/src/engine.h
index d0afe37fb209d9926a2d8cc7575434042c4092bf..0de4031568dee18e580ab74a6e668bd9bce15ce5 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -458,6 +458,9 @@ struct engine {
/* Properties of the sink model */
const struct sink_props *sink_properties;
+ /* Properties of the neutrino model */
+ const struct neutrino_props *neutrino_properties;
+
/* Properties of the self-gravity scheme */
struct gravity_props *gravity_properties;
@@ -589,7 +592,8 @@ void engine_init(
const struct entropy_floor_properties *entropy_floor,
struct gravity_props *gravity, const struct stars_props *stars,
const struct black_holes_props *black_holes, const struct sink_props *sinks,
- struct feedback_props *feedback, struct rt_props *rt, struct pm_mesh *mesh,
+ const struct neutrino_props *neutrinos, struct feedback_props *feedback,
+ struct rt_props *rt, struct pm_mesh *mesh,
const struct external_potential *potential,
struct cooling_function_data *cooling_func,
const struct star_formation *starform,
diff --git a/src/engine_maketasks.c b/src/engine_maketasks.c
index 3fd36a54f703d01a53367ca61f2c147b3ed753f1..6e3ba60a5af82f06786e8dd5dae0004695e7d4ec 100644
--- a/src/engine_maketasks.c
+++ b/src/engine_maketasks.c
@@ -50,6 +50,7 @@
#include "debug.h"
#include "error.h"
#include "feedback.h"
+#include "neutrino_properties.h"
#include "proxy.h"
#include "rt_properties.h"
#include "timers.h"
@@ -1004,6 +1005,13 @@ void engine_make_hierarchical_tasks_common(struct engine *e, struct cell *c) {
c->kick2 = scheduler_addtask(s, task_type_kick2, task_subtype_none, 0, 0,
c, NULL);
+ /* Weighting task for neutrinos after the last kick */
+ if (e->neutrino_properties->use_delta_f) {
+ c->grav.neutrino_weight = scheduler_addtask(
+ s, task_type_neutrino_weight, task_subtype_none, 0, 0, c, NULL);
+ scheduler_addunlock(s, c->kick1, c->grav.neutrino_weight);
+ }
+
#if defined(WITH_LOGGER)
struct task *kick2_or_logger;
if (with_logger) {
diff --git a/src/engine_marktasks.c b/src/engine_marktasks.c
index 0ba9e7dc1b763d258cddd67e97ca014af148d264..ee192fe0f27eeac8bf2b4656476c70f9250ead9a 100644
--- a/src/engine_marktasks.c
+++ b/src/engine_marktasks.c
@@ -1243,6 +1243,13 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
if (cell_is_active_gravity(t->ci, e)) scheduler_activate(s, t);
}
+ /* Activate the weighting task for neutrinos */
+ else if (t_type == task_type_neutrino_weight) {
+ if (cell_is_active_gravity(t->ci, e)) {
+ scheduler_activate(s, t);
+ }
+ }
+
/* Kick ? */
else if (t_type == task_type_kick1 || t_type == task_type_kick2) {
diff --git a/src/log.h b/src/log.h
new file mode 100644
index 0000000000000000000000000000000000000000..6bee9643371323ec9984a11af5f7bf3df589f1e5
--- /dev/null
+++ b/src/log.h
@@ -0,0 +1,76 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2021 Willem Elbers (whe@willemelbers.com)
+ *
+ * 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_OPTIMIZED_LOG_H
+#define SWIFT_OPTIMIZED_LOG_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Local headers. */
+#include "inline.h"
+
+/* Standard headers */
+#include <math.h>
+#include <stdint.h>
+
+/**
+ * @brief Compute the natural logarithm of a number.
+ *
+ * This function has a maximum absolute error of 5e-3 over the domain
+ * [exp(-32.), exp(32.)] and a maximum relative error of 5e-3 over the two
+ * intervals [exp(-32.), exp(-0.3)] and [exp(0.3), exp(32.)]
+ *
+ * @param x The argument of the log
+ */
+__attribute__((always_inline, const)) INLINE static float optimized_logf(
+ float val) {
+ union {
+ int32_t i;
+ float f;
+ } e;
+ e.f = val;
+
+ /* Isolate the exponent */
+ float log_2 = (float)(((e.i >> 23) & 255) - 128);
+ e.i &= ~(255 << 23);
+ e.i += 127 << 23;
+ /* Approximation based on https://stackoverflow.com/a/28730362 comment */
+ log_2 += ((-0.34484843f) * e.f + 2.02466578f) * e.f - 0.67487759f;
+ /* Return the natural log */
+ return log_2 * 0.693147181f;
+}
+
+/**
+ * @brief Compute the base-10 logarithm of a number.
+ *
+ * This function has a maximum absolute error of 5e-3 over the domain
+ * [exp(-32.), exp(32.)] and a maximum relative error of 5e-3 over the two
+ * intervals [exp(-32.), exp(-0.3)] and [exp(0.3), exp(32.)]
+ *
+ * @param x The argument of the log
+ */
+__attribute__((always_inline, const)) INLINE static float optimized_log10f(
+ float val) {
+ /* Compute the natural log */
+ float log_e = optimized_logf(val);
+ /* Return the base-10 log */
+ return log_e * 0.434294482f;
+}
+
+#endif /* SWIFT_OPTIMIZED_LOG_H */
diff --git a/src/neutrino.h b/src/neutrino.h
new file mode 100644
index 0000000000000000000000000000000000000000..a19261b3a619a62fb1aae8892aaf36e83ca1eebb
--- /dev/null
+++ b/src/neutrino.h
@@ -0,0 +1,28 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2021 Willem Elbers (whe@willemelbers.com)
+ *
+ * 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_NEUTRINO_H
+#define SWIFT_NEUTRINO_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Select the correct neutrino model */
+#include "./neutrino/Default/neutrino.h"
+
+#endif
diff --git a/src/neutrino/Default/fermi_dirac.c b/src/neutrino/Default/fermi_dirac.c
new file mode 100644
index 0000000000000000000000000000000000000000..dfd8e34d21a4de49124b6b8e8b5785824cc2a00f
--- /dev/null
+++ b/src/neutrino/Default/fermi_dirac.c
@@ -0,0 +1,327 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2021 Willem Elbers (whe@willemelbers.com)
+ *
+ * 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/>.
+ *
+ ******************************************************************************/
+
+/**
+ * @file fermi_dirac.c
+ *
+ * @brief Fast generation of pseudo-random numbers from a Fermi-Dirac
+ * distribution using numerical inversion (Hormann & Leydold, 2003). The
+ * spline tables were generated with AnyRNG. For more details, refer
+ * to https://github.com/wullm/AnyRNG.
+ */
+
+/* This object's header. */
+#include "fermi_dirac.h"
+
+/* Standard headers */
+#include <math.h>
+
+/* Fast optimized logarithm */
+#include "../../log.h"
+
+/* Cubic spline coefficients */
+struct spline {
+ double a0, a1, a2, a3;
+};
+
+/**
+ * @brief Interpolation and search tables for the Fermi-Dirac distribution
+ */
+struct anyrng {
+ /*! Number of intervals on which the interpolation is defined */
+ int intervalN;
+
+ /*! Endpoints of the intervals */
+ double *endpoints;
+
+ /*! Cubic splines of the Fermi-Dirac quantile function */
+ struct spline *splines;
+
+ /*! Length of the look up tables */
+ int tablelen;
+
+ /*! Search table to look up enclosing intervals for small u */
+ int *index_table_a;
+
+ /*! Search table to look up enclosing intervals for intermediate u */
+ int *index_table_b;
+
+ /*! Search table to look up enclosing intervals for large u */
+ int *index_table_c;
+};
+
+static double endpoints[118] = {
+ 0.000000e+00, 8.812220e-16, 3.490090e-15, 1.764874e-14, 1.080345e-13,
+ 7.475047e-13, 5.544411e-12, 4.267185e-11, 3.346975e-10, 2.649997e-09,
+ 2.107211e-08, 7.092068e-08, 1.677779e-07, 5.644860e-07, 1.334411e-06,
+ 2.599613e-06, 4.480999e-06, 1.057018e-05, 2.054614e-05, 3.533463e-05,
+ 5.584298e-05, 8.296016e-05, 1.175567e-04, 1.604850e-04, 2.125789e-04,
+ 2.746542e-04, 4.319210e-04, 6.384483e-04, 9.001106e-04, 1.222496e-03,
+ 1.610910e-03, 2.070375e-03, 2.605637e-03, 3.221165e-03, 3.921158e-03,
+ 4.709548e-03, 5.590000e-03, 6.565923e-03, 8.816531e-03, 1.148360e-02,
+ 1.458549e-02, 1.813688e-02, 2.214891e-02, 2.662937e-02, 3.158292e-02,
+ 3.701122e-02, 4.291319e-02, 4.928517e-02, 5.612111e-02, 6.341282e-02,
+ 7.115008e-02, 8.791173e-02, 1.062920e-01, 1.261570e-01, 1.473582e-01,
+ 1.697373e-01, 1.931305e-01, 2.173724e-01, 2.422994e-01, 2.677516e-01,
+ 2.935760e-01, 3.196273e-01, 3.457697e-01, 3.718776e-01, 3.978360e-01,
+ 4.235411e-01, 4.488999e-01, 4.738305e-01, 5.221310e-01, 5.679885e-01,
+ 6.110917e-01, 6.512521e-01, 6.883826e-01, 7.224780e-01, 7.535962e-01,
+ 7.818432e-01, 8.073587e-01, 8.303052e-01, 8.508587e-01, 8.692016e-01,
+ 8.855172e-01, 8.999850e-01, 9.127781e-01, 9.240609e-01, 9.339876e-01,
+ 9.427015e-01, 9.503348e-01, 9.570084e-01, 9.628322e-01, 9.679056e-01,
+ 9.723182e-01, 9.761502e-01, 9.794730e-01, 9.823505e-01, 9.848391e-01,
+ 9.869886e-01, 9.888431e-01, 9.904414e-01, 9.918172e-01, 9.930005e-01,
+ 9.948897e-01, 9.962793e-01, 9.972980e-01, 9.980425e-01, 9.985850e-01,
+ 9.989794e-01, 9.992652e-01, 9.994720e-01, 9.996213e-01, 9.998061e-01,
+ 9.999013e-01, 9.999500e-01, 9.999748e-01, 9.999937e-01, 9.999984e-01,
+ 9.999996e-01, 9.999999e-01, 1.000000e+00};
+static struct spline splines[118] = {
+ {1.000000e-05, 3.177853e-05, -3.440759e-05, 1.454999e-05},
+ {2.192092e-05, 1.957877e-05, -1.160959e-05, 3.951739e-06},
+ {3.384185e-05, 4.458279e-05, -3.298529e-05, 1.224435e-05},
+ {5.768370e-05, 9.796089e-05, -8.223074e-05, 3.195355e-05},
+ {1.053674e-04, 2.077194e-04, -1.865719e-04, 7.421998e-05},
+ {2.007348e-04, 4.293444e-04, -3.993851e-04, 1.607755e-04},
+ {3.914696e-04, 8.738366e-04, -8.274556e-04, 3.350886e-04},
+ {7.729391e-04, 1.763374e-03, -1.684705e-03, 6.842694e-04},
+ {1.535878e-03, 3.542204e-03, -3.398798e-03, 1.382473e-03},
+ {3.061757e-03, 7.097573e-03, -6.822804e-03, 2.776987e-03},
+ {6.113513e-03, 4.824439e-03, -2.643433e-03, 8.707511e-04},
+ {9.165270e-03, 4.177172e-03, -1.553626e-03, 4.282101e-04},
+ {1.221703e-02, 9.643782e-03, -5.278698e-03, 1.738429e-03},
+ {1.832054e-02, 8.348620e-03, -3.098591e-03, 8.534846e-04},
+ {2.442405e-02, 7.742955e-03, -2.147031e-03, 5.075889e-04},
+ {3.052757e-02, 7.392976e-03, -1.625943e-03, 3.364794e-04},
+ {3.663108e-02, 1.666989e-02, -6.155448e-03, 1.692581e-03},
+ {4.883811e-02, 1.546030e-02, -4.257164e-03, 1.003894e-03},
+ {6.104513e-02, 1.476135e-02, -3.217990e-03, 6.636676e-04},
+ {7.325216e-02, 1.430652e-02, -2.570199e-03, 4.707042e-04},
+ {8.545918e-02, 1.398707e-02, -2.130803e-03, 3.507572e-04},
+ {9.766621e-02, 1.375043e-02, -1.814549e-03, 2.711427e-04},
+ {1.098732e-01, 1.356811e-02, -1.576708e-03, 2.156231e-04},
+ {1.220803e-01, 1.342333e-02, -1.391687e-03, 1.753831e-04},
+ {1.342873e-01, 1.330557e-02, -1.243846e-03, 1.452992e-04},
+ {1.464943e-01, 2.851127e-02, -5.005621e-03, 9.084052e-04},
+ {1.709084e-01, 2.787360e-02, -4.132288e-03, 6.727367e-04},
+ {1.953224e-01, 2.740090e-02, -3.503576e-03, 5.167249e-04},
+ {2.197365e-01, 2.703641e-02, -3.030574e-03, 4.082211e-04},
+ {2.441505e-01, 2.674669e-02, -2.662433e-03, 3.297933e-04},
+ {2.685646e-01, 2.651082e-02, -2.368087e-03, 2.713228e-04},
+ {2.929786e-01, 2.631498e-02, -2.127539e-03, 2.266089e-04},
+ {3.173927e-01, 2.614973e-02, -1.927362e-03, 1.916808e-04},
+ {3.418067e-01, 2.600838e-02, -1.758224e-03, 1.639006e-04},
+ {3.662208e-01, 2.588603e-02, -1.613441e-03, 1.414606e-04},
+ {3.906348e-01, 2.577907e-02, -1.488110e-03, 1.230891e-04},
+ {4.150489e-01, 2.568473e-02, -1.378552e-03, 1.078702e-04},
+ {4.394629e-01, 5.362050e-02, -5.511167e-03, 7.187774e-04},
+ {4.882911e-01, 5.303602e-02, -4.778531e-03, 5.706129e-04},
+ {5.371192e-01, 5.255849e-02, -4.191326e-03, 4.609411e-04},
+ {5.859473e-01, 5.216054e-02, -3.710147e-03, 3.777099e-04},
+ {6.347754e-01, 5.182344e-02, -3.308556e-03, 3.132206e-04},
+ {6.836035e-01, 5.153391e-02, -2.968181e-03, 2.623720e-04},
+ {7.324316e-01, 5.128230e-02, -2.675876e-03, 2.216775e-04},
+ {7.812597e-01, 5.106141e-02, -2.421999e-03, 1.886903e-04},
+ {8.300878e-01, 5.086577e-02, -2.199319e-03, 1.616529e-04},
+ {8.789159e-01, 5.069114e-02, -2.002315e-03, 1.392774e-04},
+ {9.277440e-01, 5.053420e-02, -1.826699e-03, 1.206032e-04},
+ {9.765721e-01, 5.039230e-02, -1.669092e-03, 1.049023e-04},
+ {1.025400e+00, 5.026329e-02, -1.526795e-03, 9.161505e-05},
+ {1.074228e+00, 1.028683e-01, -5.815751e-03, 6.036140e-04},
+ {1.171885e+00, 1.020332e-01, -4.846724e-03, 4.697008e-04},
+ {1.269541e+00, 1.013216e-01, -4.035237e-03, 3.698061e-04},
+ {1.367197e+00, 1.007075e-01, -3.345636e-03, 2.943847e-04},
+ {1.464853e+00, 1.001716e-01, -2.752379e-03, 2.369418e-04},
+ {1.562509e+00, 9.969996e-02, -2.236687e-03, 1.929401e-04},
+ {1.660166e+00, 9.928152e-02, -1.784441e-03, 1.591361e-04},
+ {1.757822e+00, 9.890786e-02, -1.384813e-03, 1.331654e-04},
+ {1.855478e+00, 9.857228e-02, -1.029346e-03, 1.132745e-04},
+ {1.953134e+00, 9.826940e-02, -7.113359e-04, 9.814323e-05},
+ {2.050790e+00, 9.799484e-02, -4.253906e-04, 8.676449e-05},
+ {2.148447e+00, 9.774497e-02, -1.671167e-04, 7.836135e-05},
+ {2.246103e+00, 9.751678e-02, 6.710498e-05, 7.232835e-05},
+ {2.343759e+00, 9.730773e-02, 2.802876e-04, 6.818952e-05},
+ {2.441415e+00, 9.711568e-02, 4.749628e-04, 6.556786e-05},
+ {2.539071e+00, 9.693877e-02, 6.532755e-04, 6.416282e-05},
+ {2.636728e+00, 9.677542e-02, 8.170569e-04, 6.373341e-05},
+ {2.734384e+00, 1.910291e-01, 3.767224e-03, 5.160521e-04},
+ {2.929696e+00, 1.899901e-01, 4.782956e-03, 5.393740e-04},
+ {3.125009e+00, 1.890915e-01, 5.646906e-03, 5.740339e-04},
+ {3.320321e+00, 1.883086e-01, 6.388338e-03, 6.155239e-04},
+ {3.515634e+00, 1.876218e-01, 7.029801e-03, 6.608298e-04},
+ {3.710946e+00, 1.870156e-01, 7.588922e-03, 7.079357e-04},
+ {3.906258e+00, 1.864773e-01, 8.079643e-03, 7.555006e-04},
+ {4.101571e+00, 1.859967e-01, 8.513103e-03, 8.026422e-04},
+ {4.296883e+00, 1.855653e-01, 8.898285e-03, 8.487922e-04},
+ {4.492196e+00, 1.851763e-01, 9.242495e-03, 8.935968e-04},
+ {4.687508e+00, 1.848239e-01, 9.551712e-03, 9.368494e-04},
+ {4.882821e+00, 1.845031e-01, 9.830867e-03, 9.784433e-04},
+ {5.078133e+00, 1.842100e-01, 1.008405e-02, 1.018340e-03},
+ {5.273445e+00, 1.839412e-01, 1.031468e-02, 1.056544e-03},
+ {5.468758e+00, 1.836937e-01, 1.052561e-02, 1.093093e-03},
+ {5.664070e+00, 1.834651e-01, 1.071926e-02, 1.128039e-03},
+ {5.859383e+00, 1.832533e-01, 1.089768e-02, 1.161448e-03},
+ {6.054695e+00, 1.830564e-01, 1.106262e-02, 1.193391e-03},
+ {6.250007e+00, 1.828729e-01, 1.121556e-02, 1.223940e-03},
+ {6.445320e+00, 1.827015e-01, 1.135780e-02, 1.253169e-03},
+ {6.640632e+00, 1.825408e-01, 1.149043e-02, 1.281149e-03},
+ {6.835945e+00, 1.823900e-01, 1.161443e-02, 1.307949e-03},
+ {7.031257e+00, 1.822482e-01, 1.173063e-02, 1.333635e-03},
+ {7.226570e+00, 1.821144e-01, 1.183975e-02, 1.358270e-03},
+ {7.421882e+00, 1.819881e-01, 1.194245e-02, 1.381911e-03},
+ {7.617194e+00, 1.818685e-01, 1.203929e-02, 1.404614e-03},
+ {7.812507e+00, 1.817552e-01, 1.213077e-02, 1.426431e-03},
+ {8.007819e+00, 1.816477e-01, 1.221734e-02, 1.447411e-03},
+ {8.203132e+00, 1.815454e-01, 1.229938e-02, 1.467598e-03},
+ {8.398444e+00, 1.814481e-01, 1.237726e-02, 1.487035e-03},
+ {8.593757e+00, 1.813554e-01, 1.245129e-02, 1.505762e-03},
+ {8.789069e+00, 1.812668e-01, 1.252176e-02, 1.523815e-03},
+ {8.984381e+00, 3.367030e-01, 4.147623e-02, 1.244560e-02},
+ {9.375006e+00, 3.361379e-01, 4.177920e-02, 1.270774e-02},
+ {9.765631e+00, 3.356191e-01, 4.205304e-02, 1.295269e-02},
+ {1.015626e+01, 3.351411e-01, 4.230171e-02, 1.318205e-02},
+ {1.054688e+01, 3.346991e-01, 4.252850e-02, 1.339724e-02},
+ {1.093751e+01, 3.342892e-01, 4.273612e-02, 1.359950e-02},
+ {1.132813e+01, 3.339080e-01, 4.292687e-02, 1.378995e-02},
+ {1.171876e+01, 3.335526e-01, 4.310271e-02, 1.396956e-02},
+ {1.210938e+01, 3.332203e-01, 4.326528e-02, 1.413922e-02},
+ {1.250000e+01, 5.722553e-01, 9.204483e-02, 1.169495e-01},
+ {1.328125e+01, 5.704529e-01, 9.145775e-02, 1.193390e-01},
+ {1.406250e+01, 5.688543e-01, 9.089620e-02, 1.214992e-01},
+ {1.484375e+01, 5.674266e-01, 9.036152e-02, 1.234615e-01},
+ {1.562500e+01, 8.513274e-01, -3.717295e-01, 1.082901e+00},
+ {1.718750e+01, 8.456524e-01, -3.938066e-01, 1.110653e+00},
+ {1.875000e+01, 8.409733e-01, -4.130061e-01, 1.134532e+00},
+ {2.031250e+01, 8.371720e-01, -4.305845e-01, 1.155912e+00},
+ {2.187500e+01, 8.346323e-01, -4.503617e-01, 1.178229e+00},
+ {2.343750e+01, 8.359434e-01, -4.895204e-01, 1.216076e+00}};
+static int index_table_a[165] = {
+ 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2,
+ 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4,
+ 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6,
+ 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
+ 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 9,
+ 10, 10, 10, 10, 11, 11, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 13, 13, 13,
+ 13, 14, 14, 14, 15, 15, 15, 15, 16, 16, 16, 16, 17, 17, 17, 18, 18, 18, 19,
+ 19, 20, 20, 21, 21, 22, 23, 24, 24, 24, 25, 25, 26, 26, 27, 27, 28, 29, 29,
+ 30, 31, 32, 33, 34, 35, 36, 36, 37, 38, 38, 39, 40};
+static int index_table_b[165] = {
+ 41, 42, 43, 44, 45, 46, 47, 48, 48, 49, 49, 50, 50, 50, 50, 51, 51, 51, 52,
+ 52, 52, 52, 53, 53, 53, 53, 54, 54, 54, 54, 55, 55, 55, 55, 56, 56, 56, 56,
+ 57, 57, 57, 57, 57, 58, 58, 58, 58, 59, 59, 59, 59, 59, 60, 60, 60, 60, 61,
+ 61, 61, 61, 61, 62, 62, 62, 62, 63, 63, 63, 63, 63, 64, 64, 64, 64, 65, 65,
+ 65, 65, 66, 66, 66, 66, 66, 66, 66, 66, 66, 67, 67, 67, 67, 67, 67, 67, 67,
+ 68, 68, 68, 68, 68, 68, 68, 69, 69, 69, 69, 69, 69, 69, 70, 70, 70, 70, 70,
+ 70, 70, 71, 71, 71, 71, 71, 71, 72, 72, 72, 72, 72, 73, 73, 73, 73, 73, 74,
+ 74, 74, 74, 75, 75, 75, 75, 76, 76, 76, 76, 77, 77, 77, 78, 78, 78, 79, 79,
+ 80, 80, 80, 81, 81, 82, 83, 83, 84, 85, 86, 87, 88};
+static int index_table_c[165] = {
+ 89, 90, 90, 91, 92, 92, 93, 93, 94, 94, 95, 96, 96, 97, 97,
+ 98, 98, 98, 99, 99, 99, 99, 100, 100, 100, 101, 101, 101, 101, 102,
+ 102, 102, 102, 103, 103, 103, 104, 104, 104, 104, 105, 105, 105, 105, 106,
+ 106, 106, 107, 107, 107, 107, 107, 107, 107, 107, 108, 108, 108, 108, 108,
+ 108, 108, 109, 109, 109, 109, 109, 109, 109, 109, 110, 110, 110, 110, 110,
+ 110, 110, 110, 111, 111, 111, 111, 111, 111, 111, 111, 111, 111, 111, 111,
+ 111, 111, 111, 112, 112, 112, 112, 112, 112, 112, 112, 112, 112, 112, 112,
+ 112, 112, 112, 112, 113, 113, 113, 113, 113, 113, 113, 113, 113, 113, 113,
+ 113, 113, 113, 113, 113, 114, 114, 114, 114, 114, 114, 114, 114, 114, 114,
+ 114, 114, 114, 114, 114, 114, 115, 115, 115, 115, 115, 115, 115, 115, 115,
+ 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115};
+static struct anyrng anyrng = {118, endpoints, splines, 165,
+ index_table_a, index_table_b, index_table_c};
+
+/**
+ * @brief Transform a 64-bit unsigned integer seed into a (dimensionless)
+ * Fermi-Dirac momentum (units of kb*T), using cubic spline interpolation of
+ * the quantile function.
+ *
+ * @param seed Random seed to be transformed
+ */
+double neutrino_seed_to_fermi_dirac(uint64_t seed) {
+ /* Scramble the bits with splitmix64 */
+ uint64_t A = seed;
+ A = A + 0x9E3779B97f4A7C15;
+ A = (A ^ (A >> 30)) * 0xBF58476D1CE4E5B9;
+ A = (A ^ (A >> 27)) * 0x94D049BB133111EB;
+ A = A ^ (A >> 31);
+
+ /* Map the integer to the unit open interval (0, 1) */
+ const double u = ((double)A + 0.5) / ((double)UINT64_MAX + 1);
+
+ /* Use the hash table to find an enclosing interval */
+ const int tablen = anyrng.tablelen;
+ int index, interval;
+ if (u < 0.025) {
+ index = (int)((optimized_log10f(u) + 14.5229) / 12.9208 * tablen);
+ interval = anyrng.index_table_a[index < tablen ? index : tablen - 1] + 1;
+ } else if (u < 0.975) {
+ index = (int)((u - 0.025) / 0.95 * tablen);
+ interval = anyrng.index_table_b[index < tablen ? index : tablen - 1] + 1;
+ } else {
+ index = (int)(-(optimized_log10f(1 - u) + 1.60206) / 6.39794 * tablen);
+ interval = anyrng.index_table_c[index < tablen ? index : tablen - 1] + 1;
+ }
+
+ /* Retrieve the endpoints and cubic spline coefficients of this interval */
+ const double Fl = anyrng.endpoints[interval];
+ const double Fr = anyrng.endpoints[interval + 1];
+ struct spline *iv = &anyrng.splines[interval];
+
+ /* Evaluate F^-1(u) using the Hermite approximation of F in this interval */
+ const double u_tilde = (u - Fl) / (Fr - Fl);
+ const double u_tilde2 = u_tilde * u_tilde;
+ const double u_tilde3 = u_tilde2 * u_tilde;
+ return iv->a0 + iv->a1 * u_tilde + iv->a2 * u_tilde2 + iv->a3 * u_tilde3;
+}
+
+/**
+ * @brief Transform a 64-bit unsigned integer seed into a point on the sphere.
+ *
+ * @param seed Random seed to be transformed
+ */
+void neutrino_seed_to_direction(uint64_t seed, double n[3]) {
+ /* Skip the first draw, which is used for the magnitude */
+ seed += 0x9E3779B97f4A7C15;
+
+ /* Draw three numbers for the orientation of the momentum vector */
+ double u[3];
+ for (int i = 0; i < 3; i++) {
+ /* Scramble the bits with splitmix64 */
+ uint64_t N = seed += 0x9E3779B97f4A7C15;
+ N = (N ^ (N >> 30)) * 0xBF58476D1CE4E5B9;
+ N = (N ^ (N >> 27)) * 0x94D049BB133111EB;
+ N = N ^ (N >> 31);
+ /* Map the integer to the unit open interval (0,1) */
+ u[i] = ((double)N + 0.5) / ((double)UINT64_MAX + 1);
+ }
+
+ /* Map to three independent Gaussians with Box-Mueller */
+ const double sqrt_log_u = sqrt(-2 * log(u[0]));
+ n[0] = sqrt_log_u * cos(2 * M_PI * u[1]);
+ n[1] = sqrt_log_u * sin(2 * M_PI * u[1]);
+ n[2] = sqrt_log_u * cos(2 * M_PI * u[2]);
+
+ /* Normalize the vector */
+ const double nmag = sqrtf(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]);
+ if (nmag > 0) {
+ const double inv_nmag = 1.0 / nmag;
+ n[0] *= inv_nmag;
+ n[1] *= inv_nmag;
+ n[2] *= inv_nmag;
+ }
+}
diff --git a/src/neutrino/Default/fermi_dirac.h b/src/neutrino/Default/fermi_dirac.h
new file mode 100644
index 0000000000000000000000000000000000000000..8c47a12568e6b5ea4e9913e9ebbf38ace299a485
--- /dev/null
+++ b/src/neutrino/Default/fermi_dirac.h
@@ -0,0 +1,58 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2021 Willem Elbers (whe@willemelbers.com)
+ *
+ * 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_DEFAULT_FERMI_DIRAC_H
+#define SWIFT_DEFAULT_FERMI_DIRAC_H
+
+/* Standard headers */
+#include <math.h>
+#include <stdint.h>
+
+/* Faster exponential */
+#include "exp.h"
+
+/**
+ * @brief Calculate the neutrino density at the particle's location in phase
+ * space, according to the 0th order background model: f_0(x,p,t).
+ *
+ * @param z Argument of the FD function: z = p/kbT
+ */
+INLINE static float fermi_dirac_density(float z) {
+ return 1.f / (optimized_expf(z) + 1.f);
+}
+
+/**
+ * @brief Return a microscopic neutrino mass in eV based on the particle seed.
+ * We simply cycle through the masses defined in the cosmology. TODO: consider
+ * multiplicities.
+ *
+ * @param N_nu Number of distinct neutrino masses
+ * @param m_eV_array Array of masses in eV
+ * @param seed The seed of the neutrino particle
+ */
+INLINE static double neutrino_seed_to_mass(const int N_nu,
+ const double *m_eV_array,
+ uint64_t seed) {
+ return m_eV_array[(int)(seed % N_nu)];
+}
+
+double neutrino_seed_to_fermi_dirac(uint64_t seed);
+void neutrino_seed_to_direction(uint64_t seed, double n[3]);
+
+#endif /* SWIFT_DEFAULT_FERMI_DIRAC_H */
diff --git a/src/neutrino/neutrino.h b/src/neutrino/Default/neutrino.h
similarity index 58%
rename from src/neutrino/neutrino.h
rename to src/neutrino/Default/neutrino.h
index 13e7bd1d1a05ab2cf88a46b679a63a0c8c28bde4..a84688d42b0115f2c8e0bf671ef63cd4561da2ff 100644
--- a/src/neutrino/neutrino.h
+++ b/src/neutrino/Default/neutrino.h
@@ -19,6 +19,15 @@
#ifndef SWIFT_DEFAULT_NEUTRINO_H
#define SWIFT_DEFAULT_NEUTRINO_H
+/* Config parameters. */
+#include "../config.h"
+
+/* Local includes */
+#include "../../engine.h"
+#include "fermi_dirac.h"
+#include "neutrino_properties.h"
+#include "relativity.h"
+
/* Riemann function zeta(3) */
#define M_ZETA_3 1.2020569031595942853997
@@ -63,4 +72,53 @@ INLINE static double neutrino_mass_factor(
return mass_factor_eV;
}
+/**
+ * @brief Initialises the neutrino g-particles for the first time. This is done
+ * in addition to gravity_first_init_gpart().
+ *
+ * This function is called only once just after the ICs have been read in
+ * and after IDs have been remapped (if used) by space_remap_ids().
+ *
+ * @param gp The particle to act upon
+ * @param engine The engine of the run
+ */
+__attribute__((always_inline)) INLINE static void gravity_first_init_neutrino(
+ struct gpart *gp, const struct engine *e) {
+
+ /* Do we need to do anything? */
+ if (!e->neutrino_properties->generate_ics) return;
+
+ /* Retrieve physical and cosmological constants */
+ const double c_vel = e->physical_constants->const_speed_light_c;
+ const double *m_eV_array = e->cosmology->M_nu_eV;
+ const int N_nu = e->cosmology->N_nu;
+ const double T_eV = e->cosmology->T_nu_0_eV;
+ const double inv_fac = c_vel * T_eV;
+ const long long neutrino_seed = e->neutrino_properties->neutrino_seed;
+
+ /* Use a particle id dependent seed */
+ const long long seed = gp->id_or_neg_offset + neutrino_seed;
+
+ /* Compute the initial dimensionless momentum from the seed */
+ const double pi = neutrino_seed_to_fermi_dirac(seed);
+
+ /* The neutrino mass (we cycle based on the neutrino seed) */
+ const double m_eV = neutrino_seed_to_mass(N_nu, m_eV_array, seed);
+
+ /* Compute the initial direction of the momentum vector from the seed */
+ double n[3];
+ neutrino_seed_to_direction(seed, n);
+
+ /* Set the initial velocity */
+ const double vi = pi * inv_fac / m_eV;
+ gp->v_full[0] = n[0] * vi;
+ gp->v_full[1] = n[1] * vi;
+ gp->v_full[2] = n[2] * vi;
+
+ /* If running with the delta-f method, set the weight to (almost) zero */
+ if (e->neutrino_properties->use_delta_f) {
+ gp->mass = FLT_MIN;
+ }
+}
+
#endif /* SWIFT_DEFAULT_NEUTRINO_H */
diff --git a/src/neutrino/Default/neutrino_io.h b/src/neutrino/Default/neutrino_io.h
new file mode 100644
index 0000000000000000000000000000000000000000..b1ee5a7490bbde79b2db06446c9d70b09e9d320f
--- /dev/null
+++ b/src/neutrino/Default/neutrino_io.h
@@ -0,0 +1,129 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2021 Willem Elbers (willem.h.elbers@durham.ac.uk)
+ *
+ * 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_DEFAULT_NEUTRINO_IO_H
+#define SWIFT_DEFAULT_NEUTRINO_IO_H
+
+#include "../../gravity.h"
+
+/* Local includes */
+#include "fermi_dirac.h"
+#include "neutrino.h"
+#include "neutrino_properties.h"
+
+/**
+ * @brief Recover and store the initial Fermi-Dirac speed, vi, for a neutrino
+ * particle. These are used in the weight calculation of the delta-f method.
+ *
+ * @param e The engine of the run
+ * @param gp The neutrino gpart in question
+ * @param ret Output
+ */
+INLINE static void convert_gpart_vi(const struct engine* e,
+ const struct gpart* gp, float* ret) {
+
+ /* When we are running with the delta-f method, resample the momentum */
+ if (e->neutrino_properties->use_delta_f) {
+ /* Retrieve physical constants, including the neutrino mass array */
+ const double a_scale = e->cosmology->a;
+ const double N_nu = e->cosmology->N_nu;
+ const double* m_eV_array = e->cosmology->M_nu_eV;
+ const double c_vel = e->physical_constants->const_speed_light_c;
+ const double T_eV = e->cosmology->T_nu_0_eV;
+ const double a_fac = c_vel * T_eV / a_scale;
+
+ /* Use a particle id dependent seed (sum of global seed and ID) */
+ const long long neutrino_seed = e->neutrino_properties->neutrino_seed;
+ const long long seed = gp->id_or_neg_offset + neutrino_seed;
+
+ /* Convert momentum in electronvolts to speed in internal units */
+ double pi_eV = neutrino_seed_to_fermi_dirac(seed);
+ double m_eV = neutrino_seed_to_mass(N_nu, m_eV_array, seed);
+ double vi = pi_eV / m_eV * a_fac; // scales like peculiar velocity
+
+ ret[0] = vi;
+ } else {
+ /* We don't know what the initial momentum was and we don't need it */
+ ret[0] = 0.f;
+ }
+}
+
+/**
+ * @brief Recover and store the microscopic mass of a neutrino particle
+ *
+ * @param e The engine of the run
+ * @param gp The neutrino gpart in question
+ * @param ret Output
+ */
+INLINE static void convert_gpart_mnu(const struct engine* e,
+ const struct gpart* gp, double* ret) {
+
+ /* Physical constants */
+ const double c = e->physical_constants->const_speed_light_c;
+ const double eV = e->physical_constants->const_electron_volt;
+ const double eV_mass = eV / (c * c);
+
+ double micro_mass;
+
+ /* When we are running with the delta-f method, resample the mass */
+ if (e->neutrino_properties->use_delta_f) {
+
+ /* Use a particle id dependent seed (sum of global seed and ID) */
+ const long long neutrino_seed = e->neutrino_properties->neutrino_seed;
+ const long long seed = gp->id_or_neg_offset + neutrino_seed;
+
+ /* Fetch neutrino masses defined in the cosmology */
+ const int N_nu = e->cosmology->N_nu;
+ const double* m_eV_array = e->cosmology->M_nu_eV;
+
+ micro_mass = neutrino_seed_to_mass(N_nu, m_eV_array, seed); // eV
+ } else {
+ /* Otherwise, simply use the mass implied by the conversion factor */
+ const double mass_factor = e->neutrino_mass_conversion_factor;
+
+ micro_mass = gp->mass * mass_factor; // eV
+ }
+
+ /* Convert units and store the answer */
+ ret[0] = micro_mass * eV_mass;
+}
+
+/**
+ * @brief Specifies which particle fields to write to a dataset
+ *
+ * @param gparts The particle array.
+ * @param list The list of i/o properties to write.
+ *
+ * @return Returns the number of fields to write.
+ */
+__attribute__((always_inline)) INLINE static int neutrino_write_particles(
+ const struct gpart* gparts, struct io_props* list) {
+
+ list[0] = io_make_output_field_convert_gpart(
+ "SampledSpeeds", FLOAT, 1, UNIT_CONV_SPEED, 0.f, gparts, convert_gpart_vi,
+ "Initial Fermi-Dirac speed sampled at infinity. This is a * |dx/dt| "
+ "where x is the co-moving position of the particles.");
+
+ list[1] = io_make_output_field_convert_gpart(
+ "MicroscopicMasses", DOUBLE, 1, UNIT_CONV_MASS, 0.f, gparts,
+ convert_gpart_mnu, "Microscopic masses of individual neutrino particles");
+
+ return 2;
+}
+
+#endif /* SWIFT_DEFAULT_NEUTRINO_IO_H */
diff --git a/src/neutrino/Default/neutrino_properties.h b/src/neutrino/Default/neutrino_properties.h
new file mode 100644
index 0000000000000000000000000000000000000000..8dbc4cb568086fea6e8dc6210ff4a276e2fb4416
--- /dev/null
+++ b/src/neutrino/Default/neutrino_properties.h
@@ -0,0 +1,87 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2021 Willem Elbers (willem.h.elbers@durham.ac.uk)
+ *
+ * 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_DEFAULT_NEUTRINO_PROPERTIES_H
+#define SWIFT_DEFAULT_NEUTRINO_PROPERTIES_H
+
+#include "../../restart.h"
+
+/**
+ * @brief Properties of neutrinos
+ */
+struct neutrino_props {
+
+ /* Whether to run with the delta-f method for neutrino weighting */
+ char use_delta_f;
+
+ /* Whether to generate random neutrino velocities in the initial conditions */
+ char generate_ics;
+
+ /* Random seed for the neutrino weighting task */
+ long long neutrino_seed;
+};
+
+/**
+ * @brief Initialise the neutrino properties from the parameter file.
+ *
+ * @param np The #neutrino_props.
+ * @param phys_const The physical constants in the internal unit system.
+ * @param us The internal unit system.
+ * @param params The parsed parameters.
+ * @param cosmo The cosmological model.
+ */
+INLINE static void neutrino_props_init(struct neutrino_props *np,
+ const struct phys_const *phys_const,
+ const struct unit_system *us,
+ struct swift_params *params,
+ const struct cosmology *cosmo) {
+
+ np->use_delta_f = parser_get_param_int(params, "Neutrino:use_delta_f");
+ np->generate_ics = parser_get_param_int(params, "Neutrino:generate_ics");
+ np->neutrino_seed =
+ parser_get_opt_param_longlong(params, "Neutrino:neutrino_seed", 0);
+}
+
+/**
+ * @brief Write a neutrino_props struct to the given FILE as a stream of
+ * bytes.
+ *
+ * @param props the neutrino properties struct
+ * @param stream the file stream
+ */
+INLINE static void neutrino_struct_dump(const struct neutrino_props *props,
+ FILE *stream) {
+
+ restart_write_blocks((void *)props, sizeof(struct neutrino_props), 1, stream,
+ "neutrino props", "Neutrino props");
+}
+
+/**
+ * @brief Restore a neutrino_props struct from the given FILE as a stream of
+ * bytes.
+ *
+ * @param props the neutrino properties struct
+ * @param stream the file stream
+ */
+INLINE static void neutrino_struct_restore(const struct neutrino_props *props,
+ FILE *stream) {
+ restart_read_blocks((void *)props, sizeof(struct neutrino_props), 1, stream,
+ NULL, "Neutrino props");
+}
+
+#endif /* SWIFT_DEFAULT_NEUTRINO_PROPERTIES_H */
diff --git a/src/neutrino/relativity.h b/src/neutrino/Default/relativity.h
similarity index 100%
rename from src/neutrino/relativity.h
rename to src/neutrino/Default/relativity.h
diff --git a/src/neutrino_io.h b/src/neutrino_io.h
new file mode 100644
index 0000000000000000000000000000000000000000..f23d169dc5449ce8423f80f8a037ac0a4a2b1fc4
--- /dev/null
+++ b/src/neutrino_io.h
@@ -0,0 +1,28 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2021 Willem Elbers (whe@willemelbers.com)
+ *
+ * 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_NEUTRINO_IO_H
+#define SWIFT_NEUTRINO_IO_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Select the correct neutrino model */
+#include "./neutrino/Default/neutrino_io.h"
+
+#endif /* SWIFT_NEUTRINO_PROPERTIES_H */
diff --git a/src/neutrino_properties.h b/src/neutrino_properties.h
new file mode 100644
index 0000000000000000000000000000000000000000..756d4e0c6c1a3bebb3e1849c950b5ab20aca76a1
--- /dev/null
+++ b/src/neutrino_properties.h
@@ -0,0 +1,28 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2021 Willem Elbers (whe@willemelbers.com)
+ *
+ * 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_NEUTRINO_PROPERTIES_H
+#define SWIFT_NEUTRINO_PROPERTIES_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Select the correct neutrino model */
+#include "./neutrino/Default/neutrino_properties.h"
+
+#endif /* SWIFT_NEUTRINO_PROPERTIES_H */
diff --git a/src/parallel_io.c b/src/parallel_io.c
index 07a4926f021d797edb2a8c4bf98dfa98c59bc328..d9d8a5c7cc0aa2f7383361d78d39704783552d7c 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -1299,11 +1299,15 @@ void prepare_file(struct engine* e, const char* fileName,
case swift_type_dark_matter:
case swift_type_dark_matter_background:
- case swift_type_neutrino:
io_select_dm_fields(NULL, NULL, with_fof, with_stf, e, &num_fields,
list);
break;
+ case swift_type_neutrino:
+ io_select_neutrino_fields(NULL, NULL, with_fof, with_stf, e,
+ &num_fields, list);
+ break;
+
case swift_type_sink:
io_select_sink_fields(NULL, with_cosmology, with_fof, with_stf, e,
&num_fields, list);
@@ -1773,8 +1777,8 @@ void write_output_parallel(struct engine* e,
gpart_group_data_written, Ntot, Ndm_neutrino, with_stf);
/* Select the fields to write */
- io_select_dm_fields(gparts_written, gpart_group_data_written, with_fof,
- with_stf, e, &num_fields, list);
+ io_select_neutrino_fields(gparts_written, gpart_group_data_written,
+ with_fof, with_stf, e, &num_fields, list);
} break;
diff --git a/src/runner.h b/src/runner.h
index 008ea75e6790fbed32567378fa020620706b6f38..8c71570d9c2e131e3cea2706b7941bab5c2207e4 100644
--- a/src/runner.h
+++ b/src/runner.h
@@ -145,6 +145,7 @@ void runner_do_recv_spart(struct runner *r, struct cell *c, int clear_sorts,
int timer);
void runner_do_recv_bpart(struct runner *r, struct cell *c, int clear_sorts,
int timer);
+void runner_do_neutrino_weighting(struct runner *r, struct cell *c, int timer);
void *runner_main(void *data);
#endif /* SWIFT_RUNNER_H */
diff --git a/src/runner_main.c b/src/runner_main.c
index 2cef97f9b25032834bad85c3d64de4e075599563..988e42b69f29f5b09236d8d32f2b9854b524b3a0 100644
--- a/src/runner_main.c
+++ b/src/runner_main.c
@@ -610,6 +610,9 @@ void *runner_main(void *data) {
case task_type_fof_pair:
runner_do_fof_pair(r, t->ci, t->cj, 1);
break;
+ case task_type_neutrino_weight:
+ runner_do_neutrino_weighting(r, ci, 1);
+ break;
case task_type_rt_ghost1:
runner_do_rt_ghost1(r, t->ci, 1);
break;
diff --git a/src/runner_neutrino.c b/src/runner_neutrino.c
new file mode 100644
index 0000000000000000000000000000000000000000..ba380e71875680e87059149b9ddbc90279deb235
--- /dev/null
+++ b/src/runner_neutrino.c
@@ -0,0 +1,125 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2020 Willem Elbers (whe@willemelbers.com)
+ *
+ * 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/>.
+ *
+ ******************************************************************************/
+
+/* Config parameters. */
+#include "../config.h"
+
+/* This object's header. */
+#include "runner.h"
+
+/* Phase space density functions needed */
+#include "neutrino.h"
+#include "neutrino_properties.h"
+
+/* Local headers. */
+#include "active.h"
+#include "cell.h"
+#include "engine.h"
+#include "timers.h"
+
+/* Compute the dimensionless neutrino momentum (units of kb*T).
+ *
+ * @param v The internal 3-velocity
+ * @param m_eV The neutrino mass in electron-volts
+ * @param fac Conversion factor = 1. / (speed_of_light * T_nu_eV)
+ */
+INLINE static double neutrino_momentum(float v[3], double m_eV, double fac) {
+
+ float v2 = v[0] * v[0] + v[1] * v[1] + v[2] * v[2];
+ float vmag = sqrtf(v2);
+ double p = vmag * fac * m_eV;
+ return p;
+}
+
+/**
+ * @brief Weight the active neutrino particles in a cell using the delta-f
+ * method.
+ *
+ * @param r The runner thread.
+ * @param c The cell.
+ * @param timer Are we timing this ?
+ */
+void runner_do_neutrino_weighting(struct runner *r, struct cell *c, int timer) {
+
+ const struct engine *e = r->e;
+ const int with_cosmology = (e->policy & engine_policy_cosmology);
+ struct gpart *restrict gparts = c->grav.parts;
+ const int gcount = c->grav.count;
+
+ TIMER_TIC;
+
+ /* Anything to do here? */
+ if (c->grav.count == 0) return;
+ if (!cell_is_starting_gravity(c, e) && !cell_is_active_gravity(c, e)) return;
+ if (!with_cosmology)
+ error("Phase space weighting without cosmology not implemented.");
+
+ /* Retrieve physical and cosmological constants */
+ const double c_vel = e->physical_constants->const_speed_light_c;
+ const double *m_eV_array = e->cosmology->M_nu_eV;
+ const int N_nu = e->cosmology->N_nu;
+ const double T_eV = e->cosmology->T_nu_0_eV;
+ const double fac = 1.0 / (c_vel * T_eV);
+ const double inv_mass_factor = 1. / e->neutrino_mass_conversion_factor;
+ const long long neutrino_seed = e->neutrino_properties->neutrino_seed;
+
+ /* Recurse? */
+ if (c->split) {
+ for (int k = 0; k < 8; k++)
+ if (c->progeny[k] != NULL)
+ runner_do_neutrino_weighting(r, c->progeny[k], 0);
+ } else {
+ /* Loop over the gparts in this cell. */
+ for (int k = 0; k < gcount; k++) {
+ /* Get a handle on the part. */
+ struct gpart *restrict gp = &gparts[k];
+
+ /* Only act on neutrinos */
+ if (gp->type != swift_type_neutrino) continue;
+
+ /* Use a particle id dependent seed */
+ const long long seed = gp->id_or_neg_offset + neutrino_seed;
+
+ /* Compute the initial dimensionless momentum from the seed */
+ const double pi = neutrino_seed_to_fermi_dirac(seed);
+
+ /* The neutrino mass (we cycle based on the neutrino seed) */
+ const double m_eV = neutrino_seed_to_mass(N_nu, m_eV_array, seed);
+ const double mass = m_eV * inv_mass_factor;
+
+ /* Compute the current dimensionless momentum */
+ double p = neutrino_momentum(gp->v_full, m_eV, fac);
+
+ /* Compute the initial and current background phase-space density */
+ double fi = fermi_dirac_density(pi);
+ double f = fermi_dirac_density(p);
+ double weight = 1.0 - f / fi;
+
+ /* Set the statistically weighted mass */
+ gp->mass = mass * weight;
+
+ /* Prevent degeneracies */
+ if (gp->mass == 0.) {
+ gp->mass = FLT_MIN;
+ }
+ }
+ }
+
+ if (timer) TIMER_TOC(timer_weight);
+}
diff --git a/src/serial_io.c b/src/serial_io.c
index 015d0e41558b25a5e1079fc2e01acaba9aa0d344..0e792f2a012a60303035e1190e7d01ff8467f47b 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -1410,8 +1410,8 @@ void write_output_serial(struct engine* e,
gpart_group_data_written, Ntot, Ndm_neutrino, with_stf);
/* Select the fields to write */
- io_select_dm_fields(gparts_written, gpart_group_data_written,
- with_fof, with_stf, e, &num_fields, list);
+ io_select_neutrino_fields(gparts_written, gpart_group_data_written,
+ with_fof, with_stf, e, &num_fields, list);
} break;
diff --git a/src/single_io.c b/src/single_io.c
index 5deae3e179908ae239ca19a408f39710add8fd54..39377033a1ff24ec67fb5ac3180a241a5cff8e03 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -1189,8 +1189,8 @@ void write_output_single(struct engine* e,
gpart_group_data_written, Ntot, Ndm_neutrino, with_stf);
/* Select the fields to write */
- io_select_dm_fields(gparts_written, gpart_group_data_written, with_fof,
- with_stf, e, &num_fields, list);
+ io_select_neutrino_fields(gparts_written, gpart_group_data_written,
+ with_fof, with_stf, e, &num_fields, list);
} break;
diff --git a/src/space_first_init.c b/src/space_first_init.c
index 915b3646244a31b838588beb4f155ebfee5b3792..2af63445d47217acdf4b744dc6c23a77150b9351 100644
--- a/src/space_first_init.c
+++ b/src/space_first_init.c
@@ -30,6 +30,7 @@
#include "chemistry.h"
#include "engine.h"
#include "gravity.h"
+#include "neutrino.h"
#include "particle_splitting.h"
#include "pressure_floor.h"
#include "rt.h"
@@ -206,6 +207,9 @@ void space_first_init_gparts_mapper(void *restrict map_data, int count,
gravity_first_init_gpart(&gp[k], grav_props);
+ if (gp[k].type == swift_type_neutrino)
+ gravity_first_init_neutrino(&gp[k], s->e);
+
#ifdef WITH_LOGGER
logger_part_data_init(&gp[k].logger_data);
#endif
diff --git a/src/swift.h b/src/swift.h
index 4018c9d22a9e6b9e3de971b9cf54ed3391dadd20..06741a5a96fba8191643b9cd39bc620934bc1821 100644
--- a/src/swift.h
+++ b/src/swift.h
@@ -60,6 +60,7 @@
#include "minmax.h"
#include "mpiuse.h"
#include "multipole.h"
+#include "neutrino_properties.h"
#include "output_list.h"
#include "output_options.h"
#include "parallel_io.h"
diff --git a/src/task.c b/src/task.c
index edf398ff2d9556f910539da7c6a477f333ab612d..44854e5ff2ad47e05839b36a08ed25975f095cc3 100644
--- a/src/task.c
+++ b/src/task.c
@@ -107,6 +107,7 @@ const char *taskID_names[task_type_count] = {
"bh_swallow_ghost3",
"fof_self",
"fof_pair",
+ "neutrino_weight",
"sink_in",
"sink_ghost",
"sink_out",
@@ -169,7 +170,7 @@ const char *task_category_names[task_category_count] = {
"black holes", "cooling", "star formation",
"limiter", "sync", "time integration",
"mpi", "fof", "others",
- "sink", "RT"};
+ "neutrino", "sink", "RT"};
#ifdef WITH_MPI
/* MPI communicators for the subtypes. */
@@ -1764,6 +1765,9 @@ enum task_categories task_get_category(const struct task *t) {
case task_type_rt_out:
return task_category_rt;
+ case task_type_neutrino_weight:
+ return task_category_neutrino;
+
case task_type_self:
case task_type_pair:
case task_type_sub_self:
diff --git a/src/task.h b/src/task.h
index 5c6f15f67c7a154128575fcb16a4717bd0234557..96a218a87d34f928c2e58eb71ad85e6cf6ea4b10 100644
--- a/src/task.h
+++ b/src/task.h
@@ -100,6 +100,7 @@ enum task_types {
task_type_bh_swallow_ghost3, /* Implicit */
task_type_fof_self,
task_type_fof_pair,
+ task_type_neutrino_weight,
task_type_sink_in, /* Implicit */
task_type_sink_ghost, /* Implicit */
task_type_sink_out, /* Implicit */
@@ -194,6 +195,7 @@ enum task_categories {
task_category_mpi,
task_category_fof,
task_category_others,
+ task_category_neutrino,
task_category_sink,
task_category_rt,
task_category_count
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 559ea106868673dcc3b6d51bb8e0a1a829407c57..8f72493890ffd7d0eda54dc9ebc0144b1aa23526 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -30,7 +30,8 @@ TESTS = testGreetings testMaths testReading.sh testKernel testKernelLongGrav \
testPotentialPair testEOS testUtilities testSelectOutput.sh \
testCbrt testCosmology testOutputList testFormat.sh \
test27cellsStars.sh test27cellsStarsPerturbed.sh testHydroMPIrules \
- testAtomic testGravitySpeed testNeutrinoCosmology
+ testAtomic testGravitySpeed testNeutrinoCosmology testNeutrinoFermiDirac \
+ testLog
# List of test programs to compile
check_PROGRAMS = testGreetings testReading testTimeIntegration testKernelLongGrav \
@@ -43,7 +44,8 @@ check_PROGRAMS = testGreetings testReading testTimeIntegration testKernelLongGra
testGravityDerivatives testPotentialSelf testPotentialPair testEOS testUtilities \
testSelectOutput testCbrt testCosmology testOutputList test27cellsStars \
test27cellsStars_subset testCooling testComovingCooling testFeedback testHashmap \
- testAtomic testHydroMPIrules testGravitySpeed testNeutrinoCosmology
+ testAtomic testHydroMPIrules testGravitySpeed testNeutrinoCosmology \
+ testNeutrinoFermiDirac testLog
# Rebuild tests when SWIFT is updated.
$(check_PROGRAMS): ../src/.libs/libswiftsim.a
@@ -72,6 +74,10 @@ testSymmetry_SOURCES = testSymmetry.c
# Added because of issues using memcmp on clang 4.x
testSymmetry_CFLAGS = $(AM_CFLAGS) -fno-builtin-memcmp
+testNeutrinoCosmology_SOURCES = testNeutrinoCosmology.c
+
+testNeutrinoFermiDirac_SOURCES = testNeutrinoFermiDirac.c
+
testTimeIntegration_SOURCES = testTimeIntegration.c
testActivePair_SOURCES = testActivePair.c
@@ -148,6 +154,8 @@ testFeedback_SOURCES = testFeedback.c
testHashmap_SOURCES = testHashmap.c
+testLog_SOURCES = testLog.c
+
testHydroMPIrules = testHydroMPIrules.c
# Files necessary for distribution
diff --git a/tests/testLog.c b/tests/testLog.c
new file mode 100644
index 0000000000000000000000000000000000000000..61544e93b76c8383d5c07df6df9ba8a8921c73fb
--- /dev/null
+++ b/tests/testLog.c
@@ -0,0 +1,95 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (C) 2021 Willem Elbers (whe@willemelbers.com)
+ *
+ * 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/>.
+ *
+ ******************************************************************************/
+
+#include "../config.h"
+#include "log.h"
+#include "swift.h"
+
+/* Standard includes */
+#include <fenv.h>
+#include <math.h>
+
+/**
+ * @brief Check that a and b are consistent (up to some relative error)
+ *
+ * @param a First value
+ * @param b Second value
+ * @param tol Relative tolerance
+ * @param x Value for which we tested the function (for error messages)
+ */
+void check_value(double a, double b, const double tol, const double x) {
+ if (fabs(a - b) / fabs(a + b) > tol)
+ error("Values are inconsistent: %12.15e %12.15e rel=%e (for x=%e).", a, b,
+ fabs(a - b) / fabs(a + b), x);
+}
+
+/**
+ * @brief Check that a and b are consistent (up to some absolute error)
+ *
+ * @param a First value
+ * @param b Second value
+ * @param tol Relative tolerance
+ * @param x Value for which we tested the function (for error messages)
+ */
+void check_value_abs(double a, double b, const double tol, const double x) {
+ if (fabs(a - b) > tol)
+ error("Values are inconsistent: %12.15e %12.15e abs=%e (for x=%e).", a, b,
+ fabs(a - b), x);
+}
+
+int main(int argc, char* argv[]) {
+
+ /* Initialize CPU frequency, this also starts time. */
+ unsigned long long cpufreq = 0;
+ clocks_set_cpufreq(cpufreq);
+
+/* Choke on FPEs */
+#ifdef HAVE_FE_ENABLE_EXCEPT
+ feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
+#endif
+
+ /* Loop over some values */
+ for (float x = 0.; x < 32; x += 0.000001) {
+
+ const double exact_p = x;
+ const double exact_n = -x;
+ const double swift_log_p = optimized_logf(exp(x));
+ const double swift_log_n = optimized_logf(exp(-x));
+ const double swift_log10_p = optimized_log10f(exp10(x));
+ const double swift_log10_n = optimized_log10f(exp10(-x));
+
+ /* Check the absolute error */
+ check_value_abs(exact_p, swift_log_p, 5e-3, x);
+ check_value_abs(exact_n, swift_log_n, 5e-3, x);
+ check_value_abs(exact_p, swift_log10_p, 5e-3, x);
+ check_value_abs(exact_n, swift_log10_n, 5e-3, x);
+
+ /* Check the relative error in the domain where it is small */
+ if (x >= 0.3) {
+ check_value(exact_p, swift_log_p, 5e-3, x);
+ check_value(exact_n, swift_log_n, 5e-3, x);
+ check_value(exact_p, swift_log10_p, 5e-3, x);
+ check_value(exact_n, swift_log10_n, 5e-3, x);
+ }
+ }
+
+ message("Success.");
+
+ return 0;
+}
diff --git a/tests/testNeutrinoCosmology.c b/tests/testNeutrinoCosmology.c
index fc2cff575250e21be05e9fbd9acef020445e9ca7..a271c6bc0e50ffe0acd15bf8220533f93d2dc935 100644
--- a/tests/testNeutrinoCosmology.c
+++ b/tests/testNeutrinoCosmology.c
@@ -18,13 +18,11 @@
******************************************************************************/
/* Some standard headers. */
-#include "../config.h"
-
-/* Some standard headers */
#include <fenv.h>
#include <math.h>
/* Includes. */
+#include "../config.h"
#include "swift.h"
#define N_CHECK 20
diff --git a/tests/testNeutrinoFermiDirac.c b/tests/testNeutrinoFermiDirac.c
new file mode 100644
index 0000000000000000000000000000000000000000..5d010795511aa69097eb9fdcbc830e37b1b0f3ef
--- /dev/null
+++ b/tests/testNeutrinoFermiDirac.c
@@ -0,0 +1,129 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2021 Willem Elbers (whe@willemelbers.com)
+ *
+ * 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/>.
+ *
+ ******************************************************************************/
+
+/* Some standard headers */
+#include <fenv.h>
+#include <math.h>
+
+/* Includes. */
+#include "../config.h"
+#include "neutrino/Default/fermi_dirac.h"
+#include "swift.h"
+
+/* Riemann function zeta(3) and zeta(5) */
+#define M_ZETA_3 1.2020569031595942853997
+#define M_ZETA_5 1.0369277551433699263314
+
+int main(int argc, char *argv[]) {
+ /* Exact integrals of x^n / (exp(x) + 1) on (0, infinity) */
+ double integral2 = M_ZETA_3 * 1.5;
+ double integral3 = M_PI * M_PI * M_PI * M_PI * 7.0 / 120.0;
+ double integral4 = M_ZETA_5 * 22.5;
+ double integral5 = M_PI * M_PI * M_PI * M_PI * M_PI * M_PI * 31. / 252.0;
+ /* Exact moments */
+ double mu = integral3 / integral2;
+ double mu2 = mu * mu;
+ double sigma2 = integral4 / integral2 - mu2;
+ double sigma = sqrt(sigma2);
+ double sigma3 = sigma2 * sigma;
+ double skewns = (integral5 / integral2 - 3 * mu * sigma2 - mu2 * mu) / sigma3;
+
+ /* Print FermiDirac(n) for n = 1, 2, 3 for external testing */
+ message("fermi_dirac(1) = %e", neutrino_seed_to_fermi_dirac(1));
+ message("fermi_dirac(2) = %e", neutrino_seed_to_fermi_dirac(2));
+ message("fermi_dirac(3) = %e\n", neutrino_seed_to_fermi_dirac(3));
+
+ /* Generate Fermi-Dirac numbers and compute the sum, min, and max */
+ int N = 1e7;
+ double sum = 0;
+ double min = FLT_MAX, max = 0;
+ long long seed = 290009001901;
+ for (int i = 0; i < N; i++) {
+ double x = neutrino_seed_to_fermi_dirac(seed + i);
+ sum += x;
+ if (x < min) min = x;
+ if (x > max) max = x;
+ }
+
+ /* Do a second pass for the sample variance and skewness */
+ double mean = sum / N;
+ double ss_sum = 0;
+ double sss_sum = 0;
+
+ /* We also construct a histogram */
+ int bins = 1000;
+ int *histogram1 = calloc(bins, sizeof(int));
+
+ /* Generate the same numbers again and compute statistics and histogram */
+ for (int i = 0; i < N; i++) {
+ double x = neutrino_seed_to_fermi_dirac(seed + i);
+ ss_sum += (x - mean) * (x - mean);
+ sss_sum += (x - mean) * (x - mean) * (x - mean);
+
+ int bin = (int)((x - min) / (max - min) * bins);
+ histogram1[bin]++;
+ }
+
+ /* Sample statistics */
+ double var = ss_sum / (N - 1);
+ double sdev = sqrt(var);
+ double mu3 = sss_sum / (N - 2) / (var * sdev);
+
+ /* Relative errors with the exact moments */
+ double err_mu = mean / mu - 1.;
+ double err_sig = var / sigma2 - 1.;
+ double err_skew = mu3 / skewns - 1.;
+
+ message("Sample mean is %f, exact: %f (rel. %e).", mean, mu, err_mu);
+ message("Sample variance is %f, exact: %f (rel. %e).", var, sigma2, err_sig);
+ message("Sample skewness is %f, exact: %f (rel. %e).", mu3, skewns, err_skew);
+
+ assert(fabs(err_mu) < 1e-3);
+ assert(fabs(err_sig) < 1e-2);
+ assert(fabs(err_skew) < 1e-2);
+
+ /* Construct Kolmogorov-Smirnov statistic by integrating the histogram */
+ double sum_empirical = 0.;
+ double sum_analytical = 0.;
+ double KS_statistic = 0.;
+ double dx = (max - min) / bins;
+ for (int bin = 0; bin < bins; bin++) {
+ double x = min + (bin + 0.5) * dx;
+ sum_empirical += histogram1[bin] * 1. / N;
+ sum_analytical += fermi_dirac_density(x) * x * x * dx / integral2;
+
+ double delta = fabs(sum_empirical - sum_analytical);
+ if (delta > KS_statistic) {
+ KS_statistic = delta;
+ }
+ }
+
+ /* Can we reject the hypothesis that these numbers are FD distributed? */
+ double crit_val = 5.146991e-05 * sqrt(1e9 / N); // 99% confidence level
+ message("KS statistic = %e (99%% that KS < %e)\n", KS_statistic, crit_val);
+
+ /* We should not reject this */
+ assert(KS_statistic < crit_val);
+
+ free(histogram1);
+
+ message("Success.");
+
+ return 0;
+}
diff --git a/tools/analyse_runtime.py b/tools/analyse_runtime.py
index 47f1d492a9790bff3f02d03c6da6a446654e8d03..6c40e51f0300f84b02b3d6be29eb1d6a018964ad 100755
--- a/tools/analyse_runtime.py
+++ b/tools/analyse_runtime.py
@@ -91,7 +91,8 @@ tasks = [
"fof",
"others",
"sink",
- "RT",
+ "neutrino",
+ "RT",
"total",
]