FDTD solver

cmake/Dependencies.cmake

@@ -110,3 +110,12 @@ if(IPPL_ENABLE_UNIT_TESTS)
     FetchContent_MakeAvailable(googletest)
     message(STATUS "✅ GoogleTest loaded for unit tests.")
 endif()
+
+if(IPPL_ENABLE_TESTS)
+    set(DOWNLOADED_HEADERS_DIR "${CMAKE_CURRENT_BINARY_DIR}/downloaded_headers")
+    file(DOWNLOAD
+        https://raw.githubusercontent.com/manuel5975p/stb/master/stb_image_write.h
+        "${DOWNLOADED_HEADERS_DIR}/stb_image_write.h"
+    )
+    message(STATUS "✅ stb_image_write loaded for testing FDTD solver.")
+endif()

src/MaxwellSolvers/CMakeLists.txt

@@ -13,6 +13,12 @@ target_include_directories(ippl
 # Install MaxwellSolvers header
 install(FILES
     Maxwell.h
+    FDTDSolverBase.h
+    FDTDSolverBase.hpp
+    StandardFDTDSolver.h
+    StandardFDTDSolver.hpp
+    NonStandardFDTDSolver.h
+    NonStandardFDTDSolver.hpp
     DESTINATION include/MaxwellSolvers
 )
 

src/MaxwellSolvers/FDTDSolverBase.h

@@ -0,0 +1,101 @@
+#ifndef IPPL_FDTD_H
+#define IPPL_FDTD_H
+#include <cstddef>
+using std::size_t;
+#include "Types/Vector.h"
+
+#include "FieldLayout/FieldLayout.h"
+#include "MaxwellSolvers/AbsorbingBC.h"
+#include "MaxwellSolvers/Maxwell.h"
+#include "Meshes/UniformCartesian.h"
+#include "Particle/ParticleBase.h"
+
+namespace ippl {
+    enum fdtd_bc {
+        periodic,
+        absorbing
+    };
+
+    /**
+     * @class FDTDSolverBase
+     * @brief Base class for FDTD solvers in the ippl library.
+     *
+     * @tparam EMField The type representing the electromagnetic field.
+     * @tparam SourceField The type representing the source field.
+     * @tparam boundary_conditions The boundary conditions to be applied (default is periodic).
+     */
+    template <typename EMField, typename SourceField, fdtd_bc boundary_conditions = periodic>
+    class FDTDSolverBase : public Maxwell<EMField, SourceField> {
+    public:
+        constexpr static unsigned Dim = EMField::dim;
+        using scalar                  = typename EMField::value_type::value_type;
+        using Vector_t                = Vector<typename EMField::value_type::value_type, Dim>;
+        using SourceVector_t          = typename SourceField::value_type;
+
+        /**
+         * @brief Constructor for the FDTDSolverBase class.
+         *
+         * @param source Reference to the source field.
+         * @param E Reference to the electric field.
+         * @param B Reference to the magnetic field.
+         */
+        FDTDSolverBase(SourceField& source, EMField& E, EMField& B);
+
+        /**
+         * @brief Solves the FDTD equations.
+         */
+        void solve() override;
+
+        /**
+         * @brief Sets periodic boundary conditions.
+         */
+        void setPeriodicBoundaryConditions();
+
+        /**
+         * @brief Gets the time step size.
+         *
+         * @return The time step size.
+         */
+        scalar getDt() const { return dt; }
+
+        SourceField A_n;    // Current field.
+        SourceField A_np1;  // Field at the next time step.
+        SourceField A_nm1;  // Field at the previous time step.
+
+        /**
+         * @brief Shifts the saved fields in time.
+         */
+        void timeShift();
+
+        /**
+         * @brief Steps the solver forward in time. This is a pure virtual function.
+         */
+        virtual void step() = 0;
+        /**
+         * @brief Evaluates the electric and magnetic fields.
+         */
+        void evaluate_EB();
+
+        /**
+         * @brief Initializes the solver. This is a pure virtual function.
+         */
+        virtual void initialize() = 0;
+
+    protected:
+        /**
+         * @brief Applies the boundary conditions.
+         */
+        void applyBCs();
+
+        typename SourceField::Mesh_t* mesh_mp;  // Pointer to the mesh.
+        FieldLayout<Dim>* layout_mp;            // Pointer to the layout
+        NDIndex<Dim> domain_m;                  // Domain of the mesh.
+        Vector_t hr_m;                          // Mesh spacing.
+
+        Vector<int, Dim> nr_m;  // Number of grid points in each direction.
+        scalar dt;              // Time step size.
+    };
+}  // namespace ippl
+
+#include "FDTDSolverBase.hpp"
+#endif

src/MaxwellSolvers/FDTDSolverBase.hpp

@@ -0,0 +1,141 @@
+//
+// Class FDTDSolverBase
+//  Base class for solvers for Maxwell's equations using the FDTD method
+//
+
+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,
+                                                                              EMField& B) {
+        Maxwell<EMField, SourceField>::setSources(source);
+        Maxwell<EMField, SourceField>::setEMFields(E, B);
+    }
+
+    /**
+     * @brief Solves the FDTD equations.
+     *
+     * Advances the simulation by one time step, shifts the time for the fields, and evaluates the
+     * electric and magnetic fields at the new time.
+     */
+    template <typename EMField, typename SourceField, fdtd_bc boundary_conditions>
+    void FDTDSolverBase<EMField, SourceField, boundary_conditions>::solve() {
+        step();
+        timeShift();
+        evaluate_EB();
+    }
+
+    /**
+     * @brief Sets periodic boundary conditions.
+     *
+     * Configures the solver to use periodic boundary conditions by setting the appropriate boundary
+     * conditions for the fields.
+     */
+    template <typename EMField, typename SourceField, fdtd_bc boundary_conditions>
+    void
+    FDTDSolverBase<EMField, SourceField, boundary_conditions>::setPeriodicBoundaryConditions() {
+        typename SourceField::BConds_t vector_bcs;
+        auto bcsetter_single = [&vector_bcs]<size_t Idx>(const std::index_sequence<Idx>&) {
+            vector_bcs[Idx] = std::make_shared<ippl::PeriodicFace<SourceField>>(Idx);
+            return 0;
+        };
+        auto bcsetter = [bcsetter_single]<size_t... Idx>(const std::index_sequence<Idx...>&) {
+            int x = (bcsetter_single(std::index_sequence<Idx>{}) ^ ...);
+            (void)x;
+        };
+        bcsetter(std::make_index_sequence<Dim * 2>{});
+        A_n.setFieldBC(vector_bcs);
+        A_np1.setFieldBC(vector_bcs);
+        A_nm1.setFieldBC(vector_bcs);
+    }
+
+    /**
+     * @brief Shifts the saved fields in time.
+     *
+     * Copies the current field values to the previous time step field and the next time step field
+     * values to the current field.
+     */
+    template <typename EMField, typename SourceField, fdtd_bc boundary_conditions>
+    void FDTDSolverBase<EMField, SourceField, boundary_conditions>::timeShift() {
+        // Look into this, maybe cyclic swap is better
+        Kokkos::deep_copy(this->A_nm1.getView(), this->A_n.getView());
+        Kokkos::deep_copy(this->A_n.getView(), this->A_np1.getView());
+    }
+
+    /**
+     * @brief Applies the boundary conditions.
+     *
+     * Applies the specified boundary conditions (periodic or absorbing) to the fields.
+     */
+    template <typename EMField, typename SourceField, fdtd_bc boundary_conditions>
+    void FDTDSolverBase<EMField, SourceField, boundary_conditions>::applyBCs() {
+        if constexpr (boundary_conditions == periodic) {
+            A_n.getFieldBC().apply(A_n);
+            A_nm1.getFieldBC().apply(A_nm1);
+            A_np1.getFieldBC().apply(A_np1);
+        } else {
+            Vector<uint32_t, Dim> true_nr = nr_m;
+            true_nr += (A_n.getNghost() * 2);
+            second_order_mur_boundary_conditions bcs{};
+            bcs.apply(A_n, A_nm1, A_np1, this->dt, true_nr, layout_mp->getLocalNDIndex());
+        }
+    }
+
+    /**
+     * @brief Evaluates the electric and magnetic fields.
+     *
+     * Computes the electric and magnetic fields based on the current, previous, and next time step
+     * field values, as well as the source field.
+     */
+    template <typename EMField, typename SourceField, fdtd_bc boundary_conditions>
+    void FDTDSolverBase<EMField, SourceField, boundary_conditions>::evaluate_EB() {
+        *(Maxwell<EMField, SourceField>::En_mp)   = typename EMField::value_type(0);
+        *(Maxwell<EMField, SourceField>::Bn_mp)   = typename EMField::value_type(0);
+        ippl::Vector<scalar, 3> inverse_2_spacing = ippl::Vector<scalar, 3>(0.5) / hr_m;
+        const scalar idt                          = scalar(1.0) / dt;
+        auto A_np1 = this->A_np1.getView(), A_n = this->A_n.getView(),
+             A_nm1  = this->A_nm1.getView();
+        auto source = Maxwell<EMField, SourceField>::JN_mp->getView();
+        auto Eview  = Maxwell<EMField, SourceField>::En_mp->getView();
+        auto Bview  = Maxwell<EMField, SourceField>::Bn_mp->getView();
+
+        // Calculate the electric and magnetic fields
+        // Curl and grad of ippl are not used here, because we have a 4-component vector and we need
+        // it for the last three components
+        Kokkos::parallel_for(
+            this->A_n.getFieldRangePolicy(), KOKKOS_LAMBDA(size_t i, size_t j, size_t k) {
+                ippl::Vector<scalar, 3> dAdt;
+                dAdt[0] = (A_np1(i, j, k)[1] - A_n(i, j, k)[1]) * idt;
+                dAdt[1] = (A_np1(i, j, k)[2] - A_n(i, j, k)[2]) * idt;
+                dAdt[2] = (A_np1(i, j, k)[3] - A_n(i, j, k)[3]) * idt;
+
+                ippl::Vector<scalar, 4> dAdx =
+                    (A_n(i + 1, j, k) - A_n(i - 1, j, k)) * inverse_2_spacing[0];
+                ippl::Vector<scalar, 4> dAdy =
+                    (A_n(i, j + 1, k) - A_n(i, j - 1, k)) * inverse_2_spacing[1];
+                ippl::Vector<scalar, 4> dAdz =
+                    (A_n(i, j, k + 1) - A_n(i, j, k - 1)) * inverse_2_spacing[2];
+
+                ippl::Vector<scalar, 3> grad_phi{dAdx[0], dAdy[0], dAdz[0]};
+                ippl::Vector<scalar, 3> curlA{
+                    dAdy[3] - dAdz[2],
+                    dAdz[1] - dAdx[3],
+                    dAdx[2] - dAdy[1],
+                };
+                Eview(i, j, k) = -dAdt - grad_phi;
+                Bview(i, j, k) = curlA;
+            });
+        Kokkos::fence();
+    }
+
+}  // namespace ippl

src/MaxwellSolvers/NonStandardFDTDSolver.h

@@ -0,0 +1,87 @@
+/**
+ * @file NonStandardFDTDSolver.h
+ * @brief Defines the NonStandardFDTDSolver class for solving Maxwell's equations using a
+ * non-standard Finite-Difference Time-Domain (FDTD) method. The method and derivation of the
+ * discretization are based on:
+ * - A. Taflove, *Computational Electrodynamics: The Finite-difference Time-domain Method*, Artech
+ * House, 1995
+ * - A. Fallahi, MITHRA 2.0: *A Full-Wave Simulation Tool for Free Electron Lasers*, (2020)
+ */
+
+#ifndef IPPL_NON_STANDARD_FDTD_SOLVER_H
+#define IPPL_NON_STANDARD_FDTD_SOLVER_H
+
+#include <cstddef>
+using std::size_t;
+#include "Types/Vector.h"
+
+#include "FieldLayout/FieldLayout.h"
+#include "MaxwellSolvers/AbsorbingBC.h"
+#include "MaxwellSolvers/FDTDSolverBase.h"
+#include "MaxwellSolvers/Maxwell.h"
+#include "Meshes/UniformCartesian.h"
+#include "Particle/ParticleBase.h"
+
+namespace ippl {
+
+    /**
+     * @class NonStandardFDTDSolver
+     * @brief A solver for Maxwell's equations using a non-standard Finite-Difference Time-Domain
+     * (FDTD) method.
+     *
+     * @tparam EMField The type representing the electromagnetic field.
+     * @tparam SourceField The type representing the source field.
+     * @tparam boundary_conditions The boundary conditions to be applied (default is periodic).
+     */
+
+    template <typename EMField, typename SourceField, fdtd_bc boundary_conditions = periodic>
+    class NonStandardFDTDSolver : public FDTDSolverBase<EMField, SourceField, boundary_conditions> {
+    public:
+        /**
+         * @brief Constructs a NonStandardFDTDSolver.
+         *
+         * @param source The source field.
+         * @param E The electric field.
+         * @param B The magnetic field.
+         */
+        NonStandardFDTDSolver(SourceField& source, EMField& E, EMField& B);
+
+        constexpr static unsigned Dim = EMField::dim;  // Dimension of the electromagnetic field.
+        using scalar = typename EMField::value_type::value_type;  // Scalar type used in the
+                                                                  // electromagnetic field.
+        using Vector_t = Vector<typename EMField::value_type::value_type,
+                                Dim>;  // Vector type used in the electromagnetic field.
+        using SourceVector_t =
+            typename SourceField::value_type;  // Vector type used in the source field.
+
+        /**
+         * @brief A structure representing nondispersive coefficients.
+         *
+         * This structure holds the coefficients used in the nondispersive
+         * formulation of the solver. These coefficients are used to update
+         * the electromagnetic fields during the simulation.
+         *
+         * @tparam scalar The type of the coefficients.
+         */
+        template <typename scalar>
+        struct nondispersive {
+            scalar a1;
+            scalar a2;
+            scalar a4;
+            scalar a6;
+            scalar a8;
+        };
+
+        /**
+         * @brief Advances the simulation by one time step.
+         */
+        void step() override;
+        /**
+         * @brief Initializes the solver.
+         */
+        void initialize() override;
+    };
+}  // namespace ippl
+
+#include "NonStandardFDTDSolver.hpp"
+#endif