Branch data Line data Source code
1 : : //
2 : : // Class HaloCells
3 : : // The guard / ghost cells of BareField.
4 : : //
5 : :
6 : : #include <memory>
7 : : #include <vector>
8 : :
9 : : #include "Utility/IpplException.h"
10 : :
11 : : #include "Communicate/Communicator.h"
12 : :
13 : : namespace ippl {
14 : : namespace detail {
15 : : template <typename T, unsigned Dim, class... ViewArgs>
16 : 736 : HaloCells<T, Dim, ViewArgs...>::HaloCells() {}
17 : :
18 : : template <typename T, unsigned Dim, class... ViewArgs>
19 : 0 : void HaloCells<T, Dim, ViewArgs...>::accumulateHalo(view_type& view, Layout_t* layout) {
20 : 0 : exchangeBoundaries<lhs_plus_assign>(view, layout, HALO_TO_INTERNAL);
21 : 0 : }
22 : :
23 : : template <typename T, unsigned Dim, class... ViewArgs>
24 : 0 : void HaloCells<T, Dim, ViewArgs...>::fillHalo(view_type& view, Layout_t* layout) {
25 : 0 : exchangeBoundaries<assign>(view, layout, INTERNAL_TO_HALO);
26 : 0 : }
27 : :
28 : : template <typename T, unsigned Dim, class... ViewArgs>
29 : : template <class Op>
30 : 0 : void HaloCells<T, Dim, ViewArgs...>::exchangeBoundaries(view_type& view, Layout_t* layout,
31 : : SendOrder order) {
32 : : using neighbor_list = typename Layout_t::neighbor_list;
33 : : using range_list = typename Layout_t::neighbor_range_list;
34 : :
35 : 0 : auto& comm = layout->comm;
36 : :
37 : 0 : const neighbor_list& neighbors = layout->getNeighbors();
38 : 0 : const range_list &sendRanges = layout->getNeighborsSendRange(),
39 : 0 : &recvRanges = layout->getNeighborsRecvRange();
40 : :
41 : 0 : size_t totalRequests = 0;
42 [ # # ]: 0 : for (const auto& componentNeighbors : neighbors) {
43 : 0 : totalRequests += componentNeighbors.size();
44 : : }
45 : :
46 : : using memory_space = typename view_type::memory_space;
47 : : using buffer_type = mpi::Communicator::buffer_type<memory_space>;
48 [ # # ]: 0 : std::vector<MPI_Request> requests(totalRequests);
49 : :
50 : : // sending loop
51 : 0 : constexpr size_t cubeCount = detail::countHypercubes(Dim) - 1;
52 : 0 : size_t requestIndex = 0;
53 [ # # ]: 0 : for (size_t index = 0; index < cubeCount; index++) {
54 : 0 : int tag = mpi::tag::HALO + index;
55 : 0 : const auto& componentNeighbors = neighbors[index];
56 [ # # ]: 0 : for (size_t i = 0; i < componentNeighbors.size(); i++) {
57 : 0 : int targetRank = componentNeighbors[i];
58 : :
59 : : bound_type range;
60 [ # # ]: 0 : if (order == INTERNAL_TO_HALO) {
61 : : /*We store only the sending and receiving ranges
62 : : * of INTERNAL_TO_HALO and use the fact that the
63 : : * sending range of HALO_TO_INTERNAL is the receiving
64 : : * range of INTERNAL_TO_HALO and vice versa
65 : : */
66 : 0 : range = sendRanges[index][i];
67 : : } else {
68 : 0 : range = recvRanges[index][i];
69 : : }
70 : :
71 : : size_type nsends;
72 [ # # ]: 0 : pack(range, view, haloData_m, nsends);
73 : :
74 [ # # ]: 0 : buffer_type buf = comm.template getBuffer<memory_space, T>(nsends);
75 : :
76 [ # # ]: 0 : comm.isend(targetRank, tag, haloData_m, *buf, requests[requestIndex++], nsends);
77 : 0 : buf->resetWritePos();
78 : : }
79 : : }
80 : :
81 : : // receiving loop
82 [ # # ]: 0 : for (size_t index = 0; index < cubeCount; index++) {
83 : 0 : int tag = mpi::tag::HALO + Layout_t::getMatchingIndex(index);
84 : 0 : const auto& componentNeighbors = neighbors[index];
85 [ # # ]: 0 : for (size_t i = 0; i < componentNeighbors.size(); i++) {
86 : 0 : int sourceRank = componentNeighbors[i];
87 : :
88 : : bound_type range;
89 [ # # ]: 0 : if (order == INTERNAL_TO_HALO) {
90 : 0 : range = recvRanges[index][i];
91 : : } else {
92 : 0 : range = sendRanges[index][i];
93 : : }
94 : :
95 : 0 : size_type nrecvs = range.size();
96 : :
97 [ # # ]: 0 : buffer_type buf = comm.template getBuffer<memory_space, T>(nrecvs);
98 : :
99 [ # # ]: 0 : comm.recv(sourceRank, tag, haloData_m, *buf, nrecvs * sizeof(T), nrecvs);
100 : 0 : buf->resetReadPos();
101 : :
102 [ # # ]: 0 : unpack<Op>(range, view, haloData_m);
103 : : }
104 : : }
105 : :
106 [ # # ]: 0 : if (totalRequests > 0) {
107 [ # # ]: 0 : MPI_Waitall(totalRequests, requests.data(), MPI_STATUSES_IGNORE);
108 : : }
109 [ # # ]: 0 : comm.freeAllBuffers();
110 : 0 : }
111 : :
112 : : template <typename T, unsigned Dim, class... ViewArgs>
113 : 0 : void HaloCells<T, Dim, ViewArgs...>::pack(const bound_type& range, const view_type& view,
114 : : databuffer_type& fd, size_type& nsends) {
115 [ # # ]: 0 : auto subview = makeSubview(view, range);
116 : :
117 : 0 : auto& buffer = fd.buffer;
118 : :
119 [ # # ]: 0 : size_t size = subview.size();
120 : 0 : nsends = size;
121 [ # # # # ]: 0 : if (buffer.size() < size) {
[ # # ]
122 : 0 : int overalloc = Comm->getDefaultOverallocation();
123 [ # # ]: 0 : Kokkos::realloc(buffer, size * overalloc);
124 : : }
125 : :
126 : : using index_array_type =
127 : : typename RangePolicy<Dim, typename view_type::execution_space>::index_array_type;
128 [ # # # # ]: 0 : ippl::parallel_for(
129 : 0 : "HaloCells::pack()", getRangePolicy(subview),
130 [ # # # # : 0 : KOKKOS_LAMBDA(const index_array_type& args) {
# # # # ]
131 : 0 : int l = 0;
132 : :
133 [ # # # # : 0 : for (unsigned d1 = 0; d1 < Dim; d1++) {
# # # # #
# # # # #
# # # # #
# # # #
# ][ # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # #
# ][ # # #
# # # # #
# # # # ]
[ # # ]
[ # # # # ]
134 : 0 : int next = args[d1];
135 [ # # # # : 0 : for (unsigned d2 = 0; d2 < d1; d2++) {
# # # # #
# # # # #
# # # # #
# # # #
# ][ # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # #
# ][ # # #
# # # # #
# # # # ]
[ # # ]
[ # # # # ]
136 : 0 : next *= subview.extent(d2);
137 : : }
138 : 0 : l += next;
139 : : }
140 : :
141 : 0 : buffer(l) = apply(subview, args);
142 : : });
143 [ # # # # ]: 0 : Kokkos::fence();
144 : 0 : }
145 : :
146 : : template <typename T, unsigned Dim, class... ViewArgs>
147 : : template <typename Op>
148 : 0 : void HaloCells<T, Dim, ViewArgs...>::unpack(const bound_type& range, const view_type& view,
149 : : databuffer_type& fd) {
150 [ # # ]: 0 : auto subview = makeSubview(view, range);
151 [ # # ]: 0 : auto buffer = fd.buffer;
152 : :
153 : : // 29. November 2020
154 : : // https://stackoverflow.com/questions/3735398/operator-as-template-parameter
155 : : Op op;
156 : :
157 : : using index_array_type =
158 : : typename RangePolicy<Dim, typename view_type::execution_space>::index_array_type;
159 [ # # # # ]: 0 : ippl::parallel_for(
160 : 0 : "HaloCells::unpack()", getRangePolicy(subview),
161 [ # # # # : 0 : KOKKOS_LAMBDA(const index_array_type& args) {
# # # # ]
162 : 0 : int l = 0;
163 : :
164 [ # # # # : 0 : for (unsigned d1 = 0; d1 < Dim; d1++) {
# # # # #
# # # # #
# # # # #
# # # #
# ][ # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # #
# ][ # # #
# # # # #
# # # # ]
[ # # ]
[ # # # # ]
165 : 0 : int next = args[d1];
166 [ # # # # : 0 : for (unsigned d2 = 0; d2 < d1; d2++) {
# # # # #
# # # # #
# # # # #
# # # #
# ][ # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # #
# ][ # # #
# # # # #
# # # # ]
[ # # ]
[ # # # # ]
167 : 0 : next *= subview.extent(d2);
168 : : }
169 : 0 : l += next;
170 : : }
171 : :
172 : 0 : op(apply(subview, args), buffer(l));
173 : : });
174 [ # # # # ]: 0 : Kokkos::fence();
175 : 0 : }
176 : :
177 : : template <typename T, unsigned Dim, class... ViewArgs>
178 : 0 : auto HaloCells<T, Dim, ViewArgs...>::makeSubview(const view_type& view,
179 : : const bound_type& intersect) {
180 : 0 : auto makeSub = [&]<size_t... Idx>(const std::index_sequence<Idx...>&) {
181 : : return Kokkos::subview(view,
182 : 0 : Kokkos::make_pair(intersect.lo[Idx], intersect.hi[Idx])...);
183 : : };
184 [ # # ]: 0 : return makeSub(std::make_index_sequence<Dim>{});
185 : : }
186 : :
187 : : template <typename T, unsigned Dim, class... ViewArgs>
188 : : template <typename Op>
189 : 36 : void HaloCells<T, Dim, ViewArgs...>::applyPeriodicSerialDim(view_type& view,
190 : : const Layout_t* layout,
191 : : const int nghost) {
192 : 36 : int myRank = layout->comm.rank();
193 [ + - ]: 36 : const auto& lDomains = layout->getHostLocalDomains();
194 : 36 : const auto& domain = layout->getDomain();
195 : :
196 : : using exec_space = typename view_type::execution_space;
197 : : using index_type = typename RangePolicy<Dim, exec_space>::index_type;
198 : :
199 : : Kokkos::Array<index_type, Dim> ext, begin, end;
200 : :
201 [ + + ]: 162 : for (size_t i = 0; i < Dim; ++i) {
202 : 126 : ext[i] = view.extent(i);
203 : 126 : begin[i] = 0;
204 : : }
205 : :
206 : : Op op;
207 : :
208 [ + + ]: 162 : for (unsigned d = 0; d < Dim; ++d) {
209 : 126 : end = ext;
210 [ + - ]: 126 : end[d] = nghost;
211 : :
212 [ + - ]: 126 : if (lDomains[myRank][d].length() == domain[d].length()) {
213 : 126 : int N = view.extent(d) - 1;
214 : :
215 : : using index_array_type =
216 : : typename RangePolicy<Dim,
217 : : typename view_type::execution_space>::index_array_type;
218 [ + - + - ]: 126 : ippl::parallel_for(
219 : 126 : "applyPeriodicSerialDim", createRangePolicy<Dim, exec_space>(begin, end),
220 [ + - + - : 7234056 : KOKKOS_LAMBDA(index_array_type & coords) {
- - ]
221 : : // The ghosts are filled starting from the inside
222 : : // of the domain proceeding outwards for both lower
223 : : // and upper faces. The extra brackets and explicit
224 : : // mention
225 : :
226 : : // nghost + i
227 : 3616902 : coords[d] += nghost;
228 : 3616902 : auto&& left = apply(view, coords);
229 : :
230 : : // N - nghost - i
231 : 3616902 : coords[d] = N - coords[d];
232 : 3616902 : auto&& right = apply(view, coords);
233 : :
234 : : // nghost - 1 - i
235 : 3616902 : coords[d] += 2 * nghost - 1 - N;
236 : 3616902 : op(apply(view, coords), right);
237 : :
238 : : // N - (nghost - 1 - i) = N - (nghost - 1) + i
239 : 3616902 : coords[d] = N - coords[d];
240 : 3616902 : op(apply(view, coords), left);
241 : : });
242 : :
243 [ + - + - ]: 252 : Kokkos::fence();
244 : : }
245 : : }
246 : 36 : }
247 : : } // namespace detail
248 : : } // namespace ippl
|
