From 03f7a60205d0f5518ef9afef35b69a6025eabeec Mon Sep 17 00:00:00 2001 From: Anna Date: Tue, 17 Jun 2025 15:37:41 +0200 Subject: [PATCH 1/6] Change CMake to download additional dependencie and add possibility to link tests with external files --- cmake/Dependencies.cmake | 9 +++++++++ src/MaxwellSolvers/CMakeLists.txt | 6 ++++++ test/cmake/AddIpplIntegrationTest.cmake | 15 +++++++++++---- test/solver/CMakeLists.txt | 4 ++++ 4 files changed, 30 insertions(+), 4 deletions(-) diff --git a/cmake/Dependencies.cmake b/cmake/Dependencies.cmake index 955d292e2..d800a402f 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 free electron laster.") +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/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/solver/CMakeLists.txt b/test/solver/CMakeLists.txt index d095eeb01..323d6467c 100644 --- a/test/solver/CMakeLists.txt +++ b/test/solver/CMakeLists.txt @@ -19,4 +19,8 @@ if(IPPL_ENABLE_FFT) add_ippl_integration_test(TestGaussian_biharmonic LABELS solver integration) add_ippl_integration_test(TestGaussian_hessian LABELS solver integration) add_ippl_integration_test(TestP3MSolver LABELS solver integration) + 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() From 5654d23d47e5cd2c0535a342eb7fdf83b90aea65 Mon Sep 17 00:00:00 2001 From: Anna Date: Tue, 17 Jun 2025 15:39:41 +0200 Subject: [PATCH 2/6] Add Standard and NonStandard FDTD Solver with Base class and Absorbing BC (Absorbing BC are not implemented very nicely, this could be improved, base code from Manuel) --- src/MaxwellSolvers/AbsorbingBC.h | 573 +++++++++++++++++++ src/MaxwellSolvers/FDTDSolverBase.h | 104 ++++ src/MaxwellSolvers/FDTDSolverBase.hpp | 136 +++++ src/MaxwellSolvers/NonStandardFDTDSolver.h | 78 +++ src/MaxwellSolvers/NonStandardFDTDSolver.hpp | 129 +++++ src/MaxwellSolvers/StandardFDTDSolver.h | 59 ++ src/MaxwellSolvers/StandardFDTDSolver.hpp | 113 ++++ 7 files changed, 1192 insertions(+) create mode 100644 src/MaxwellSolvers/AbsorbingBC.h create mode 100644 src/MaxwellSolvers/FDTDSolverBase.h create mode 100644 src/MaxwellSolvers/FDTDSolverBase.hpp create mode 100644 src/MaxwellSolvers/NonStandardFDTDSolver.h create mode 100644 src/MaxwellSolvers/NonStandardFDTDSolver.hpp create mode 100644 src/MaxwellSolvers/StandardFDTDSolver.h create mode 100644 src/MaxwellSolvers/StandardFDTDSolver.hpp diff --git a/src/MaxwellSolvers/AbsorbingBC.h b/src/MaxwellSolvers/AbsorbingBC.h new file mode 100644 index 000000000..3678e3f32 --- /dev/null +++ b/src/MaxwellSolvers/AbsorbingBC.h @@ -0,0 +1,573 @@ +#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...}; + 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}; + 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) + 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/FDTDSolverBase.h b/src/MaxwellSolvers/FDTDSolverBase.h new file mode 100644 index 000000000..0f479b25a --- /dev/null +++ b/src/MaxwellSolvers/FDTDSolverBase.h @@ -0,0 +1,104 @@ +#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. + + // TODO make this private again + /** + * @brief Shifts the time for the fields. + */ + void timeShift(); + protected: + /** + * @brief Applies the boundary conditions. + */ + void applyBCs(); + + public: + /** + * @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(); + + bool periodic_bc; // Flag indicating if periodic boundary conditions are used. + + 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. + + /** + * @brief Initializes the solver. This is a pure virtual function. + */ + virtual void initialize() = 0; + }; +} // 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..93386a1cc --- /dev/null +++ b/src/MaxwellSolvers/FDTDSolverBase.hpp @@ -0,0 +1,136 @@ +// +// 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() { + periodic_bc = true; + 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 time for the fields. + * + * 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() { + // TODO: should this comment stay? + // 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(); + + 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; + // TODO make this work with tail!!! + //ippl::Vector dAdt = + //(A_n(i, j, k).template tail<3>() - A_nm1(i, j, k).template tail<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..08ae1ed4f --- /dev/null +++ b/src/MaxwellSolvers/NonStandardFDTDSolver.h @@ -0,0 +1,78 @@ +/** + * @file NonStandardFDTDSolver.h + * @brief Defines the NonStandardFDTDSolver class for solving Maxwell's equations using a non-standard Finite-Difference Time-Domain (FDTD) method. + */ + +#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/Maxwell.h" +#include "MaxwellSolvers/FDTDSolverBase.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..88892d31a --- /dev/null +++ b/src/MaxwellSolvers/NonStandardFDTDSolver.hpp @@ -0,0 +1,129 @@ +/** + * @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) { + std::cerr << "Dispersion-free CFL condition not satisfiable\n"; + std::abort(); + } + 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; + }; +} \ No newline at end of file diff --git a/src/MaxwellSolvers/StandardFDTDSolver.h b/src/MaxwellSolvers/StandardFDTDSolver.h new file mode 100644 index 000000000..7a2f42389 --- /dev/null +++ b/src/MaxwellSolvers/StandardFDTDSolver.h @@ -0,0 +1,59 @@ +/** + * @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/Maxwell.h" +#include "MaxwellSolvers/FDTDSolverBase.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..b74b55841 --- /dev/null +++ b/src/MaxwellSolvers/StandardFDTDSolver.hpp @@ -0,0 +1,113 @@ +/** + * @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; + } +} \ No newline at end of file From d5b8783c0b2d4ab8904e65503b6b0a4f28835c0a Mon Sep 17 00:00:00 2001 From: Anna Date: Tue, 17 Jun 2025 15:42:12 +0200 Subject: [PATCH 3/6] Add tests for Standard and NonStandard FDTD solver and README.md with pythonscript to recreate convergence plot --- test/solver/README.md | 92 ++++++++++ test/solver/TestNonStandardFDTDSolver.cpp | 162 +++++++++++++++++ .../TestNonStandardFDTDSolver_convergence.cpp | 166 ++++++++++++++++++ test/solver/TestStandardFDTDSolver.cpp | 162 +++++++++++++++++ .../TestStandardFDTDSolver_convergence.cpp | 166 ++++++++++++++++++ 5 files changed, 748 insertions(+) create mode 100644 test/solver/README.md create mode 100644 test/solver/TestNonStandardFDTDSolver.cpp create mode 100644 test/solver/TestNonStandardFDTDSolver_convergence.cpp create mode 100644 test/solver/TestStandardFDTDSolver.cpp create mode 100644 test/solver/TestStandardFDTDSolver_convergence.cpp diff --git a/test/solver/README.md b/test/solver/README.md new file mode 100644 index 000000000..94b4cac77 --- /dev/null +++ b/test/solver/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/solver/TestNonStandardFDTDSolver.cpp b/test/solver/TestNonStandardFDTDSolver.cpp new file mode 100644 index 000000000..5aac6a7ad --- /dev/null +++ b/test/solver/TestNonStandardFDTDSolver.cpp @@ -0,0 +1,162 @@ +#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]; + + } + #ifndef __CUDA_ARCH__ + using std::exp; + #endif + return 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>; + constexpr size_t n = 100; + ippl::Vector nr{n / 2, n / 2, 2 * n}; + ippl::NDIndex<3> owned(nr[0], nr[1], nr[2]); + ippl::Vector extents{scalar(1), scalar(1), scalar(1)}; + std::array isParallel; + isParallel.fill(false); + isParallel[2] = true; + + // all parallel layout, standard domain, normal axis order + ippl::FieldLayout<3> layout(MPI_COMM_WORLD, owned, isParallel); + + 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); + + // 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); + + // Set boundary conditions + sfdsolver.setPeriodicBoundaryConditions(); + + // 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(); + Kokkos::parallel_for(ippl::getRangePolicy(aview, 1), KOKKOS_LAMBDA(size_t i, size_t j, size_t k){ + 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; + const scalar magnitude = gauss(scalar(0.5), scalar(0.05), z); + aview (i,j,k) = vector4_type{scalar(0), scalar(0), magnitude, scalar(0)}; + am1view(i,j,k) = vector4_type{scalar(0), scalar(0), magnitude, scalar(0)}; + }); + 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.dt ;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; + 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; + ref += Kokkos::abs(gauss(scalar(0.5), scalar(0.05), z) - aview(i,j,k)[2]); + }, sum_error); + std::cout << "Sum after evaluating boundaryconditions: " << sum_error / double(nr[0] * nr[1] * nr[2]) << "\n"; + + } + ippl::finalize(); +} \ No newline at end of file diff --git a/test/solver/TestNonStandardFDTDSolver_convergence.cpp b/test/solver/TestNonStandardFDTDSolver_convergence.cpp new file mode 100644 index 000000000..dd848d25b --- /dev/null +++ b/test/solver/TestNonStandardFDTDSolver_convergence.cpp @@ -0,0 +1,166 @@ +#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]; + + } + #ifndef __CUDA_ARCH__ + using std::exp; + #endif + return 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') ? 0 : (direction == 'y') ? 1 : 2; + + // 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; 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 NonStandardFDTDSolver object + ippl::NonStandardFDTDSolver sfdsolver(source, E, B); + + // Set boundary conditions + sfdsolver.setPeriodicBoundaryConditions(); + + // 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(); + Kokkos::parallel_for(ippl::getRangePolicy(aview, 1), KOKKOS_LAMBDA(size_t i, size_t j, size_t k){ + // 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 == 0) ? x : (dir == 1) ? 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; + }); + 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.dt ;s++){ + // Solve the FDTD equations + sfdsolver.solve(); + } + + // Calculate the L2 norm error between the computed and expected values + double sum_error = 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 == 0) ? x : (dir == 1) ? y : z; + + // Calculate error in given direction at this point + aview(i,j,k)[dir] -= gauss(scalar(0.5), scalar(0.05), coord); + ref += Kokkos::sqrt(aview(i,j,k)[0] * aview(i,j,k)[0] + aview(i,j,k)[1] * aview(i,j,k)[1] + aview(i,j,k)[2] * aview(i,j,k)[2]); + }, 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 << "," << sum_error / double(nr[1] * nr[1] * nr[2]) << 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/solver/TestStandardFDTDSolver.cpp b/test/solver/TestStandardFDTDSolver.cpp new file mode 100644 index 000000000..2a36fd55d --- /dev/null +++ b/test/solver/TestStandardFDTDSolver.cpp @@ -0,0 +1,162 @@ +#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]; + + } + #ifndef __CUDA_ARCH__ + using std::exp; + #endif + return 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>; + constexpr size_t n = 100; + ippl::Vector nr{n, n, n}; + ippl::NDIndex<3> owned(nr[0], nr[1], nr[2]); + ippl::Vector extents{scalar(1), scalar(1), scalar(1)}; + std::array isParallel; + isParallel.fill(false); + isParallel[2] = true; + + // all parallel layout, standard domain, normal axis order + ippl::FieldLayout<3> layout(MPI_COMM_WORLD, owned, isParallel); + + 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); + + // 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); + + // Set boundary conditions + sfdsolver.setPeriodicBoundaryConditions(); + + // 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(); + Kokkos::parallel_for(ippl::getRangePolicy(aview, 1), KOKKOS_LAMBDA(size_t i, size_t j, size_t k){ + 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; + const scalar magnitude = gauss(scalar(0.5), scalar(0.05), z); + aview (i,j,k) = vector4_type{scalar(0), scalar(0), magnitude, scalar(0)}; + am1view(i,j,k) = vector4_type{scalar(0), scalar(0), magnitude, scalar(0)}; + }); + 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.dt;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; + 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; + ref += Kokkos::abs(gauss(scalar(0.5), scalar(0.05), z) - aview(i,j,k)[2]); + }, sum_error); + std::cout << "Sum after evaluating boundaryconditions: " << sum_error / double(nr[0] * nr[1] * nr[2]) << "\n"; + + } + ippl::finalize(); +} \ No newline at end of file diff --git a/test/solver/TestStandardFDTDSolver_convergence.cpp b/test/solver/TestStandardFDTDSolver_convergence.cpp new file mode 100644 index 000000000..b1bc760e3 --- /dev/null +++ b/test/solver/TestStandardFDTDSolver_convergence.cpp @@ -0,0 +1,166 @@ +#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]; + + } + #ifndef __CUDA_ARCH__ + using std::exp; + #endif + return 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') ? 0 : (direction == 'y') ? 1 : 2; + + // 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); + + // Set boundary conditions + sfdsolver.setPeriodicBoundaryConditions(); + + // 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(); + Kokkos::parallel_for(ippl::getRangePolicy(aview, 1), KOKKOS_LAMBDA(size_t i, size_t j, size_t k){ + // 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 == 0) ? x : (dir == 1) ? 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; + }); + 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.dt ;s++){ + // Solve the FDTD equations + sfdsolver.solve(); + } + + // Calculate the L2 norm error between the computed and expected values + double sum_error = 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 == 0) ? x : (dir == 1) ? y : z; + + // Calculate error in given direction at this point + aview(i,j,k)[dir] -= gauss(scalar(0.5), scalar(0.05), coord); + ref += Kokkos::sqrt(aview(i,j,k)[0] * aview(i,j,k)[0] + aview(i,j,k)[1] * aview(i,j,k)[1] + aview(i,j,k)[2] * aview(i,j,k)[2]); + }, 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 << "," << sum_error / double(nr[1] * nr[1] * nr[2]) << 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 From d6c23f4303dde4002e6f21e04cf5810e2221c1cc Mon Sep 17 00:00:00 2001 From: Anna Date: Thu, 19 Jun 2025 16:19:39 +0200 Subject: [PATCH 4/6] Adress comments in pullrequest --- cmake/Dependencies.cmake | 2 +- src/MaxwellSolvers/FDTDSolverBase.h | 1 - src/MaxwellSolvers/FDTDSolverBase.hpp | 5 +---- src/MaxwellSolvers/NonStandardFDTDSolver.h | 3 +++ test/solver/TestNonStandardFDTDSolver.cpp | 2 ++ test/solver/TestNonStandardFDTDSolver_convergence.cpp | 2 ++ test/solver/TestStandardFDTDSolver.cpp | 2 ++ test/solver/TestStandardFDTDSolver_convergence.cpp | 2 ++ 8 files changed, 13 insertions(+), 6 deletions(-) diff --git a/cmake/Dependencies.cmake b/cmake/Dependencies.cmake index d800a402f..a87953f83 100644 --- a/cmake/Dependencies.cmake +++ b/cmake/Dependencies.cmake @@ -117,5 +117,5 @@ if(IPPL_ENABLE_TESTS) 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 free electron laster.") + message(STATUS "✅ stb_image_write loaded for testing FDTD solver.") endif() \ No newline at end of file diff --git a/src/MaxwellSolvers/FDTDSolverBase.h b/src/MaxwellSolvers/FDTDSolverBase.h index 0f479b25a..197ae00e3 100644 --- a/src/MaxwellSolvers/FDTDSolverBase.h +++ b/src/MaxwellSolvers/FDTDSolverBase.h @@ -62,7 +62,6 @@ namespace ippl { SourceField A_np1; // Field at the next time step. SourceField A_nm1; // Field at the previous time step. - // TODO make this private again /** * @brief Shifts the time for the fields. */ diff --git a/src/MaxwellSolvers/FDTDSolverBase.hpp b/src/MaxwellSolvers/FDTDSolverBase.hpp index 93386a1cc..1505de56b 100644 --- a/src/MaxwellSolvers/FDTDSolverBase.hpp +++ b/src/MaxwellSolvers/FDTDSolverBase.hpp @@ -62,7 +62,6 @@ namespace ippl { */ template void FDTDSolverBase::timeShift() { - // TODO: should this comment stay? // 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()); @@ -110,9 +109,7 @@ namespace ippl { 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; - // TODO make this work with tail!!! - //ippl::Vector dAdt = - //(A_n(i, j, k).template tail<3>() - A_nm1(i, j, k).template tail<3>()) * idt; + ippl::Vector dAdx = (A_n(i + 1, j, k) - A_n(i - 1, j, k)) * inverse_2_spacing[0]; ippl::Vector dAdy = diff --git a/src/MaxwellSolvers/NonStandardFDTDSolver.h b/src/MaxwellSolvers/NonStandardFDTDSolver.h index 08ae1ed4f..023f14f0d 100644 --- a/src/MaxwellSolvers/NonStandardFDTDSolver.h +++ b/src/MaxwellSolvers/NonStandardFDTDSolver.h @@ -1,6 +1,9 @@ /** * @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 diff --git a/test/solver/TestNonStandardFDTDSolver.cpp b/test/solver/TestNonStandardFDTDSolver.cpp index 5aac6a7ad..ee8be3cf6 100644 --- a/test/solver/TestNonStandardFDTDSolver.cpp +++ b/test/solver/TestNonStandardFDTDSolver.cpp @@ -1,3 +1,5 @@ +// Check the README.md file in this directory for information about the test and how to run it. + #include using std::size_t; #include diff --git a/test/solver/TestNonStandardFDTDSolver_convergence.cpp b/test/solver/TestNonStandardFDTDSolver_convergence.cpp index dd848d25b..842cfce95 100644 --- a/test/solver/TestNonStandardFDTDSolver_convergence.cpp +++ b/test/solver/TestNonStandardFDTDSolver_convergence.cpp @@ -1,3 +1,5 @@ +// Check the README.md file in this directory for information about the test and how to run it. + #include using std::size_t; #include diff --git a/test/solver/TestStandardFDTDSolver.cpp b/test/solver/TestStandardFDTDSolver.cpp index 2a36fd55d..69db6a5e9 100644 --- a/test/solver/TestStandardFDTDSolver.cpp +++ b/test/solver/TestStandardFDTDSolver.cpp @@ -1,3 +1,5 @@ +// Check the README.md file in this directory for information about the test and how to run it. + #include using std::size_t; #include diff --git a/test/solver/TestStandardFDTDSolver_convergence.cpp b/test/solver/TestStandardFDTDSolver_convergence.cpp index b1bc760e3..2eaed6a7f 100644 --- a/test/solver/TestStandardFDTDSolver_convergence.cpp +++ b/test/solver/TestStandardFDTDSolver_convergence.cpp @@ -1,3 +1,5 @@ +// Check the README.md file in this directory for information about the test and how to run it. + #include using std::size_t; #include From 6d6c01e8396a9d5db84b1dd097a2df8f067c31fa Mon Sep 17 00:00:00 2001 From: Anna Date: Tue, 24 Jun 2025 15:08:30 +0200 Subject: [PATCH 5/6] Adress Sonalis comment to pull request --- src/MaxwellSolvers/AbsorbingBC.h | 768 ++++++++++-------- src/MaxwellSolvers/FDTDSolverBase.h | 49 +- src/MaxwellSolvers/FDTDSolverBase.hpp | 48 +- src/MaxwellSolvers/NonStandardFDTDSolver.h | 34 +- src/MaxwellSolvers/NonStandardFDTDSolver.hpp | 85 +- src/MaxwellSolvers/StandardFDTDSolver.h | 25 +- src/MaxwellSolvers/StandardFDTDSolver.hpp | 42 +- test/solver/TestNonStandardFDTDSolver.cpp | 212 +++-- .../TestNonStandardFDTDSolver_convergence.cpp | 170 ++-- test/solver/TestStandardFDTDSolver.cpp | 211 +++-- .../TestStandardFDTDSolver_convergence.cpp | 166 ++-- 11 files changed, 1030 insertions(+), 780 deletions(-) diff --git a/src/MaxwellSolvers/AbsorbingBC.h b/src/MaxwellSolvers/AbsorbingBC.h index 3678e3f32..b4f8156f1 100644 --- a/src/MaxwellSolvers/AbsorbingBC.h +++ b/src/MaxwellSolvers/AbsorbingBC.h @@ -1,21 +1,23 @@ #ifndef ABC_H #define ABC_H #include + #include "Types/Vector.h" + #include "Field/Field.h" -template -constexpr KOKKOS_INLINE_FUNCTION auto first(){ +template +constexpr KOKKOS_INLINE_FUNCTION auto first() { return a; } -template -constexpr KOKKOS_INLINE_FUNCTION auto second(){ +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. + * @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 @@ -28,13 +30,14 @@ constexpr KOKKOS_INLINE_FUNCTION auto second(){ * @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) - +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. @@ -42,31 +45,34 @@ struct second_order_abc_face{ * @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){ + 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 + // 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) - ); + (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); + 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; + 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; } /*! @@ -77,56 +83,66 @@ struct second_order_abc_face{ * @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{ + 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...}; - 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}; - ippl::Vector mainaxis_off = ippl::Vector{(main_axis == 0) * sign, (main_axis == 1) * sign, (main_axis == 2) * sign}; + + // 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 direciton 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 + // 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) - ); + 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; + 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. + * @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). + * 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. @@ -135,40 +151,54 @@ struct second_order_abc_face{ * @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). + * @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 - - +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 + 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)); + 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 ); + 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; + 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; } /*! @@ -179,103 +209,125 @@ struct second_order_abc_edge{ * @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{ + 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) - 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); + // 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) - ); + 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; - } + 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. + * @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. + * 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). + * @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 +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){ + 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[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); } @@ -287,8 +339,10 @@ struct second_order_abc_corner{ * @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{ + 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); @@ -297,72 +351,60 @@ struct second_order_abc_corner{ // 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}, + 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]) - ); + 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]; - } + 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. + * @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. + * 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. + * 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{ - +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. @@ -372,8 +414,9 @@ struct second_order_mur_boundary_conditions{ * @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){ + 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 @@ -385,188 +428,199 @@ struct second_order_mur_boundary_conditions{ 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)) - }; + 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 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::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}); + 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); + } } - 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::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(); } }; diff --git a/src/MaxwellSolvers/FDTDSolverBase.h b/src/MaxwellSolvers/FDTDSolverBase.h index 197ae00e3..7c9f7451c 100644 --- a/src/MaxwellSolvers/FDTDSolverBase.h +++ b/src/MaxwellSolvers/FDTDSolverBase.h @@ -19,7 +19,7 @@ namespace ippl { /** * @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). @@ -30,17 +30,17 @@ namespace ippl { 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; + 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. */ @@ -53,26 +53,27 @@ namespace ippl { /** * @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. - + + 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 time for the fields. + * @brief Shifts the saved fields in time. */ void timeShift(); - protected: + + protected: /** * @brief Applies the boundary conditions. */ void applyBCs(); - - public: + + public: /** * @brief Steps the solver forward in time. This is a pure virtual function. */ @@ -81,23 +82,21 @@ namespace ippl { * @brief Evaluates the electric and magnetic fields. */ void evaluate_EB(); - - bool periodic_bc; // Flag indicating if periodic boundary conditions are used. - - 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. /** * @brief Initializes the solver. This is a pure virtual function. */ virtual void initialize() = 0; + + 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 +} // namespace ippl #include "FDTDSolverBase.hpp" #endif \ No newline at end of file diff --git a/src/MaxwellSolvers/FDTDSolverBase.hpp b/src/MaxwellSolvers/FDTDSolverBase.hpp index 1505de56b..6be31a55e 100644 --- a/src/MaxwellSolvers/FDTDSolverBase.hpp +++ b/src/MaxwellSolvers/FDTDSolverBase.hpp @@ -7,23 +7,26 @@ 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) { + 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. + * + * 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() { @@ -34,12 +37,13 @@ namespace ippl { /** * @brief Sets periodic boundary conditions. - * - * Configures the solver to use periodic boundary conditions by setting the appropriate boundary conditions for the fields. + * + * Configures the solver to use periodic boundary conditions by setting the appropriate boundary + * conditions for the fields. */ template - void FDTDSolverBase::setPeriodicBoundaryConditions() { - periodic_bc = true; + 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); @@ -56,9 +60,10 @@ namespace ippl { } /** - * @brief Shifts the time for the fields. - * - * Copies the current field values to the previous time step field and the next time step field values to the current field. + * @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() { @@ -69,7 +74,7 @@ namespace ippl { /** * @brief Applies the boundary conditions. - * + * * Applies the specified boundary conditions (periodic or absorbing) to the fields. */ template @@ -88,21 +93,25 @@ namespace ippl { /** * @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. + * + * 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); + *(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(); + 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; @@ -129,5 +138,4 @@ namespace ippl { Kokkos::fence(); } - -} // namespace ippl \ No newline at end of file +} // namespace ippl \ No newline at end of file diff --git a/src/MaxwellSolvers/NonStandardFDTDSolver.h b/src/MaxwellSolvers/NonStandardFDTDSolver.h index 023f14f0d..f33046b84 100644 --- a/src/MaxwellSolvers/NonStandardFDTDSolver.h +++ b/src/MaxwellSolvers/NonStandardFDTDSolver.h @@ -1,8 +1,10 @@ /** * @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 + * @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) */ @@ -15,8 +17,8 @@ using std::size_t; #include "FieldLayout/FieldLayout.h" #include "MaxwellSolvers/AbsorbingBC.h" -#include "MaxwellSolvers/Maxwell.h" #include "MaxwellSolvers/FDTDSolverBase.h" +#include "MaxwellSolvers/Maxwell.h" #include "Meshes/UniformCartesian.h" #include "Particle/ParticleBase.h" @@ -24,8 +26,9 @@ namespace ippl { /** * @class NonStandardFDTDSolver - * @brief A solver for Maxwell's equations using a non-standard Finite-Difference Time-Domain (FDTD) method. - * + * @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). @@ -36,25 +39,28 @@ namespace ippl { 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. + 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 @@ -75,7 +81,7 @@ namespace ippl { */ void initialize() override; }; -} // namespace ippl +} // namespace ippl #include "NonStandardFDTDSolver.hpp" #endif \ No newline at end of file diff --git a/src/MaxwellSolvers/NonStandardFDTDSolver.hpp b/src/MaxwellSolvers/NonStandardFDTDSolver.hpp index 88892d31a..a5b565ef9 100644 --- a/src/MaxwellSolvers/NonStandardFDTDSolver.hpp +++ b/src/MaxwellSolvers/NonStandardFDTDSolver.hpp @@ -4,23 +4,26 @@ */ 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. + * 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) { + 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) { - std::cerr << "Dispersion-free CFL condition not satisfiable\n"; - std::abort(); + throw IpplException("NonStandardFDTDSolver Constructor", + "Dispersion-free CFL condition not satisfiable\n"); } initialize(); } @@ -41,22 +44,33 @@ namespace ippl { 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)))); + 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))), + .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)), + .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 + 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) { + 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]; @@ -71,26 +85,22 @@ namespace ippl { + (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); + 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(); @@ -125,5 +135,10 @@ namespace ippl { 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(); + } }; -} \ No newline at end of file +} // namespace ippl \ No newline at end of file diff --git a/src/MaxwellSolvers/StandardFDTDSolver.h b/src/MaxwellSolvers/StandardFDTDSolver.h index 7a2f42389..60fd6bb37 100644 --- a/src/MaxwellSolvers/StandardFDTDSolver.h +++ b/src/MaxwellSolvers/StandardFDTDSolver.h @@ -1,6 +1,7 @@ /** * @file StandardFDTDSolver.h - * @brief Defines the StandardFDTDSolver class for solving Maxwell's equations using the FDTD method. + * @brief Defines the StandardFDTDSolver class for solving Maxwell's equations using the FDTD + * method. */ #ifndef IPPL_STANDARD_FDTD_SOLVER_H @@ -12,8 +13,8 @@ using std::size_t; #include "FieldLayout/FieldLayout.h" #include "MaxwellSolvers/AbsorbingBC.h" -#include "MaxwellSolvers/Maxwell.h" #include "MaxwellSolvers/FDTDSolverBase.h" +#include "MaxwellSolvers/Maxwell.h" #include "Meshes/UniformCartesian.h" #include "Particle/ParticleBase.h" @@ -21,8 +22,9 @@ namespace ippl { /** * @class StandardFDTDSolver - * @brief A solver for Maxwell's equations using the Finite-Difference Time-Domain (FDTD) method. - * + * @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). @@ -32,17 +34,20 @@ namespace ippl { 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. + 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. @@ -53,7 +58,7 @@ namespace ippl { */ virtual void initialize() override; }; -} // namespace ippl +} // namespace ippl #include "StandardFDTDSolver.hpp" #endif \ No newline at end of file diff --git a/src/MaxwellSolvers/StandardFDTDSolver.hpp b/src/MaxwellSolvers/StandardFDTDSolver.hpp index b74b55841..beefb3660 100644 --- a/src/MaxwellSolvers/StandardFDTDSolver.hpp +++ b/src/MaxwellSolvers/StandardFDTDSolver.hpp @@ -7,21 +7,24 @@ 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. - * + * + * 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) { + 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. @@ -36,8 +39,10 @@ namespace ippl { 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 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); @@ -45,9 +50,11 @@ namespace ippl { 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 + boundary_conditions == absorbing + ? 1 + : 0; // 1 if absorbing, 0 otherwise, indicates the start index of the field - // Update the field values + // 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) { @@ -66,10 +73,10 @@ namespace ippl { // 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); + + 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; } }); @@ -80,7 +87,7 @@ namespace ippl { /** * @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. */ @@ -109,5 +116,10 @@ namespace ippl { 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(); + } } -} \ No newline at end of file +} // namespace ippl \ No newline at end of file diff --git a/test/solver/TestNonStandardFDTDSolver.cpp b/test/solver/TestNonStandardFDTDSolver.cpp index ee8be3cf6..a6ff6b99a 100644 --- a/test/solver/TestNonStandardFDTDSolver.cpp +++ b/test/solver/TestNonStandardFDTDSolver.cpp @@ -1,57 +1,72 @@ +// 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 +template requires((std::is_floating_point_v)) -KOKKOS_INLINE_FUNCTION float gauss(scalar1 mean, scalar1 stddev, scalar... x){ +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++){ + for (unsigned d = 0; d < dim; d++) { vecsum += vec[d] * vec[d]; - } - #ifndef __CUDA_ARCH__ - using std::exp; - #endif - return exp(-(vecsum) / (stddev * stddev)); + return Kokkos::exp(-(vecsum) / (stddev * stddev)); } -int main(int argc, char* argv[]){ +int main(int argc, char* argv[]) { ippl::initialize(argc, argv); { - using scalar = float; + using scalar = float; const unsigned dim = 3; - using vector_type = ippl::Vector; + 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>; + + 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]); - ippl::Vector extents{scalar(1), scalar(1), scalar(1)}; - std::array isParallel; - isParallel.fill(false); - isParallel[2] = true; - // all parallel layout, standard domain, normal axis order - ippl::FieldLayout<3> layout(MPI_COMM_WORLD, owned, isParallel); + // 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++){ + for (unsigned d = 0; d < 3; d++) { hx[d] = extents[d] / (scalar)nr[d]; } - ippl::Vector origin = {0,0,0}; + 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); @@ -61,29 +76,45 @@ int main(int argc, char* argv[]){ // Create the NonStandardFDTDSolver object ippl::NonStandardFDTDSolver sfdsolver(source, E, B); - // Set boundary conditions - sfdsolver.setPeriodicBoundaryConditions(); - - // 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(); + // 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(); - Kokkos::parallel_for(ippl::getRangePolicy(aview, 1), KOKKOS_LAMBDA(size_t i, size_t j, size_t k){ - 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; - const scalar magnitude = gauss(scalar(0.5), scalar(0.05), z); - aview (i,j,k) = vector4_type{scalar(0), scalar(0), magnitude, scalar(0)}; - am1view(i,j,k) = vector4_type{scalar(0), scalar(0), magnitude, scalar(0)}; - }); + + // 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 @@ -94,49 +125,56 @@ int main(int argc, char* argv[]){ sfdsolver.A_nm1.fillHalo(); // Create a 500x500 image with 3 channels (RGB) - int img_width = 500; - int img_height = 500; + 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.dt ;s++){ + + // 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.dt; 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++){ - + 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 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]){ - + 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]); - + 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; + 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}; + 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) + if (s % 4 == 0) stbi_write_bmp(output, img_width, img_height, 3, imagedata_final); // Solve the FDTD equations @@ -144,21 +182,37 @@ int main(int argc, char* argv[]){ } // Calculate the error between the computed and expected values - double sum_error = 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; - ref += Kokkos::abs(gauss(scalar(0.5), scalar(0.05), z) - aview(i,j,k)[2]); - }, sum_error); - std::cout << "Sum after evaluating boundaryconditions: " << sum_error / double(nr[0] * nr[1] * nr[2]) << "\n"; - + 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/solver/TestNonStandardFDTDSolver_convergence.cpp b/test/solver/TestNonStandardFDTDSolver_convergence.cpp index 842cfce95..f03798312 100644 --- a/test/solver/TestNonStandardFDTDSolver_convergence.cpp +++ b/test/solver/TestNonStandardFDTDSolver_convergence.cpp @@ -1,51 +1,55 @@ +// 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 +template requires((std::is_floating_point_v)) -KOKKOS_INLINE_FUNCTION float gauss(scalar1 mean, scalar1 stddev, scalar... x){ +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++){ + for (unsigned d = 0; d < dim; d++) { vecsum += vec[d] * vec[d]; - } - #ifndef __CUDA_ARCH__ - using std::exp; - #endif - return exp(-(vecsum) / (stddev * stddev)); + return Kokkos::exp(-(vecsum) / (stddev * stddev)); } void compute_convergence(char direction, unsigned int np, std::string fname) { - using scalar = float; + using scalar = float; const unsigned dim = 3; - using vector_type = ippl::Vector; + 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') ? 0 : (direction == 'y') ? 1 : 2; + 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) + // 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; here all dimensions are parallel + // specifies decomposition; std::array isParallel; isParallel.fill(true); // unit box ippl::Vector hx; - for(unsigned d = 0;d < 3;d++){ + for (unsigned d = 0; d < 3; d++) { hx[d] = 1 / (scalar)nr[d]; } ippl::Vector origin = {0.0, 0.0, 0.0}; @@ -63,40 +67,47 @@ void compute_convergence(char direction, unsigned int np, std::string fname) { // Create the NonStandardFDTDSolver object ippl::NonStandardFDTDSolver sfdsolver(source, E, B); - // Set boundary conditions - sfdsolver.setPeriodicBoundaryConditions(); - - // 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(); + // 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(); - Kokkos::parallel_for(ippl::getRangePolicy(aview, 1), KOKKOS_LAMBDA(size_t i, size_t j, size_t k){ - // 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 == 0) ? x : (dir == 1) ? 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; - }); - Kokkos::fence(); + // 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); @@ -105,45 +116,56 @@ void compute_convergence(char direction, unsigned int np, std::string fname) { 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.dt ;s++){ + // 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.dt; s++) { // Solve the FDTD equations sfdsolver.solve(); } // Calculate the L2 norm error between the computed and expected values - double sum_error = 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 == 0) ? x : (dir == 1) ? y : z; - - // Calculate error in given direction at this point - aview(i,j,k)[dir] -= gauss(scalar(0.5), scalar(0.05), coord); - ref += Kokkos::sqrt(aview(i,j,k)[0] * aview(i,j,k)[0] + aview(i,j,k)[1] * aview(i,j,k)[1] + aview(i,j,k)[2] * aview(i,j,k)[2]); - }, sum_error); + 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 << "," << sum_error / double(nr[1] * nr[1] * nr[2]) << endl; + csvout << direction << "," << np << "," << Kokkos::sqrt(sum_error / sum_norm) << endl; } -int main(int argc, char* argv[]){ +int main(int argc, char* argv[]) { ippl::initialize(argc, argv); - { + { // gridsizes to iterate over std::array N = {4, 8, 16, 32, 64, 128, 256}; diff --git a/test/solver/TestStandardFDTDSolver.cpp b/test/solver/TestStandardFDTDSolver.cpp index 69db6a5e9..e02165abb 100644 --- a/test/solver/TestStandardFDTDSolver.cpp +++ b/test/solver/TestStandardFDTDSolver.cpp @@ -1,57 +1,72 @@ +// 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 +template requires((std::is_floating_point_v)) -KOKKOS_INLINE_FUNCTION float gauss(scalar1 mean, scalar1 stddev, scalar... x){ +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++){ + for (unsigned d = 0; d < dim; d++) { vecsum += vec[d] * vec[d]; - } - #ifndef __CUDA_ARCH__ - using std::exp; - #endif - return exp(-(vecsum) / (stddev * stddev)); + return Kokkos::exp(-(vecsum) / (stddev * stddev)); } -int main(int argc, char* argv[]){ +int main(int argc, char* argv[]) { ippl::initialize(argc, argv); { - using scalar = float; + using scalar = float; const unsigned dim = 3; - using vector_type = ippl::Vector; + 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>; + + 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]); - ippl::Vector extents{scalar(1), scalar(1), scalar(1)}; - std::array isParallel; - isParallel.fill(false); - isParallel[2] = true; - // all parallel layout, standard domain, normal axis order - ippl::FieldLayout<3> layout(MPI_COMM_WORLD, owned, isParallel); + // 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++){ + for (unsigned d = 0; d < 3; d++) { hx[d] = extents[d] / (scalar)nr[d]; } - ippl::Vector origin = {0,0,0}; + 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); @@ -61,29 +76,45 @@ int main(int argc, char* argv[]){ // Create the StandardFDTDSolver object ippl::StandardFDTDSolver sfdsolver(source, E, B); - // Set boundary conditions - sfdsolver.setPeriodicBoundaryConditions(); - - // 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(); + // 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(); - Kokkos::parallel_for(ippl::getRangePolicy(aview, 1), KOKKOS_LAMBDA(size_t i, size_t j, size_t k){ - 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; - const scalar magnitude = gauss(scalar(0.5), scalar(0.05), z); - aview (i,j,k) = vector4_type{scalar(0), scalar(0), magnitude, scalar(0)}; - am1view(i,j,k) = vector4_type{scalar(0), scalar(0), magnitude, scalar(0)}; - }); + + // 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 @@ -94,49 +125,56 @@ int main(int argc, char* argv[]){ sfdsolver.A_nm1.fillHalo(); // Create a 500x500 image with 3 channels (RGB) - int img_width = 500; - int img_height = 500; + 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.dt;s++){ + + // 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.dt; 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++){ - + 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 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]){ - + 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]); - + 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; + 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}; + 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) + if (s % 4 == 0) stbi_write_bmp(output, img_width, img_height, 3, imagedata_final); // Solve the FDTD equations @@ -144,21 +182,36 @@ int main(int argc, char* argv[]){ } // Calculate the error between the computed and expected values - double sum_error = 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; - ref += Kokkos::abs(gauss(scalar(0.5), scalar(0.05), z) - aview(i,j,k)[2]); - }, sum_error); - std::cout << "Sum after evaluating boundaryconditions: " << sum_error / double(nr[0] * nr[1] * nr[2]) << "\n"; - + 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/solver/TestStandardFDTDSolver_convergence.cpp b/test/solver/TestStandardFDTDSolver_convergence.cpp index 2eaed6a7f..ce50bc2f1 100644 --- a/test/solver/TestStandardFDTDSolver_convergence.cpp +++ b/test/solver/TestStandardFDTDSolver_convergence.cpp @@ -1,41 +1,45 @@ +// 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 +template requires((std::is_floating_point_v)) -KOKKOS_INLINE_FUNCTION float gauss(scalar1 mean, scalar1 stddev, scalar... x){ +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++){ + for (unsigned d = 0; d < dim; d++) { vecsum += vec[d] * vec[d]; - } - #ifndef __CUDA_ARCH__ - using std::exp; - #endif - return exp(-(vecsum) / (stddev * stddev)); + return Kokkos::exp(-(vecsum) / (stddev * stddev)); } void compute_convergence(char direction, unsigned int np, std::string fname) { - using scalar = float; + using scalar = float; const unsigned dim = 3; - using vector_type = ippl::Vector; + 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>; + + 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') ? 0 : (direction == 'y') ? 1 : 2; + 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) + // 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]); @@ -45,7 +49,7 @@ void compute_convergence(char direction, unsigned int np, std::string fname) { // unit box ippl::Vector hx; - for(unsigned d = 0;d < 3;d++){ + for (unsigned d = 0; d < 3; d++) { hx[d] = 1 / (scalar)nr[d]; } ippl::Vector origin = {0.0, 0.0, 0.0}; @@ -63,40 +67,47 @@ void compute_convergence(char direction, unsigned int np, std::string fname) { // Create the StandardFDTDSolver object ippl::StandardFDTDSolver sfdsolver(source, E, B); - // Set boundary conditions - sfdsolver.setPeriodicBoundaryConditions(); - - // 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(); + // 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(); - Kokkos::parallel_for(ippl::getRangePolicy(aview, 1), KOKKOS_LAMBDA(size_t i, size_t j, size_t k){ - // 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 == 0) ? x : (dir == 1) ? 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; - }); - Kokkos::fence(); + // 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); @@ -105,45 +116,56 @@ void compute_convergence(char direction, unsigned int np, std::string fname) { 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.dt ;s++){ + // 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.dt; s++) { // Solve the FDTD equations sfdsolver.solve(); } // Calculate the L2 norm error between the computed and expected values - double sum_error = 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 == 0) ? x : (dir == 1) ? y : z; - - // Calculate error in given direction at this point - aview(i,j,k)[dir] -= gauss(scalar(0.5), scalar(0.05), coord); - ref += Kokkos::sqrt(aview(i,j,k)[0] * aview(i,j,k)[0] + aview(i,j,k)[1] * aview(i,j,k)[1] + aview(i,j,k)[2] * aview(i,j,k)[2]); - }, sum_error); + 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 << "," << sum_error / double(nr[1] * nr[1] * nr[2]) << endl; + csvout << direction << "," << np << "," << Kokkos::sqrt(sum_error / sum_norm) << endl; } -int main(int argc, char* argv[]){ +int main(int argc, char* argv[]) { ippl::initialize(argc, argv); - { + { // gridsizes to iterate over std::array N = {4, 8, 16, 32, 64, 128, 256}; From 1ec7686bedcbe04877549c9a4b0cbf92b1fab108 Mon Sep 17 00:00:00 2001 From: Anna Date: Tue, 24 Jun 2025 17:38:33 +0200 Subject: [PATCH 6/6] Adress comments in pullrequest --- src/MaxwellSolvers/AbsorbingBC.h | 2 +- src/MaxwellSolvers/FDTDSolverBase.h | 15 +++++++-------- test/CMakeLists.txt | 1 + test/maxwell/CMakeLists.txt | 15 +++++++++++++++ test/{solver => maxwell}/README.md | 2 +- .../TestNonStandardFDTDSolver.cpp | 2 +- .../TestNonStandardFDTDSolver_convergence.cpp | 2 +- .../TestStandardFDTDSolver.cpp | 2 +- .../TestStandardFDTDSolver_convergence.cpp | 2 +- test/solver/CMakeLists.txt | 4 ---- 10 files changed, 29 insertions(+), 18 deletions(-) create mode 100644 test/maxwell/CMakeLists.txt rename test/{solver => maxwell}/README.md (99%) rename test/{solver => maxwell}/TestNonStandardFDTDSolver.cpp (99%) rename test/{solver => maxwell}/TestNonStandardFDTDSolver_convergence.cpp (99%) rename test/{solver => maxwell}/TestStandardFDTDSolver.cpp (99%) rename test/{solver => maxwell}/TestStandardFDTDSolver_convergence.cpp (99%) diff --git a/src/MaxwellSolvers/AbsorbingBC.h b/src/MaxwellSolvers/AbsorbingBC.h index b4f8156f1..81474241d 100644 --- a/src/MaxwellSolvers/AbsorbingBC.h +++ b/src/MaxwellSolvers/AbsorbingBC.h @@ -98,7 +98,7 @@ struct second_order_abc_face { 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 direciton of the main axis, the sign is used to determine the + // 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}; diff --git a/src/MaxwellSolvers/FDTDSolverBase.h b/src/MaxwellSolvers/FDTDSolverBase.h index 7c9f7451c..0a16eb3e3 100644 --- a/src/MaxwellSolvers/FDTDSolverBase.h +++ b/src/MaxwellSolvers/FDTDSolverBase.h @@ -67,13 +67,6 @@ namespace ippl { */ void timeShift(); - protected: - /** - * @brief Applies the boundary conditions. - */ - void applyBCs(); - - public: /** * @brief Steps the solver forward in time. This is a pure virtual function. */ @@ -82,12 +75,18 @@ namespace ippl { * @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. 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/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/solver/README.md b/test/maxwell/README.md similarity index 99% rename from test/solver/README.md rename to test/maxwell/README.md index 94b4cac77..9b7299fcc 100644 --- a/test/solver/README.md +++ b/test/maxwell/README.md @@ -75,7 +75,7 @@ for filename in filenames: 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']["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) diff --git a/test/solver/TestNonStandardFDTDSolver.cpp b/test/maxwell/TestNonStandardFDTDSolver.cpp similarity index 99% rename from test/solver/TestNonStandardFDTDSolver.cpp rename to test/maxwell/TestNonStandardFDTDSolver.cpp index a6ff6b99a..ccb9f413e 100644 --- a/test/solver/TestNonStandardFDTDSolver.cpp +++ b/test/maxwell/TestNonStandardFDTDSolver.cpp @@ -135,7 +135,7 @@ int main(int argc, char* argv[]) { // 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.dt; s++) { + 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++) { diff --git a/test/solver/TestNonStandardFDTDSolver_convergence.cpp b/test/maxwell/TestNonStandardFDTDSolver_convergence.cpp similarity index 99% rename from test/solver/TestNonStandardFDTDSolver_convergence.cpp rename to test/maxwell/TestNonStandardFDTDSolver_convergence.cpp index f03798312..ae9305037 100644 --- a/test/solver/TestNonStandardFDTDSolver_convergence.cpp +++ b/test/maxwell/TestNonStandardFDTDSolver_convergence.cpp @@ -118,7 +118,7 @@ void compute_convergence(char direction, unsigned int np, std::string fname) { // 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.dt; s++) { + for (size_t s = 0; s < 1. / sfdsolver.getDt(); s++) { // Solve the FDTD equations sfdsolver.solve(); } diff --git a/test/solver/TestStandardFDTDSolver.cpp b/test/maxwell/TestStandardFDTDSolver.cpp similarity index 99% rename from test/solver/TestStandardFDTDSolver.cpp rename to test/maxwell/TestStandardFDTDSolver.cpp index e02165abb..1b85f50e9 100644 --- a/test/solver/TestStandardFDTDSolver.cpp +++ b/test/maxwell/TestStandardFDTDSolver.cpp @@ -135,7 +135,7 @@ int main(int argc, char* argv[]) { // 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.dt; s++) { + 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++) { diff --git a/test/solver/TestStandardFDTDSolver_convergence.cpp b/test/maxwell/TestStandardFDTDSolver_convergence.cpp similarity index 99% rename from test/solver/TestStandardFDTDSolver_convergence.cpp rename to test/maxwell/TestStandardFDTDSolver_convergence.cpp index ce50bc2f1..e4a259758 100644 --- a/test/solver/TestStandardFDTDSolver_convergence.cpp +++ b/test/maxwell/TestStandardFDTDSolver_convergence.cpp @@ -118,7 +118,7 @@ void compute_convergence(char direction, unsigned int np, std::string fname) { // 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.dt; s++) { + for (size_t s = 0; s < 1. / sfdsolver.getDt(); s++) { // Solve the FDTD equations sfdsolver.solve(); } diff --git a/test/solver/CMakeLists.txt b/test/solver/CMakeLists.txt index 323d6467c..d095eeb01 100644 --- a/test/solver/CMakeLists.txt +++ b/test/solver/CMakeLists.txt @@ -19,8 +19,4 @@ if(IPPL_ENABLE_FFT) add_ippl_integration_test(TestGaussian_biharmonic LABELS solver integration) add_ippl_integration_test(TestGaussian_hessian LABELS solver integration) add_ippl_integration_test(TestP3MSolver LABELS solver integration) - 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()