improved code documentation, added doxygen target

CMakeLists.txt

@@ -183,19 +183,42 @@ list (APPEND SOURCES
 endif()
 endif()
 
+########################################################################################################################
+# DOXYGEN
+# This will generate Doxygen documentation when you run make doc_doxygen in your build directory.
+option(BUILD_DOCUMENTATION "Build doxygen documentation" OFF)
+if(BUILD_DOCUMENTATION)
+    find_package(Doxygen)
+    if(DOXYGEN_FOUND)
+        set(DOXYGEN_IN ${CMAKE_CURRENT_SOURCE_DIR}/Doxyfile.in)
+        set(DOXYGEN_OUT ${CMAKE_CURRENT_BINARY_DIR}/Doxyfile)
+
+        configure_file(${DOXYGEN_IN} ${DOXYGEN_OUT} @ONLY)
+
+        add_custom_target( doc_doxygen ALL
+            COMMAND ${DOXYGEN_EXECUTABLE} ${DOXYGEN_OUT}
+            WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
+            COMMENT "Generating API documentation with Doxygen"
+            VERBATIM )
+    endif(DOXYGEN_FOUND)
+endif(BUILD_DOCUMENTATION)
+
+########################################################################################################################
 # project configuration header
 configure_file(
   ${PROJECT_SOURCE_DIR}/include/cmake_config.hh.in
   ${PROJECT_SOURCE_DIR}/include/cmake_config.hh
 )
 
+########################################################################################################################
+# executable and linking
 add_executable(${PRGNAME} ${SOURCES} ${PLUGINS})
 # target_setup_class(${PRGNAME})
 
 set_target_properties(${PRGNAME} PROPERTIES CXX_STANDARD 14)
 
-
+########################################################################################################################
-# mpi flags
+# mpi flags and precision set up
 if(MPI_CXX_FOUND)
   if(CODE_PRECISION STREQUAL "FLOAT")
     if(FFTW3_SINGLE_MPI_FOUND)

include/convolution.hh

@@ -22,6 +22,9 @@
 #include <general.hh>
 #include <grid_fft.hh>
 
+/// @brief base class for convolutions of two or three fields
+/// @tparam data_t 
+/// @tparam derived_t 
 template <typename data_t, typename derived_t>
 class BaseConvolver
 {
@@ -30,23 +33,44 @@ protected:
     std::array<real_t, 3> length_;
 
 public:
+
+    /// @brief Construct a new Base Convolver object
+    /// @param N linear grid size
+    /// @param L physical box size
     BaseConvolver(const std::array<size_t, 3> &N, const std::array<real_t, 3> &L)
         : np_(N), length_(L) {}
 
+
+    /// @brief Construct a new Base Convolver object [deleted copy constructor]  
     BaseConvolver( const BaseConvolver& ) = delete;
     
+    /// @brief destructor (virtual)
     virtual ~BaseConvolver() {}
 
-    // implements convolution of two Fourier-space fields
+    /// @brief implements convolution of two Fourier-space fields
+    /// @tparam kfunc1 field 1
+    /// @tparam kfunc2 field 2
+    /// @tparam opp output operator
     template <typename kfunc1, typename kfunc2, typename opp>
     void convolve2(kfunc1 kf1, kfunc2 kf2, opp op) {}
 
-    // implements convolution of three Fourier-space fields
+    /// @brief implements convolution of three Fourier-space fields
+    /// @tparam kfunc1 field 1
+    /// @tparam kfunc2 field 2
+    /// @tparam kfunc3 field 3
+    /// @tparam opp output operator
     template <typename kfunc1, typename kfunc2, typename kfunc3, typename opp>
     void convolve3(kfunc1 kf1, kfunc2 kf2, kfunc3 kf3, opp op) {}
 
 public:
 
+    /// @brief convolve two gradient fields in Fourier space a_{,i} * b_{,j}
+    /// @tparam opp output operator type
+    /// @param inl left input field a
+    /// @param d1l direction of first gradient (,i)
+    /// @param inr right input field b
+    /// @param d1r direction of second gradient (,j)
+    /// @param output_op output operator
     template <typename opp>
     void convolve_Gradients(Grid_FFT<data_t> &inl, const std::array<int, 1> &d1l,
                             Grid_FFT<data_t> &inr, const std::array<int, 1> &d1r,
@@ -55,19 +79,29 @@ public:
         // transform to FS in case fields are not
         inl.FourierTransformForward();
         inr.FourierTransformForward();
-        // perform convolution of Hessians
+        // perform convolution of two gradients
         static_cast<derived_t &>(*this).convolve2(
+            // first gradient
             [&inl,&d1l](size_t i, size_t j, size_t k) -> ccomplex_t {
                 auto grad1 = inl.gradient(d1l[0],{i,j,k});
                 return grad1*inl.kelem(i, j, k);
             },
+            // second gradient
             [&inr,&d1r](size_t i, size_t j, size_t k) -> ccomplex_t {
                 auto grad1 = inr.gradient(d1r[0],{i,j,k});
                 return grad1*inr.kelem(i, j, k);
             },
+            // -> output operator
             output_op);
     }
 
+    /// @brief convolve a gradient and a Hessian field in Fourier space a_{,i} * b_{,jk}
+    /// @tparam opp output operator type
+    /// @param inl left input field a
+    /// @param d1l direction of gradient (,i)
+    /// @param inr right input field b
+    /// @param d2r directions of Hessian (,jk)
+    /// @param output_op output operator
     template <typename opp>
     void convolve_Gradient_and_Hessian(Grid_FFT<data_t> &inl, const std::array<int, 1> &d1l,
                                        Grid_FFT<data_t> &inr, const std::array<int, 2> &d2r,
@@ -76,19 +110,29 @@ public:
         // transform to FS in case fields are not
         inl.FourierTransformForward();
         inr.FourierTransformForward();
-        // perform convolution of Hessians
+        // perform convolution of gradient and Hessian
         static_cast<derived_t &>(*this).convolve2(
+            // gradient
             [&](size_t i, size_t j, size_t k) -> ccomplex_t {
                 auto kk = inl.template get_k<real_t>(i, j, k);
                 return ccomplex_t(0.0, -kk[d1l[0]]) * inl.kelem(i, j, k);
             },
+            // Hessian
             [&](size_t i, size_t j, size_t k) -> ccomplex_t {
                 auto kk = inr.template get_k<real_t>(i, j, k);
                 return -kk[d2r[0]] * kk[d2r[1]] * inr.kelem(i, j, k);
             },
+            // -> output operator
             output_op);
     }
 
+    /// @brief convolve two Hessian fields in Fourier space a_{,ij} * b_{,kl}
+    /// @tparam opp output operator type
+    /// @param inl left input field a
+    /// @param d2l directions of first Hessian (,ij)
+    /// @param inr right input field b
+    /// @param d2r directions of second Hessian (,kl)
+    /// @param output_op output operator
     template <typename opp>
     void convolve_Hessians(Grid_FFT<data_t> &inl, const std::array<int, 2> &d2l,
                            Grid_FFT<data_t> &inr, const std::array<int, 2> &d2r,
@@ -99,19 +143,31 @@ public:
         inr.FourierTransformForward();
         // perform convolution of Hessians
         static_cast<derived_t &>(*this).convolve2(
+            // first Hessian
             [&inl,&d2l](size_t i, size_t j, size_t k) -> ccomplex_t {
                 auto grad1 = inl.gradient(d2l[0],{i,j,k});
                 auto grad2 = inl.gradient(d2l[1],{i,j,k});
                 return grad1*grad2*inl.kelem(i, j, k);
             },
+            // second Hessian
             [&inr,&d2r](size_t i, size_t j, size_t k) -> ccomplex_t {
                 auto grad1 = inr.gradient(d2r[0],{i,j,k});
                 auto grad2 = inr.gradient(d2r[1],{i,j,k});
                 return grad1*grad2*inr.kelem(i, j, k);
             },
+            // -> output operator
             output_op);
     }
 
+    /// @brief convolve three Hessian fields in Fourier space a_{,ij} * b_{,kl} * c_{,mn}
+    /// @tparam opp output operator
+    /// @param inl first input field a
+    /// @param d2l directions of first Hessian (,ij)
+    /// @param inm second input field b
+    /// @param d2m directions of second Hessian (,kl)
+    /// @param inr third input field c
+    /// @param d2r directions of third Hessian (,mn)
+    /// @param output_op output operator
     template <typename opp>
     void convolve_Hessians(Grid_FFT<data_t> &inl, const std::array<int, 2> &d2l,
                            Grid_FFT<data_t> &inm, const std::array<int, 2> &d2m,
@@ -124,24 +180,36 @@ public:
         inr.FourierTransformForward();
         // perform convolution of Hessians
         static_cast<derived_t &>(*this).convolve3(
+            // first Hessian
             [&inl, &d2l](size_t i, size_t j, size_t k) -> ccomplex_t {
                 auto grad1 = inl.gradient(d2l[0],{i,j,k});
                 auto grad2 = inl.gradient(d2l[1],{i,j,k});
                 return grad1*grad2*inl.kelem(i, j, k);
             },
+            // second Hessian
             [&inm, &d2m](size_t i, size_t j, size_t k) -> ccomplex_t {
                 auto grad1 = inm.gradient(d2m[0],{i,j,k});
                 auto grad2 = inm.gradient(d2m[1],{i,j,k});
                 return grad1*grad2*inm.kelem(i, j, k);
             },
+            // third Hessian
             [&inr, &d2r](size_t i, size_t j, size_t k) -> ccomplex_t {
                 auto grad1 = inr.gradient(d2r[0],{i,j,k});
                 auto grad2 = inr.gradient(d2r[1],{i,j,k});
                 return grad1*grad2*inr.kelem(i, j, k);
             },
+            // -> output operator
             output_op);
     }
 
+    /// @brief convolve Hessian field with sum of two Hessian fields in Fourier space a_{,ij} * (b_{,kl} + c_{,mn})
+    /// @tparam opp output operator type
+    /// @param inl left input field a
+    /// @param d2l directions of first Hessian (,ij)
+    /// @param inr right input field b
+    /// @param d2r1 directions of second Hessian (,kl)
+    /// @param d2r2 directions of third Hessian (,mn)
+    /// @param output_op output operator
     template <typename opp>
     void convolve_SumOfHessians(Grid_FFT<data_t> &inl, const std::array<int, 2> &d2l,
                                 Grid_FFT<data_t> &inr, const std::array<int, 2> &d2r1, const std::array<int, 2> &d2r2,
@@ -152,11 +220,13 @@ public:
         inr.FourierTransformForward();
         // perform convolution of Hessians
         static_cast<derived_t &>(*this).convolve2(
+            // first Hessian
             [&inl, &d2l](size_t i, size_t j, size_t k) -> ccomplex_t {
                 auto grad1 = inl.gradient(d2l[0],{i,j,k});
                 auto grad2 = inl.gradient(d2l[1],{i,j,k});
                 return grad1*grad2*inl.kelem(i, j, k);
             },
+            // second two Hessian and sum
             [&inr, &d2r1, &d2r2](size_t i, size_t j, size_t k) -> ccomplex_t {
                 auto grad11 = inr.gradient(d2r1[0],{i,j,k});
                 auto grad12 = inr.gradient(d2r1[1],{i,j,k});
@@ -164,9 +234,18 @@ public:
                 auto grad22 = inr.gradient(d2r2[1],{i,j,k});
                 return (grad11*grad12+grad21*grad22)*inr.kelem(i, j, k);
             },
+            // -> output operator
             output_op);
     }
 
+    /// @brief convolve Hessian field with difference of two Hessian fields in Fourier space a_{,ij} * (b_{,kl} - c_{,mn})
+    /// @tparam opp output operator type
+    /// @param inl left input field a
+    /// @param d2l directions of first Hessian (,ij)
+    /// @param inr right input field b
+    /// @param d2r1 directions of second Hessian (,kl)
+    /// @param d2r2 directions of third Hessian (,mn)
+    /// @param output_op output operator
     template <typename opp>
     void convolve_DifferenceOfHessians(Grid_FFT<data_t> &inl, const std::array<int, 2> &d2l,
                                        Grid_FFT<data_t> &inr, const std::array<int, 2> &d2r1, const std::array<int, 2> &d2r2,
@@ -177,11 +256,13 @@ public:
         inr.FourierTransformForward();
         // perform convolution of Hessians
         static_cast<derived_t &>(*this).convolve2(
+            // first Hessian
             [&inl, &d2l](size_t i, size_t j, size_t k) -> ccomplex_t {
                 auto grad1 = inl.gradient(d2l[0],{i,j,k});
                 auto grad2 = inl.gradient(d2l[1],{i,j,k});
                 return grad1*grad2*inl.kelem(i, j, k);
             },
+            // second two Hessian and difference
             [&inr, &d2r1, &d2r2](size_t i, size_t j, size_t k) -> ccomplex_t {
                 auto grad11 = inr.gradient(d2r1[0],{i,j,k});
                 auto grad12 = inr.gradient(d2r1[1],{i,j,k});
@@ -189,21 +270,29 @@ public:
                 auto grad22 = inr.gradient(d2r2[1],{i,j,k});
                 return (grad11*grad12-grad21*grad22)*inr.kelem(i, j, k);
             },
+            // -> output operator
             output_op);
     }
 };
 
-//! naive convolution class, disrespecting aliasing
+//! low-level implementation of convolutions -- naive convolution class, ignoring aliasing (no padding)
 template <typename data_t>
 class NaiveConvolver : public BaseConvolver<data_t, NaiveConvolver<data_t>>
 {
 protected:
+    /// @brief buffer for Fourier transformed fields
     Grid_FFT<data_t> *fbuf1_, *fbuf2_;
 
+    /// @brief number of points in each direction
     using BaseConvolver<data_t, NaiveConvolver<data_t>>::np_;
+
+    /// @brief length of each direction
     using BaseConvolver<data_t, NaiveConvolver<data_t>>::length_;
 
 public:
+    /// @brief constructor
+    /// @param N number of points in each direction
+    /// @param L length of each direction
     NaiveConvolver(const std::array<size_t, 3> &N, const std::array<real_t, 3> &L)
         : BaseConvolver<data_t, NaiveConvolver<data_t>>(N, L)
     {
@@ -211,12 +300,14 @@ public:
         fbuf2_ = new Grid_FFT<data_t>(N, length_, true, kspace_id);
     }
 
+    /// @brief destructor
     ~NaiveConvolver()
     {
         delete fbuf1_;
         delete fbuf2_;
     }
 
+    /// @brief convolution of two fields
     template <typename kfunc1, typename kfunc2, typename opp>
     void convolve2(kfunc1 kf1, kfunc2 kf2, opp output_op)
     {
@@ -249,6 +340,7 @@ public:
 
     }
 
+    /// @brief convolution of three fields
     template <typename kfunc1, typename kfunc2, typename kfunc3, typename opp>
     void convolve3(kfunc1 kf1, kfunc2 kf2, kfunc3 kf3, opp output_op)
     {
@@ -292,7 +384,13 @@ public:
         }
     }
 
+//--------------------------------------------------------------------------------------------------------
 private:
+
+    /// @brief copy data into a grid
+    /// @tparam kfunc abstract function type generating data
+    /// @param kf abstract function generating data
+    /// @param g grid to copy data into
     template <typename kfunc>
     void copy_in(kfunc kf, Grid_FFT<data_t> &g)
     {
@@ -310,23 +408,28 @@ private:
     }
 };
 
-//! convolution class, respecting Orszag's 3/2 rule
+//! convolution class, respecting Orszag's 3/2 rule (padding in Fourier space to avoid aliasing)
 template <typename data_t>
 class OrszagConvolver : public BaseConvolver<data_t, OrszagConvolver<data_t>>
 {
 private:
-    Grid_FFT<data_t> *f1p_, *f2p_;
+    /// @brief buffer for Fourier transformed fields
-    Grid_FFT<data_t> *fbuf_;
+    Grid_FFT<data_t> *f1p_, *f2p_, *fbuf_;
 
     using BaseConvolver<data_t, OrszagConvolver<data_t>>::np_;
     using BaseConvolver<data_t, OrszagConvolver<data_t>>::length_;
 
-    ccomplex_t *crecvbuf_;
+    ccomplex_t *crecvbuf_; //!< receive buffer for MPI (complex)
-    real_t *recvbuf_;
+    real_t *recvbuf_;     //!< receive buffer for MPI (real)
-    size_t maxslicesz_;
+    size_t maxslicesz_;  //!< maximum size of a slice
-    std::vector<ptrdiff_t> offsets_, offsetsp_;
+    std::vector<ptrdiff_t> offsets_, offsetsp_; //!< offsets for MPI
-    std::vector<size_t> sizes_, sizesp_;
+    std::vector<size_t> sizes_, sizesp_;       //!< sizes for MPI
 
+    /// @brief get task index for a given index
+    /// @param index index
+    /// @param offsets offsets
+    /// @param sizes sizes
+    /// @param ntasks number of tasks
     int get_task(ptrdiff_t index, const std::vector<ptrdiff_t> &offsets, const std::vector<size_t> &sizes, const int ntasks)
     {
         int itask = 0;
@@ -336,6 +439,10 @@ private:
     }
 
 public:
+
+    /// @brief constructor
+    /// @param N grid size
+    /// @param L grid length
     OrszagConvolver(const std::array<size_t, 3> &N, const std::array<real_t, 3> &L)
         : BaseConvolver<data_t, OrszagConvolver<data_t>>({3 * N[0] / 2, 3 * N[1] / 2, 3 * N[2] / 2}, L)
     {
@@ -370,6 +477,7 @@ public:
 #endif
     }
 
+    /// @brief destructor
     ~OrszagConvolver()
     {
         delete f1p_;
@@ -380,6 +488,10 @@ public:
 #endif
     }
 
+    /// @brief convolve two fields
+    /// @tparam kfunc1 abstract function type generating data for the first field
+    /// @tparam kfunc2 abstract function type generating data for the second field
+    /// @tparam opp abstract function type for the output operation
     template <typename kfunc1, typename kfunc2, typename opp>
     void convolve2(kfunc1 kf1, kfunc2 kf2, opp output_op)
     {
@@ -405,6 +517,11 @@ public:
         unpad(*f2p_, output_op);
     }
 
+    /// @brief convolve three fields
+    /// @tparam kfunc1 abstract function type generating data for the first field
+    /// @tparam kfunc2 abstract function type generating data for the second field
+    /// @tparam kfunc3 abstract function type generating data for the third field
+    /// @tparam opp abstract function type for the output operation
     template <typename kfunc1, typename kfunc2, typename kfunc3, typename opp>
     void convolve3(kfunc1 kf1, kfunc2 kf2, kfunc3 kf3, opp output_op)
     {
@@ -428,6 +545,10 @@ public:
 
 private:
 
+    /// @brief unpad the result of a convolution and copy it to a grid
+    /// @tparam kdep_functor abstract function type generating data for the result
+    /// @param kfunc abstract function generating data for the result
+    /// @param fp grid to copy the result to
     template <typename kdep_functor>
     void pad_insert( kdep_functor kfunc, Grid_FFT<data_t> &fp)
     {
@@ -472,6 +593,10 @@ private:
 #endif //defined(USE_MPI)
     }
 
+    /// @brief unpad the result of a convolution and write it to an output operator
+    /// @tparam operator_t abstract function type for the output operation
+    /// @param fp grid to copy the result from
+    /// @param output_op abstract function to write the result to
     template <typename operator_t>
     void unpad( Grid_FFT<data_t> &fp, operator_t output_op)
     {

include/cosmology_calculator.hh

@@ -152,8 +152,10 @@ private:
 
 public:
     
+    //! default constructor [deleted]
     calculator() = delete;
     
+    //! copy constructor [deleted]
     calculator(const calculator& c) = delete;
 
     //! constructor for a cosmology calculator object
@@ -204,9 +206,12 @@ public:
         m_sqrtpnorm_ = cosmo_param_["sqrtpnorm"];
     }
 
+    //! destructor
     ~calculator() { }
 
-    //! Write out a correctly scaled power spectrum at time a
+    /// @brief Write out a correctly scaled power spectrum at time a
+    /// @param a scale factor
+    /// @param fname file name
     void write_powerspectrum(real_t a, std::string fname) const
     {
         // const real_t Dplus0 = this->get_growth_factor(a);
@@ -253,7 +258,8 @@ public:
         music::ilog << "Wrote power spectrum at a=" << a << " to file \'" << fname << "\'" << std::endl;
     }
 
-    //! Write out a correctly scaled power spectrum at starting time
+    /// @brief Write out a correctly scaled transfer function at time a
+    /// @param[in] fname filename to write to
     void write_transfer( std::string fname ) const
     {
         // const real_t Dplus0 = this->get_growth_factor(a);
@@ -299,12 +305,16 @@ public:
         music::ilog << "Wrote input transfer functions at a=" << astart_ << " to file \'" << fname << "\'" << std::endl;
     }
 
+    /// @brief return the cosmological parameter object
+    /// @return cosmological parameter object
     const cosmology::parameters &get_parameters(void) const noexcept
     {
         return cosmo_param_;
     }
 
-    //! return the value of the Hubble function H(a) = dloga/dt 
+    /// @brief return the value of the Hubble function H(a) = dloga/dt 
+    /// @param[in] a scale factor
+    /// @return H(a)
     inline double H_of_a(double a) const noexcept
     {
         double HH2 = 0.0;
@@ -315,13 +325,17 @@ public:
         return cosmo_param_["H0"] * std::sqrt(HH2);
     }
 
-    //! Computes the linear theory growth factor D+, normalised to D+(a=1)=1
+    /// @brief Computes the linear theory growth factor D+, normalised to D+(a=1)=1
+    /// @param[in] a scale factor
+    /// @return D+(a)
     real_t get_growth_factor(real_t a) const noexcept
     {
         return D_of_a_(a) / Dnow_;
     }
 
-    //! Computes the inverse of get_growth_factor
+    /// @brief Computes the inverse of get_growth_factor, i.e. a(D+)
+    /// @param[in] Dplus growth factor
+    /// @return a(D+)
     real_t get_a( real_t Dplus ) const noexcept
     {
         return a_of_D_( Dplus * Dnow_ );

include/cosmology_parameters.hh

@@ -29,6 +29,7 @@ namespace cosmology
     class parameters
     {
     public:
+        //! type for default parameter maps
         using defaultmmap_t = std::map<std::string,std::map<std::string,real_t>>;
 
     private:

include/grid_fft.hh

@@ -25,26 +25,27 @@
 #include <bounding_box.hh>
 #include <typeinfo>
 
-enum space_t
+/// @brief enum to indicate whether a grid is currently in real or k-space
-{
+enum space_t { kspace_id, rspace_id };
-    kspace_id,
-    rspace_id
-};
-
 
 #ifdef USE_MPI
-template <typename data_t_, bool bdistributed=true>
+#define GRID_FFT_DISTRIBUTED true
 #else
-template <typename data_t_, bool bdistributed=false>
+#define GRID_FFT_DISTRIBUTED false
 #endif
+
+/// @brief class for FFTable grids
+/// @tparam data_t_ data type
+/// @tparam bdistributed flag to indicate whether this grid is distributed in memory
+template <typename data_t_, bool bdistributed=GRID_FFT_DISTRIBUTED>
 class Grid_FFT
 {
 public:
-    using data_t = data_t_;
+    using data_t = data_t_; ///< data type
-    static constexpr bool is_distributed_trait{bdistributed};
+    static constexpr bool is_distributed_trait{bdistributed}; ///< flag to indicate whether this grid is distributed in memory
 
 protected:
-    using grid_fft_t = Grid_FFT<data_t,bdistributed>;
+    using grid_fft_t = Grid_FFT<data_t,bdistributed>; ///< type of this grid
     
 public:
     std::array<size_t, 3> n_, nhalf_;
@@ -68,7 +69,11 @@ public:
     ptrdiff_t local_0_start_, local_1_start_;
     ptrdiff_t local_0_size_, local_1_size_;
 
-    //! constructor for FTable grid object
+    /// @brief constructor for FTable grid object
+    /// @param N number of grid points in each dimension
+    /// @param L physical size of the grid in each dimension
+    /// @param allocate flag to indicate whether to allocate memory for the grid
+    /// @param initialspace flag to indicate whether the grid is initially in real or k-space
     Grid_FFT(const std::array<size_t, 3> &N, const std::array<real_t, 3> &L, bool allocate = true, space_t initialspace = rspace_id)
         : n_(N), length_(L), space_(initialspace), data_(nullptr), cdata_(nullptr), plan_(nullptr), iplan_(nullptr), ballocated_( false )
     {
@@ -77,11 +82,16 @@ public:
         }
     }
 
-    // avoid implicit copying of data
+    /// @brief copy constructor [deleted] -- to avoid implicit copying of data
     Grid_FFT(const grid_fft_t &g) = delete;
 
+    /// @brief assignment operator [deleted] -- to avoid implicit copying of data
+    grid_fft_t &operator=(const grid_fft_t &g) = delete;
+
+    /// @brief destructor
     ~Grid_FFT() { reset(); }
 
+    /// @brief reset grid object (free memory, etc.)
     void reset()
     {
         if (data_ != nullptr)  { FFTW_API(free)(data_); data_ = nullptr; }
@@ -90,13 +100,18 @@ public:
         ballocated_ = false;
     }
 
+    /// @brief return the grid object for a given refinement level [dummy implementation for backward compatibility with MUSIC1]
     const grid_fft_t *get_grid(size_t ilevel) const { return this; }
 
-    //! return if grid object is 
+    /// @brief return if grid object is distributed in memory
+    /// @return true if grid object is distributed in memory
     bool is_distributed( void ) const noexcept { return bdistributed; }
 
+    /// @brief allocate memory for grid object
     void allocate();
 
+    /// @brief return if grid object is allocated
+    /// @return true if grid object is allocated
     bool is_allocated( void ) const noexcept { return ballocated_; }
 
     //! return the number of data_t elements that we store in the container

include/grid_ghosts.hh

@@ -18,12 +18,16 @@
 
 #include <array>
 #include <vector>
-
 #include <numeric>
 
 #include <general.hh>
 #include <math/vec3.hh>
 
+/// @brief implements a wrapper class for grids with ghost zones for MPI communication (slab decomposition)
+/// @tparam numghosts number of ghost zones on each side
+/// @tparam haveleft flag whether to have left ghost zone
+/// @tparam haveright flag whether to have right ghost zone
+/// @tparam grid_t grid type to wrap
 template <int numghosts, bool haveleft, bool haveright, typename grid_t>
 struct grid_with_ghosts
 {
@@ -43,6 +47,9 @@ struct grid_with_ghosts
   //... determine communication offsets
   std::vector<ptrdiff_t> offsets_, sizes_;
 
+  /// @brief get task index for a given index
+  /// @param index index
+  /// @return task index
   int get_task(ptrdiff_t index) const
   {
     int itask = 0;
@@ -51,6 +58,8 @@ struct grid_with_ghosts
     return itask;
   }
 
+  /// @brief constructor for grid with ghosts
+  /// @param g grid to wrap
   explicit grid_with_ghosts(const grid_t &g)
   : gridref(g), nx_(g.n_[0]), ny_(g.n_[1]), nz_(g.n_[2]), nzp_(g.n_[2]+2)
   {
@@ -74,6 +83,8 @@ struct grid_with_ghosts
     }
   }
 
+  /// @brief update ghost zones via MPI communication
+  /// @param g grid to wrap
   void update_ghosts_allow_multiple( const grid_t &g )
   {
   #if defined(USE_MPI)
@@ -181,11 +192,19 @@ struct grid_with_ghosts
     #endif
   }
 
+  /// @brief return the element at position (i,j,k) in the grid
+  /// @param i index in x direction
+  /// @param j index in y direction
+  /// @param k index in z direction
+  /// @return grid element at position (i,j,k)
   data_t relem(const ptrdiff_t& i, const ptrdiff_t& j, const ptrdiff_t&k ) const noexcept
   {
     return this->relem({i,j,k});
   }
 
+  /// @brief return the element at position (i,j,k) in the grid
+  /// @param pos position in the grid (array of size 3: {i,j,k})
+  /// @return grid element at position (i,j,k)
   data_t relem(const std::array<ptrdiff_t, 3> &pos) const noexcept
   {
     const ptrdiff_t ix = pos[0];

include/math/interpolate.hh

@@ -22,16 +22,22 @@
 #include <gsl/gsl_spline.h>
 #include <gsl/gsl_errno.h>
 
+
+/// @brief 1D interpolation class
+/// @tparam logx static flag to indicate logarithmic interpolation in x
+/// @tparam logy static flag to indicate logarithmic interpolation in y
+/// @tparam periodic static flag to indicate periodic interpolation in x
 template <bool logx, bool logy, bool periodic>
 class interpolated_function_1d
 {
 
 private:
-  bool isinit_;
+  bool isinit_; ///< flag to indicate whether the interpolation has been initialized
-  std::vector<double> data_x_, data_y_;
+  std::vector<double> data_x_, data_y_; ///< data vectors
-  gsl_interp_accel *gsl_ia_;
+  gsl_interp_accel *gsl_ia_; ///< GSL interpolation accelerator
-  gsl_spline *gsl_sp_;
+  gsl_spline *gsl_sp_; ///< GSL spline object
 
+  /// @brief deallocate GSL objects
   void deallocate()
   {
     gsl_spline_free(gsl_sp_);
@@ -39,10 +45,16 @@ private:
   }
 
 public:
+
+  /// @brief default copy constructor (deleted)
   interpolated_function_1d(const interpolated_function_1d &) = delete;
 
+  /// @brief empty constructor (without data)
   interpolated_function_1d() : isinit_(false){}
 
+  /// @brief constructor with data
+  /// @param data_x x data vector
+  /// @param data_y y data vector
   interpolated_function_1d(const std::vector<double> &data_x, const std::vector<double> &data_y)
   : isinit_(false)
   {
@@ -50,11 +62,15 @@ public:
     this->set_data( data_x, data_y );
   }
 
+  /// @brief destructor
   ~interpolated_function_1d()
   {
     if (isinit_) this->deallocate();
   }
 
+  /// @brief set data
+  /// @param data_x x data vector
+  /// @param data_y y data vector
   void set_data(const std::vector<double> &data_x, const std::vector<double> &data_y)
   {
     data_x_ = data_x;
@@ -75,6 +91,9 @@ public:
     isinit_ = true;
   }
 
+  /// @brief evaluate the interpolation
+  /// @param x x value
+  /// @return y value
   double operator()(double x) const noexcept
   {
     assert( isinit_ && !(logx&&x<=0.0) );

include/math/mat3.hh

@@ -22,18 +22,20 @@
 
 #include <math/vec3.hh>
 
-//! class for 3x3 matrix calculations
+/// @brief class for 3x3 matrix calculations
+/// @tparam T type of matrix elements
 template<typename T>
 class mat3_t{
 protected:
-    std::array<T,9> data_;
+    std::array<T,9> data_; //< data array
-    std::array<double,9> data_double_;
+    std::array<double,9> data_double_; //< data array for GSL operations
-    gsl_matrix_view m_;
+    gsl_matrix_view m_; //< GSL matrix view
-    gsl_vector *eval_;
+    gsl_vector *eval_; //< GSL eigenvalue vector
-    gsl_matrix *evec_;
+    gsl_matrix *evec_; //< GSL eigenvector matrix
-	gsl_eigen_symmv_workspace * wsp_;
+	gsl_eigen_symmv_workspace * wsp_; //< GSL workspace
-    bool bdid_alloc_gsl_;
+    bool bdid_alloc_gsl_; //< flag to indicate whether GSL memory has been allocated
 						
+    /// @brief initialize GSL memory
     void init_gsl(){
         // allocate memory for GSL operations if we haven't done so yet
         if( !bdid_alloc_gsl_ )
@@ -54,6 +56,7 @@ protected:
         }
     }
 
+    /// @brief free GSL memory
     void free_gsl(){
         // free memory for GSL operations if it was allocated
         if( bdid_alloc_gsl_ )
@@ -66,54 +69,76 @@ protected:
 
 public:
 
+    /// @brief default constructor
     mat3_t()
     : bdid_alloc_gsl_(false) 
     {}
 
-    //! copy constructor
+    /// @brief copy constructor
+    /// @param m matrix to copy
     mat3_t( const mat3_t<T> &m)
     : data_(m.data_), bdid_alloc_gsl_(false) 
     {}
     
-    //! move constructor
+    /// @brief move constructor
+    /// @param m matrix to move
     mat3_t( mat3_t<T> &&m)
     : data_(std::move(m.data_)), bdid_alloc_gsl_(false) 
     {}
 
-    //! construct mat3_t from initializer list
+    /// @brief construct mat3_t from initializer list
+    /// @param e initializer list
     template<typename ...E>
     mat3_t(E&&...e) 
     : data_{{std::forward<E>(e)...}}, bdid_alloc_gsl_(false)
     {}
 
+    /// @brief assignment operator
+    /// @param m matrix to copy
+    /// @return reference to this
     mat3_t<T>& operator=(const mat3_t<T>& m) noexcept{
         data_ = m.data_;
         return *this;
     }
 
+    /// @brief move assignment operator
+    /// @param m matrix to move
+    /// @return reference to this
     mat3_t<T>& operator=(const mat3_t<T>&& m) noexcept{
         data_ = std::move(m.data_);
         return *this;
     }
 
-    //! destructor
+    /// @brief destructor
     ~mat3_t(){
         this->free_gsl();
     }
     
-    //! bracket index access to vector components
+    /// @brief bracket index access to flattened matrix components
+    /// @param i index
+    /// @return reference to i-th component
     T &operator[](size_t i) noexcept { return data_[i];}
     
-    //! const bracket index access to vector components
+    /// @brief const bracket index access to flattened matrix components
+    /// @param i index
+    /// @return const reference to i-th component
     const T &operator[](size_t i) const noexcept { return data_[i]; }
 
-    //! matrix 2d index access
+    /// @brief matrix 2d index access
+    /// @param i row index
+    /// @param j column index
+    /// @return reference to (i,j)-th component
     T &operator()(size_t i, size_t j) noexcept { return data_[3*i+j]; }
 
-    //! const matrix 2d index access
+    /// @brief const matrix 2d index access
+    /// @param i row index
+    /// @param j column index
+    /// @return const reference to (i,j)-th component
     const T &operator()(size_t i, size_t j) const noexcept { return data_[3*i+j]; }
 
-    //! in-place addition
+    /// @brief in-place addition
+    /// @param rhs matrix to add
+    /// @return reference to this
     mat3_t<T>& operator+=( const mat3_t<T>& rhs ) noexcept{
         for (size_t i = 0; i < 9; ++i) {
            (*this)[i] += rhs[i];
@@ -121,7 +146,9 @@ public:
         return *this;
     }
 
-    //! in-place subtraction
+    /// @brief in-place subtraction
+    /// @param rhs matrix to subtract
+    /// @return reference to this
     mat3_t<T>& operator-=( const mat3_t<T>& rhs ) noexcept{
         for (size_t i = 0; i < 9; ++i) {
            (*this)[i] -= rhs[i];
@@ -129,10 +156,16 @@ public:
         return *this;
     }
 
+    /// @brief zeroing of matrix
     void zero() noexcept{
         for (size_t i = 0; i < 9; ++i) data_[i]=0;
     }
 
+    /// @brief compute eigenvalues and eigenvectors
+    /// @param evals eigenvalues
+    /// @param evec1 first eigenvector
+    /// @param evec2 second eigenvector
+    /// @param evec3 third eigenvector
     void eigen( vec3_t<T>& evals, vec3_t<T>& evec1, vec3_t<T>& evec2, vec3_t<T>& evec3_t )
     {
         this->init_gsl();
@@ -149,6 +182,11 @@ public:
     }
 };
 
+/// @brief matrix addition
+/// @tparam T type of matrix components
+/// @param lhs left hand side matrix
+/// @param rhs right hand side matrix
+/// @return matrix result = lhs + rhs
 template<typename T>
 constexpr const mat3_t<T> operator+(const mat3_t<T> &lhs, const mat3_t<T> &rhs) noexcept
 {
@@ -159,7 +197,11 @@ constexpr const mat3_t<T> operator+(const mat3_t<T> &lhs, const mat3_t<T> &rhs)
     return result;
 }
 
-// matrix - vector multiplication
+/// @brief matrix - vector multiplication
+/// @tparam T type of matrix and vector components
+/// @param A matrix
+/// @param v vector
+/// @return vector result = A*v
 template<typename T>
 inline vec3_t<T> operator*( const mat3_t<T> &A, const vec3_t<T> &v ) noexcept
 {

include/math/vec3.hh

@@ -17,7 +17,8 @@
 
 #pragma once
 
-//! implements a simple class of 3-vectors of arbitrary scalar type
+/// @brief implements a simple class of 3-vectors of arbitrary scalar type
+/// @tparam T scalar type
 template< typename T >
 class vec3_t{
 private:
@@ -28,19 +29,22 @@ public:
     //! expose access to elements via references
     T &x,&y,&z;
 
-    //! empty constructor
+    /// @brief empty constructor
     vec3_t()
     : data_{{T(0),T(0),T(0)}},x(data_[0]),y(data_[1]),z(data_[2]){}
 
-    //! copy constructor
+    /// @brief copy constructor
+    /// @param v vector to copy from
     vec3_t( const vec3_t<T> &v)
     : data_(v.data_), x(data_[0]),y(data_[1]),z(data_[2]){}
 
-    //! copy constructor for non-const reference, needed to avoid variadic template being called for non-const reference
+    /// @brief copy constructor for non-const reference, needed to avoid variadic template being called for non-const reference
+    /// @param v vector to copy from
     vec3_t( vec3_t<T>& v)
     : data_(v.data_), x(data_[0]),y(data_[1]),z(data_[2]){}
 
-    //! move constructor
+    /// @brief move constructor
+    /// @param v vector to move from
     vec3_t( vec3_t<T> &&v)
     : data_(std::move(v.data_)), x(data_[0]), y(data_[1]), z(data_[2]){}