Branch data Line data Source code
1 : : /**
2 : : * @file StandardFDTDSolver.hpp
3 : : * @brief Impelmentation of the StandardFDTDSolver class functions.
4 : : */
5 : :
6 : : namespace ippl {
7 : :
8 : : /**
9 : : * @brief Constructor for the StandardFDTDSolver class.
10 : : *
11 : : * This constructor initializes the StandardFDTDSolver with the given source field and
12 : : * electromagnetic fields and initializes the solver.
13 : : *
14 : : * @param source The source field.
15 : : * @param E The electric field.
16 : : * @param B The magnetic field.
17 : : */
18 : : template <typename EMField, typename SourceField, fdtd_bc boundary_conditions>
19 : 0 : StandardFDTDSolver<EMField, SourceField, boundary_conditions>::StandardFDTDSolver(
20 : : SourceField& source, EMField& E, EMField& B)
21 : 0 : : FDTDSolverBase<EMField, SourceField, boundary_conditions>(source, E, B) {
22 [ # # ]: 0 : initialize();
23 : 0 : }
24 : :
25 : : /**
26 : : * @brief Advances the simulation by one time step.
27 : : *
28 : : * This function updates the fields by one time step using the FDTD method.
29 : : * It calculates the new field values based on the current and previous field values,
30 : : * as well as the source field. The boundary conditions are applied after the update.
31 : : */
32 : : template <typename EMField, typename SourceField, fdtd_bc boundary_conditions>
33 : 0 : void StandardFDTDSolver<EMField, SourceField, boundary_conditions>::step() {
34 [ # # ]: 0 : const auto& ldom = this->layout_mp->getLocalNDIndex();
35 : 0 : const int nghost = this->A_n.getNghost();
36 [ # # ]: 0 : const auto aview = this->A_n.getView();
37 [ # # ]: 0 : const auto anp1view = this->A_np1.getView();
38 [ # # ]: 0 : const auto anm1view = this->A_nm1.getView();
39 [ # # ]: 0 : const auto source_view = Maxwell<EMField, SourceField>::JN_mp->getView();
40 : :
41 : : // Calculate the coefficients for the nondispersive update
42 : 0 : const scalar a1 = scalar(2)
43 : 0 : * (scalar(1) - Kokkos::pow(this->dt / this->hr_m[0], 2)
44 : 0 : - Kokkos::pow(this->dt / this->hr_m[1], 2)
45 : 0 : - Kokkos::pow(this->dt / this->hr_m[2], 2));
46 : 0 : const scalar a2 = Kokkos::pow(this->dt / this->hr_m[0], 2);
47 : 0 : const scalar a4 = Kokkos::pow(this->dt / this->hr_m[1], 2);
48 : 0 : const scalar a6 = Kokkos::pow(this->dt / this->hr_m[2], 2);
49 : 0 : const scalar a8 = Kokkos::pow(this->dt, 2);
50 [ # # ]: 0 : Vector<uint32_t, Dim> true_nr = this->nr_m;
51 : 0 : true_nr += (nghost * 2);
52 : 0 : constexpr uint32_t one_if_absorbing_otherwise_0 =
53 : : boundary_conditions == absorbing
54 : : ? 1
55 : : : 0; // 1 if absorbing, 0 otherwise, indicates the start index of the field
56 : :
57 : : // Update the field values
58 [ # # # # ]: 0 : Kokkos::parallel_for(
59 : 0 : "Source field update", ippl::getRangePolicy(aview, nghost),
60 [ # # # # : 0 : KOKKOS_LAMBDA(const size_t i, const size_t j, const size_t k) {
# # # # #
# # # # #
# # # # ]
61 : : // global indices
62 [ # # ]: 0 : const uint32_t ig = i + ldom.first()[0];
63 [ # # ]: 0 : const uint32_t jg = j + ldom.first()[1];
64 [ # # ]: 0 : const uint32_t kg = k + ldom.first()[2];
65 : : // check if at a boundary of the field
66 : 0 : uint32_t val =
67 : 0 : uint32_t(ig == one_if_absorbing_otherwise_0)
68 [ # # ]: 0 : + (uint32_t(jg == one_if_absorbing_otherwise_0) << 1)
69 [ # # ]: 0 : + (uint32_t(kg == one_if_absorbing_otherwise_0) << 2)
70 [ # # ]: 0 : + (uint32_t(ig == true_nr[0] - one_if_absorbing_otherwise_0 - 1) << 3)
71 [ # # ]: 0 : + (uint32_t(jg == true_nr[1] - one_if_absorbing_otherwise_0 - 1) << 4)
72 [ # # ]: 0 : + (uint32_t(kg == true_nr[2] - one_if_absorbing_otherwise_0 - 1) << 5);
73 : : // update the interior field values
74 [ # # ]: 0 : if (val == 0) {
75 [ # # # # : 0 : SourceVector_t interior = -anm1view(i, j, k) + a1 * aview(i, j, k)
# # ]
76 [ # # # # : 0 : + a2 * (aview(i + 1, j, k) + aview(i - 1, j, k))
# # ]
77 [ # # # # : 0 : + a4 * (aview(i, j + 1, k) + aview(i, j - 1, k))
# # ]
78 [ # # # # : 0 : + a6 * (aview(i, j, k + 1) + aview(i, j, k - 1))
# # ]
79 [ # # # # ]: 0 : + a8 * source_view(i, j, k);
80 : 0 : anp1view(i, j, k) = interior;
81 : 0 : }
82 : : });
83 [ # # # # ]: 0 : Kokkos::fence();
84 [ # # ]: 0 : this->A_np1.fillHalo();
85 [ # # ]: 0 : this->applyBCs();
86 : 0 : }
87 : :
88 : : /**
89 : : * @brief Initializes the solver.
90 : : *
91 : : * This function initializes the solver by setting up the layout, mesh, and fields.
92 : : * It retrieves the mesh spacing, domain, and mesh size, and initializes the fields to zero.
93 : : */
94 : : template <typename EMField, typename SourceField, fdtd_bc boundary_conditions>
95 : 0 : void StandardFDTDSolver<EMField, SourceField, boundary_conditions>::initialize() {
96 : : // get layout and mesh
97 : 0 : this->layout_mp = &(this->JN_mp->getLayout());
98 : 0 : this->mesh_mp = &(this->JN_mp->get_mesh());
99 : :
100 : : // get mesh spacing, timestep, domain, and mesh size
101 : 0 : this->hr_m = this->mesh_mp->getMeshSpacing();
102 : 0 : this->dt = this->hr_m[0] / 2;
103 [ # # ]: 0 : for (unsigned int i = 0; i < Dim; ++i) {
104 : 0 : this->dt = std::min(this->dt, this->hr_m[i] / 2);
105 : : }
106 : 0 : this->domain_m = this->layout_mp->getDomain();
107 [ # # ]: 0 : for (unsigned int i = 0; i < Dim; ++i)
108 : 0 : this->nr_m[i] = this->domain_m[i].length();
109 : :
110 : : // initialize fields
111 : 0 : this->A_nm1.initialize(*(this->mesh_mp), *(this->layout_mp));
112 : 0 : this->A_n.initialize(*(this->mesh_mp), *(this->layout_mp));
113 : 0 : this->A_np1.initialize(*(this->mesh_mp), *(this->layout_mp));
114 : :
115 : : // Initialize fields to zero
116 [ # # ]: 0 : this->A_nm1 = 0.0;
117 [ # # ]: 0 : this->A_n = 0.0;
118 [ # # ]: 0 : this->A_np1 = 0.0;
119 : :
120 : : // set the mesh domain
121 : : if constexpr (boundary_conditions == periodic) {
122 : 0 : this->setPeriodicBoundaryConditions();
123 : : }
124 : 0 : }
125 : : } // namespace ippl
|
