diff --git a/cmake/Dependencies.cmake b/cmake/Dependencies.cmake index 955d292e2..a87953f83 100644 --- a/cmake/Dependencies.cmake +++ b/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() \ No newline at end of file diff --git a/src/MaxwellSolvers/AbsorbingBC.h b/src/MaxwellSolvers/AbsorbingBC.h new file mode 100644 index 000000000..81474241d --- /dev/null +++ b/src/MaxwellSolvers/AbsorbingBC.h @@ -0,0 +1,627 @@ +#ifndef ABC_H +#define ABC_H +#include + +#include "Types/Vector.h" + +#include "Field/Field.h" + +template +constexpr KOKKOS_INLINE_FUNCTION auto first() { + return a; +} +template +constexpr KOKKOS_INLINE_FUNCTION auto second() { + return b; +} + +/* @struct second_order_abc_face + * @brief Implements a second-order absorbing boundary condition (ABC) for a face of a 3D + * computational domain. + * + * This struct computes and applies the second-order Mur ABC on a specified face of the domain, + * minimizing reflections for electromagnetic field solvers. The face is defined by its main axis + * (the normal direction) and two side axes (the tangential directions). + * + * The struct precomputes weights based on mesh spacing and time step, and provides an operator() + * to update the field at the boundary using current, previous, and next time-step values. + * + * @tparam _scalar Scalar type (e.g., float or double). + * @tparam _main_axis The main axis (0=x, 1=y, 2=z) normal to the face. + * @tparam _side_axes The two side axes (tangential to the face). + */ +template +struct second_order_abc_face { + using scalar = _scalar; // Define the scalar type for the weights and calculations + scalar Cweights[5]; // Weights for the second-order ABC, precomputed based on mesh spacing and + // time step + int sign; // Sign of the main axis (1 for lower boundary, -1 for upper). + constexpr static unsigned main_axis = + _main_axis; // Main axis normal to the face (0=x, 1=y, 2=z) + + /*! + * @brief Constructor for the second-order ABC face. + * @param hr Mesh spacing in each dimension. + * @param dt Time step size. + * @param _sign Sign of the main axis (1 or -1). + */ + KOKKOS_FUNCTION second_order_abc_face(ippl::Vector hr, scalar dt, int _sign) + : sign(_sign) { + // Speed of light in the medium + constexpr scalar c = 1; + + // Check that the main axis is one of the three axes and the side axes are different from + // each other and from the main axis + constexpr unsigned side_axes[2] = {_side_axes...}; + static_assert( + (main_axis == 0 && first<_side_axes...>() == 1 && second<_side_axes...>() == 2) + || (main_axis == 1 && first<_side_axes...>() == 0 && second<_side_axes...>() == 2) + || (main_axis == 2 && first<_side_axes...>() == 0 && second<_side_axes...>() == 1)); + assert(_main_axis != side_axes[0]); + assert(_main_axis != side_axes[1]); + assert(side_axes[0] != side_axes[1]); + + // Calculate prefactor for the weights used to calculate A_np1 at the boundary + scalar d = 1.0 / (2.0 * dt * hr[main_axis]) + 1.0 / (2.0 * c * dt * dt); + + // Calculate the weights for the boundary condition + Cweights[0] = (1.0 / (2.0 * dt * hr[main_axis]) - 1.0 / (2.0 * c * dt * dt)) / d; + Cweights[1] = (-1.0 / (2.0 * dt * hr[main_axis]) - 1.0 / (2.0 * c * dt * dt)) / d; + assert(abs(Cweights[1] + 1) < 1e-6); // Cweight[1] should always be -1 + Cweights[2] = (1.0 / (c * dt * dt) - c / (2.0 * hr[side_axes[0]] * hr[side_axes[0]]) + - c / (2.0 * hr[side_axes[1]] * hr[side_axes[1]])) + / d; + Cweights[3] = (c / (4.0 * hr[side_axes[0]] * hr[side_axes[0]])) / d; + Cweights[4] = (c / (4.0 * hr[side_axes[1]] * hr[side_axes[1]])) / d; + } + + /*! + * @brief Applies the second-order ABC to the boundary of the field. + * @param A_n Current time step field. + * @param A_nm1 Previous time step field. + * @param A_np1 Next time step field. + * @param c Coordinates of the current point in the field. + * @return Updated value at the boundary. + */ + template + KOKKOS_INLINE_FUNCTION auto operator()(const view_type& A_n, const view_type& A_nm1, + const view_type& A_np1, const Coords& c) const -> + typename view_type::value_type { + uint32_t i = c[0]; + uint32_t j = c[1]; + uint32_t k = c[2]; + using ippl::apply; + constexpr unsigned side_axes[2] = {_side_axes...}; + + // Create vectors with 1 in the direction of the specified axes and 0 in the others + ippl::Vector side_axis1_onehot = + ippl::Vector{side_axes[0] == 0, side_axes[0] == 1, side_axes[0] == 2}; + ippl::Vector side_axis2_onehot = + ippl::Vector{side_axes[1] == 0, side_axes[1] == 1, side_axes[1] == 2}; + // Create a vector in the direction of the main axis, the sign is used to determine the + // direction of the offset depending on which side of the boundary we are at + ippl::Vector mainaxis_off = ippl::Vector{ + (main_axis == 0) * sign, (main_axis == 1) * sign, (main_axis == 2) * sign}; + + // Apply the boundary condition + // The main axis is the one that is normal to the face, so we need to apply the boundary + // condition in the direction of the main axis, this means we need to offset the coordinates + // by the sign of the main axis + return advanceBoundaryS( + A_nm1(i, j, k), A_n(i, j, k), apply(A_nm1, c + mainaxis_off), + apply(A_n, c + mainaxis_off), apply(A_np1, c + mainaxis_off), + apply(A_n, c + side_axis1_onehot + mainaxis_off), + apply(A_n, c - side_axis1_onehot + mainaxis_off), + apply(A_n, c + side_axis2_onehot + mainaxis_off), + apply(A_n, c - side_axis2_onehot + mainaxis_off), apply(A_n, c + side_axis1_onehot), + apply(A_n, c - side_axis1_onehot), apply(A_n, c + side_axis2_onehot), + apply(A_n, c - side_axis2_onehot)); + } + + /*! + * @brief Advances the boundary condition using the precomputed weights. + */ + template + KOKKOS_FUNCTION value_type advanceBoundaryS(const value_type& v1, const value_type& v2, + const value_type& v3, const value_type& v4, + const value_type& v5, const value_type& v6, + const value_type& v7, const value_type& v8, + const value_type& v9, const value_type& v10, + const value_type& v11, const value_type& v12, + const value_type& v13) const noexcept { + value_type v0 = Cweights[0] * (v1 + v5) + (Cweights[1]) * v3 + (Cweights[2]) * (v2 + v4) + + (Cweights[3]) * (v6 + v7 + v10 + v11) + + (Cweights[4]) * (v8 + v9 + v12 + v13); + return v0; + } +}; + +/* @struct second_order_abc_edge + * @brief Implements a second-order absorbing boundary condition (ABC) for an edge of a 3D + * computational domain. + * + * This struct computes and applies the second-order Mur ABC on a specified edge of the domain, + * minimizing reflections for electromagnetic field solvers. The edge is defined by its axis (the + * edge direction) and two normal axes (the directions normal to the edge). + * + * The struct precomputes weights based on mesh spacing and time step, and provides an operator() + * to update the field at the boundary using current, previous, and next time-step values. + * + * @tparam _scalar Scalar type (e.g., float or double). + * @tparam edge_axis The axis (0=x, 1=y, 2=z) along the edge. + * @tparam normal_axis1 The first normal axis to the edge. + * @tparam normal_axis2 The second normal axis to the edge. + * @tparam na1_zero Boolean indicating direction for normal_axis1 (true for lower boundary, + * false for upper). + * @tparam na2_zero Boolean indicating direction for normal_axis2 (true for lower boundary, + * false for upper). + */ +template +struct second_order_abc_edge { + using scalar = _scalar; // Define the scalar type for the weights and calculations + scalar Eweights[5]; // Weights for the second-order ABC, precomputed based on mesh spacing and + // time step + + /*! + * @brief Constructor for the second-order ABC edge. + * @param hr Mesh spacing in each dimension. + * @param dt Time step size. + */ + KOKKOS_FUNCTION second_order_abc_edge(ippl::Vector hr, scalar dt) { + // Check that the edge axis is one of the three axes and the normal axes are different from + // each other and from the edge axis + static_assert(normal_axis1 != normal_axis2); + static_assert(edge_axis != normal_axis2); + static_assert(edge_axis != normal_axis1); + // Check that the axes are in the correct order + static_assert((edge_axis == 2 && normal_axis1 == 0 && normal_axis2 == 1) + || (edge_axis == 0 && normal_axis1 == 1 && normal_axis2 == 2) + || (edge_axis == 1 && normal_axis1 == 2 && normal_axis2 == 0)); + + // Speed of light in the medium + constexpr scalar c0_ = scalar(1); + + // Calculate prefactor for the weights used to calculate A_np1 at the boundary + scalar d = (1.0 / hr[normal_axis1] + 1.0 / hr[normal_axis2]) / (4.0 * dt) + + 3.0 / (8.0 * c0_ * dt * dt); + + // Calculate the weights for the boundary condition + Eweights[0] = (-(1.0 / hr[normal_axis2] - 1.0 / hr[normal_axis1]) / (4.0 * dt) + - 3.0 / (8.0 * c0_ * dt * dt)) + / d; + Eweights[1] = ((1.0 / hr[normal_axis2] - 1.0 / hr[normal_axis1]) / (4.0 * dt) + - 3.0 / (8.0 * c0_ * dt * dt)) + / d; + Eweights[2] = ((1.0 / hr[normal_axis2] + 1.0 / hr[normal_axis1]) / (4.0 * dt) + - 3.0 / (8.0 * c0_ * dt * dt)) + / d; + Eweights[3] = + (3.0 / (4.0 * c0_ * dt * dt) - c0_ / (4.0 * hr[edge_axis] * hr[edge_axis])) / d; + Eweights[4] = c0_ / (8.0 * hr[edge_axis] * hr[edge_axis]) / d; + } + + /*! + * @brief Applies the second-order ABC to the edge of the field. + * @param A_n Current time step field. + * @param A_nm1 Previous time step field. + * @param A_np1 Next time step field. + * @param c Coordinates of the current point in the field. + * @return Updated value at the edge. + */ + template + KOKKOS_INLINE_FUNCTION auto operator()(const view_type& A_n, const view_type& A_nm1, + const view_type& A_np1, const Coords& c) const -> + typename view_type::value_type { + uint32_t i = c[0]; + uint32_t j = c[1]; + uint32_t k = c[2]; + using ippl::apply; + + // Get the opposite direction of the normal vectors for the two normal axes (direction to + // the interior of the domain), 1 or -1 in the direction of the normal axis and 0 in the others + ippl::Vector normal_axis1_onehot = + ippl::Vector{normal_axis1 == 0, normal_axis1 == 1, normal_axis1 == 2} + * int32_t(na1_zero ? 1 : -1); + ippl::Vector normal_axis2_onehot = + ippl::Vector{normal_axis2 == 0, normal_axis2 == 1, normal_axis2 == 2} + * int32_t(na2_zero ? 1 : -1); + + // Calculate the coordinates of the neighboring points in the interior of the domain + ippl::Vector acc0 = {i, j, k}; + ippl::Vector acc1 = acc0 + normal_axis1_onehot; + ippl::Vector acc2 = acc0 + normal_axis2_onehot; + ippl::Vector acc3 = acc0 + normal_axis1_onehot + normal_axis2_onehot; + + // Create a vector in the direction of the edge axis + ippl::Vector axisp{edge_axis == 0, edge_axis == 1, edge_axis == 2}; + + // Apply the boundary condition + return advanceEdgeS(A_n(i, j, k), A_nm1(i, j, k), apply(A_np1, acc1), apply(A_n, acc1), + apply(A_nm1, acc1), apply(A_np1, acc2), apply(A_n, acc2), + apply(A_nm1, acc2), apply(A_np1, acc3), apply(A_n, acc3), + apply(A_nm1, acc3), apply(A_n, acc0 - axisp), apply(A_n, acc1 - axisp), + apply(A_n, acc2 - axisp), apply(A_n, acc3 - axisp), + apply(A_n, acc0 + axisp), apply(A_n, acc1 + axisp), + apply(A_n, acc2 + axisp), apply(A_n, acc3 + axisp)); + } + + /*! + * @brief Advances the edge boundary condition using the precomputed weights. + */ + template + KOKKOS_INLINE_FUNCTION value_type advanceEdgeS(value_type v1, value_type v2, value_type v3, + value_type v4, value_type v5, value_type v6, + value_type v7, value_type v8, value_type v9, + value_type v10, value_type v11, value_type v12, + value_type v13, value_type v14, value_type v15, + value_type v16, value_type v17, value_type v18, + value_type v19) const noexcept { + value_type v0 = Eweights[0] * (v3 + v8) + Eweights[1] * (v5 + v6) + Eweights[2] * (v2 + v9) + + Eweights[3] * (v1 + v4 + v7 + v10) + + Eweights[4] * (v12 + v13 + v14 + v15 + v16 + v17 + v18 + v19) - v11; + return v0; + } +}; + +/* @struct second_order_abc_corner + * @brief Implements a second-order absorbing boundary condition (ABC) for a corner of a 3D + * computational domain. + * + * This struct computes and applies the second-order Mur ABC on a specified corner of the domain, + * minimizing reflections for electromagnetic field solvers. The corner is defined by its three + * axes, with each axis having a boolean indicating whether it is at the lower or upper boundary. + * + * The struct precomputes weights based on mesh spacing and time step, and provides an operator() + * to update the field at the boundary using current, previous, and next time-step values. + * + * @tparam _scalar Scalar type (e.g., float or double). + * @tparam x0 Boolean indicating direction for x-axis (true for lower boundary, false for + * upper). + * @tparam y0 Boolean indicating direction for y-axis (true for lower boundary, false for + * upper). + * @tparam z0 Boolean indicating direction for z-axis (true for lower boundary, false for + * upper). + */ +template +struct second_order_abc_corner { + using scalar = _scalar; // Define the scalar type for the weights and calculations + scalar Cweights[17]; // Weights for the second-order ABC, precomputed based on mesh spacing and + // time step + + /*! + * @brief Constructor for the second-order ABC corner. + * @param hr Mesh spacing in each dimension. + * @param dt Time step size. + */ + KOKKOS_FUNCTION second_order_abc_corner(ippl::Vector hr, scalar dt) { + constexpr scalar c0_ = scalar(1); + Cweights[0] = + (-1.0 / hr[0] - 1.0 / hr[1] - 1.0 / hr[2]) / (8.0 * dt) - 1.0 / (4.0 * c0_ * dt * dt); + Cweights[1] = + (1.0 / hr[0] - 1.0 / hr[1] - 1.0 / hr[2]) / (8.0 * dt) - 1.0 / (4.0 * c0_ * dt * dt); + Cweights[2] = + (-1.0 / hr[0] + 1.0 / hr[1] - 1.0 / hr[2]) / (8.0 * dt) - 1.0 / (4.0 * c0_ * dt * dt); + Cweights[3] = + (-1.0 / hr[0] - 1.0 / hr[1] + 1.0 / hr[2]) / (8.0 * dt) - 1.0 / (4.0 * c0_ * dt * dt); + Cweights[4] = + (1.0 / hr[0] + 1.0 / hr[1] - 1.0 / hr[2]) / (8.0 * dt) - 1.0 / (4.0 * c0_ * dt * dt); + Cweights[5] = + (1.0 / hr[0] - 1.0 / hr[1] + 1.0 / hr[2]) / (8.0 * dt) - 1.0 / (4.0 * c0_ * dt * dt); + Cweights[6] = + (-1.0 / hr[0] + 1.0 / hr[1] + 1.0 / hr[2]) / (8.0 * dt) - 1.0 / (4.0 * c0_ * dt * dt); + Cweights[7] = + (1.0 / hr[0] + 1.0 / hr[1] + 1.0 / hr[2]) / (8.0 * dt) - 1.0 / (4.0 * c0_ * dt * dt); + Cweights[8] = + -(-1.0 / hr[0] - 1.0 / hr[1] - 1.0 / hr[2]) / (8.0 * dt) - 1.0 / (4.0 * c0_ * dt * dt); + Cweights[9] = + -(1.0 / hr[0] - 1.0 / hr[1] - 1.0 / hr[2]) / (8.0 * dt) - 1.0 / (4.0 * c0_ * dt * dt); + Cweights[10] = + -(-1.0 / hr[0] + 1.0 / hr[1] - 1.0 / hr[2]) / (8.0 * dt) - 1.0 / (4.0 * c0_ * dt * dt); + Cweights[11] = + -(-1.0 / hr[0] - 1.0 / hr[1] + 1.0 / hr[2]) / (8.0 * dt) - 1.0 / (4.0 * c0_ * dt * dt); + Cweights[12] = + -(1.0 / hr[0] + 1.0 / hr[1] - 1.0 / hr[2]) / (8.0 * dt) - 1.0 / (4.0 * c0_ * dt * dt); + Cweights[13] = + -(1.0 / hr[0] - 1.0 / hr[1] + 1.0 / hr[2]) / (8.0 * dt) - 1.0 / (4.0 * c0_ * dt * dt); + Cweights[14] = + -(-1.0 / hr[0] + 1.0 / hr[1] + 1.0 / hr[2]) / (8.0 * dt) - 1.0 / (4.0 * c0_ * dt * dt); + Cweights[15] = + -(1.0 / hr[0] + 1.0 / hr[1] + 1.0 / hr[2]) / (8.0 * dt) - 1.0 / (4.0 * c0_ * dt * dt); + Cweights[16] = 1.0 / (2.0 * c0_ * dt * dt); + } + + /*! + * @brief Applies the second-order ABC to the corner of the field. + * @param A_n Current time step field. + * @param A_nm1 Previous time step field. + * @param A_np1 Next time step field. + * @param c Coordinates of the current point in the field. + * @return Updated value at the corner. + */ + template + KOKKOS_INLINE_FUNCTION auto operator()(const view_type& A_n, const view_type& A_nm1, + const view_type& A_np1, const Coords& c) const -> + typename view_type::value_type { + // Get direction of the vector to the interior of the domain + constexpr uint32_t xoff = (x0) ? 1 : uint32_t(-1); + constexpr uint32_t yoff = (y0) ? 1 : uint32_t(-1); + constexpr uint32_t zoff = (z0) ? 1 : uint32_t(-1); + using ippl::apply; + + // Calculate the offset vectors for the neighboring points in the interior of the domain + const ippl::Vector offsets[8] = { + ippl::Vector{0, 0, 0}, ippl::Vector{xoff, 0, 0}, + ippl::Vector{0, yoff, 0}, ippl::Vector{0, 0, zoff}, + ippl::Vector{xoff, yoff, 0}, ippl::Vector{xoff, 0, zoff}, + ippl::Vector{0, yoff, zoff}, ippl::Vector{xoff, yoff, zoff}, + }; + + // Apply the boundary condition + return advanceCornerS( + apply(A_n, c), apply(A_nm1, c), apply(A_np1, c + offsets[1]), + apply(A_n, c + offsets[1]), apply(A_nm1, c + offsets[1]), apply(A_np1, c + offsets[2]), + apply(A_n, c + offsets[2]), apply(A_nm1, c + offsets[2]), apply(A_np1, c + offsets[3]), + apply(A_n, c + offsets[3]), apply(A_nm1, c + offsets[3]), apply(A_np1, c + offsets[4]), + apply(A_n, c + offsets[4]), apply(A_nm1, c + offsets[4]), apply(A_np1, c + offsets[5]), + apply(A_n, c + offsets[5]), apply(A_nm1, c + offsets[5]), apply(A_np1, c + offsets[6]), + apply(A_n, c + offsets[6]), apply(A_nm1, c + offsets[6]), apply(A_np1, c + offsets[7]), + apply(A_n, c + offsets[7]), apply(A_nm1, c + offsets[7])); + } + + /*! + * @brief Advances the corner boundary condition using the precomputed weights. + */ + template + KOKKOS_INLINE_FUNCTION value_type + advanceCornerS(value_type v1, value_type v2, value_type v3, value_type v4, value_type v5, + value_type v6, value_type v7, value_type v8, value_type v9, value_type v10, + value_type v11, value_type v12, value_type v13, value_type v14, value_type v15, + value_type v16, value_type v17, value_type v18, value_type v19, value_type v20, + value_type v21, value_type v22, value_type v23) const noexcept { + return -(v1 * (Cweights[16]) + v2 * (Cweights[8]) + v3 * Cweights[1] + v4 * Cweights[16] + + v5 * Cweights[9] + v6 * Cweights[2] + v7 * Cweights[16] + v8 * Cweights[10] + + v9 * Cweights[3] + v10 * Cweights[16] + v11 * Cweights[11] + v12 * Cweights[4] + + v13 * Cweights[16] + v14 * Cweights[12] + v15 * Cweights[5] + v16 * Cweights[16] + + v17 * Cweights[13] + v18 * Cweights[6] + v19 * Cweights[16] + v20 * Cweights[14] + + v21 * Cweights[7] + v22 * Cweights[16] + v23 * Cweights[15]) + / Cweights[0]; + } +}; + +/* @struct second_order_mur_boundary_conditions + * @brief Applies second-order Mur absorbing boundary conditions (ABC) to all boundaries of a 3D + * computational domain. + * + * This struct coordinates the application of second-order Mur ABCs on faces, edges, and corners of + * the domain, minimizing reflections for electromagnetic field solvers. It uses the appropriate + * face, edge, or corner ABC implementation depending on the location of each boundary point. + * + * The apply() method detects the type of boundary (face, edge, or corner) for each point and + * applies the corresponding second-order ABC using precomputed weights based on mesh spacing and + * time step. + * + * @tparam field_type The field type (must provide getView() and get_mesh().getMeshSpacing()). + * @tparam dt_type The type of the time step. + */ +struct second_order_mur_boundary_conditions { + /*! + * @brief Applies second-order Mur ABC to the boundaries of the field. + * @param FA_n Field at the current time step. + * @param FA_nm1 Field at the previous time step. + * @param FA_np1 Field at the next time step. + * @param dt Time step size. + * @param true_nr True number of grid points in each dimension. + * @param lDom Local domain index for parallel processing. + */ + template + void apply(field_type& FA_n, field_type& FA_nm1, field_type& FA_np1, dt_type dt, + ippl::Vector true_nr, ippl::NDIndex<3> lDom) { + using scalar = decltype(dt); + + // Get the mesh spacing + const ippl::Vector hr = FA_n.get_mesh().getMeshSpacing(); + + // Get views of the fields + auto A_n = FA_n.getView(); + auto A_np1 = FA_np1.getView(); + auto A_nm1 = FA_nm1.getView(); + + // Get the local number of grid points in each dimension + ippl::Vector local_nr{uint32_t(A_n.extent(0)), uint32_t(A_n.extent(1)), + uint32_t(A_n.extent(2))}; + + constexpr uint32_t min_abc_boundary = 1; + constexpr uint32_t max_abc_boundary_sub = min_abc_boundary + 1; + + // Apply second-order ABC to faces + Kokkos::parallel_for( + ippl::getRangePolicy(A_n, 1), KOKKOS_LAMBDA(uint32_t i, uint32_t j, uint32_t k) { + // Get the global indices of the current point + uint32_t ig = i + lDom.first()[0]; + uint32_t jg = j + lDom.first()[1]; + uint32_t kg = k + lDom.first()[2]; + + // Check on how many boundaries of the local domain the current point is located + uint32_t lval = uint32_t(i == 0) + (uint32_t(j == 0) << 1) + (uint32_t(k == 0) << 2) + + (uint32_t(i == local_nr[0] - 1) << 3) + + (uint32_t(j == local_nr[1] - 1) << 4) + + (uint32_t(k == local_nr[2] - 1) << 5); + + // Exits current iteration if the point is on an edge or corner of the local domain + if (Kokkos::popcount(lval) > 1) + return; + + // Check on how many boundaries of the global domain the current point is located + uint32_t val = uint32_t(ig == min_abc_boundary) + + (uint32_t(jg == min_abc_boundary) << 1) + + (uint32_t(kg == min_abc_boundary) << 2) + + (uint32_t(ig == true_nr[0] - max_abc_boundary_sub) << 3) + + (uint32_t(jg == true_nr[1] - max_abc_boundary_sub) << 4) + + (uint32_t(kg == true_nr[2] - max_abc_boundary_sub) << 5); + + // If the point is on a global face, apply the second-order ABC + if (Kokkos::popcount(val) == 1) { + if (ig == min_abc_boundary) { + second_order_abc_face soa(hr, dt, 1); + A_np1(i, j, k) = soa(A_n, A_nm1, A_np1, ippl::Vector{i, j, k}); + } + if (jg == min_abc_boundary) { + second_order_abc_face soa(hr, dt, 1); + A_np1(i, j, k) = soa(A_n, A_nm1, A_np1, ippl::Vector{i, j, k}); + } + if (kg == min_abc_boundary) { + second_order_abc_face soa(hr, dt, 1); + A_np1(i, j, k) = soa(A_n, A_nm1, A_np1, ippl::Vector{i, j, k}); + } + if (ig == true_nr[0] - max_abc_boundary_sub) { + second_order_abc_face soa(hr, dt, -1); + A_np1(i, j, k) = soa(A_n, A_nm1, A_np1, ippl::Vector{i, j, k}); + } + if (jg == true_nr[1] - max_abc_boundary_sub) { + second_order_abc_face soa(hr, dt, -1); + A_np1(i, j, k) = soa(A_n, A_nm1, A_np1, ippl::Vector{i, j, k}); + } + if (kg == true_nr[2] - max_abc_boundary_sub) { + second_order_abc_face soa(hr, dt, -1); + A_np1(i, j, k) = soa(A_n, A_nm1, A_np1, ippl::Vector{i, j, k}); + } + } + }); + Kokkos::fence(); + + // Apply second-order ABC to edges + Kokkos::parallel_for( + ippl::getRangePolicy(A_n, 1), KOKKOS_LAMBDA(uint32_t i, uint32_t j, uint32_t k) { + // Get the global indices of the current point + uint32_t ig = i + lDom.first()[0]; + uint32_t jg = j + lDom.first()[1]; + uint32_t kg = k + lDom.first()[2]; + + // Check on how many boundaries of the local domain the current point is located + uint32_t lval = uint32_t(i == 0) + (uint32_t(j == 0) << 1) + (uint32_t(k == 0) << 2) + + (uint32_t(i == local_nr[0] - 1) << 3) + + (uint32_t(j == local_nr[1] - 1) << 4) + + (uint32_t(k == local_nr[2] - 1) << 5); + + // Exits current iteration if the point is on a corner of the local domain + if (Kokkos::popcount(lval) > 2) + return; + + // Check on how many boundaries of the global domain the current point is located + uint32_t val = uint32_t(ig == min_abc_boundary) + + (uint32_t(jg == min_abc_boundary) << 1) + + (uint32_t(kg == min_abc_boundary) << 2) + + (uint32_t(ig == true_nr[0] - max_abc_boundary_sub) << 3) + + (uint32_t(jg == true_nr[1] - max_abc_boundary_sub) << 4) + + (uint32_t(kg == true_nr[2] - max_abc_boundary_sub) << 5); + + // If the point is on a global edge, apply the second-order ABC + if (Kokkos::popcount(val) == 2) { + if (ig == min_abc_boundary && kg == min_abc_boundary) { + second_order_abc_edge soa(hr, dt); + A_np1(i, j, k) = soa(A_n, A_nm1, A_np1, ippl::Vector{i, j, k}); + } else if (ig == min_abc_boundary && jg == min_abc_boundary) { + second_order_abc_edge soa(hr, dt); + A_np1(i, j, k) = soa(A_n, A_nm1, A_np1, ippl::Vector{i, j, k}); + } else if (jg == min_abc_boundary && kg == min_abc_boundary) { + second_order_abc_edge soa(hr, dt); + A_np1(i, j, k) = soa(A_n, A_nm1, A_np1, ippl::Vector{i, j, k}); + } else if (ig == min_abc_boundary && kg == true_nr[2] - max_abc_boundary_sub) { + second_order_abc_edge soa(hr, dt); + A_np1(i, j, k) = soa(A_n, A_nm1, A_np1, ippl::Vector{i, j, k}); + } else if (ig == min_abc_boundary && jg == true_nr[1] - max_abc_boundary_sub) { + second_order_abc_edge soa(hr, dt); + A_np1(i, j, k) = soa(A_n, A_nm1, A_np1, ippl::Vector{i, j, k}); + } else if (jg == min_abc_boundary && kg == true_nr[2] - max_abc_boundary_sub) { + second_order_abc_edge soa(hr, dt); + A_np1(i, j, k) = soa(A_n, A_nm1, A_np1, ippl::Vector{i, j, k}); + } else if (ig == true_nr[0] - max_abc_boundary_sub && kg == min_abc_boundary) { + second_order_abc_edge soa(hr, dt); + A_np1(i, j, k) = soa(A_n, A_nm1, A_np1, ippl::Vector{i, j, k}); + } else if (ig == true_nr[0] - max_abc_boundary_sub && jg == min_abc_boundary) { + second_order_abc_edge soa(hr, dt); + A_np1(i, j, k) = soa(A_n, A_nm1, A_np1, ippl::Vector{i, j, k}); + } else if (jg == true_nr[1] - max_abc_boundary_sub && kg == min_abc_boundary) { + second_order_abc_edge soa(hr, dt); + A_np1(i, j, k) = soa(A_n, A_nm1, A_np1, ippl::Vector{i, j, k}); + } else if (ig == true_nr[0] - max_abc_boundary_sub + && kg == true_nr[2] - max_abc_boundary_sub) { + second_order_abc_edge soa(hr, dt); + A_np1(i, j, k) = soa(A_n, A_nm1, A_np1, ippl::Vector{i, j, k}); + } else if (ig == true_nr[0] - max_abc_boundary_sub + && jg == true_nr[1] - max_abc_boundary_sub) { + second_order_abc_edge soa(hr, dt); + A_np1(i, j, k) = soa(A_n, A_nm1, A_np1, ippl::Vector{i, j, k}); + } else if (jg == true_nr[1] - max_abc_boundary_sub + && kg == true_nr[2] - max_abc_boundary_sub) { + second_order_abc_edge soa(hr, dt); + A_np1(i, j, k) = soa(A_n, A_nm1, A_np1, ippl::Vector{i, j, k}); + } else { + assert(false); + } + } + }); + Kokkos::fence(); + + // Apply second-order ABC to corners + Kokkos::parallel_for( + ippl::getRangePolicy(A_n, 1), KOKKOS_LAMBDA(uint32_t i, uint32_t j, uint32_t k) { + // Get the global indices of the current point + uint32_t ig = i + lDom.first()[0]; + uint32_t jg = j + lDom.first()[1]; + uint32_t kg = k + lDom.first()[2]; + + // Check on how many boundaries of the global domain the current point is located + uint32_t val = uint32_t(ig == min_abc_boundary) + + (uint32_t(jg == min_abc_boundary) << 1) + + (uint32_t(kg == min_abc_boundary) << 2) + + (uint32_t(ig == true_nr[0] - max_abc_boundary_sub) << 3) + + (uint32_t(jg == true_nr[1] - max_abc_boundary_sub) << 4) + + (uint32_t(kg == true_nr[2] - max_abc_boundary_sub) << 5); + + // If the point is on a global corner, apply the second-order ABC + if (Kokkos::popcount(val) == 3) { + if (ig == min_abc_boundary && jg == min_abc_boundary + && kg == min_abc_boundary) { + second_order_abc_corner coa(hr, dt); + A_np1(i, j, k) = coa(A_n, A_nm1, A_np1, ippl::Vector{i, j, k}); + } else if (ig == true_nr[0] - max_abc_boundary_sub && jg == min_abc_boundary + && kg == min_abc_boundary) { + second_order_abc_corner coa(hr, dt); + A_np1(i, j, k) = coa(A_n, A_nm1, A_np1, ippl::Vector{i, j, k}); + } else if (ig == min_abc_boundary && jg == true_nr[1] - max_abc_boundary_sub + && kg == min_abc_boundary) { + second_order_abc_corner coa(hr, dt); + A_np1(i, j, k) = coa(A_n, A_nm1, A_np1, ippl::Vector{i, j, k}); + } else if (ig == true_nr[0] - max_abc_boundary_sub + && jg == true_nr[1] - max_abc_boundary_sub + && kg == min_abc_boundary) { + second_order_abc_corner coa(hr, dt); + A_np1(i, j, k) = coa(A_n, A_nm1, A_np1, ippl::Vector{i, j, k}); + } else if (ig == min_abc_boundary && jg == min_abc_boundary + && kg == true_nr[2] - max_abc_boundary_sub) { + second_order_abc_corner coa(hr, dt); + A_np1(i, j, k) = coa(A_n, A_nm1, A_np1, ippl::Vector{i, j, k}); + } else if (ig == true_nr[0] - max_abc_boundary_sub && jg == min_abc_boundary + && kg == true_nr[2] - max_abc_boundary_sub) { + second_order_abc_corner coa(hr, dt); + A_np1(i, j, k) = coa(A_n, A_nm1, A_np1, ippl::Vector{i, j, k}); + } else if (ig == min_abc_boundary && jg == true_nr[1] - max_abc_boundary_sub + && kg == true_nr[2] - max_abc_boundary_sub) { + second_order_abc_corner coa(hr, dt); + A_np1(i, j, k) = coa(A_n, A_nm1, A_np1, ippl::Vector{i, j, k}); + } else if (ig == true_nr[0] - max_abc_boundary_sub + && jg == true_nr[1] - max_abc_boundary_sub + && kg == true_nr[2] - max_abc_boundary_sub) { + second_order_abc_corner coa(hr, dt); + A_np1(i, j, k) = coa(A_n, A_nm1, A_np1, ippl::Vector{i, j, k}); + } else { + assert(false); + } + } + }); + Kokkos::fence(); + } +}; +#endif \ No newline at end of file diff --git a/src/MaxwellSolvers/CMakeLists.txt b/src/MaxwellSolvers/CMakeLists.txt index 441f723ee..758c07106 100644 --- a/src/MaxwellSolvers/CMakeLists.txt +++ b/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 ) diff --git a/src/MaxwellSolvers/FDTDSolverBase.h b/src/MaxwellSolvers/FDTDSolverBase.h new file mode 100644 index 000000000..0a16eb3e3 --- /dev/null +++ b/src/MaxwellSolvers/FDTDSolverBase.h @@ -0,0 +1,101 @@ +#ifndef IPPL_FDTD_H +#define IPPL_FDTD_H +#include +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 + class FDTDSolverBase : public Maxwell { + public: + constexpr static unsigned Dim = EMField::dim; + using scalar = typename EMField::value_type::value_type; + using Vector_t = Vector; + 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* layout_mp; // Pointer to the layout + NDIndex domain_m; // Domain of the mesh. + Vector_t hr_m; // Mesh spacing. + + Vector nr_m; // Number of grid points in each direction. + scalar dt; // Time step size. + }; +} // namespace ippl + +#include "FDTDSolverBase.hpp" +#endif \ No newline at end of file diff --git a/src/MaxwellSolvers/FDTDSolverBase.hpp b/src/MaxwellSolvers/FDTDSolverBase.hpp new file mode 100644 index 000000000..6be31a55e --- /dev/null +++ b/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 + FDTDSolverBase::FDTDSolverBase(SourceField& source, + EMField& E, + EMField& B) { + Maxwell::setSources(source); + Maxwell::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 + void FDTDSolverBase::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 + void + FDTDSolverBase::setPeriodicBoundaryConditions() { + typename SourceField::BConds_t vector_bcs; + auto bcsetter_single = [&vector_bcs](const std::index_sequence&) { + vector_bcs[Idx] = std::make_shared>(Idx); + return 0; + }; + auto bcsetter = [bcsetter_single](const std::index_sequence&) { + int x = (bcsetter_single(std::index_sequence{}) ^ ...); + (void)x; + }; + bcsetter(std::make_index_sequence{}); + 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 + void FDTDSolverBase::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 + void FDTDSolverBase::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 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 + void FDTDSolverBase::evaluate_EB() { + *(Maxwell::En_mp) = typename EMField::value_type(0); + *(Maxwell::Bn_mp) = typename EMField::value_type(0); + ippl::Vector inverse_2_spacing = ippl::Vector(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::JN_mp->getView(); + auto Eview = Maxwell::En_mp->getView(); + auto Bview = Maxwell::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 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 dAdx = + (A_n(i + 1, j, k) - A_n(i - 1, j, k)) * inverse_2_spacing[0]; + ippl::Vector dAdy = + (A_n(i, j + 1, k) - A_n(i, j - 1, k)) * inverse_2_spacing[1]; + ippl::Vector dAdz = + (A_n(i, j, k + 1) - A_n(i, j, k - 1)) * inverse_2_spacing[2]; + + ippl::Vector grad_phi{dAdx[0], dAdy[0], dAdz[0]}; + ippl::Vector 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 \ No newline at end of file diff --git a/src/MaxwellSolvers/NonStandardFDTDSolver.h b/src/MaxwellSolvers/NonStandardFDTDSolver.h new file mode 100644 index 000000000..f33046b84 --- /dev/null +++ b/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 +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 + class NonStandardFDTDSolver : public FDTDSolverBase { + 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; // 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 + 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 \ No newline at end of file diff --git a/src/MaxwellSolvers/NonStandardFDTDSolver.hpp b/src/MaxwellSolvers/NonStandardFDTDSolver.hpp new file mode 100644 index 000000000..a5b565ef9 --- /dev/null +++ b/src/MaxwellSolvers/NonStandardFDTDSolver.hpp @@ -0,0 +1,144 @@ +/** + * @file NonStandardFDTDSolver.hpp + * @brief Impelmentation of the NonStandardFDTDSolver class functions. + */ + +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 + NonStandardFDTDSolver::NonStandardFDTDSolver( + SourceField& source, EMField& E, EMField& B) + : FDTDSolverBase(source, E, B) { + auto hx = source.get_mesh().getMeshSpacing(); + if ((hx[2] / hx[0]) * (hx[2] / hx[0]) + (hx[2] / hx[1]) * (hx[2] / hx[1]) >= 1) { + throw IpplException("NonStandardFDTDSolver Constructor", + "Dispersion-free CFL condition not satisfiable\n"); + } + initialize(); + } + + /** + * @brief Advances the simulation by one time step. + * + * This function updates the fields by one time step using non-standard FDTD method. + * It calculates the new field values based on the current and previous field values, + * as well as the source field. The boundary conditions are applied after the update. + */ + template + void NonStandardFDTDSolver::step() { + const auto& ldom = this->layout_mp->getLocalNDIndex(); + const int nghost = this->A_n.getNghost(); + const auto aview = this->A_n.getView(); + const auto anp1view = this->A_np1.getView(); + const auto anm1view = this->A_nm1.getView(); + const auto source_view = Maxwell::JN_mp->getView(); + + const scalar calA = + 0.25 + * (1 + + 0.02 + / (static_cast(Kokkos::pow(this->hr_m[2] / this->hr_m[0], 2)) + + static_cast(Kokkos::pow(this->hr_m[2] / this->hr_m[1], 2)))); + nondispersive ndisp{ + .a1 = + 2 + * (1 + - (1 - 2 * calA) * static_cast(Kokkos::pow(this->dt / this->hr_m[0], 2)) + - (1 - 2 * calA) * static_cast(Kokkos::pow(this->dt / this->hr_m[1], 2)) + - static_cast(Kokkos::pow(this->dt / this->hr_m[2], 2))), + .a2 = static_cast(Kokkos::pow(this->dt / this->hr_m[0], 2)), + .a4 = static_cast(Kokkos::pow(this->dt / this->hr_m[1], 2)), + .a6 = static_cast(Kokkos::pow(this->dt / this->hr_m[2], 2)) + - 2 * calA * static_cast(Kokkos::pow(this->dt / this->hr_m[0], 2)) + - 2 * calA * static_cast(Kokkos::pow(this->dt / this->hr_m[1], 2)), + .a8 = static_cast(Kokkos::pow(this->dt, 2))}; + Vector true_nr = this->nr_m; + true_nr += (nghost * 2); + constexpr uint32_t one_if_absorbing_otherwise_0 = + boundary_conditions == absorbing + ? 1 + : 0; // 1 if absorbing, 0 otherwise, indicates the start index of the field + Kokkos::parallel_for( + ippl::getRangePolicy(aview, nghost), KOKKOS_LAMBDA(size_t i, size_t j, size_t k) { + // global indices + uint32_t ig = i + ldom.first()[0]; + uint32_t jg = j + ldom.first()[1]; + uint32_t kg = k + ldom.first()[2]; + // check if at a boundary of the field + uint32_t val = + uint32_t(ig == one_if_absorbing_otherwise_0) + + (uint32_t(jg == one_if_absorbing_otherwise_0) << 1) + + (uint32_t(kg == one_if_absorbing_otherwise_0) << 2) + + (uint32_t(ig == true_nr[0] - one_if_absorbing_otherwise_0 - 1) << 3) + + (uint32_t(jg == true_nr[1] - one_if_absorbing_otherwise_0 - 1) << 4) + + (uint32_t(kg == true_nr[2] - one_if_absorbing_otherwise_0 - 1) << 5); + // update the interior field values + if (!val) { + anp1view(i, j, k) = + -anm1view(i, j, k) + ndisp.a1 * aview(i, j, k) + + ndisp.a2 + * (calA * aview(i + 1, j, k - 1) + (1 - 2 * calA) * aview(i + 1, j, k) + + calA * aview(i + 1, j, k + 1)) + + ndisp.a2 + * (calA * aview(i - 1, j, k - 1) + (1 - 2 * calA) * aview(i - 1, j, k) + + calA * aview(i - 1, j, k + 1)) + + ndisp.a4 + * (calA * aview(i, j + 1, k - 1) + (1 - 2 * calA) * aview(i, j + 1, k) + + calA * aview(i, j + 1, k + 1)) + + ndisp.a4 + * (calA * aview(i, j - 1, k - 1) + (1 - 2 * calA) * aview(i, j - 1, k) + + calA * aview(i, j - 1, k + 1)) + + ndisp.a6 * aview(i, j, k + 1) + ndisp.a6 * aview(i, j, k - 1) + + ndisp.a8 * source_view(i, j, k); + } + }); + Kokkos::fence(); + this->A_np1.fillHalo(); + this->applyBCs(); + } + + /** + * @brief Initializes the solver. + * + * This function initializes the solver by setting up the layout, mesh, and fields. + * It retrieves the mesh spacing, domain, and mesh size, and initializes the fields to zero. + */ + template + void NonStandardFDTDSolver::initialize() { + // get layout and mesh + this->layout_mp = &(this->JN_mp->getLayout()); + this->mesh_mp = &(this->JN_mp->get_mesh()); + + // get mesh spacing, domain, and mesh size + this->hr_m = this->mesh_mp->getMeshSpacing(); + this->dt = this->hr_m[2]; + this->domain_m = this->layout_mp->getDomain(); + for (unsigned int i = 0; i < Dim; ++i) + this->nr_m[i] = this->domain_m[i].length(); + + // initialize fields + this->A_nm1.initialize(*this->mesh_mp, *this->layout_mp); + this->A_n.initialize(*this->mesh_mp, *this->layout_mp); + this->A_np1.initialize(*this->mesh_mp, *this->layout_mp); + + this->A_nm1 = 0.0; + this->A_n = 0.0; + this->A_np1 = 0.0; + + // set the mesh domain + if constexpr (boundary_conditions == periodic) { + this->setPeriodicBoundaryConditions(); + } + }; +} // namespace ippl \ No newline at end of file diff --git a/src/MaxwellSolvers/StandardFDTDSolver.h b/src/MaxwellSolvers/StandardFDTDSolver.h new file mode 100644 index 000000000..60fd6bb37 --- /dev/null +++ b/src/MaxwellSolvers/StandardFDTDSolver.h @@ -0,0 +1,64 @@ +/** + * @file StandardFDTDSolver.h + * @brief Defines the StandardFDTDSolver class for solving Maxwell's equations using the FDTD + * method. + */ + +#ifndef IPPL_STANDARD_FDTD_SOLVER_H +#define IPPL_STANDARD_FDTD_SOLVER_H + +#include +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 StandardFDTDSolver + * @brief A solver for Maxwell's equations using the 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 + class StandardFDTDSolver : public FDTDSolverBase { + public: + /** + * @brief Constructs a StandardFDTDSolver. + * + * @param source The source field. + * @param E The electric field. + * @param B The magnetic field. + */ + StandardFDTDSolver(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; // Vector type used in the electromagnetic field. + using SourceVector_t = + typename SourceField::value_type; // Vector type used in the source field. + + /** + * @brief Advances the simulation by one time step. + */ + virtual void step() override; + /** + * @brief Initializes the solver. + */ + virtual void initialize() override; + }; +} // namespace ippl + +#include "StandardFDTDSolver.hpp" +#endif \ No newline at end of file diff --git a/src/MaxwellSolvers/StandardFDTDSolver.hpp b/src/MaxwellSolvers/StandardFDTDSolver.hpp new file mode 100644 index 000000000..beefb3660 --- /dev/null +++ b/src/MaxwellSolvers/StandardFDTDSolver.hpp @@ -0,0 +1,125 @@ +/** + * @file StandardFDTDSolver.hpp + * @brief Impelmentation of the StandardFDTDSolver class functions. + */ + +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 + StandardFDTDSolver::StandardFDTDSolver( + SourceField& source, EMField& E, EMField& B) + : FDTDSolverBase(source, E, B) { + initialize(); + } + + /** + * @brief Advances the simulation by one time step. + * + * This function updates the fields by one time step using the FDTD method. + * It calculates the new field values based on the current and previous field values, + * as well as the source field. The boundary conditions are applied after the update. + */ + template + void StandardFDTDSolver::step() { + const auto& ldom = this->layout_mp->getLocalNDIndex(); + const int nghost = this->A_n.getNghost(); + const auto aview = this->A_n.getView(); + const auto anp1view = this->A_np1.getView(); + const auto anm1view = this->A_nm1.getView(); + const auto source_view = Maxwell::JN_mp->getView(); + + // Calculate the coefficients for the nondispersive update + const scalar a1 = scalar(2) + * (scalar(1) - Kokkos::pow(this->dt / this->hr_m[0], 2) + - Kokkos::pow(this->dt / this->hr_m[1], 2) + - Kokkos::pow(this->dt / this->hr_m[2], 2)); + const scalar a2 = Kokkos::pow(this->dt / this->hr_m[0], 2); + const scalar a4 = Kokkos::pow(this->dt / this->hr_m[1], 2); + const scalar a6 = Kokkos::pow(this->dt / this->hr_m[2], 2); + const scalar a8 = Kokkos::pow(this->dt, 2); + Vector true_nr = this->nr_m; + true_nr += (nghost * 2); + constexpr uint32_t one_if_absorbing_otherwise_0 = + boundary_conditions == absorbing + ? 1 + : 0; // 1 if absorbing, 0 otherwise, indicates the start index of the field + + // Update the field values + Kokkos::parallel_for( + "Source field update", ippl::getRangePolicy(aview, nghost), + KOKKOS_LAMBDA(const size_t i, const size_t j, const size_t k) { + // global indices + const uint32_t ig = i + ldom.first()[0]; + const uint32_t jg = j + ldom.first()[1]; + const uint32_t kg = k + ldom.first()[2]; + // check if at a boundary of the field + uint32_t val = + uint32_t(ig == one_if_absorbing_otherwise_0) + + (uint32_t(jg == one_if_absorbing_otherwise_0) << 1) + + (uint32_t(kg == one_if_absorbing_otherwise_0) << 2) + + (uint32_t(ig == true_nr[0] - one_if_absorbing_otherwise_0 - 1) << 3) + + (uint32_t(jg == true_nr[1] - one_if_absorbing_otherwise_0 - 1) << 4) + + (uint32_t(kg == true_nr[2] - one_if_absorbing_otherwise_0 - 1) << 5); + // update the interior field values + if (val == 0) { + SourceVector_t interior = -anm1view(i, j, k) + a1 * aview(i, j, k) + + a2 * (aview(i + 1, j, k) + aview(i - 1, j, k)) + + a4 * (aview(i, j + 1, k) + aview(i, j - 1, k)) + + a6 * (aview(i, j, k + 1) + aview(i, j, k - 1)) + + a8 * source_view(i, j, k); + anp1view(i, j, k) = interior; + } + }); + Kokkos::fence(); + this->A_np1.fillHalo(); + this->applyBCs(); + } + + /** + * @brief Initializes the solver. + * + * This function initializes the solver by setting up the layout, mesh, and fields. + * It retrieves the mesh spacing, domain, and mesh size, and initializes the fields to zero. + */ + template + void StandardFDTDSolver::initialize() { + // get layout and mesh + this->layout_mp = &(this->JN_mp->getLayout()); + this->mesh_mp = &(this->JN_mp->get_mesh()); + + // get mesh spacing, timestep, domain, and mesh size + this->hr_m = this->mesh_mp->getMeshSpacing(); + this->dt = this->hr_m[0] / 2; + for (unsigned int i = 0; i < Dim; ++i) { + this->dt = std::min(this->dt, this->hr_m[i] / 2); + } + this->domain_m = this->layout_mp->getDomain(); + for (unsigned int i = 0; i < Dim; ++i) + this->nr_m[i] = this->domain_m[i].length(); + + // initialize fields + this->A_nm1.initialize(*(this->mesh_mp), *(this->layout_mp)); + this->A_n.initialize(*(this->mesh_mp), *(this->layout_mp)); + this->A_np1.initialize(*(this->mesh_mp), *(this->layout_mp)); + + // Initialize fields to zero + this->A_nm1 = 0.0; + this->A_n = 0.0; + this->A_np1 = 0.0; + + // set the mesh domain + if constexpr (boundary_conditions == periodic) { + this->setPeriodicBoundaryConditions(); + } + } +} // namespace ippl \ No newline at end of file diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 51b1eb206..dd4b72791 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -20,6 +20,7 @@ endif() if (IPPL_ENABLE_SOLVERS) add_subdirectory(solver) add_subdirectory(solver/fem) + add_subdirectory(maxwell) endif() add_subdirectory(particle) diff --git a/test/cmake/AddIpplIntegrationTest.cmake b/test/cmake/AddIpplIntegrationTest.cmake index f12592d4c..dd204cc30 100644 --- a/test/cmake/AddIpplIntegrationTest.cmake +++ b/test/cmake/AddIpplIntegrationTest.cmake @@ -8,7 +8,7 @@ function(add_ippl_integration_test TEST_NAME) set(options) - set(oneValueArgs COMMAND) + set(oneValueArgs COMMAND LINK_DIRS) set(multiValueArgs LABELS ARGS) cmake_parse_arguments(TEST "${options}" "${oneValueArgs}" "${multiValueArgs}" ${ARGN}) @@ -18,9 +18,16 @@ function(add_ippl_integration_test TEST_NAME) PRIVATE ippl ${MPI_CXX_LIBRARIES} ) - target_include_directories(${TEST_NAME} - PRIVATE $ - ) + if(TEST_LINK_DIRS) + target_include_directories(${TEST_NAME} + PRIVATE $ + ${TEST_LINK_DIRS} + ) + else() + target_include_directories(${TEST_NAME} + PRIVATE $ + ) + endif() if(TEST_COMMAND) add_test(NAME Integration.${TEST_NAME} COMMAND ${TEST_COMMAND}) diff --git a/test/maxwell/CMakeLists.txt b/test/maxwell/CMakeLists.txt new file mode 100644 index 000000000..e368d649c --- /dev/null +++ b/test/maxwell/CMakeLists.txt @@ -0,0 +1,15 @@ +# ----------------------------------------------------------------------------- +# test/maxwell/CMakeLists.txt +# +# Integration tests for IPPL maxwell solvers. +# ----------------------------------------------------------------------------- + +file(RELATIVE_PATH _relPath "${CMAKE_SOURCE_DIR}" "${CMAKE_CURRENT_SOURCE_DIR}") +message(STATUS "🔧 Adding maxwell integration tests from: ${_relPath}") + +if(IPPL_ENABLE_FFT) + add_ippl_integration_test(TestStandardFDTDSolver LINK_DIRS "${DOWNLOADED_HEADERS_DIR}" LABELS solver integration) + add_ippl_integration_test(TestNonStandardFDTDSolver LINK_DIRS "${DOWNLOADED_HEADERS_DIR}" LABELS solver integration) + add_ippl_integration_test(TestStandardFDTDSolver_convergence LINK_DIRS LABELS solver integration) + add_ippl_integration_test(TestNonStandardFDTDSolver_convergence LINK_DIRS LABELS solver integration) +endif() \ No newline at end of file diff --git a/test/maxwell/README.md b/test/maxwell/README.md new file mode 100644 index 000000000..9b7299fcc --- /dev/null +++ b/test/maxwell/README.md @@ -0,0 +1,92 @@ +## FDTD Solver Testcases + +This directory contains testcases for both the Standard and Non-Standard FDTD solvers. There are two types of tests: + +- **Image Output Testcases**: Visualize the field evolution of a gaussian pulse over a time increment of one. With periodic boundary conditions the wave is expected to return to the initial state. +- **Convergence Testcases**: Compute and record the convergence error for different grid sizes and directions, and plot the results. The L2 error between the initial field and the field after a time of one is calculated. + +--- + +### 1. Running the Image Output Testcases + +1. **Build the executables** (in existing IPPL build directory): + ```bash + cd ippl/build/test/solver/ + make TestStandardFDTDSolver + make TestNonStandardFDTDSolver + ``` + +2. **Create the required output folders**: + ```bash + mkdir -p renderdataStandard + mkdir -p renderdataNonStandard + ``` + +3. **Run the testcases**: + ```bash + ./TestStandardFDTDSolver + ./TestNonStandardFDTDSolver + ``` + +4. **Image Output**: + Every 4th timestep an images of the filed will be written to: + - `renderdataStandard/outimageXXXXX.bmp` (Standard solver) + - `renderdataNonStandard/outimageXXXXX.bmp` (Non-Standard solver) + +--- + +### 2. Running the Convergence Testcases + +1. **Build the executables**: + ```bash + make TestStandardFDTDSolver_convergence + make TestNonStandardFDTDSolver_convergence + ``` + +2. **Run the testcases**: + ```bash + ./TestStandardFDTDSolver_convergence + ./TestNonStandardFDTDSolver_convergence + ``` + +3. **Output**: + This will generate CSV files: + - `StandardFDTDSolver_convergence.csv` + - `NonStandardFDTDSolver_convergence.csv` + +### 3. Plotting the Convergence Results + +The following python code can be used to create the convergance plot: + +```bash +import pandas as pd +import matplotlib.pyplot as plt + +# Read the data from the CSV file +filenames = ["StandardFDTDSolver_convergence.csv", "NonStandardFDTDSolver_convergence.csv"] + +for filename in filenames: + df = pd.read_csv(filename) + + solvername = filename.split("FDTD")[0] + + # Plot figure + plt.figure(figsize=(8,6)) + + plt.plot(df.loc[df["GaussianPulseDir"] == 'x']["NGridpoints"], df.loc[df["GaussianPulseDir"] == 'x']["ConverganceError"], 'b-', label='L_2 error for x wave') + plt.plot(df.loc[df["GaussianPulseDir"] == 'y']["NGridpoints"], df.loc[df["GaussianPulseDir"] == 'y']["ConverganceError"], 'c--', label='L_2 error for y wave') + plt.plot(df.loc[df["GaussianPulseDir"] == 'z']["NGridpoints"], df.loc[df["GaussianPulseDir"] == 'z']["ConverganceError"], 'm:', label='L_2 error for z wave') + plt.plot(df.loc[df["GaussianPulseDir"] == 'z']["NGridpoints"], df.loc[df["GaussianPulseDir"] == 'z']["NGridpoints"].astype(float)**(-2) * 2*max(df["ConverganceError"]), 'k-', linewidth=0.5, label=r'O($\Delta x^{2}$)') + plt.plot(df.loc[df["GaussianPulseDir"] == 'z']["NGridpoints"], df.loc[df["GaussianPulseDir"] == 'z']["NGridpoints"].astype(float)**(-1) * 2*max(df["ConverganceError"]), 'k--', linewidth=0.5, label=r'O($\Delta x$)') + plt.yscale('log', base=2) + #plt.xscale('log', base=2) + plt.xlabel('Number of Gridpoints') + plt.ylabel('Error') + plt.title(f'{solvername} Finite Differences Solver Error') + plt.grid(True, which="both", linestyle="--", linewidth=0.5) + plt.legend() + plt.tight_layout() + plt.savefig(f'{solvername}FDTDSolver_convergence') +``` + +This will generate the convergence plots for both solvers and save them as image files in the current directory. diff --git a/test/maxwell/TestNonStandardFDTDSolver.cpp b/test/maxwell/TestNonStandardFDTDSolver.cpp new file mode 100644 index 000000000..ccb9f413e --- /dev/null +++ b/test/maxwell/TestNonStandardFDTDSolver.cpp @@ -0,0 +1,218 @@ +// TestNonStandardFDTDSolver +// Check the README.md file in this directory for information about the test and how to run it. + +#include +using std::size_t; +#include +#include "Ippl.h" + +#include "Types/Vector.h" + +#include "Field/Field.h" + +#include "MaxwellSolvers/NonStandardFDTDSolver.h" +#define STB_IMAGE_WRITE_IMPLEMENTATION +#include +template + requires((std::is_floating_point_v)) +KOKKOS_INLINE_FUNCTION float gauss(scalar1 mean, scalar1 stddev, scalar... x) { + uint32_t dim = sizeof...(scalar); + ippl::Vector vec{scalar1(x - mean)...}; + scalar1 vecsum(0); + for (unsigned d = 0; d < dim; d++) { + vecsum += vec[d] * vec[d]; + } + return Kokkos::exp(-(vecsum) / (stddev * stddev)); +} + +int main(int argc, char* argv[]) { + ippl::initialize(argc, argv); + { + using scalar = float; + const unsigned dim = 3; + using vector_type = ippl::Vector; + using vector4_type = ippl::Vector; + + using SourceField = + ippl::Field, + typename ippl::UniformCartesian::DefaultCentering>; + using EMField = ippl::Field, + typename ippl::UniformCartesian::DefaultCentering>; + + // Direction of the Gaussian pulse (can be 'x', 'y', or 'z') + char direction = 'z'; + + // Get variable for direction (1 for x, 2 for y and 3 for z) + const int dir = (direction == 'x') ? 1 : (direction == 'y') ? 2 : 3; + + // Specifie number of gridpoints in each direction (more gridpoints in z needed for CFL + // condition) + constexpr size_t n = 100; + ippl::Vector nr{n / 2, n / 2, 2 * n}; + ippl::NDIndex<3> owned(nr[0], nr[1], nr[2]); + + // specifies decomposition; + std::array isParallel; + isParallel.fill(true); + + // unit box + ippl::Vector extents{scalar(1), scalar(1), scalar(1)}; + ippl::Vector hx; + for (unsigned d = 0; d < 3; d++) { + hx[d] = extents[d] / (scalar)nr[d]; + } + ippl::Vector origin = {0, 0, 0}; + ippl::UniformCartesian mesh(owned, hx, origin); + + // all parallel layout, standard domain, normal axis order + ippl::FieldLayout<3> layout(MPI_COMM_WORLD, owned, isParallel); + + // Define the source and field types + SourceField source(mesh, layout); + EMField E(mesh, layout); + EMField B(mesh, layout); + source = vector4_type(0); + + // Create the NonStandardFDTDSolver object + ippl::NonStandardFDTDSolver sfdsolver(source, E, B); + + // Initialize the source field with a Gaussian distribution in the z-direction and zeros in + // the x and y directions + auto aview = sfdsolver.A_n.getView(); + auto am1view = sfdsolver.A_nm1.getView(); + auto ap1view = sfdsolver.A_np1.getView(); + auto lDom = layout.getLocalNDIndex(); + size_t nghost = sfdsolver.A_n.getNghost(); + + // Initialize the fields and calculate the sum of the squared magnitudes for error + // calculation later + double sum_norm = 0.0; + Kokkos::parallel_reduce( + ippl::getRangePolicy(aview, 1), + KOKKOS_LAMBDA(size_t i, size_t j, size_t k, double& ref) { + const size_t ig = i + lDom.first()[0] - nghost; + const size_t jg = j + lDom.first()[1] - nghost; + const size_t kg = k + lDom.first()[2] - nghost; + scalar x = scalar(ig) / nr[0]; + scalar y = scalar(jg) / nr[1]; + scalar z = scalar(kg) / nr[2]; + (void)x; + (void)y; + (void)z; + + // Calculate gaussian pules in direction dir + const scalar coord = (dir == 1) ? x : (dir == 2) ? y : z; + const scalar magnitude = gauss(scalar(0.5), scalar(0.05), coord); + + // Initialize fields + aview(i, j, k) = vector4_type{scalar(0), scalar(0), scalar(0), scalar(0)}; + aview(i, j, k)[dir] = magnitude; + + am1view(i, j, k) = vector4_type{scalar(0), scalar(0), scalar(0), scalar(0)}; + am1view(i, j, k)[dir] = magnitude; + + ref += magnitude + * magnitude; // Accumulate the squared magnitude for error calculation + }, + sum_norm); + Kokkos::fence(); + + // Apply the boundary conditions to the initialized fields + sfdsolver.A_n.getFieldBC().apply(sfdsolver.A_n); + sfdsolver.A_np1.getFieldBC().apply(sfdsolver.A_np1); + sfdsolver.A_nm1.getFieldBC().apply(sfdsolver.A_nm1); + sfdsolver.A_n.fillHalo(); + sfdsolver.A_nm1.fillHalo(); + + // Create a 500x500 image with 3 channels (RGB) + int img_width = 500; + int img_height = 500; + float* imagedata = new float[img_width * img_height * 3]; + + // Initialize the image data to black + std::fill(imagedata, imagedata + img_width * img_height * 3, 0.0f); + uint8_t* imagedata_final = new uint8_t[img_width * img_height * 3]; + + // Run the simulation for 1s, with periodic boundary conditions this should be the same + // state as at time 0 + for (size_t s = 0; s < 1. / sfdsolver.getDt(); s++) { + auto ebh = sfdsolver.A_n.getHostMirror(); + Kokkos::deep_copy(ebh, sfdsolver.A_n.getView()); + for (int i = 1; i < img_width; i++) { + for (int j = 1; j < img_height; j++) { + // Map the pixel coordinates to the simulation domain + int i_remap = (double(i) / (img_width - 1)) * (nr[2] - 4) + 2; + int j_remap = (double(j) / (img_height - 1)) * (nr[0] - 4) + 2; + + // Check if the remapped coordinates are within the local domain + if (i_remap >= lDom.first()[2] && i_remap <= lDom.last()[2]) { + if (j_remap >= lDom.first()[0] && j_remap <= lDom.last()[0]) { + // Get the corresponding field vector + ippl::Vector acc = + ebh(j_remap + 1 - lDom.first()[0], nr[1] / 2, + i_remap + 1 - lDom.first()[2]); + + // Normalize field vector and set pixel color + float normalized_colorscale_value = acc.Pnorm(); + imagedata[(j * img_width + i) * 3 + 0] = + normalized_colorscale_value * 255.0f; + imagedata[(j * img_width + i) * 3 + 1] = + normalized_colorscale_value * 255.0f; + imagedata[(j * img_width + i) * 3 + 2] = + normalized_colorscale_value * 255.0f; + } + } + } + } + + // Convert the float image data to unsigned char + std::transform(imagedata, imagedata + img_height * img_width * 3, imagedata_final, + [](float x) { + return (unsigned char)std::min(255.0f, std::max(0.0f, x)); + }); + char output[1024] = {0}; + snprintf(output, 1023, "%soutimage%.05lu.bmp", "renderdataNonStandard/", s); + + // Save the image every 4th step + if (s % 4 == 0) + stbi_write_bmp(output, img_width, img_height, 3, imagedata_final); + + // Solve the FDTD equations + sfdsolver.solve(); + } + + // Calculate the error between the computed and expected values + double sum_error = 0.0; + Kokkos::parallel_reduce( + ippl::getRangePolicy(aview, 1), + KOKKOS_LAMBDA(size_t i, size_t j, size_t k, double& ref) { + const size_t ig = i + lDom.first()[0] - nghost; + const size_t jg = j + lDom.first()[1] - nghost; + const size_t kg = k + lDom.first()[2] - nghost; + scalar x = scalar(ig) / nr[0]; + scalar y = scalar(jg) / nr[1]; + scalar z = scalar(kg) / nr[2]; + (void)x; + (void)y; + (void)z; + + // Get coordinate in the right direction + const scalar coord = (dir == 1) ? x : (dir == 2) ? y : z; + + // Calculate error in given direction at this point + double original_value = gauss(scalar(0.5), scalar(0.05), coord); + + // Calculate the difference between the computed and expected values + ippl::Vector diff = aview(i, j, k); + diff[dir] -= original_value; + + // Accumulate the squared differences for L2 norm and the original value for + // normalization + ref += ippl::dot(diff, diff).apply(); + }, + sum_error); + std::cout << "Sum after evaluating boundaryconditions: " + << Kokkos::sqrt(sum_error / sum_norm) << "\n"; + } + ippl::finalize(); +} \ No newline at end of file diff --git a/test/maxwell/TestNonStandardFDTDSolver_convergence.cpp b/test/maxwell/TestNonStandardFDTDSolver_convergence.cpp new file mode 100644 index 000000000..ae9305037 --- /dev/null +++ b/test/maxwell/TestNonStandardFDTDSolver_convergence.cpp @@ -0,0 +1,190 @@ +// TestNonStandardFDTDSolver_convergence +// Check the README.md file in this directory for information about the test and how to run it. + +#include +using std::size_t; +#include +#include "Ippl.h" + +#include "Types/Vector.h" + +#include "Field/Field.h" + +#include "MaxwellSolvers/NonStandardFDTDSolver.h" + +template + requires((std::is_floating_point_v)) +KOKKOS_INLINE_FUNCTION float gauss(scalar1 mean, scalar1 stddev, scalar... x) { + uint32_t dim = sizeof...(scalar); + ippl::Vector vec{scalar1(x - mean)...}; + scalar1 vecsum(0); + for (unsigned d = 0; d < dim; d++) { + vecsum += vec[d] * vec[d]; + } + return Kokkos::exp(-(vecsum) / (stddev * stddev)); +} + +void compute_convergence(char direction, unsigned int np, std::string fname) { + using scalar = float; + const unsigned dim = 3; + using vector_type = ippl::Vector; + using vector4_type = ippl::Vector; + + using SourceField = ippl::Field, + typename ippl::UniformCartesian::DefaultCentering>; + using EMField = ippl::Field, + typename ippl::UniformCartesian::DefaultCentering>; + + // Get variable for direction (1 for x, 2 for y and 3 for z) + const int dir = (direction == 'x') ? 1 : (direction == 'y') ? 2 : 3; + + // Specifie number of gridpoints in each direction (more gridpoints in z needed for CFL + // condition) + ippl::Vector nr{np / 2, np / 2, 2 * np}; + ippl::NDIndex<3> owned(nr[0], nr[1], nr[2]); + + // specifies decomposition; + std::array isParallel; + isParallel.fill(true); + + // unit box + ippl::Vector hx; + for (unsigned d = 0; d < 3; d++) { + hx[d] = 1 / (scalar)nr[d]; + } + ippl::Vector origin = {0.0, 0.0, 0.0}; + ippl::UniformCartesian mesh(owned, hx, origin); + + // all parallel layout, standard domain, normal axis order + ippl::FieldLayout<3> layout(MPI_COMM_WORLD, owned, isParallel); + + // Define the source and field types + SourceField source(mesh, layout); + EMField E(mesh, layout); + EMField B(mesh, layout); + source = vector4_type(0); + + // Create the NonStandardFDTDSolver object + ippl::NonStandardFDTDSolver sfdsolver(source, E, B); + + // Initialize the source field with a Gaussian distribution in the given direction and zeros in + // the other directions + auto aview = sfdsolver.A_n.getView(); + auto am1view = sfdsolver.A_nm1.getView(); + auto ap1view = sfdsolver.A_np1.getView(); + auto lDom = layout.getLocalNDIndex(); + size_t nghost = sfdsolver.A_n.getNghost(); + + // Initialize the fields and calculate the sum of the squared magnitudes for error calculation + // later + double sum_norm = 0.0; + Kokkos::parallel_reduce( + ippl::getRangePolicy(aview, 1), + KOKKOS_LAMBDA(size_t i, size_t j, size_t k, double& ref) { + // Calculate position in x, y and z + const size_t ig = i + lDom.first()[0] - nghost; + const size_t jg = j + lDom.first()[1] - nghost; + const size_t kg = k + lDom.first()[2] - nghost; + scalar x = scalar(ig) / nr[0]; + scalar y = scalar(jg) / nr[1]; + scalar z = scalar(kg) / nr[2]; + (void)x; + (void)y; + (void)z; + + // Calculate gaussian pules in direction dir + const scalar coord = (dir == 1) ? x : (dir == 2) ? y : z; + const scalar magnitude = gauss(scalar(0.5), scalar(0.05), coord); + + // Initialize fields + aview(i, j, k) = vector4_type{scalar(0), scalar(0), scalar(0), scalar(0)}; + aview(i, j, k)[dir] = magnitude; + + am1view(i, j, k) = vector4_type{scalar(0), scalar(0), scalar(0), scalar(0)}; + am1view(i, j, k)[dir] = magnitude; + + // Calculate the sum of the squared magnitudes for error calculation later + ref += magnitude * magnitude; + }, + sum_norm); + Kokkos::fence(); + + // Apply the boundary conditions to the initialized fields + sfdsolver.A_n.getFieldBC().apply(sfdsolver.A_n); + sfdsolver.A_np1.getFieldBC().apply(sfdsolver.A_np1); + sfdsolver.A_nm1.getFieldBC().apply(sfdsolver.A_nm1); + sfdsolver.A_n.fillHalo(); + sfdsolver.A_nm1.fillHalo(); + + // Run the simulation for 1s, with periodic boundary conditions this should be the same state as + // at time 0 + for (size_t s = 0; s < 1. / sfdsolver.getDt(); s++) { + // Solve the FDTD equations + sfdsolver.solve(); + } + + // Calculate the L2 norm error between the computed and expected values + double sum_error = 0.0; + Kokkos::parallel_reduce( + ippl::getRangePolicy(aview, 1), + KOKKOS_LAMBDA(size_t i, size_t j, size_t k, double& ref) { + // Calculate position in x, y and z + const size_t ig = i + lDom.first()[0] - nghost; + const size_t jg = j + lDom.first()[1] - nghost; + const size_t kg = k + lDom.first()[2] - nghost; + scalar x = scalar(ig) / nr[0]; + scalar y = scalar(jg) / nr[1]; + scalar z = scalar(kg) / nr[2]; + (void)x; + (void)y; + (void)z; + + // Get coordinate in the right direction + const scalar coord = (dir == 1) ? x : (dir == 2) ? y : z; + + // Calculate error in given direction at this point + double original_value = gauss(scalar(0.5), scalar(0.05), coord); + + // Calculate the difference between the computed and expected values + ippl::Vector diff = aview(i, j, k); + diff[dir] -= original_value; + + // Accumulate the squared differences for L2 norm and the original value for + // normalization + ref += ippl::dot(diff, diff).apply(); + }, + sum_error); + + // Write results to output file + Inform csvout(NULL, fname.c_str(), Inform::APPEND); + csvout.precision(16); + csvout.setf(std::ios::scientific, std::ios::floatfield); + + csvout << direction << "," << np << "," << Kokkos::sqrt(sum_error / sum_norm) << endl; +} + +int main(int argc, char* argv[]) { + ippl::initialize(argc, argv); + { + // gridsizes to iterate over + std::array N = {4, 8, 16, 32, 64, 128, 256}; + + // directions to iterate over + std::array directions = {'x', 'y', 'z'}; + + // Create outputfile and plot header + std::string fname = "NonStandardFDTDSolver_convergence.csv"; + + Inform csvout(NULL, fname.c_str(), Inform::OVERWRITE); + csvout.precision(16); + csvout.setf(std::ios::scientific, std::ios::floatfield); + csvout << "GaussianPulseDir,NGridpoints,ConverganceError" << endl; + + for (char dir : directions) { + for (int pt : N) { + compute_convergence(dir, pt, fname); + } + } + } + ippl::finalize(); +} \ No newline at end of file diff --git a/test/maxwell/TestStandardFDTDSolver.cpp b/test/maxwell/TestStandardFDTDSolver.cpp new file mode 100644 index 000000000..1b85f50e9 --- /dev/null +++ b/test/maxwell/TestStandardFDTDSolver.cpp @@ -0,0 +1,217 @@ +// TestStandardFDTDSolver +// Check the README.md file in this directory for information about the test and how to run it. + +#include +using std::size_t; +#include +#include "Ippl.h" + +#include "Types/Vector.h" + +#include "Field/Field.h" + +#include "MaxwellSolvers/StandardFDTDSolver.h" +#define STB_IMAGE_WRITE_IMPLEMENTATION +#include +template + requires((std::is_floating_point_v)) +KOKKOS_INLINE_FUNCTION float gauss(scalar1 mean, scalar1 stddev, scalar... x) { + uint32_t dim = sizeof...(scalar); + ippl::Vector vec{scalar1(x - mean)...}; + scalar1 vecsum(0); + for (unsigned d = 0; d < dim; d++) { + vecsum += vec[d] * vec[d]; + } + return Kokkos::exp(-(vecsum) / (stddev * stddev)); +} + +int main(int argc, char* argv[]) { + ippl::initialize(argc, argv); + { + using scalar = float; + const unsigned dim = 3; + using vector_type = ippl::Vector; + using vector4_type = ippl::Vector; + + using SourceField = + ippl::Field, + typename ippl::UniformCartesian::DefaultCentering>; + using EMField = ippl::Field, + typename ippl::UniformCartesian::DefaultCentering>; + + // Direction of the Gaussian pulse (can be 'x', 'y', or 'z') + char direction = 'z'; + + // Get variable for direction (1 for x, 2 for y and 3 for z) + const int dir = (direction == 'x') ? 1 : (direction == 'y') ? 2 : 3; + + // Specifie number of gridpoints in each direction (more gridpoints in z needed for CFL + // condition) + constexpr size_t n = 100; + ippl::Vector nr{n, n, n}; + ippl::NDIndex<3> owned(nr[0], nr[1], nr[2]); + + // specifies decomposition; + std::array isParallel; + isParallel.fill(true); + + // unit box + ippl::Vector extents{scalar(1), scalar(1), scalar(1)}; + ippl::Vector hx; + for (unsigned d = 0; d < 3; d++) { + hx[d] = extents[d] / (scalar)nr[d]; + } + ippl::Vector origin = {0, 0, 0}; + ippl::UniformCartesian mesh(owned, hx, origin); + + // all parallel layout, standard domain, normal axis order + ippl::FieldLayout<3> layout(MPI_COMM_WORLD, owned, isParallel); + + // Define the source and field types + SourceField source(mesh, layout); + EMField E(mesh, layout); + EMField B(mesh, layout); + source = vector4_type(0); + + // Create the StandardFDTDSolver object + ippl::StandardFDTDSolver sfdsolver(source, E, B); + + // Initialize the source field with a Gaussian distribution in the z-direction and zeros in + // the x and y directions + auto aview = sfdsolver.A_n.getView(); + auto am1view = sfdsolver.A_nm1.getView(); + auto ap1view = sfdsolver.A_np1.getView(); + auto lDom = layout.getLocalNDIndex(); + size_t nghost = sfdsolver.A_n.getNghost(); + + // Initialize the fields and calculate the sum of the squared magnitudes for error + // calculation later + double sum_norm = 0.0; + Kokkos::parallel_reduce( + ippl::getRangePolicy(aview, 1), + KOKKOS_LAMBDA(size_t i, size_t j, size_t k, double& ref) { + const size_t ig = i + lDom.first()[0] - nghost; + const size_t jg = j + lDom.first()[1] - nghost; + const size_t kg = k + lDom.first()[2] - nghost; + scalar x = scalar(ig) / nr[0]; + scalar y = scalar(jg) / nr[1]; + scalar z = scalar(kg) / nr[2]; + (void)x; + (void)y; + (void)z; + + // Calculate gaussian pules in direction dir + const scalar coord = (dir == 1) ? x : (dir == 2) ? y : z; + const scalar magnitude = gauss(scalar(0.5), scalar(0.05), coord); + + // Initialize fields + aview(i, j, k) = vector4_type{scalar(0), scalar(0), scalar(0), scalar(0)}; + aview(i, j, k)[dir] = magnitude; + + am1view(i, j, k) = vector4_type{scalar(0), scalar(0), scalar(0), scalar(0)}; + am1view(i, j, k)[dir] = magnitude; + + // Accumulate the squared magnitude for error calculation later + ref += magnitude * magnitude; + }, + sum_norm); + Kokkos::fence(); + + // Apply the boundary conditions to the initialized fields + sfdsolver.A_n.getFieldBC().apply(sfdsolver.A_n); + sfdsolver.A_np1.getFieldBC().apply(sfdsolver.A_np1); + sfdsolver.A_nm1.getFieldBC().apply(sfdsolver.A_nm1); + sfdsolver.A_n.fillHalo(); + sfdsolver.A_nm1.fillHalo(); + + // Create a 500x500 image with 3 channels (RGB) + int img_width = 500; + int img_height = 500; + float* imagedata = new float[img_width * img_height * 3]; + + // Initialize the image data to black + std::fill(imagedata, imagedata + img_width * img_height * 3, 0.0f); + uint8_t* imagedata_final = new uint8_t[img_width * img_height * 3]; + + // Run the simulation for a duration of 1, with periodic boundary conditions the wave will + // return to the initial state + for (size_t s = 0; s < 1. / sfdsolver.getDt(); s++) { + auto ebh = sfdsolver.A_n.getHostMirror(); + Kokkos::deep_copy(ebh, sfdsolver.A_n.getView()); + for (int i = 1; i < img_width; i++) { + for (int j = 1; j < img_height; j++) { + // Map the pixel coordinates to the simulation domain + int i_remap = (double(i) / (img_width - 1)) * (nr[2] - 4) + 2; + int j_remap = (double(j) / (img_height - 1)) * (nr[0] - 4) + 2; + + // Check if the remapped coordinates are within the local domain + if (i_remap >= lDom.first()[2] && i_remap <= lDom.last()[2]) { + if (j_remap >= lDom.first()[0] && j_remap <= lDom.last()[0]) { + // Get the corresponding field vector + ippl::Vector acc = + ebh(j_remap + 1 - lDom.first()[0], nr[1] / 2, + i_remap + 1 - lDom.first()[2]); + + // Normalize field vector and set pixel color + float normalized_colorscale_value = acc.Pnorm(); + imagedata[(j * img_width + i) * 3 + 0] = + normalized_colorscale_value * 255.0f; + imagedata[(j * img_width + i) * 3 + 1] = + normalized_colorscale_value * 255.0f; + imagedata[(j * img_width + i) * 3 + 2] = + normalized_colorscale_value * 255.0f; + } + } + } + } + + // Convert the float image data to unsigned char + std::transform(imagedata, imagedata + img_height * img_width * 3, imagedata_final, + [](float x) { + return (unsigned char)std::min(255.0f, std::max(0.0f, x)); + }); + char output[1024] = {0}; + snprintf(output, 1023, "%soutimage%.05lu.bmp", "renderdataStandard/", s); + + // Save the image every 4th step + if (s % 4 == 0) + stbi_write_bmp(output, img_width, img_height, 3, imagedata_final); + + // Solve the FDTD equations + sfdsolver.solve(); + } + + // Calculate the error between the computed and expected values + double sum_error = 0.0; + Kokkos::parallel_reduce( + ippl::getRangePolicy(aview, 1), + KOKKOS_LAMBDA(size_t i, size_t j, size_t k, double& ref) { + const size_t ig = i + lDom.first()[0] - nghost; + const size_t jg = j + lDom.first()[1] - nghost; + const size_t kg = k + lDom.first()[2] - nghost; + scalar x = scalar(ig) / nr[0]; + scalar y = scalar(jg) / nr[1]; + scalar z = scalar(kg) / nr[2]; + (void)x; + (void)y; + (void)z; + // Get coordinate in the right direction + const scalar coord = (dir == 1) ? x : (dir == 2) ? y : z; + + // Calculate error in given direction at this point + double original_value = gauss(scalar(0.5), scalar(0.05), coord); + + // Calculate the difference between the computed and expected values + ippl::Vector diff = aview(i, j, k); + diff[dir] -= original_value; + + // Accumulate the squared differences for L2 norm and the original value for + // normalization + ref += ippl::dot(diff, diff).apply(); + }, + sum_error); + std::cout << "Sum after evaluating boundaryconditions: " + << Kokkos::sqrt(sum_error / sum_norm) << "\n"; + } + ippl::finalize(); +} \ No newline at end of file diff --git a/test/maxwell/TestStandardFDTDSolver_convergence.cpp b/test/maxwell/TestStandardFDTDSolver_convergence.cpp new file mode 100644 index 000000000..e4a259758 --- /dev/null +++ b/test/maxwell/TestStandardFDTDSolver_convergence.cpp @@ -0,0 +1,190 @@ +// TestStandardFDTDSolver_convergence +// Check the README.md file in this directory for information about the test and how to run it. + +#include +using std::size_t; +#include +#include "Ippl.h" + +#include "Types/Vector.h" + +#include "Field/Field.h" + +#include "MaxwellSolvers/StandardFDTDSolver.h" + +template + requires((std::is_floating_point_v)) +KOKKOS_INLINE_FUNCTION float gauss(scalar1 mean, scalar1 stddev, scalar... x) { + uint32_t dim = sizeof...(scalar); + ippl::Vector vec{scalar1(x - mean)...}; + scalar1 vecsum(0); + for (unsigned d = 0; d < dim; d++) { + vecsum += vec[d] * vec[d]; + } + return Kokkos::exp(-(vecsum) / (stddev * stddev)); +} + +void compute_convergence(char direction, unsigned int np, std::string fname) { + using scalar = float; + const unsigned dim = 3; + using vector_type = ippl::Vector; + using vector4_type = ippl::Vector; + + using SourceField = ippl::Field, + typename ippl::UniformCartesian::DefaultCentering>; + using EMField = ippl::Field, + typename ippl::UniformCartesian::DefaultCentering>; + + // Get variable for direction (0 for x, 1 for y and 2 for z) + const int dir = (direction == 'x') ? 1 : (direction == 'y') ? 2 : 3; + + // Specifie number of gridpoints in each direction (more gridpoints in z needed for CFL + // condition) + ippl::Vector nr{np, np, np}; + ippl::NDIndex<3> owned(nr[0], nr[1], nr[2]); + + // specifies decomposition; here all dimensions are parallel + std::array isParallel; + isParallel.fill(true); + + // unit box + ippl::Vector hx; + for (unsigned d = 0; d < 3; d++) { + hx[d] = 1 / (scalar)nr[d]; + } + ippl::Vector origin = {0.0, 0.0, 0.0}; + ippl::UniformCartesian mesh(owned, hx, origin); + + // all parallel layout, standard domain, normal axis order + ippl::FieldLayout<3> layout(MPI_COMM_WORLD, owned, isParallel); + + // Define the source and field types + SourceField source(mesh, layout); + EMField E(mesh, layout); + EMField B(mesh, layout); + source = vector4_type(0); + + // Create the StandardFDTDSolver object + ippl::StandardFDTDSolver sfdsolver(source, E, B); + + // Initialize the source field with a Gaussian distribution in the given direction and zeros in + // the other directions + auto aview = sfdsolver.A_n.getView(); + auto am1view = sfdsolver.A_nm1.getView(); + auto ap1view = sfdsolver.A_np1.getView(); + auto lDom = layout.getLocalNDIndex(); + size_t nghost = sfdsolver.A_n.getNghost(); + + // Initialize the fields and calculate the sum of the squared magnitudes for error calculation + // later + double sum_norm = 0.0; + Kokkos::parallel_reduce( + ippl::getRangePolicy(aview, 1), + KOKKOS_LAMBDA(size_t i, size_t j, size_t k, double& ref) { + // Calculate position in x, y and z + const size_t ig = i + lDom.first()[0] - nghost; + const size_t jg = j + lDom.first()[1] - nghost; + const size_t kg = k + lDom.first()[2] - nghost; + scalar x = scalar(ig) / nr[0]; + scalar y = scalar(jg) / nr[1]; + scalar z = scalar(kg) / nr[2]; + (void)x; + (void)y; + (void)z; + + // Calculate gaussian pules in direction dir + const scalar coord = (dir == 1) ? x : (dir == 2) ? y : z; + const scalar magnitude = gauss(scalar(0.5), scalar(0.05), coord); + + // Initialize fields + aview(i, j, k) = vector4_type{scalar(0), scalar(0), scalar(0), scalar(0)}; + aview(i, j, k)[dir] = magnitude; + + am1view(i, j, k) = vector4_type{scalar(0), scalar(0), scalar(0), scalar(0)}; + am1view(i, j, k)[dir] = magnitude; + + // Calculate the sum of the squared magnitudes for L2 norm error calculation later + ref += magnitude * magnitude; + }, + sum_norm); + Kokkos::fence(); + + // Apply the boundary conditions to the initialized fields + sfdsolver.A_n.getFieldBC().apply(sfdsolver.A_n); + sfdsolver.A_np1.getFieldBC().apply(sfdsolver.A_np1); + sfdsolver.A_nm1.getFieldBC().apply(sfdsolver.A_nm1); + sfdsolver.A_n.fillHalo(); + sfdsolver.A_nm1.fillHalo(); + + // Run the simulation for 1s, with periodic boundary conditions this should be the same state as + // at time 0 + for (size_t s = 0; s < 1. / sfdsolver.getDt(); s++) { + // Solve the FDTD equations + sfdsolver.solve(); + } + + // Calculate the L2 norm error between the computed and expected values + double sum_error = 0.0; + Kokkos::parallel_reduce( + ippl::getRangePolicy(aview, 1), + KOKKOS_LAMBDA(size_t i, size_t j, size_t k, double& ref) { + // Calculate position in x, y and z + const size_t ig = i + lDom.first()[0] - nghost; + const size_t jg = j + lDom.first()[1] - nghost; + const size_t kg = k + lDom.first()[2] - nghost; + scalar x = scalar(ig) / nr[0]; + scalar y = scalar(jg) / nr[1]; + scalar z = scalar(kg) / nr[2]; + (void)x; + (void)y; + (void)z; + + // Get coordinate in the right direction + const scalar coord = (dir == 1) ? x : (dir == 2) ? y : z; + + // Calculate error in given direction at this point + double original_value = gauss(scalar(0.5), scalar(0.05), coord); + + // Calculate the difference between the computed and expected values + ippl::Vector diff = aview(i, j, k); + diff[dir] -= original_value; + + // Accumulate the squared differences for L2 norm and the original value for + // normalization + ref += ippl::dot(diff, diff).apply(); + }, + sum_error); + + // Write results to output file + Inform csvout(NULL, fname.c_str(), Inform::APPEND); + csvout.precision(16); + csvout.setf(std::ios::scientific, std::ios::floatfield); + + csvout << direction << "," << np << "," << Kokkos::sqrt(sum_error / sum_norm) << endl; +} + +int main(int argc, char* argv[]) { + ippl::initialize(argc, argv); + { + // gridsizes to iterate over + std::array N = {4, 8, 16, 32, 64, 128, 256}; + + // directions to iterate over + std::array directions = {'x', 'y', 'z'}; + + // Create outputfile and plot header + std::string fname = "StandardFDTDSolver_convergence.csv"; + + Inform csvout(NULL, fname.c_str(), Inform::OVERWRITE); + csvout.precision(16); + csvout.setf(std::ios::scientific, std::ios::floatfield); + csvout << "GaussianPulseDir,NGridpoints,ConverganceError" << endl; + + for (char dir : directions) { + for (int pt : N) { + compute_convergence(dir, pt, fname); + } + } + } + ippl::finalize(); +} \ No newline at end of file