Branch data Line data Source code
1 : : //
2 : : // Class FDTDSolverBase
3 : : // Base class for solvers for Maxwell's equations using the FDTD method
4 : : //
5 : :
6 : : namespace ippl {
7 : :
8 : : /**
9 : : * @brief Constructor for the FDTDSolverBase class.
10 : : *
11 : : * Initializes the solver by setting the source and electromagnetic fields.
12 : : *
13 : : * @param source Reference to the source field.
14 : : * @param E Reference to the electric field.
15 : : * @param B Reference to the magnetic field.
16 : : */
17 : : template <typename EMField, typename SourceField, fdtd_bc boundary_conditions>
18 : 0 : FDTDSolverBase<EMField, SourceField, boundary_conditions>::FDTDSolverBase(SourceField& source,
19 : : EMField& E,
20 [ # # # # : 0 : EMField& B) {
# # # # #
# ]
21 : 0 : Maxwell<EMField, SourceField>::setSources(source);
22 : 0 : Maxwell<EMField, SourceField>::setEMFields(E, B);
23 : 0 : }
24 : :
25 : : /**
26 : : * @brief Solves the FDTD equations.
27 : : *
28 : : * Advances the simulation by one time step, shifts the time for the fields, and evaluates the
29 : : * electric and magnetic fields at the new time.
30 : : */
31 : : template <typename EMField, typename SourceField, fdtd_bc boundary_conditions>
32 : 0 : void FDTDSolverBase<EMField, SourceField, boundary_conditions>::solve() {
33 : 0 : step();
34 : 0 : timeShift();
35 : 0 : evaluate_EB();
36 : 0 : }
37 : :
38 : : /**
39 : : * @brief Sets periodic boundary conditions.
40 : : *
41 : : * Configures the solver to use periodic boundary conditions by setting the appropriate boundary
42 : : * conditions for the fields.
43 : : */
44 : : template <typename EMField, typename SourceField, fdtd_bc boundary_conditions>
45 : : void
46 : 0 : FDTDSolverBase<EMField, SourceField, boundary_conditions>::setPeriodicBoundaryConditions() {
47 : 0 : typename SourceField::BConds_t vector_bcs;
48 : 0 : auto bcsetter_single = [&vector_bcs]<size_t Idx>(const std::index_sequence<Idx>&) {
49 [ # # # # : 0 : vector_bcs[Idx] = std::make_shared<ippl::PeriodicFace<SourceField>>(Idx);
# # # # #
# # # ]
50 : 0 : return 0;
51 : : };
52 : 0 : auto bcsetter = [bcsetter_single]<size_t... Idx>(const std::index_sequence<Idx...>&) {
53 [ # # # # : 0 : int x = (bcsetter_single(std::index_sequence<Idx>{}) ^ ...);
# # # # #
# # # ]
54 : : (void)x;
55 : : };
56 [ # # ]: 0 : bcsetter(std::make_index_sequence<Dim * 2>{});
57 [ # # ]: 0 : A_n.setFieldBC(vector_bcs);
58 [ # # ]: 0 : A_np1.setFieldBC(vector_bcs);
59 [ # # ]: 0 : A_nm1.setFieldBC(vector_bcs);
60 : 0 : }
61 : :
62 : : /**
63 : : * @brief Shifts the saved fields in time.
64 : : *
65 : : * Copies the current field values to the previous time step field and the next time step field
66 : : * values to the current field.
67 : : */
68 : : template <typename EMField, typename SourceField, fdtd_bc boundary_conditions>
69 : 0 : void FDTDSolverBase<EMField, SourceField, boundary_conditions>::timeShift() {
70 : : // Look into this, maybe cyclic swap is better
71 : 0 : Kokkos::deep_copy(this->A_nm1.getView(), this->A_n.getView());
72 : 0 : Kokkos::deep_copy(this->A_n.getView(), this->A_np1.getView());
73 : 0 : }
74 : :
75 : : /**
76 : : * @brief Applies the boundary conditions.
77 : : *
78 : : * Applies the specified boundary conditions (periodic or absorbing) to the fields.
79 : : */
80 : : template <typename EMField, typename SourceField, fdtd_bc boundary_conditions>
81 : 0 : void FDTDSolverBase<EMField, SourceField, boundary_conditions>::applyBCs() {
82 : : if constexpr (boundary_conditions == periodic) {
83 : 0 : A_n.getFieldBC().apply(A_n);
84 : 0 : A_nm1.getFieldBC().apply(A_nm1);
85 : 0 : A_np1.getFieldBC().apply(A_np1);
86 : : } else {
87 : : Vector<uint32_t, Dim> true_nr = nr_m;
88 : : true_nr += (A_n.getNghost() * 2);
89 : : second_order_mur_boundary_conditions bcs{};
90 : : bcs.apply(A_n, A_nm1, A_np1, this->dt, true_nr, layout_mp->getLocalNDIndex());
91 : : }
92 : 0 : }
93 : :
94 : : /**
95 : : * @brief Evaluates the electric and magnetic fields.
96 : : *
97 : : * Computes the electric and magnetic fields based on the current, previous, and next time step
98 : : * field values, as well as the source field.
99 : : */
100 : : template <typename EMField, typename SourceField, fdtd_bc boundary_conditions>
101 : 0 : void FDTDSolverBase<EMField, SourceField, boundary_conditions>::evaluate_EB() {
102 [ # # ]: 0 : *(Maxwell<EMField, SourceField>::En_mp) = typename EMField::value_type(0);
103 [ # # ]: 0 : *(Maxwell<EMField, SourceField>::Bn_mp) = typename EMField::value_type(0);
104 [ # # # # ]: 0 : ippl::Vector<scalar, 3> inverse_2_spacing = ippl::Vector<scalar, 3>(0.5) / hr_m;
105 : 0 : const scalar idt = scalar(1.0) / dt;
106 [ # # # # ]: 0 : auto A_np1 = this->A_np1.getView(), A_n = this->A_n.getView(),
107 [ # # ]: 0 : A_nm1 = this->A_nm1.getView();
108 [ # # ]: 0 : auto source = Maxwell<EMField, SourceField>::JN_mp->getView();
109 [ # # ]: 0 : auto Eview = Maxwell<EMField, SourceField>::En_mp->getView();
110 [ # # ]: 0 : auto Bview = Maxwell<EMField, SourceField>::Bn_mp->getView();
111 : :
112 : : // Calculate the electric and magnetic fields
113 : : // Curl and grad of ippl are not used here, because we have a 4-component vector and we need
114 : : // it for the last three components
115 [ # # ]: 0 : Kokkos::parallel_for(
116 [ # # # # : 0 : this->A_n.getFieldRangePolicy(), KOKKOS_LAMBDA(size_t i, size_t j, size_t k) {
# # # # #
# # # # #
# # # # ]
117 [ # # ]: 0 : ippl::Vector<scalar, 3> dAdt;
118 : 0 : dAdt[0] = (A_np1(i, j, k)[1] - A_n(i, j, k)[1]) * idt;
119 : 0 : dAdt[1] = (A_np1(i, j, k)[2] - A_n(i, j, k)[2]) * idt;
120 : 0 : dAdt[2] = (A_np1(i, j, k)[3] - A_n(i, j, k)[3]) * idt;
121 : :
122 [ # # ]: 0 : ippl::Vector<scalar, 4> dAdx =
123 [ # # # # ]: 0 : (A_n(i + 1, j, k) - A_n(i - 1, j, k)) * inverse_2_spacing[0];
124 [ # # ]: 0 : ippl::Vector<scalar, 4> dAdy =
125 [ # # # # ]: 0 : (A_n(i, j + 1, k) - A_n(i, j - 1, k)) * inverse_2_spacing[1];
126 [ # # ]: 0 : ippl::Vector<scalar, 4> dAdz =
127 [ # # # # ]: 0 : (A_n(i, j, k + 1) - A_n(i, j, k - 1)) * inverse_2_spacing[2];
128 : :
129 : 0 : ippl::Vector<scalar, 3> grad_phi{dAdx[0], dAdy[0], dAdz[0]};
130 : 0 : ippl::Vector<scalar, 3> curlA{
131 : 0 : dAdy[3] - dAdz[2],
132 : 0 : dAdz[1] - dAdx[3],
133 : 0 : dAdx[2] - dAdy[1],
134 : : };
135 [ # # # # ]: 0 : Eview(i, j, k) = -dAdt - grad_phi;
136 : 0 : Bview(i, j, k) = curlA;
137 : 0 : });
138 [ # # # # ]: 0 : Kokkos::fence();
139 : 0 : }
140 : :
141 : : } // namespace ippl
|
