Formatting

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));

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.",

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;

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);