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
|
