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