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
|
