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
|
