diff --git a/examples/main.c b/examples/main.c
index 081374f4fda3f14ebaab7035a093163b59b9db63..a16464e658e13fb724fcd96800cbfb2bdaee806a 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -1303,8 +1303,8 @@ int main(int argc, char *argv[]) {
/* Initialize the neutrino mesh if used */
bzero(&neutrino_mesh, sizeof(struct neutrino_mesh));
if (neutrino_properties.use_linear_response)
- neutrino_mesh_init(params, &us, s.dim, &cosmo, &neutrino_properties,
- &gravity_properties, &neutrino_mesh, verbose);
+ neutrino_mesh_init(&neutrino_mesh, params, &us, s.dim, &cosmo,
+ &neutrino_properties, &gravity_properties, verbose);
/* Initialise the external potential properties */
bzero(&potential, sizeof(struct external_potential));
diff --git a/src/mesh_gravity.c b/src/mesh_gravity.c
index 4b1293ddf32d9a1caafbfce749b0e40257f1878c..a6202ceb9add2584925c949932db60034af93e5f 100644
--- a/src/mesh_gravity.c
+++ b/src/mesh_gravity.c
@@ -759,7 +759,7 @@ void compute_potential_distributed(struct pm_mesh* mesh, const struct space* s,
/* If using linear response neutrinos, apply to local slice of the MPI mesh */
if (s->e->neutrino_properties->use_linear_response)
- neutrino_mesh_compute(s, mesh, tp, frho, verbose, local_0_start, local_n0);
+ neutrino_mesh_compute(s, mesh, tp, frho, local_0_start, local_n0, verbose);
if (verbose)
message("Applying neutrino response took %.3f %s.",
@@ -957,8 +957,8 @@ void compute_potential_global(struct pm_mesh* mesh, const struct space* s,
/* If using linear response neutrinos, apply the response to the mesh */
if (s->e->neutrino_properties->use_linear_response)
- neutrino_mesh_compute(s, mesh, tp, frho, verbose, /*slice_offset=*/0,
- /*slice_width=*/N);
+ neutrino_mesh_compute(s, mesh, tp, frho, /*slice_offset=*/0,
+ /*slice_width=*/N, verbose);
if (verbose)
message("Applying neutrino response took %.3f %s.",
diff --git a/src/neutrino/Default/neutrino_mesh_linres.c b/src/neutrino/Default/neutrino_mesh_linres.c
index df24add2f57c8544ed7e9e0653a4e4cec4820e7c..81ae5330bce06d2ced345d2b55f2334311b68086 100644
--- a/src/neutrino/Default/neutrino_mesh_linres.c
+++ b/src/neutrino/Default/neutrino_mesh_linres.c
@@ -35,6 +35,13 @@ const hsize_t timestep_length = 10000;
/*! Number of wavenumbers in the neutrino density interpolation table */
const hsize_t wavenumber_length = 10000;
+/**
+ * @brief Determine the length of an HDF5 dataset
+ *
+ * @param h_file The file object
+ * @param title The title of the dataset
+ * @param length (Output) The vector length
+ */
void read_vector_length(hid_t h_file, char title[PARSER_MAX_LINE_SIZE],
hsize_t *length) {
@@ -60,6 +67,15 @@ void read_vector_length(hid_t h_file, char title[PARSER_MAX_LINE_SIZE],
H5Dclose(h_data);
}
+/**
+ * @brief Read transfer funtion from an HDF5 dataset with expected size N_z*N_k
+ *
+ * @param h_file The file object
+ * @param title The title of the dataset
+ * @param length (Output) Memory fot the transfer function data
+ * @param N_z Expected number of timesteps
+ * @param N_k Expected number of wavenumbers
+ */
void read_transfer_function(hid_t h_file, char title[PARSER_MAX_LINE_SIZE],
double *dest, hsize_t N_z, hsize_t N_k) {
@@ -91,12 +107,26 @@ void read_transfer_function(hid_t h_file, char title[PARSER_MAX_LINE_SIZE],
H5Dclose(h_data);
}
-void neutrino_mesh_init(struct swift_params *params,
+/**
+ * @brief Initialize the #neutrino_mesh object by reading an HDF5 file with
+ * transfer functions and pre-computing a transfer function ratio.
+ *
+ * @param numesh The #neutrino_mesh to be initialized
+ * @param params The parsed parameter file.
+ * @param us The system of units used internally.
+ * @param dim Spatial dimensions of the domain.
+ * @param c The #cosmology used for this run.
+ * @param np The #neutrino_props used for this run.
+ * @param gp The #gravity_props used for this run.
+ * @param verbose Are we talkative ?
+ */
+void neutrino_mesh_init(struct neutrino_mesh *numesh,
+ struct swift_params *params,
const struct unit_system *us, const double dim[3],
const struct cosmology *c,
const struct neutrino_props *np,
const struct gravity_props *gp,
- struct neutrino_mesh *numesh, int verbose) {
+ int verbose) {
/* Do we need to do anything? */
if (!np->use_linear_response) return;
@@ -265,7 +295,7 @@ void neutrino_mesh_init(struct swift_params *params,
message("(a_min, a_max) = (%g, %g)", a_min, a_max);
}
- /* We will remap the data such that it exactly covers the required domain
+ /* We will remap the data such that it just covers the required domain
* with a constant log spacing, allowing for faster interpolation.
* We only use the slower GSL interpolation for the remapping. */
@@ -373,7 +403,7 @@ struct neutrino_mesh_tp_data {
/* Background and perturbed density ratios */
double bg_density_ratio;
- double *ncdm_over_cb_arr;
+ const double *ncdm_over_cb_arr;
};
/**
@@ -454,10 +484,24 @@ void neutrino_mesh_apply_neutrino_response_mapper(void *map_data, const int num,
}
}
+/**
+ * @brief Apply the linear neutrino response to the Fourier transform of the
+ * gravitational potential.
+ *
+ * @param s The current #space
+ * @param mesh The #pm_mesh used to store the potential
+ * @param tp The #threadpool object used for parallelisation
+ * @param frho The NxNx(N/2) complex array of the Fourier transform of the
+ * density field
+ * @param slice_offset The x coordinate of the start of the slice on this MPI
+ * rank
+ * @param slice_width The width of the local slice on this MPI rank
+ * @param verbose Are we talkative?
+ */
void neutrino_mesh_compute(const struct space *s, struct pm_mesh *mesh,
struct threadpool *tp, fftw_complex *frho,
- int verbose, const int slice_offset,
- const int slice_width) {
+ const int slice_offset, const int slice_width,
+ int verbose) {
#ifdef HAVE_FFTW
const struct cosmology *c = s->e->cosmology;
diff --git a/src/neutrino/Default/neutrino_mesh_linres.h b/src/neutrino/Default/neutrino_mesh_linres.h
index 70f00c760d49af66684018a2a179a989fa9d4703..11b2e11864ab4a13beb3fdaaddd43c2a9b84bfb5 100644
--- a/src/neutrino/Default/neutrino_mesh_linres.h
+++ b/src/neutrino/Default/neutrino_mesh_linres.h
@@ -56,17 +56,18 @@ struct neutrino_mesh {
double *ncdm_over_cb;
};
-void neutrino_mesh_init(struct swift_params *params,
+void neutrino_mesh_init(struct neutrino_mesh *numesh,
+ struct swift_params *params,
const struct unit_system *us, const double dim[3],
const struct cosmology *c,
const struct neutrino_props *np,
const struct gravity_props *gp,
- struct neutrino_mesh *numesh, int verbose);
+ int verbose);
void neutrino_mesh_clean(struct neutrino_mesh *numesh);
void neutrino_mesh_compute(const struct space *s, struct pm_mesh *mesh,
struct threadpool *tp, fftw_complex *frho,
- int verbose, const int slice_offset,
- const int slice_width);
+ const int slice_offset, const int slice_width,
+ int verbose);
void neutrino_mesh_struct_dump(const struct neutrino_mesh *numesh,
FILE *stream);
void neutrino_mesh_struct_restore(struct neutrino_mesh *numesh, FILE *stream);