diff --git a/AGENTS.md b/AGENTS.md
new file mode 100644
index 000000000..88bdddc35
--- /dev/null
+++ b/AGENTS.md
@@ -0,0 +1,22 @@
+# AGENTS.md
+
+## Repository purpose
+IPPL is a C++ Particle and fields Framework. Preserve physical correctness over stylistic cleanup.
+
+## Build
+- Configure with CMake using the project’s normal toolchain and dependency prefixes.
+- Build in the local build or build_openmp directory.
+- build instructions can be found online.
+- Run the relevant test subset after changes.
+- External dependencies are in the _deps directory fetched with cmake
+
+## External dependencies
+- heFFte and Kokkos are external dependencies.
+- Prefer changing IPPL adapter/wrapper code over patching upstream dependencies.
+- When a change touches parallel kernels, explain execution space, memory space, and data movement implications.
+
+## Physics / numerics rules
+- For algorithmic changes, state expected impact on conservation, stability, and reproducibility.
+- Flag any change that may alter floating-point behavior or MPI/GPU execution order.
+- Add or update at least one regression/sanity test for physics-facing changes.
+- Keep am eye on parallel efficincy at least on the openmp level
\ No newline at end of file
diff --git a/scripts/pre-commit b/scripts/pre-commit
index 1517055b5..a05e84c01 100755
--- a/scripts/pre-commit
+++ b/scripts/pre-commit
@@ -14,7 +14,16 @@
 # from the .git/hooks directory, so the symlink has to be redirected)
 # ln -s -f ../../scripts/pre-commit .git/hooks/pre-commit
 
-CLANG_FORMAT_VERSION=clang-format-21
+CLANG_FORMAT_VERSION=${CLANG_FORMAT_VERSION:-clang-format-21}
+
+if ! command -v "$CLANG_FORMAT_VERSION" >/dev/null 2>&1; then
+    if command -v clang-format >/dev/null 2>&1; then
+        CLANG_FORMAT_VERSION=clang-format
+    else
+        printf "clang-format executable not found. Install clang-format-21 or set CLANG_FORMAT_VERSION.\n"
+        exit 1
+    fi
+fi
 
 red=$(tput setaf 1)
 green=$(tput setaf 2)
@@ -25,9 +34,33 @@ normal=$(tput sgr0)
 cxx_match="c|cc|cp|cxx|cpp|c\+\+|h|hh|hpp|hxx|h\+\+|ipp|tpp|tcc|txx|inl|inc"
 cmake_match="CMakeLists\.txt|\.cmake"
 
+relative_to_git_prefix() {
+    local path="$1"
+    local prefix="${GIT_PREFIX:-}"
+
+    if [ -n "$prefix" ]; then
+        case "$path" in
+            "$prefix"*) printf "%s\n" "${path#$prefix}" ;;
+            *) printf "%s\n" "$path" ;;
+        esac
+    else
+        printf "%s\n" "$path"
+    fi
+}
+
+clang_assume_filename() {
+    local path="$1"
+
+    case "$path" in
+        *.h) printf "%s.cpp\n" "$path" ;;
+        *) printf "%s\n" "$path" ;;
+    esac
+}
+
 cxxfiles=()
 for file in `git diff --cached --name-only --diff-filter=ACMRT | grep -E "\.(${cxx_match})$"`; do
-    if ! cmp -s <(git show :${file}) <(git show :${file}|$CLANG_FORMAT_VERSION -style=file --assume-filename="${file}"); then
+    assume_filename=$(clang_assume_filename "$file")
+    if ! cmp -s <(git show :${file}) <(git show :${file}|$CLANG_FORMAT_VERSION -style=file --assume-filename="${assume_filename}"); then
         cxxfiles+=("${file}")
     fi
 done
@@ -49,7 +82,7 @@ full_list=
 if [ -n "${cxxfiles}" ]; then
     printf "# ${blue}clang-format ${red}error pre-commit${normal} : To fix run the following (use git commit ${yellow}--no-verify${normal} to bypass)\n"
     for f in "${cxxfiles[@]}" ; do
-        rel=$(realpath --relative-to "./$GIT_PREFIX" "$f")
+        rel=$(relative_to_git_prefix "$f")
         printf "$CLANG_FORMAT_VERSION -style=file -i %s\n" "$rel"
         full_list="${rel} ${full_list}"
     done
@@ -59,7 +92,7 @@ fi
 if [ -n "${cmakefiles}" ]; then
     printf "# ${green}cmake-format ${red}error pre-commit${normal} : To fix run the following (use git commit ${yellow}--no-verify${normal} to bypass)\n"
     for f in "${cmakefiles[@]}" ; do
-        rel=$(realpath --relative-to "./$GIT_PREFIX" "$f")
+        rel=$(relative_to_git_prefix "$f")
         printf "cmake-format -i %s\n" "$rel"
         full_list="${rel} ${full_list}"
     done
diff --git a/src/Communicate/BufferHandler.h b/src/Communicate/BufferHandler.h
index 117c12f22..9975a3631 100644
--- a/src/Communicate/BufferHandler.h
+++ b/src/Communicate/BufferHandler.h
@@ -121,22 +121,22 @@ namespace ippl {
         void freeBuffer(buffer_type buffer) override;
 
         /**
-         * @copydoc BufferHandler::freeBuffer
+         * @copydoc BufferHandler::freeAllBuffers
          */
         void freeAllBuffers() override;
 
         /**
-         * @copydoc BufferHandler::freeBuffer
+         * @copydoc BufferHandler::deleteAllBuffers
          */
         void deleteAllBuffers() override;
 
         /**
-         * @copydoc BufferHandler::freeBuffer
+         * @copydoc BufferHandler::getUsedSize
          */
         size_type getUsedSize() const override;
 
         /**
-         * @copydoc BufferHandler::freeBuffer
+         * @copydoc BufferHandler::getFreeSize
          */
         size_type getFreeSize() const override;
 
diff --git a/src/Communicate/Communicator.cpp b/src/Communicate/Communicator.cpp
index 55d601715..2c33f939f 100644
--- a/src/Communicate/Communicator.cpp
+++ b/src/Communicate/Communicator.cpp
@@ -25,7 +25,7 @@ namespace ippl::mpi {
         return *this;
     }
 
-    Communicator Communicator::Communicator::split(int color, int key) const {
+    Communicator Communicator::split(int color, int key) const {
         MPI_Comm newcomm;
         MPI_Comm_split(*comm_m, color, key, &newcomm);
         return Communicator(newcomm);
diff --git a/src/Communicate/Operations.h b/src/Communicate/Operations.h
index 0096b70ca..98b9c2be3 100644
--- a/src/Communicate/Operations.h
+++ b/src/Communicate/Operations.h
@@ -91,14 +91,9 @@ namespace ippl {
                 constexpr binaryOperationKind opKind = extractBinaryOperationKind<Op>::value;
                 MPI_Op ret;
                 MPI_Op_create(
-                    /**
-                     * @brief Construct a new lambda object without captures, therefore convertible
-                     * to a function pointer
-                     *
-                     * @param inputBuffer pointing to a Type object
-                     * @param outputBuffer pointing to a Type object
-                     * @param len Amount of _Type objects_! NOT amount of bytes!
-                     */
+                    // Captureless lambda, therefore convertible to a function pointer.
+                    // inputBuffer and outputBuffer point to Type objects; len is the number of
+                    // Type objects, not bytes.
                     [](void* inputBuffer, void* outputBuffer, int* len, MPI_Datatype*) {
                         Type* input = (Type*)inputBuffer;
 
diff --git a/src/FEM/Elements/EdgeElement.h b/src/FEM/Elements/EdgeElement.h
index 15dece5f2..99cc9a58f 100644
--- a/src/FEM/Elements/EdgeElement.h
+++ b/src/FEM/Elements/EdgeElement.h
@@ -59,7 +59,8 @@ namespace ippl {
          *
          * @return point_t
          */
-        KOKKOS_FUNCTION point_t globalToLocal(const vertex_points_t&, const point_t&) const;
+        KOKKOS_FUNCTION point_t globalToLocal(const vertex_points_t& global_vertices,
+                                              const point_t& point) const;
 
         /**
          * @brief Transforms a point from local to global coordinates.
diff --git a/src/FEM/Elements/HexahedralElement.h b/src/FEM/Elements/HexahedralElement.h
index 1a659f124..4c08d27a2 100644
--- a/src/FEM/Elements/HexahedralElement.h
+++ b/src/FEM/Elements/HexahedralElement.h
@@ -58,7 +58,8 @@ namespace ippl {
          *
          * @return point_t
          */
-        KOKKOS_FUNCTION point_t globalToLocal(const vertex_points_t&, const point_t&) const;
+        KOKKOS_FUNCTION point_t globalToLocal(const vertex_points_t& global_vertices,
+                                              const point_t& point) const;
 
         /**
          * @brief Transforms a point from local to global coordinates.
diff --git a/src/FEM/Elements/QuadrilateralElement.h b/src/FEM/Elements/QuadrilateralElement.h
index ebd1d7085..d7f7f5b18 100644
--- a/src/FEM/Elements/QuadrilateralElement.h
+++ b/src/FEM/Elements/QuadrilateralElement.h
@@ -59,7 +59,8 @@ namespace ippl {
          *
          * @return point_t
          */
-        KOKKOS_FUNCTION point_t globalToLocal(const vertex_points_t&, const point_t&) const;
+        KOKKOS_FUNCTION point_t globalToLocal(const vertex_points_t& global_vertices,
+                                              const point_t& point) const;
 
         /**
          * @brief Transforms a point from local to global coordinates.
diff --git a/src/FEM/LagrangeSpace.h b/src/FEM/LagrangeSpace.h
index ef980a48c..dcbfdd559 100644
--- a/src/FEM/LagrangeSpace.h
+++ b/src/FEM/LagrangeSpace.h
@@ -112,7 +112,7 @@ namespace ippl {
 
         ///////////////////////////////////////////////////////////////////////
         /**
-         * @brief Function to update the element partition and the layout of 
+         * @brief Function to update the element partition and the layout of
          * fields in the LagrangeSpace if the layout has been changed during
          * the simulation (for example by the load balancer).
          */
@@ -138,7 +138,7 @@ namespace ippl {
          * @return size_t - The local DOF index
          */
         KOKKOS_FUNCTION size_t getLocalDOFIndex(const size_t& elementIndex,
-                                 const size_t& globalDOFIndex) const override;
+                                                const size_t& globalDOFIndex) const override;
 
         /**
          * @brief Get the global DOF index from the element index and local DOF
@@ -163,7 +163,7 @@ namespace ippl {
         /**
          * @brief Get the global DOF indices (vector of global DOF indices) of an element
          *
-         * @param elementIndex size_t - The index of the element
+         * @param element_index size_t - The index of the element
          *
          * @return Vector<size_t, NumElementDOFs> - The global DOF indices
          */
@@ -203,7 +203,8 @@ namespace ippl {
         /// Functions to access element info from outside /////////////////////
         ///////////////////////////////////////////////////////////////////////
 
-        KOKKOS_FUNCTION point_t getInverseTransposeTransformationJacobian(vertex_points_t pt) const {
+        KOKKOS_FUNCTION point_t
+        getInverseTransposeTransformationJacobian(vertex_points_t pt) const {
             return this->ref_element_m.getInverseTransposeTransformationJacobian(pt);
         }
 
@@ -238,8 +239,8 @@ namespace ippl {
         FieldLHS evaluateAx_diag(FieldLHS& field, F& evalFunction);
 
         /**
-         * @brief Assemble the left stiffness matrix A of the system 
-         * but only for the boundary values, so that they can be 
+         * @brief Assemble the left stiffness matrix A of the system
+         * but only for the boundary values, so that they can be
          * subtracted from the RHS for treatment of Dirichlet BCs
          *
          * @param field The field to assemble the matrix for
@@ -253,9 +254,7 @@ namespace ippl {
         /**
          * @brief Assemble the load vector b of the system Ax = b
          *
-         * @param rhs_field The field to set with the load vector
-         *
-         * @return FieldRHS - The RHS field containing b
+         * @param field The field to set with the load vector
          */
         void evaluateLoadVector(FieldRHS& field) const;
         void evaluateLumpedMass(FieldRHS& field) const;
@@ -298,17 +297,17 @@ namespace ippl {
             KOKKOS_FUNCTION indices_t getMeshVertexNDIndex(const size_t& vertex_index) const;
 
             KOKKOS_FUNCTION size_t getLocalDOFIndex(const indices_t& elementNDIndex,
-                const size_t& globalDOFIndex) const;
+                                                    const size_t& globalDOFIndex) const;
             KOKKOS_FUNCTION Vector<size_t, numElementDOFs> getGlobalDOFIndices(
                 const indices_t& elementNDIndex) const;
 
             KOKKOS_FUNCTION T evaluateRefElementShapeFunction(const size_t& localDOF,
-                const point_t& localPoint) const;
+                                                              const point_t& localPoint) const;
             KOKKOS_FUNCTION point_t evaluateRefElementShapeFunctionGradient(
                 const size_t& localDOF, const point_t& localPoint) const;
         };
 
-        DeviceStruct getDeviceMirror() const; 
+        DeviceStruct getDeviceMirror() const;
 
     private:
         /**
diff --git a/src/FEM/NedelecSpace.h b/src/FEM/NedelecSpace.h
index d272b43a0..c023cbf34 100644
--- a/src/FEM/NedelecSpace.h
+++ b/src/FEM/NedelecSpace.h
@@ -7,13 +7,13 @@
 
 #include <cmath>
 
-#include "FEM/FiniteElementSpace.h"
 #include "FEM/FEMVector.h"
+#include "FEM/FiniteElementSpace.h"
 
 constexpr unsigned getNedelecNumElementDOFs(unsigned Dim, [[maybe_unused]] unsigned Order) {
     // needs to be constexpr pow function to work at compile time. Kokkos::pow
     // doesn't work.
-    return static_cast<unsigned>(static_cast<int>(Dim)*power(2, static_cast<int>(Dim-1)));
+    return static_cast<unsigned>(static_cast<int>(Dim) * power(2, static_cast<int>(Dim - 1)));
 }
 
 namespace ippl {
@@ -31,22 +31,25 @@ namespace ippl {
     template <typename T, unsigned Dim, unsigned Order, typename ElementType,
               typename QuadratureType, typename FieldType>
     // requires IsQuadrature<QuadratureType>
-    class NedelecSpace : public FiniteElementSpace<T, Dim, getNedelecNumElementDOFs(Dim, Order),
-            ElementType, QuadratureType, FEMVector<T>, FEMVector<T>> {
+    class NedelecSpace
+        : public FiniteElementSpace<T, Dim, getNedelecNumElementDOFs(Dim, Order), ElementType,
+                                    QuadratureType, FEMVector<T>, FEMVector<T>> {
     public:
         // The number of degrees of freedom per element
         static constexpr unsigned numElementDOFs = getNedelecNumElementDOFs(Dim, Order);
 
         // The dimension of the mesh
-        static constexpr unsigned dim = FiniteElementSpace<T, Dim, numElementDOFs, ElementType,
-                                            QuadratureType, FEMVector<T>, FEMVector<T>>::dim;
+        static constexpr unsigned dim =
+            FiniteElementSpace<T, Dim, numElementDOFs, ElementType, QuadratureType, FEMVector<T>,
+                               FEMVector<T>>::dim;
 
         // The order of the Nedelec space
         static constexpr unsigned order = Order;
 
         // The number of mesh vertices per element
-        static constexpr unsigned numElementVertices = FiniteElementSpace<T, Dim, numElementDOFs,
-                        ElementType, QuadratureType, FEMVector<T>, FEMVector<T>>::numElementVertices;
+        static constexpr unsigned numElementVertices =
+            FiniteElementSpace<T, Dim, numElementDOFs, ElementType, QuadratureType, FEMVector<T>,
+                               FEMVector<T>>::numElementVertices;
 
         // A vector with the position of the element in the mesh in each dimension
         typedef typename FiniteElementSpace<T, Dim, numElementDOFs, ElementType, QuadratureType,
@@ -57,17 +60,16 @@ namespace ippl {
                                             FEMVector<T>, FEMVector<T>>::point_t point_t;
 
         typedef typename FiniteElementSpace<T, Dim, numElementDOFs, ElementType, QuadratureType,
-                                            FEMVector<T>, FEMVector<T>>::vertex_points_t vertex_points_t;
+                                            FEMVector<T>, FEMVector<T>>::vertex_points_t
+            vertex_points_t;
 
         // Field layout type for domain decomposition info
         typedef FieldLayout<Dim> Layout_t;
 
         // View types
         typedef typename detail::ViewType<T, 1>::view_type ViewType;
-        typedef typename detail::ViewType<T, 1,
-                                Kokkos::MemoryTraits<Kokkos::Atomic>>::view_type AtomicViewType;
-
-
+        typedef typename detail::ViewType<T, 1, Kokkos::MemoryTraits<Kokkos::Atomic>>::view_type
+            AtomicViewType;
 
         ///////////////////////////////////////////////////////////////////////
         // Constructors ///////////////////////////////////////////////////////
@@ -82,7 +84,7 @@ namespace ippl {
          * @param layout Reference to the field layout
          */
         NedelecSpace(UniformCartesian<T, Dim>& mesh, ElementType& ref_element,
-                                    const QuadratureType& quadrature, const Layout_t& layout);
+                     const QuadratureType& quadrature, const Layout_t& layout);
 
         /**
          * @brief Construct a new NedelecSpace object (without layout)
@@ -94,7 +96,7 @@ namespace ippl {
          * @param quadrature Reference to the quadrature rule
          */
         NedelecSpace(UniformCartesian<T, Dim>& mesh, ElementType& ref_element,
-                                    const QuadratureType& quadrature);
+                     const QuadratureType& quadrature);
 
         /**
          * @brief Initialize a NedelecSpace object created with the default
@@ -140,7 +142,7 @@ namespace ippl {
          * @return size_t - The local DOF index
          */
         KOKKOS_FUNCTION size_t getLocalDOFIndex(const size_t& elementIndex,
-                                    const size_t& globalDOFIndex) const override;
+                                                const size_t& globalDOFIndex) const override;
 
         /**
          * @brief Get the global DOF index from the element index and local DOF.
@@ -151,7 +153,7 @@ namespace ippl {
          * @return size_t - The global DOF index
          */
         KOKKOS_FUNCTION size_t getGlobalDOFIndex(const size_t& elementIndex,
-                                    const size_t& localDOFIndex) const override;
+                                                 const size_t& localDOFIndex) const override;
 
         /**
          * @brief Get the local DOF indices (vector of local DOF indices)
@@ -171,8 +173,8 @@ namespace ippl {
          * @return Vector<size_t, NumElementDOFs> - The global DOF indices
          */
         KOKKOS_FUNCTION Vector<size_t, numElementDOFs> getGlobalDOFIndices(
-                                    const size_t& elementIndex) const override;
-        
+            const size_t& elementIndex) const override;
+
         /**
          * @brief Get the global DOF indices (vector of global DOF indices) of
          * an element.
@@ -182,38 +184,39 @@ namespace ippl {
          * @return Vector<size_t, NumElementDOFs> - The global DOF indices
          */
         KOKKOS_FUNCTION Vector<size_t, numElementDOFs> getGlobalDOFIndices(
-                                    const indices_t& elementIndex) const;
-        
+            const indices_t& elementIndex) const;
+
         /**
          * @brief Get the DOF indices (vector of indices) corresponding to the
          * position inside the FEMVector of an element
          *
          * @param elementIndex size_t - The index of the element
+         * @param ldom local domain used for FEMVector indexing
          *
          * @return Vector<size_t, NumElementDOFs> - The DOF indices
          */
         KOKKOS_FUNCTION Vector<size_t, numElementDOFs> getFEMVectorDOFIndices(
-                                    const size_t& elementIndex, NDIndex<Dim> ldom) const;
+            const size_t& elementIndex, NDIndex<Dim> ldom) const;
 
         /**
          * @brief Get the DOF indices (vector of indices) corresponding to the
          * position inside the FEMVector of an element
          *
          * @param elementIndex indices_t - The index of the element
+         * @param ldom local domain used for FEMVector indexing
          *
          * @return Vector<size_t, NumElementDOFs> - The DOF indices
          */
         KOKKOS_FUNCTION Vector<size_t, numElementDOFs> getFEMVectorDOFIndices(
-                                    indices_t elementIndex, NDIndex<Dim> ldom) const;
+            indices_t elementIndex, NDIndex<Dim> ldom) const;
 
-        
         /**
          * @brief Get the cartesion position of a local DOF in the reference
          * element.
-         * 
+         *
          * Given the local DOF index this function will return the cartesian
          * position of this DOF with respect to the reference element.
-         * 
+         *
          */
         KOKKOS_FUNCTION point_t getLocalDOFPosition(size_t localDOFIndex) const;
 
@@ -221,7 +224,6 @@ namespace ippl {
         /// Basis functions and gradients /////////////////////////////////////
         ///////////////////////////////////////////////////////////////////////
 
-
         /**
          * @brief Evaluate the shape function of a local degree of freedom at a
          * given point in the reference element
@@ -233,8 +235,7 @@ namespace ippl {
          * @return T - The value of the shape function at the given point
          */
         KOKKOS_FUNCTION point_t evaluateRefElementShapeFunction(const size_t& localDOF,
-                                    const point_t& localPoint) const;
-        
+                                                                const point_t& localPoint) const;
 
         /**
          * @brief Evaluate the curl of the shape function of a local degree of
@@ -247,9 +248,9 @@ namespace ippl {
          * @return point_t (Vector<T, Dim>) - The curl of the shape function
          * at the given point
          */
-        KOKKOS_FUNCTION point_t evaluateRefElementShapeFunctionCurl(const size_t& localDOF,
-            const point_t& localPoint) const;
-        
+        KOKKOS_FUNCTION point_t evaluateRefElementShapeFunctionCurl(
+            const size_t& localDOF, const point_t& localPoint) const;
+
         ///////////////////////////////////////////////////////////////////////
         /// Assembly operations ///////////////////////////////////////////////
         ///////////////////////////////////////////////////////////////////////
@@ -258,6 +259,7 @@ namespace ippl {
          * @brief Assemble the left stiffness matrix A of the system Ax = b
          *
          * @param x The vector which we want to multiply
+         * @param evalFunction The element-wise operator evaluation functor
          *
          * @return The vector containing A*x
          */
@@ -267,7 +269,7 @@ namespace ippl {
         /**
          * @brief Assemble the load vector b of the system Ax = b, given a field
          * of the right hand side defined at the Nédélec DOF positions. If a
-         * functor instead of a field should be used, use the function 
+         * functor instead of a field should be used, use the function
          * \c NedelecSpace::evaluateLoadVectorFunctor.
          *
          * @param f The source field defined at the Nédélec degrees fo freedom.
@@ -275,7 +277,6 @@ namespace ippl {
          * @return The resulting rhs b of the Galerkin discretization.
          */
         FEMVector<T> evaluateLoadVector(const FEMVector<point_t>& f) const;
-        
 
         /**
          * @brief Assemble the load vector b of the system Ax = b, given a
@@ -284,41 +285,38 @@ namespace ippl {
          *
          * @param f The source function, which can be evaluated at arbitrary
          * points.
-         * 
+         *
          * @tparam F The functor type.
-         * 
+         *
          * @return The resulting rhs b of the Galerkin discretization.
          */
         template <typename F>
         FEMVector<T> evaluateLoadVectorFunctor(const F& f) const;
 
-
-
         ///////////////////////////////////////////////////////////////////////
         /// FEMVector conversion and creation//////////////////////////////////
         ///////////////////////////////////////////////////////////////////////
-        
+
         /**
          * @brief Creates and empty FEMVector.
-         * 
+         *
          * Creates and empty FEMVector which corresponds to the domain this MPI
          * rank owns (according to the ippl layout created for this mesh). To
          * this extend it will also setup all the information needed to exchange
          * halo cells.
-         * 
+         *
          * @returns An empty FEMVector for this domain.
          */
         FEMVector<T> createFEMVector() const;
-        
 
         /**
          * @brief Reconstructs function values at arbitrary points in the mesh
          * given the Nedelec DOF coefficients.
-         * 
+         *
          * This function can be used to retrieve the values of a solution
          * function at arbitrary points inside of the mesh given the Nedelec
          * DOF coefficients which solved the problem using FEM.
-         * 
+         *
          * @note Currently the function is able to handle both cases, where we
          * have that \p positions only contains positions which are inside of
          * local domain of this MPI rank (i.e. each rank gets its own unique
@@ -327,23 +325,22 @@ namespace ippl {
          * it can be guaranteed, that each rank will get its own \p positions
          * then certain parts of the function implementation can be removed.
          * Instructions for this are given in the implementation itself.
-         * 
+         *
          * @param positions The points at which the function should be
-         * evaluated. A \c Kokkos::View which stores in each element a 2D/3D 
+         * evaluated. A \c Kokkos::View which stores in each element a 2D/3D
          * point.
          * @param coef The basis function coefficients obtained via FEM.
-         * 
+         *
          * @return The function evaluated at the given points, stored inside of
          * \c Kokkos::View where each element corresponts to the function value
          * at the point described by the same element inside of \p positions.
          */
         Kokkos::View<point_t*> reconstructToPoints(const Kokkos::View<point_t*>& positions,
-            const FEMVector<T>& coef) const;
+                                                   const FEMVector<T>& coef) const;
 
         ///////////////////////////////////////////////////////////////////////
         /// Error norm computations ///////////////////////////////////////////
         ///////////////////////////////////////////////////////////////////////
-        
 
         /**
          * @brief Given the Nedelec space DoF coefficients and an analytical
@@ -354,8 +351,8 @@ namespace ippl {
          *
          * @return error - The error ||u_h - u_sol||_L2
          */
-        template <typename F> T computeError(const FEMVector<T>& u_h, const F& u_sol) const;
-
+        template <typename F>
+        T computeError(const FEMVector<T>& u_h, const F& u_sol) const;
 
         /**
          * @brief Check if a DOF is on the boundary of the mesh
@@ -363,7 +360,7 @@ namespace ippl {
          * This function takes as input the global index of a DoF and returns
          * if this DoF is on. If one would like to know which boundary this is
          * the function \c NedelecSpace::getBoundarySide can be used.
-         * 
+         *
          * @param dofIdx The global DoF index for which should be checked if it
          * is on the boundary.
          *
@@ -372,30 +369,28 @@ namespace ippl {
          */
         KOKKOS_FUNCTION bool isDOFOnBoundary(const size_t& dofIdx) const;
 
-        /** 
-        * @brief Returns which side the boundary is on.
-        * 
-        * This function takes as input the global index of a DoF and then
-        * returns on which side of the domain boundary it is on, in 2d that
-        * would be either south, north, west, east and in 3d space and ground is
-        * added. The mapping is as follows:
-        * 0 = south
-        * 1 = west
-        * 2 = north
-        * 3 = east
-        * 4 = ground
-        * 5 = space
-        * -1 = not on a boundary.
-        * 
-        * @param dofIdx the global DoF index for which the boundary side should
-        * be retrieved.
-        * 
-        * @returns Which boundary side the DoF is on or -1 if on no boundary.
-        */
+        /**
+         * @brief Returns which side the boundary is on.
+         *
+         * This function takes as input the global index of a DoF and then
+         * returns on which side of the domain boundary it is on, in 2d that
+         * would be either south, north, west, east and in 3d space and ground is
+         * added. The mapping is as follows:
+         * 0 = south
+         * 1 = west
+         * 2 = north
+         * 3 = east
+         * 4 = ground
+         * 5 = space
+         * -1 = not on a boundary.
+         *
+         * @param dofIdx the global DoF index for which the boundary side should
+         * be retrieved.
+         *
+         * @returns Which boundary side the DoF is on or -1 if on no boundary.
+         */
         KOKKOS_FUNCTION int getBoundarySide(const size_t& dofIdx) const;
 
-
-    
     private:
         /**
          * @brief Implementation of the \c NedelecSpace::createFEMVector
@@ -408,7 +403,7 @@ namespace ippl {
          * function for 3d.
          */
         FEMVector<T> createFEMVector3d() const;
-        
+
         /**
          * @brief Stores which elements (squares or cubes) belong to the current
          * MPI rank.
@@ -416,23 +411,22 @@ namespace ippl {
         Kokkos::View<size_t*> elementIndices;
 
         /**
-         * @brief Stores the positions of the local Degrees of Freedoms on the 
+         * @brief Stores the positions of the local Degrees of Freedoms on the
          * reference elements.
-         * 
+         *
          * We are saying that the local degree of freedom positions are simply
-         * the centers of the edges. 
+         * the centers of the edges.
          */
         Vector<point_t, 12> localDofPositions_m;
-        
+
         /**
          * @brief The layout of the MPI ranks over the mesh.
-         * 
+         *
          * Standart ippl layout which dictates how the MPI ranks are layed out
          * over the mesh. It is used in order to be able to create FEMVectors,
          * retreive correct DOF indices and intitalize the elementIndices.
          */
         Layout_t layout_m;
-
     };
 
 }  // namespace ippl
diff --git a/src/MaxwellSolvers/FDTDSolverBase.hpp b/src/MaxwellSolvers/FDTDSolverBase.hpp
index f929e2159..cafa81363 100644
--- a/src/MaxwellSolvers/FDTDSolverBase.hpp
+++ b/src/MaxwellSolvers/FDTDSolverBase.hpp
@@ -5,15 +5,6 @@
 
 namespace ippl {
 
-    /**
-     * @brief Constructor for the FDTDSolverBase class.
-     *
-     * Initializes the solver by setting the source and electromagnetic fields.
-     *
-     * @param source Reference to the source field.
-     * @param E Reference to the electric field.
-     * @param B Reference to the magnetic field.
-     */
     template <typename EMField, typename SourceField, fdtd_bc boundary_conditions>
     FDTDSolverBase<EMField, SourceField, boundary_conditions>::FDTDSolverBase(SourceField& source,
                                                                               EMField& E,
@@ -138,4 +129,4 @@ namespace ippl {
         Kokkos::fence();
     }
 
-}  // namespace ippl
\ No newline at end of file
+}  // namespace ippl
diff --git a/src/MaxwellSolvers/FEMMaxwellDiffusionSolver.h b/src/MaxwellSolvers/FEMMaxwellDiffusionSolver.h
index 6c627a3ff..0319ce8c4 100644
--- a/src/MaxwellSolvers/FEMMaxwellDiffusionSolver.h
+++ b/src/MaxwellSolvers/FEMMaxwellDiffusionSolver.h
@@ -1,24 +1,25 @@
 // Class FEMMaxwellDifussionSolver
-//   Solves the electric diffusion probelm given by curl(curl(E)) + E = f in 
+//   Solves the electric diffusion probelm given by curl(curl(E)) + E = f in
 //   the domain and n x E = 0 on the boundary.
 
 #ifndef IPPL_FEM_MAXWELL_DIFFUSION_SOLVER_H
 #define IPPL_FEM_MAXWELL_DIFFUSION_SOLVER_H
 
-#include "LinearSolvers/PCG.h"
-#include "Maxwell.h"
+#include <fstream>
 #include <iomanip>
 #include <iostream>
-#include <fstream>
+
+#include "LinearSolvers/PCG.h"
+#include "Maxwell.h"
 
 namespace ippl {
 
     /**
      * @brief Representation of the lhs of the problem we are trying to solve.
-     * 
-     * In our case this corresponds to the variational formulation of the 
+     *
+     * In our case this corresponds to the variational formulation of the
      * curl(curl(E)) + E and is curl(b_i)*curl(b_j) + b_i*b_j.
-     * 
+     *
      * @tparam T The type we are working with.
      * @tparam Dim the dimension of the space.
      * @tparam numElementDOFs the number of DOFs per element that we have.
@@ -27,7 +28,7 @@ namespace ippl {
     struct EvalFunctor {
         /**
          * @brief The inverse transpose Jacobian.
-         * 
+         *
          * As we have a unirectangular grid it is the same for all the differnt
          * Elements and we therefore have to store it only once.
          */
@@ -35,7 +36,7 @@ namespace ippl {
 
         /**
          * @brief The determinant of the Jacobian.
-         * 
+         *
          * As we have a unirectangular grid it is the same for all the differnt
          * Elements and we therefore have to store it only once.
          */
@@ -51,28 +52,28 @@ namespace ippl {
         /**
          * @brief Returns the evaluation of
          * (curl(b_i)*curl(b_j) + b_i*b_j)*absDetDPhi.
-         * 
+         *
          * This function takes as input the basis function values and their curl
          * for the different DOFs and returns the evaluation of the inner part
          * of the integral of the variational formuation, which corresponds to
          * (curl(b_i)*curl(b_j) + b_i*b_j), but note that we additionally also
          * multiply this with absDetDPhi, which is required by the quadrature
          * rule. In theroy this could also be done outside of this.
-         * 
+         *
          * @param i The first DOF index.
          * @param j The second DOF index.
          * @param curl_b_q_k The curl of the DOFs.
          * @param val_b_q_k The values of the DOFs.
-         * 
+         *
          * @returns (curl(b_i)*curl(b_j) + b_i*b_j)*absDetDPhi
          */
-        KOKKOS_FUNCTION auto operator()(size_t i, size_t j,
+        KOKKOS_FUNCTION auto operator()(
+            size_t i, size_t j,
             const ippl::Vector<ippl::Vector<T, Dim>, numElementDOFs>& curl_b_q_k,
             const ippl::Vector<ippl::Vector<T, Dim>, numElementDOFs>& val_b_q_k) const {
-            
-            T curlTerm = dot(DPhiInvT*curl_b_q_k[j], DPhiInvT*curl_b_q_k[i]).apply();
+            T curlTerm = dot(DPhiInvT * curl_b_q_k[j], DPhiInvT * curl_b_q_k[i]).apply();
             T massTerm = dot(val_b_q_k[j], val_b_q_k[i]).apply();
-            return (curlTerm + massTerm)*absDetDPhi;
+            return (curlTerm + massTerm) * absDetDPhi;
         }
     };
 
@@ -92,41 +93,41 @@ namespace ippl {
         // vector data represented by an ippl::Vector
         using T = typename FieldType::value_type::value_type;
 
-        typedef Vector<T,Dim> point_t;
+        typedef Vector<T, Dim> point_t;
 
     public:
-        using Base = Maxwell<FieldType, FieldType>;
+        using Base     = Maxwell<FieldType, FieldType>;
         using MeshType = typename FieldType::Mesh_t;
 
         // PCG (Preconditioned Conjugate Gradient) is the solver algorithm used
         using PCGSolverAlgorithm_t = CG<FEMVector<T>, FEMVector<T>, FEMVector<T>, FEMVector<T>,
-            FEMVector<T>, FEMVector<T>, FEMVector<T>>;
+                                        FEMVector<T>, FEMVector<T>, FEMVector<T>>;
 
         // FEM Space types
-        using ElementType = std::conditional_t<Dim == 2, ippl::QuadrilateralElement<T>, 
-                                ippl::HexahedralElement<T>>;
+        using ElementType =
+            std::conditional_t<Dim == 2, ippl::QuadrilateralElement<T>, ippl::HexahedralElement<T>>;
 
         using QuadratureType = GaussJacobiQuadrature<T, 5, ElementType>;
 
         using NedelecType = NedelecSpace<T, Dim, 1, ElementType, QuadratureType, FieldType>;
 
         // default constructor (compatibility with Alpine)
-        FEMMaxwellDiffusionSolver() 
+        FEMMaxwellDiffusionSolver()
             : Base()
             , rhsVector_m(nullptr)
             , refElement_m()
             , quadrature_m(refElement_m, 0.0, 0.0)
-            , nedelecSpace_m(*(new MeshType(NDIndex<Dim>(Vector<unsigned, Dim>(0)), Vector<T, Dim>(0),
-                                Vector<T, Dim>(0))), refElement_m, quadrature_m)
-        {}
+            , nedelecSpace_m(*(new MeshType(NDIndex<Dim>(Vector<unsigned, Dim>(0)),
+                                            Vector<T, Dim>(0), Vector<T, Dim>(0))),
+                             refElement_m, quadrature_m) {}
 
-        FEMMaxwellDiffusionSolver(FieldType& lhs, FieldType& rhs, const FEMVector<point_t>& rhsVector)
+        FEMMaxwellDiffusionSolver(FieldType& lhs, FieldType& rhs,
+                                  const FEMVector<point_t>& rhsVector)
             : Base(lhs, lhs, rhs)
             , rhsVector_m(nullptr)
             , refElement_m()
             , quadrature_m(refElement_m, 0.0, 0.0)
             , nedelecSpace_m(rhs.get_mesh(), refElement_m, quadrature_m, rhs.getLayout()) {
-            
             static_assert(std::is_floating_point<T>::value, "Not a floating point type");
             setDefaultParameters();
 
@@ -136,11 +137,9 @@ namespace ippl {
 
             rhsVector_m->accumulateHalo();
             rhsVector_m->fillHalo();
-            
         }
 
         void setRhs(FieldType& rhs, const FEMVector<point_t>& rhsVector) {
-            
             Base::setRhs(rhs);
 
             // Calcualte the rhs, using the Nedelec space
@@ -155,7 +154,6 @@ namespace ippl {
          * @brief Solve the equation using finite element methods.
          */
         void solve() override {
-
             const Vector<size_t, Dim> zeroNdIndex = Vector<size_t, Dim>(0);
 
             // We can pass the zeroNdIndex here, since the transformation
@@ -174,16 +172,16 @@ namespace ippl {
 
             // Create the functor object which stores the function we have to
             // solve for the lhs
-            EvalFunctor<T, Dim, NedelecType::numElementDOFs> maxwellDiffusionEval(
-                DPhiInvT, absDetDPhi);
-            
-            // The Ax operator
-            const auto algoOperator = [maxwellDiffusionEval, this](FEMVector<T> vector)
-                                                                            -> FEMVector<T> {
+            EvalFunctor<T, Dim, NedelecType::numElementDOFs> maxwellDiffusionEval(DPhiInvT,
+                                                                                  absDetDPhi);
 
+            // The Ax operator
+            const auto algoOperator = [maxwellDiffusionEval,
+                                       this](FEMVector<T> vector) -> FEMVector<T> {
                 vector.fillHalo();
 
-                FEMVector<T> return_vector = nedelecSpace_m.evaluateAx(vector,maxwellDiffusionEval);
+                FEMVector<T> return_vector =
+                    nedelecSpace_m.evaluateAx(vector, maxwellDiffusionEval);
 
                 return_vector.accumulateHalo();
 
@@ -192,10 +190,10 @@ namespace ippl {
 
             // setup the CG solver
             pcg_algo_m.setOperator(algoOperator);
-            
+
             // Create the coefficient vector for the solution
             FEMVector<T> lhsVector = rhsVector_m->deepCopy();
-            
+
             // Solve the system using CG
             try {
                 pcg_algo_m(lhsVector, *rhsVector_m, this->params_m);
@@ -203,13 +201,12 @@ namespace ippl {
                 std::string msg = e.where() + ": " + e.what() + "\n";
                 Kokkos::abort(msg.c_str());
             }
-            
+
             // store solution.
             lhsVector_m = std::make_unique<FEMVector<T>>(lhsVector);
 
             // set the boundary values to the correct values.
             lhsVector.fillHalo();
-
         }
 
         /**
@@ -227,10 +224,10 @@ namespace ippl {
 
         /**
          * @brief Reconstructs function values at arbitrary points in the mesh.
-         * 
+         *
          * This function can be used to retrieve the values of a solution
          * function at arbitrary points inside of the mesh.
-         * 
+         *
          * @note Currently the function is able to handle both cases, where we
          * have that \p positions only contains positions which are inside of
          * local domain of this MPI rank (i.e. each rank gets its own unique
@@ -239,11 +236,11 @@ namespace ippl {
          * it can be guaranteed, that each rank will get its own \p positions
          * then certain parts of the function implementation can be removed.
          * Instructions for this are given in the implementation itself.
-         * 
+         *
          * @param positions The points at which the function should be
-         * evaluated. A \c Kokkos::View which stores in each element a 2D/3D 
+         * evaluated. A \c Kokkos::View which stores in each element a 2D/3D
          * point.
-         * 
+         *
          * @return The function evaluated at the given points, stored inside of
          * \c Kokkos::View where each element corresponts to the function value
          * at the point described by the same element inside of \p positions.
@@ -252,12 +249,10 @@ namespace ippl {
             return this->nedelecSpace_m.reconstructToPoints(positions, *lhsVector_m);
         }
 
-
-        
         /**
          * @brief Given an analytical solution computes the L2 norm error.
          *
-         * @param analytical The analytical solution (functor)
+         * @param analytic The analytical solution functor
          *
          * @return error - The error ||u - analytical||_L2
          */
@@ -267,14 +262,12 @@ namespace ippl {
             return error_norm;
         }
 
-
     protected:
-
         /**
          * @brief The CG Solver we use
          */
         PCGSolverAlgorithm_t pcg_algo_m;
-        
+
         /**
          * @brief Sets the default values for the CG solver.
          * Defaults are: max Iterations = 10, tolerance = 1e-13
@@ -310,7 +303,7 @@ namespace ippl {
 
         /**
          * @brief The Nedelec Space object.
-         * 
+         *
          * This is the representation of the Nedelec space that we have and
          * which we use to interact with all the Nedelec stuff.
          */
@@ -319,6 +312,4 @@ namespace ippl {
 
 }  // namespace ippl
 
-
-
-#endif // IPPL_FEM_MAXWELL_DIFFUSION_SOLVER_H
+#endif  // IPPL_FEM_MAXWELL_DIFFUSION_SOLVER_H
diff --git a/src/MaxwellSolvers/NonStandardFDTDSolver.hpp b/src/MaxwellSolvers/NonStandardFDTDSolver.hpp
index a5b565ef9..5bd5fee20 100644
--- a/src/MaxwellSolvers/NonStandardFDTDSolver.hpp
+++ b/src/MaxwellSolvers/NonStandardFDTDSolver.hpp
@@ -5,17 +5,6 @@
 
 namespace ippl {
 
-    /**
-     * @brief Constructor for the NonStandardFDTDSolver class.
-     *
-     * This constructor initializes the NonStandardFDTDSolver with the given source field and
-     * electromagnetic fields. It checks the dispersion-free CFL condition and initializes the
-     * solver.
-     *
-     * @param source The source field.
-     * @param E The electric field.
-     * @param B The magnetic field.
-     */
     template <typename EMField, typename SourceField, fdtd_bc boundary_conditions>
     NonStandardFDTDSolver<EMField, SourceField, boundary_conditions>::NonStandardFDTDSolver(
         SourceField& source, EMField& E, EMField& B)
@@ -141,4 +130,4 @@ namespace ippl {
             this->setPeriodicBoundaryConditions();
         }
     };
-}  // namespace ippl
\ No newline at end of file
+}  // namespace ippl
diff --git a/src/MaxwellSolvers/StandardFDTDSolver.hpp b/src/MaxwellSolvers/StandardFDTDSolver.hpp
index beefb3660..e61250245 100644
--- a/src/MaxwellSolvers/StandardFDTDSolver.hpp
+++ b/src/MaxwellSolvers/StandardFDTDSolver.hpp
@@ -5,16 +5,6 @@
 
 namespace ippl {
 
-    /**
-     * @brief Constructor for the StandardFDTDSolver class.
-     *
-     * This constructor initializes the StandardFDTDSolver with the given source field and
-     * electromagnetic fields and initializes the solver.
-     *
-     * @param source The source field.
-     * @param E The electric field.
-     * @param B The magnetic field.
-     */
     template <typename EMField, typename SourceField, fdtd_bc boundary_conditions>
     StandardFDTDSolver<EMField, SourceField, boundary_conditions>::StandardFDTDSolver(
         SourceField& source, EMField& E, EMField& B)
@@ -122,4 +112,4 @@ namespace ippl {
             this->setPeriodicBoundaryConditions();
         }
     }
-}  // namespace ippl
\ No newline at end of file
+}  // namespace ippl
diff --git a/src/Particle/ParticleBase.h b/src/Particle/ParticleBase.h
index afb8e32d2..9fbed2faa 100644
--- a/src/Particle/ParticleBase.h
+++ b/src/Particle/ParticleBase.h
@@ -68,7 +68,7 @@ namespace ippl {
      *
      *  Minimal empty base class for all ParticleBase specializations.
      *  Needed for e.g: c++20 constraints and concepts using std::derived_from
-     * 
+     *
      */
     class ParticleBaseBase {
     public:
@@ -84,7 +84,7 @@ namespace ippl {
      * IDs will be disabled for the bunch)
      */
     template <class PLayout, typename... IDProperties>
-    class ParticleBase: public ParticleBaseBase {
+    class ParticleBase : public ParticleBaseBase {
         constexpr static bool EnableIDs = sizeof...(IDProperties) > 0;
 
     public:
@@ -173,7 +173,7 @@ namespace ippl {
 
         /*!
          * Set all boundary conditions
-         * @param bc the boundary conditions
+         * @param bcs the boundary conditions
          */
         void setParticleBC(const bc_container_type& bcs) { layout_m->setParticleBC(bcs); }
 
@@ -308,7 +308,6 @@ namespace ippl {
          * @tparam HashType the hash view type
          * @param rank the destination rank
          * @param tag the MPI tag
-         * @param sendNum the number of messages already sent (to distinguish the buffers)
          * @param requests destination vector in which to store the MPI requests for polling
          * purposes
          * @param hash a hash view indicating which particles need to be sent to which rank
@@ -321,7 +320,6 @@ namespace ippl {
          * Receives particles from another rank
          * @param rank the source rank
          * @param tag the MPI tag
-         * @param recvNum the number of messages already received (to distinguish the buffers)
          * @param nRecvs the number of particles to receive
          */
         void recvFromRank(int rank, int tag, size_type nRecvs);
@@ -329,6 +327,7 @@ namespace ippl {
         /*!
          * Serialize to do MPI calls.
          * @param ar archive
+         * @param nsends number of particles to serialize
          */
         template <typename Archive>
         void serialize(Archive& ar, size_type nsends);
@@ -336,6 +335,7 @@ namespace ippl {
         /*!
          * Deserialize to do MPI calls.
          * @param ar archive
+         * @param nrecvs number of particles to deserialize
          */
         template <typename Archive>
         void deserialize(Archive& ar, size_type nrecvs);
@@ -352,14 +352,13 @@ namespace ippl {
     protected:
         /*!
          * Fill attributes of buffer.
-         * @param buffer to send
          * @param hash function to access index.
          */
         void pack(const hash_container_type& hash);
 
         /*!
          * Fill my attributes.
-         * @param buffer received
+         * @param nrecvs number of particles to unpack
          */
         void unpack(size_type nrecvs);
 
diff --git a/src/Particle/ParticleLayout.h b/src/Particle/ParticleLayout.h
index 15914322a..4d2cb60fd 100644
--- a/src/Particle/ParticleLayout.h
+++ b/src/Particle/ParticleLayout.h
@@ -73,22 +73,20 @@ namespace ippl {
             void setParticleBC(bc_container_type bcs) { bcs_m = bcs; }
 
             /*!
-             * Copy over the given boundary conditions.
-             * @param bcs are the boundary conditions
+             * Return the current boundary conditions.
              */
             const bc_container_type& getParticleBC() const { return bcs_m; }
 
             /*!
              * Use the same boundary condition on each face
-             * @param bcs are the boundary conditions
+             * @param bc is the boundary condition to apply on all faces
              */
             void setParticleBC(BC bc) { bcs_m.fill(bc); }
 
             /*!
              * Apply the given boundary conditions to the current particle positions.
-             * @tparam R is the particle position attribute
-             * @tparam nr is the NDRegion
-             * @param
+             * @param R is the particle position attribute
+             * @param nr is the particle domain
              */
             void applyBC(const particle_position_type& R, const NDRegion<T, Dim>& nr);
 
diff --git a/src/Particle/ParticleSpatialLayout.h b/src/Particle/ParticleSpatialLayout.h
index 9f82c870b..568b30753 100644
--- a/src/Particle/ParticleSpatialLayout.h
+++ b/src/Particle/ParticleSpatialLayout.h
@@ -23,16 +23,16 @@
 #ifndef IPPL_PARTICLE_SPATIAL_LAYOUT_H
 #define IPPL_PARTICLE_SPATIAL_LAYOUT_H
 
+#include <vector>
+
 #include "Types/IpplTypes.h"
 
+#include "Communicate/Window.h"
 #include "FieldLayout/FieldLayout.h"
 #include "Particle/ParticleBase.h"
 #include "Particle/ParticleLayout.h"
 #include "Region/RegionLayout.h"
 
-#include "Communicate/Window.h"
-#include <vector>
-
 namespace ippl {
 
     /*!
@@ -81,14 +81,13 @@ namespace ippl {
 
         //! The FieldLayout containing information on nearest neighbors
         FieldLayout_t& flayout_m;
-        
-        
+
         // Vector keeping track of the recieves from all ranks
         std::vector<size_type> nRecvs_m;
-  
+
         // MPI RMA window for one-sided communication
         mpi::rma::Window<mpi::rma::Active> window_m;
-        
+
         //! Type of the Kokkos view containing the local regions.
         using region_view_type = typename RegionLayout_t::view_type;
         //! Type of a single Region object.
@@ -117,11 +116,15 @@ namespace ippl {
          * @param ranks the integer view in which to store the destination ranks
          * @param invalid the boolean view in which to store whether each particle
          * is invalidated, i.e. needs to be sent to another rank
+         * @param nSends_dview the view storing per-rank send counts
+         * @param sends_dview the view storing per-rank send offsets
          * @return The total number of invalidated particles
          */
         template <typename ParticleContainer>
-        std::pair<size_type,size_type> locateParticles(const ParticleContainer& pc, locate_type& ranks,
-                                  bool_type& invalid, locate_type& nSends_dview, locate_type& sends_dview) const;
+        std::pair<size_type, size_type> locateParticles(const ParticleContainer& pc,
+                                                        locate_type& ranks, bool_type& invalid,
+                                                        locate_type& nSends_dview,
+                                                        locate_type& sends_dview) const;
 
         /*!
          * @param rank we sent to
@@ -141,4 +144,3 @@ namespace ippl {
 #include "Particle/ParticleSpatialLayout.hpp"
 
 #endif
-
diff --git a/src/Particle/ParticleSpatialOverlapLayout.h b/src/Particle/ParticleSpatialOverlapLayout.h
index cae305e6d..34e341fa2 100644
--- a/src/Particle/ParticleSpatialOverlapLayout.h
+++ b/src/Particle/ParticleSpatialOverlapLayout.h
@@ -315,7 +315,7 @@ namespace ippl {
             hash_type cellPermutationForward);
 
         /*!
-         * @brieg compute the nd-cell-index from a flattened (non-permuted) index
+         * @brief compute the nd-cell-index from a flattened (non-permuted) index
          * @param nonPermutedIndex the flat index to transform
          * @param numCells in each dimension
          */
diff --git a/src/PoissonSolvers/FEMPoissonSolver.h b/src/PoissonSolvers/FEMPoissonSolver.h
index 2a1bdf862..8c183a6ea 100644
--- a/src/PoissonSolvers/FEMPoissonSolver.h
+++ b/src/PoissonSolvers/FEMPoissonSolver.h
@@ -5,9 +5,9 @@
 #ifndef IPPL_FEMPOISSONSOLVER_H
 #define IPPL_FEMPOISSONSOLVER_H
 
+#include "EvalFunctor.h"
 #include "LinearSolvers/PCG.h"
 #include "Poisson.h"
-#include "EvalFunctor.h"
 
 namespace ippl {
 
@@ -18,7 +18,8 @@ namespace ippl {
      * @tparam FieldLHS field type for the left hand side
      * @tparam FieldRHS field type for the right hand side
      */
-    template <typename FieldLHS, typename FieldRHS = FieldLHS, unsigned Order = 1, unsigned QuadNumNodes = 5>
+    template <typename FieldLHS, typename FieldRHS = FieldLHS, unsigned Order = 1,
+              unsigned QuadNumNodes = 5>
     class FEMPoissonSolver : public Poisson<FieldLHS, FieldRHS> {
         constexpr static unsigned Dim = FieldLHS::dim;
         using Tlhs                    = typename FieldLHS::value_type;
@@ -40,16 +41,17 @@ namespace ippl {
 
         using QuadratureType = GaussJacobiQuadrature<Tlhs, QuadNumNodes, ElementType>;
 
-        using LagrangeType = LagrangeSpace<Tlhs, Dim, Order, ElementType, QuadratureType, FieldLHS, FieldRHS>;
+        using LagrangeType =
+            LagrangeSpace<Tlhs, Dim, Order, ElementType, QuadratureType, FieldLHS, FieldRHS>;
 
         // default constructor (compatibility with Alpine)
-        FEMPoissonSolver() 
+        FEMPoissonSolver()
             : Base()
             , refElement_m()
             , quadrature_m(refElement_m, 0.0, 0.0)
-            , lagrangeSpace_m(*(new MeshType(NDIndex<Dim>(Vector<unsigned, Dim>(0)), Vector<Tlhs, Dim>(0),
-                                Vector<Tlhs, Dim>(0))), refElement_m, quadrature_m)
-        {
+            , lagrangeSpace_m(*(new MeshType(NDIndex<Dim>(Vector<unsigned, Dim>(0)),
+                                             Vector<Tlhs, Dim>(0), Vector<Tlhs, Dim>(0))),
+                              refElement_m, quadrature_m) {
             setDefaultParameters();
         }
 
@@ -57,8 +59,7 @@ namespace ippl {
             : Base(lhs, rhs)
             , refElement_m()
             , quadrature_m(refElement_m, 0.0, 0.0)
-            , lagrangeSpace_m(rhs.get_mesh(), refElement_m, quadrature_m, rhs.getLayout())
-        {
+            , lagrangeSpace_m(rhs.get_mesh(), refElement_m, quadrature_m, rhs.getLayout()) {
             static_assert(std::is_floating_point<Tlhs>::value, "Not a floating point type");
             setDefaultParameters();
             pcg_algo_m.initializeFields(rhs.get_mesh(), rhs.getLayout());
@@ -74,9 +75,7 @@ namespace ippl {
         /**
          * @brief Return the LagrangeSpace object.
          */
-        LagrangeType& getSpace() {
-            return lagrangeSpace_m;
-        }
+        LagrangeType& getSpace() { return lagrangeSpace_m; }
 
         /**
          * @brief Solve the poisson equation using finite element methods.
@@ -103,14 +102,15 @@ namespace ippl {
             const Tlhs absDetDPhi = Kokkos::abs(
                 refElement_m.getDeterminantOfTransformationJacobian(firstElementVertexPoints));
 
-            EvalFunctor<Tlhs, Dim, LagrangeType::numElementDOFs> poissonEquationEval(
-                DPhiInvT, absDetDPhi);
+            EvalFunctor<Tlhs, Dim, LagrangeType::numElementDOFs> poissonEquationEval(DPhiInvT,
+                                                                                     absDetDPhi);
 
             // get BC type of our RHS
             BConds<FieldRHS, Dim>& bcField = (this->rhs_mp)->getFieldBC();
-            FieldBC bcType = bcField[0]->getBCType();
+            FieldBC bcType                 = bcField[0]->getBCType();
 
-            const auto algoOperator = [poissonEquationEval, &bcField, this](rhs_type field) -> lhs_type {
+            const auto algoOperator = [poissonEquationEval, &bcField,
+                                       this](rhs_type field) -> lhs_type {
                 // set appropriate BCs for the field as the info gets lost in the CG iteration
                 field.setFieldBC(bcField);
 
@@ -125,8 +125,9 @@ namespace ippl {
 
             // send boundary values to RHS (load vector) i.e. lifting (Dirichlet BCs)
             if (bcType == CONSTANT_FACE) {
-                *(this->rhs_mp) = *(this->rhs_mp) -
-                    lagrangeSpace_m.evaluateAx_lift(*(this->rhs_mp), poissonEquationEval);
+                *(this->rhs_mp) =
+                    *(this->rhs_mp)
+                    - lagrangeSpace_m.evaluateAx_lift(*(this->rhs_mp), poissonEquationEval);
             }
 
             // start a timer
@@ -167,16 +168,16 @@ namespace ippl {
 
         /**
          * Query the average of the solution
-         * @param vol Boolean indicating whether we divide by volume or not
+         * @param Vol Boolean indicating whether we divide by volume or not
          * @return avg (offset for null space test cases if divided by volume)
          */
         Tlhs getAvg(bool Vol = false) {
             Tlhs avg = this->lagrangeSpace_m.computeAvg(*(this->lhs_mp));
             if (Vol) {
                 lhs_type unit((this->lhs_mp)->get_mesh(), (this->lhs_mp)->getLayout());
-                unit = 1.0;
+                unit     = 1.0;
                 Tlhs vol = this->lagrangeSpace_m.computeAvg(unit);
-                return avg/vol;
+                return avg / vol;
             } else {
                 return avg;
             }
diff --git a/src/PoissonSolvers/PreconditionedFEMPoissonSolver.h b/src/PoissonSolvers/PreconditionedFEMPoissonSolver.h
index c085478b1..64a917a60 100644
--- a/src/PoissonSolvers/PreconditionedFEMPoissonSolver.h
+++ b/src/PoissonSolvers/PreconditionedFEMPoissonSolver.h
@@ -6,10 +6,10 @@
 #define IPPL_PRECONFEMPOISSONSOLVER_H
 
 // #include "FEM/FiniteElementSpace.h"
+#include "EvalFunctor.h"
 #include "LaplaceHelpers.h"
 #include "LinearSolvers/PCG.h"
 #include "Poisson.h"
-#include "EvalFunctor.h"
 
 namespace ippl {
     /**
@@ -41,16 +41,17 @@ namespace ippl {
 
         using QuadratureType = GaussJacobiQuadrature<Tlhs, 5, ElementType>;
 
-        using LagrangeType = LagrangeSpace<Tlhs, Dim, 1, ElementType, QuadratureType, FieldLHS, FieldRHS>;
+        using LagrangeType =
+            LagrangeSpace<Tlhs, Dim, 1, ElementType, QuadratureType, FieldLHS, FieldRHS>;
 
         // default constructor (compatibility with Alpine)
-        PreconditionedFEMPoissonSolver() 
+        PreconditionedFEMPoissonSolver()
             : Base()
             , refElement_m()
             , quadrature_m(refElement_m, 0.0, 0.0)
-            , lagrangeSpace_m(*(new MeshType(NDIndex<Dim>(Vector<unsigned, Dim>(0)), Vector<Tlhs, Dim>(0),
-                                Vector<Tlhs, Dim>(0))), refElement_m, quadrature_m)
-        {
+            , lagrangeSpace_m(*(new MeshType(NDIndex<Dim>(Vector<unsigned, Dim>(0)),
+                                             Vector<Tlhs, Dim>(0), Vector<Tlhs, Dim>(0))),
+                              refElement_m, quadrature_m) {
             setDefaultParameters();
         }
 
@@ -58,8 +59,7 @@ namespace ippl {
             : Base(lhs, rhs)
             , refElement_m()
             , quadrature_m(refElement_m, 0.0, 0.0)
-            , lagrangeSpace_m(rhs.get_mesh(), refElement_m, quadrature_m, rhs.getLayout())
-        {
+            , lagrangeSpace_m(rhs.get_mesh(), refElement_m, quadrature_m, rhs.getLayout()) {
             static_assert(std::is_floating_point<Tlhs>::value, "Not a floating point type");
             setDefaultParameters();
         }
@@ -73,9 +73,7 @@ namespace ippl {
         /**
          * @brief Return the LagrangeSpace object.
          */
-        LagrangeType& getSpace() {
-            return lagrangeSpace_m;
-        }
+        LagrangeType& getSpace() { return lagrangeSpace_m; }
 
         /**
          * @brief Solve the poisson equation using finite element methods.
@@ -102,14 +100,15 @@ namespace ippl {
             const Tlhs absDetDPhi = Kokkos::abs(
                 refElement_m.getDeterminantOfTransformationJacobian(firstElementVertexPoints));
 
-            EvalFunctor<Tlhs, Dim, LagrangeType::numElementDOFs> poissonEquationEval(
-                DPhiInvT, absDetDPhi);
+            EvalFunctor<Tlhs, Dim, LagrangeType::numElementDOFs> poissonEquationEval(DPhiInvT,
+                                                                                     absDetDPhi);
 
             // get BC type of our RHS
             BConds<FieldRHS, Dim>& bcField = (this->rhs_mp)->getFieldBC();
-            FieldBC bcType = bcField[0]->getBCType();
+            FieldBC bcType                 = bcField[0]->getBCType();
 
-            const auto algoOperator = [poissonEquationEval, &bcField, this](rhs_type field) -> lhs_type {
+            const auto algoOperator = [poissonEquationEval, &bcField,
+                                       this](rhs_type field) -> lhs_type {
                 // set appropriate BCs for the field as the info gets lost in the CG iteration
                 field.setFieldBC(bcField);
 
@@ -120,7 +119,8 @@ namespace ippl {
                 return return_field;
             };
 
-            const auto algoOperatorL = [poissonEquationEval, &bcField, this](lhs_type field) -> lhs_type {
+            const auto algoOperatorL = [poissonEquationEval, &bcField,
+                                        this](lhs_type field) -> lhs_type {
                 // set appropriate BCs for the field as the info gets lost in the CG iteration
                 field.setFieldBC(bcField);
 
@@ -131,7 +131,8 @@ namespace ippl {
                 return return_field;
             };
 
-            const auto algoOperatorU = [poissonEquationEval, &bcField, this](lhs_type field) -> lhs_type {
+            const auto algoOperatorU = [poissonEquationEval, &bcField,
+                                        this](lhs_type field) -> lhs_type {
                 // set appropriate BCs for the field as the info gets lost in the CG iteration
                 field.setFieldBC(bcField);
 
@@ -142,29 +143,34 @@ namespace ippl {
                 return return_field;
             };
 
-            const auto algoOperatorUL = [poissonEquationEval, &bcField, this](lhs_type field) -> lhs_type {
+            const auto algoOperatorUL = [poissonEquationEval, &bcField,
+                                         this](lhs_type field) -> lhs_type {
                 // set appropriate BCs for the field as the info gets lost in the CG iteration
                 field.setFieldBC(bcField);
 
                 field.fillHalo();
 
-                auto return_field = lagrangeSpace_m.evaluateAx_upperlower(field, poissonEquationEval);
+                auto return_field =
+                    lagrangeSpace_m.evaluateAx_upperlower(field, poissonEquationEval);
 
                 return return_field;
             };
 
-            const auto algoOperatorInvD = [poissonEquationEval, &bcField, this](lhs_type field) -> lhs_type {
+            const auto algoOperatorInvD = [poissonEquationEval, &bcField,
+                                           this](lhs_type field) -> lhs_type {
                 // set appropriate BCs for the field as the info gets lost in the CG iteration
                 field.setFieldBC(bcField);
 
                 field.fillHalo();
 
-                auto return_field = lagrangeSpace_m.evaluateAx_inversediag(field, poissonEquationEval);
+                auto return_field =
+                    lagrangeSpace_m.evaluateAx_inversediag(field, poissonEquationEval);
 
                 return return_field;
             };
 
-            const auto algoOperatorD = [poissonEquationEval, &bcField, this](lhs_type field) -> lhs_type {
+            const auto algoOperatorD = [poissonEquationEval, &bcField,
+                                        this](lhs_type field) -> lhs_type {
                 // set appropriate BCs for the field as the info gets lost in the CG iteration
                 field.setFieldBC(bcField);
 
@@ -183,19 +189,19 @@ namespace ippl {
             int inner    = this->params_m.template get<int>("gauss_seidel_inner_iterations");
             int outer    = this->params_m.template get<int>("gauss_seidel_outer_iterations");
             double omega = this->params_m.template get<double>("ssor_omega");
-            int richardson_iterations =
-                this->params_m.template get<int>("richardson_iterations");
+            int richardson_iterations = this->params_m.template get<int>("richardson_iterations");
 
             pcg_algo_m.setPreconditioner(algoOperator, algoOperatorL, algoOperatorU, algoOperatorUL,
-                                     algoOperatorInvD, algoOperatorD, 0, 0, preconditioner_type,
-                                     level, degree, richardson_iterations, inner, outer, omega);
+                                         algoOperatorInvD, algoOperatorD, 0, 0, preconditioner_type,
+                                         level, degree, richardson_iterations, inner, outer, omega);
 
             pcg_algo_m.setOperator(algoOperator);
 
             // send boundary values to RHS (load vector) i.e. lifting (Dirichlet BCs)
             if (bcType == CONSTANT_FACE) {
-                *(this->rhs_mp) = *(this->rhs_mp) -
-                    lagrangeSpace_m.evaluateAx_lift(*(this->rhs_mp), poissonEquationEval);
+                *(this->rhs_mp) =
+                    *(this->rhs_mp)
+                    - lagrangeSpace_m.evaluateAx_lift(*(this->rhs_mp), poissonEquationEval);
             }
 
             // start a timer
@@ -240,16 +246,16 @@ namespace ippl {
 
         /**
          * Query the average of the solution
-         * @param vol Boolean indicating whether we divide by volume or not
+         * @param Vol Boolean indicating whether we divide by volume or not
          * @return avg (offset for null space test cases if divided by volume)
          */
         Tlhs getAvg(bool Vol = false) {
             Tlhs avg = this->lagrangeSpace_m.computeAvg(*(this->lhs_mp));
             if (Vol) {
                 lhs_type unit((this->lhs_mp)->get_mesh(), (this->lhs_mp)->getLayout());
-                unit = 1.0;
+                unit     = 1.0;
                 Tlhs vol = this->lagrangeSpace_m.computeAvg(unit);
-                return avg/vol;
+                return avg / vol;
             } else {
                 return avg;
             }
diff --git a/src/Random/NormalDistribution.h b/src/Random/NormalDistribution.h
index 763afbe54..5dedecca6 100644
--- a/src/Random/NormalDistribution.h
+++ b/src/Random/NormalDistribution.h
@@ -43,7 +43,7 @@ namespace ippl {
          * @brief An estimator for the initial guess that is used in Newton-Raphson method of
          * Inverste Transfrom Sampling
          *
-         * @param x Input value.
+         * @param u Input value.
          * @param mean Mean of the distribution.
          * @param stddev Standard deviation of the distribution.
          * @return The estimate value.
@@ -156,7 +156,7 @@ namespace ippl {
              */
             KOKKOS_INLINE_FUNCTION NormalDistribution(const T* par_p)
                 : ippl::random::Distribution<T, Dim, 2 * Dim, NormalDistributionFunctions<T>>(
-                    par_p) {}
+                      par_p) {}
         };
 
     }  // namespace random