Line data Source code
1 : // This file contains the abstract base class for
2 : // field boundary conditions and other child classes
3 : // which represent specific BCs. At the moment the
4 : // following field BCs are supported
5 : //
6 : // 1. Periodic BC
7 : // 2. Zero BC
8 : // 3. Specifying a constant BC
9 : // 4. No BC (default option)
10 : // 5. Constant extrapolation BC
11 : // Only cell-centered field BCs are implemented
12 : // at the moment.
13 : //
14 :
15 : #include "Utility/IpplException.h"
16 :
17 : #include "Field/HaloCells.h"
18 :
19 : namespace ippl {
20 : namespace detail {
21 :
22 : template <typename Field>
23 16440 : BCondBase<Field>::BCondBase(unsigned int face)
24 16440 : : face_m(face)
25 16440 : , changePhysical_m(false) {}
26 :
27 : template <typename Field>
28 : inline std::ostream& operator<<(std::ostream& os, const BCondBase<Field>& bc) {
29 : bc.write(os);
30 : return os;
31 : }
32 :
33 : } // namespace detail
34 :
35 : template <typename Field>
36 1344 : void ExtrapolateFace<Field>::apply(Field& field) {
37 : // We only support constant extrapolation for the moment, other
38 : // higher order extrapolation stuffs need to be added.
39 :
40 1344 : unsigned int face = this->face_m;
41 1344 : unsigned d = face / 2;
42 1344 : if (field.getCommunicator().size() > 1) {
43 1344 : const Layout_t& layout = field.getLayout();
44 1344 : const auto& lDomains = layout.getHostLocalDomains();
45 1344 : const auto& domain = layout.getDomain();
46 1344 : int myRank = field.getCommunicator().rank();
47 :
48 1344 : bool isBoundary = (lDomains[myRank][d].max() == domain[d].max())
49 1536 : || (lDomains[myRank][d].min() == domain[d].min());
50 :
51 1344 : if (!isBoundary) {
52 0 : return;
53 : }
54 1344 : }
55 :
56 : // If we are here then it is a processor with the face on the physical
57 : // boundary or it is the single core case. Then the following code is same
58 : // irrespective of either it is a single core or multi-core case as the
59 : // non-periodic BC is local to apply.
60 1344 : typename Field::view_type& view = field.getView();
61 1344 : const int nghost = field.getNghost();
62 : int src, dest;
63 :
64 : // It is not clear what it exactly means to do extrapolate
65 : // BC for nghost >1
66 1344 : if (nghost > 1) {
67 0 : throw IpplException("ExtrapolateFace::apply", "nghost > 1 not supported");
68 : }
69 :
70 1344 : if (d >= Dim) {
71 0 : throw IpplException("ExtrapolateFace::apply", "face number wrong");
72 : }
73 :
74 : // If face & 1 is true, then it is an upper BC
75 1344 : if (face & 1) {
76 672 : src = view.extent(d) - 2;
77 672 : dest = src + 1;
78 : } else {
79 672 : src = 1;
80 672 : dest = src - 1;
81 : }
82 :
83 : using exec_space = typename Field::execution_space;
84 : using index_type = typename RangePolicy<Dim, exec_space>::index_type;
85 : Kokkos::Array<index_type, Dim> begin, end;
86 7168 : for (unsigned i = 0; i < Dim; i++) {
87 5824 : begin[i] = nghost;
88 5824 : end[i] = view.extent(i) - nghost;
89 : }
90 1344 : begin[d] = src;
91 1344 : end[d] = src + 1;
92 : using index_array_type = typename RangePolicy<Dim, exec_space>::index_array_type;
93 1344 : ippl::parallel_for(
94 1344 : "Assign extrapolate BC", createRangePolicy<Dim, exec_space>(begin, end),
95 8985344 : KOKKOS_CLASS_LAMBDA(index_array_type & args) {
96 : // to avoid ambiguity with the member function
97 : using ippl::apply;
98 :
99 4491328 : T value = apply(view, args);
100 :
101 4491328 : args[d] = dest;
102 :
103 4491328 : apply(view, args) = slope_m * value + offset_m;
104 0 : });
105 : }
106 :
107 : template <typename Field>
108 0 : void ExtrapolateFace<Field>::write(std::ostream& out) const {
109 : out << "Constant Extrapolation Face"
110 0 : << ", Face = " << this->face_m;
111 0 : }
112 :
113 : template <typename Field>
114 0 : void NoBcFace<Field>::write(std::ostream& out) const {
115 : out << "NoBcFace"
116 0 : << ", Face = " << this->face_m;
117 0 : }
118 :
119 : template <typename Field>
120 0 : void ConstantFace<Field>::write(std::ostream& out) const {
121 : out << "ConstantFace"
122 0 : << ", Face = " << this->face_m << ", Constant = " << this->offset_m;
123 0 : }
124 :
125 : template <typename Field>
126 0 : void ZeroFace<Field>::write(std::ostream& out) const {
127 : out << "ZeroFace"
128 0 : << ", Face = " << this->face_m;
129 0 : }
130 :
131 : template <typename Field>
132 0 : void PeriodicFace<Field>::write(std::ostream& out) const {
133 : out << "PeriodicFace"
134 0 : << ", Face = " << this->face_m;
135 0 : }
136 :
137 : template <typename Field>
138 3360 : void PeriodicFace<Field>::findBCNeighbors(Field& field) {
139 3360 : auto& comm = field.getCommunicator();
140 : // For cell centering only face neighbors are needed
141 3360 : unsigned int face = this->face_m;
142 3360 : unsigned int d = face / 2;
143 3360 : const int nghost = field.getNghost();
144 3360 : int myRank = comm.rank();
145 3360 : const Layout_t& layout = field.getLayout();
146 3360 : const auto& lDomains = layout.getHostLocalDomains();
147 3360 : const auto& domain = layout.getDomain();
148 :
149 32480 : for (auto& neighbors : faceNeighbors_m) {
150 29120 : neighbors.clear();
151 : }
152 :
153 3360 : if (lDomains[myRank][d].length() < domain[d].length()) {
154 : // Only along this dimension we need communication.
155 :
156 960 : bool isBoundary = (lDomains[myRank][d].max() == domain[d].max())
157 1440 : || (lDomains[myRank][d].min() == domain[d].min());
158 :
159 960 : if (isBoundary) {
160 : // this face is on mesh/physical boundary
161 : // get my local box
162 960 : auto& nd = lDomains[myRank];
163 :
164 : // grow the box by nghost cells in dimension d of face
165 960 : auto gnd = nd.grow(nghost, d);
166 :
167 : int offset;
168 960 : if (face & 1) {
169 : // upper face
170 480 : offset = -domain[d].length();
171 : } else {
172 : // lower face
173 480 : offset = domain[d].length();
174 : }
175 : // shift by offset
176 960 : gnd[d] = gnd[d] + offset;
177 :
178 : // Now, we are ready to intersect
179 2880 : for (int rank = 0; rank < comm.size(); ++rank) {
180 1920 : if (rank == myRank) {
181 960 : continue;
182 : }
183 :
184 1920 : if (gnd.touches(lDomains[rank])) {
185 480 : faceNeighbors_m[face].push_back(rank);
186 : }
187 : }
188 : }
189 : }
190 3360 : }
191 :
192 : template <typename Field>
193 1344 : void PeriodicFace<Field>::apply(Field& field) {
194 1344 : auto& comm = field.getCommunicator();
195 1344 : unsigned int face = this->face_m;
196 1344 : unsigned int d = face / 2;
197 1344 : typename Field::view_type& view = field.getView();
198 1344 : const Layout_t& layout = field.getLayout();
199 1344 : const int nghost = field.getNghost();
200 1344 : int myRank = comm.rank();
201 1344 : const auto& lDomains = layout.getHostLocalDomains();
202 1344 : const auto& domain = layout.getDomain();
203 :
204 : // We have to put tag here so that the matchtag inside
205 : // the if is proper.
206 1344 : int tag = comm.next_tag(mpi::tag::BC_PARALLEL_PERIODIC, mpi::tag::BC_CYCLE);
207 :
208 1344 : if (lDomains[myRank][d].length() < domain[d].length()) {
209 : // Only along this dimension we need communication.
210 :
211 384 : bool isBoundary = (lDomains[myRank][d].max() == domain[d].max())
212 576 : || (lDomains[myRank][d].min() == domain[d].min());
213 :
214 384 : if (isBoundary) {
215 : // this face is on mesh/physical boundary
216 : // get my local box
217 384 : auto& nd = lDomains[myRank];
218 :
219 : int offset, offsetRecv, matchtag;
220 384 : if (face & 1) {
221 : // upper face
222 192 : offset = -domain[d].length();
223 192 : offsetRecv = nghost;
224 192 : matchtag = comm.preceding_tag(mpi::tag::BC_PARALLEL_PERIODIC);
225 : } else {
226 : // lower face
227 192 : offset = domain[d].length();
228 192 : offsetRecv = -nghost;
229 192 : matchtag = comm.following_tag(mpi::tag::BC_PARALLEL_PERIODIC);
230 : }
231 :
232 384 : auto& neighbors = faceNeighbors_m[face];
233 :
234 : using memory_space = typename Field::memory_space;
235 : using buffer_type = mpi::Communicator::buffer_type<memory_space>;
236 384 : std::vector<MPI_Request> requests(neighbors.size());
237 :
238 : using HaloCells_t = typename Field::halo_type;
239 : using range_t = typename HaloCells_t::bound_type;
240 384 : HaloCells_t& halo = field.getHalo();
241 384 : std::vector<range_t> rangeNeighbors;
242 :
243 576 : for (size_t i = 0; i < neighbors.size(); ++i) {
244 192 : int rank = neighbors[i];
245 :
246 192 : auto ndNeighbor = lDomains[rank];
247 192 : ndNeighbor[d] = ndNeighbor[d] - offset;
248 :
249 192 : NDIndex<Dim> gndNeighbor = ndNeighbor.grow(nghost, d);
250 :
251 192 : NDIndex<Dim> overlap = gndNeighbor.intersect(nd);
252 :
253 : range_t range;
254 :
255 864 : for (size_t j = 0; j < Dim; ++j) {
256 672 : range.lo[j] = overlap[j].first() - nd[j].first() + nghost;
257 672 : range.hi[j] = overlap[j].last() - nd[j].first() + nghost + 1;
258 : }
259 :
260 192 : rangeNeighbors.push_back(range);
261 :
262 : detail::size_type nSends;
263 192 : halo.pack(range, view, haloData_m, nSends);
264 :
265 192 : buffer_type buf = comm.template getBuffer<memory_space, T>(nSends);
266 :
267 192 : comm.isend(rank, tag, haloData_m, *buf, requests[i], nSends);
268 192 : buf->resetWritePos();
269 : }
270 :
271 576 : for (size_t i = 0; i < neighbors.size(); ++i) {
272 192 : int rank = neighbors[i];
273 :
274 192 : range_t range = rangeNeighbors[i];
275 :
276 192 : range.lo[d] = range.lo[d] + offsetRecv;
277 192 : range.hi[d] = range.hi[d] + offsetRecv;
278 :
279 192 : detail::size_type nRecvs = range.size();
280 :
281 192 : buffer_type buf = comm.template getBuffer<memory_space, T>(nRecvs);
282 192 : comm.recv(rank, matchtag, haloData_m, *buf, nRecvs * sizeof(T), nRecvs);
283 192 : buf->resetReadPos();
284 :
285 : using assign_t = typename HaloCells_t::assign;
286 192 : halo.template unpack<assign_t>(range, view, haloData_m);
287 : }
288 384 : if (!requests.empty()) {
289 192 : MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
290 : }
291 384 : comm.freeAllBuffers();
292 384 : }
293 : // For all other processors do nothing
294 : } else {
295 960 : if (d >= Dim) {
296 0 : throw IpplException("PeriodicFace::apply", "face number wrong");
297 : }
298 :
299 960 : auto N = view.extent(d) - 1;
300 :
301 : using exec_space = typename Field::execution_space;
302 : using index_type = typename RangePolicy<Dim, exec_space>::index_type;
303 : Kokkos::Array<index_type, Dim> begin, end;
304 :
305 : // For the axis along which BCs are being applied, iterate
306 : // through only the ghost cells. For all other axes, iterate
307 : // through all internal cells.
308 5440 : for (size_t i = 0; i < Dim; ++i) {
309 4480 : end[i] = view.extent(i) - nghost;
310 4480 : begin[i] = nghost;
311 : }
312 960 : begin[d] = 0;
313 960 : end[d] = nghost;
314 :
315 : using index_array_type = typename RangePolicy<Dim, exec_space>::index_array_type;
316 960 : ippl::parallel_for(
317 960 : "Assign periodic field BC", createRangePolicy<Dim, exec_space>(begin, end),
318 8507264 : KOKKOS_CLASS_LAMBDA(index_array_type & coords) {
319 : // The ghosts are filled starting from the inside of
320 : // the domain proceeding outwards for both lower and
321 : // upper faces.
322 :
323 : // to avoid ambiguity with the member function
324 : using ippl::apply;
325 :
326 : // x -> nghost + x
327 0 : coords[d] += nghost;
328 0 : auto&& left = apply(view, coords);
329 :
330 : // nghost + x -> N - (nghost + x) = N - nghost - x
331 0 : coords[d] = N - coords[d];
332 0 : auto&& right = apply(view, coords);
333 :
334 : // N - nghost - x -> nghost - 1 - x
335 0 : coords[d] += 2 * nghost - 1 - N;
336 0 : apply(view, coords) = right;
337 :
338 : // nghost - 1 - x -> N - (nghost - 1 - x)
339 : // = N - (nghost - 1) + x
340 0 : coords[d] = N - coords[d];
341 0 : apply(view, coords) = left;
342 : });
343 : }
344 1344 : }
345 :
346 : template <typename Field>
347 0 : void PeriodicFace<Field>::assignGhostToPhysical(Field& field) {
348 0 : unsigned int face = this->face_m;
349 0 : unsigned int d = face / 2;
350 0 : typename Field::view_type& view = field.getView();
351 0 : const Layout_t& layout = field.getLayout();
352 0 : const int nghost = field.getNghost();
353 0 : const auto& ldom = layout.getLocalNDIndex();
354 0 : const auto& domain = layout.getDomain();
355 :
356 0 : if (d >= Dim) {
357 0 : throw IpplException("PeriodicFace::apply", "face number wrong");
358 : }
359 :
360 0 : bool upperFace = (face & 1);
361 0 : bool isBoundary = ((ldom[d].max() == domain[d].max()) && upperFace)
362 0 : || ((ldom[d].min() == domain[d].min()) && !(upperFace));
363 :
364 0 : if (isBoundary) {
365 :
366 0 : auto N = view.extent(d) - 1;
367 :
368 : using exec_space = typename Field::execution_space;
369 : using index_type = typename RangePolicy<Dim, exec_space>::index_type;
370 : Kokkos::Array<index_type, Dim> begin, end;
371 :
372 : // For the axis along which BCs are being applied, iterate
373 : // through only the ghost cells. For all other axes, iterate
374 : // through all internal cells.
375 0 : bool isCorner = (d != 0);
376 0 : for (size_t i = 0; i < Dim; ++i) {
377 0 : bool upperFace_i = (ldom[i].max() == domain[i].max());
378 0 : bool lowerFace_i = (ldom[i].min() == domain[i].min());
379 0 : end[i] = view.extent(i) - nghost - (upperFace_i)*(isCorner);
380 0 : begin[i] = nghost + (lowerFace_i)*(isCorner);
381 : }
382 0 : begin[d] = ((0 + nghost - 1) * (1 - upperFace)) + (N * upperFace);
383 0 : end[d] = begin[d] + 1;
384 :
385 : using index_array_type = typename RangePolicy<Dim, exec_space>::index_array_type;
386 0 : ippl::parallel_for(
387 0 : "Assign periodic field BC", createRangePolicy<Dim, exec_space>(begin, end),
388 0 : KOKKOS_CLASS_LAMBDA(index_array_type & coords) {
389 : // we add the ghost cell values to the appropriate
390 : // neighbouring physical boundary cell
391 :
392 : // to avoid ambiguity with the member function
393 : using ippl::apply;
394 :
395 : // get the value at ghost cells
396 0 : auto&& right = apply(view, coords);
397 :
398 : // apply to the last physical cells (boundary)
399 0 : int shift = 1 - (2 * upperFace);
400 0 : coords[d] += shift;
401 :
402 0 : apply(view, coords) += right;
403 : });
404 : }
405 0 : }
406 :
407 : template <typename Field>
408 0 : void ExtrapolateFace<Field>::assignGhostToPhysical(Field& field) {
409 0 : unsigned int face = this->face_m;
410 0 : unsigned int d = face / 2;
411 0 : typename Field::view_type& view = field.getView();
412 0 : const Layout_t& layout = field.getLayout();
413 0 : const int nghost = field.getNghost();
414 0 : const auto& ldom = layout.getLocalNDIndex();
415 0 : const auto& domain = layout.getDomain();
416 :
417 0 : if (d >= Dim) {
418 0 : throw IpplException("ExtrapolateFace::apply", "face number wrong");
419 : }
420 :
421 0 : bool upperFace = (face & 1);
422 0 : bool isBoundary = ((ldom[d].max() == domain[d].max()) && upperFace)
423 0 : || ((ldom[d].min() == domain[d].min()) && !(upperFace));
424 :
425 0 : if (isBoundary) {
426 0 : auto N = view.extent(d) - 1;
427 :
428 : using exec_space = typename Field::execution_space;
429 : using index_type = typename RangePolicy<Dim, exec_space>::index_type;
430 : Kokkos::Array<index_type, Dim> begin, end;
431 :
432 : // For the axis along which BCs are being applied, iterate
433 : // through only the ghost cells. For all other axes, iterate
434 : // through all internal cells.
435 0 : for (size_t i = 0; i < Dim; ++i) {
436 0 : end[i] = view.extent(i) - nghost;
437 0 : begin[i] = nghost;
438 : }
439 0 : begin[d] = ((0 + nghost - 1) * (1 - upperFace)) + (N * upperFace);
440 0 : end[d] = begin[d] + 1;
441 :
442 : using index_array_type = typename RangePolicy<Dim, exec_space>::index_array_type;
443 0 : ippl::parallel_for(
444 0 : "Assign field BC", createRangePolicy<Dim, exec_space>(begin, end),
445 0 : KOKKOS_CLASS_LAMBDA(index_array_type & coords) {
446 : // we assign the ghost cell values to the appropriate
447 : // neighbouring physical boundary cell
448 :
449 : // to avoid ambiguity with the member function
450 : using ippl::apply;
451 :
452 : // get the value at ghost cells
453 0 : auto&& right = apply(view, coords);
454 :
455 : // apply to the last physical cells (boundary)
456 0 : int shift = 1 - (2 * upperFace);
457 :
458 0 : coords[d] += shift;
459 :
460 0 : apply(view, coords) = right;
461 : });
462 : }
463 0 : }
464 : } // namespace ippl
|
